#include <kelly_error_estimator.h>

Inheritance diagram for libMesh::KellyErrorEstimator:

Public Types
typedef std::map< std::pair< const System , unsigned int >, ErrorVector >	ErrorMap

Public Member Functions
	KellyErrorEstimator ()

	KellyErrorEstimator (const KellyErrorEstimator &)=delete

KellyErrorEstimator &	operator= (const KellyErrorEstimator &)=delete

	KellyErrorEstimator (KellyErrorEstimator &&)=default

KellyErrorEstimator &	operator= (KellyErrorEstimator &&)=default

virtual	~KellyErrorEstimator ()=default

void	attach_flux_bc_function (std::pair< bool, Real > fptr(const System &system, const Point &p, const std::string &var_name))

virtual ErrorEstimatorType	type () const override

virtual void	estimate_error (const System &system, ErrorVector &error_per_cell, const NumericVector< Number > *solution_vector=nullptr, bool estimate_parent_error=false) override

virtual void	estimate_errors (const EquationSystems &equation_systems, ErrorVector &error_per_cell, const std::map< const System , SystemNorm > &error_norms, const std::map< const System , const NumericVector< Number > > solution_vectors=nullptr, bool estimate_parent_error=false)

virtual void	estimate_errors (const EquationSystems &equation_systems, ErrorMap &errors_per_cell, const std::map< const System , const NumericVector< Number > > *solution_vectors=nullptr, bool estimate_parent_error=false)

Public Attributes
bool	scale_by_n_flux_faces

SystemNorm	error_norm

Protected Member Functions
virtual void	init_context (FEMContext &c) override

virtual void	internal_side_integration () override

virtual bool	boundary_side_integration () override

void	reinit_sides ()

float	coarse_n_flux_faces_increment ()

void	reduce_error (std::vector< ErrorVectorReal > &error_per_cell, const Parallel::Communicator &comm) const

Protected Attributes
std::pair< bool, Real >(*	_bc_function )(const System &system, const Point &p, const std::string &var_name)

bool	integrate_boundary_sides

std::unique_ptr< FEMContext >	fine_context

std::unique_ptr< FEMContext >	coarse_context

Real	fine_error

Real	coarse_error

unsigned int	var

Detailed Description

This class implements the Kelly error indicator which is based on the flux jumps between elements. See the JumpErrorEstimator class for most user APIs

Full BibteX reference:

* @Article{Kelly83error,
* author = {D.~W.~Kelly and J.~P.~Gago and O.~C.~Zienkiewicz and I.~Babuska},
* title  = {{A posteriori error analysis and adaptive
*            processes in the finite element method: Part I Error analysis}},
* journal = {Int. J. Num. Meth. Engng.},
* volume  = {19},
* pages   = {1593--1619},
* year    = {1983}
* }
*

Author: Benjamin S. Kirk

Date: 2003

Definition at line 59 of file kelly_error_estimator.h.

Member Typedef Documentation

◆ ErrorMap

typedef std::map<std::pair<const System *, unsigned int>, ErrorVector *> libMesh::ErrorEstimator::ErrorMap

inherited

When calculating many error vectors at once, we need a data structure to hold them all

Definition at line 124 of file error_estimator.h.

Constructor & Destructor Documentation

◆ KellyErrorEstimator() [1/3]

libMesh::KellyErrorEstimator::KellyErrorEstimator ( )

Constructor. Responsible for initializing the _bc_function function pointer to nullptr. Defaults to H1 seminorm; changes to system norm are ignored.

Definition at line 41 of file kelly_error_estimator.C.

References libMesh::ErrorEstimator::error_norm, and libMesh::H1_SEMINORM.

                                          :
   JumpErrorEstimator(),
   _bc_function(nullptr)
 {
   error_norm = H1_SEMINORM;
 }

◆ KellyErrorEstimator() [2/3]

libMesh::KellyErrorEstimator::KellyErrorEstimator ( const KellyErrorEstimator & )

delete

This class cannot be (default) copy constructed/assigned because its base class has unique_ptr members.

◆ KellyErrorEstimator() [3/3]

libMesh::KellyErrorEstimator::KellyErrorEstimator ( KellyErrorEstimator && )

default

Defaulted move ctor, move assignment operator, and destructor.

◆ ~KellyErrorEstimator()

virtual libMesh::KellyErrorEstimator::~KellyErrorEstimator ( )

virtualdefault

Member Function Documentation

◆ attach_flux_bc_function()

void libMesh::KellyErrorEstimator::attach_flux_bc_function ( std::pair< bool, Real > fptrconst System &system, const Point &p, const std::string &var_name )

Register a user function to use in computing the flux BCs.

Definition at line 207 of file kelly_error_estimator.C.

References _bc_function, and libMesh::JumpErrorEstimator::integrate_boundary_sides.

 {
   _bc_function = fptr;
 
   // We may be turning boundary side integration on or off
   if (fptr)
     integrate_boundary_sides = true;
   else
     integrate_boundary_sides = false;
 }

◆ boundary_side_integration()

bool libMesh::KellyErrorEstimator::boundary_side_integration ( )

overrideprotectedvirtual

The function which calculates a normal derivative jump based error term on a boundary side.

