libMesh
uniform_refinement_estimator.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 // Local Includes
20 #include "libmesh/dof_map.h"
21 #include "libmesh/elem.h"
22 #include "libmesh/equation_systems.h"
23 #include "libmesh/error_vector.h"
24 #include "libmesh/fe.h"
25 #include "libmesh/libmesh_common.h"
26 #include "libmesh/libmesh_logging.h"
27 #include "libmesh/mesh_base.h"
28 #include "libmesh/mesh_refinement.h"
29 #include "libmesh/numeric_vector.h"
30 #include "libmesh/quadrature.h"
31 #include "libmesh/system.h"
32 #include "libmesh/uniform_refinement_estimator.h"
33 #include "libmesh/partitioner.h"
34 #include "libmesh/tensor_tools.h"
35 #include "libmesh/enum_error_estimator_type.h"
36 #include "libmesh/enum_norm_type.h"
37 #include "libmesh/int_range.h"
38 
39 // C++ includes
40 #include <algorithm> // for std::fill
41 #include <sstream>
42 #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
43 #include <cmath> // for sqrt
44 #include <memory>
45 
46 #ifdef LIBMESH_ENABLE_AMR
47 
48 namespace libMesh
49 {
50 
51 //-----------------------------------------------------------------
52 // ErrorEstimator implementations
53 
56  number_h_refinements(1),
57  number_p_refinements(0),
58  _extra_order(1)
59 {
60  error_norm = H1;
61 }
62 
63 
64 
66 {
67  return UNIFORM_REFINEMENT;
68 }
69 
70 
72  ErrorVector & error_per_cell,
73  const NumericVector<Number> * solution_vector,
74  bool estimate_parent_error)
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 }
87 
89  ErrorVector & error_per_cell,
90  const std::map<const System *, SystemNorm> & error_norms,
91  const std::map<const System *, const NumericVector<Number> *> * solution_vectors,
92  bool estimate_parent_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 }
103 
105  ErrorMap & errors_per_cell,
106  const std::map<const System *, const NumericVector<Number> *> * solution_vectors,
107  bool estimate_parent_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 }
118 
120  const System * _system,
121  ErrorVector * error_per_cell,
122  ErrorMap * errors_per_cell,
123  const std::map<const System *, SystemNorm> * _error_norms,
124  const std::map<const System *, const NumericVector<Number> *> * solution_vectors,
125  bool)
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 }
742 
743 } // namespace libMesh
744 
745 #endif // #ifdef LIBMESH_ENABLE_AMR
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:196
std::map< std::string, std::unique_ptr< NumericVector< Number > >, std::less<> >::iterator vectors_iterator
Vector iterator typedefs.
Definition: system.h:757
vectors_iterator vectors_end()
End of vectors container.
Definition: system.h:2576
unsigned char number_h_refinements
How many h refinements to perform to get the fine grid.
void uniformly_p_refine(unsigned int n=1)
Uniformly p refines the mesh n times.
This is the EquationSystems class.
const Elem * parent() const
Definition: elem.h:3030
unsigned int n_systems() const
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:310
SystemNorm error_norm
When estimating the error in a single system, the error_norm is used to control the scaling and norm ...
virtual void get(const std::vector< numeric_index_type > &index, T *values) const
Access multiple components at once.
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Definition: dof_map.C:2164
unsigned int dim
unsigned int n_qois() const
Number of currently active quantities of interest.
Definition: system.h:2621
The ErrorVector is a specialization of the StatisticsVector for error data computed on a finite eleme...
Definition: error_vector.h:50
virtual void disable_cache()
Avoids use of any cached data that might affect any solve result.
Definition: system.h:2618
ErrorEstimatorType
Defines an enum for the different types of error estimators which are available.
const FEType & variable_type(const unsigned int c) const
Definition: dof_map.h:2220
const EquationSystems & get_equation_systems() const
Definition: system.h:721
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
auto norm_sq() const -> decltype(std::norm(T()))
Definition: type_tensor.h:1348
MeshBase & mesh
void uniformly_p_coarsen(unsigned int n=1)
Attempts to uniformly p coarsen the mesh n times.
const Parallel::Communicator & comm() const
This class defines a norm/seminorm to be applied to a NumericVector which contains coefficients in a ...
Definition: system_norm.h:49
This class defines a vector in LIBMESH_DIM dimensional Real or Complex space.
virtual void adjoint_solve(const QoISet &qoi_indices=QoISet())
Call adjoint_solve on all the individual equation systems.
The libMesh namespace provides an interface to certain functionality in the library.
vectors_iterator vectors_begin()
Beginning of vectors container.
Definition: system.h:2564
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...
const T_sys & get_system(std::string_view name) const
DIE A HORRIBLE DEATH HERE typedef float ErrorVectorReal
This is the MeshBase class.
Definition: mesh_base.h:75
Number current_solution(const dof_id_type global_dof_number) const
Definition: system.C:165
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:179
std::unique_ptr< QBase > default_quadrature_rule(const unsigned int dim, const int extraorder=0) const
Definition: fe_type.C:34
auto norm_sq() const -> decltype(std::norm(T()))
Definition: type_vector.h:926
Implements (adaptive) mesh refinement algorithms for a MeshBase.
int _extra_order
Extra order to use for quadrature rule.
dof_id_type id() const
Definition: dof_object.h:828
unsigned int n_vars
Manages consistently variables, degrees of freedom, and coefficient vectors.
Definition: system.h:96
virtual void reinit()
Handle any mesh changes and reinitialize all the systems on the updated mesh.
FEMNormType type(unsigned int var) const
Definition: system_norm.C:112
void uniformly_coarsen(unsigned int n=1)
Attempts to uniformly coarsen the mesh n times.
std::unique_ptr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1593
static std::unique_ptr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.
This class holds functions that will estimate the error in a finite element solution on a given mesh...
libmesh_assert(ctx)
virtual std::pair< unsigned int, Real > adjoint_solve(const QoISet &qoi_indices=QoISet())
Solves the adjoint system, for the specified qoi indices, or for every qoi if qoi_indices is nullptr...
Definition: system.h:2647
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...
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...
unsigned char number_p_refinements
How many p refinements to perform to get the fine grid.
virtual void solve()
Solves the system.
Definition: system.h:347
const NumericVector< Number > * request_vector(std::string_view vec_name) const
Definition: system.C:891
virtual void solve()
Call solve on all the individual equation systems.
Real weight_sq(unsigned int var) const
Definition: system_norm.C:177
virtual void update()
Update the local values to reflect the solution on neighboring processors.
Definition: system.C:510
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual void swap(NumericVector< T > &v)
Swaps the contents of this with v.
auto norm_sq(const T &a) -> decltype(std::norm(a))
Definition: tensor_tools.h:104
bool & project_solution_on_reinit(void)
Tells the System whether or not to project the solution vector onto new grids when the system is rein...
Definition: system.h:838
const MeshBase & get_mesh() const
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...
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:140
static std::unique_ptr< NumericVector< T > > build(const Parallel::Communicator &comm, SolverPackage solver_package=libMesh::default_solver_package(), ParallelType parallel_type=AUTOMATIC)
Builds a NumericVector on the processors in communicator comm using the linear solver package specifi...
NumericVector< Number > & get_adjoint_solution(unsigned int i=0)
Definition: system.C:1245
std::unique_ptr< NumericVector< Number > > current_local_solution
All the values I need to compute my contribution to the simulation at hand.
Definition: system.h:1605
virtual void clear()
Restores the NumericVector<T> to a pristine state.
unsigned int n_vars() const
Definition: system.h:2430
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_...
const DofMap & get_dof_map() const
Definition: system.h:2374
virtual ErrorEstimatorType type() const override
const std::vector< dof_id_type > & get_send_list() const
Definition: dof_map.h:526
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:117
const NumericVector< Number > & get_vector(std::string_view vec_name) const
Definition: system.C:943
This class defines a tensor in LIBMESH_DIM dimensional Real or Complex space.
uint8_t dof_id_type
Definition: id_types.h:67
void uniformly_refine(unsigned int n=1)
Uniformly refines the mesh n times.