libMesh
jump_error_estimator.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 // C++ includes
20 #include <algorithm> // for std::fill
21 #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
22 #include <cmath> // for sqrt
23 
24 
25 // Local Includes
26 #include "libmesh/libmesh_common.h"
27 #include "libmesh/jump_error_estimator.h"
28 #include "libmesh/dof_map.h"
29 #include "libmesh/elem.h"
30 #include "libmesh/error_vector.h"
31 #include "libmesh/fe_base.h"
32 #include "libmesh/fe_interface.h"
33 #include "libmesh/fem_context.h"
34 #include "libmesh/libmesh_logging.h"
35 #include "libmesh/mesh_base.h"
36 #include "libmesh/quadrature_gauss.h"
37 #include "libmesh/system.h"
38 
39 #include "libmesh/dense_vector.h"
40 #include "libmesh/numeric_vector.h"
41 
42 namespace libMesh
43 {
44 
45 //-----------------------------------------------------------------
46 // JumpErrorEstimator implementations
48 {
49 }
50 
51 
52 
54  ErrorVector & error_per_cell,
55  const NumericVector<Number> * solution_vector,
56  bool estimate_parent_error)
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  for (const auto & e : mesh.active_local_element_ptr_range())
177  {
178  const dof_id_type e_id = e->id();
179 
180 #ifdef LIBMESH_ENABLE_AMR
181  // See if the parent of element e has been examined yet;
182  // if not, we may want to compute the estimator on it
183  const Elem * parent = e->parent();
184 
185  // We only can compute and only need to compute on
186  // parents with all active children
187  bool compute_on_parent = true;
188  if (!parent || !estimate_parent_error)
189  compute_on_parent = false;
190  else
191  for (auto & child : parent->child_ref_range())
192  if (!child.active())
193  compute_on_parent = false;
194 
195  if (compute_on_parent &&
196  !error_per_cell[parent->id()])
197  {
198  // Compute a projection onto the parent
199  DenseVector<Number> Uparent;
201  (*(system.solution), dof_map, parent, Uparent, false);
202 
203  // Loop over the neighbors of the parent
204  for (auto n_p : parent->side_index_range())
205  {
206  if (parent->neighbor_ptr(n_p) != libmesh_nullptr) // parent has a neighbor here
207  {
208  // Find the active neighbors in this direction
209  std::vector<const Elem *> active_neighbors;
210  parent->neighbor_ptr(n_p)->
211  active_family_tree_by_neighbor(active_neighbors,
212  parent);
213  // Compute the flux to each active neighbor
214  for (unsigned int a=0;
215  a != active_neighbors.size(); ++a)
216  {
217  const Elem * f = active_neighbors[a];
218  // FIXME - what about when f->level <
219  // parent->level()??
220  if (f->level() >= parent->level())
221  {
222  fine_context->pre_fe_reinit(system, f);
223  coarse_context->pre_fe_reinit(system, parent);
224  libmesh_assert_equal_to
225  (coarse_context->get_elem_solution().size(),
226  Uparent.size());
227  coarse_context->get_elem_solution() = Uparent;
228 
229  this->reinit_sides();
230 
231  // Loop over all significant variables in the system
232  for (var=0; var<n_vars; var++)
233  if (error_norm.weight(var) != 0.0)
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 = 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  {
274  if (this->boundary_side_integration())
275  {
276  error_per_cell[fine_context->get_elem().id()] +=
277  static_cast<ErrorVectorReal>(fine_error);
278  found_boundary_flux = true;
279  }
280  }
281 
282  if (scale_by_n_flux_faces && found_boundary_flux)
283  n_flux_faces[fine_context->get_elem().id()]++;
284  }
285  }
286  }
287 #endif // #ifdef LIBMESH_ENABLE_AMR
288 
289  // If we do any more flux integration, e will be the fine element
290  fine_context->pre_fe_reinit(system, e);
291 
292  // Loop over the neighbors of element e
293  for (auto n_e : e->side_index_range())
294  {
295  if ((e->neighbor_ptr(n_e) != libmesh_nullptr) ||
297  {
298  fine_context->side = n_e;
299  fine_context->side_fe_reinit();
300  }
301 
302  if (e->neighbor_ptr(n_e) != libmesh_nullptr) // e is not on the boundary
303  {
304  const Elem * f = e->neighbor_ptr(n_e);
305  const dof_id_type f_id = f->id();
306 
307  // Compute flux jumps if we are in case 1 or case 2.
308  if ((f->active() && (f->level() == e->level()) && (e_id < f_id))
309  || (f->level() < e->level()))
310  {
311  // f is now the coarse element
312  coarse_context->pre_fe_reinit(system, f);
313 
314  this->reinit_sides();
315 
316  // Loop over all significant variables in the system
317  for (var=0; var<n_vars; var++)
318  if (error_norm.weight(var) != 0.0)
319  {
321 
322  error_per_cell[fine_context->get_elem().id()] +=
323  static_cast<ErrorVectorReal>(fine_error);
324  error_per_cell[coarse_context->get_elem().id()] +=
325  static_cast<ErrorVectorReal>(coarse_error);
326  }
327 
328  // Keep track of the number of internal flux
329  // sides found on each element
331  {
332  n_flux_faces[fine_context->get_elem().id()]++;
333  n_flux_faces[coarse_context->get_elem().id()] +=
335  }
336  } // end if (case1 || case2)
337  } // if (e->neighbor(n_e) != libmesh_nullptr)
338 
339  // Otherwise, e is on the boundary. If it happens to
340  // be on a Dirichlet boundary, we need not do anything.
341  // On the other hand, if e is on a Neumann (flux) boundary
342  // with grad(u).n = g, we need to compute the additional residual
343  // (h * \int |g - grad(u_h).n|^2 dS)^(1/2).
344  // We can only do this with some knowledge of the boundary
345  // conditions, i.e. the user must have attached an appropriate
346  // BC function.
347  else if (integrate_boundary_sides)
348  {
349  bool found_boundary_flux = false;
350 
351  for (var=0; var<n_vars; var++)
352  if (error_norm.weight(var) != 0.0)
353  if (this->boundary_side_integration())
354  {
355  error_per_cell[fine_context->get_elem().id()] +=
356  static_cast<ErrorVectorReal>(fine_error);
357  found_boundary_flux = true;
358  }
359 
360  if (scale_by_n_flux_faces && found_boundary_flux)
361  n_flux_faces[fine_context->get_elem().id()]++;
362  } // end if (e->neighbor_ptr(n_e) == libmesh_nullptr)
363  } // end loop over neighbors
364  } // End loop over active local elements
365 
366 
367  // Each processor has now computed the error contributions
368  // for its local elements. We need to sum the vector
369  // and then take the square-root of each component. Note
370  // that we only need to sum if we are running on multiple
371  // processors, and we only need to take the square-root
372  // if the value is nonzero. There will in general be many
373  // zeros for the inactive elements.
374 
375  // First sum the vector of estimated error values
376  this->reduce_error(error_per_cell, system.comm());
377 
378  // Compute the square-root of each component.
379  for (std::size_t i=0; i<error_per_cell.size(); i++)
380  if (error_per_cell[i] != 0.)
381  error_per_cell[i] = std::sqrt(error_per_cell[i]);
382 
383 
384  if (this->scale_by_n_flux_faces)
385  {
386  // Sum the vector of flux face counts
387  this->reduce_error(n_flux_faces, system.comm());
388 
389  // Sanity check: Make sure the number of flux faces is
390  // always an integer value
391 #ifdef DEBUG
392  for (std::size_t i=0; i<n_flux_faces.size(); ++i)
393  libmesh_assert_equal_to (n_flux_faces[i], static_cast<float>(static_cast<unsigned int>(n_flux_faces[i])) );
394 #endif
395 
396  // Scale the error by the number of flux faces for each element
397  for (std::size_t i=0; i<n_flux_faces.size(); ++i)
398  {
399  if (n_flux_faces[i] == 0.0) // inactive or non-local element
400  continue;
401 
402  //libMesh::out << "Element " << i << " has " << n_flux_faces[i] << " flux faces." << std::endl;
403  error_per_cell[i] /= static_cast<ErrorVectorReal>(n_flux_faces[i]);
404  }
405  }
406 
407  // If we used a non-standard solution before, now is the time to fix
408  // the current_local_solution
409  if (solution_vector && solution_vector != system.solution.get())
410  {
411  NumericVector<Number> * newsol =
412  const_cast<NumericVector<Number> *>(solution_vector);
413  System & sys = const_cast<System &>(system);
414  newsol->swap(*sys.solution);
415  sys.update();
416  }
417 }
418 
419 
420 
421 void
423 {
424  fine_context->side_fe_reinit();
425 
426  unsigned int dim = fine_context->get_elem().dim();
427  libmesh_assert_equal_to(dim, coarse_context->get_elem().dim());
428 
429  FEBase * fe_fine = libmesh_nullptr;
430  fine_context->get_side_fe( 0, fe_fine, dim );
431 
432  // Get the physical locations of the fine element quadrature points
433  std::vector<Point> qface_point = fe_fine->get_xyz();
434 
435  // Find the master quadrature point locations on the coarse element
436  FEBase * fe_coarse = libmesh_nullptr;
437  coarse_context->get_side_fe( 0, fe_coarse, dim );
438 
439  std::vector<Point> qp_coarse;
440 
442  (coarse_context->get_elem().dim(), fe_coarse->get_fe_type(),
443  &coarse_context->get_elem(), qface_point, qp_coarse);
444 
445  // The number of variables in the system
446  const unsigned int n_vars = fine_context->n_vars();
447 
448  // Calculate all coarse element shape functions at those locations
449  for (unsigned int v=0; v<n_vars; v++)
450  if (error_norm.weight(v) != 0.0)
451  {
452  coarse_context->get_side_fe( v, fe_coarse, dim );
453  fe_coarse->reinit (&coarse_context->get_elem(), &qp_coarse);
454  }
455 }
456 
457 
458 
460 {
461  // Keep track of the number of internal flux sides found on each
462  // element
463  unsigned int dim = coarse_context->get_elem().dim();
464 
465  const unsigned int divisor =
466  1 << (dim-1)*(fine_context->get_elem().level() -
467  coarse_context->get_elem().level());
468 
469  // With a difference of n levels between fine and coarse elements,
470  // we compute a fractional flux face for the coarse element by adding:
471  // 1/2^n in 2D
472  // 1/4^n in 3D
473  // each time. This code will get hit 2^n times in 2D and 4^n
474  // times in 3D so that the final flux face count for the coarse
475  // element will be an integer value.
476 
477  return 1.0f / static_cast<float>(divisor);
478 }
479 
480 } // namespace libMesh
FEFamily family
The type of finite element.
Definition: fe_type.h:203
bool active() const
Definition: elem.h:2257
SystemNorm error_norm
When estimating the error in a single system, the error_norm is used to control the scaling and norm ...
unsigned int dim
The ErrorVector is a specialization of the StatisticsVector for error data computed on a finite eleme...
Definition: error_vector.h:50
ImplicitSystem & sys
IntRange< unsigned short > side_index_range() const
Definition: elem.h:2083
Real weight(unsigned int var) const
Definition: system_norm.h:292
const Elem * parent() const
Definition: elem.h:2346
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
MeshBase & mesh
const class libmesh_nullptr_t libmesh_nullptr
virtual void internal_side_integration()=0
The function, to be implemented by derived classes, which calculates an error term on an internal sid...
UniquePtr< FEMContext > fine_context
Context objects for integrating on the fine and coarse elements sharing a face.
const unsigned int n_vars
Definition: tecplot_io.C:68
The libMesh namespace provides an interface to certain functionality in the library.
virtual unsigned int size() const libmesh_override
Definition: dense_vector.h:87
FEType get_fe_type() const
Definition: fe_abstract.h:421
DIE A HORRIBLE DEATH HERE typedef float ErrorVectorReal
This is the MeshBase class.
Definition: mesh_base.h:68
SimpleRange< ChildRefIter > child_ref_range()
Returns a range with all children of a parent element, usable in range-based for loops.
Definition: elem.h:1699
const Elem * neighbor_ptr(unsigned int i) const
Definition: elem.h:1967
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:802
virtual SimpleRange< element_iterator > active_local_element_ptr_range()=0
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:167
void reinit_sides()
A utility function to reinit the finite element data on elements sharing a side.
virtual dof_id_type max_elem_id() const =0
const MeshBase & get_mesh() const
Definition: system.h:2014
virtual bool boundary_side_integration()
The function, to be implemented by derived classes, which calculates an error term on a boundary side...
const DofMap & get_dof_map() const
Definition: system.h:2030
This is the base class for classes which contain information related to any physical process that mig...
Definition: system.h:76
This class provides all data required for a physics package (e.g.
Definition: fem_context.h:61
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:569
UniquePtr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1523
void reduce_error(std::vector< ErrorVectorReal > &error_per_cell, const Parallel::Communicator &comm LIBMESH_CAN_DEFAULT_TO_COMMWORLD) const
This method takes the local error contributions in error_per_cell from each processor and combines th...
Real fine_error
The fine and coarse error values to be set by each side_integration();.
void libmesh_ignore(const T &)
virtual void swap(NumericVector< T > &v)
Swaps the contents of this with v.
const Parallel::Communicator & comm() const
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
This function uses the derived class&#39;s jump error estimate formula to estimate the error on each cell...
float coarse_n_flux_faces_increment()
A utility function to correctly increase n_flux_faces for the coarse element.
unsigned int level() const
Definition: elem.h:2388
unsigned int var
The variable number currently being evaluated.
unsigned int n_vars() const
Definition: system.h:2086
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...
dof_id_type id() const
Definition: dof_object.h:632
bool integrate_boundary_sides
A boolean flag, by default false, to be set to true if integrations with boundary_side_integration() ...
UniquePtr< FEMContext > coarse_context
const std::vector< Point > & get_xyz() const
Definition: fe_abstract.h:230
virtual void init_context(FEMContext &c)
An initialization function, to give derived classes a chance to request specific data from the FE obj...
This class forms the foundation from which generic finite elements may be derived.
virtual void reinit(const Elem *elem, const std::vector< Point > *const pts=libmesh_nullptr, const std::vector< Real > *const weights=libmesh_nullptr)=0
This is at the core of this class.
uint8_t dof_id_type
Definition: id_types.h:64