Returns: true if the flux bc function is in fact defined on the current side.

Reimplemented from libMesh::JumpErrorEstimator.

Definition at line 135 of file kelly_error_estimator.C.

References _bc_function, libMesh::Elem::dim(), libMesh::ErrorEstimator::error_norm, libMesh::JumpErrorEstimator::fine_context, libMesh::JumpErrorEstimator::fine_error, libMesh::FEGenericBase< OutputType >::get_dphi(), libMesh::FEAbstract::get_JxW(), libMesh::FEAbstract::get_normals(), libMesh::FEAbstract::get_xyz(), libMesh::Elem::hmax(), libMesh::FEAbstract::n_quadrature_points(), libMesh::TensorTools::norm_sq(), libMesh::Real, libMesh::JumpErrorEstimator::var, and libMesh::SystemNorm::weight().

 {
   const Elem & fine_elem = fine_context->get_elem();
 
   FEBase * fe_fine = nullptr;
   fine_context->get_side_fe( var, fe_fine, fine_elem.dim() );
 
   const std::string & var_name =
     fine_context->get_system().variable_name(var);
 
   std::vector<std::vector<RealGradient>> dphi_fine = fe_fine->get_dphi();
   std::vector<Point> face_normals = fe_fine->get_normals();
   std::vector<Real> JxW_face = fe_fine->get_JxW();
   std::vector<Point> qface_point = fe_fine->get_xyz();
 
   // The reinitialization also recomputes the locations of
   // the quadrature points on the side.  By checking if the
   // first quadrature point on the side is on a flux boundary
   // for a particular variable, we will determine if the whole
   // element is on a flux boundary (assuming quadrature points
   // are strictly contained in the side).
   if (this->_bc_function(fine_context->get_system(),
                          qface_point[0], var_name).first)
     {
       const Real h = fine_elem.hmax();
 
       // The number of quadrature points
       const unsigned int n_qp = fe_fine->n_quadrature_points();
 
       // The error contribution from this face
       Real error = 1.e-30;
 
       // loop over the integration points on the face.
       for (unsigned int qp=0; qp<n_qp; qp++)
         {
           // Value of the imposed flux BC at this quadrature point.
           const std::pair<bool,Real> flux_bc =
             this->_bc_function(fine_context->get_system(),
                                qface_point[qp], var_name);
 
           // Be sure the BC function still thinks we're on the
           // flux boundary.
           libmesh_assert_equal_to (flux_bc.first, true);
 
           // The solution gradient from each element
           Gradient grad_fine = fine_context->side_gradient(var, qp);
 
           // The difference between the desired BC and the approximate solution.
           const Number jump = flux_bc.second - grad_fine*face_normals[qp];
 
           // The flux jump squared.  If using complex numbers,
           // TensorTools::norm_sq(z) returns |z|^2, where |z| is the modulus of z.
           const Real jump2 = TensorTools::norm_sq(jump);
 
           // Integrate the error on the face.  The error is
           // scaled by an additional power of h, where h is
           // the maximum side length for the element.  This
           // arises in the definition of the indicator.
           error += JxW_face[qp]*jump2;
 
         } // End quadrature point loop
 
       fine_error = error*h*error_norm.weight(var);
 
       return true;
     } // end if side on flux boundary
   return false;
 }

◆ coarse_n_flux_faces_increment()

float libMesh::JumpErrorEstimator::coarse_n_flux_faces_increment ( )

protectedinherited

A utility function to correctly increase n_flux_faces for the coarse element

Definition at line 462 of file jump_error_estimator.C.

References libMesh::JumpErrorEstimator::coarse_context, and libMesh::JumpErrorEstimator::fine_context.

Referenced by libMesh::JumpErrorEstimator::estimate_error().

 {
   // Keep track of the number of internal flux sides found on each
   // element
   unsigned short dim = coarse_context->get_elem().dim();
 
   const unsigned int divisor =
     1 << (dim-1)*(fine_context->get_elem().level() -
                   coarse_context->get_elem().level());
 
   // With a difference of n levels between fine and coarse elements,
   // we compute a fractional flux face for the coarse element by adding:
   // 1/2^n in 2D
   // 1/4^n in 3D
   // each time.  This code will get hit 2^n times in 2D and 4^n
   // times in 3D so that the final flux face count for the coarse
   // element will be an integer value.
 
   return 1.0f / static_cast<float>(divisor);
 }

◆ estimate_error()

void libMesh::JumpErrorEstimator::estimate_error	(	const System &	system,
		ErrorVector &	error_per_cell,
		const NumericVector< Number > *	solution_vector = `nullptr`,
		bool	estimate_parent_error = `false`
	)

overridevirtualinherited

This function uses the derived class's jump error estimate formula to estimate the error on each cell. The estimated error is output in the vector error_per_cell

Conventions for assigning the direction of the normal:

e & f are global element ids

Case (1.) Elements are at the same level, e<f Compute the flux jump on the face and add it as a contribution to error_per_cell[e] and error_per_cell[f]

| | |

	f
e	—> n

Case (2.) The neighbor is at a higher level. Compute the flux jump on e's face and add it as a contribution to error_per_cell[e] and error_per_cell[f]

| | | | | | e |—> n | | | | | |--------—| f |

Implements libMesh::ErrorEstimator.

Definition at line 53 of file jump_error_estimator.C.

 {
   LOG_SCOPE("estimate_error()", "JumpErrorEstimator");
 
   // This parameter is not used when !LIBMESH_ENABLE_AMR.
   libmesh_ignore(estimate_parent_error);
 
   // The current mesh
   const MeshBase & mesh = system.get_mesh();
 
   // The number of variables in the system
   const unsigned int n_vars = system.n_vars();
 
   // The DofMap for this system
 #ifdef LIBMESH_ENABLE_AMR
   const DofMap & dof_map = system.get_dof_map();
 #endif
 
   // Resize the error_per_cell vector to be
   // the number of elements, initialize it to 0.
   error_per_cell.resize (mesh.max_elem_id());
   std::fill (error_per_cell.begin(), error_per_cell.end(), 0.);
 
   // Declare a vector of floats which is as long as
   // error_per_cell above, and fill with zeros.  This vector will be
   // used to keep track of the number of edges (faces) on each active
   // element which are either:
   // 1) an internal edge
   // 2) an edge on a Neumann boundary for which a boundary condition
   //    function has been specified.
   // The error estimator can be scaled by the number of flux edges (faces)
   // which the element actually has to obtain a more uniform measure
   // of the error.  Use floats instead of ints since in case 2 (above)
   // f gets 1/2 of a flux face contribution from each of his
   // neighbors
   std::vector<float> n_flux_faces;
   if (scale_by_n_flux_faces)
     n_flux_faces.resize(error_per_cell.size(), 0);
 
   // Prepare current_local_solution to localize a non-standard
   // solution vector if necessary
   if (solution_vector && solution_vector != system.solution.get())
     {
       NumericVector<Number> * newsol =
         const_cast<NumericVector<Number> *>(solution_vector);
       System & sys = const_cast<System &>(system);
       newsol->swap(*sys.solution);
       sys.update();
     }
 
   fine_context.reset(new FEMContext(system));
   coarse_context.reset(new FEMContext(system));
 
   // Loop over all the variables we've been requested to find jumps in, to
   // pre-request
   for (var=0; var<n_vars; var++)
     {
       // Skip variables which aren't part of our norm,
       // as well as SCALAR variables, which have no jumps
       if (error_norm.weight(var) == 0.0 ||
           system.variable_type(var).family == SCALAR)
         continue;
 
       // FIXME: Need to generalize this to vector-valued elements. [PB]
       FEBase * side_fe = nullptr;
 
       const std::set<unsigned char> & elem_dims =
         fine_context->elem_dimensions();
 
       for (const auto & dim : elem_dims)
         {
           fine_context->get_side_fe( var, side_fe, dim );
 
           side_fe->get_xyz();
         }
     }
 
   this->init_context(*fine_context);
   this->init_context(*coarse_context);
 
   // Iterate over all the active elements in the mesh
   // that live on this processor.
   for (const auto & e : mesh.active_local_element_ptr_range())
     {
       const dof_id_type e_id = e->id();
 
 #ifdef LIBMESH_ENABLE_AMR
       // See if the parent of element e has been examined yet;
       // if not, we may want to compute the estimator on it
       const Elem * parent = e->parent();
 
       // We only can compute and only need to compute on
       // parents with all active children
       bool compute_on_parent = true;
       if (!parent || !estimate_parent_error)
         compute_on_parent = false;
       else
         for (auto & child : parent->child_ref_range())
           if (!child.active())
             compute_on_parent = false;
 
       if (compute_on_parent &&
           !error_per_cell[parent->id()])
         {
           // Compute a projection onto the parent
           DenseVector<Number> Uparent;
           FEBase::coarsened_dof_values
             (*(system.solution), dof_map, parent, Uparent, false);
 
           // Loop over the neighbors of the parent
           for (auto n_p : parent->side_index_range())
             {
               if (parent->neighbor_ptr(n_p) != nullptr) // parent has a neighbor here
                 {
                   // Find the active neighbors in this direction
                   std::vector<const Elem *> active_neighbors;
                   parent->neighbor_ptr(n_p)->
                     active_family_tree_by_neighbor(active_neighbors,
                                                    parent);
                   // Compute the flux to each active neighbor
                   for (std::size_t a=0,
                         n_active_neighbors = active_neighbors.size();
                        a != n_active_neighbors; ++a)
                     {
                       const Elem * f = active_neighbors[a];
                       // FIXME - what about when f->level <
                       // parent->level()??
                       if (f->level() >= parent->level())
                         {
                           fine_context->pre_fe_reinit(system, f);
                           coarse_context->pre_fe_reinit(system, parent);
                           libmesh_assert_equal_to
                             (coarse_context->get_elem_solution().size(),
                              Uparent.size());
                           coarse_context->get_elem_solution() = Uparent;
 
                           this->reinit_sides();
 
                           // Loop over all significant variables in the system
                           for (var=0; var<n_vars; var++)
                             if (error_norm.weight(var) != 0.0 &&
                                 system.variable_type(var).family != SCALAR)
                               {
                                 this->internal_side_integration();
 
                                 error_per_cell[fine_context->get_elem().id()] +=
                                   static_cast<ErrorVectorReal>(fine_error);
                                 error_per_cell[coarse_context->get_elem().id()] +=
                                   static_cast<ErrorVectorReal>(coarse_error);
                               }
 
                           // Keep track of the number of internal flux
                           // sides found on each element
                           if (scale_by_n_flux_faces)
                             {
                               n_flux_faces[fine_context->get_elem().id()]++;
                               n_flux_faces[coarse_context->get_elem().id()] +=
                                 this->coarse_n_flux_faces_increment();
                             }
                         }
                     }
                 }
               else if (integrate_boundary_sides)
                 {
                   fine_context->pre_fe_reinit(system, parent);
                   libmesh_assert_equal_to
                     (fine_context->get_elem_solution().size(),
                      Uparent.size());
                   fine_context->get_elem_solution() = Uparent;
                   fine_context->side = cast_int<unsigned char>(n_p);
                   fine_context->side_fe_reinit();
 
                   // If we find a boundary flux for any variable,
                   // let's just count it as a flux face for all
                   // variables.  Otherwise we'd need to keep track of
                   // a separate n_flux_faces and error_per_cell for
                   // every single var.
                   bool found_boundary_flux = false;
 
                   for (var=0; var<n_vars; var++)
                     if (error_norm.weight(var) != 0.0 &&
                         system.variable_type(var).family != SCALAR)
                       {
                         if (this->boundary_side_integration())
                           {
                             error_per_cell[fine_context->get_elem().id()] +=
                               static_cast<ErrorVectorReal>(fine_error);
                             found_boundary_flux = true;
                           }
                       }
 
                   if (scale_by_n_flux_faces && found_boundary_flux)
                     n_flux_faces[fine_context->get_elem().id()]++;
                 }
             }
         }
 #endif // #ifdef LIBMESH_ENABLE_AMR
 
       // If we do any more flux integration, e will be the fine element
       fine_context->pre_fe_reinit(system, e);
 
       // Loop over the neighbors of element e
       for (auto n_e : e->side_index_range())
         {
           if ((e->neighbor_ptr(n_e) != nullptr) ||
               integrate_boundary_sides)
             {
               fine_context->side = cast_int<unsigned char>(n_e);
               fine_context->side_fe_reinit();
             }
 
           if (e->neighbor_ptr(n_e) != nullptr) // e is not on the boundary
             {
               const Elem * f           = e->neighbor_ptr(n_e);
               const dof_id_type f_id = f->id();
 
               // Compute flux jumps if we are in case 1 or case 2.
               if ((f->active() && (f->level() == e->level()) && (e_id < f_id))
                   || (f->level() < e->level()))
                 {
                   // f is now the coarse element
                   coarse_context->pre_fe_reinit(system, f);
 
                   this->reinit_sides();
 
                   // Loop over all significant variables in the system
                   for (var=0; var<n_vars; var++)
                     if (error_norm.weight(var) != 0.0 &&
                         system.variable_type(var).family != SCALAR)
                       {
                         this->internal_side_integration();
 
                         error_per_cell[fine_context->get_elem().id()] +=
                           static_cast<ErrorVectorReal>(fine_error);
                         error_per_cell[coarse_context->get_elem().id()] +=
                           static_cast<ErrorVectorReal>(coarse_error);
                       }
 
                   // Keep track of the number of internal flux
                   // sides found on each element
                   if (scale_by_n_flux_faces)
                     {
                       n_flux_faces[fine_context->get_elem().id()]++;
                       n_flux_faces[coarse_context->get_elem().id()] +=
                         this->coarse_n_flux_faces_increment();
                     }
                 } // end if (case1 || case2)
             } // if (e->neighbor(n_e) != nullptr)
 
           // Otherwise, e is on the boundary.  If it happens to
           // be on a Dirichlet boundary, we need not do anything.
           // On the other hand, if e is on a Neumann (flux) boundary
           // with grad(u).n = g, we need to compute the additional residual
           // (h * \int |g - grad(u_h).n|^2 dS)^(1/2).
           // We can only do this with some knowledge of the boundary
           // conditions, i.e. the user must have attached an appropriate
           // BC function.
           else if (integrate_boundary_sides)
             {
               bool found_boundary_flux = false;
 
               for (var=0; var<n_vars; var++)
                 if (error_norm.weight(var) != 0.0 &&
                     system.variable_type(var).family != SCALAR)
                   if (this->boundary_side_integration())
                     {
                       error_per_cell[fine_context->get_elem().id()] +=
                         static_cast<ErrorVectorReal>(fine_error);
                       found_boundary_flux = true;
                     }
 
               if (scale_by_n_flux_faces && found_boundary_flux)
                 n_flux_faces[fine_context->get_elem().id()]++;
             } // end if (e->neighbor_ptr(n_e) == nullptr)
         } // end loop over neighbors
     } // End loop over active local elements
 
 
   // Each processor has now computed the error contributions
   // for its local elements.  We need to sum the vector
   // and then take the square-root of each component.  Note
   // that we only need to sum if we are running on multiple
   // processors, and we only need to take the square-root
   // if the value is nonzero.  There will in general be many
   // zeros for the inactive elements.
 
   // First sum the vector of estimated error values
   this->reduce_error(error_per_cell, system.comm());
 
   // Compute the square-root of each component.
   for (auto i : index_range(error_per_cell))
     if (error_per_cell[i] != 0.)
       error_per_cell[i] = std::sqrt(error_per_cell[i]);
 
 
   if (this->scale_by_n_flux_faces)
     {
       // Sum the vector of flux face counts
       this->reduce_error(n_flux_faces, system.comm());
 
       // Sanity check: Make sure the number of flux faces is
       // always an integer value
 #ifdef DEBUG
       for (const auto & val : n_flux_faces)
         libmesh_assert_equal_to (val, static_cast<float>(static_cast<unsigned int>(val)));
 #endif
 
       // Scale the error by the number of flux faces for each element
       for (auto i : index_range(n_flux_faces))
         {
           if (n_flux_faces[i] == 0.0) // inactive or non-local element
             continue;
 
           error_per_cell[i] /= static_cast<ErrorVectorReal>(n_flux_faces[i]);
         }
     }
 
   // If we used a non-standard solution before, now is the time to fix
   // the current_local_solution
   if (solution_vector && solution_vector != system.solution.get())
     {
       NumericVector<Number> * newsol =
         const_cast<NumericVector<Number> *>(solution_vector);
       System & sys = const_cast<System &>(system);
       newsol->swap(*sys.solution);
       sys.update();
     }
 }

