libMesh::DiscontinuityMeasure Class Reference

#include <discontinuity_measure.h>

Inheritance diagram for libMesh::DiscontinuityMeasure:

Public Types

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

Public Member Functions

 DiscontinuityMeasure ()
 
 ~DiscontinuityMeasure ()
 
void attach_essential_bc_function (std::pair< bool, Real > fptr(const System &system, const Point &p, const std::string &var_name))
 
virtual ErrorEstimatorType type () const libmesh_override
 
virtual void estimate_error (const System &system, ErrorVector &error_per_cell, const NumericVector< Number > *solution_vector=libmesh_nullptr, bool estimate_parent_error=false) libmesh_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=libmesh_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=libmesh_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) libmesh_override
 
virtual void internal_side_integration () libmesh_override
 
virtual bool boundary_side_integration () libmesh_override
 
void reinit_sides ()
 
float coarse_n_flux_faces_increment ()
 
void reduce_error (std::vector< ErrorVectorReal > &error_per_cell, const Parallel::Communicator &comm LIBMESH_CAN_DEFAULT_TO_COMMWORLD) const
 

Protected Attributes

std::pair< bool, Real >(* _bc_function )(const System &system, const Point &p, const std::string &var_name)
 
bool integrate_boundary_sides
 
UniquePtr< FEMContextfine_context
 
UniquePtr< FEMContextcoarse_context
 
Real fine_error
 
Real coarse_error
 
unsigned int var
 

Detailed Description

This class measures discontinuities between elements for debugging purposes. It derives from ErrorEstimator just in case someone finds it useful in a DG framework.

Author
Roy H. Stogner
Date
2006

Definition at line 49 of file discontinuity_measure.h.

Member Typedef Documentation

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 113 of file error_estimator.h.

Constructor & Destructor Documentation

libMesh::DiscontinuityMeasure::DiscontinuityMeasure ( )
inline

Constructor. Responsible for initializing the _bc_function function pointer to libmesh_nullptr. Defaults to L2 norm; changes to system norm are ignored.

Definition at line 58 of file discontinuity_measure.h.

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

58  :
61  { error_norm = L2; }
const class libmesh_nullptr_t libmesh_nullptr
std::pair< bool, Real >(* _bc_function)(const System &system, const Point &p, const std::string &var_name)
libMesh::DiscontinuityMeasure::~DiscontinuityMeasure ( )
inline

Destructor.

Definition at line 66 of file discontinuity_measure.h.

References attach_essential_bc_function().

66 {}

Member Function Documentation

void libMesh::DiscontinuityMeasure::attach_essential_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 essential BCs.

Definition at line 186 of file discontinuity_measure.C.

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

Referenced by ~DiscontinuityMeasure().

189 {
190  _bc_function = fptr;
191 
192  // We may be turning boundary side integration on or off
193  if (fptr)
195  else
196  integrate_boundary_sides = false;
197 }
std::pair< bool, Real >(* _bc_function)(const System &system, const Point &p, const std::string &var_name)
bool libMesh::DiscontinuityMeasure::boundary_side_integration ( )
protectedvirtual

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 115 of file discontinuity_measure.C.

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

Referenced by type().

