libMesh
Public Types | Public Member Functions | Public Attributes | Protected Member Functions | List of all members
libMesh::UniformRefinementEstimator Class Reference

This class implements a ``brute force'' error estimator which integrates differences between the current solution and the solution on a uniformly refined (in h and/or p, for an arbitrary number of levels) grid. More...

#include <uniform_refinement_estimator.h>

Inheritance diagram for libMesh::UniformRefinementEstimator:
[legend]

Public Types

typedef std::map< std::pair< const System *, unsigned int >, ErrorVector * > ErrorMap
 When calculating many error vectors at once, we need a data structure to hold them all. More...
 

Public Member Functions

 UniformRefinementEstimator ()
 Constructor. More...
 
 ~UniformRefinementEstimator ()
 Destructor. More...
 
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 does uniform refinements and a solve to get an improved solution on each cell, then estimates the error by integrating differences between the coarse and fine solutions. More...
 
virtual void estimate_errors (const EquationSystems &equation_systems, ErrorVector &error_per_cell, const std::map< const System *, SystemNorm > &error_norms, const std::map< const System *, const NumericVector< Number > * > *solution_vectors=libmesh_nullptr, bool estimate_parent_error=false) libmesh_override
 Currently this function ignores the error_norm member variable, and uses the function argument error_norms instead. More...
 
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) libmesh_override
 Currently this function ignores the component_scale member variable, because it calculates each error individually, unscaled. More...
 
virtual ErrorEstimatorType type () const libmesh_override
 

Public Attributes

unsigned char number_h_refinements
 How many h refinements to perform to get the fine grid. More...
 
unsigned char number_p_refinements
 How many p refinements to perform to get the fine grid. More...
 
SystemNorm error_norm
 When estimating the error in a single system, the error_norm is used to control the scaling and norm choice for each variable. More...
 

Protected Member Functions

virtual void _estimate_error (const EquationSystems *equation_systems, const System *system, ErrorVector *error_per_cell, std::map< std::pair< const System *, unsigned int >, ErrorVector * > *errors_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)
 The code for estimate_error and both estimate_errors versions is very similar, so we use the same function for all three. More...
 
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 them to get the global error vector. More...
 

Detailed Description

This class implements a ``brute force'' error estimator which integrates differences between the current solution and the solution on a uniformly refined (in h and/or p, for an arbitrary number of levels) grid.

Author
Roy H. Stogner
Date
2006

Definition at line 45 of file uniform_refinement_estimator.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::UniformRefinementEstimator::UniformRefinementEstimator ( )

Constructor.

Sets the most common default parameter values.

Definition at line 52 of file uniform_refinement_estimator.h.

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

52  :
56  { error_norm = H1; }
unsigned char number_h_refinements
How many h refinements to perform to get the fine grid.
SystemNorm error_norm
When estimating the error in a single system, the error_norm is used to control the scaling and norm ...
ErrorEstimator()
Constructor.
unsigned char number_p_refinements
How many p refinements to perform to get the fine grid.
libMesh::UniformRefinementEstimator::~UniformRefinementEstimator ( )

Destructor.

Definition at line 61 of file uniform_refinement_estimator.h.

References estimate_error(), estimate_errors(), and libmesh_nullptr.

61 {}

Member Function Documentation

void libMesh::UniformRefinementEstimator::_estimate_error ( const EquationSystems equation_systems,
const System system,
ErrorVector error_per_cell,
std::map< std::pair< const System *, unsigned int >, ErrorVector * > *  errors_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 
)
protectedvirtual

The code for estimate_error and both estimate_errors versions is very similar, so we use the same function for all three.

Definition at line 97 of file uniform_refinement_estimator.C.

