libMesh::KellyErrorEstimator Class Reference

#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
 
KellyErrorEstimatoroperator= (const KellyErrorEstimator &)=delete
 
 KellyErrorEstimator (KellyErrorEstimator &&)=default
 
KellyErrorEstimatoroperator= (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< FEMContextfine_context
 
std::unique_ptr< FEMContextcoarse_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.

41  :
43  _bc_function(nullptr)
44 {
46 }
std::pair< bool, Real >(* _bc_function)(const System &system, const Point &p, const std::string &var_name)

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

210 {
211  _bc_function = fptr;
212 
213  // We may be turning boundary side integration on or off
214  if (fptr)
216  else
217  integrate_boundary_sides = false;
218 }
std::pair< bool, Real >(* _bc_function)(const System &system, const Point &p, const std::string &var_name)

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

136 {
137  const Elem & fine_elem = fine_context->get_elem();
138 
139  FEBase * fe_fine = nullptr;
140  fine_context->get_side_fe( var, fe_fine, fine_elem.dim() );
141 
142  const std::string & var_name =
143  fine_context->get_system().variable_name(var);
144 
145  std::vector<std::vector<RealGradient>> dphi_fine = fe_fine->get_dphi();
146  std::vector<Point> face_normals = fe_fine->get_normals();
147  std::vector<Real> JxW_face = fe_fine->get_JxW();
148  std::vector<Point> qface_point = fe_fine->get_xyz();
149 
150  // The reinitialization also recomputes the locations of
151  // the quadrature points on the side. By checking if the
152  // first quadrature point on the side is on a flux boundary
153  // for a particular variable, we will determine if the whole
154  // element is on a flux boundary (assuming quadrature points
155  // are strictly contained in the side).
156  if (this->_bc_function(fine_context->get_system(),
157  qface_point[0], var_name).first)
158  {
159  const Real h = fine_elem.hmax();
160 
161  // The number of quadrature points
162  const unsigned int n_qp = fe_fine->n_quadrature_points();
163 
164  // The error contribution from this face
165  Real error = 1.e-30;
166 
167  // loop over the integration points on the face.
168  for (unsigned int qp=0; qp<n_qp; qp++)
169  {
170  // Value of the imposed flux BC at this quadrature point.
171  const std::pair<bool,Real> flux_bc =
172  this->_bc_function(fine_context->get_system(),
173  qface_point[qp], var_name);
174 
175  // Be sure the BC function still thinks we're on the
176  // flux boundary.
177  libmesh_assert_equal_to (flux_bc.first, true);
178 
179  // The solution gradient from each element
180  Gradient grad_fine = fine_context->side_gradient(var, qp);
181 
182  // The difference between the desired BC and the approximate solution.
183  const Number jump = flux_bc.second - grad_fine*face_normals[qp];
184 
185  // The flux jump squared. If using complex numbers,
186  // TensorTools::norm_sq(z) returns |z|^2, where |z| is the modulus of z.
187  const Real jump2 = TensorTools::norm_sq(jump);
188 
189  // Integrate the error on the face. The error is
190  // scaled by an additional power of h, where h is
191  // the maximum side length for the element. This
192  // arises in the definition of the indicator.
193  error += JxW_face[qp]*jump2;
194 
195  } // End quadrature point loop
196 
197  fine_error = error*h*error_norm.weight(var);
198 
199  return true;
200  } // end if side on flux boundary
201  return false;
202 }
std::pair< bool, Real >(* _bc_function)(const System &system, const Point &p, const std::string &var_name)
std::unique_ptr< FEMContext > fine_context
NumberVectorValue Gradient
FEGenericBase< Real > FEBase
Real weight(unsigned int var) const
Definition: system_norm.C:132
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real

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

463 {
464  // Keep track of the number of internal flux sides found on each
465  // element
466  unsigned short dim = coarse_context->get_elem().dim();
467 
468  const unsigned int divisor =
469  1 << (dim-1)*(fine_context->get_elem().level() -
470  coarse_context->get_elem().level());
471 
472  // With a difference of n levels between fine and coarse elements,
473  // we compute a fractional flux face for the coarse element by adding:
474  // 1/2^n in 2D
475  // 1/4^n in 3D
476  // each time. This code will get hit 2^n times in 2D and 4^n
477  // times in 3D so that the final flux face count for the coarse
478  // element will be an integer value.
479 
480  return 1.0f / static_cast<float>(divisor);
481 }
std::unique_ptr< FEMContext > fine_context
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 
)
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.

References libMesh::Elem::active(), libMesh::JumpErrorEstimator::boundary_side_integration(), libMesh::Elem::child_ref_range(), libMesh::JumpErrorEstimator::coarse_context, libMesh::JumpErrorEstimator::coarse_error, libMesh::JumpErrorEstimator::coarse_n_flux_faces_increment(), libMesh::FEGenericBase< OutputType >::coarsened_dof_values(), libMesh::ParallelObject::comm(), libMesh::ErrorEstimator::error_norm, libMesh::ErrorVectorReal, libMesh::FEType::family, libMesh::JumpErrorEstimator::fine_context, libMesh::JumpErrorEstimator::fine_error, libMesh::System::get_dof_map(), libMesh::System::get_mesh(), libMesh::FEAbstract::get_xyz(), libMesh::DofObject::id(), libMesh::index_range(), libMesh::JumpErrorEstimator::init_context(), libMesh::JumpErrorEstimator::integrate_boundary_sides, libMesh::JumpErrorEstimator::internal_side_integration(), libMesh::Elem::level(), libMesh::libmesh_ignore(), mesh, n_vars, libMesh::System::n_vars(), libMesh::Elem::neighbor_ptr(), libMesh::Elem::parent(), libMesh::ErrorEstimator::reduce_error(), libMesh::JumpErrorEstimator::reinit_sides(), libMesh::SCALAR, libMesh::JumpErrorEstimator::scale_by_n_flux_faces, libMesh::Elem::side_index_range(), libMesh::DenseVector< T >::size(), libMesh::System::solution, libMesh::NumericVector< T >::swap(), libMesh::JumpErrorEstimator::var, libMesh::System::variable_type(), and libMesh::SystemNorm::weight().

57 {
58  LOG_SCOPE("estimate_error()", "JumpErrorEstimator");
59 
96  // This parameter is not used when !LIBMESH_ENABLE_AMR.
97  libmesh_ignore(estimate_parent_error);
98 
99  // The current mesh
100  const MeshBase & mesh = system.get_mesh();
101 
102  // The number of variables in the system
103  const unsigned int n_vars = system.n_vars();
104 
105  // The DofMap for this system
106 #ifdef LIBMESH_ENABLE_AMR
107  const DofMap & dof_map = system.get_dof_map();
108 #endif
109 
110  // Resize the error_per_cell vector to be
111  // the number of elements, initialize it to 0.
112  error_per_cell.resize (mesh.max_elem_id());
113  std::fill (error_per_cell.begin(), error_per_cell.end(), 0.);
114 
115  // Declare a vector of floats which is as long as
116  // error_per_cell above, and fill with zeros. This vector will be
117  // used to keep track of the number of edges (faces) on each active
118  // element which are either:
119  // 1) an internal edge
120  // 2) an edge on a Neumann boundary for which a boundary condition
121  // function has been specified.
122  // The error estimator can be scaled by the number of flux edges (faces)
123  // which the element actually has to obtain a more uniform measure
124  // of the error. Use floats instead of ints since in case 2 (above)
125  // f gets 1/2 of a flux face contribution from each of his
126  // neighbors
127  std::vector<float> n_flux_faces;
129  n_flux_faces.resize(error_per_cell.size(), 0);
130 
131  // Prepare current_local_solution to localize a non-standard
132  // solution vector if necessary
133  if (solution_vector && solution_vector != system.solution.get())
134  {
135  NumericVector<Number> * newsol =
136  const_cast<NumericVector<Number> *>(solution_vector);
137  System & sys = const_cast<System &>(system);
138  newsol->swap(*sys.solution);
139  sys.update();
140  }
141 
142  fine_context.reset(new FEMContext(system));
143  coarse_context.reset(new FEMContext(system));
144 
145  // Loop over all the variables we've been requested to find jumps in, to
146  // pre-request
147  for (var=0; var<n_vars; var++)
148  {
149  // Skip variables which aren't part of our norm,
150  // as well as SCALAR variables, which have no jumps
151  if (error_norm.weight(var) == 0.0 ||
152  system.variable_type(var).family == SCALAR)
153  continue;
154 
155  // FIXME: Need to generalize this to vector-valued elements. [PB]
156  FEBase * side_fe = nullptr;
157 
158  const std::set<unsigned char> & elem_dims =
159  fine_context->elem_dimensions();
160 
161  for (const auto & dim : elem_dims)
162  {
163  fine_context->get_side_fe( var, side_fe, dim );
164 
165  side_fe->get_xyz();
166  }
167  }
168 
169  this->init_context(*fine_context);
171 
172  // Iterate over all the active elements in the mesh
173  // that live on this processor.
174  for (const auto & e : mesh.active_local_element_ptr_range())
175  {
176  const dof_id_type e_id = e->id();
177 
178 #ifdef LIBMESH_ENABLE_AMR
179  // See if the parent of element e has been examined yet;
180  // if not, we may want to compute the estimator on it
181  const Elem * parent = e->parent();
182 
183  // We only can compute and only need to compute on
184  // parents with all active children
185  bool compute_on_parent = true;
186  if (!parent || !estimate_parent_error)
187  compute_on_parent = false;
188  else
189  for (auto & child : parent->child_ref_range())
190  if (!child.active())
191  compute_on_parent = false;
192 
193  if (compute_on_parent &&
194  !error_per_cell[parent->id()])
195  {
196  // Compute a projection onto the parent
197  DenseVector<Number> Uparent;
199  (*(system.solution), dof_map, parent, Uparent, false);
200 
201  // Loop over the neighbors of the parent
202  for (auto n_p : parent->side_index_range())
203  {
204  if (parent->neighbor_ptr(n_p) != nullptr) // parent has a neighbor here
205  {
206  // Find the active neighbors in this direction
207  std::vector<const Elem *> active_neighbors;
208  parent->neighbor_ptr(n_p)->
209  active_family_tree_by_neighbor(active_neighbors,
210  parent);
211  // Compute the flux to each active neighbor
212  for (std::size_t a=0,
213  n_active_neighbors = active_neighbors.size();
214  a != n_active_neighbors; ++a)
215  {
216  const Elem * f = active_neighbors[a];
217  // FIXME - what about when f->level <
218  // parent->level()??
219  if (f->level() >= parent->level())
220  {
221  fine_context->pre_fe_reinit(system, f);
222  coarse_context->pre_fe_reinit(system, parent);
223  libmesh_assert_equal_to
224  (coarse_context->get_elem_solution().size(),
225  Uparent.size());
226  coarse_context->get_elem_solution() = Uparent;
227 
228  this->reinit_sides();
229 
230  // Loop over all significant variables in the system
231  for (var=0; var<n_vars; var++)
232  if (error_norm.weight(var) != 0.0 &&
233  system.variable_type(var).family != SCALAR)
234  {
236 
237  error_per_cell[fine_context->get_elem().id()] +=
238  static_cast<ErrorVectorReal>(fine_error);
239  error_per_cell[coarse_context->get_elem().id()] +=
240  static_cast<ErrorVectorReal>(coarse_error);
241  }
242 
243  // Keep track of the number of internal flux
244  // sides found on each element
246  {
247  n_flux_faces[fine_context->get_elem().id()]++;
248  n_flux_faces[coarse_context->get_elem().id()] +=
250  }
251  }
252  }
253  }
254  else if (integrate_boundary_sides)
255  {
256  fine_context->pre_fe_reinit(system, parent);
257  libmesh_assert_equal_to
258  (fine_context->get_elem_solution().size(),
259  Uparent.size());
260  fine_context->get_elem_solution() = Uparent;
261  fine_context->side = cast_int<unsigned char>(n_p);
262  fine_context->side_fe_reinit();
263 
264  // If we find a boundary flux for any variable,
265  // let's just count it as a flux face for all
266  // variables. Otherwise we'd need to keep track of
267  // a separate n_flux_faces and error_per_cell for
268  // every single var.
269  bool found_boundary_flux = false;
270 
271  for (var=0; var<n_vars; var++)
272  if (error_norm.weight(var) != 0.0 &&
273  system.variable_type(var).family != SCALAR)
274  {
275  if (this->boundary_side_integration())
276  {
277  error_per_cell[fine_context->get_elem().id()] +=
278  static_cast<ErrorVectorReal>(fine_error);
279  found_boundary_flux = true;
280  }
281  }
282 
283  if (scale_by_n_flux_faces && found_boundary_flux)
284  n_flux_faces[fine_context->get_elem().id()]++;
285  }
286  }
287  }
288 #endif // #ifdef LIBMESH_ENABLE_AMR
289 
290  // If we do any more flux integration, e will be the fine element
291  fine_context->pre_fe_reinit(system, e);
292 
293  // Loop over the neighbors of element e
294  for (auto n_e : e->side_index_range())
295  {
296  if ((e->neighbor_ptr(n_e) != nullptr) ||
298  {
299  fine_context->side = cast_int<unsigned char>(n_e);
300  fine_context->side_fe_reinit();
301  }
302 
303  if (e->neighbor_ptr(n_e) != nullptr) // e is not on the boundary
304  {
305  const Elem * f = e->neighbor_ptr(n_e);
306  const dof_id_type f_id = f->id();
307 
308  // Compute flux jumps if we are in case 1 or case 2.
309  if ((f->active() && (f->level() == e->level()) && (e_id < f_id))
310  || (f->level() < e->level()))
311  {
312  // f is now the coarse element
313  coarse_context->pre_fe_reinit(system, f);
314 
315  this->reinit_sides();
316 
317  // Loop over all significant variables in the system
318  for (var=0; var<n_vars; var++)
319  if (error_norm.weight(var) != 0.0 &&
320  system.variable_type(var).family != SCALAR)
321  {
323 
324  error_per_cell[fine_context->get_elem().id()] +=
325  static_cast<ErrorVectorReal>(fine_error);
326  error_per_cell[coarse_context->get_elem().id()] +=
327  static_cast<ErrorVectorReal>(coarse_error);
328  }
329 
330  // Keep track of the number of internal flux
331  // sides found on each element
333  {
334  n_flux_faces[fine_context->get_elem().id()]++;
335  n_flux_faces[coarse_context->get_elem().id()] +=
337  }
338  } // end if (case1 || case2)
339  } // if (e->neighbor(n_e) != nullptr)
340 
341  // Otherwise, e is on the boundary. If it happens to
342  // be on a Dirichlet boundary, we need not do anything.
343  // On the other hand, if e is on a Neumann (flux) boundary
344  // with grad(u).n = g, we need to compute the additional residual
345  // (h * \int |g - grad(u_h).n|^2 dS)^(1/2).
346  // We can only do this with some knowledge of the boundary
347  // conditions, i.e. the user must have attached an appropriate
348  // BC function.
349  else if (integrate_boundary_sides)
350  {
351  bool found_boundary_flux = false;
352 
353  for (var=0; var<n_vars; var++)
354  if (error_norm.weight(var) != 0.0 &&
355  system.variable_type(var).family != SCALAR)
356  if (this->boundary_side_integration())
357  {
358  error_per_cell[fine_context->get_elem().id()] +=
359  static_cast<ErrorVectorReal>(fine_error);
360  found_boundary_flux = true;
361  }
362 
363  if (scale_by_n_flux_faces && found_boundary_flux)
364  n_flux_faces[fine_context->get_elem().id()]++;
365  } // end if (e->neighbor_ptr(n_e) == nullptr)
366  } // end loop over neighbors
367  } // End loop over active local elements
368 
369 
370  // Each processor has now computed the error contributions
371  // for its local elements. We need to sum the vector
372  // and then take the square-root of each component. Note
373  // that we only need to sum if we are running on multiple
374  // processors, and we only need to take the square-root
375  // if the value is nonzero. There will in general be many
376  // zeros for the inactive elements.
377 
378  // First sum the vector of estimated error values
379  this->reduce_error(error_per_cell, system.comm());
380 
381  // Compute the square-root of each component.
382  for (auto i : index_range(error_per_cell))
383  if (error_per_cell[i] != 0.)
384  error_per_cell[i] = std::sqrt(error_per_cell[i]);
385 
386 
387  if (this->scale_by_n_flux_faces)
388  {
389  // Sum the vector of flux face counts
390  this->reduce_error(n_flux_faces, system.comm());
391 
392  // Sanity check: Make sure the number of flux faces is
393  // always an integer value
394 #ifdef DEBUG
395  for (const auto & val : n_flux_faces)
396  libmesh_assert_equal_to (val, static_cast<float>(static_cast<unsigned int>(val)));
397 #endif
398 
399  // Scale the error by the number of flux faces for each element
400  for (auto i : index_range(n_flux_faces))
401  {
402  if (n_flux_faces[i] == 0.0) // inactive or non-local element
403  continue;
404 
405  error_per_cell[i] /= static_cast<ErrorVectorReal>(n_flux_faces[i]);
406  }
407  }
408 
409  // If we used a non-standard solution before, now is the time to fix
410  // the current_local_solution
411  if (solution_vector && solution_vector != system.solution.get())
412  {
413  NumericVector<Number> * newsol =
414  const_cast<NumericVector<Number> *>(solution_vector);
415  System & sys = const_cast<System &>(system);
416  newsol->swap(*sys.solution);
417  sys.update();
418  }
419 }
std::unique_ptr< FEMContext > fine_context
MeshBase & mesh
IntRange< std::size_t > index_range(const std::vector< T > &vec)
Definition: int_range.h:104
virtual void internal_side_integration()=0
const unsigned int n_vars
Definition: tecplot_io.C:69
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)
Definition: fe_base.C:791
void libmesh_ignore(const Args &...)
void reduce_error(std::vector< ErrorVectorReal > &error_per_cell, const Parallel::Communicator &comm) const
std::unique_ptr< FEMContext > coarse_context
FEGenericBase< Real > FEBase
Real weight(unsigned int var) const
Definition: system_norm.C:132
virtual void init_context(FEMContext &c)
uint8_t dof_id_type
Definition: id_types.h:64

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