116 {
117  const Elem & fine_elem = fine_context->get_elem();
118 
119  FEBase * fe_fine = libmesh_nullptr;
120  fine_context->get_side_fe( var, fe_fine, fine_elem.dim() );
121 
122  const std::string & var_name =
123  fine_context->get_system().variable_name(var);
124 
125  std::vector<std::vector<Real> > phi_fine = fe_fine->get_phi();
126  std::vector<Real> JxW_face = fe_fine->get_JxW();
127  std::vector<Point> qface_point = fe_fine->get_xyz();
128 
129  // The reinitialization also recomputes the locations of
130  // the quadrature points on the side. By checking if the
131  // first quadrature point on the side is on an essential boundary
132  // for a particular variable, we will determine if the whole
133  // element is on an essential boundary (assuming quadrature points
134  // are strictly contained in the side).
135  if (this->_bc_function(fine_context->get_system(),
136  qface_point[0], var_name).first)
137  {
138  const Real h = fine_elem.hmax();
139 
140  // The number of quadrature points
141  const unsigned int n_qp = fe_fine->n_quadrature_points();
142 
143  // The error contribution from this face
144  Real error = 1.e-30;
145 
146  // loop over the integration points on the face.
147  for (unsigned int qp=0; qp<n_qp; qp++)
148  {
149  // Value of the imposed essential BC at this quadrature point.
150  const std::pair<bool,Real> essential_bc =
151  this->_bc_function(fine_context->get_system(), qface_point[qp],
152  var_name);
153 
154  // Be sure the BC function still thinks we're on the
155  // essential boundary.
156  libmesh_assert_equal_to (essential_bc.first, true);
157 
158  // The solution value on each point
159  Number u_fine = fine_context->side_value(var, qp);
160 
161  // The difference between the desired BC and the approximate solution.
162  const Number jump = essential_bc.second - u_fine;
163 
164  // The flux jump squared. If using complex numbers,
165  // norm_sq(z) returns |z|^2, where |z| is the modulus of z.
166  const Real jump2 = TensorTools::norm_sq(jump);
167 
168  // Integrate the error on the face. The error is
169  // scaled by an additional power of h, where h is
170  // the maximum side length for the element. This
171  // arises in the definition of the indicator.
172  error += JxW_face[qp]*jump2;
173 
174  } // End quadrature point loop
175 
176  fine_error = error*h*error_norm.weight(var);
177 
178  return true;
179  } // end if side on flux boundary
180  return false;
181 }
Real weight(unsigned int var) const
Definition: system_norm.h:292
const class libmesh_nullptr_t libmesh_nullptr
UniquePtr< FEMContext > fine_context
std::pair< bool, Real >(* _bc_function)(const System &system, const Point &p, const std::string &var_name)
FEGenericBase< Real > FEBase
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
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 464 of file jump_error_estimator.C.

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

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

