adjoint_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 <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
21 #include <cmath> // for sqrt
22 #include <set>
23 #include <sstream> // for ostringstream
24 #include <unordered_map>
25 #include <unordered_set>
26 
27 // Local Includes
28 #include "libmesh/dof_map.h"
29 #include "libmesh/elem.h"
31 #include "libmesh/error_vector.h"
32 #include "libmesh/fe.h"
33 #include "libmesh/fe_interface.h"
34 #include "libmesh/libmesh_common.h"
36 #include "libmesh/mesh_base.h"
38 #include "libmesh/numeric_vector.h"
39 #include "libmesh/quadrature.h"
40 #include "libmesh/system.h"
41 #include "libmesh/diff_physics.h"
42 #include "libmesh/fem_system.h"
44 #include "libmesh/partitioner.h"
47 #include "libmesh/enum_norm_type.h"
48 #include "libmesh/utility.h"
49 
50 #ifdef LIBMESH_ENABLE_AMR
51 
52 namespace libMesh
53 {
54 
55 //-----------------------------------------------------------------
56 // AdjointRefinementErrorEstimator
57 
58 // As of 10/2/2012, this function implements a 'brute-force' adjoint based QoI
59 // error estimator, using the following relationship:
60 // Q(u) - Q(u_h) \approx - R( (u_h)_(h/2), z_(h/2) ) .
61 // where: Q(u) is the true QoI
62 // u_h is the approximate primal solution on the current FE space
63 // Q(u_h) is the approximate QoI
64 // z_(h/2) is the adjoint corresponding to Q, on a richer FE space
65 // (u_h)_(h/2) is a projection of the primal solution on the richer FE space
66 // By richer FE space, we mean a grid that has been refined once and a polynomial order
67 // that has been increased once, i.e. one h and one p refinement
68 
69 // Both a global QoI error estimate and element wise error indicators are included
70 // Note that the element wise error indicators slightly over estimate the error in
71 // each element
72 
75  number_h_refinements(1),
76  number_p_refinements(0),
77  _residual_evaluation_physics(nullptr),
78  _qoi_set(QoISet())
79 {
80  // We're not actually going to use error_norm; our norms are
81  // absolute values of QoI error.
83 }
84 
86 {
87  return ADJOINT_REFINEMENT;
88 }
89 
91  ErrorVector & error_per_cell,
92  const NumericVector<Number> * solution_vector,
93  bool /*estimate_parent_error*/)
94 {
95  // We have to break the rules here, because we can't refine a const System
96  System & system = const_cast<System &>(_system);
97 
98  // An EquationSystems reference will be convenient.
99  EquationSystems & es = system.get_equation_systems();
100 
101  // The current mesh
102  MeshBase & mesh = es.get_mesh();
103 
104  // Get coarse grid adjoint solutions. This should be a relatively
105  // quick (especially with preconditioner reuse) way to get a good
106  // initial guess for the fine grid adjoint solutions. More
107  // importantly, subtracting off a coarse adjoint approximation gives
108  // us better local error indication, and subtracting off *some* lift
109  // function is necessary for correctness if we have heterogeneous
110  // adjoint Dirichlet conditions.
111 
112  // Solve the adjoint problem(s) on the coarse FE space
113  // Only if the user didn't already solve it for us
114  if (!system.is_adjoint_already_solved())
115  system.adjoint_solve(_qoi_set);
116 
117  // Loop over all the adjoint problems and, if any have heterogenous
118  // Dirichlet conditions, get the corresponding coarse lift
119  // function(s)
120  for (unsigned int j=0,
121  n_qois = system.n_qois(); j != n_qois; j++)
122  {
123  // Skip this QoI if it is not in the QoI Set or if there are no
124  // heterogeneous Dirichlet boundaries for it
125  if (_qoi_set.has_index(j) &&
127  {
128  // Next, we are going to build up the residual for evaluating the flux QoI
129  NumericVector<Number> * coarse_residual = nullptr;
130 
131  // The definition of a flux QoI is R(u^h, L) where R is the residual as defined
132  // by a conservation law. Therefore, if we are using stabilization, the
133  // R should be specified by the user via the residual_evaluation_physics
134 
135  // If the residual physics pointer is not null, use it when
136  // evaluating here.
137  {
138  const bool swapping_physics = _residual_evaluation_physics;
139  if (swapping_physics)
140  dynamic_cast<DifferentiableSystem &>(system).swap_physics(_residual_evaluation_physics);
141 
142  // Assemble without applying constraints, to capture the solution values on the boundary
143  (dynamic_cast<ImplicitSystem &>(system)).assembly(true, false, false, true);
144 
145  // Get the residual vector (no constraints applied on boundary, so we wont blow away the lift)
146  coarse_residual = &(dynamic_cast<ExplicitSystem &>(system)).get_vector("RHS Vector");
147  coarse_residual->close();
148 
149  // Now build the lift function and add it to the system vectors
150  std::ostringstream liftfunc_name;
151  liftfunc_name << "adjoint_lift_function" << j;
152 
153  system.add_vector(liftfunc_name.str());
154 
155  // Initialize lift with coarse adjoint solve associate with this flux QoI to begin with
156  system.get_vector(liftfunc_name.str()).init(system.get_adjoint_solution(j), false);
157 
158  // Build the actual lift using adjoint dirichlet conditions
160  (system.get_vector(liftfunc_name.str()), static_cast<unsigned int>(j));
161 
162  // Compute the flux R(u^h, L)
163  std::cout<<"The flux QoI "<<static_cast<unsigned int>(j)<<" is: "<<coarse_residual->dot(system.get_vector(liftfunc_name.str()))<<std::endl<<std::endl;
164 
165  // Swap back if needed
166  if (swapping_physics)
167  dynamic_cast<DifferentiableSystem &>(system).swap_physics(_residual_evaluation_physics);
168  }
169  } // End if QoI in set and flux/dirichlet boundary QoI
170  } // End loop over QoIs
171 
172  // We'll want to back up all coarse grid vectors
173  std::map<std::string, std::unique_ptr<NumericVector<Number>>> coarse_vectors;
174  for (const auto & pr : as_range(system.vectors_begin(), system.vectors_end()))
175  {
176  // The (string) name of this vector
177  const std::string & var_name = pr.first;
178 
179  coarse_vectors[var_name] = pr.second->clone();
180  }
181 
182  // Back up the coarse solution and coarse local solution
183  std::unique_ptr<NumericVector<Number>> coarse_solution = system.solution->clone();
184  std::unique_ptr<NumericVector<Number>> coarse_local_solution = system.current_local_solution->clone();
185 
186  // And we'll need to temporarily change solution projection settings
187  bool old_projection_setting;
188  old_projection_setting = system.project_solution_on_reinit();
189 
190  // Make sure the solution is projected when we refine the mesh
191  system.project_solution_on_reinit() = true;
192 
193  // And it'll be best to avoid any repartitioning
194  std::unique_ptr<Partitioner> old_partitioner(mesh.partitioner().release());
195 
196  // And we can't allow any renumbering
197  const bool old_renumbering_setting = mesh.allow_renumbering();
198  mesh.allow_renumbering(false);
199 
200  // Use a non-standard solution vector if necessary
201  if (solution_vector && solution_vector != system.solution.get())
202  {
203  NumericVector<Number> * newsol =
204  const_cast<NumericVector<Number> *> (solution_vector);
205  newsol->swap(*system.solution);
206  system.update();
207  }
208 
209  // Resize the error_per_cell vector to be
210  // the number of elements, initialized to 0.
211  error_per_cell.clear();
212  error_per_cell.resize (mesh.max_elem_id(), 0.);
213 
214 #ifndef NDEBUG
215  // These variables are only used in assertions later so
216  // avoid declaring them unless asserts are active.
217  const dof_id_type n_coarse_elem = mesh.n_active_elem();
218 
219  dof_id_type local_dof_bearing_nodes = 0;
220  const unsigned int sysnum = system.number();
221  for (const auto * node : mesh.local_node_ptr_range())
222  for (unsigned int v=0, nvars = node->n_vars(sysnum);
223  v != nvars; ++v)
224  if (node->n_comp(sysnum, v))
225  {
226  local_dof_bearing_nodes++;
227  break;
228  }
229 #endif // NDEBUG
230 
231  // Uniformly refine the mesh
232  MeshRefinement mesh_refinement(mesh);
233 
234  libmesh_assert (number_h_refinements > 0 || number_p_refinements > 0);
235 
236  // FIXME: this may break if there is more than one System
237  // on this mesh but estimate_error was still called instead of
238  // estimate_errors
239  for (unsigned int i = 0; i != number_h_refinements; ++i)
240  {
241  mesh_refinement.uniformly_refine(1);
242  es.reinit();
243  }
244 
245  for (unsigned int i = 0; i != number_p_refinements; ++i)
246  {
247  mesh_refinement.uniformly_p_refine(1);
248  es.reinit();
249  }
250 
251  // Copy the projected coarse grid solutions, which will be
252  // overwritten by solve()
253  std::vector<std::unique_ptr<NumericVector<Number>>> coarse_adjoints;
254  for (unsigned int j=0; j != system.n_qois(); j++)
255  {
256  if (_qoi_set.has_index(j))
257  {
258  auto coarse_adjoint = NumericVector<Number>::build(mesh.comm());
259 
260  // Can do "fast" init since we're overwriting this in a sec
261  coarse_adjoint->init(system.get_adjoint_solution(j),
262  /* fast = */ true);
263 
264  *coarse_adjoint = system.get_adjoint_solution(j);
265 
266  coarse_adjoints.emplace_back(std::move(coarse_adjoint));
267  }
268  else
269  coarse_adjoints.emplace_back(nullptr);
270  }
271 
272  // Next, we are going to build up the residual for evaluating the
273  // error estimate.
274 
275  // If the residual physics pointer is not null, use it when
276  // evaluating here.
277  {
278  const bool swapping_physics = _residual_evaluation_physics;
279  if (swapping_physics)
280  dynamic_cast<DifferentiableSystem &>(system).swap_physics(_residual_evaluation_physics);
281 
282  // Rebuild the rhs with the projected primal solution, constraints
283  // have to be applied to get the correct error estimate since error
284  // on the Dirichlet boundary is zero
285  (dynamic_cast<ImplicitSystem &>(system)).assembly(true, false);
286 
287  // Swap back if needed
288  if (swapping_physics)
289  dynamic_cast<DifferentiableSystem &>(system).swap_physics(_residual_evaluation_physics);
290  }
291 
292  NumericVector<Number> & projected_residual = (dynamic_cast<ExplicitSystem &>(system)).get_vector("RHS Vector");
293  projected_residual.close();
294 
295  // Solve the adjoint problem(s) on the refined FE space
296  system.adjoint_solve(_qoi_set);
297 
298  // Now that we have the refined adjoint solution and the projected primal solution,
299  // we first compute the global QoI error estimate
300 
301  // Resize the computed_global_QoI_errors vector to hold the error estimates for each QoI
302  computed_global_QoI_errors.resize(system.n_qois());
303 
304  // Loop over all the adjoint solutions and get the QoI error
305  // contributions from all of them. While we're looping anyway we'll
306  // pull off the coarse adjoints
307  for (unsigned int j=0; j != system.n_qois(); j++)
308  {
309  // Skip this QoI if not in the QoI Set
310  if (_qoi_set.has_index(j))
311  {
312 
313  // If the adjoint solution has heterogeneous dirichlet
314  // values, then to get a proper error estimate here we need
315  // to subtract off a coarse grid lift function.
316  // |Q(u) - Q(u^h)| = |R([u^h]+, z^h+ - [L]+)| + HOT
318  {
319  // Need to create a string with current loop index to retrieve
320  // the correct vector from the liftvectors map
321  std::ostringstream liftfunc_name;
322  liftfunc_name << "adjoint_lift_function" << j;
323 
324  // Subtract off the corresponding lift vector from the adjoint solution
325  system.get_adjoint_solution(j) -= system.get_vector(liftfunc_name.str());
326 
327  // Now evaluate R(u^h, z^h+ - lift)
328  computed_global_QoI_errors[j] = projected_residual.dot(system.get_adjoint_solution(j));
329 
330  // Add the lift back to get back the adjoint
331  system.get_adjoint_solution(j) += system.get_vector(liftfunc_name.str());
332 
333  }
334  else
335  {
336  // Usual dual weighted residual error estimate
337  // |Q(u) - Q(u^h)| = |R([u^h]+, z^h+)| + HOT
338  computed_global_QoI_errors[j] = projected_residual.dot(system.get_adjoint_solution(j));
339  }
340 
341  }
342  }
343 
344  // Done with the global error estimates, now construct the element wise error indicators
345 
346  // To get a better element wise breakdown of the error estimate,
347  // we subtract off a coarse representation of the adjoint solution.
348  // |Q(u) - Q(u^h)| = |R([u^h]+, z^h+ - [z^h]+)|
349  // This remains valid for all combinations of heterogenous adjoint bcs and
350  // stabilized/non-stabilized formulations, except for the case where we not using a
351  // heterogenous adjoint bc and have a stabilized formulation.
352  // Then, R(u^h_s, z^h_s) != 0 (no Galerkin orthogonality w.r.t the non-stabilized residual)
353  for (unsigned int j=0; j != system.n_qois(); j++)
354  {
355  // Skip this QoI if not in the QoI Set
356  if (_qoi_set.has_index(j))
357  {
358  // If we have a nullptr residual evaluation physics pointer, we
359  // assume the user's formulation is consistent from mesh to
360  // mesh, so we have Galerkin orthogonality and we can get
361  // better indicator results by subtracting a coarse adjoint.
362 
363  // If we have a residual evaluation physics pointer, but we
364  // also have heterogeneous adjoint dirichlet boundaries,
365  // then we have to subtract off *some* lift function for
366  // consistency, and we choose the coarse adjoint in lieu of
367  // having any better ideas.
368 
369  // If we have a residual evaluation physics pointer and we
370  // have homogeneous adjoint dirichlet boundaries, then we
371  // don't have to subtract off anything, and with stabilized
372  // formulations we get the best results if we don't.
375  {
376  // z^h+ -> z^h+ - [z^h]+
377  system.get_adjoint_solution(j) -= *coarse_adjoints[j];
378  }
379  }
380  }
381 
382  // We ought to account for 'spill-over' effects while computing the
383  // element error indicators This happens because the same dof is
384  // shared by multiple elements, one way of mitigating this is to
385  // scale the contribution from each dof by the number of elements it
386  // belongs to We first obtain the number of elements each node
387  // belongs to
388 
389  // A map that relates a node id to an int that will tell us how many elements it is a node of
390  std::unordered_map<dof_id_type, unsigned int> shared_element_count;
391 
392  // To fill this map, we will loop over elements, and then in each element, we will loop
393  // over the nodes each element contains, and then query it for the number of coarse
394  // grid elements it was a node of
395 
396  // Keep track of which nodes we have already dealt with
397  std::unordered_set<dof_id_type> processed_node_ids;
398 
399  // We will be iterating over all the active elements in the fine mesh that live on
400  // this processor.
401  for (const auto & elem : mesh.active_local_element_ptr_range())
402  for (unsigned int n=0; n != elem->n_nodes(); ++n)
403  {
404  // Get a reference to the current node
405  const Node & node = elem->node_ref(n);
406 
407  // Get the id of this node
408  dof_id_type node_id = node.id();
409 
410  // If we havent already processed this node, do so now
411  if (processed_node_ids.find(node_id) == processed_node_ids.end())
412  {
413  // Declare a neighbor_set to be filled by the find_point_neighbors
414  std::set<const Elem *> fine_grid_neighbor_set;
415 
416  // Call find_point_neighbors to fill the neighbor_set
417  elem->find_point_neighbors(node, fine_grid_neighbor_set);
418 
419  // A vector to hold the coarse grid parents neighbors
420  std::vector<dof_id_type> coarse_grid_neighbors;
421 
422  // Loop over all the fine neighbors of this node
423  for (const auto & fine_elem : fine_grid_neighbor_set)
424  {
425  // Find the element id for the corresponding coarse grid element
426  const Elem * coarse_elem = fine_elem;
427  for (unsigned int j = 0; j != number_h_refinements; ++j)
428  {
429  libmesh_assert (coarse_elem->parent());
430 
431  coarse_elem = coarse_elem->parent();
432  }
433 
434  // Loop over the existing coarse neighbors and check if this one is
435  // already in there
436  const dof_id_type coarse_id = coarse_elem->id();
437  std::size_t j = 0;
438 
439  // If the set already contains this element break out of the loop
440  for (; j != coarse_grid_neighbors.size(); j++)
441  if (coarse_grid_neighbors[j] == coarse_id)
442  break;
443 
444  // If we didn't leave the loop even at the last element,
445  // this is a new neighbour, put in the coarse_grid_neighbor_set
446  if (j == coarse_grid_neighbors.size())
447  coarse_grid_neighbors.push_back(coarse_id);
448  } // End loop over fine neighbors
449 
450  // Set the shared_neighbour index for this node to the
451  // size of the coarse grid neighbor set
452  shared_element_count[node_id] =
453  cast_int<unsigned int>(coarse_grid_neighbors.size());
454 
455  // Add this node to processed_node_ids vector
456  processed_node_ids.insert(node_id);
457  } // End if not processed node
458  } // End loop over nodes
459 
460  // Get a DoF map, will be used to get the nodal dof_indices for each element
461  DofMap & dof_map = system.get_dof_map();
462 
463  // The global DOF indices, we will use these later on when we compute the element wise indicators
464  std::vector<dof_id_type> dof_indices;
465 
466  // Localize the global rhs and adjoint solution vectors (which might be shared on multiple processors) onto a
467  // local ghosted vector, this ensures each processor has all the dof_indices to compute an error indicator for
468  // an element it owns
469  std::unique_ptr<NumericVector<Number>> localized_projected_residual = NumericVector<Number>::build(system.comm());
470  localized_projected_residual->init(system.n_dofs(), system.n_local_dofs(), system.get_dof_map().get_send_list(), false, GHOSTED);
471  projected_residual.localize(*localized_projected_residual, system.get_dof_map().get_send_list());
472 
473  // Each adjoint solution will also require ghosting; for efficiency we'll reuse the same memory
474  std::unique_ptr<NumericVector<Number>> localized_adjoint_solution = NumericVector<Number>::build(system.comm());
475  localized_adjoint_solution->init(system.n_dofs(), system.n_local_dofs(), system.get_dof_map().get_send_list(), false, GHOSTED);
476 
477  // We will loop over each adjoint solution, localize that adjoint
478  // solution and then loop over local elements
479  for (unsigned int i=0; i != system.n_qois(); i++)
480  {
481  // Skip this QoI if not in the QoI Set
482  if (_qoi_set.has_index(i))
483  {
484  // Get the weight for the current QoI
485  Real error_weight = _qoi_set.weight(i);
486 
487  (system.get_adjoint_solution(i)).localize(*localized_adjoint_solution, system.get_dof_map().get_send_list());
488 
489  // Loop over elements
490  for (const auto & elem : mesh.active_local_element_ptr_range())
491  {
492  // Go up number_h_refinements levels up to find the coarse parent
493  const Elem * coarse = elem;
494 
495  for (unsigned int j = 0; j != number_h_refinements; ++j)
496  {
497  libmesh_assert (coarse->parent());
498 
499  coarse = coarse->parent();
500  }
501 
502  const dof_id_type e_id = coarse->id();
503 
504  // Get the local to global degree of freedom maps for this element
505  dof_map.dof_indices (elem, dof_indices);
506 
507  // We will have to manually do the dot products.
508  Number local_contribution = 0.;
509 
510  // Sum the contribution to the error indicator for each element from the current QoI
511  for (const auto & dof : dof_indices)
512  local_contribution += (*localized_projected_residual)(dof) * (*localized_adjoint_solution)(dof);
513 
514  // Multiply by the error weight for this QoI
515  local_contribution *= error_weight;
516 
517  // FIXME: we're throwing away information in the
518  // --enable-complex case
519  error_per_cell[e_id] += static_cast<ErrorVectorReal>
520  (std::abs(local_contribution));
521 
522  } // End loop over elements
523  } // End if belong to QoI set
524  } // End loop over QoIs
525 
526  // Don't bother projecting the solution; we'll restore from backup
527  // after coarsening
528  system.project_solution_on_reinit() = false;
529 
530  // Uniformly coarsen the mesh, without projecting the solution
531  libmesh_assert (number_h_refinements > 0 || number_p_refinements > 0);
532 
533  for (unsigned int i = 0; i != number_h_refinements; ++i)
534  {
535  mesh_refinement.uniformly_coarsen(1);
536  // FIXME - should the reinits here be necessary? - RHS
537  es.reinit();
538  }
539 
540  for (unsigned int i = 0; i != number_p_refinements; ++i)
541  {
542  mesh_refinement.uniformly_p_coarsen(1);
543  es.reinit();
544  }
545 
546  // We should have the same number of active elements as when we started,
547  // but not necessarily the same number of elements since reinit() doesn't
548  // always call contract()
549  libmesh_assert_equal_to (n_coarse_elem, mesh.n_active_elem());
550 
551  // We should have the same number of dof-bearing nodes as when we
552  // started
553 #ifndef NDEBUG
554  dof_id_type final_local_dof_bearing_nodes = 0;
555  for (const auto * node : mesh.local_node_ptr_range())
556  for (unsigned int v=0, nvars = node->n_vars(sysnum);
557  v != nvars; ++v)
558  if (node->n_comp(sysnum, v))
559  {
560  final_local_dof_bearing_nodes++;
561  break;
562  }
563  libmesh_assert_equal_to (local_dof_bearing_nodes,
564  final_local_dof_bearing_nodes);
565 #endif // NDEBUG
566 
567  // Restore old solutions and clean up the heap
568  system.project_solution_on_reinit() = old_projection_setting;
569 
570  // Restore the coarse solution vectors and delete their copies
571  *system.solution = *coarse_solution;
572  *system.current_local_solution = *coarse_local_solution;
573 
574  for (const auto & pr : as_range(system.vectors_begin(), system.vectors_end()))
575  {
576  // The (string) name of this vector
577  const std::string & var_name = pr.first;
578 
579  // If it's a vector we already had (and not a newly created
580  // vector like an adjoint rhs), we need to restore it.
581  auto it = coarse_vectors.find(var_name);
582  if (it != coarse_vectors.end())
583  {
584  NumericVector<Number> * coarsevec = it->second.get();
585  system.get_vector(var_name) = *coarsevec;
586 
587  coarsevec->clear();
588  }
589  }
590 
591  // Restore old partitioner and renumbering settings
592  mesh.partitioner().reset(old_partitioner.release());
593  mesh.allow_renumbering(old_renumbering_setting);
594 
595  // Finally sum the vector of estimated error values.
596  this->reduce_error(error_per_cell, system.comm());
597 
598  // We don't take a square root here; this is a goal-oriented
599  // estimate not a Hilbert norm estimate.
600 } // end estimate_error function
601 
602 }// namespace libMesh
603 
604 #endif // #ifdef LIBMESH_ENABLE_AMR
vectors_iterator vectors_end()
Definition: system.h:2257
virtual ErrorEstimatorType type() const
void uniformly_p_refine(unsigned int n=1)
double abs(double a)
Manages multiples systems of equations.
const Elem * parent() const
Definition: elem.h:2479
A geometric point in (x,y,z) space associated with a DOF.
Definition: node.h:52
virtual void assembly(bool get_residual, bool get_jacobian, bool apply_heterogeneous_constraints=false, bool apply_no_constraints=false) override=0
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
Used to specify quantities of interest in a simulation.
Definition: qoi_set.h:45
unsigned int n_qois() const
Definition: system.h:2278
const EquationSystems & get_equation_systems() const
Definition: system.h:712
The base class for all geometric element types.
Definition: elem.h:100
MeshBase & mesh
void uniformly_p_coarsen(unsigned int n=1)
const Parallel::Communicator & comm() const
vectors_iterator vectors_begin()
Definition: system.h:2245
virtual T dot(const NumericVector< T > &v) const =0
dof_id_type n_local_dofs() const
Definition: system.C:187
DIE A HORRIBLE DEATH HERE typedef float ErrorVectorReal
dof_id_type n_dofs() const
Definition: system.C:150
Base class for Mesh.
Definition: mesh_base.h:77
void enforce_adjoint_constraints_exactly(NumericVector< Number > &v, unsigned int q) const
Manages the degrees of freedom (DOFs) in a simulation.
Definition: dof_map.h:176
bool has_index(std::size_t) const
Definition: qoi_set.h:221
unsigned int number() const
Definition: system.h:2025
NumericVector< Number > & add_vector(const std::string &vec_name, const bool projections=true, const ParallelType type=PARALLEL)
Definition: system.C:661
Responsible for mesh refinement algorithms and data.
bool has_adjoint_dirichlet_boundaries(unsigned int q) const
dof_id_type id() const
Definition: dof_object.h:655
Manages consistently variables, degrees of freedom, and coefficient vectors.
Definition: system.h:92
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
void init(triangulateio &t)
SimpleRange< I > as_range(const std::pair< I, I > &p)
Definition: simple_range.h:57
virtual std::pair< unsigned int, Real > adjoint_solve(const QoISet &qoi_indices=QoISet())
Definition: system.h:2300
void reduce_error(std::vector< ErrorVectorReal > &error_per_cell, const Parallel::Communicator &comm) const
bool is_adjoint_already_solved() const
Definition: system.h:388
virtual void close()=0
static std::unique_ptr< NumericVector< T > > build(const Parallel::Communicator &comm, const SolverPackage solver_package=libMesh::default_solver_package())
Real weight(std::size_t) const
Definition: qoi_set.h:240
virtual void update()
Definition: system.C:408
virtual void estimate_error(const System &system, ErrorVector &error_per_cell, const NumericVector< Number > *solution_vector=nullptr, bool estimate_parent_error=false)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual void swap(NumericVector< T > &v)
virtual std::unique_ptr< DifferentiableQoI > clone() override
Definition: diff_system.h:169
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
const DofMap & get_dof_map() const
Definition: system.h:2049
const std::vector< dof_id_type > & get_send_list() const
Definition: dof_map.h:450
Manages consistently variables, degrees of freedom, and coefficient vectors for explicit systems...
uint8_t dof_id_type
Definition: id_types.h:64
void uniformly_refine(unsigned int n=1)
Manages consistently variables, degrees of freedom, coefficient vectors, and matrices for implicit sy...
virtual void localize(std::vector< T > &v_local) const =0