52 {
53  SystemNorm old_error_norm = this->error_norm;
54 
55  // Sum the error values from each system
56  for (unsigned int s = 0; s != equation_systems.n_systems(); ++s)
57  {
58  ErrorVector system_error_per_cell;
59  const System & sys = equation_systems.get_system(s);
60  if (error_norms.find(&sys) == error_norms.end())
61  this->error_norm = old_error_norm;
62  else
63  this->error_norm = error_norms.find(&sys)->second;
64 
65  const NumericVector<Number> * solution_vector = nullptr;
66  if (solution_vectors &&
67  solution_vectors->find(&sys) != solution_vectors->end())
68  solution_vector = solution_vectors->find(&sys)->second;
69 
70  this->estimate_error(sys, system_error_per_cell,
71  solution_vector, estimate_parent_error);
72 
73  if (s)
74  {
75  libmesh_assert_equal_to (error_per_cell.size(), system_error_per_cell.size());
76  for (auto i : index_range(error_per_cell))
77  error_per_cell[i] += system_error_per_cell[i];
78  }
79  else
80  error_per_cell = system_error_per_cell;
81  }
82 
83  // Restore our old state before returning
84  this->error_norm = old_error_norm;
85 }
IntRange< std::size_t > index_range(const std::vector< T > &vec)
Definition: int_range.h:104
virtual void estimate_error(const System &system, ErrorVector &error_per_cell, const NumericVector< Number > *solution_vector=nullptr, bool estimate_parent_error=false)=0

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