465 {
466  // Keep track of the number of internal flux sides found on each
467  // element
468  unsigned int dim = coarse_context->get_elem().dim();
469 
470  const unsigned int divisor =
471  1 << (dim-1)*(fine_context->get_elem().level() -
472  coarse_context->get_elem().level());
473 
474  // With a difference of n levels between fine and coarse elements,
475  // we compute a fractional flux face for the coarse element by adding:
476  // 1/2^n in 2D
477  // 1/4^n in 3D
478  // each time. This code will get hit 2^n times in 2D and 4^n
479  // times in 3D so that the final flux face count for the coarse
480  // element will be an integer value.
481 
482  return 1.0f / static_cast<float>(divisor);
483 }
UniquePtr< FEMContext > fine_context
UniquePtr< FEMContext > coarse_context
void libMesh::JumpErrorEstimator::estimate_error ( const System system,
ErrorVector error_per_cell,
const NumericVector< Number > *  solution_vector = libmesh_nullptr,
bool  estimate_parent_error = false 
)
virtualinherited

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::MeshBase::active_local_elements_begin(), libMesh::MeshBase::active_local_elements_end(), libMesh::JumpErrorEstimator::boundary_side_integration(), libMesh::Elem::child_ptr(), 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::FEAbstract::get_fe_type(), libMesh::System::get_mesh(), libMesh::FEAbstract::get_xyz(), libMesh::DofObject::id(), libMesh::JumpErrorEstimator::init_context(), libMesh::JumpErrorEstimator::integrate_boundary_sides, libMesh::JumpErrorEstimator::internal_side_integration(), libMesh::Elem::level(), libMesh::libmesh_ignore(), libmesh_nullptr, libMesh::MeshBase::max_elem_id(), mesh, libMesh::Elem::n_children(), libMesh::Elem::n_neighbors(), 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::DenseVector< T >::size(), libMesh::System::solution, libMesh::NumericVector< T >::swap(), libMesh::sys, libMesh::JumpErrorEstimator::var, and libMesh::SystemNorm::weight().

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

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  // Possibly skip this variable
150  if (error_norm.weight(var) == 0.0) continue;
151 
152  // FIXME: Need to generalize this to vector-valued elements. [PB]
153  FEBase * side_fe = libmesh_nullptr;
154 
155  const std::set<unsigned char> & elem_dims =
156  fine_context->elem_dimensions();
157 
158  for (std::set<unsigned char>::const_iterator dim_it =
159  elem_dims.begin(); dim_it != elem_dims.end(); ++dim_it)
160  {
161  const unsigned char dim = *dim_it;
162 
163  fine_context->get_side_fe( var, side_fe, dim );
164 
165  libmesh_assert_not_equal_to(side_fe->get_fe_type().family, SCALAR);
166 
167  side_fe->get_xyz();
168  }
169  }
170 
171  this->init_context(*fine_context);
173 
174  // Iterate over all the active elements in the mesh
175  // that live on this processor.
176  MeshBase::const_element_iterator elem_it = mesh.active_local_elements_begin();
177  const MeshBase::const_element_iterator elem_end = mesh.active_local_elements_end();
178 
179  for (; elem_it != elem_end; ++elem_it)
180  {
181  // e is necessarily an active element on the local processor
182  const Elem * e = *elem_it;
183  const dof_id_type e_id = e->id();
184 
185 #ifdef LIBMESH_ENABLE_AMR
186  // See if the parent of element e has been examined yet;
187  // if not, we may want to compute the estimator on it
188  const Elem * parent = e->parent();
189 
190  // We only can compute and only need to compute on
191  // parents with all active children
192  bool compute_on_parent = true;
193  if (!parent || !estimate_parent_error)
194  compute_on_parent = false;
195  else
196  for (unsigned int c=0; c != parent->n_children(); ++c)
197  if (!parent->child_ptr(c)->active())
198  compute_on_parent = false;
199 
200  if (compute_on_parent &&
201  !error_per_cell[parent->id()])
202  {
203  // Compute a projection onto the parent
204  DenseVector<Number> Uparent;
206  (*(system.solution), dof_map, parent, Uparent, false);
207 
208  // Loop over the neighbors of the parent
209  for (unsigned int n_p=0; n_p<parent->n_neighbors(); n_p++)
210  {
211  if (parent->neighbor_ptr(n_p) != libmesh_nullptr) // parent has a neighbor here
212  {
213  // Find the active neighbors in this direction
214  std::vector<const Elem *> active_neighbors;
215  parent->neighbor_ptr(n_p)->
216  active_family_tree_by_neighbor(active_neighbors,
217  parent);
218  // Compute the flux to each active neighbor
219  for (unsigned int a=0;
220  a != active_neighbors.size(); ++a)
221  {
222  const Elem * f = active_neighbors[a];
223  // FIXME - what about when f->level <
224  // parent->level()??
225  if (f->level() >= parent->level())
226  {
227  fine_context->pre_fe_reinit(system, f);
228  coarse_context->pre_fe_reinit(system, parent);
229  libmesh_assert_equal_to
230  (coarse_context->get_elem_solution().size(),
231  Uparent.size());
232  coarse_context->get_elem_solution() = Uparent;
233 
234  this->reinit_sides();
235 
236  // Loop over all significant variables in the system
237  for (var=0; var<n_vars; var++)
238  if (error_norm.weight(var) != 0.0)
239  {
241 
242  error_per_cell[fine_context->get_elem().id()] +=
243  static_cast<ErrorVectorReal>(fine_error);
244  error_per_cell[coarse_context->get_elem().id()] +=
245  static_cast<ErrorVectorReal>(coarse_error);
246  }
247 
248  // Keep track of the number of internal flux
249  // sides found on each element
251  {
252  n_flux_faces[fine_context->get_elem().id()]++;
253  n_flux_faces[coarse_context->get_elem().id()] +=
255  }
256  }
257  }
258  }
259  else if (integrate_boundary_sides)
260  {
261  fine_context->pre_fe_reinit(system, parent);
262  libmesh_assert_equal_to
263  (fine_context->get_elem_solution().size(),
264  Uparent.size());
265  fine_context->get_elem_solution() = Uparent;
266  fine_context->side = n_p;
267  fine_context->side_fe_reinit();
268 
269  // If we find a boundary flux for any variable,
270  // let's just count it as a flux face for all
271  // variables. Otherwise we'd need to keep track of
272  // a separate n_flux_faces and error_per_cell for
273  // every single var.
274  bool found_boundary_flux = false;
275 
276  for (var=0; var<n_vars; var++)
277  if (error_norm.weight(var) != 0.0)
278  {
279  if (this->boundary_side_integration())
280  {
281  error_per_cell[fine_context->get_elem().id()] +=
282  static_cast<ErrorVectorReal>(fine_error);
283  found_boundary_flux = true;
284  }
285  }
286 
287  if (scale_by_n_flux_faces && found_boundary_flux)
288  n_flux_faces[fine_context->get_elem().id()]++;
289  }
290  }
291  }
292 #endif // #ifdef LIBMESH_ENABLE_AMR
293 
294  // If we do any more flux integration, e will be the fine element
295  fine_context->pre_fe_reinit(system, e);
296 
297  // Loop over the neighbors of element e
298  for (unsigned int n_e=0; n_e<e->n_neighbors(); n_e++)
299  {
300  if ((e->neighbor_ptr(n_e) != libmesh_nullptr) ||
302  {
303  fine_context->side = n_e;
304  fine_context->side_fe_reinit();
305  }
306 
307  if (e->neighbor_ptr(n_e) != libmesh_nullptr) // e is not on the boundary
308  {
309  const Elem * f = e->neighbor_ptr(n_e);
310  const dof_id_type f_id = f->id();
311 
312  // Compute flux jumps if we are in case 1 or case 2.
313  if ((f->active() && (f->level() == e->level()) && (e_id < f_id))
314  || (f->level() < e->level()))
315  {
316  // f is now the coarse element
317  coarse_context->pre_fe_reinit(system, f);
318 
319  this->reinit_sides();
320 
321  // Loop over all significant variables in the system
322  for (var=0; var<n_vars; var++)
323  if (error_norm.weight(var) != 0.0)
324  {
326 
327  error_per_cell[fine_context->get_elem().id()] +=
328  static_cast<ErrorVectorReal>(fine_error);
329  error_per_cell[coarse_context->get_elem().id()] +=
330  static_cast<ErrorVectorReal>(coarse_error);
331  }
332 
333  // Keep track of the number of internal flux
334  // sides found on each element
336  {
337  n_flux_faces[fine_context->get_elem().id()]++;
338  n_flux_faces[coarse_context->get_elem().id()] +=
340  }
341  } // end if (case1 || case2)
342  } // if (e->neighbor(n_e) != libmesh_nullptr)
343 
344  // Otherwise, e is on the boundary. If it happens to
345  // be on a Dirichlet boundary, we need not do anything.
346  // On the other hand, if e is on a Neumann (flux) boundary
347  // with grad(u).n = g, we need to compute the additional residual
348  // (h * \int |g - grad(u_h).n|^2 dS)^(1/2).
349  // We can only do this with some knowledge of the boundary
350  // conditions, i.e. the user must have attached an appropriate
351  // BC function.
352  else if (integrate_boundary_sides)
353  {
354  bool found_boundary_flux = false;
355 
356  for (var=0; var<n_vars; var++)
357  if (error_norm.weight(var) != 0.0)
358  if (this->boundary_side_integration())
359  {
360  error_per_cell[fine_context->get_elem().id()] +=
361  static_cast<ErrorVectorReal>(fine_error);
362  found_boundary_flux = true;
363  }
364 
365  if (scale_by_n_flux_faces && found_boundary_flux)
366  n_flux_faces[fine_context->get_elem().id()]++;
367  } // end if (e->neighbor_ptr(n_e) == libmesh_nullptr)
368  } // end loop over neighbors
369  } // End loop over active local elements
370 
371 
372  // Each processor has now computed the error contributions
373  // for its local elements. We need to sum the vector
374  // and then take the square-root of each component. Note
375  // that we only need to sum if we are running on multiple
376  // processors, and we only need to take the square-root
377  // if the value is nonzero. There will in general be many
378  // zeros for the inactive elements.
379 
380  // First sum the vector of estimated error values
381  this->reduce_error(error_per_cell, system.comm());
382 
383  // Compute the square-root of each component.
384  for (std::size_t i=0; i<error_per_cell.size(); i++)
385  if (error_per_cell[i] != 0.)
386  error_per_cell[i] = std::sqrt(error_per_cell[i]);
387 
388 
389  if (this->scale_by_n_flux_faces)
390  {
391  // Sum the vector of flux face counts
392  this->reduce_error(n_flux_faces, system.comm());
393 
394  // Sanity check: Make sure the number of flux faces is
395  // always an integer value
396 #ifdef DEBUG
397  for (std::size_t i=0; i<n_flux_faces.size(); ++i)
398  libmesh_assert_equal_to (n_flux_faces[i], static_cast<float>(static_cast<unsigned int>(n_flux_faces[i])) );
399 #endif
400 
401  // Scale the error by the number of flux faces for each element
402  for (std::size_t i=0; i<n_flux_faces.size(); ++i)
403  {
404  if (n_flux_faces[i] == 0.0) // inactive or non-local element
405  continue;
406 
407  //libMesh::out << "Element " << i << " has " << n_flux_faces[i] << " flux faces." << std::endl;
408  error_per_cell[i] /= static_cast<ErrorVectorReal>(n_flux_faces[i]);
409  }
410  }
411 
412  // If we used a non-standard solution before, now is the time to fix
413  // the current_local_solution
414  if (solution_vector && solution_vector != system.solution.get())
415  {
416  NumericVector<Number> * newsol =
417  const_cast<NumericVector<Number> *>(solution_vector);
418  System & sys = const_cast<System &>(system);
419  newsol->swap(*sys.solution);
420  sys.update();
421  }
422 }
ImplicitSystem & sys
Real weight(unsigned int var) const
Definition: system_norm.h:292
MeshBase & mesh
const class libmesh_nullptr_t libmesh_nullptr
virtual void internal_side_integration()=0
UniquePtr< FEMContext > fine_context
const unsigned int n_vars
Definition: tecplot_io.C:68
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:804
FEGenericBase< Real > FEBase
void reduce_error(std::vector< ErrorVectorReal > &error_per_cell, const Parallel::Communicator &comm LIBMESH_CAN_DEFAULT_TO_COMMWORLD) const
void libmesh_ignore(const T &)
virtual void update()
Definition: system.C:423
virtual void get(const std::vector< numeric_index_type > &index, T *values) const
UniquePtr< FEMContext > coarse_context
virtual void init_context(FEMContext &c)
uint8_t dof_id_type
Definition: id_types.h:64
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 = libmesh_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 48 of file error_estimator.C.

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

Referenced by libMesh::ErrorEstimator::~ErrorEstimator().

53 {
54  SystemNorm old_error_norm = this->error_norm;
55 
56  // Sum the error values from each system
57  for (unsigned int s = 0; s != equation_systems.n_systems(); ++s)
58  {
59  ErrorVector system_error_per_cell;
60  const System & sys = equation_systems.get_system(s);
61  if (error_norms.find(&sys) == error_norms.end())
62  this->error_norm = old_error_norm;
63  else
64  this->error_norm = error_norms.find(&sys)->second;
65 
66  const NumericVector<Number> * solution_vector = libmesh_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 (std::size_t i=0; i != error_per_cell.size(); ++i)
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 }
ImplicitSystem & sys
const class libmesh_nullptr_t libmesh_nullptr
virtual void estimate_error(const System &system, ErrorVector &error_per_cell, const NumericVector< Number > *solution_vector=libmesh_nullptr, bool estimate_parent_error=false)=0
void libMesh::ErrorEstimator::estimate_errors ( const EquationSystems equation_systems,
ErrorMap errors_per_cell,
const std::map< const System *, const NumericVector< Number > * > *  solution_vectors = libmesh_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 94 of file error_estimator.C.

References libMesh::ErrorEstimator::error_norm, libMesh::ErrorEstimator::estimate_error(), libMesh::EquationSystems::get_system(), libmesh_nullptr, libMesh::EquationSystems::n_systems(), n_vars, libMesh::System::n_vars(), libMesh::sys, 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 (unsigned int s = 0; s != equation_systems.n_systems(); ++s)
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.find(std::make_pair(&sys, v)) ==
112  errors_per_cell.end())
113  continue;
114 
115  // Calculate error in only one variable
116  std::vector<Real> weights(n_vars, 0.0);
117  weights[v] = 1.0;
118  this->error_norm =
119  SystemNorm(std::vector<FEMNormType>(n_vars, old_error_norm.type(v)),
120  weights);
121 
122  const NumericVector<Number> * solution_vector = libmesh_nullptr;
123  if (solution_vectors &&
124  solution_vectors->find(&sys) != solution_vectors->end())
125  solution_vector = solution_vectors->find(&sys)->second;
126 
127  this->estimate_error
128  (sys, *errors_per_cell[std::make_pair(&sys, v)],
129  solution_vector, estimate_parent_error);
130  }
131  }
132 
133  // Restore our old state before returning
134  this->error_norm = old_error_norm;
135 }
ImplicitSystem & sys
const class libmesh_nullptr_t libmesh_nullptr
const unsigned int n_vars
Definition: tecplot_io.C:68
virtual void estimate_error(const System &system, ErrorVector &error_per_cell, const NumericVector< Number > *solution_vector=libmesh_nullptr, bool estimate_parent_error=false)=0
void libMesh::DiscontinuityMeasure::init_context ( FEMContext c)
protectedvirtual

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

