libMesh::JumpErrorEstimator Class Reference

This abstract base class implements utility functions for error estimators which are based on integrated jumps between elements. More...

#include <jump_error_estimator.h>

Inheritance diagram for libMesh::JumpErrorEstimator:

Public Types
typedef std::map< std::pair< const System *, unsigned int >, std::unique_ptr< ErrorVector > >	ErrorMap
	When calculating many error vectors at once, we need a data structure to hold them all. More...

Public Member Functions
	JumpErrorEstimator ()
	Constructor. More...
	JumpErrorEstimator (const JumpErrorEstimator &)=delete
	This class cannot be (default) copy constructed/assigned because it has unique_ptr members. More...
JumpErrorEstimator &	operator= (const JumpErrorEstimator &)=delete
	JumpErrorEstimator (JumpErrorEstimator &&)=default
	Defaulted move ctor, move assignment operator, and destructor. More...
JumpErrorEstimator &	operator= (JumpErrorEstimator &&)=default
virtual	~JumpErrorEstimator ()=default
virtual void	estimate_error (const System &system, ErrorVector &error_per_cell, const NumericVector< Number > *solution_vector=nullptr, bool estimate_parent_error=false) override
	This function uses the derived class's jump error estimate formula to estimate the error on each cell. More...
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)
	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. More...
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)
	This virtual function can be redefined in derived classes, but by default it calls estimate_error repeatedly to calculate the requested error vectors. More...
virtual ErrorEstimatorType	type () const =0

Public Attributes
bool	scale_by_n_flux_faces
	This boolean flag allows you to scale the error indicator result for each element by the number of "flux faces" the element actually has. More...
bool	use_unweighted_quadrature_rules
	This boolean flag allows you to use "unweighted" quadrature rules (sized to exactly integrate unweighted shape functions in master element space) rather than "default" quadrature rules (sized to exactly integrate polynomials of one higher degree than mass matrix terms). More...
bool	integrate_slits
	A boolean flag, by default false, to be set to true if integrations should be performed on "slits" where two elements' faces overlap even if those elements are not connected by neighbor links. More...
SystemNorm	error_norm
	When estimating the error in a single system, the error_norm is used to control the scaling and norm choice for each variable. More...

Protected Member Functions
void	reinit_sides ()
	A utility function to reinit the finite element data on elements sharing a side. More...
float	coarse_n_flux_faces_increment ()
	A utility function to correctly increase n_flux_faces for the coarse element. More...
virtual void	init_context (FEMContext &c)
	An initialization function, to give derived classes a chance to request specific data from the FE objects. More...
virtual void	internal_side_integration ()=0
	The function, to be implemented by derived classes, which calculates an error term on an internal side. More...
virtual bool	boundary_side_integration ()
	The function, to be implemented by derived classes, which calculates an error term on a boundary side. More...
void	reduce_error (std::vector< ErrorVectorReal > &error_per_cell, const Parallel::Communicator &comm) const
	This method takes the local error contributions in error_per_cell from each processor and combines them to get the global error vector. More...

Protected Attributes
bool	integrate_boundary_sides
	A boolean flag, by default false, to be set to true if integrations with boundary_side_integration() should be performed. More...
std::unique_ptr< FEMContext >	fine_context
	Context objects for integrating on the fine and coarse elements sharing a face. More...
std::unique_ptr< FEMContext >	coarse_context
Real	fine_error
	The fine and coarse error values to be set by each side_integration();. More...
Real	coarse_error
unsigned int	var
	The variable number currently being evaluated. More...

Detailed Description

This abstract base class implements utility functions for error estimators which are based on integrated jumps between elements.

Author

Roy H. Stogner

Date

2006

Definition at line 48 of file jump_error_estimator.h.

Member Typedef Documentation

ErrorMap

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

inherited

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

Definition at line 121 of file error_estimator.h.

Constructor & Destructor Documentation

JumpErrorEstimator() [1/3]

libMesh::JumpErrorEstimator::JumpErrorEstimator ( )

inline

Constructor.

Definition at line 55 of file jump_error_estimator.h.

    : ErrorEstimator(),
      scale_by_n_flux_faces(false),
      use_unweighted_quadrature_rules(false),
      integrate_slits(false),
      integrate_boundary_sides(false),
      fine_context(),
      coarse_context(),
      fine_error(0),
      coarse_error(0) {}

JumpErrorEstimator() [2/3]

libMesh::JumpErrorEstimator::JumpErrorEstimator ( const JumpErrorEstimator &  )

delete

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

Explicitly deleting these functions is the best way to document this fact.

JumpErrorEstimator() [3/3]

libMesh::JumpErrorEstimator::JumpErrorEstimator ( JumpErrorEstimator &&  )

