libMesh
Public Types | Public Member Functions | Public Attributes | Protected Member Functions | Protected Attributes | 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 >, std::unique_ptr< 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 (const UniformRefinementEstimator &)=default
 Copy/move ctor, copy/move assignment operator, and destructor are all explicitly defaulted for this simple class. More...
 
 UniformRefinementEstimator (UniformRefinementEstimator &&)=default
 
UniformRefinementEstimatoroperator= (const UniformRefinementEstimator &)=default
 
UniformRefinementEstimatoroperator= (UniformRefinementEstimator &&)=default
 
virtual ~UniformRefinementEstimator ()=default
 
virtual void estimate_error (const System &system, ErrorVector &error_per_cell, const NumericVector< Number > *solution_vector=nullptr, bool estimate_parent_error=false) 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=nullptr, bool estimate_parent_error=false) 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=nullptr, bool estimate_parent_error=false) override
 Currently this function ignores the component_scale member variable, because it calculates each error individually, unscaled. More...
 
void extra_quadrature_order (const int extraorder)
 Increases or decreases the order of the quadrature rule used for numerical integration. More...
 
virtual ErrorEstimatorType type () const 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, ErrorMap *errors_per_cell, const std::map< const System *, SystemNorm > *error_norms, const std::map< const System *, const NumericVector< Number > *> *solution_vectors=nullptr, bool estimate_parent_error=false)
 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) 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...
 

Protected Attributes

int _extra_order
 Extra order to use for quadrature rule. 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

◆ ErrorMap

typedef std::map<std::pair<const System *, unsigned int>, std::unique_ptr<ErrorVector> > libMesh::ErrorEstimator::ErrorMap
inherited

When calculating many error vectors at once, we need a data structure to hold them all.

Definition at line 121 of file error_estimator.h.

Constructor & Destructor Documentation

◆ UniformRefinementEstimator() [1/3]

libMesh::UniformRefinementEstimator::UniformRefinementEstimator ( )

Constructor.

Sets the most common default parameter values.

Definition at line 54 of file uniform_refinement_estimator.C.

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

54  :
58  _extra_order(1)
59 {
60  error_norm = H1;
61 }
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 ...
int _extra_order
Extra order to use for quadrature rule.
ErrorEstimator()=default
Constructor.
unsigned char number_p_refinements
How many p refinements to perform to get the fine grid.

◆ UniformRefinementEstimator() [2/3]

libMesh::UniformRefinementEstimator::UniformRefinementEstimator ( const UniformRefinementEstimator )
default

Copy/move ctor, copy/move assignment operator, and destructor are all explicitly defaulted for this simple class.

◆ UniformRefinementEstimator() [3/3]

libMesh::UniformRefinementEstimator::UniformRefinementEstimator ( UniformRefinementEstimator &&  )
default

◆ ~UniformRefinementEstimator()

virtual libMesh::UniformRefinementEstimator::~UniformRefinementEstimator ( )
virtualdefault

Member Function Documentation

◆ _estimate_error()

void libMesh::UniformRefinementEstimator::_estimate_error ( const EquationSystems equation_systems,
const System system,
ErrorVector error_per_cell,
ErrorMap errors_per_cell,
const std::map< const System *, SystemNorm > *  error_norms,
const std::map< const System *, const NumericVector< Number > *> *  solution_vectors = nullptr,
bool  estimate_parent_error = false 
)
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 119 of file uniform_refinement_estimator.C.

References _extra_order, 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::NumericVector< T >::get(), 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::GHOSTED, libMesh::H1, libMesh::H1_SEMINORM, libMesh::H2, libMesh::H2_SEMINORM, libMesh::DofObject::id(), libMesh::index_range(), libMesh::invalid_uint, libMesh::L2, libMesh::libmesh_assert(), libMesh::make_range(), mesh, libMesh::System::n_qois(), 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::System::project_solution_on_reinit(), libMesh::Real, libMesh::ErrorEstimator::reduce_error(), libMesh::EquationSystems::reinit(), libMesh::System::request_vector(), libMesh::System::solution, libMesh::EquationSystems::solve(), libMesh::System::solve(), std::sqrt(), libMesh::NumericVector< T >::swap(), 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().

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

◆ estimate_error()

void libMesh::UniformRefinementEstimator::estimate_error ( const System system,
ErrorVector error_per_cell,
const NumericVector< Number > *  solution_vector = nullptr,
bool  estimate_parent_error = false 
)
overridevirtual

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 71 of file uniform_refinement_estimator.C.

References _estimate_error().

Referenced by main().

75 {
76  LOG_SCOPE("estimate_error()", "UniformRefinementEstimator");
77  std::map<const System *, const NumericVector<Number> *> solution_vectors;
78  solution_vectors[&_system] = solution_vector;
79  this->_estimate_error (nullptr,
80  &_system,
81  &error_per_cell,
82  nullptr,
83  nullptr,
84  &solution_vectors,
85  estimate_parent_error);
86 }
virtual void _estimate_error(const EquationSystems *equation_systems, const System *system, ErrorVector *error_per_cell, ErrorMap *errors_per_cell, const std::map< const System *, SystemNorm > *error_norms, const std::map< const System *, const NumericVector< Number > *> *solution_vectors=nullptr, bool estimate_parent_error=false)
The code for estimate_error and both estimate_errors versions is very similar, so we use the same fun...