Reimplemented from libMesh::JumpErrorEstimator.

Definition at line 42 of file discontinuity_measure.C.

References libMesh::FEMContext::elem_dimensions(), libMesh::ErrorEstimator::error_norm, libMesh::JumpErrorEstimator::fine_context, libMesh::FEGenericBase< OutputType >::get_phi(), libmesh_nullptr, n_vars, libMesh::DiffContext::n_vars(), and libMesh::SystemNorm::weight().

Referenced by type().

43 {
44  const unsigned int n_vars = c.n_vars();
45  for (unsigned int v=0; v<n_vars; v++)
46  {
47  // Possibly skip this variable
48  if (error_norm.weight(v) == 0.0) continue;
49 
50  // FIXME: Need to generalize this to vector-valued elements. [PB]
51  FEBase * side_fe = libmesh_nullptr;
52 
53  const std::set<unsigned char> & elem_dims =
54  c.elem_dimensions();
55 
56  for (std::set<unsigned char>::const_iterator dim_it =
57  elem_dims.begin(); dim_it != elem_dims.end(); ++dim_it)
58  {
59  const unsigned char dim = *dim_it;
60 
61  fine_context->get_side_fe( v, side_fe, dim );
62 
63  // We'll need values on both sides for discontinuity computation
64  side_fe->get_phi();
65  }
66  }
67 }
Real weight(unsigned int var) const
Definition: system_norm.h:292
const class libmesh_nullptr_t libmesh_nullptr
UniquePtr< FEMContext > fine_context
const unsigned int n_vars
Definition: tecplot_io.C:68
FEGenericBase< Real > FEBase
void libMesh::DiscontinuityMeasure::internal_side_integration ( )
protectedvirtual

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