References libMesh::MeshBase::active_local_element_ptr_range(), libMesh::EquationSystems::adjoint_solve(), libMesh::System::adjoint_solve(), libMesh::FEGenericBase< OutputType >::build(), libMesh::NumericVector< T >::build(), libMesh::NumericVector< T >::clear(), libMesh::ParallelObject::comm(), libMesh::System::current_local_solution, libMesh::System::current_solution(), libMesh::FEType::default_quadrature_rule(), dim, libMesh::System::disable_cache(), libMesh::DofMap::dof_indices(), libMesh::ErrorEstimator::error_norm, libMesh::ErrorVectorReal, libMesh::System::get_adjoint_solution(), libMesh::System::get_dof_map(), libMesh::System::get_equation_systems(), libMesh::EquationSystems::get_mesh(), libMesh::DofMap::get_send_list(), libMesh::EquationSystems::get_system(), libMesh::System::get_vector(), libMesh::H1, libMesh::H1_SEMINORM, libMesh::H2, libMesh::H2_SEMINORM, libMesh::DofObject::id(), libMesh::invalid_uint, libMesh::L2, libMesh::libmesh_assert(), libMesh::MeshBase::max_elem_id(), mesh, libMesh::MeshBase::mesh_dimension(), libMesh::MeshBase::n_elem(), libMesh::EquationSystems::n_systems(), n_vars, libMesh::System::n_vars(), libMesh::TensorTools::norm_sq(), libMesh::TypeVector< T >::norm_sq(), libMesh::TypeTensor< T >::norm_sq(), number_h_refinements, number_p_refinements, libMesh::Elem::parent(), libMesh::MeshBase::partitioner(), libMesh::System::project_solution_on_reinit(), libMesh::System::qoi, libMesh::Real, libMesh::ErrorEstimator::reduce_error(), libMesh::EquationSystems::reinit(), libMesh::System::request_vector(), libMesh::SERIAL, libMesh::System::solution, libMesh::EquationSystems::solve(), libMesh::System::solve(), libMesh::NumericVector< T >::swap(), libMesh::sys, libMesh::SystemNorm::type(), libMesh::MeshRefinement::uniformly_coarsen(), libMesh::MeshRefinement::uniformly_p_coarsen(), libMesh::MeshRefinement::uniformly_p_refine(), libMesh::MeshRefinement::uniformly_refine(), libMesh::System::update(), libMesh::DofMap::variable_type(), libMesh::System::vectors_begin(), libMesh::System::vectors_end(), and libMesh::SystemNorm::weight_sq().

Referenced by estimate_error(), and estimate_errors().

104 {
105  // Get a vector of the Systems we're going to work on,
106  // and set up a error_norms map if necessary
107  std::vector<System *> system_list;
108  UniquePtr<std::map<const System *, SystemNorm>> error_norms =
109  UniquePtr<std::map<const System *, SystemNorm>>
110  (new std::map<const System *, SystemNorm>);
111 
112  if (_es)
113  {
114  libmesh_assert(!_system);
115  libmesh_assert(_es->n_systems());
116  _system = &(_es->get_system(0));
117  libmesh_assert_equal_to (&(_system->get_equation_systems()), _es);
118 
119  libmesh_assert(_es->n_systems());
120  for (unsigned int i=0; i != _es->n_systems(); ++i)
121  // We have to break the rules here, because we can't refine a const System
122  system_list.push_back(const_cast<System *>(&(_es->get_system(i))));
123 
124  // If we're computing one vector, we need to know how to scale
125  // each variable's contributions to it.
126  if (_error_norms)
127  {
128  libmesh_assert(!errors_per_cell);
129  }
130  else
131  // If we're computing many vectors, we just need to know which
132  // variables to skip
133  {
134  libmesh_assert (errors_per_cell);
135 
136  _error_norms = error_norms.get();
137 
138  for (unsigned int i=0; i!= _es->n_systems(); ++i)
139  {
140  const System & sys = _es->get_system(i);
141  unsigned int n_vars = sys.n_vars();
142 
143  std::vector<Real> weights(n_vars, 0.0);
144  for (unsigned int v = 0; v != n_vars; ++v)
145  {
146  if (errors_per_cell->find(std::make_pair(&sys, v)) ==
147  errors_per_cell->end())
148  continue;
149 
150  weights[v] = 1.0;
151  }
152  (*error_norms)[&sys] =
153  SystemNorm(std::vector<FEMNormType>(n_vars, error_norm.type(0)),
154  weights);
155  }
156  }
157  }
158  else
159  {
160  libmesh_assert(_system);
161  // We have to break the rules here, because we can't refine a const System
162  system_list.push_back(const_cast<System *>(_system));
163 
164  libmesh_assert(!_error_norms);
165  (*error_norms)[_system] = error_norm;
166  _error_norms = error_norms.get();
167  }
168 
169  // An EquationSystems reference will be convenient.
170  // We have to break the rules here, because we can't refine a const System
171  EquationSystems & es =
172  const_cast<EquationSystems &>(_system->get_equation_systems());
173 
174  // The current mesh
175  MeshBase & mesh = es.get_mesh();
176 
177  // The dimensionality of the mesh
178  const unsigned int dim = mesh.mesh_dimension();
179 
180  // Resize the error_per_cell vectors to be
181  // the number of elements, initialize them to 0.
182  if (error_per_cell)
183  {
184  error_per_cell->clear();
185  error_per_cell->resize (mesh.max_elem_id(), 0.);
186  }
187  else
188  {
189  libmesh_assert(errors_per_cell);
190  for (ErrorMap::iterator i = errors_per_cell->begin();
191  i != errors_per_cell->end(); ++i)
192  {
193  ErrorVector * e = i->second;
194  e->clear();
195  e->resize(mesh.max_elem_id(), 0.);
196  }
197  }
198 
199  // We'll want to back up all coarse grid vectors
200  std::vector<std::map<std::string, NumericVector<Number> *>>
201  coarse_vectors(system_list.size());
202  std::vector<NumericVector<Number> *>
203  coarse_solutions(system_list.size());
204  std::vector<NumericVector<Number> *>
205  coarse_local_solutions(system_list.size());
206  // And make copies of projected solutions
207  std::vector<NumericVector<Number> *>
208  projected_solutions(system_list.size());
209 
210  // And we'll need to temporarily change solution projection settings
211  std::vector<bool> old_projection_settings(system_list.size());
212 
213  // And it'll be best to avoid any repartitioning
214  UniquePtr<Partitioner> old_partitioner(mesh.partitioner().release());
215 
216  for (std::size_t i=0; i != system_list.size(); ++i)
217  {
218  System & system = *system_list[i];
219 
220  // Check for valid error_norms
221  libmesh_assert (_error_norms->find(&system) !=
222  _error_norms->end());
223 
224  // Back up the solution vector
225  coarse_solutions[i] = system.solution->clone().release();
226  coarse_local_solutions[i] =
227  system.current_local_solution->clone().release();
228 
229  // Back up all other coarse grid vectors
230  for (System::vectors_iterator vec = system.vectors_begin(); vec !=
231  system.vectors_end(); ++vec)
232  {
233  // The (string) name of this vector
234  const std::string & var_name = vec->first;
235 
236  coarse_vectors[i][var_name] = vec->second->clone().release();
237  }
238 
239  // Use a non-standard solution vector if necessary
240  if (solution_vectors &&
241  solution_vectors->find(&system) != solution_vectors->end() &&
242  solution_vectors->find(&system)->second &&
243  solution_vectors->find(&system)->second != system.solution.get())
244  {
245  NumericVector<Number> * newsol =
246  const_cast<NumericVector<Number> *>
247  (solution_vectors->find(&system)->second);
248  newsol->swap(*system.solution);
249  system.update();
250  }
251 
252  // Make sure the solution is projected when we refine the mesh
253  old_projection_settings[i] = system.project_solution_on_reinit();
254  system.project_solution_on_reinit() = true;
255  }
256 
257  // Find the number of coarse mesh elements, to make it possible
258  // to find correct coarse elem ids later
259  const dof_id_type max_coarse_elem_id = mesh.max_elem_id();
260 #ifndef NDEBUG
261  // n_coarse_elem is only used in an assertion later so
262  // avoid declaring it unless asserts are active.
263  const dof_id_type n_coarse_elem = mesh.n_elem();
264 #endif
265 
266  // Uniformly refine the mesh
267  MeshRefinement mesh_refinement(mesh);
268 
270 
271  // FIXME: this may break if there is more than one System
272  // on this mesh but estimate_error was still called instead of
273  // estimate_errors
274  for (unsigned int i = 0; i != number_h_refinements; ++i)
275  {
276  mesh_refinement.uniformly_refine(1);
277  es.reinit();
278  }
279 
280  for (unsigned int i = 0; i != number_p_refinements; ++i)
281  {
282  mesh_refinement.uniformly_p_refine(1);
283  es.reinit();
284  }
285 
286  for (std::size_t i=0; i != system_list.size(); ++i)
287  {
288  System & system = *system_list[i];
289 
290  // Copy the projected coarse grid solutions, which will be
291  // overwritten by solve()
292  // projected_solutions[i] = system.solution->clone().release();
293  projected_solutions[i] = NumericVector<Number>::build(system.comm()).release();
294  projected_solutions[i]->init(system.solution->size(), true, SERIAL);
295  system.solution->localize(*projected_solutions[i],
296  system.get_dof_map().get_send_list());
297  }
298 
299  // Are we doing a forward or an adjoint solve?
300  bool solve_adjoint = false;
301  if (solution_vectors)
302  {
303  System * sys = system_list[0];
304  libmesh_assert (solution_vectors->find(sys) !=
305  solution_vectors->end());
306  const NumericVector<Number> * vec = solution_vectors->find(sys)->second;
307  for (std::size_t j=0; j != sys->qoi.size(); ++j)
308  {
309  std::ostringstream adjoint_name;
310  adjoint_name << "adjoint_solution" << j;
311 
312  if (vec == sys->request_vector(adjoint_name.str()))
313  {
314  solve_adjoint = true;
315  break;
316  }
317  }
318  }
319 
320  // Get the uniformly refined solution.
321 
322  if (_es)
323  {
324  // Even if we had a decent preconditioner, valid matrix etc. before
325  // refinement, we don't any more.
326  for (unsigned int i=0; i != es.n_systems(); ++i)
327  es.get_system(i).disable_cache();
328 
329  // No specified vectors == forward solve
330  if (!solution_vectors)
331  es.solve();
332  else
333  {
334  libmesh_assert_equal_to (solution_vectors->size(), es.n_systems());
335  libmesh_assert (solution_vectors->find(system_list[0]) !=
336  solution_vectors->end());
337  libmesh_assert(solve_adjoint ||
338  (solution_vectors->find(system_list[0])->second ==
339  system_list[0]->solution.get()) ||
340  !solution_vectors->find(system_list[0])->second);
341 
342 #ifdef DEBUG
343  for (std::size_t i=0; i != system_list.size(); ++i)
344  {
345  System * sys = system_list[i];
346  libmesh_assert (solution_vectors->find(sys) !=
347  solution_vectors->end());
348  const NumericVector<Number> * vec = solution_vectors->find(sys)->second;
349  if (solve_adjoint)
350  {
351  bool found_vec = false;
352  for (std::size_t j=0; j != sys->qoi.size(); ++j)
353  {
354  std::ostringstream adjoint_name;
355  adjoint_name << "adjoint_solution" << j;
356 
357  if (vec == sys->request_vector(adjoint_name.str()))
358  {
359  found_vec = true;
360  break;
361  }
362  }
363  libmesh_assert(found_vec);
364  }
365  else
366  libmesh_assert(vec == sys->solution.get() || !vec);
367  }
368 #endif
369 
370  if (solve_adjoint)
371  {
372  std::vector<unsigned int> adjs(system_list.size(),
374  // Set up proper initial guesses
375  for (std::size_t i=0; i != system_list.size(); ++i)
376  {
377  System * sys = system_list[i];
378  libmesh_assert (solution_vectors->find(sys) !=
379  solution_vectors->end());
380  const NumericVector<Number> * vec = solution_vectors->find(sys)->second;
381  for (std::size_t j=0; j != sys->qoi.size(); ++j)
382  {
383  std::ostringstream adjoint_name;
384  adjoint_name << "adjoint_solution" << j;
385 
386  if (vec == sys->request_vector(adjoint_name.str()))
387  {
388  adjs[i] = j;
389  break;
390  }
391  }
392  libmesh_assert_not_equal_to (adjs[i], libMesh::invalid_uint);
393  system_list[i]->get_adjoint_solution(adjs[i]) =
394  *system_list[i]->solution;
395  }
396 
397  es.adjoint_solve();
398 
399  // Put the adjoint_solution into solution for
400  // comparisons
401  for (std::size_t i=0; i != system_list.size(); ++i)
402  {
403  system_list[i]->get_adjoint_solution(adjs[i]).swap(*system_list[i]->solution);
404  system_list[i]->update();
405  }
406  }
407  else
408  es.solve();
409  }
410  }
411  else
412  {
413  System * sys = system_list[0];
414 
415  // Even if we had a decent preconditioner, valid matrix etc. before
416  // refinement, we don't any more.
417  sys->disable_cache();
418 
419  // No specified vectors == forward solve
420  if (!solution_vectors)
421  sys->solve();
422  else
423  {
424  libmesh_assert (solution_vectors->find(sys) !=
425  solution_vectors->end());
426 
427  const NumericVector<Number> * vec = solution_vectors->find(sys)->second;
428 
429  libmesh_assert(solve_adjoint ||
430  (solution_vectors->find(sys)->second ==
431  sys->solution.get()) ||
432  !solution_vectors->find(sys)->second);
433 
434  if (solve_adjoint)
435  {
436  unsigned int adj = libMesh::invalid_uint;
437  for (std::size_t j=0; j != sys->qoi.size(); ++j)
438  {
439  std::ostringstream adjoint_name;
440  adjoint_name << "adjoint_solution" << j;
441 
442  if (vec == sys->request_vector(adjoint_name.str()))
443  {
444  adj = j;
445  break;
446  }
447  }
448  libmesh_assert_not_equal_to (adj, libMesh::invalid_uint);
449 
450  // Set up proper initial guess
451  sys->get_adjoint_solution(adj) = *sys->solution;
452  sys->adjoint_solve();
453  // Put the adjoint_solution into solution for
454  // comparisons
455  sys->get_adjoint_solution(adj).swap(*sys->solution);
456  sys->update();
457  }
458  else
459  sys->solve();
460  }
461  }
462 
463  // Get the error in the uniformly refined solution(s).
464 
465  for (std::size_t sysnum=0; sysnum != system_list.size(); ++sysnum)
466  {
467  System & system = *system_list[sysnum];
468 
469  unsigned int n_vars = system.n_vars();
470 
471  DofMap & dof_map = system.get_dof_map();
472 
473  const SystemNorm & system_i_norm =
474  _error_norms->find(&system)->second;
475 
476  NumericVector<Number> * projected_solution = projected_solutions[sysnum];
477 
478  // Loop over all the variables in the system
479  for (unsigned int var=0; var<n_vars; var++)
480  {
481  // Get the error vector to fill for this system and variable
482  ErrorVector * err_vec = error_per_cell;
483  if (!err_vec)
484  {
485  libmesh_assert(errors_per_cell);
486  err_vec =
487  (*errors_per_cell)[std::make_pair(&system,var)];
488  }
489 
490  // The type of finite element to use for this variable
491  const FEType & fe_type = dof_map.variable_type (var);
492 
493  // Finite element object for each fine element
494  UniquePtr<FEBase> fe (FEBase::build (dim, fe_type));
495 
496  // Build and attach an appropriate quadrature rule
497  UniquePtr<QBase> qrule = fe_type.default_quadrature_rule(dim);
498  fe->attach_quadrature_rule (qrule.get());
499 
500  const std::vector<Real> & JxW = fe->get_JxW();
501  const std::vector<std::vector<Real>> & phi = fe->get_phi();
502  const std::vector<std::vector<RealGradient>> & dphi =
503  fe->get_dphi();
504 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
505  const std::vector<std::vector<RealTensor>> & d2phi =
506  fe->get_d2phi();
507 #endif
508 
509  // The global DOF indices for the fine element
510  std::vector<dof_id_type> dof_indices;
511 
512  // Iterate over all the active elements in the fine mesh
513  // that live on this processor.
514  for (const auto & elem : mesh.active_local_element_ptr_range())
515  {
516  // Find the element id for the corresponding coarse grid element
517  const Elem * coarse = elem;
518  dof_id_type e_id = coarse->id();
519  while (e_id >= max_coarse_elem_id)
520  {
521  libmesh_assert (coarse->parent());
522  coarse = coarse->parent();
523  e_id = coarse->id();
524  }
525 
526  Real L2normsq = 0., H1seminormsq = 0.;
527 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
528  Real H2seminormsq = 0.;
529 #endif
530 
531  // reinitialize the element-specific data
532  // for the current element
533  fe->reinit (elem);
534 
535  // Get the local to global degree of freedom maps
536  dof_map.dof_indices (elem, dof_indices, var);
537 
538  // The number of quadrature points
539  const unsigned int n_qp = qrule->n_points();
540 
541  // The number of shape functions
542  const unsigned int n_sf =
543  cast_int<unsigned int>(dof_indices.size());
544 
545  //
546  // Begin the loop over the Quadrature points.
547  //
548  for (unsigned int qp=0; qp<n_qp; qp++)
549  {
550  Number u_fine = 0., u_coarse = 0.;
551 
552  Gradient grad_u_fine, grad_u_coarse;
553 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
554  Tensor grad2_u_fine, grad2_u_coarse;
555 #endif
556 
557  // Compute solution values at the current
558  // quadrature point. This requires a sum
559  // over all the shape functions evaluated
560  // at the quadrature point.
561  for (unsigned int i=0; i<n_sf; i++)
562  {
563  u_fine += phi[i][qp]*system.current_solution (dof_indices[i]);
564  u_coarse += phi[i][qp]*(*projected_solution) (dof_indices[i]);
565  grad_u_fine += dphi[i][qp]*system.current_solution (dof_indices[i]);
566  grad_u_coarse += dphi[i][qp]*(*projected_solution) (dof_indices[i]);
567 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
568  grad2_u_fine += d2phi[i][qp]*system.current_solution (dof_indices[i]);
569  grad2_u_coarse += d2phi[i][qp]*(*projected_solution) (dof_indices[i]);
570 #endif
571  }
572 
573  // Compute the value of the error at this quadrature point
574  const Number val_error = u_fine - u_coarse;
575 
576  // Add the squares of the error to each contribution
577  if (system_i_norm.type(var) == L2 ||
578  system_i_norm.type(var) == H1 ||
579  system_i_norm.type(var) == H2)
580  {
581  L2normsq += JxW[qp] * system_i_norm.weight_sq(var) *
582  TensorTools::norm_sq(val_error);
583  libmesh_assert_greater_equal (L2normsq, 0.);
584  }
585 
586 
587  // Compute the value of the error in the gradient at this
588  // quadrature point
589  if (system_i_norm.type(var) == H1 ||
590  system_i_norm.type(var) == H2 ||
591  system_i_norm.type(var) == H1_SEMINORM)
592  {
593  Gradient grad_error = grad_u_fine - grad_u_coarse;
594 
595  H1seminormsq += JxW[qp] * system_i_norm.weight_sq(var) *
596  grad_error.norm_sq();
597  libmesh_assert_greater_equal (H1seminormsq, 0.);
598  }
599 
600  // Compute the value of the error in the hessian at this
601  // quadrature point
602  if (system_i_norm.type(var) == H2 ||
603  system_i_norm.type(var) == H2_SEMINORM)
604  {
605 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
606  Tensor grad2_error = grad2_u_fine - grad2_u_coarse;
607 
608  H2seminormsq += JxW[qp] * system_i_norm.weight_sq(var) *
609  grad2_error.norm_sq();
610  libmesh_assert_greater_equal (H2seminormsq, 0.);
611 #else
612  libmesh_error_msg
613  ("libMesh was not configured with --enable-second");
614 #endif
615  }
616  } // end qp loop
617 
618  if (system_i_norm.type(var) == L2 ||
619  system_i_norm.type(var) == H1 ||
620  system_i_norm.type(var) == H2)
621  (*err_vec)[e_id] +=
622  static_cast<ErrorVectorReal>(L2normsq);
623  if (system_i_norm.type(var) == H1 ||
624  system_i_norm.type(var) == H2 ||
625  system_i_norm.type(var) == H1_SEMINORM)
626  (*err_vec)[e_id] +=
627  static_cast<ErrorVectorReal>(H1seminormsq);
628 
629 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
630  if (system_i_norm.type(var) == H2 ||
631  system_i_norm.type(var) == H2_SEMINORM)
632  (*err_vec)[e_id] +=
633  static_cast<ErrorVectorReal>(H2seminormsq);
634 #endif
635  } // End loop over active local elements
636  } // End loop over variables
637 
638  // Don't bother projecting the solution; we'll restore from backup
639  // after coarsening
640  system.project_solution_on_reinit() = false;
641  }
642 
643 
644  // Uniformly coarsen the mesh, without projecting the solution
646 
647  for (unsigned int i = 0; i != number_h_refinements; ++i)
648  {
649  mesh_refinement.uniformly_coarsen(1);
650  // FIXME - should the reinits here be necessary? - RHS
651  es.reinit();
652  }
653 
654  for (unsigned int i = 0; i != number_p_refinements; ++i)
655  {
656  mesh_refinement.uniformly_p_coarsen(1);
657  es.reinit();
658  }
659 
660  // We should be back where we started
661  libmesh_assert_equal_to (n_coarse_elem, mesh.n_elem());
662 
663  // Each processor has now computed the error contributions
664  // for its local elements. We need to sum the vector
665  // and then take the square-root of each component. Note
666  // that we only need to sum if we are running on multiple
667  // processors, and we only need to take the square-root
668  // if the value is nonzero. There will in general be many
669  // zeros for the inactive elements.
670 
671  if (error_per_cell)
672  {
673  // First sum the vector of estimated error values
674  this->reduce_error(*error_per_cell, es.comm());
675 
676  // Compute the square-root of each component.
677  LOG_SCOPE("std::sqrt()", "UniformRefinementEstimator");
678  for (std::size_t i=0; i<error_per_cell->size(); i++)
679  if ((*error_per_cell)[i] != 0.)
680  (*error_per_cell)[i] = std::sqrt((*error_per_cell)[i]);
681  }
682  else
683  {
684  for (ErrorMap::iterator it = errors_per_cell->begin();
685  it != errors_per_cell->end(); ++it)
686  {
687  ErrorVector * e = it->second;
688  // First sum the vector of estimated error values
689  this->reduce_error(*e, es.comm());
690 
691  // Compute the square-root of each component.
692  LOG_SCOPE("std::sqrt()", "UniformRefinementEstimator");
693  for (std::size_t i=0; i<e->size(); i++)
694  if ((*e)[i] != 0.)
695  (*e)[i] = std::sqrt((*e)[i]);
696  }
697  }
698 
699  // Restore old solutions and clean up the heap
700  for (std::size_t i=0; i != system_list.size(); ++i)
701  {
702  System & system = *system_list[i];
703 
704  system.project_solution_on_reinit() = old_projection_settings[i];
705 
706  // Restore the coarse solution vectors and delete their copies
707  *system.solution = *coarse_solutions[i];
708  delete coarse_solutions[i];
709  *system.current_local_solution = *coarse_local_solutions[i];
710  delete coarse_local_solutions[i];
711  delete projected_solutions[i];
712 
713  for (System::vectors_iterator vec = system.vectors_begin(); vec !=
714  system.vectors_end(); ++vec)
715  {
716  // The (string) name of this vector
717  const std::string & var_name = vec->first;
718 
719  system.get_vector(var_name) = *coarse_vectors[i][var_name];
720 
721  coarse_vectors[i][var_name]->clear();
722  delete coarse_vectors[i][var_name];
723  }
724  }
725 
726  // Restore old partitioner settings
727  mesh.partitioner().reset(old_partitioner.release());
728 }
unsigned char number_h_refinements
How many h refinements to perform to get the fine grid.
const unsigned int invalid_uint
A number which is used quite often to represent an invalid or uninitialized value.
Definition: libmesh.h:184
SystemNorm error_norm
When estimating the error in a single system, the error_norm is used to control the scaling and norm ...
Real norm_sq() const
Definition: type_tensor.h:1270
unsigned int dim
ImplicitSystem & sys
MeshBase & mesh
const unsigned int n_vars
Definition: tecplot_io.C:68
DIE A HORRIBLE DEATH HERE typedef float ErrorVectorReal
libmesh_assert(j)
static UniquePtr< NumericVector< T > > build(const Parallel::Communicator &comm, const SolverPackage solver_package=libMesh::default_solver_package())
Builds a NumericVector on the processors in communicator comm using the linear solver package specifi...
NumberVectorValue Gradient
Real norm_sq() const
Definition: type_vector.h:940
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...
unsigned char number_p_refinements
How many p refinements to perform to get the fine grid.
std::map< std::string, NumericVector< Number > * >::iterator vectors_iterator
Vector iterator typedefs.
Definition: system.h:748
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
NumberTensorValue Tensor
FEMNormType type(unsigned int var) const
Definition: system_norm.h:268
uint8_t dof_id_type
Definition: id_types.h:64
static UniquePtr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.
void libMesh::UniformRefinementEstimator::estimate_error ( const System system,
ErrorVector error_per_cell,
const NumericVector< Number > *  solution_vector = libmesh_nullptr,
bool  estimate_parent_error = false 
)
virtual

This function does uniform refinements and a solve to get an improved solution on each cell, then estimates the error by integrating differences between the coarse and fine solutions.

system.solve() must be called, and so should have no side effects.

Only the provided system is solved on the refined mesh; for problems decoupled into multiple systems, use of estimate_errors() should be more reliable.

The estimated error is output in the vector error_per_cell

Implements libMesh::ErrorEstimator.

Definition at line 49 of file uniform_refinement_estimator.C.

References _estimate_error(), and libmesh_nullptr.

Referenced by main(), and ~UniformRefinementEstimator().

53 {
54  LOG_SCOPE("estimate_error()", "UniformRefinementEstimator");
55  std::map<const System *, const NumericVector<Number> *> solution_vectors;
56  solution_vectors[&_system] = solution_vector;
58  &_system,
59  &error_per_cell,
62  &solution_vectors,
63  estimate_parent_error);
64 }
const class libmesh_nullptr_t libmesh_nullptr
virtual void _estimate_error(const EquationSystems *equation_systems, const System *system, ErrorVector *error_per_cell, std::map< std::pair< const System *, unsigned int >, ErrorVector * > *errors_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)
The code for estimate_error and both estimate_errors versions is very similar, so we use the same fun...
void libMesh::UniformRefinementEstimator::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

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 from libMesh::ErrorEstimator.

Definition at line 66 of file uniform_refinement_estimator.C.

References _estimate_error(), and libmesh_nullptr.

Referenced by ~UniformRefinementEstimator().

71 {
72  LOG_SCOPE("estimate_errors()", "UniformRefinementEstimator");
73  this->_estimate_error (&_es,
75  &error_per_cell,
77  &error_norms,
78  solution_vectors,
79  estimate_parent_error);
80 }
const class libmesh_nullptr_t libmesh_nullptr
virtual void _estimate_error(const EquationSystems *equation_systems, const System *system, ErrorVector *error_per_cell, std::map< std::pair< const System *, unsigned int >, ErrorVector * > *errors_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)
The code for estimate_error and both estimate_errors versions is very similar, so we use the same fun...
void libMesh::UniformRefinementEstimator::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 
)
virtual

Currently this function ignores the component_scale member variable, because it calculates each error individually, unscaled.

The user selects which errors get computed by filling a map with error vectors: If errors_per_cell[&system][v] exists, it will be filled with the error values in variable v of system

Reimplemented from libMesh::ErrorEstimator.

Definition at line 82 of file uniform_refinement_estimator.C.

References _estimate_error(), and libmesh_nullptr.

86 {
87  LOG_SCOPE("estimate_errors()", "UniformRefinementEstimator");
88  this->_estimate_error (&_es,
91  &errors_per_cell,
93  solution_vectors,
94  estimate_parent_error);
95 }
const class libmesh_nullptr_t libmesh_nullptr
virtual void _estimate_error(const EquationSystems *equation_systems, const System *system, ErrorVector *error_per_cell, std::map< std::pair< const System *, unsigned int >, ErrorVector * > *errors_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)
The code for estimate_error and both estimate_errors versions is very similar, so we use the same fun...
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 _estimate_error(), libMesh::WeightedPatchRecoveryErrorEstimator::estimate_error(), libMesh::JumpErrorEstimator::estimate_error(), libMesh::PatchRecoveryErrorEstimator::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 }
virtual ErrorEstimatorType libMesh::UniformRefinementEstimator::type ( ) const
virtual
Returns
The type for the ErrorEstimator subclass.

Implements libMesh::ErrorEstimator.

Definition at line 111 of file uniform_refinement_estimator.h.

References libMesh::UNIFORM_REFINEMENT.

Member Data Documentation

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

unsigned char libMesh::UniformRefinementEstimator::number_h_refinements

How many h refinements to perform to get the fine grid.

Definition at line 117 of file uniform_refinement_estimator.h.

Referenced by _estimate_error().

unsigned char libMesh::UniformRefinementEstimator::number_p_refinements

How many p refinements to perform to get the fine grid.

Definition at line 122 of file uniform_refinement_estimator.h.

Referenced by _estimate_error().


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