◆ estimate_errors() [1/2]

void libMesh::ErrorEstimator::estimate_errors	(	const EquationSystems &	equation_systems,
		ErrorVector &	error_per_cell,
		const std::map< const System *, SystemNorm > &	error_norms,
		const std::map< const System , const NumericVector< Number > > *	solution_vectors = `nullptr`,
		bool	estimate_parent_error = `false`
	)

virtualinherited

This virtual function can be redefined in derived classes, but by default computes the sum of the error_per_cell for each system in the equation_systems.

Currently this function ignores the error_norm member variable, and uses the function argument error_norms instead.

This function is named estimate_errors instead of estimate_error because otherwise C++ can get confused.

Reimplemented in libMesh::UniformRefinementEstimator.

Definition at line 47 of file error_estimator.C.

References libMesh::ErrorEstimator::error_norm, libMesh::ErrorEstimator::estimate_error(), libMesh::EquationSystems::get_system(), libMesh::index_range(), and libMesh::EquationSystems::n_systems().

 {
   SystemNorm old_error_norm = this->error_norm;
 
   // Sum the error values from each system
   for (unsigned int s = 0; s != equation_systems.n_systems(); ++s)
     {
       ErrorVector system_error_per_cell;
       const System & sys = equation_systems.get_system(s);
       if (error_norms.find(&sys) == error_norms.end())
         this->error_norm = old_error_norm;
       else
         this->error_norm = error_norms.find(&sys)->second;
 
       const NumericVector<Number> * solution_vector = nullptr;
       if (solution_vectors &&
           solution_vectors->find(&sys) != solution_vectors->end())
         solution_vector = solution_vectors->find(&sys)->second;
 
       this->estimate_error(sys, system_error_per_cell,
                            solution_vector, estimate_parent_error);
 
       if (s)
         {
           libmesh_assert_equal_to (error_per_cell.size(), system_error_per_cell.size());
           for (auto i : index_range(error_per_cell))
             error_per_cell[i] += system_error_per_cell[i];
         }
       else
         error_per_cell = system_error_per_cell;
     }
 
   // Restore our old state before returning
   this->error_norm = old_error_norm;
 }

◆ estimate_errors() [2/2]

void libMesh::ErrorEstimator::estimate_errors	(	const EquationSystems &	equation_systems,
		ErrorMap &	errors_per_cell,
		const std::map< const System , const NumericVector< Number > > *	solution_vectors = `nullptr`,
		bool	estimate_parent_error = `false`
	)

virtualinherited

This virtual function can be redefined in derived classes, but by default it calls estimate_error repeatedly to calculate the requested error vectors.

Currently this function ignores the error_norm.weight() values because it calculates each variable's error individually, unscaled.

The user selects which errors get computed by filling a map with error vectors: If errors_per_cell[&system][v] exists, it will be filled with the error values in variable v of system

FIXME: This is a default implementation - derived classes should reimplement it for efficiency.

Reimplemented in libMesh::UniformRefinementEstimator.

Definition at line 93 of file error_estimator.C.

References libMesh::ErrorEstimator::error_norm, libMesh::ErrorEstimator::estimate_error(), libMesh::EquationSystems::get_system(), libMesh::EquationSystems::n_systems(), n_vars, libMesh::System::n_vars(), and libMesh::SystemNorm::type().

 {
   SystemNorm old_error_norm = this->error_norm;
 
   // Find the requested error values from each system
   for (unsigned int s = 0; s != equation_systems.n_systems(); ++s)
     {
       const System & sys = equation_systems.get_system(s);
 
       unsigned int n_vars = sys.n_vars();
 
       for (unsigned int v = 0; v != n_vars; ++v)
         {
           // Only fill in ErrorVectors the user asks for
           if (errors_per_cell.find(std::make_pair(&sys, v)) ==
               errors_per_cell.end())
             continue;
 
           // Calculate error in only one variable
           std::vector<Real> weights(n_vars, 0.0);
           weights[v] = 1.0;
           this->error_norm =
             SystemNorm(std::vector<FEMNormType>(n_vars, old_error_norm.type(v)),
                        weights);
 
           const NumericVector<Number> * solution_vector = nullptr;
           if (solution_vectors &&
               solution_vectors->find(&sys) != solution_vectors->end())
             solution_vector = solution_vectors->find(&sys)->second;
 
           this->estimate_error
             (sys, *errors_per_cell[std::make_pair(&sys, v)],
              solution_vector, estimate_parent_error);
         }
     }
 
   // Restore our old state before returning
   this->error_norm = old_error_norm;
 }