Implements libMesh::JumpErrorEstimator.

Definition at line 72 of file discontinuity_measure.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::FEAbstract::get_JxW(), libMesh::FEGenericBase< OutputType >::get_phi(), libMesh::Elem::hmax(), libmesh_nullptr, libMesh::FEAbstract::n_quadrature_points(), libMesh::TensorTools::norm_sq(), libMesh::Real, libMesh::JumpErrorEstimator::var, and libMesh::SystemNorm::weight().

Referenced by type().

73 {
74  const Elem & coarse_elem = coarse_context->get_elem();
75  const Elem & fine_elem = fine_context->get_elem();
76 
77  FEBase * fe_fine = libmesh_nullptr;
78  fine_context->get_side_fe( var, fe_fine, fine_elem.dim() );
79 
80  FEBase * fe_coarse = libmesh_nullptr;
81  coarse_context->get_side_fe( var, fe_coarse, fine_elem.dim() );
82 
83  Real error = 1.e-30;
84  unsigned int n_qp = fe_fine->n_quadrature_points();
85 
86  std::vector<std::vector<Real> > phi_coarse = fe_coarse->get_phi();
87  std::vector<std::vector<Real> > phi_fine = fe_fine->get_phi();
88  std::vector<Real> JxW_face = fe_fine->get_JxW();
89 
90  for (unsigned int qp=0; qp != n_qp; ++qp)
91  {
92  // Calculate solution values on fine and coarse elements
93  // at this quadrature point
94  Number
95  u_fine = fine_context->side_value(var, qp),
96  u_coarse = coarse_context->side_value(var, qp);
97 
98  // Find the jump in the value
99  // at this quadrature point
100  const Number jump = u_fine - u_coarse;
101  const Real jump2 = TensorTools::norm_sq(jump);
102  // Accumulate the jump integral
103  error += JxW_face[qp] * jump2;
104  }
105 
106  // Add the h-weighted jump integral to each error term
107  fine_error =
108  error * fine_elem.hmax() * error_norm.weight(var);
109  coarse_error =
110  error * coarse_elem.hmax() * error_norm.weight(var);
111 }
Real weight(unsigned int var) const
Definition: system_norm.h:292
const class libmesh_nullptr_t libmesh_nullptr
UniquePtr< FEMContext > fine_context
FEGenericBase< Real > FEBase
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
UniquePtr< FEMContext > coarse_context
void libMesh::ErrorEstimator::reduce_error ( std::vector< ErrorVectorReal > &  error_per_cell,
const Parallel::Communicator &comm  LIBMESH_CAN_DEFAULT_TO_COMMWORLD 
) 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 33 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().

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

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

