uniform_refinement_estimator.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2018 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 // C++ includes
19 #include <algorithm> // for std::fill
20 #include <sstream>
21 #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
22 #include <cmath> // for sqrt
23 
24 
25 // Local Includes
26 #include "libmesh/dof_map.h"
27 #include "libmesh/elem.h"
29 #include "libmesh/error_vector.h"
30 #include "libmesh/fe.h"
31 #include "libmesh/libmesh_common.h"
33 #include "libmesh/mesh_base.h"
35 #include "libmesh/numeric_vector.h"
36 #include "libmesh/quadrature.h"
37 #include "libmesh/system.h"
39 #include "libmesh/partitioner.h"
40 #include "libmesh/tensor_tools.h"
42 #include "libmesh/enum_norm_type.h"
43 #include "libmesh/int_range.h"
44 
45 #ifdef LIBMESH_ENABLE_AMR
46 
47 namespace libMesh
48 {
49 
50 //-----------------------------------------------------------------
51 // ErrorEstimator implementations
52 
55  number_h_refinements(1),
56  number_p_refinements(0)
57 {
58  error_norm = H1;
59 }
60 
61 
62 
64 {
65  return UNIFORM_REFINEMENT;
66 }
67 
68 
70  ErrorVector & error_per_cell,
71  const NumericVector<Number> * solution_vector,
72  bool estimate_parent_error)
73 {
74  LOG_SCOPE("estimate_error()", "UniformRefinementEstimator");
75  std::map<const System *, const NumericVector<Number> *> solution_vectors;
76  solution_vectors[&_system] = solution_vector;
77  this->_estimate_error (nullptr,
78  &_system,
79  &error_per_cell,
80  nullptr,
81  nullptr,
82  &solution_vectors,
83  estimate_parent_error);
84 }
85 
87  ErrorVector & error_per_cell,
88  const std::map<const System *, SystemNorm> & error_norms,
89  const std::map<const System *, const NumericVector<Number> *> * solution_vectors,
90  bool estimate_parent_error)
91 {
92  LOG_SCOPE("estimate_errors()", "UniformRefinementEstimator");
93  this->_estimate_error (&_es,
94  nullptr,
95  &error_per_cell,
96  nullptr,
97  &error_norms,
98  solution_vectors,
99  estimate_parent_error);
100 }
101 
103  ErrorMap & errors_per_cell,
104  const std::map<const System *, const NumericVector<Number> *> * solution_vectors,
105  bool estimate_parent_error)
106 {
107  LOG_SCOPE("estimate_errors()", "UniformRefinementEstimator");
108  this->_estimate_error (&_es,
109  nullptr,
110  nullptr,
111  &errors_per_cell,
112  nullptr,
113  solution_vectors,
114  estimate_parent_error);
115 }
116 
118  const System * _system,
119  ErrorVector * error_per_cell,
120  ErrorMap * errors_per_cell,
121  const std::map<const System *, SystemNorm> * _error_norms,
122  const std::map<const System *, const NumericVector<Number> *> * solution_vectors,
123  bool)
124 {
125  // Get a vector of the Systems we're going to work on,
126  // and set up a error_norms map if necessary
127  std::vector<System *> system_list;
128  std::unique_ptr<std::map<const System *, SystemNorm>> error_norms =
129  std::unique_ptr<std::map<const System *, SystemNorm>>
130  (new 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 (unsigned int i=0; i != _es->n_systems(); ++i)
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 (unsigned int i=0; i!= _es->n_systems(); ++i)
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->find(std::make_pair(&sys, v)) ==
167  errors_per_cell->end())
168  continue;
169 
170  weights[v] = 1.0;
171  }
172  (*error_norms)[&sys] =
173  SystemNorm(std::vector<FEMNormType>(n_vars, error_norm.type(0)),
174  weights);
175  }
176  }
177  }
178  else
179  {
180  libmesh_assert(_system);
181  // We have to break the rules here, because we can't refine a const System
182  system_list.push_back(const_cast<System *>(_system));
183 
184  libmesh_assert(!_error_norms);
185  (*error_norms)[_system] = error_norm;
186  _error_norms = error_norms.get();
187  }
188 
189  // An EquationSystems reference will be convenient.
190  // We have to break the rules here, because we can't refine a const System
191  EquationSystems & es =
192  const_cast<EquationSystems &>(_system->get_equation_systems());
193 
194  // The current mesh
195  MeshBase & mesh = es.get_mesh();
196 
197  // The dimensionality of the mesh
198  const unsigned int dim = mesh.mesh_dimension();
199 
200  // Resize the error_per_cell vectors to be
201  // the number of elements, initialize them to 0.
202  if (error_per_cell)
203  {
204  error_per_cell->clear();
205  error_per_cell->resize (mesh.max_elem_id(), 0.);
206  }
207  else
208  {
209  libmesh_assert(errors_per_cell);
210  for (const auto & pr : *errors_per_cell)
211  {
212  ErrorVector * e = pr.second;
213  e->clear();
214  e->resize(mesh.max_elem_id(), 0.);
215  }
216  }
217 
218  // We'll want to back up all coarse grid vectors
219  std::vector<std::map<std::string, std::unique_ptr<NumericVector<Number>>>> coarse_vectors(system_list.size());
220  std::vector<std::unique_ptr<NumericVector<Number>>> coarse_solutions(system_list.size());
221  std::vector<std::unique_ptr<NumericVector<Number>>> coarse_local_solutions(system_list.size());
222  // And make copies of projected solutions
223  std::vector<std::unique_ptr<NumericVector<Number>>> projected_solutions(system_list.size());
224 
225  // And we'll need to temporarily change solution projection settings
226  std::vector<bool> old_projection_settings(system_list.size());
227 
228  // And it'll be best to avoid any repartitioning
229  std::unique_ptr<Partitioner> old_partitioner(mesh.partitioner().release());
230 
231  // And we can't allow any renumbering
232  const bool old_renumbering_setting = mesh.allow_renumbering();
233  mesh.allow_renumbering(false);
234 
235  for (auto i : index_range(system_list))
236  {
237  System & system = *system_list[i];
238 
239  // Check for valid error_norms
240  libmesh_assert (_error_norms->find(&system) !=
241  _error_norms->end());
242 
243  // Back up the solution vector
244  coarse_solutions[i] = system.solution->clone();
245  coarse_local_solutions[i] = system.current_local_solution->clone();
246 
247  // Back up all other coarse grid vectors
248  for (System::vectors_iterator vec = system.vectors_begin(); vec !=
249  system.vectors_end(); ++vec)
250  {
251  // The (string) name of this vector
252  const std::string & var_name = vec->first;
253 
254  coarse_vectors[i][var_name] = vec->second->clone();
255  }
256 
257  // Use a non-standard solution vector if necessary
258  if (solution_vectors &&
259  solution_vectors->find(&system) != solution_vectors->end() &&
260  solution_vectors->find(&system)->second &&
261  solution_vectors->find(&system)->second != system.solution.get())
262  {
263  NumericVector<Number> * newsol =
264  const_cast<NumericVector<Number> *>
265  (solution_vectors->find(&system)->second);
266  newsol->swap(*system.solution);
267  system.update();
268  }
269 
270  // Make sure the solution is projected when we refine the mesh
271  old_projection_settings[i] = system.project_solution_on_reinit();
272  system.project_solution_on_reinit() = true;
273  }
274 
275  // Find the number of coarse mesh elements, to make it possible
276  // to find correct coarse elem ids later
277  const dof_id_type max_coarse_elem_id = mesh.max_elem_id();
278 #ifndef NDEBUG
279  // n_coarse_elem is only used in an assertion later so
280  // avoid declaring it unless asserts are active.
281  const dof_id_type n_coarse_elem = mesh.n_elem();
282 #endif
283 
284  // Uniformly refine the mesh
285  MeshRefinement mesh_refinement(mesh);
286 
287  libmesh_assert (number_h_refinements > 0 || number_p_refinements > 0);
288 
289  // FIXME: this may break if there is more than one System
290  // on this mesh but estimate_error was still called instead of
291  // estimate_errors
292  for (unsigned int i = 0; i != number_h_refinements; ++i)
293  {
294  mesh_refinement.uniformly_refine(1);
295  es.reinit();
296  }
297 
298  for (unsigned int i = 0; i != number_p_refinements; ++i)
299  {
300  mesh_refinement.uniformly_p_refine(1);
301  es.reinit();
302  }
303 
304  for (auto i : index_range(system_list))
305  {
306  System & system = *system_list[i];
307 
308  // Copy the projected coarse grid solutions, which will be
309  // overwritten by solve()
310  projected_solutions[i] = NumericVector<Number>::build(system.comm());
311  projected_solutions[i]->init(system.solution->size(), true, SERIAL);
312  system.solution->localize(*projected_solutions[i],
313  system.get_dof_map().get_send_list());
314  }
315 
316  // Are we doing a forward or an adjoint solve?
317  bool solve_adjoint = false;
318  if (solution_vectors)
319  {
320  System * sys = system_list[0];
321  libmesh_assert (solution_vectors->find(sys) !=
322  solution_vectors->end());
323  const NumericVector<Number> * vec = solution_vectors->find(sys)->second;
324  for (unsigned int j=0; j != sys->n_qois(); ++j)
325  {
326  std::ostringstream adjoint_name;
327  adjoint_name << "adjoint_solution" << j;
328 
329  if (vec == sys->request_vector(adjoint_name.str()))
330  {
331  solve_adjoint = true;
332  break;
333  }
334  }
335  }
336 
337  // Get the uniformly refined solution.
338 
339  if (_es)
340  {
341  // Even if we had a decent preconditioner, valid matrix etc. before
342  // refinement, we don't any more.
343  for (unsigned int i=0; i != es.n_systems(); ++i)
344  es.get_system(i).disable_cache();
345 
346  // No specified vectors == forward solve
347  if (!solution_vectors)
348  es.solve();
349  else
350  {
351  libmesh_assert_equal_to (solution_vectors->size(), es.n_systems());
352  libmesh_assert (solution_vectors->find(system_list[0]) !=
353  solution_vectors->end());
354  libmesh_assert(solve_adjoint ||
355  (solution_vectors->find(system_list[0])->second ==
356  system_list[0]->solution.get()) ||
357  !solution_vectors->find(system_list[0])->second);
358 
359 #ifdef DEBUG
360  for (const auto & sys : system_list)
361  {
362  libmesh_assert (solution_vectors->find(sys) !=
363  solution_vectors->end());
364  const NumericVector<Number> * vec = solution_vectors->find(sys)->second;
365  if (solve_adjoint)
366  {
367  bool found_vec = false;
368  for (unsigned int j=0; j != sys->n_qois(); ++j)
369  {
370  std::ostringstream adjoint_name;
371  adjoint_name << "adjoint_solution" << j;
372 
373  if (vec == sys->request_vector(adjoint_name.str()))
374  {
375  found_vec = true;
376  break;
377  }
378  }
379  libmesh_assert(found_vec);
380  }
381  else
382  libmesh_assert(vec == sys->solution.get() || !vec);
383  }
384 #endif
385 
386  if (solve_adjoint)
387  {
388  std::vector<unsigned int> adjs(system_list.size(),
390  // Set up proper initial guesses
391  for (auto i : index_range(system_list))
392  {
393  System * sys = system_list[i];
394  libmesh_assert (solution_vectors->find(sys) !=
395  solution_vectors->end());
396  const NumericVector<Number> * vec = solution_vectors->find(sys)->second;
397  for (unsigned int j=0, n_qois = sys->n_qois();
398  j != n_qois; ++j)
399  {
400  std::ostringstream adjoint_name;
401  adjoint_name << "adjoint_solution" << j;
402 
403  if (vec == sys->request_vector(adjoint_name.str()))
404  {
405  adjs[i] = j;
406  break;
407  }
408  }
409  libmesh_assert_not_equal_to (adjs[i], libMesh::invalid_uint);
410  sys->get_adjoint_solution(adjs[i]) = *sys->solution;
411  }
412 
413  es.adjoint_solve();
414 
415  // Put the adjoint_solution into solution for
416  // comparisons
417  for (auto i : index_range(system_list))
418  {
419  system_list[i]->get_adjoint_solution(adjs[i]).swap(*system_list[i]->solution);
420  system_list[i]->update();
421  }
422  }
423  else
424  es.solve();
425  }
426  }
427  else
428  {
429  System * sys = system_list[0];
430 
431  // Even if we had a decent preconditioner, valid matrix etc. before
432  // refinement, we don't any more.
433  sys->disable_cache();
434 
435  // No specified vectors == forward solve
436  if (!solution_vectors)
437  sys->solve();
438  else
439  {
440  libmesh_assert (solution_vectors->find(sys) !=
441  solution_vectors->end());
442 
443  const NumericVector<Number> * vec = solution_vectors->find(sys)->second;
444 
445  libmesh_assert(solve_adjoint ||
446  (solution_vectors->find(sys)->second ==
447  sys->solution.get()) ||
448  !solution_vectors->find(sys)->second);
449 
450  if (solve_adjoint)
451  {
452  unsigned int adj = libMesh::invalid_uint;
453  for (unsigned int j=0, n_qois = sys->n_qois();
454  j != n_qois; ++j)
455  {
456  std::ostringstream adjoint_name;
457  adjoint_name << "adjoint_solution" << j;
458 
459  if (vec == sys->request_vector(adjoint_name.str()))
460  {
461  adj = j;
462  break;
463  }
464  }
465  libmesh_assert_not_equal_to (adj, libMesh::invalid_uint);
466 
467  // Set up proper initial guess
468  sys->get_adjoint_solution(adj) = *sys->solution;
469  sys->adjoint_solve();
470  // Put the adjoint_solution into solution for
471  // comparisons
472  sys->get_adjoint_solution(adj).swap(*sys->solution);
473  sys->update();
474  }
475  else
476  sys->solve();
477  }
478  }
479 
480  // Get the error in the uniformly refined solution(s).
481  for (auto sysnum : index_range(system_list))
482  {
483  System & system = *system_list[sysnum];
484 
485  unsigned int n_vars = system.n_vars();
486 
487  DofMap & dof_map = system.get_dof_map();
488 
489  const SystemNorm & system_i_norm =
490  _error_norms->find(&system)->second;
491 
492  NumericVector<Number> * projected_solution = projected_solutions[sysnum].get();
493 
494  // Loop over all the variables in the system
495  for (unsigned int var=0; var<n_vars; var++)
496  {
497  // Get the error vector to fill for this system and variable
498  ErrorVector * err_vec = error_per_cell;
499  if (!err_vec)
500  {
501  libmesh_assert(errors_per_cell);
502  err_vec =
503  (*errors_per_cell)[std::make_pair(&system,var)];
504  }
505 
506  // The type of finite element to use for this variable
507  const FEType & fe_type = dof_map.variable_type (var);
508 
509  // Finite element object for each fine element
510  std::unique_ptr<FEBase> fe (FEBase::build (dim, fe_type));
511 
512  // Build and attach an appropriate quadrature rule
513  std::unique_ptr<QBase> qrule = fe_type.default_quadrature_rule(dim);
514  fe->attach_quadrature_rule (qrule.get());
515 
516  const std::vector<Real> & JxW = fe->get_JxW();
517  const std::vector<std::vector<Real>> & phi = fe->get_phi();
518  const std::vector<std::vector<RealGradient>> & dphi =
519  fe->get_dphi();
520 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
521  const std::vector<std::vector<RealTensor>> & d2phi =
522  fe->get_d2phi();
523 #endif
524 
525  // The global DOF indices for the fine element
526  std::vector<dof_id_type> dof_indices;
527 
528  // Iterate over all the active elements in the fine mesh
529  // that live on this processor.
530  for (const auto & elem : mesh.active_local_element_ptr_range())
531  {
532  // Find the element id for the corresponding coarse grid element
533  const Elem * coarse = elem;
534  dof_id_type e_id = coarse->id();
535  while (e_id >= max_coarse_elem_id)
536  {
537  libmesh_assert (coarse->parent());
538  coarse = coarse->parent();
539  e_id = coarse->id();
540  }
541 
542  Real L2normsq = 0., H1seminormsq = 0.;
543 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
544  Real H2seminormsq = 0.;
545 #endif
546 
547  // reinitialize the element-specific data
548  // for the current element
549  fe->reinit (elem);
550 
551  // Get the local to global degree of freedom maps
552  dof_map.dof_indices (elem, dof_indices, var);
553 
554  // The number of quadrature points
555  const unsigned int n_qp = qrule->n_points();
556 
557  // The number of shape functions
558  const unsigned int n_sf =
559  cast_int<unsigned int>(dof_indices.size());
560 
561  //
562  // Begin the loop over the Quadrature points.
563  //
564  for (unsigned int qp=0; qp<n_qp; qp++)
565  {
566  Number u_fine = 0., u_coarse = 0.;
567 
568  Gradient grad_u_fine, grad_u_coarse;
569 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
570  Tensor grad2_u_fine, grad2_u_coarse;
571 #endif
572 
573  // Compute solution values at the current
574  // quadrature point. This requires a sum
575  // over all the shape functions evaluated
576  // at the quadrature point.
577  for (unsigned int i=0; i<n_sf; i++)
578  {
579  u_fine += phi[i][qp]*system.current_solution (dof_indices[i]);
580  u_coarse += phi[i][qp]*(*projected_solution) (dof_indices[i]);
581  grad_u_fine += dphi[i][qp]*system.current_solution (dof_indices[i]);
582  grad_u_coarse += dphi[i][qp]*(*projected_solution) (dof_indices[i]);
583 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
584  grad2_u_fine += d2phi[i][qp]*system.current_solution (dof_indices[i]);
585  grad2_u_coarse += d2phi[i][qp]*(*projected_solution) (dof_indices[i]);
586 #endif
587  }
588 
589  // Compute the value of the error at this quadrature point
590  const Number val_error = u_fine - u_coarse;
591 
592  // Add the squares of the error to each contribution
593  if (system_i_norm.type(var) == L2 ||
594  system_i_norm.type(var) == H1 ||
595  system_i_norm.type(var) == H2)
596  {
597  L2normsq += JxW[qp] * system_i_norm.weight_sq(var) *
598  TensorTools::norm_sq(val_error);
599  libmesh_assert_greater_equal (L2normsq, 0.);
600  }
601 
602 
603  // Compute the value of the error in the gradient at this
604  // quadrature point
605  if (system_i_norm.type(var) == H1 ||
606  system_i_norm.type(var) == H2 ||
607  system_i_norm.type(var) == H1_SEMINORM)
608  {
609  Gradient grad_error = grad_u_fine - grad_u_coarse;
610 
611  H1seminormsq += JxW[qp] * system_i_norm.weight_sq(var) *
612  grad_error.norm_sq();
613  libmesh_assert_greater_equal (H1seminormsq, 0.);
614  }
615 
616  // Compute the value of the error in the hessian at this
617  // quadrature point
618  if (system_i_norm.type(var) == H2 ||
619  system_i_norm.type(var) == H2_SEMINORM)
620  {
621 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
622  Tensor grad2_error = grad2_u_fine - grad2_u_coarse;
623 
624  H2seminormsq += JxW[qp] * system_i_norm.weight_sq(var) *
625  grad2_error.norm_sq();
626  libmesh_assert_greater_equal (H2seminormsq, 0.);
627 #else
628  libmesh_error_msg
629  ("libMesh was not configured with --enable-second");
630 #endif
631  }
632  } // end qp loop
633 
634  if (system_i_norm.type(var) == L2 ||
635  system_i_norm.type(var) == H1 ||
636  system_i_norm.type(var) == H2)
637  (*err_vec)[e_id] +=
638  static_cast<ErrorVectorReal>(L2normsq);
639  if (system_i_norm.type(var) == H1 ||
640  system_i_norm.type(var) == H2 ||
641  system_i_norm.type(var) == H1_SEMINORM)
642  (*err_vec)[e_id] +=
643  static_cast<ErrorVectorReal>(H1seminormsq);
644 
645 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
646  if (system_i_norm.type(var) == H2 ||
647  system_i_norm.type(var) == H2_SEMINORM)
648  (*err_vec)[e_id] +=
649  static_cast<ErrorVectorReal>(H2seminormsq);
650 #endif
651  } // End loop over active local elements
652  } // End loop over variables
653 
654  // Don't bother projecting the solution; we'll restore from backup
655  // after coarsening
656  system.project_solution_on_reinit() = false;
657  }
658 
659 
660  // Uniformly coarsen the mesh, without projecting the solution
661  libmesh_assert (number_h_refinements > 0 || number_p_refinements > 0);
662 
663  for (unsigned int i = 0; i != number_h_refinements; ++i)
664  {
665  mesh_refinement.uniformly_coarsen(1);
666  // FIXME - should the reinits here be necessary? - RHS
667  es.reinit();
668  }
669 
670  for (unsigned int i = 0; i != number_p_refinements; ++i)
671  {
672  mesh_refinement.uniformly_p_coarsen(1);
673  es.reinit();
674  }
675 
676  // We should be back where we started
677  libmesh_assert_equal_to (n_coarse_elem, mesh.n_elem());
678 
679  // Each processor has now computed the error contributions
680  // for its local elements. We need to sum the vector
681  // and then take the square-root of each component. Note
682  // that we only need to sum if we are running on multiple
683  // processors, and we only need to take the square-root
684  // if the value is nonzero. There will in general be many
685  // zeros for the inactive elements.
686 
687  if (error_per_cell)
688  {
689  // First sum the vector of estimated error values
690  this->reduce_error(*error_per_cell, es.comm());
691 
692  // Compute the square-root of each component.
693  LOG_SCOPE("std::sqrt()", "UniformRefinementEstimator");
694  for (auto & val : *error_per_cell)
695  if (val != 0.)
696  val = std::sqrt(val);
697  }
698  else
699  {
700  for (const auto & pr : *errors_per_cell)
701  {
702  ErrorVector & e = *(pr.second);
703  // First sum the vector of estimated error values
704  this->reduce_error(e, es.comm());
705 
706  // Compute the square-root of each component.
707  LOG_SCOPE("std::sqrt()", "UniformRefinementEstimator");
708  for (auto & val : e)
709  if (val != 0.)
710  val = std::sqrt(val);
711  }
712  }
713 
714  // Restore old solutions and clean up the heap
715  for (auto i : index_range(system_list))
716  {
717  System & system = *system_list[i];
718 
719  system.project_solution_on_reinit() = old_projection_settings[i];
720 
721  // Restore the coarse solution vectors and delete their copies
722  *system.solution = *coarse_solutions[i];
723  *system.current_local_solution = *coarse_local_solutions[i];
724 
725  for (System::vectors_iterator vec = system.vectors_begin(); vec !=
726  system.vectors_end(); ++vec)
727  {
728  // The (string) name of this vector
729  const std::string & var_name = vec->first;
730 
731  system.get_vector(var_name) = *coarse_vectors[i][var_name];
732 
733  coarse_vectors[i][var_name]->clear();
734  }
735  }
736 
737  // Restore old partitioner and renumbering settings
738  mesh.partitioner().reset(old_partitioner.release());
739  mesh.allow_renumbering(old_renumbering_setting);
740 }
741 
742 } // namespace libMesh
743 
744 #endif // #ifdef LIBMESH_ENABLE_AMR
Manages the family, order, etc. parameters for a given FE.
Definition: fe_type.h:179
vectors_iterator vectors_end()
Definition: system.h:2257
void uniformly_p_refine(unsigned int n=1)
Manages multiples systems of equations.
const Elem * parent() const
Definition: elem.h:2479
unsigned int n_systems() const
const unsigned int invalid_uint
Definition: libmesh.h:245
virtual void get(const std::vector< numeric_index_type > &index, T *values) const
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Definition: dof_map.C:1930
unsigned int n_qois() const
Definition: system.h:2278
virtual void disable_cache()
Definition: system.h:2275
const T_sys & get_system(const std::string &name) const
const FEType & variable_type(const unsigned int c) const
Definition: dof_map.h:1792
const EquationSystems & get_equation_systems() const
Definition: system.h:712
The base class for all geometric element types.
Definition: elem.h:100
MeshBase & mesh
IntRange< std::size_t > index_range(const std::vector< T > &vec)
Definition: int_range.h:104
void uniformly_p_coarsen(unsigned int n=1)
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=nullptr, bool estimate_parent_error=false)
const Parallel::Communicator & comm() const
virtual void adjoint_solve(const QoISet &qoi_indices=QoISet())
const unsigned int n_vars
Definition: tecplot_io.C:69
vectors_iterator vectors_begin()
Definition: system.h:2245
DIE A HORRIBLE DEATH HERE typedef float ErrorVectorReal
Base class for Mesh.
Definition: mesh_base.h:77
Number current_solution(const dof_id_type global_dof_number) const
Definition: system.C:194
Manages the degrees of freedom (DOFs) in a simulation.
Definition: dof_map.h:176
std::unique_ptr< QBase > default_quadrature_rule(const unsigned int dim, const int extraorder=0) const
Definition: fe_type.C:31
Responsible for mesh refinement algorithms and data.
dof_id_type id() const
Definition: dof_object.h:655
Manages consistently variables, degrees of freedom, and coefficient vectors.
Definition: system.h:92
FEMNormType type(unsigned int var) const
Definition: system_norm.C:110
void uniformly_coarsen(unsigned int n=1)
const NumericVector< Number > & get_vector(const std::string &vec_name) const
Definition: system.C:774
std::unique_ptr< NumericVector< Number > > solution
Definition: system.h:1523
static std::unique_ptr< FEGenericBase > build(const unsigned int dim, const FEType &type)
virtual std::pair< unsigned int, Real > adjoint_solve(const QoISet &qoi_indices=QoISet())
Definition: system.h:2300
std::map< std::pair< const System *, unsigned int >, ErrorVector * > ErrorMap
void reduce_error(std::vector< ErrorVectorReal > &error_per_cell, const Parallel::Communicator &comm) const
virtual void estimate_error(const System &system, ErrorVector &error_per_cell, const NumericVector< Number > *solution_vector=nullptr, bool estimate_parent_error=false) override
virtual void solve()
Definition: system.h:326
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:748
Real weight_sq(unsigned int var) const
Definition: system_norm.C:175
Real norm_sq() const
Definition: type_vector.h:943
virtual void update()
Definition: system.C:408
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual void swap(NumericVector< T > &v)
const NumericVector< Number > * request_vector(const std::string &vec_name) const
Definition: system.C:716
bool & project_solution_on_reinit(void)
Definition: system.h:794
const MeshBase & get_mesh() const
NumericVector< Number > & get_adjoint_solution(unsigned int i=0)
Definition: system.C:969
std::unique_ptr< NumericVector< Number > > current_local_solution
Definition: system.h:1535
unsigned int n_vars() const
Definition: system.h:2105
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
const DofMap & get_dof_map() const
Definition: system.h:2049
virtual ErrorEstimatorType type() const override
const std::vector< dof_id_type > & get_send_list() const
Definition: dof_map.h:450
Real norm_sq() const
Definition: type_tensor.h:1299
uint8_t dof_id_type
Definition: id_types.h:64
void uniformly_refine(unsigned int n=1)