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