jump_error_estimator.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2018 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"
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"
35 #include "libmesh/mesh_base.h"
37 #include "libmesh/system.h"
38 #include "libmesh/dense_vector.h"
39 #include "libmesh/numeric_vector.h"
40 #include "libmesh/int_range.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  // 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 }
420 
421 
422 
423 void
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 }
459 
460 
461 
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 }
482 
483 } // namespace libMesh
virtual unsigned int size() const override
Definition: dense_vector.h:92
FEFamily family
Definition: fe_type.h:204
const Elem * parent() const
Definition: elem.h:2479
const std::vector< Point > & get_xyz() const
Definition: fe_abstract.h:238
IntRange< unsigned short > side_index_range() const
Definition: elem.h:2166
std::unique_ptr< FEMContext > fine_context
The base class for all geometric element types.
Definition: elem.h:100
MeshBase & mesh
IntRange< std::size_t > index_range(const std::vector< T > &vec)
Definition: int_range.h:104
const Parallel::Communicator & comm() const
virtual void internal_side_integration()=0
const unsigned int n_vars
Definition: tecplot_io.C:69
const MeshBase & get_mesh() const
Definition: system.h:2033
DIE A HORRIBLE DEATH HERE typedef float ErrorVectorReal
Base class for Mesh.
Definition: mesh_base.h:77
SimpleRange< ChildRefIter > child_ref_range()
Definition: elem.h:1779
FEType get_fe_type() const
Definition: fe_abstract.h:429
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
Manages the degrees of freedom (DOFs) in a simulation.
Definition: dof_map.h:176
void libmesh_ignore(const Args &...)
dof_id_type id() const
Definition: dof_object.h:655
Manages consistently variables, degrees of freedom, and coefficient vectors.
Definition: system.h:92
std::unique_ptr< NumericVector< Number > > solution
Definition: system.h:1523
virtual void reinit(const Elem *elem, const std::vector< Point > *const pts=nullptr, const std::vector< Real > *const weights=nullptr)=0
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
void reduce_error(std::vector< ErrorVectorReal > &error_per_cell, const Parallel::Communicator &comm) const
std::unique_ptr< FEMContext > coarse_context
const Elem * neighbor_ptr(unsigned int i) const
Definition: elem.h:2050
Real weight(unsigned int var) const
Definition: system_norm.C:132
unsigned int level() const
Definition: elem.h:2521
const FEType & variable_type(const unsigned int i) const
Definition: system.h:2183
virtual void swap(NumericVector< T > &v)
virtual void estimate_error(const System &system, ErrorVector &error_per_cell, const NumericVector< Number > *solution_vector=nullptr, bool estimate_parent_error=false) override
unsigned int n_vars() const
Definition: system.h:2105
bool active() const
Definition: elem.h:2390
const DofMap & get_dof_map() const
Definition: system.h:2049
virtual void init_context(FEMContext &c)
uint8_t dof_id_type
Definition: id_types.h:64