◆ init_context()

void libMesh::KellyErrorEstimator::init_context ( FEMContext & c )

overrideprotectedvirtual

An initialization function, for requesting specific data from the FE objects.

Reimplemented from libMesh::JumpErrorEstimator.

Definition at line 59 of file kelly_error_estimator.C.

References libMesh::JumpErrorEstimator::coarse_context, libMesh::FEMContext::elem_dimensions(), libMesh::ErrorEstimator::error_norm, libMesh::JumpErrorEstimator::fine_context, libMesh::FEGenericBase< OutputType >::get_dphi(), libMesh::FEAbstract::get_normals(), n_vars, libMesh::DiffContext::n_vars(), and libMesh::SystemNorm::weight().

 {
   const unsigned int n_vars = c.n_vars();
   for (unsigned int v=0; v<n_vars; v++)
     {
       // Possibly skip this variable
       if (error_norm.weight(v) == 0.0) continue;
 
       // FIXME: Need to generalize this to vector-valued elements. [PB]
       FEBase * side_fe = nullptr;
 
       const std::set<unsigned char> & elem_dims =
         c.elem_dimensions();
 
       for (const auto & dim : elem_dims)
         {
           fine_context->get_side_fe( v, side_fe, dim );
 
           // We'll need gradients on both sides for flux jump computation
           side_fe->get_dphi();
 
           // But we only need normal vectors from one side
           if (&c != coarse_context.get())
             side_fe->get_normals();
         }
     }
 }

◆ internal_side_integration()

void libMesh::KellyErrorEstimator::internal_side_integration ( )

overrideprotectedvirtual

The function which calculates a normal derivative jump based error term on an internal side.

Implements libMesh::JumpErrorEstimator.

Definition at line 90 of file kelly_error_estimator.C.

References libMesh::JumpErrorEstimator::coarse_context, libMesh::JumpErrorEstimator::coarse_error, libMesh::Elem::dim(), libMesh::ErrorEstimator::error_norm, libMesh::JumpErrorEstimator::fine_context, libMesh::JumpErrorEstimator::fine_error, libMesh::FEGenericBase< OutputType >::get_dphi(), libMesh::FEAbstract::get_JxW(), libMesh::FEAbstract::get_normals(), libMesh::Elem::hmax(), libMesh::FEAbstract::n_quadrature_points(), libMesh::TensorTools::norm_sq(), libMesh::Real, libMesh::JumpErrorEstimator::var, and libMesh::SystemNorm::weight().

 {
   const Elem & coarse_elem = coarse_context->get_elem();
   const Elem & fine_elem = fine_context->get_elem();
 
   FEBase * fe_fine = nullptr;
   fine_context->get_side_fe( var, fe_fine, fine_elem.dim() );
 
   FEBase * fe_coarse = nullptr;
   coarse_context->get_side_fe( var, fe_coarse, fine_elem.dim() );
 
   Real error = 1.e-30;
   unsigned int n_qp = fe_fine->n_quadrature_points();
 
   std::vector<std::vector<RealGradient>> dphi_coarse = fe_coarse->get_dphi();
   std::vector<std::vector<RealGradient>> dphi_fine = fe_fine->get_dphi();
   std::vector<Point> face_normals = fe_fine->get_normals();
   std::vector<Real> JxW_face = fe_fine->get_JxW();
 
   for (unsigned int qp=0; qp != n_qp; ++qp)
     {
       // Calculate solution gradients on fine and coarse elements
       // at this quadrature point
       Gradient
         grad_fine   = fine_context->side_gradient(var, qp),
         grad_coarse = coarse_context->side_gradient(var, qp);
 
       // Find the jump in the normal derivative
       // at this quadrature point
       const Number jump = (grad_fine - grad_coarse)*face_normals[qp];
       const Real jump2 = TensorTools::norm_sq(jump);
 
       // Accumulate the jump integral
       error += JxW_face[qp] * jump2;
     }
 
   // Add the h-weighted jump integral to each error term
   fine_error =
     error * fine_elem.hmax() * error_norm.weight(var);
   coarse_error =
     error * coarse_elem.hmax() * error_norm.weight(var);
 }

◆ operator=() [1/2]

KellyErrorEstimator& libMesh::KellyErrorEstimator::operator= ( const KellyErrorEstimator & )

delete

◆ operator=() [2/2]

KellyErrorEstimator& libMesh::KellyErrorEstimator::operator= ( KellyErrorEstimator && )

default

◆ reduce_error()

