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_elements_begin(), libMesh::MeshBase::active_local_elements_end(), 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::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  MeshBase::const_element_iterator elem_it = mesh.active_local_elements_begin();
515  const MeshBase::const_element_iterator elem_end = mesh.active_local_elements_end();
516 
517  for (; elem_it != elem_end; ++elem_it)
518  {
519  // e is necessarily an active element on the local processor
520  const Elem * elem = *elem_it;
521 
522  // Find the element id for the corresponding coarse grid element
523  const Elem * coarse = elem;
524  dof_id_type e_id = coarse->id();
525  while (e_id >= max_coarse_elem_id)
526  {
527  libmesh_assert (coarse->parent());
528  coarse = coarse->parent();
529  e_id = coarse->id();
530  }
531 
532  Real L2normsq = 0., H1seminormsq = 0.;
533 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
534  Real H2seminormsq = 0.;
535 #endif
536 
537  // reinitialize the element-specific data
538  // for the current element
539  fe->reinit (elem);
540 
541  // Get the local to global degree of freedom maps
542  dof_map.dof_indices (elem, dof_indices, var);
543 
544  // The number of quadrature points
545  const unsigned int n_qp = qrule->n_points();
546 
547  // The number of shape functions
548  const unsigned int n_sf =
549  cast_int<unsigned int>(dof_indices.size());
550 
551  //
552  // Begin the loop over the Quadrature points.
553  //
554  for (unsigned int qp=0; qp<n_qp; qp++)
555  {
556  Number u_fine = 0., u_coarse = 0.;
557 
558  Gradient grad_u_fine, grad_u_coarse;
559 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
560  Tensor grad2_u_fine, grad2_u_coarse;
561 #endif
562 
563  // Compute solution values at the current
564  // quadrature point. This requires a sum
565  // over all the shape functions evaluated
566  // at the quadrature point.
567  for (unsigned int i=0; i<n_sf; i++)
568  {
569  u_fine += phi[i][qp]*system.current_solution (dof_indices[i]);
570  u_coarse += phi[i][qp]*(*projected_solution) (dof_indices[i]);
571  grad_u_fine += dphi[i][qp]*system.current_solution (dof_indices[i]);
572  grad_u_coarse += dphi[i][qp]*(*projected_solution) (dof_indices[i]);
573 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
574  grad2_u_fine += d2phi[i][qp]*system.current_solution (dof_indices[i]);
575  grad2_u_coarse += d2phi[i][qp]*(*projected_solution) (dof_indices[i]);
576 #endif
577  }
578 
579  // Compute the value of the error at this quadrature point
580  const Number val_error = u_fine - u_coarse;
581 
582  // Add the squares of the error to each contribution
583  if (system_i_norm.type(var) == L2 ||
584  system_i_norm.type(var) == H1 ||
585  system_i_norm.type(var) == H2)
586  {
587  L2normsq += JxW[qp] * system_i_norm.weight_sq(var) *
588  TensorTools::norm_sq(val_error);
589  libmesh_assert_greater_equal (L2normsq, 0.);
590  }
591 
592 
593  // Compute the value of the error in the gradient at this
594  // quadrature point
595  if (system_i_norm.type(var) == H1 ||
596  system_i_norm.type(var) == H2 ||
597  system_i_norm.type(var) == H1_SEMINORM)
598  {
599  Gradient grad_error = grad_u_fine - grad_u_coarse;
600 
601  H1seminormsq += JxW[qp] * system_i_norm.weight_sq(var) *
602  grad_error.norm_sq();
603  libmesh_assert_greater_equal (H1seminormsq, 0.);
604  }
605 
606  // Compute the value of the error in the hessian at this
607  // quadrature point
608  if (system_i_norm.type(var) == H2 ||
609  system_i_norm.type(var) == H2_SEMINORM)
610  {
611 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
612  Tensor grad2_error = grad2_u_fine - grad2_u_coarse;
613 
614  H2seminormsq += JxW[qp] * system_i_norm.weight_sq(var) *
615  grad2_error.norm_sq();
616  libmesh_assert_greater_equal (H2seminormsq, 0.);
617 #else
618  libmesh_error_msg
619  ("libMesh was not configured with --enable-second");
620 #endif
621  }
622  } // end qp loop
623 
624  if (system_i_norm.type(var) == L2 ||
625  system_i_norm.type(var) == H1 ||
626  system_i_norm.type(var) == H2)
627  (*err_vec)[e_id] +=
628  static_cast<ErrorVectorReal>(L2normsq);
629  if (system_i_norm.type(var) == H1 ||
630  system_i_norm.type(var) == H2 ||
631  system_i_norm.type(var) == H1_SEMINORM)
632  (*err_vec)[e_id] +=
633  static_cast<ErrorVectorReal>(H1seminormsq);
634 
635 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
636  if (system_i_norm.type(var) == H2 ||
637  system_i_norm.type(var) == H2_SEMINORM)
638  (*err_vec)[e_id] +=
639  static_cast<ErrorVectorReal>(H2seminormsq);
640 #endif
641  } // End loop over active local elements
642  } // End loop over variables
643 
644  // Don't bother projecting the solution; we'll restore from backup
645  // after coarsening
646  system.project_solution_on_reinit() = false;
647  }
648 
649 
650  // Uniformly coarsen the mesh, without projecting the solution
652 
653  for (unsigned int i = 0; i != number_h_refinements; ++i)
654  {
655  mesh_refinement.uniformly_coarsen(1);
656  // FIXME - should the reinits here be necessary? - RHS
657  es.reinit();
658  }
659 
660  for (unsigned int i = 0; i != number_p_refinements; ++i)
661  {
662  mesh_refinement.uniformly_p_coarsen(1);
663  es.reinit();
664  }
665 
666  // We should be back where we started
667  libmesh_assert_equal_to (n_coarse_elem, mesh.n_elem());
668 
669  // Each processor has now computed the error contributions
670  // for its local elements. We need to sum the vector
671  // and then take the square-root of each component. Note
672  // that we only need to sum if we are running on multiple
673  // processors, and we only need to take the square-root
674  // if the value is nonzero. There will in general be many
675  // zeros for the inactive elements.
676 
677  if (error_per_cell)
678  {
679  // First sum the vector of estimated error values
680  this->reduce_error(*error_per_cell, es.comm());
681 
682  // Compute the square-root of each component.
683  LOG_SCOPE("std::sqrt()", "UniformRefinementEstimator");
684  for (std::size_t i=0; i<error_per_cell->size(); i++)
685  if ((*error_per_cell)[i] != 0.)
686  (*error_per_cell)[i] = std::sqrt((*error_per_cell)[i]);
687  }
688  else
689  {
690  for (ErrorMap::iterator it = errors_per_cell->begin();
691  it != errors_per_cell->end(); ++it)
692  {
693  ErrorVector * e = it->second;
694  // First sum the vector of estimated error values
695  this->reduce_error(*e, es.comm());
696 
697  // Compute the square-root of each component.
698  LOG_SCOPE("std::sqrt()", "UniformRefinementEstimator");
699  for (std::size_t i=0; i<e->size(); i++)
700  if ((*e)[i] != 0.)
701  (*e)[i] = std::sqrt((*e)[i]);
702  }
703  }
704 
705  // Restore old solutions and clean up the heap
706  for (std::size_t i=0; i != system_list.size(); ++i)
707  {
708  System & system = *system_list[i];
709 
710  system.project_solution_on_reinit() = old_projection_settings[i];
711 
712  // Restore the coarse solution vectors and delete their copies
713  *system.solution = *coarse_solutions[i];
714  delete coarse_solutions[i];
715  *system.current_local_solution = *coarse_local_solutions[i];
716  delete coarse_local_solutions[i];
717  delete projected_solutions[i];
718 
719  for (System::vectors_iterator vec = system.vectors_begin(); vec !=
720  system.vectors_end(); ++vec)
721  {
722  // The (string) name of this vector
723  const std::string & var_name = vec->first;
724 
725  system.get_vector(var_name) = *coarse_vectors[i][var_name];
726 
727  coarse_vectors[i][var_name]->clear();
728  delete coarse_vectors[i][var_name];
729  }
730  }
731 
732  // Restore old partitioner settings
733  mesh.partitioner().reset(old_partitioner.release());
734 }
const unsigned int invalid_uint
Definition: libmesh.h:184
Real norm_sq() const
Definition: type_tensor.h:1262
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())
NumberVectorValue Gradient
Real norm_sq() const
Definition: type_vector.h:932
void reduce_error(std::vector< ErrorVectorReal > &error_per_cell, const Parallel::Communicator &comm LIBMESH_CAN_DEFAULT_TO_COMMWORLD) const
std::map< std::string, NumericVector< Number > * >::iterator vectors_iterator
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)
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::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
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: