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             : #include "libmesh/libmesh_config.h"
      20             : 
      21             : #include "libmesh/fourth_error_estimators.h"
      22             : 
      23             : #include "libmesh/enum_error_estimator_type.h"
      24             : 
      25             : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
      26             : 
      27             : // Local Includes
      28             : #include "libmesh/libmesh_common.h"
      29             : #include "libmesh/error_vector.h"
      30             : #include "libmesh/fe_base.h"
      31             : #include "libmesh/libmesh_logging.h"
      32             : #include "libmesh/elem.h"
      33             : #include "libmesh/system.h"
      34             : #include "libmesh/dense_vector.h"
      35             : #include "libmesh/tensor_tools.h"
      36             : #include "libmesh/enum_norm_type.h"
      37             : 
      38             : // C++ includes
      39             : #include <algorithm> // for std::fill
      40             : #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
      41             : #include <cmath>    // for sqrt
      42             : 
      43             : namespace libMesh
      44             : {
      45             : 
      46             : 
      47        2982 : LaplacianErrorEstimator::LaplacianErrorEstimator() :
      48        2982 :   JumpErrorEstimator()
      49             : {
      50        2982 :   error_norm = H2_SEMINORM;
      51        2982 : }
      52             : 
      53             : 
      54             : 
      55             : ErrorEstimatorType
      56           0 : LaplacianErrorEstimator::type() const
      57             : {
      58           0 :   return LAPLACIAN;
      59             : }
      60             : 
      61             : 
      62             : 
      63             : void
      64        5964 : LaplacianErrorEstimator::init_context(FEMContext & c)
      65             : {
      66         168 :   const unsigned int n_vars = c.n_vars();
      67       11928 :   for (unsigned int v=0; v<n_vars; v++)
      68             :     {
      69             :       // Possibly skip this variable
      70        5964 :       if (error_norm.weight(v) == 0.0) continue;
      71             : 
      72             :       // FIXME: Need to generalize this to vector-valued elements. [PB]
      73         168 :       FEBase * side_fe = nullptr;
      74             : 
      75             :       const std::set<unsigned char> & elem_dims =
      76         168 :         c.elem_dimensions();
      77             : 
      78       11928 :       for (const auto & dim : elem_dims)
      79             :         {
      80        5964 :           c.get_side_fe( v, side_fe, dim );
      81             : 
      82             :           // We'll need hessians on both sides for flux jump computation
      83         168 :           side_fe->get_d2phi();
      84             :         }
      85             :     }
      86        5964 : }
      87             : 
      88             : 
      89             : 
      90             : void
      91       36809 : LaplacianErrorEstimator::internal_side_integration ()
      92             : {
      93        6108 :   const Elem & coarse_elem = coarse_context->get_elem();
      94        6108 :   const Elem & fine_elem   = fine_context->get_elem();
      95             : 
      96        3054 :   const DenseVector<Number> & Ucoarse = coarse_context->get_elem_solution();
      97        3054 :   const DenseVector<Number> & Ufine   = fine_context->get_elem_solution();
      98             : 
      99       36809 :   unsigned short dim = fine_elem.dim();
     100             : 
     101        3054 :   FEBase * fe_fine = nullptr;
     102       36809 :   fine_context->get_side_fe( var, fe_fine, dim );
     103             : 
     104        3054 :   FEBase * fe_coarse = nullptr;
     105        3054 :   coarse_context->get_side_fe( var, fe_coarse, dim );
     106             : 
     107        3054 :   Real error = 1.e-30;
     108       36809 :   unsigned int n_qp = fe_fine->n_quadrature_points();
     109             : 
     110        3054 :   const std::vector<std::vector<RealTensor>> & d2phi_coarse = fe_coarse->get_d2phi();
     111        3054 :   const std::vector<std::vector<RealTensor>> & d2phi_fine = fe_fine->get_d2phi();
     112        9359 :   const std::vector<Real> & JxW_face = fe_fine->get_JxW();
     113             : 
     114      110259 :   for (unsigned int qp=0; qp != n_qp; ++qp)
     115             :     {
     116             :       // Calculate solution gradients on fine and coarse elements
     117             :       // at this quadrature point
     118        6094 :       Number laplacian_fine = 0., laplacian_coarse = 0.;
     119             : 
     120        6094 :       const unsigned int n_coarse_dofs = Ucoarse.size();
     121     1070162 :       for (unsigned int i=0; i != n_coarse_dofs; ++i)
     122             :         {
     123     1161972 :           laplacian_coarse += d2phi_coarse[i][qp](0,0) * Ucoarse(i);
     124      996712 :           if (dim > 1)
     125      907008 :             laplacian_coarse += d2phi_coarse[i][qp](1,1) * Ucoarse(i);
     126      996012 :           if (dim > 2)
     127           0 :             laplacian_coarse += d2phi_coarse[i][qp](2,2) * Ucoarse(i);
     128             :         }
     129             : 
     130        6094 :       const unsigned int n_fine_dofs = Ufine.size();
     131     1070162 :       for (unsigned int i=0; i != n_fine_dofs; ++i)
     132             :         {
     133     1161972 :           laplacian_fine += d2phi_fine[i][qp](0,0) * Ufine(i);
     134      996712 :           if (dim > 1)
     135      907008 :             laplacian_fine += d2phi_fine[i][qp](1,1) * Ufine(i);
     136      996012 :           if (dim > 2)
     137           0 :             laplacian_fine += d2phi_fine[i][qp](2,2) * Ufine(i);
     138             :         }
     139             : 
     140             : 
     141             :       // Find the jump in the Laplacian
     142             :       // at this quadrature point
     143       66962 :       const Number jump = laplacian_fine - laplacian_coarse;
     144        6094 :       const Real jump2 = TensorTools::norm_sq(jump);
     145             : 
     146             :       // Accumulate the jump integral
     147       79544 :       error += JxW_face[qp] * jump2;
     148             :     }
     149             : 
     150             :   // Add the h-weighted jump integral to each error term
     151       36809 :   fine_error =
     152       36809 :     error * fine_elem.hmax() * error_norm.weight(var);
     153       36809 :   coarse_error =
     154       36809 :     error * coarse_elem.hmax() * error_norm.weight(var);
     155       36809 : }
     156             : 
     157             : } // namespace libMesh
     158             : 
     159             : #else // defined (LIBMESH_ENABLE_SECOND_DERIVATIVES)
     160             : 
     161             : namespace libMesh
     162             : {
     163             : 
     164             : LaplacianErrorEstimator::LaplacianErrorEstimator() :
     165             :   JumpErrorEstimator()
     166             : {
     167             :   libmesh_error_msg("Error: LaplacianErrorEstimator requires second " \
     168             :                     << "derivative support; try configuring libmesh with " \
     169             :                     << "--enable-second");
     170             : }
     171             : 
     172             : ErrorEstimatorType
     173             : LaplacianErrorEstimator::type() const
     174             : {
     175             :   return LAPLACIAN;
     176             : }
     177             : 
     178             : 
     179             : 
     180             : void LaplacianErrorEstimator::init_context (FEMContext &) {}
     181             : 
     182             : void LaplacianErrorEstimator::internal_side_integration () {}
     183             : 
     184             : } // namespace libMesh
     185             : 
     186             : #endif // defined (LIBMESH_ENABLE_SECOND_DERIVATIVES)