Definition at line 427 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(), libmesh_nullptr, n_vars, libMesh::FEAbstract::reinit(), and libMesh::SystemNorm::weight().

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

428 {
429  fine_context->side_fe_reinit();
430 
431  unsigned int dim = fine_context->get_elem().dim();
432  libmesh_assert_equal_to(dim, coarse_context->get_elem().dim());
433 
434  FEBase * fe_fine = libmesh_nullptr;
435  fine_context->get_side_fe( 0, fe_fine, dim );
436 
437  // Get the physical locations of the fine element quadrature points
438  std::vector<Point> qface_point = fe_fine->get_xyz();
439 
440  // Find the master quadrature point locations on the coarse element
441  FEBase * fe_coarse = libmesh_nullptr;
442  coarse_context->get_side_fe( 0, fe_coarse, dim );
443 
444  std::vector<Point> qp_coarse;
445 
447  (coarse_context->get_elem().dim(), fe_coarse->get_fe_type(),
448  &coarse_context->get_elem(), qface_point, qp_coarse);
449 
450  // The number of variables in the system
451  const unsigned int n_vars = fine_context->n_vars();
452 
453  // Calculate all coarse element shape functions at those locations
454  for (unsigned int v=0; v<n_vars; v++)
455  if (error_norm.weight(v) != 0.0)
456  {
457  coarse_context->get_side_fe( v, fe_coarse, dim );
458  fe_coarse->reinit (&coarse_context->get_elem(), &qp_coarse);
459  }
460 }
Real weight(unsigned int var) const
Definition: system_norm.h:292
const class libmesh_nullptr_t libmesh_nullptr
UniquePtr< FEMContext > fine_context
const unsigned int n_vars
Definition: tecplot_io.C:68
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:556
FEGenericBase< Real > FEBase
UniquePtr< FEMContext > coarse_context
virtual ErrorEstimatorType libMesh::DiscontinuityMeasure::type ( ) const
inlinevirtual

Member Data Documentation

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

Pointer to function that provides BC information.

Definition at line 103 of file discontinuity_measure.h.

Referenced by attach_essential_bc_function(), and boundary_side_integration().

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 150 of file error_estimator.h.

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

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

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

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

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


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