libMesh
Public Types | Public Member Functions | Public Attributes | Protected Member Functions | Protected Attributes | List of all members
libMesh::JumpErrorEstimator Class Referenceabstract

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:
[legend]

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...
 
JumpErrorEstimatoroperator= (const JumpErrorEstimator &)=delete
 
 JumpErrorEstimator (JumpErrorEstimator &&)=default
 Defaulted move ctor, move assignment operator, and destructor. More...
 
JumpErrorEstimatoroperator= (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< FEMContextfine_context
 Context objects for integrating on the fine and coarse elements sharing a face. More...
 
std::unique_ptr< FEMContextcoarse_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.

56  : ErrorEstimator(),
57  scale_by_n_flux_faces(false),
59  integrate_slits(false),
61  fine_context(),
63  fine_error(0),
64  coarse_error(0) {}
std::unique_ptr< FEMContext > fine_context
Context objects for integrating on the fine and coarse elements sharing a face.
bool integrate_slits
A boolean flag, by default false, to be set to true if integrations should be performed on "slits" wh...
std::unique_ptr< FEMContext > coarse_context
ErrorEstimator()=default
Constructor.
Real fine_error
The fine and coarse error values to be set by each side_integration();.
bool use_unweighted_quadrature_rules
This boolean flag allows you to use "unweighted" quadrature rules (sized to exactly integrate unweigh...
bool scale_by_n_flux_faces
This boolean flag allows you to scale the error indicator result for each element by the number of "f...
bool integrate_boundary_sides
A boolean flag, by default false, to be set to true if integrations with boundary_side_integration() ...

◆ 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().

163 { 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().

592 {
593  // Keep track of the number of internal flux sides found on each
594  // element
595  unsigned short dim = coarse_context->get_elem().dim();
596 
597  const unsigned int divisor =
598  1 << (dim-1)*(fine_context->get_elem().level() -
599  coarse_context->get_elem().level());
600 
601  // With a difference of n levels between fine and coarse elements,
602  // we compute a fractional flux face for the coarse element by adding:
603  // 1/2^n in 2D
604  // 1/4^n in 3D
605  // each time. This code will get hit 2^n times in 2D and 4^n
606  // times in 3D so that the final flux face count for the coarse
607  // element will be an integer value.
608 
609  return 1.0f / static_cast<float>(divisor);
610 }
std::unique_ptr< FEMContext > fine_context
Context objects for integrating on the fine and coarse elements sharing a face.
unsigned int dim
std::unique_ptr< FEMContext > coarse_context

◆ 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().

60 {
61  LOG_SCOPE("estimate_error()", "JumpErrorEstimator");
62 
99  // This parameter is not used when !LIBMESH_ENABLE_AMR.
100  libmesh_ignore(estimate_parent_error);
101 
102  // The current mesh
103  const MeshBase & mesh = system.get_mesh();
104 
105  // The number of variables in the system
106  const unsigned int n_vars = system.n_vars();
107 
108  // The DofMap for this system
109 #ifdef LIBMESH_ENABLE_AMR
110  const DofMap & dof_map = system.get_dof_map();
111 #endif
112 
113  // Resize the error_per_cell vector to be
114  // the number of elements, initialize it to 0.
115  error_per_cell.resize (mesh.max_elem_id());
116  std::fill (error_per_cell.begin(), error_per_cell.end(), 0.);
117 
118  // Declare a vector of floats which is as long as
119  // error_per_cell above, and fill with zeros. This vector will be
120  // used to keep track of the number of edges (faces) on each active
121  // element which are either:
122  // 1) an internal edge
123  // 2) an edge on a Neumann boundary for which a boundary condition
124  // function has been specified.
125  // The error estimator can be scaled by the number of flux edges (faces)
126  // which the element actually has to obtain a more uniform measure
127  // of the error. Use floats instead of ints since in case 2 (above)
128  // f gets 1/2 of a flux face contribution from each of his
129  // neighbors
130  std::vector<float> n_flux_faces;
132  n_flux_faces.resize(error_per_cell.size(), 0);
133 
134  // Prepare current_local_solution to localize a non-standard
135  // solution vector if necessary
136  if (solution_vector && solution_vector != system.solution.get())
137  {
138  NumericVector<Number> * newsol =
139  const_cast<NumericVector<Number> *>(solution_vector);
140  System & sys = const_cast<System &>(system);
141  newsol->swap(*sys.solution);
142  sys.update();
143  }
144 
145  // We don't use full elem_jacobian or subjacobians here.
146  fine_context = std::make_unique<FEMContext>
147  (system, nullptr, /* allocate_local_matrices = */ false);
148  coarse_context = std::make_unique<FEMContext>
149  (system, nullptr, /* allocate_local_matrices = */ false);
150 
151  // Don't overintegrate - we're evaluating differences of FE values,
152  // not products of them.
154  fine_context->use_unweighted_quadrature_rules(system.extra_quadrature_order);
155 
156  // Loop over all the variables we've been requested to find jumps in, to
157  // pre-request
158  for (var=0; var<n_vars; var++)
159  {
160  // Skip variables which aren't part of our norm,
161  // as well as SCALAR variables, which have no jumps
162  if (error_norm.weight(var) == 0.0 ||
163  system.variable_type(var).family == SCALAR)
164  continue;
165 
166  // FIXME: Need to generalize this to vector-valued elements. [PB]
167  FEBase * side_fe = nullptr;
168 
169  const std::set<unsigned char> & elem_dims =
170  fine_context->elem_dimensions();
171 
172  for (const auto & dim : elem_dims)
173  {
174  fine_context->get_side_fe( var, side_fe, dim );
175 
176  side_fe->get_xyz();
177  }
178  }
179 
180  this->init_context(*fine_context);
181  this->init_context(*coarse_context);
182 
183  // If we're integrating jumps across mesh slits, then we'll need a
184  // point locator to find slits, and we'll need to integrate point by
185  // point on sides.
186  std::unique_ptr<PointLocatorBase> point_locator;
187  std::unique_ptr<const Elem> side_ptr;
188 
189  if (integrate_slits)
190  point_locator = mesh.sub_point_locator();
191 
192  // Iterate over all the active elements in the mesh
193  // that live on this processor.
194  for (const auto & e : mesh.active_local_element_ptr_range())
195  {
196  const dof_id_type e_id = e->id();
197 
198 #ifdef LIBMESH_ENABLE_AMR
199 
200  if (e->infinite())
201  {
202  libmesh_warning("Warning: Jumps on the border of infinite elements are ignored."
203  << std::endl);
204  continue;
205  }
206 
207  // See if the parent of element e has been examined yet;
208  // if not, we may want to compute the estimator on it
209  const Elem * parent = e->parent();
210 
211  // We only can compute and only need to compute on
212  // parents with all active children
213  bool compute_on_parent = true;
214  if (!parent || !estimate_parent_error)
215  compute_on_parent = false;
216  else
217  for (auto & child : parent->child_ref_range())
218  if (!child.active())
219  compute_on_parent = false;
220 
221  if (compute_on_parent &&
222  !error_per_cell[parent->id()])
223  {
224  // Compute a projection onto the parent
225  DenseVector<Number> Uparent;
227  (*(system.solution), dof_map, parent, Uparent, false);
228 
229  // Loop over the neighbors of the parent
230  for (auto n_p : parent->side_index_range())
231  {
232  if (parent->neighbor_ptr(n_p) != nullptr) // parent has a neighbor here
233  {
234  // Find the active neighbors in this direction
235  std::vector<const Elem *> active_neighbors;
236  parent->neighbor_ptr(n_p)->
237  active_family_tree_by_neighbor(active_neighbors,
238  parent);
239  // Compute the flux to each active neighbor
240  for (std::size_t a=0,
241  n_active_neighbors = active_neighbors.size();
242  a != n_active_neighbors; ++a)
243  {
244  const Elem * f = active_neighbors[a];
245 
246  if (f ->infinite()) // don't take infinite elements into account
247  continue;
248 
249  // FIXME - what about when f->level <
250  // parent->level()??
251  if (f->level() >= parent->level())
252  {
253  fine_context->pre_fe_reinit(system, f);
254  coarse_context->pre_fe_reinit(system, parent);
255  libmesh_assert_equal_to
256  (coarse_context->get_elem_solution().size(),
257  Uparent.size());
258  coarse_context->get_elem_solution() = Uparent;
259 
260  this->reinit_sides();
261 
262  // Loop over all significant variables in the system
263  for (var=0; var<n_vars; var++)
264  if (error_norm.weight(var) != 0.0 &&
265  system.variable_type(var).family != SCALAR)
266  {
268 
269  error_per_cell[fine_context->get_elem().id()] +=
270  static_cast<ErrorVectorReal>(fine_error);
271  error_per_cell[coarse_context->get_elem().id()] +=
272  static_cast<ErrorVectorReal>(coarse_error);
273  }
274 
275  // Keep track of the number of internal flux
276  // sides found on each element
278  {
279  n_flux_faces[fine_context->get_elem().id()]++;
280  n_flux_faces[coarse_context->get_elem().id()] +=
282  }
283  }
284  }
285  }
286  else if (integrate_boundary_sides)
287  {
288  fine_context->pre_fe_reinit(system, parent);
289  libmesh_assert_equal_to
290  (fine_context->get_elem_solution().size(),
291  Uparent.size());
292  fine_context->get_elem_solution() = Uparent;
293  fine_context->side = cast_int<unsigned char>(n_p);
294  fine_context->side_fe_reinit();
295 
296  // If we find a boundary flux for any variable,
297  // let's just count it as a flux face for all
298  // variables. Otherwise we'd need to keep track of
299  // a separate n_flux_faces and error_per_cell for
300  // every single var.
301  bool found_boundary_flux = false;
302 
303  for (var=0; var<n_vars; var++)
304  if (error_norm.weight(var) != 0.0 &&
305  system.variable_type(var).family != SCALAR)
306  {
307  if (this->boundary_side_integration())
308  {
309  error_per_cell[fine_context->get_elem().id()] +=
310  static_cast<ErrorVectorReal>(fine_error);
311  found_boundary_flux = true;
312  }
313  }
314 
315  if (scale_by_n_flux_faces && found_boundary_flux)
316  n_flux_faces[fine_context->get_elem().id()]++;
317  }
318  }
319  }
320 #endif // #ifdef LIBMESH_ENABLE_AMR
321 
322  // If we do any more flux integration, e will be the fine element
323  fine_context->pre_fe_reinit(system, e);
324 
325  // Loop over the neighbors of element e
326  for (auto n_e : e->side_index_range())
327  {
328  if ((e->neighbor_ptr(n_e) != nullptr) ||
330  {
331  fine_context->side = cast_int<unsigned char>(n_e);
332  fine_context->side_fe_reinit();
333  }
334 
335  // e is not on the boundary (infinite elements are treated as boundary)
336  if (e->neighbor_ptr(n_e) != nullptr
337  && !e->neighbor_ptr(n_e) ->infinite())
338  {
339 
340  const Elem * f = e->neighbor_ptr(n_e);
341  const dof_id_type f_id = f->id();
342 
343  // Compute flux jumps if we are in case 1 or case 2.
344  if ((f->active() && (f->level() == e->level()) && (e_id < f_id))
345  || (f->level() < e->level()))
346  {
347  // f is now the coarse element
348  coarse_context->pre_fe_reinit(system, f);
349 
350  this->reinit_sides();
351 
352  // Loop over all significant variables in the system
353  for (var=0; var<n_vars; var++)
354  if (error_norm.weight(var) != 0.0 &&
355  system.variable_type(var).family != SCALAR)
356  {
358 
359  error_per_cell[fine_context->get_elem().id()] +=
360  static_cast<ErrorVectorReal>(fine_error);
361  error_per_cell[coarse_context->get_elem().id()] +=
362  static_cast<ErrorVectorReal>(coarse_error);
363  }
364 
365  // Keep track of the number of internal flux
366  // sides found on each element
368  {
369  n_flux_faces[fine_context->get_elem().id()]++;
370  n_flux_faces[coarse_context->get_elem().id()] +=
372  }
373  } // end if (case1 || case2)
374  } // if (e->neighbor(n_e) != nullptr)
375 
376  // e might not have a neighbor_ptr, but might still have
377  // another element sharing its side. This can happen in a
378  // mesh where solution continuity is maintained via nodal
379  // constraint rows.
380  else if (integrate_slits)
381  {
382  side_ptr = e->build_side_ptr(n_e);
383  std::set<const Elem *> candidate_elements;
384  (*point_locator)(side_ptr->vertex_average(), candidate_elements);
385 
386  // We should have at least found ourselves...
387  libmesh_assert(candidate_elements.count(e));
388 
389  // If we only found ourselves, this probably isn't a
390  // slit; we don't yet support meshes so non-conforming
391  // as to have overlap of part of an element side without
392  // overlap of its center.
393  if (candidate_elements.size() < 2)
394  continue;
395 
396  FEType hardest_fe_type = fine_context->find_hardest_fe_type();
397 
398  auto dim = e->dim();
399 
400  auto side_qrule =
401  hardest_fe_type.unweighted_quadrature_rule
402  (dim-1, system.extra_quadrature_order);
403  auto side_fe = FEAbstract::build(dim, hardest_fe_type);
404  side_fe->attach_quadrature_rule(side_qrule.get());
405  const std::vector<Point> & qface_point = side_fe->get_xyz();
406  side_fe->reinit(e, n_e);
407 
408  for (auto qp : make_range(side_qrule->n_points()))
409  {
410  const Point p = qface_point[qp];
411  const std::vector<Point> qp_pointvec(1, p);
412  std::set<const Elem *> side_elements;
413  (*point_locator)(side_ptr->vertex_average(), side_elements);
414 
415  // If we have multiple neighbors meeting here we'll just
416  // take weighted jumps from all of them.
417  //
418  // We'll also do integrations from both sides of slits,
419  // rather than try to figure out a disambiguation rule
420  // that makes sense for non-conforming slits in general.
421  // This means we want an extra factor of 0.5 on the
422  // integrals to compensate for doubling them.
423  const std::size_t n_neighbors = side_elements.size() - 1;
424  const Real neighbor_frac = Real(1)/n_neighbors;
425 
426  const std::vector<Real>
427  qp_weightvec(1, neighbor_frac * side_qrule->w(qp));
428 
429  for (const Elem * f : side_elements)
430  {
431  if (f == e)
432  continue;
433 
434  coarse_context->pre_fe_reinit(system, f);
435  fine_context->pre_fe_reinit(system, e);
436  std::vector<Point> qp_coarse, qp_fine;
437  for (unsigned int v=0; v<n_vars; v++)
438  if (error_norm.weight(v) != 0.0 &&
439  fine_context->get_system().variable_type(v).family != SCALAR)
440  {
441  FEBase * coarse_fe = coarse_context->get_side_fe(v, dim);
442  if (qp_coarse.empty())
443  FEMap::inverse_map (dim, f, qp_pointvec, qp_coarse);
444  coarse_fe->reinit(f, &qp_coarse, &qp_weightvec);
445  FEBase * fine_fe = fine_context->get_side_fe(v, dim);
446  if (qp_fine.empty())
447  FEMap::inverse_map (dim, e, qp_pointvec, qp_fine);
448  fine_fe->reinit(e, &qp_fine, &qp_weightvec);
449  }
450 
451  // Loop over all significant variables in the system
452  for (var=0; var<n_vars; var++)
453  if (error_norm.weight(var) != 0.0 &&
454  system.variable_type(var).family != SCALAR)
455  {
457 
458  error_per_cell[fine_context->get_elem().id()] +=
459  static_cast<ErrorVectorReal>(fine_error);
460  error_per_cell[coarse_context->get_elem().id()] +=
461  static_cast<ErrorVectorReal>(coarse_error);
462  }
463  }
464  }
465  }
466 
467  // Otherwise, e is on the boundary. If it happens to
468  // be on a Dirichlet boundary, we need not do anything.
469  // On the other hand, if e is on a Neumann (flux) boundary
470  // with grad(u).n = g, we need to compute the additional residual
471  // (h * \int |g - grad(u_h).n|^2 dS)^(1/2).
472  // We can only do this with some knowledge of the boundary
473  // conditions, i.e. the user must have attached an appropriate
474  // BC function.
475  else if (integrate_boundary_sides)
476  {
477  if (integrate_slits)
478  libmesh_not_implemented();
479 
480  bool found_boundary_flux = false;
481 
482  for (var=0; var<n_vars; var++)
483  if (error_norm.weight(var) != 0.0 &&
484  system.variable_type(var).family != SCALAR)
485  if (this->boundary_side_integration())
486  {
487  error_per_cell[fine_context->get_elem().id()] +=
488  static_cast<ErrorVectorReal>(fine_error);
489  found_boundary_flux = true;
490  }
491 
492  if (scale_by_n_flux_faces && found_boundary_flux)
493  n_flux_faces[fine_context->get_elem().id()]++;
494  } // end if (e->neighbor_ptr(n_e) == nullptr)
495  } // end loop over neighbors
496  } // End loop over active local elements
497 
498 
499  // Each processor has now computed the error contributions
500  // for its local elements. We need to sum the vector
501  // and then take the square-root of each component. Note
502  // that we only need to sum if we are running on multiple
503  // processors, and we only need to take the square-root
504  // if the value is nonzero. There will in general be many
505  // zeros for the inactive elements.
506 
507  // First sum the vector of estimated error values
508  this->reduce_error(error_per_cell, system.comm());
509 
510  // Compute the square-root of each component.
511  for (auto i : index_range(error_per_cell))
512  if (error_per_cell[i] != 0.)
513  error_per_cell[i] = std::sqrt(error_per_cell[i]);
514 
515 
516  if (this->scale_by_n_flux_faces)
517  {
518  // Sum the vector of flux face counts
519  this->reduce_error(n_flux_faces, system.comm());
520 
521  // Sanity check: Make sure the number of flux faces is
522  // always an integer value
523 #ifdef DEBUG
524  for (const auto & val : n_flux_faces)
525  libmesh_assert_equal_to (val, static_cast<float>(static_cast<unsigned int>(val)));
526 #endif
527 
528  // Scale the error by the number of flux faces for each element
529  for (auto i : index_range(n_flux_faces))
530  {
531  if (n_flux_faces[i] == 0.0) // inactive or non-local element
532  continue;
533 
534  error_per_cell[i] /= static_cast<ErrorVectorReal>(n_flux_faces[i]);
535  }
536  }
537 
538  // If we used a non-standard solution before, now is the time to fix
539  // the current_local_solution
540  if (solution_vector && solution_vector != system.solution.get())
541  {
542  NumericVector<Number> * newsol =
543  const_cast<NumericVector<Number> *>(solution_vector);
544  System & sys = const_cast<System &>(system);
545  newsol->swap(*sys.solution);
546  sys.update();
547  }
548 }
SystemNorm error_norm
When estimating the error in a single system, the error_norm is used to control the scaling and norm ...
void active_family_tree_by_neighbor(T elem, std::vector< T > &family, T neighbor_in, bool reset=true)
std::unique_ptr< FEMContext > fine_context
Context objects for integrating on the fine and coarse elements sharing a face.
unsigned int dim
static Point inverse_map(const unsigned int dim, const Elem *elem, const Point &p, const Real tolerance=TOLERANCE, const bool secure=true, const bool extra_checks=true)
Definition: fe_map.C:1628
MeshBase & mesh
virtual void internal_side_integration()=0
The function, to be implemented by derived classes, which calculates an error term on an internal sid...
bool integrate_slits
A boolean flag, by default false, to be set to true if integrations should be performed on "slits" wh...
DIE A HORRIBLE DEATH HERE typedef float ErrorVectorReal
static void coarsened_dof_values(const NumericVector< Number > &global_vector, const DofMap &dof_map, const Elem *coarse_elem, DenseVector< Number > &coarse_dofs, const unsigned int var, const bool use_old_dof_indices=false)
Creates a local projection on coarse_elem, based on the DoF values in global_vector for it&#39;s children...
Definition: fe_base.C:1012
void reinit_sides()
A utility function to reinit the finite element data on elements sharing a side.
void libmesh_ignore(const Args &...)
virtual bool boundary_side_integration()
The function, to be implemented by derived classes, which calculates an error term on a boundary side...
unsigned int n_vars
libmesh_assert(ctx)
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 th...
std::unique_ptr< FEMContext > coarse_context
FEGenericBase< Real > FEBase
Real fine_error
The fine and coarse error values to be set by each side_integration();.
Real weight(unsigned int var) const
Definition: system_norm.C:134
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
static std::unique_ptr< FEAbstract > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.
Definition: fe_abstract.C:78
float coarse_n_flux_faces_increment()
A utility function to correctly increase n_flux_faces for the coarse element.
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:140
bool use_unweighted_quadrature_rules
This boolean flag allows you to use "unweighted" quadrature rules (sized to exactly integrate unweigh...
unsigned int var
The variable number currently being evaluated.
bool scale_by_n_flux_faces
This boolean flag allows you to scale the error indicator result for each element by the number of "f...
bool integrate_boundary_sides
A boolean flag, by default false, to be set to true if integrations with boundary_side_integration() ...
template class LIBMESH_EXPORT NumericVector< Number >
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
Definition: int_range.h:117
virtual void init_context(FEMContext &c)
An initialization function, to give derived classes a chance to request specific data from the FE obj...
uint8_t dof_id_type
Definition: id_types.h:67

◆ 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().

52 {
53  SystemNorm old_error_norm = this->error_norm;
54 
55  // Sum the error values from each system
56  for (auto s : make_range(equation_systems.n_systems()))
57  {
58  ErrorVector system_error_per_cell;
59  const System & sys = equation_systems.get_system(s);
60  if (const auto it = error_norms.find(&sys);
61  it == error_norms.end())
62  this->error_norm = old_error_norm;
63  else
64  this->error_norm = it->second;
65 
66  const NumericVector<Number> * solution_vector = nullptr;
67  if (solution_vectors &&
68  solution_vectors->find(&sys) != solution_vectors->end())
69  solution_vector = solution_vectors->find(&sys)->second;
70 
71  this->estimate_error(sys, system_error_per_cell,
72  solution_vector, estimate_parent_error);
73 
74  if (s)
75  {
76  libmesh_assert_equal_to (error_per_cell.size(), system_error_per_cell.size());
77  for (auto i : index_range(error_per_cell))
78  error_per_cell[i] += system_error_per_cell[i];
79  }
80  else
81  error_per_cell = system_error_per_cell;
82  }
83 
84  // Restore our old state before returning
85  this->error_norm = old_error_norm;
86 }
SystemNorm error_norm
When estimating the error in a single system, the error_norm is used to control the scaling and norm ...
virtual void estimate_error(const System &system, ErrorVector &error_per_cell, const NumericVector< Number > *solution_vector=nullptr, bool estimate_parent_error=false)=0
This pure virtual function must be redefined in derived classes to compute the error for each active ...
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:140
template class LIBMESH_EXPORT NumericVector< Number >
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
Definition: int_range.h:117

◆ 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().

98 {
99  SystemNorm old_error_norm = this->error_norm;
100 
101  // Find the requested error values from each system
102  for (auto s : make_range(equation_systems.n_systems()))
103  {
104  const System & sys = equation_systems.get_system(s);
105 
106  unsigned int n_vars = sys.n_vars();
107 
108  for (unsigned int v = 0; v != n_vars; ++v)
109  {
110  // Only fill in ErrorVectors the user asks for
111  if (!errors_per_cell.count(std::make_pair(&sys, v)))
112  continue;
113 
114  // Calculate error in only one variable
115  std::vector<Real> weights(n_vars, 0.0);
116  weights[v] = 1.0;
117  this->error_norm =
118  SystemNorm(std::vector<FEMNormType>(n_vars, old_error_norm.type(v)),
119  weights);
120 
121  const NumericVector<Number> * solution_vector = nullptr;
122  if (solution_vectors &&
123  solution_vectors->find(&sys) != solution_vectors->end())
124  solution_vector = solution_vectors->find(&sys)->second;
125 
126  this->estimate_error
127  (sys, *errors_per_cell[std::make_pair(&sys, v)],
128  solution_vector, estimate_parent_error);
129  }
130  }
131 
132  // Restore our old state before returning
133  this->error_norm = old_error_norm;
134 }
SystemNorm error_norm
When estimating the error in a single system, the error_norm is used to control the scaling and norm ...
virtual void estimate_error(const System &system, ErrorVector &error_per_cell, const NumericVector< Number > *solution_vector=nullptr, bool estimate_parent_error=false)=0
This pure virtual function must be redefined in derived classes to compute the error for each active ...
unsigned int n_vars
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:140
template class LIBMESH_EXPORT NumericVector< Number >

◆ 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().

49 {
50  // Derived classes are *supposed* to rederive this
51  libmesh_deprecated();
52 }

◆ 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().

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

◆ 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().

554 {
555  fine_context->side_fe_reinit();
556 
557  unsigned short dim = fine_context->get_elem().dim();
558  libmesh_assert_equal_to(dim, coarse_context->get_elem().dim());
559 
560  FEBase * fe_fine = nullptr;
561  fine_context->get_side_fe( 0, fe_fine, dim );
562 
563  // Get the physical locations of the fine element quadrature points
564  std::vector<Point> qface_point = fe_fine->get_xyz();
565 
566  // Find the master quadrature point locations on the coarse element
567  FEBase * fe_coarse = nullptr;
568  coarse_context->get_side_fe( 0, fe_coarse, dim );
569 
570  std::vector<Point> qp_coarse;
571 
572  FEMap::inverse_map (coarse_context->get_elem().dim(),
573  &coarse_context->get_elem(), qface_point,
574  qp_coarse);
575 
576  // The number of variables in the system
577  const unsigned int n_vars = fine_context->n_vars();
578 
579  // Calculate all coarse element shape functions at those locations
580  for (unsigned int v=0; v<n_vars; v++)
581  if (error_norm.weight(v) != 0.0 &&
582  fine_context->get_system().variable_type(v).family != SCALAR)
583  {
584  coarse_context->get_side_fe( v, fe_coarse, dim );
585  fe_coarse->reinit (&coarse_context->get_elem(), &qp_coarse);
586  }
587 }
SystemNorm error_norm
When estimating the error in a single system, the error_norm is used to control the scaling and norm ...
std::unique_ptr< FEMContext > fine_context
Context objects for integrating on the fine and coarse elements sharing a face.
unsigned int dim
static Point inverse_map(const unsigned int dim, const Elem *elem, const Point &p, const Real tolerance=TOLERANCE, const bool secure=true, const bool extra_checks=true)
Definition: fe_map.C:1628
unsigned int n_vars
std::unique_ptr< FEMContext > coarse_context
FEGenericBase< Real > FEBase
Real weight(unsigned int var) const
Definition: system_norm.C:134

◆ type()

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

Member Data Documentation

◆ coarse_context

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

◆ coarse_error

Real libMesh::JumpErrorEstimator::coarse_error
protected

◆ 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

◆ fine_error

Real libMesh::JumpErrorEstimator::fine_error
protected

◆ 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 documentation for this class was generated from the following files: