libMesh::UniformRefinementEstimator Class Reference

#include <uniform_refinement_estimator.h>

Inheritance diagram for libMesh::UniformRefinementEstimator:

Public Types

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

Public Member Functions

 UniformRefinementEstimator ()
 
 ~UniformRefinementEstimator ()
 
virtual void estimate_error (const System &system, ErrorVector &error_per_cell, const NumericVector< Number > *solution_vector=libmesh_nullptr, bool estimate_parent_error=false) libmesh_override
 
virtual void estimate_errors (const EquationSystems &equation_systems, ErrorVector &error_per_cell, const std::map< const System *, SystemNorm > &error_norms, const std::map< const System *, const NumericVector< Number > * > *solution_vectors=libmesh_nullptr, bool estimate_parent_error=false) libmesh_override
 
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
 
virtual ErrorEstimatorType type () const libmesh_override
 

Public Attributes

unsigned char number_h_refinements
 
unsigned char number_p_refinements
 
SystemNorm error_norm
 

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)
 
void reduce_error (std::vector< ErrorVectorReal > &error_per_cell, const Parallel::Communicator &comm LIBMESH_CAN_DEFAULT_TO_COMMWORLD) const
 

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 ( )
inline

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.

libMesh::UniformRefinementEstimator::~UniformRefinementEstimator ( )
inline

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

35 {
36  // This function must be run on all processors at once
37  // parallel_object_only();
38 
39  // Each processor has now computed the error contributions
40  // for its local elements. We may need to sum the vector to
41  // recover the error for each element.
42 
43  comm.sum(error_per_cell);
44 }
virtual ErrorEstimatorType libMesh::UniformRefinementEstimator::type ( ) const
inlinevirtual
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(), 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: