libMesh
uniform_refinement_estimator.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 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"
28 #include "libmesh/equation_systems.h"
29 #include "libmesh/error_vector.h"
30 #include "libmesh/fe.h"
31 #include "libmesh/libmesh_common.h"
32 #include "libmesh/libmesh_logging.h"
33 #include "libmesh/mesh_base.h"
34 #include "libmesh/mesh_refinement.h"
35 #include "libmesh/numeric_vector.h"
36 #include "libmesh/quadrature.h"
37 #include "libmesh/system.h"
38 #include "libmesh/uniform_refinement_estimator.h"
39 #include "libmesh/partitioner.h"
40 #include "libmesh/tensor_tools.h"
41 
42 #ifdef LIBMESH_ENABLE_AMR
43 
44 namespace libMesh
45 {
46 
47 //-----------------------------------------------------------------
48 // ErrorEstimator implementations
50  ErrorVector & error_per_cell,
51  const NumericVector<Number> * solution_vector,
52  bool estimate_parent_error)
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 }
65 
67  ErrorVector & error_per_cell,
68  const std::map<const System *, SystemNorm> & error_norms,
69  const std::map<const System *, const NumericVector<Number> *> * solution_vectors,
70  bool estimate_parent_error)
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 }
81 
83  ErrorMap & errors_per_cell,
84  const std::map<const System *, const NumericVector<Number> *> * solution_vectors,
85  bool estimate_parent_error)
86 {
87  LOG_SCOPE("estimate_errors()", "UniformRefinementEstimator");
88  this->_estimate_error (&_es,
91  &errors_per_cell,
93  solution_vectors,
94  estimate_parent_error);
95 }
96 
98  const System * _system,
99  ErrorVector * error_per_cell,
100  ErrorMap * errors_per_cell,
101  const std::map<const System *, SystemNorm> * _error_norms,
102  const std::map<const System *, const NumericVector<Number> *> * solution_vectors,
103  bool)
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;
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  for (const auto & elem : mesh.active_local_element_ptr_range())
515  {
516  // Find the element id for the corresponding coarse grid element
517  const Elem * coarse = elem;
518  dof_id_type e_id = coarse->id();
519  while (e_id >= max_coarse_elem_id)
520  {
521  libmesh_assert (coarse->parent());
522  coarse = coarse->parent();
523  e_id = coarse->id();
524  }
525 
526  Real L2normsq = 0., H1seminormsq = 0.;
527 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
528  Real H2seminormsq = 0.;
529 #endif
530 
531  // reinitialize the element-specific data
532  // for the current element
533  fe->reinit (elem);
534 
535  // Get the local to global degree of freedom maps
536  dof_map.dof_indices (elem, dof_indices, var);
537 
538  // The number of quadrature points
539  const unsigned int n_qp = qrule->n_points();
540 
541  // The number of shape functions
542  const unsigned int n_sf =
543  cast_int<unsigned int>(dof_indices.size());
544 
545  //
546  // Begin the loop over the Quadrature points.
547  //
548  for (unsigned int qp=0; qp<n_qp; qp++)
549  {
550  Number u_fine = 0., u_coarse = 0.;
551 
552  Gradient grad_u_fine, grad_u_coarse;
553 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
554  Tensor grad2_u_fine, grad2_u_coarse;
555 #endif
556 
557  // Compute solution values at the current
558  // quadrature point. This requires a sum
559  // over all the shape functions evaluated
560  // at the quadrature point.
561  for (unsigned int i=0; i<n_sf; i++)
562  {
563  u_fine += phi[i][qp]*system.current_solution (dof_indices[i]);
564  u_coarse += phi[i][qp]*(*projected_solution) (dof_indices[i]);
565  grad_u_fine += dphi[i][qp]*system.current_solution (dof_indices[i]);
566  grad_u_coarse += dphi[i][qp]*(*projected_solution) (dof_indices[i]);
567 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
568  grad2_u_fine += d2phi[i][qp]*system.current_solution (dof_indices[i]);
569  grad2_u_coarse += d2phi[i][qp]*(*projected_solution) (dof_indices[i]);
570 #endif
571  }
572 
573  // Compute the value of the error at this quadrature point
574  const Number val_error = u_fine - u_coarse;
575 
576  // Add the squares of the error to each contribution
577  if (system_i_norm.type(var) == L2 ||
578  system_i_norm.type(var) == H1 ||
579  system_i_norm.type(var) == H2)
580  {
581  L2normsq += JxW[qp] * system_i_norm.weight_sq(var) *
582  TensorTools::norm_sq(val_error);
583  libmesh_assert_greater_equal (L2normsq, 0.);
584  }
585 
586 
587  // Compute the value of the error in the gradient at this
588  // quadrature point
589  if (system_i_norm.type(var) == H1 ||
590  system_i_norm.type(var) == H2 ||
591  system_i_norm.type(var) == H1_SEMINORM)
592  {
593  Gradient grad_error = grad_u_fine - grad_u_coarse;
594 
595  H1seminormsq += JxW[qp] * system_i_norm.weight_sq(var) *
596  grad_error.norm_sq();
597  libmesh_assert_greater_equal (H1seminormsq, 0.);
598  }
599 
600  // Compute the value of the error in the hessian at this
601  // quadrature point
602  if (system_i_norm.type(var) == H2 ||
603  system_i_norm.type(var) == H2_SEMINORM)
604  {
605 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
606  Tensor grad2_error = grad2_u_fine - grad2_u_coarse;
607 
608  H2seminormsq += JxW[qp] * system_i_norm.weight_sq(var) *
609  grad2_error.norm_sq();
610  libmesh_assert_greater_equal (H2seminormsq, 0.);
611 #else
612  libmesh_error_msg
613  ("libMesh was not configured with --enable-second");
614 #endif
615  }
616  } // end qp loop
617 
618  if (system_i_norm.type(var) == L2 ||
619  system_i_norm.type(var) == H1 ||
620  system_i_norm.type(var) == H2)
621  (*err_vec)[e_id] +=
622  static_cast<ErrorVectorReal>(L2normsq);
623  if (system_i_norm.type(var) == H1 ||
624  system_i_norm.type(var) == H2 ||
625  system_i_norm.type(var) == H1_SEMINORM)
626  (*err_vec)[e_id] +=
627  static_cast<ErrorVectorReal>(H1seminormsq);
628 
629 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
630  if (system_i_norm.type(var) == H2 ||
631  system_i_norm.type(var) == H2_SEMINORM)
632  (*err_vec)[e_id] +=
633  static_cast<ErrorVectorReal>(H2seminormsq);
634 #endif
635  } // End loop over active local elements
636  } // End loop over variables
637 
638  // Don't bother projecting the solution; we'll restore from backup
639  // after coarsening
640  system.project_solution_on_reinit() = false;
641  }
642 
643 
644  // Uniformly coarsen the mesh, without projecting the solution
646 
647  for (unsigned int i = 0; i != number_h_refinements; ++i)
648  {
649  mesh_refinement.uniformly_coarsen(1);
650  // FIXME - should the reinits here be necessary? - RHS
651  es.reinit();
652  }
653 
654  for (unsigned int i = 0; i != number_p_refinements; ++i)
655  {
656  mesh_refinement.uniformly_p_coarsen(1);
657  es.reinit();
658  }
659 
660  // We should be back where we started
661  libmesh_assert_equal_to (n_coarse_elem, mesh.n_elem());
662 
663  // Each processor has now computed the error contributions
664  // for its local elements. We need to sum the vector
665  // and then take the square-root of each component. Note
666  // that we only need to sum if we are running on multiple
667  // processors, and we only need to take the square-root
668  // if the value is nonzero. There will in general be many
669  // zeros for the inactive elements.
670 
671  if (error_per_cell)
672  {
673  // First sum the vector of estimated error values
674  this->reduce_error(*error_per_cell, es.comm());
675 
676  // Compute the square-root of each component.
677  LOG_SCOPE("std::sqrt()", "UniformRefinementEstimator");
678  for (std::size_t i=0; i<error_per_cell->size(); i++)
679  if ((*error_per_cell)[i] != 0.)
680  (*error_per_cell)[i] = std::sqrt((*error_per_cell)[i]);
681  }
682  else
683  {
684  for (ErrorMap::iterator it = errors_per_cell->begin();
685  it != errors_per_cell->end(); ++it)
686  {
687  ErrorVector * e = it->second;
688  // First sum the vector of estimated error values
689  this->reduce_error(*e, es.comm());
690 
691  // Compute the square-root of each component.
692  LOG_SCOPE("std::sqrt()", "UniformRefinementEstimator");
693  for (std::size_t i=0; i<e->size(); i++)
694  if ((*e)[i] != 0.)
695  (*e)[i] = std::sqrt((*e)[i]);
696  }
697  }
698 
699  // Restore old solutions and clean up the heap
700  for (std::size_t i=0; i != system_list.size(); ++i)
701  {
702  System & system = *system_list[i];
703 
704  system.project_solution_on_reinit() = old_projection_settings[i];
705 
706  // Restore the coarse solution vectors and delete their copies
707  *system.solution = *coarse_solutions[i];
708  delete coarse_solutions[i];
709  *system.current_local_solution = *coarse_local_solutions[i];
710  delete coarse_local_solutions[i];
711  delete projected_solutions[i];
712 
713  for (System::vectors_iterator vec = system.vectors_begin(); vec !=
714  system.vectors_end(); ++vec)
715  {
716  // The (string) name of this vector
717  const std::string & var_name = vec->first;
718 
719  system.get_vector(var_name) = *coarse_vectors[i][var_name];
720 
721  coarse_vectors[i][var_name]->clear();
722  delete coarse_vectors[i][var_name];
723  }
724  }
725 
726  // Restore old partitioner settings
727  mesh.partitioner().reset(old_partitioner.release());
728 }
729 
730 } // namespace libMesh
731 
732 #endif // #ifdef LIBMESH_ENABLE_AMR
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
vectors_iterator vectors_end()
End of vectors container.
Definition: system.h:2238
unsigned char number_h_refinements
How many h refinements to perform to get the fine grid.
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
Currently this function ignores the error_norm member variable, and uses the function argument error_...
void uniformly_p_refine(unsigned int n=1)
Uniformly p refines the mesh n times.
This is the EquationSystems class.
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
This function does uniform refinements and a solve to get an improved solution on each cell...
const unsigned int invalid_uint
A number which is used quite often to represent an invalid or uninitialized value.
Definition: libmesh.h:184
SystemNorm error_norm
When estimating the error in a single system, the error_norm is used to control the scaling and norm ...
const FEType & variable_type(const unsigned int c) const
Definition: dof_map.h:1697
Real norm_sq() const
Definition: type_tensor.h:1270
unsigned int dim
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:2256
ImplicitSystem & sys
const Elem * parent() const
Definition: elem.h:2346
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
MeshBase & mesh
UniquePtr< NumericVector< Number > > current_local_solution
All the values I need to compute my contribution to the simulation at hand.
Definition: system.h:1535
void uniformly_p_coarsen(unsigned int n=1)
Attempts to uniformly p coarsen the mesh n times.
const class libmesh_nullptr_t libmesh_nullptr
This class defines a norm/seminorm to be applied to a NumericVector which contains coefficients in a ...
Definition: system_norm.h:43
virtual void adjoint_solve(const QoISet &qoi_indices=QoISet())
Call adjoint_solve on all the individual equation systems.
const unsigned int n_vars
Definition: tecplot_io.C:68
The libMesh namespace provides an interface to certain functionality in the library.
vectors_iterator vectors_begin()
Beginning of vectors container.
Definition: system.h:2226
DIE A HORRIBLE DEATH HERE typedef float ErrorVectorReal
This is the MeshBase class.
Definition: mesh_base.h:68
libmesh_assert(j)
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)
The code for estimate_error and both estimate_errors versions is very similar, so we use the same fun...
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
static UniquePtr< NumericVector< T > > build(const Parallel::Communicator &comm, const SolverPackage solver_package=libMesh::default_solver_package())
Builds a NumericVector on the processors in communicator comm using the linear solver package specifi...
virtual SimpleRange< element_iterator > active_local_element_ptr_range()=0
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:167
const NumericVector< Number > * request_vector(const std::string &vec_name) const
Definition: system.C:736
virtual dof_id_type max_elem_id() const =0
const DofMap & get_dof_map() const
Definition: system.h:2030
This is the MeshRefinement class.
std::vector< Number > qoi
Values of the quantities of interest.
Definition: system.h:1553
This is the base class for classes which contain information related to any physical process that mig...
Definition: system.h:76
virtual void reinit()
Reinitialize all the systems.
void uniformly_coarsen(unsigned int n=1)
Attempts to uniformly coarsen the mesh n times.
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 NULL...
Definition: system.h:2275
std::map< std::pair< const System *, unsigned int >, ErrorVector * > ErrorMap
When calculating many error vectors at once, we need a data structure to hold them all...
UniquePtr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1523
Real norm_sq() const
Definition: type_vector.h:940
const NumericVector< Number > & get_vector(const std::string &vec_name) const
Definition: system.C:794
void reduce_error(std::vector< ErrorVectorReal > &error_per_cell, const Parallel::Communicator &comm LIBMESH_CAN_DEFAULT_TO_COMMWORLD) const
This method takes the local error contributions in error_per_cell from each processor and combines th...
unsigned char number_p_refinements
How many p refinements to perform to get the fine grid.
virtual void solve()
Solves the system.
Definition: system.h:310
std::map< std::string, NumericVector< Number > * >::iterator vectors_iterator
Vector iterator typedefs.
Definition: system.h:748
unsigned int n_systems() const
virtual void solve()
Call solve on all the individual equation systems.
Real weight_sq(unsigned int var) const
Definition: system_norm.h:338
virtual void update()
Update the local values to reflect the solution on neighboring processors.
Definition: system.C:425
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Number current_solution(const dof_id_type global_dof_number) const
Definition: system.C:192
virtual UniquePtr< Partitioner > & partitioner()
A partitioner to use at each prepare_for_use()
Definition: mesh_base.h:112
virtual void swap(NumericVector< T > &v)
Swaps the contents of this with v.
UniquePtr< QBase > default_quadrature_rule(const unsigned int dim, const int extraorder=0) const
Definition: fe_type.C:30
const Parallel::Communicator & comm() const
unsigned int mesh_dimension() const
Definition: mesh_base.C:148
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:794
NumericVector< Number > & get_adjoint_solution(unsigned int i=0)
Definition: system.C:989
const EquationSystems & get_equation_systems() const
Definition: system.h:712
const MeshBase & get_mesh() const
virtual void clear()
Restores the NumericVector<T> to a pristine state.
const T_sys & get_system(const std::string &name) const
unsigned int n_vars() const
Definition: system.h:2086
FEMNormType type(unsigned int var) const
Definition: system_norm.h:268
dof_id_type id() const
Definition: dof_object.h:632
virtual dof_id_type n_elem() const =0
const std::vector< dof_id_type > & get_send_list() const
Definition: dof_map.h:394
This class defines a tensor in LIBMESH_DIM dimensional Real or Complex space.
uint8_t dof_id_type
Definition: id_types.h:64
void uniformly_refine(unsigned int n=1)
Uniformly refines the mesh n times.
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Fills the vector di with the global degree of freedom indices for the element.
Definition: dof_map.C:1917
static UniquePtr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.