97 {
98  SystemNorm old_error_norm = this->error_norm;
99 
100  // Find the requested error values from each system
101  for (unsigned int s = 0; s != equation_systems.n_systems(); ++s)
102  {
103  const System & sys = equation_systems.get_system(s);
104 
105  unsigned int n_vars = sys.n_vars();
106 
107  for (unsigned int v = 0; v != n_vars; ++v)
108  {
109  // Only fill in ErrorVectors the user asks for
110  if (errors_per_cell.find(std::make_pair(&sys, v)) ==
111  errors_per_cell.end())
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 }
virtual void estimate_error(const System &system, ErrorVector &error_per_cell, const NumericVector< Number > *solution_vector=nullptr, bool estimate_parent_error=false)=0
const unsigned int n_vars
Definition: tecplot_io.C:69

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

60 {
61  const unsigned int n_vars = c.n_vars();
62  for (unsigned int v=0; v<n_vars; v++)
63  {
64  // Possibly skip this variable
65  if (error_norm.weight(v) == 0.0) continue;
66 
67  // FIXME: Need to generalize this to vector-valued elements. [PB]
68  FEBase * side_fe = nullptr;
69 
70  const std::set<unsigned char> & elem_dims =
71  c.elem_dimensions();
72 
73  for (const auto & dim : elem_dims)
74  {
75  fine_context->get_side_fe( v, side_fe, dim );
76 
77  // We'll need gradients on both sides for flux jump computation
78  side_fe->get_dphi();
79 
80  // But we only need normal vectors from one side
81  if (&c != coarse_context.get())
82  side_fe->get_normals();
83  }
84  }
85 }
std::unique_ptr< FEMContext > fine_context
const unsigned int n_vars
Definition: tecplot_io.C:69
std::unique_ptr< FEMContext > coarse_context
FEGenericBase< Real > FEBase
Real weight(unsigned int var) const
Definition: system_norm.C:132

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

91 {
92  const Elem & coarse_elem = coarse_context->get_elem();
93  const Elem & fine_elem = fine_context->get_elem();
94 
95  FEBase * fe_fine = nullptr;
96  fine_context->get_side_fe( var, fe_fine, fine_elem.dim() );
97 
98  FEBase * fe_coarse = nullptr;
99  coarse_context->get_side_fe( var, fe_coarse, fine_elem.dim() );
100 
101  Real error = 1.e-30;
102  unsigned int n_qp = fe_fine->n_quadrature_points();
103 
104  std::vector<std::vector<RealGradient>> dphi_coarse = fe_coarse->get_dphi();
105  std::vector<std::vector<RealGradient>> dphi_fine = fe_fine->get_dphi();
106  std::vector<Point> face_normals = fe_fine->get_normals();
107  std::vector<Real> JxW_face = fe_fine->get_JxW();
108 
109  for (unsigned int qp=0; qp != n_qp; ++qp)
110  {
111  // Calculate solution gradients on fine and coarse elements
112  // at this quadrature point
113  Gradient
114  grad_fine = fine_context->side_gradient(var, qp),
115  grad_coarse = coarse_context->side_gradient(var, qp);
116 
117  // Find the jump in the normal derivative
118  // at this quadrature point
119  const Number jump = (grad_fine - grad_coarse)*face_normals[qp];
120  const Real jump2 = TensorTools::norm_sq(jump);
121 
122  // Accumulate the jump integral
123  error += JxW_face[qp] * jump2;
124  }
125 
126  // Add the h-weighted jump integral to each error term
127  fine_error =
128  error * fine_elem.hmax() * error_norm.weight(var);
129  coarse_error =
130  error * coarse_elem.hmax() * error_norm.weight(var);
131 }
std::unique_ptr< FEMContext > fine_context
NumberVectorValue Gradient
std::unique_ptr< FEMContext > coarse_context
FEGenericBase< Real > FEBase
Real weight(unsigned int var) const
Definition: system_norm.C:132
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real

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

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

425 {
426  fine_context->side_fe_reinit();
427 
428  unsigned short dim = fine_context->get_elem().dim();
429  libmesh_assert_equal_to(dim, coarse_context->get_elem().dim());
430 
431  FEBase * fe_fine = nullptr;
432  fine_context->get_side_fe( 0, fe_fine, dim );
433 
434  // Get the physical locations of the fine element quadrature points
435  std::vector<Point> qface_point = fe_fine->get_xyz();
436 
437  // Find the master quadrature point locations on the coarse element
438  FEBase * fe_coarse = nullptr;
439  coarse_context->get_side_fe( 0, fe_coarse, dim );
440 
441  std::vector<Point> qp_coarse;
442 
444  (coarse_context->get_elem().dim(), fe_coarse->get_fe_type(),
445  &coarse_context->get_elem(), qface_point, qp_coarse);
446 
447  // The number of variables in the system
448  const unsigned int n_vars = fine_context->n_vars();
449 
450  // Calculate all coarse element shape functions at those locations
451  for (unsigned int v=0; v<n_vars; v++)
452  if (error_norm.weight(v) != 0.0 &&
453  fine_context->get_system().variable_type(v).family != SCALAR)
454  {
455  coarse_context->get_side_fe( v, fe_coarse, dim );
456  fe_coarse->reinit (&coarse_context->get_elem(), &qp_coarse);
457  }
458 }
std::unique_ptr< FEMContext > fine_context
const unsigned int n_vars
Definition: tecplot_io.C:69
static Point inverse_map(const unsigned int dim, const FEType &fe_t, const Elem *elem, const Point &p, const Real tolerance=TOLERANCE, const bool secure=true)
Definition: fe_interface.C:590
std::unique_ptr< FEMContext > coarse_context
FEGenericBase< Real > FEBase
Real weight(unsigned int var) const
Definition: system_norm.C:132

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

52 {
53  return KELLY;
54 }

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

◆ coarse_error

◆ 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

◆ fine_error

◆ 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


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