void libMesh::ErrorEstimator::reduce_error	(	std::vector< ErrorVectorReal > &	error_per_cell,
		const Parallel::Communicator &	comm
	)		const

protectedinherited

This method takes the local error contributions in error_per_cell from each processor and combines them to get the global error vector.

Definition at line 32 of file error_estimator.C.

References libMesh::Parallel::Communicator::sum().

Referenced by libMesh::UniformRefinementEstimator::_estimate_error(), libMesh::WeightedPatchRecoveryErrorEstimator::estimate_error(), libMesh::PatchRecoveryErrorEstimator::estimate_error(), libMesh::JumpErrorEstimator::estimate_error(), libMesh::AdjointRefinementEstimator::estimate_error(), and libMesh::ExactErrorEstimator::estimate_error().

 {
   // This function must be run on all processors at once
   // parallel_object_only();
 
   // Each processor has now computed the error contributions
   // for its local elements.  We may need to sum the vector to
   // recover the error for each element.
 
   comm.sum(error_per_cell);
 }

◆ reinit_sides()

void libMesh::JumpErrorEstimator::reinit_sides ( )

protectedinherited

A utility function to reinit the finite element data on elements sharing a side

Definition at line 424 of file jump_error_estimator.C.

References libMesh::JumpErrorEstimator::coarse_context, libMesh::ErrorEstimator::error_norm, libMesh::JumpErrorEstimator::fine_context, libMesh::FEAbstract::get_fe_type(), libMesh::FEInterface::inverse_map(), n_vars, libMesh::FEAbstract::reinit(), libMesh::SCALAR, and libMesh::SystemNorm::weight().

Referenced by libMesh::JumpErrorEstimator::estimate_error().

 {
   fine_context->side_fe_reinit();
 
   unsigned short dim = fine_context->get_elem().dim();
   libmesh_assert_equal_to(dim, coarse_context->get_elem().dim());
 
   FEBase * fe_fine = nullptr;
   fine_context->get_side_fe( 0, fe_fine, dim );
 
   // Get the physical locations of the fine element quadrature points
   std::vector<Point> qface_point = fe_fine->get_xyz();
 
   // Find the master quadrature point locations on the coarse element
   FEBase * fe_coarse = nullptr;
   coarse_context->get_side_fe( 0, fe_coarse, dim );
 
   std::vector<Point> qp_coarse;
 
   FEInterface::inverse_map
     (coarse_context->get_elem().dim(), fe_coarse->get_fe_type(),
      &coarse_context->get_elem(), qface_point, qp_coarse);
 
   // The number of variables in the system
   const unsigned int n_vars = fine_context->n_vars();
 
   // Calculate all coarse element shape functions at those locations
   for (unsigned int v=0; v<n_vars; v++)
     if (error_norm.weight(v) != 0.0 &&
         fine_context->get_system().variable_type(v).family != SCALAR)
       {
         coarse_context->get_side_fe( v, fe_coarse, dim );
         fe_coarse->reinit (&coarse_context->get_elem(), &qp_coarse);
       }
 }

◆ type()

ErrorEstimatorType libMesh::KellyErrorEstimator::type ( ) const

overridevirtual

Returns: The type for the ErrorEstimator subclass.

Implements libMesh::ErrorEstimator.

Definition at line 51 of file kelly_error_estimator.C.

References libMesh::KELLY.

 {
   return KELLY;
 }

Member Data Documentation

◆ _bc_function

std::pair<bool,Real>(* libMesh::KellyErrorEstimator::_bc_function) (const System &system, const Point &p, const std::string &var_name)

protected

Pointer to function that provides BC information.

Definition at line 119 of file kelly_error_estimator.h.

Referenced by attach_flux_bc_function(), and boundary_side_integration().

◆ coarse_context

std::unique_ptr<FEMContext> libMesh::JumpErrorEstimator::coarse_context

protectedinherited

Definition at line 144 of file jump_error_estimator.h.

Referenced by libMesh::JumpErrorEstimator::coarse_n_flux_faces_increment(), libMesh::JumpErrorEstimator::estimate_error(), init_context(), libMesh::LaplacianErrorEstimator::internal_side_integration(), libMesh::DiscontinuityMeasure::internal_side_integration(), internal_side_integration(), and libMesh::JumpErrorEstimator::reinit_sides().

◆ coarse_error

Real libMesh::JumpErrorEstimator::coarse_error

protectedinherited

Definition at line 149 of file jump_error_estimator.h.

Referenced by libMesh::JumpErrorEstimator::estimate_error(), libMesh::LaplacianErrorEstimator::internal_side_integration(), libMesh::DiscontinuityMeasure::internal_side_integration(), and internal_side_integration().

◆ error_norm

SystemNorm libMesh::ErrorEstimator::error_norm

inherited

When estimating the error in a single system, the error_norm is used to control the scaling and norm choice for each variable. Not all estimators will support all norm choices. The default scaling is for all variables to be weighted equally. The default norm choice depends on the error estimator.

Part of this functionality was supported via component_scale and sobolev_order in older libMesh versions, and a small part was supported via component_mask in even older versions. Hopefully the encapsulation here will allow us to avoid changing this API again.

Definition at line 161 of file error_estimator.h.