default

Defaulted move ctor, move assignment operator, and destructor.

~JumpErrorEstimator()

virtual libMesh::JumpErrorEstimator::~JumpErrorEstimator ( )

virtualdefault

Member Function Documentation

boundary_side_integration()

virtual bool libMesh::JumpErrorEstimator::boundary_side_integration ( )

inlineprotectedvirtual

The function, to be implemented by derived classes, which calculates an error term on a boundary side.

Returns

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

Reimplemented in libMesh::KellyErrorEstimator, and libMesh::DiscontinuityMeasure.

Definition at line 163 of file jump_error_estimator.h.

Referenced by estimate_error().

{ return false; }

coarse_n_flux_faces_increment()

float libMesh::JumpErrorEstimator::coarse_n_flux_faces_increment ( )

protected

A utility function to correctly increase n_flux_faces for the coarse element.

Definition at line 591 of file jump_error_estimator.C.

References coarse_context, dim, and fine_context.

Referenced by 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  )

overridevirtual

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 56 of file jump_error_estimator.C.

References libMesh::Elem::active(), libMesh::ElemInternal::active_family_tree_by_neighbor(), boundary_side_integration(), libMesh::FEAbstract::build(), libMesh::Elem::child_ref_range(), coarse_context, coarse_error, coarse_n_flux_faces_increment(), libMesh::FEGenericBase< OutputType >::coarsened_dof_values(), dim, libMesh::ErrorEstimator::error_norm, libMesh::ErrorVectorReal, fine_context, fine_error, libMesh::NumericVector< T >::get(), libMesh::System::get_dof_map(), libMesh::System::get_mesh(), libMesh::FEAbstract::get_xyz(), libMesh::DofObject::id(), libMesh::index_range(), init_context(), integrate_boundary_sides, integrate_slits, internal_side_integration(), libMesh::FEMap::inverse_map(), libMesh::Elem::level(), libMesh::libmesh_assert(), libMesh::libmesh_ignore(), libMesh::make_range(), mesh, n_vars, libMesh::System::n_vars(), libMesh::Elem::neighbor_ptr(), libMesh::Elem::parent(), libMesh::Real, libMesh::ErrorEstimator::reduce_error(), libMesh::FEAbstract::reinit(), reinit_sides(), libMesh::SCALAR, scale_by_n_flux_faces, libMesh::Elem::side_index_range(), libMesh::DenseVector< T >::size(), libMesh::System::solution, libMesh::NumericVector< T >::swap(), libMesh::FEType::unweighted_quadrature_rule(), use_unweighted_quadrature_rules, var, and libMesh::SystemNorm::weight().

Referenced by assemble_and_solve(), main(), and SlitMeshRefinedSystemTest::testSystem().

