Line data    Source code

       1             : // The libMesh Finite Element Library.
       2             : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
       3             : 
       4             : // This library is free software; you can redistribute it and/or
       5             : // modify it under the terms of the GNU Lesser General Public
       6             : // License as published by the Free Software Foundation; either
       7             : // version 2.1 of the License, or (at your option) any later version.
       8             : 
       9             : // This library is distributed in the hope that it will be useful,
      10             : // but WITHOUT ANY WARRANTY; without even the implied warranty of
      11             : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
      12             : // Lesser General Public License for more details.
      13             : 
      14             : // You should have received a copy of the GNU Lesser General Public
      15             : // License along with this library; if not, write to the Free Software
      16             : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
      17             : 
      18             : 
      19             : // C++ includes
      20             : #include <algorithm> // for std::fill
      21             : #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
      22             : #include <cmath>    // for sqrt
      23             : 
      24             : 
      25             : // Local Includes
      26             : #include "libmesh/libmesh_common.h"
      27             : #include "libmesh/kelly_error_estimator.h"
      28             : #include "libmesh/error_vector.h"
      29             : #include "libmesh/fe_base.h"
      30             : #include "libmesh/libmesh_logging.h"
      31             : #include "libmesh/elem.h"
      32             : #include "libmesh/system.h"
      33             : #include "libmesh/dense_vector.h"
      34             : #include "libmesh/tensor_tools.h"
      35             : #include "libmesh/enum_error_estimator_type.h"
      36             : #include "libmesh/enum_norm_type.h"
      37             : 
      38             : namespace libMesh
      39             : {
      40             : 
      41       10695 : KellyErrorEstimator::KellyErrorEstimator() :
      42             :   JumpErrorEstimator(),
      43       10695 :   _bc_function(nullptr)
      44             : {
      45       10695 :   error_norm = H1_SEMINORM;
      46       10695 : }
      47             : 
      48             : 
      49             : 
      50             : ErrorEstimatorType
      51           0 : KellyErrorEstimator::type() const
      52             : {
      53           0 :   return KELLY;
      54             : }
      55             : 
      56             : 
      57             : 
      58             : void
      59       21410 : KellyErrorEstimator::init_context(FEMContext & c)
      60             : {
      61         804 :   const unsigned int n_vars = c.n_vars();
      62       42820 :   for (unsigned int v=0; v<n_vars; v++)
      63             :     {
      64             :       // Possibly skip this variable
      65       21410 :       if (error_norm.weight(v) == 0.0) continue;
      66             : 
      67             :       // FIXME: Need to generalize this to vector-valued elements. [PB]
      68         804 :       FEBase * side_fe = nullptr;
      69             : 
      70             :       const std::set<unsigned char> & elem_dims =
      71         804 :         c.elem_dimensions();
      72             : 
      73       43104 :       for (const auto & dim : elem_dims)
      74             :         {
      75       21694 :           c.get_side_fe( v, side_fe, dim );
      76             : 
      77             :           // We'll need gradients on both sides for flux jump computation
      78         812 :           side_fe->get_dphi();
      79             : 
      80             :           // But we only need normal vectors from one side
      81       21694 :           if (&c != coarse_context.get())
      82        1013 :             side_fe->get_normals();
      83             :         }
      84             :     }
      85       21410 : }
      86             : 
      87             : 
      88             : 
      89             : void
      90     1676423 : KellyErrorEstimator::internal_side_integration ()
      91             : {
      92      397608 :   const Elem & coarse_elem = coarse_context->get_elem();
      93      397608 :   const Elem & fine_elem = fine_context->get_elem();
      94             : 
      95      198804 :   FEBase * fe_fine = nullptr;
      96     1676423 :   fine_context->get_side_fe( var, fe_fine, fine_elem.dim() );
      97             : 
      98      198804 :   FEBase * fe_coarse = nullptr;
      99     1676423 :   coarse_context->get_side_fe( var, fe_coarse, fine_elem.dim() );
     100             : 
     101      198804 :   Real error = 1.e-30;
     102     1676423 :   unsigned int n_qp = fe_fine->n_quadrature_points();
     103             : 
     104      596068 :   const std::vector<Point> & face_normals = fe_fine->get_normals();
     105      596068 :   const std::vector<Real> & JxW_face = fe_fine->get_JxW();
     106             : 
     107     3536054 :   for (unsigned int qp=0; qp != n_qp; ++qp)
     108             :     {
     109             :       // Calculate solution gradients on fine and coarse elements
     110             :       // at this quadrature point
     111             :       Gradient
     112     1859631 :         grad_fine   = fine_context->side_gradient(var, qp),
     113     1859631 :         grad_coarse = coarse_context->side_gradient(var, qp);
     114             : 
     115             :       // Find the jump in the normal derivative
     116             :       // at this quadrature point
     117     1859631 :       const Number jump = (grad_fine - grad_coarse)*face_normals[qp];
     118      224355 :       const Real jump2 = TensorTools::norm_sq(jump);
     119             : 
     120             :       // Accumulate the jump integral
     121     2083986 :       error += JxW_face[qp] * jump2;
     122             :     }
     123             : 
     124             :   // Add the h-weighted jump integral to each error term
     125     1676423 :   fine_error =
     126     1676423 :     error * fine_elem.hmax() * error_norm.weight(var);
     127     1676423 :   coarse_error =
     128     1676423 :     error * coarse_elem.hmax() * error_norm.weight(var);
     129     1676423 : }
     130             : 
     131             : 
     132             : bool
     133           0 : KellyErrorEstimator::boundary_side_integration ()
     134             : {
     135           0 :   const Elem & fine_elem = fine_context->get_elem();
     136             : 
     137           0 :   FEBase * fe_fine = nullptr;
     138           0 :   fine_context->get_side_fe( var, fe_fine, fine_elem.dim() );
     139             : 
     140             :   const std::string & var_name =
     141           0 :     fine_context->get_system().variable_name(var);
     142             : 
     143           0 :   const std::vector<Point> & face_normals = fe_fine->get_normals();
     144           0 :   const std::vector<Real> & JxW_face = fe_fine->get_JxW();
     145           0 :   const std::vector<Point> & qface_point = fe_fine->get_xyz();
     146             : 
     147             :   // The reinitialization also recomputes the locations of
     148             :   // the quadrature points on the side.  By checking if the
     149             :   // first quadrature point on the side is on a flux boundary
     150             :   // for a particular variable, we will determine if the whole
     151             :   // element is on a flux boundary (assuming quadrature points
     152             :   // are strictly contained in the side).
     153           0 :   if (this->_bc_function(fine_context->get_system(),
     154           0 :                          qface_point[0], var_name).first)
     155             :     {
     156           0 :       const Real h = fine_elem.hmax();
     157             : 
     158             :       // The number of quadrature points
     159           0 :       const unsigned int n_qp = fe_fine->n_quadrature_points();
     160             : 
     161             :       // The error contribution from this face
     162           0 :       Real error = 1.e-30;
     163             : 
     164             :       // loop over the integration points on the face.
     165           0 :       for (unsigned int qp=0; qp<n_qp; qp++)
     166             :         {
     167             :           // Value of the imposed flux BC at this quadrature point.
     168             :           const std::pair<bool,Real> flux_bc =
     169           0 :             this->_bc_function(fine_context->get_system(),
     170           0 :                                qface_point[qp], var_name);
     171             : 
     172             :           // Be sure the BC function still thinks we're on the
     173             :           // flux boundary.
     174           0 :           libmesh_assert_equal_to (flux_bc.first, true);
     175             : 
     176             :           // The solution gradient from each element
     177           0 :           Gradient grad_fine = fine_context->side_gradient(var, qp);
     178             : 
     179             :           // The difference between the desired BC and the approximate solution.
     180           0 :           const Number jump = flux_bc.second - grad_fine*face_normals[qp];
     181             : 
     182             :           // The flux jump squared.  If using complex numbers,
     183             :           // TensorTools::norm_sq(z) returns |z|^2, where |z| is the modulus of z.
     184           0 :           const Real jump2 = TensorTools::norm_sq(jump);
     185             : 
     186             :           // Integrate the error on the face.  The error is
     187             :           // scaled by an additional power of h, where h is
     188             :           // the maximum side length for the element.  This
     189             :           // arises in the definition of the indicator.
     190           0 :           error += JxW_face[qp]*jump2;
     191             : 
     192             :         } // End quadrature point loop
     193             : 
     194           0 :       fine_error = error*h*error_norm.weight(var);
     195             : 
     196           0 :       return true;
     197             :     } // end if side on flux boundary
     198           0 :   return false;
     199             : }
     200             : 
     201             : 
     202             : 
     203             : void
     204           0 : KellyErrorEstimator::attach_flux_bc_function (std::pair<bool,Real> fptr(const System & system,
     205             :                                                                         const Point & p,
     206             :                                                                         const std::string & var_name))
     207             : {
     208           0 :   _bc_function = fptr;
     209             : 
     210             :   // We may be turning boundary side integration on or off
     211           0 :   if (fptr)
     212           0 :     integrate_boundary_sides = true;
     213             :   else
     214           0 :     integrate_boundary_sides = false;
     215           0 : }
     216             : 
     217             : } // namespace libMesh