20 #include "libmesh/dof_map.h" 21 #include "libmesh/elem.h" 22 #include "libmesh/equation_systems.h" 23 #include "libmesh/error_vector.h" 24 #include "libmesh/fe.h" 25 #include "libmesh/libmesh_common.h" 26 #include "libmesh/libmesh_logging.h" 27 #include "libmesh/mesh_base.h" 28 #include "libmesh/mesh_refinement.h" 29 #include "libmesh/numeric_vector.h" 30 #include "libmesh/quadrature.h" 31 #include "libmesh/system.h" 32 #include "libmesh/uniform_refinement_estimator.h" 33 #include "libmesh/partitioner.h" 34 #include "libmesh/tensor_tools.h" 35 #include "libmesh/enum_error_estimator_type.h" 36 #include "libmesh/enum_norm_type.h" 37 #include "libmesh/int_range.h" 46 #ifdef LIBMESH_ENABLE_AMR 56 number_h_refinements(1),
57 number_p_refinements(0),
74 bool estimate_parent_error)
76 LOG_SCOPE(
"estimate_error()",
"UniformRefinementEstimator");
77 std::map<const System *, const NumericVector<Number> *> solution_vectors;
78 solution_vectors[&_system] = solution_vector;
85 estimate_parent_error);
90 const std::map<const System *, SystemNorm> & error_norms,
92 bool estimate_parent_error)
94 LOG_SCOPE(
"estimate_errors()",
"UniformRefinementEstimator");
101 estimate_parent_error);
107 bool estimate_parent_error)
109 LOG_SCOPE(
"estimate_errors()",
"UniformRefinementEstimator");
116 estimate_parent_error);
123 const std::map<const System *, SystemNorm> * _error_norms,
129 std::vector<System *> system_list;
130 auto error_norms = std::make_unique<std::map<const System *, SystemNorm>>();
137 libmesh_assert_equal_to (&(_system->get_equation_systems()), _es);
142 system_list.push_back(const_cast<System *>(&(_es->
get_system(i))));
156 _error_norms = error_norms.get();
163 std::vector<Real> weights(
n_vars, 0.0);
164 for (
unsigned int v = 0; v !=
n_vars; ++v)
166 if (!errors_per_cell->count(std::make_pair(&sys, v)))
171 (*error_norms)[&sys] =
181 system_list.push_back(const_cast<System *>(_system));
185 _error_norms = error_norms.get();
197 const unsigned int dim =
mesh.mesh_dimension();
203 error_per_cell->clear();
204 error_per_cell->resize (
mesh.max_elem_id(), 0.);
209 for (
const auto & pr : *errors_per_cell)
213 e->resize(
mesh.max_elem_id(), 0.);
218 std::vector<std::map<std::string, std::unique_ptr<NumericVector<Number>>>> coarse_vectors(system_list.size());
219 std::vector<std::unique_ptr<NumericVector<Number>>> coarse_solutions(system_list.size());
220 std::vector<std::unique_ptr<NumericVector<Number>>> coarse_local_solutions(system_list.size());
222 std::vector<std::unique_ptr<NumericVector<Number>>> projected_solutions(system_list.size());
225 std::vector<bool> old_projection_settings(system_list.size());
228 std::unique_ptr<Partitioner> old_partitioner(
mesh.partitioner().release());
231 const bool old_renumbering_setting =
mesh.allow_renumbering();
232 mesh.allow_renumbering(
false);
236 System & system = *system_list[i];
240 _error_norms->end());
243 coarse_solutions[i] = system.
solution->clone();
251 const std::string & var_name = vec->first;
253 coarse_vectors[i][var_name] = vec->second->clone();
257 if (solution_vectors &&
258 solution_vectors->find(&system) != solution_vectors->end() &&
259 solution_vectors->find(&system)->second &&
260 solution_vectors->find(&system)->second != system.
solution.get())
264 (solution_vectors->find(&system)->second);
305 System & system = *system_list[i];
310 projected_solutions[i]->init(system.
solution->size(),
314 system.
solution->localize(*projected_solutions[i],
319 bool solve_adjoint =
false;
320 if (solution_vectors)
322 System * sys = system_list[0];
324 solution_vectors->end());
328 std::ostringstream adjoint_name;
329 adjoint_name <<
"adjoint_solution" << j;
333 solve_adjoint =
true;
349 if (!solution_vectors)
353 libmesh_assert_equal_to (solution_vectors->size(), es.
n_systems());
355 solution_vectors->end());
357 (solution_vectors->find(system_list[0])->second ==
358 system_list[0]->solution.get()) ||
359 !solution_vectors->find(system_list[0])->second);
362 for (
const auto & sys : system_list)
365 solution_vectors->end());
369 bool found_vec =
false;
372 std::ostringstream adjoint_name;
373 adjoint_name <<
"adjoint_solution" << j;
375 if (vec == sys->request_vector(adjoint_name.str()))
390 std::vector<unsigned int> adjs(system_list.size(),
395 System * sys = system_list[i];
397 solution_vectors->end());
401 std::ostringstream adjoint_name;
402 adjoint_name <<
"adjoint_solution" << j;
420 system_list[i]->get_adjoint_solution(adjs[i]).swap(*system_list[i]->solution);
421 system_list[i]->update();
430 System * sys = system_list[0];
437 if (!solution_vectors)
442 solution_vectors->end());
447 (solution_vectors->find(sys)->second ==
449 !solution_vectors->find(sys)->second);
454 for (
unsigned int j=0, n_qois = sys->
n_qois();
457 std::ostringstream adjoint_name;
458 adjoint_name <<
"adjoint_solution" << j;
484 System & system = *system_list[sysnum];
491 _error_norms->find(&system)->second;
496 for (
unsigned int var=0; var<
n_vars; var++)
504 (*errors_per_cell)[std::make_pair(&system,var)].get();
515 fe->attach_quadrature_rule (qrule.get());
517 const std::vector<Real> & JxW = fe->get_JxW();
518 const std::vector<std::vector<Real>> & phi = fe->get_phi();
519 const std::vector<std::vector<RealGradient>> & dphi =
521 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 522 const std::vector<std::vector<RealTensor>> & d2phi =
527 std::vector<dof_id_type> dof_indices;
531 for (
const auto & elem :
mesh.active_local_element_ptr_range())
534 const Elem * coarse = elem;
536 while (e_id >= max_coarse_elem_id)
539 coarse = coarse->
parent();
543 Real L2normsq = 0., H1seminormsq = 0.;
544 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 545 Real H2seminormsq = 0.;
556 const unsigned int n_qp = qrule->n_points();
559 const unsigned int n_sf =
560 cast_int<unsigned int>(dof_indices.size());
565 for (
unsigned int qp=0; qp<n_qp; qp++)
567 Number u_fine = 0., u_coarse = 0.;
569 Gradient grad_u_fine, grad_u_coarse;
570 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 571 Tensor grad2_u_fine, grad2_u_coarse;
578 for (
unsigned int i=0; i<n_sf; i++)
581 u_coarse += phi[i][qp]*(*projected_solution) (dof_indices[i]);
583 grad_u_coarse += dphi[i][qp]*(*projected_solution) (dof_indices[i]);
584 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 586 grad2_u_coarse += d2phi[i][qp]*(*projected_solution) (dof_indices[i]);
591 const Number val_error = u_fine - u_coarse;
594 if (system_i_norm.
type(var) ==
L2 ||
595 system_i_norm.
type(var) ==
H1 ||
596 system_i_norm.
type(var) ==
H2)
598 L2normsq += JxW[qp] * system_i_norm.
weight_sq(var) *
600 libmesh_assert_greater_equal (L2normsq, 0.);
606 if (system_i_norm.
type(var) ==
H1 ||
607 system_i_norm.
type(var) ==
H2 ||
610 Gradient grad_error = grad_u_fine - grad_u_coarse;
612 H1seminormsq += JxW[qp] * system_i_norm.
weight_sq(var) *
614 libmesh_assert_greater_equal (H1seminormsq, 0.);
619 if (system_i_norm.
type(var) ==
H2 ||
622 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 623 Tensor grad2_error = grad2_u_fine - grad2_u_coarse;
625 H2seminormsq += JxW[qp] * system_i_norm.
weight_sq(var) *
627 libmesh_assert_greater_equal (H2seminormsq, 0.);
630 (
"libMesh was not configured with --enable-second");
635 if (system_i_norm.
type(var) ==
L2 ||
636 system_i_norm.
type(var) ==
H1 ||
637 system_i_norm.
type(var) ==
H2)
640 if (system_i_norm.
type(var) ==
H1 ||
641 system_i_norm.
type(var) ==
H2 ||
646 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 647 if (system_i_norm.
type(var) ==
H2 ||
678 libmesh_assert_equal_to (n_coarse_elem,
mesh.n_elem());
694 LOG_SCOPE(
"std::sqrt()",
"UniformRefinementEstimator");
695 for (
auto & val : *error_per_cell)
697 val = std::sqrt(val);
701 for (
const auto & pr : *errors_per_cell)
708 LOG_SCOPE(
"std::sqrt()",
"UniformRefinementEstimator");
711 val = std::sqrt(val);
718 System & system = *system_list[i];
723 *system.
solution = *coarse_solutions[i];
730 const std::string & var_name = vec->first;
732 system.
get_vector(var_name) = *coarse_vectors[i][var_name];
734 coarse_vectors[i][var_name]->
clear();
739 mesh.partitioner().reset(old_partitioner.release());
740 mesh.allow_renumbering(old_renumbering_setting);
745 #endif // #ifdef LIBMESH_ENABLE_AMR class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
std::map< std::string, std::unique_ptr< NumericVector< Number > >, std::less<> >::iterator vectors_iterator
Vector iterator typedefs.
vectors_iterator vectors_end()
End of vectors container.
void uniformly_p_refine(unsigned int n=1)
Uniformly p refines the mesh n times.
This is the EquationSystems class.
const Elem * parent() const
unsigned int n_systems() const
const unsigned int invalid_uint
A number which is used quite often to represent an invalid or uninitialized value for an unsigned int...
SystemNorm error_norm
When estimating the error in a single system, the error_norm is used to control the scaling and norm ...
virtual void get(const std::vector< numeric_index_type > &index, T *values) const
Access multiple components at once.
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
unsigned int n_qois() const
Number of currently active quantities of interest.
The ErrorVector is a specialization of the StatisticsVector for error data computed on a finite eleme...
virtual void disable_cache()
Avoids use of any cached data that might affect any solve result.
ErrorEstimatorType
Defines an enum for the different types of error estimators which are available.
const FEType & variable_type(const unsigned int c) const
const EquationSystems & get_equation_systems() const
This is the base class from which all geometric element types are derived.
auto norm_sq() const -> decltype(std::norm(T()))
void uniformly_p_coarsen(unsigned int n=1)
Attempts to uniformly p coarsen the mesh n times.
const Parallel::Communicator & comm() const
This class defines a norm/seminorm to be applied to a NumericVector which contains coefficients in a ...
This class defines a vector in LIBMESH_DIM dimensional Real or Complex space.
virtual void adjoint_solve(const QoISet &qoi_indices=QoISet())
Call adjoint_solve on all the individual equation systems.
The libMesh namespace provides an interface to certain functionality in the library.
vectors_iterator vectors_begin()
Beginning of vectors container.
std::map< std::pair< const System *, unsigned int >, std::unique_ptr< ErrorVector > > ErrorMap
When calculating many error vectors at once, we need a data structure to hold them all...
const T_sys & get_system(std::string_view name) const
DIE A HORRIBLE DEATH HERE typedef float ErrorVectorReal
This is the MeshBase class.
Number current_solution(const dof_id_type global_dof_number) const
This class handles the numbering of degrees of freedom on a mesh.
std::unique_ptr< QBase > default_quadrature_rule(const unsigned int dim, const int extraorder=0) const
auto norm_sq() const -> decltype(std::norm(T()))
Implements (adaptive) mesh refinement algorithms for a MeshBase.
Manages consistently variables, degrees of freedom, and coefficient vectors.
virtual void reinit()
Handle any mesh changes and reinitialize all the systems on the updated mesh.
FEMNormType type(unsigned int var) const
void uniformly_coarsen(unsigned int n=1)
Attempts to uniformly coarsen the mesh n times.
std::unique_ptr< NumericVector< Number > > solution
Data structure to hold solution values.
static std::unique_ptr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.
This class holds functions that will estimate the error in a finite element solution on a given mesh...
virtual std::pair< unsigned int, Real > adjoint_solve(const QoISet &qoi_indices=QoISet())
Solves the adjoint system, for the specified qoi indices, or for every qoi if qoi_indices is nullptr...
void reduce_error(std::vector< ErrorVectorReal > &error_per_cell, const Parallel::Communicator &comm) const
This method takes the local error contributions in error_per_cell from each processor and combines th...
virtual void solve()
Solves the system.
const NumericVector< Number > * request_vector(std::string_view vec_name) const
virtual void solve()
Call solve on all the individual equation systems.
Real weight_sq(unsigned int var) const
virtual void update()
Update the local values to reflect the solution on neighboring processors.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual void swap(NumericVector< T > &v)
Swaps the contents of this with v.
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...
const MeshBase & get_mesh() const
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
static std::unique_ptr< NumericVector< T > > build(const Parallel::Communicator &comm, SolverPackage solver_package=libMesh::default_solver_package(), ParallelType parallel_type=AUTOMATIC)
Builds a NumericVector on the processors in communicator comm using the linear solver package specifi...
NumericVector< Number > & get_adjoint_solution(unsigned int i=0)
std::unique_ptr< NumericVector< Number > > current_local_solution
All the values I need to compute my contribution to the simulation at hand.
virtual void clear()
Restores the NumericVector<T> to a pristine state.
unsigned int n_vars() const
const DofMap & get_dof_map() const
const std::vector< dof_id_type > & get_send_list() const
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
const NumericVector< Number > & get_vector(std::string_view vec_name) const
This class defines a tensor in LIBMESH_DIM dimensional Real or Complex space.
void uniformly_refine(unsigned int n=1)
Uniformly refines the mesh n times.