◆ estimate_errors() [1/2]

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 = nullptr,
bool  estimate_parent_error = false 
)
overridevirtual

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 88 of file uniform_refinement_estimator.C.

References _estimate_error().

93 {
94  LOG_SCOPE("estimate_errors()", "UniformRefinementEstimator");
95  this->_estimate_error (&_es,
96  nullptr,
97  &error_per_cell,
98  nullptr,
99  &error_norms,
100  solution_vectors,
101  estimate_parent_error);
102 }
virtual void _estimate_error(const EquationSystems *equation_systems, const System *system, ErrorVector *error_per_cell, ErrorMap *errors_per_cell, const std::map< const System *, SystemNorm > *error_norms, const std::map< const System *, const NumericVector< Number > *> *solution_vectors=nullptr, bool estimate_parent_error=false)
The code for estimate_error and both estimate_errors versions is very similar, so we use the same fun...

◆ estimate_errors() [2/2]

void libMesh::UniformRefinementEstimator::estimate_errors ( const EquationSystems equation_systems,
ErrorMap errors_per_cell,
const std::map< const System *, const NumericVector< Number > *> *  solution_vectors = nullptr,
bool  estimate_parent_error = false 
)
overridevirtual

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 104 of file uniform_refinement_estimator.C.

References _estimate_error().

108 {
109  LOG_SCOPE("estimate_errors()", "UniformRefinementEstimator");
110  this->_estimate_error (&_es,
111  nullptr,
112  nullptr,
113  &errors_per_cell,
114  nullptr,
115  solution_vectors,
116  estimate_parent_error);
117 }
virtual void _estimate_error(const EquationSystems *equation_systems, const System *system, ErrorVector *error_per_cell, ErrorMap *errors_per_cell, const std::map< const System *, SystemNorm > *error_norms, const std::map< const System *, const NumericVector< Number > *> *solution_vectors=nullptr, bool estimate_parent_error=false)
The code for estimate_error and both estimate_errors versions is very similar, so we use the same fun...

◆ extra_quadrature_order()

void libMesh::UniformRefinementEstimator::extra_quadrature_order ( const int  extraorder)
inline

Increases or decreases the order of the quadrature rule used for numerical integration.

The default extraorder is 1, because properly integrating L2 error requires integrating the squares of terms with order p+1, and 2p+2 is 1 higher than what we default to using for reasonable mass matrix integration.

Definition at line 118 of file uniform_refinement_estimator.h.

References _extra_order.

119  { _extra_order = extraorder; }
int _extra_order
Extra order to use for quadrature rule.

◆ operator=() [1/2]

UniformRefinementEstimator& libMesh::UniformRefinementEstimator::operator= ( const UniformRefinementEstimator )
default

◆ operator=() [2/2]

UniformRefinementEstimator& libMesh::UniformRefinementEstimator::operator= ( UniformRefinementEstimator &&  )
default

◆ reduce_error()

void libMesh::ErrorEstimator::reduce_error ( std::vector< ErrorVectorReal > &  error_per_cell,
const Parallel::Communicator comm 
) const
protectedinherited

This method takes the local error contributions in error_per_cell from each processor and combines them to get the global error vector.

Definition at line 32 of file error_estimator.C.

References TIMPI::Communicator::sum().

Referenced by _estimate_error(), libMesh::PatchRecoveryErrorEstimator::estimate_error(), libMesh::WeightedPatchRecoveryErrorEstimator::estimate_error(), libMesh::JumpErrorEstimator::estimate_error(), libMesh::AdjointRefinementEstimator::estimate_error(), and libMesh::ExactErrorEstimator::estimate_error().

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

◆ type()

ErrorEstimatorType libMesh::UniformRefinementEstimator::type ( ) const
overridevirtual
Returns
The type for the ErrorEstimator subclass.

Implements libMesh::ErrorEstimator.

Definition at line 65 of file uniform_refinement_estimator.C.

References libMesh::UNIFORM_REFINEMENT.

Member Data Documentation

◆ _extra_order

int libMesh::UniformRefinementEstimator::_extra_order
protected

Extra order to use for quadrature rule.

Definition at line 138 of file uniform_refinement_estimator.h.

Referenced by _estimate_error(), and extra_quadrature_order().

◆ error_norm

SystemNorm libMesh::ErrorEstimator::error_norm
inherited

When estimating the error in a single system, the error_norm is used to control the scaling and norm choice for each variable.

Not all estimators will support all norm choices. The default scaling is for all variables to be weighted equally. The default norm choice depends on the error estimator.

Part of this functionality was supported via component_scale and sobolev_order in older libMesh versions, and a small part was supported via component_mask in even older versions. Hopefully the encapsulation here will allow us to avoid changing this API again.

Definition at line 158 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(), libMesh::WeightedPatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::PatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::PatchRecoveryErrorEstimator::PatchRecoveryErrorEstimator(), libMesh::JumpErrorEstimator::reinit_sides(), and UniformRefinementEstimator().

◆ number_h_refinements

unsigned char libMesh::UniformRefinementEstimator::number_h_refinements

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

Definition at line 126 of file uniform_refinement_estimator.h.

Referenced by _estimate_error().

◆ number_p_refinements

unsigned char libMesh::UniformRefinementEstimator::number_p_refinements

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

Definition at line 131 of file uniform_refinement_estimator.h.

Referenced by _estimate_error().


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