{
  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();
    }

  // We don't use full elem_jacobian or subjacobians here.
  fine_context = std::make_unique<FEMContext>
    (system, nullptr, /* allocate_local_matrices = */ false);
  coarse_context = std::make_unique<FEMContext>
    (system, nullptr, /* allocate_local_matrices = */ false);

  // Don't overintegrate - we're evaluating differences of FE values,
  // not products of them.
  if (this->use_unweighted_quadrature_rules)
    fine_context->use_unweighted_quadrature_rules(system.extra_quadrature_order);

  // 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);

  // If we're integrating jumps across mesh slits, then we'll need a
  // point locator to find slits, and we'll need to integrate point by
  // point on sides.
  std::unique_ptr<PointLocatorBase> point_locator;
  std::unique_ptr<const Elem> side_ptr;

  if (integrate_slits)
    point_locator = mesh.sub_point_locator();

  // 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

      if (e->infinite())
      {
         libmesh_warning("Warning: Jumps on the border of infinite elements are ignored."
                         << std::endl);
         continue;
      }

      // 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];

                      if (f ->infinite()) // don't take infinite elements into account
                         continue;

                      // 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();
            }

          // e is not on the boundary (infinite elements are treated as boundary)
          if (e->neighbor_ptr(n_e) != nullptr
              && !e->neighbor_ptr(n_e) ->infinite())
            {

              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)

          // e might not have a neighbor_ptr, but might still have
          // another element sharing its side.  This can happen in a
          // mesh where solution continuity is maintained via nodal
          // constraint rows.
          else if (integrate_slits)
            {
              side_ptr = e->build_side_ptr(n_e);
              std::set<const Elem *> candidate_elements;
              (*point_locator)(side_ptr->vertex_average(), candidate_elements);

              // We should have at least found ourselves...
              libmesh_assert(candidate_elements.count(e));

              // If we only found ourselves, this probably isn't a
              // slit; we don't yet support meshes so non-conforming
              // as to have overlap of part of an element side without
              // overlap of its center.
              if (candidate_elements.size() < 2)
                continue;

              FEType hardest_fe_type = fine_context->find_hardest_fe_type();

              auto dim = e->dim();

              auto side_qrule =
                hardest_fe_type.unweighted_quadrature_rule
                (dim-1, system.extra_quadrature_order);
              auto side_fe = FEAbstract::build(dim, hardest_fe_type);
              side_fe->attach_quadrature_rule(side_qrule.get());
              const std::vector<Point> & qface_point = side_fe->get_xyz();
              side_fe->reinit(e, n_e);

              for (auto qp : make_range(side_qrule->n_points()))
                {
                  const Point p = qface_point[qp];
                  const std::vector<Point> qp_pointvec(1, p);
                  std::set<const Elem *> side_elements;
                  (*point_locator)(side_ptr->vertex_average(), side_elements);

                  // If we have multiple neighbors meeting here we'll just
                  // take weighted jumps from all of them.
                  //
                  // We'll also do integrations from both sides of slits,
                  // rather than try to figure out a disambiguation rule
                  // that makes sense for non-conforming slits in general.
                  // This means we want an extra factor of 0.5 on the
                  // integrals to compensate for doubling them.
                  const std::size_t n_neighbors = side_elements.size() - 1;
                  const Real neighbor_frac = Real(1)/n_neighbors;

                  const std::vector<Real>
                    qp_weightvec(1, neighbor_frac * side_qrule->w(qp));

                  for (const Elem * f : side_elements)
                    {
                      if (f == e)
                        continue;

                      coarse_context->pre_fe_reinit(system, f);
                      fine_context->pre_fe_reinit(system, e);
                      std::vector<Point> qp_coarse, qp_fine;
                      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)
                          {
                            FEBase * coarse_fe = coarse_context->get_side_fe(v, dim);
                            if (qp_coarse.empty())
                              FEMap::inverse_map (dim, f, qp_pointvec, qp_coarse);
                            coarse_fe->reinit(f, &qp_coarse, &qp_weightvec);
                            FEBase * fine_fe = fine_context->get_side_fe(v, dim);
                            if (qp_fine.empty())
                              FEMap::inverse_map (dim, e, qp_pointvec, qp_fine);
                            fine_fe->reinit(e, &qp_fine, &qp_weightvec);
                          }

                      // 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);
                          }
                    }
                }
            }

          // 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)
            {
              if (integrate_slits)
                libmesh_not_implemented();

              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(), libMesh::make_range(), and libMesh::EquationSystems::n_systems().

{
  SystemNorm old_error_norm = this->error_norm;

  // Sum the error values from each system
  for (auto s : make_range(equation_systems.n_systems()))
    {
      ErrorVector system_error_per_cell;
      const System & sys = equation_systems.get_system(s);
      if (const auto it = error_norms.find(&sys);
          it == error_norms.end())
        this->error_norm = old_error_norm;
      else
        this->error_norm = it->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.

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

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

Reimplemented in libMesh::UniformRefinementEstimator.

Definition at line 94 of file error_estimator.C.

References libMesh::ErrorEstimator::error_norm, libMesh::ErrorEstimator::estimate_error(), libMesh::EquationSystems::get_system(), libMesh::make_range(), 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 (auto s : make_range(equation_systems.n_systems()))
    {
      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.count(std::make_pair(&sys, v)))
            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::JumpErrorEstimator::init_context ( FEMContext &  c )

protectedvirtual

An initialization function, to give derived classes a chance to request specific data from the FE objects.

Reimplemented in libMesh::KellyErrorEstimator, libMesh::DiscontinuityMeasure, and libMesh::LaplacianErrorEstimator.

Definition at line 48 of file jump_error_estimator.C.

Referenced by estimate_error().

{
  // Derived classes are *supposed* to rederive this
  libmesh_deprecated();
}

internal_side_integration()

virtual void libMesh::JumpErrorEstimator::internal_side_integration ( )

protectedpure virtual

The function, to be implemented by derived classes, which calculates an error term on an internal side.

Implemented in libMesh::KellyErrorEstimator, libMesh::DiscontinuityMeasure, and libMesh::LaplacianErrorEstimator.

Referenced by estimate_error().

operator=() [1/2]

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

delete

operator=() [2/2]

JumpErrorEstimator& libMesh::JumpErrorEstimator::operator= ( JumpErrorEstimator &&  )

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 TIMPI::Communicator::sum().

Referenced by libMesh::UniformRefinementEstimator::_estimate_error(), libMesh::PatchRecoveryErrorEstimator::estimate_error(), libMesh::WeightedPatchRecoveryErrorEstimator::estimate_error(), 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 ( )

protected

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

Definition at line 553 of file jump_error_estimator.C.

References coarse_context, dim, libMesh::ErrorEstimator::error_norm, fine_context, libMesh::FEMap::inverse_map(), n_vars, libMesh::FEAbstract::reinit(), libMesh::SCALAR, and libMesh::SystemNorm::weight().

Referenced by 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;

  FEMap::inverse_map (coarse_context->get_elem().dim(),
                      &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()

virtual ErrorEstimatorType libMesh::ErrorEstimator::type ( ) const

pure virtualinherited

Returns

The type for the ErrorEstimator subclass.

Implemented in libMesh::ExactErrorEstimator, libMesh::AdjointResidualErrorEstimator, libMesh::UniformRefinementEstimator, libMesh::AdjointRefinementEstimator, libMesh::PatchRecoveryErrorEstimator, libMesh::KellyErrorEstimator, libMesh::WeightedPatchRecoveryErrorEstimator, libMesh::DiscontinuityMeasure, and libMesh::LaplacianErrorEstimator.

Member Data Documentation

coarse_context

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

protected

Definition at line 175 of file jump_error_estimator.h.

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

coarse_error

Real libMesh::JumpErrorEstimator::coarse_error

protected

Definition at line 180 of file jump_error_estimator.h.

Referenced by estimate_error(), libMesh::LaplacianErrorEstimator::internal_side_integration(), libMesh::DiscontinuityMeasure::internal_side_integration(), and libMesh::KellyErrorEstimator::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 158 of file error_estimator.h.

Referenced by libMesh::UniformRefinementEstimator::_estimate_error(), libMesh::AdjointRefinementEstimator::AdjointRefinementEstimator(), libMesh::DiscontinuityMeasure::boundary_side_integration(), libMesh::KellyErrorEstimator::boundary_side_integration(), libMesh::DiscontinuityMeasure::DiscontinuityMeasure(), 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(), libMesh::KellyErrorEstimator::init_context(), libMesh::LaplacianErrorEstimator::internal_side_integration(), libMesh::DiscontinuityMeasure::internal_side_integration(), libMesh::KellyErrorEstimator::internal_side_integration(), libMesh::KellyErrorEstimator::KellyErrorEstimator(), libMesh::LaplacianErrorEstimator::LaplacianErrorEstimator(), libMesh::WeightedPatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::PatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::PatchRecoveryErrorEstimator::PatchRecoveryErrorEstimator(), reinit_sides(), and libMesh::UniformRefinementEstimator::UniformRefinementEstimator().

fine_context

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

protected

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

Definition at line 175 of file jump_error_estimator.h.

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

fine_error

Real libMesh::JumpErrorEstimator::fine_error

protected

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

Definition at line 180 of file jump_error_estimator.h.

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

integrate_boundary_sides

bool libMesh::JumpErrorEstimator::integrate_boundary_sides

protected

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

Definition at line 169 of file jump_error_estimator.h.

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

integrate_slits

bool libMesh::JumpErrorEstimator::integrate_slits

A boolean flag, by default false, to be set to true if integrations should be performed on "slits" where two elements' faces overlap even if those elements are not connected by neighbor links.

This may only be useful for flex-IGA meshes, where highly-nonconforming meshes are given pseudo-conforming solutions via nodal constraints.

Note that, to safely use this option in parallel, it is also necessary to expand the default algebraic ghosting requirements to include elements on the opposite sides of slits from local elements.

Definition at line 130 of file jump_error_estimator.h.

Referenced by estimate_error(), and SlitMeshRefinedSystemTest::testSystem().

scale_by_n_flux_faces

bool libMesh::JumpErrorEstimator::scale_by_n_flux_faces

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 100 of file jump_error_estimator.h.

Referenced by estimate_error().

use_unweighted_quadrature_rules

bool libMesh::JumpErrorEstimator::use_unweighted_quadrature_rules

This boolean flag allows you to use "unweighted" quadrature rules (sized to exactly integrate unweighted shape functions in master element space) rather than "default" quadrature rules (sized to exactly integrate polynomials of one higher degree than mass matrix terms).

The results with the former, lower-order rules will be somewhat less accurate in many cases but will be much cheaper to compute.

The value is initialized to false, simply set it to true if you want to use the feature.

Definition at line 114 of file jump_error_estimator.h.

Referenced by estimate_error(), and main().

var

unsigned int libMesh::JumpErrorEstimator::var

protected

The variable number currently being evaluated.

Definition at line 185 of file jump_error_estimator.h.

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

---

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

• include/error_estimation/jump_error_estimator.h
• src/error_estimation/jump_error_estimator.C