Referenced by libMesh::UniformRefinementEstimator::_estimate_error(), libMesh::AdjointRefinementEstimator::AdjointRefinementEstimator(), libMesh::DiscontinuityMeasure::boundary_side_integration(), boundary_side_integration(), libMesh::DiscontinuityMeasure::DiscontinuityMeasure(), libMesh::JumpErrorEstimator::estimate_error(), libMesh::AdjointResidualErrorEstimator::estimate_error(), libMesh::ExactErrorEstimator::estimate_error(), libMesh::ErrorEstimator::estimate_errors(), libMesh::ExactErrorEstimator::ExactErrorEstimator(), libMesh::ExactErrorEstimator::find_squared_element_error(), libMesh::LaplacianErrorEstimator::init_context(), libMesh::DiscontinuityMeasure::init_context(), init_context(), libMesh::LaplacianErrorEstimator::internal_side_integration(), libMesh::DiscontinuityMeasure::internal_side_integration(), internal_side_integration(), KellyErrorEstimator(), libMesh::LaplacianErrorEstimator::LaplacianErrorEstimator(), libMesh::WeightedPatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::PatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::PatchRecoveryErrorEstimator::PatchRecoveryErrorEstimator(), libMesh::JumpErrorEstimator::reinit_sides(), and libMesh::UniformRefinementEstimator::UniformRefinementEstimator().

◆ fine_context

std::unique_ptr<FEMContext> libMesh::JumpErrorEstimator::fine_context

protectedinherited

Context objects for integrating on the fine and coarse elements sharing a face

Definition at line 144 of file jump_error_estimator.h.

Referenced by libMesh::DiscontinuityMeasure::boundary_side_integration(), boundary_side_integration(), libMesh::JumpErrorEstimator::coarse_n_flux_faces_increment(), libMesh::JumpErrorEstimator::estimate_error(), libMesh::LaplacianErrorEstimator::init_context(), libMesh::DiscontinuityMeasure::init_context(), init_context(), libMesh::LaplacianErrorEstimator::internal_side_integration(), libMesh::DiscontinuityMeasure::internal_side_integration(), internal_side_integration(), and libMesh::JumpErrorEstimator::reinit_sides().

◆ fine_error

Real libMesh::JumpErrorEstimator::fine_error

protectedinherited

The fine and coarse error values to be set by each side_integration();

Definition at line 149 of file jump_error_estimator.h.

Referenced by libMesh::DiscontinuityMeasure::boundary_side_integration(), boundary_side_integration(), libMesh::JumpErrorEstimator::estimate_error(), libMesh::LaplacianErrorEstimator::internal_side_integration(), libMesh::DiscontinuityMeasure::internal_side_integration(), and internal_side_integration().

◆ integrate_boundary_sides

bool libMesh::JumpErrorEstimator::integrate_boundary_sides

protectedinherited

A boolean flag, by default false, to be set to true if integrations with boundary_side_integration() should be performed

Definition at line 138 of file jump_error_estimator.h.

Referenced by libMesh::DiscontinuityMeasure::attach_essential_bc_function(), attach_flux_bc_function(), and libMesh::JumpErrorEstimator::estimate_error().

◆ scale_by_n_flux_faces

bool libMesh::JumpErrorEstimator::scale_by_n_flux_faces

inherited

This boolean flag allows you to scale the error indicator result for each element by the number of "flux faces" the element actually has. This tends to weight more evenly cells which are on the boundaries and thus have fewer contributions to their flux. The value is initialized to false, simply set it to true if you want to use the feature.

Definition at line 99 of file jump_error_estimator.h.

Referenced by libMesh::JumpErrorEstimator::estimate_error().

◆ var

unsigned int libMesh::JumpErrorEstimator::var

protectedinherited

The variable number currently being evaluated

Definition at line 154 of file jump_error_estimator.h.

Referenced by libMesh::DiscontinuityMeasure::boundary_side_integration(), boundary_side_integration(), libMesh::JumpErrorEstimator::estimate_error(), libMesh::LaplacianErrorEstimator::internal_side_integration(), libMesh::DiscontinuityMeasure::internal_side_integration(), and internal_side_integration().

The documentation for this class was generated from the following files:

Public Types

Public Member Functions

Public Attributes

Protected Member Functions

Protected Attributes

Detailed Description

Member Typedef Documentation

◆ ErrorMap

Constructor & Destructor Documentation

◆ KellyErrorEstimator() [1/3]

◆ KellyErrorEstimator() [2/3]

◆ KellyErrorEstimator() [3/3]

◆ ~KellyErrorEstimator()

Member Function Documentation

◆ attach_flux_bc_function()

◆ boundary_side_integration()

◆ coarse_n_flux_faces_increment()

◆ estimate_error()

◆ estimate_errors() [1/2]

◆ estimate_errors() [2/2]

◆ init_context()

◆ internal_side_integration()

◆ operator=() [1/2]

◆ operator=() [2/2]

◆ reduce_error()

◆ reinit_sides()

◆ type()

Member Data Documentation

◆ _bc_function

◆ coarse_context

◆ coarse_error

◆ error_norm

◆ fine_context

◆ fine_error

◆ integrate_boundary_sides

◆ scale_by_n_flux_faces

◆ var