#include <uniform_refinement_estimator.h>
typedef std::map< std::pair< const System *, unsigned int >, ErrorVector * > ErrorMap
UniformRefinementEstimator ()
~UniformRefinementEstimator ()
virtual void estimate_error (const System &system, ErrorVector &error_per_cell, const NumericVector< Number > *solution_vector=NULL, bool estimate_parent_error=false)
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=NULL, bool estimate_parent_error=false)
virtual void estimate_errors (const EquationSystems &equation_systems, ErrorMap &errors_per_cell, const std::map< const System *, const NumericVector< Number > * > *solution_vectors=NULL, bool estimate_parent_error=false)
unsigned char number_h_refinements
unsigned char number_p_refinements
SystemNorm error_norm
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=NULL, bool estimate_parent_error=false)
void reduce_error (std::vector< float > &error_per_cell) const
This class implements a ``brute force'' error estimator which integrates differences between the current solution and the solution on a uniformly refined (in h and/or p, for an arbitrary number of levels) grid.
Author:
Definition at line 42 of file uniform_refinement_estimator.h.
Definition at line 105 of file error_estimator.h.
Definition at line 49 of file uniform_refinement_estimator.h.
References ErrorEstimator::error_norm, and libMeshEnums::H1.
: number_h_refinements(1),
number_p_refinements(0)
{ error_norm = H1; }
Definition at line 56 of file uniform_refinement_estimator.h.
{}
Definition at line 82 of file uniform_refinement_estimator.C.
References MeshBase::active_local_elements_begin(), MeshBase::active_local_elements_end(), EquationSystems::adjoint_solve(), FEBase::build(), NumericVector< T >::clear(), System::current_local_solution, System::current_solution(), FEType::default_quadrature_rule(), DofMap::dof_indices(), ErrorEstimator::error_norm, AutoPtr< Tp >::get(), System::get_dof_map(), System::get_equation_systems(), EquationSystems::get_mesh(), DofMap::get_send_list(), EquationSystems::get_system(), System::get_vector(), libMeshEnums::H1, libMeshEnums::H1_SEMINORM, libMeshEnums::H2, libMeshEnums::H2_SEMINORM, DofObject::id(), NumericVector< T >::init(), libMeshEnums::L2, libmesh_norm(), MeshBase::max_elem_id(), MeshBase::mesh_dimension(), MeshBase::n_elem(), EquationSystems::n_systems(), System::n_vars(), number_h_refinements, number_p_refinements, Elem::parent(), MeshBase::partitioner(), System::project_solution_on_reinit(), ErrorEstimator::reduce_error(), EquationSystems::reinit(), AutoPtr< Tp >::release(), AutoPtr< Tp >::reset(), libMeshEnums::SERIAL, TypeTensor< T >::size_sq(), TypeVector< T >::size_sq(), System::solution, EquationSystems::solve(), NumericVector< T >::swap(), SystemNorm::type(), MeshRefinement::uniformly_coarsen(), MeshRefinement::uniformly_p_coarsen(), MeshRefinement::uniformly_p_refine(), MeshRefinement::uniformly_refine(), System::update(), DofMap::variable_type(), System::vectors_begin(), System::vectors_end(), and SystemNorm::weight_sq().
Referenced by estimate_error(), and estimate_errors().
{
// Get a vector of the Systems we're going to work on,
// and set up a error_norms map if necessary
std::vector<System *> system_list;
AutoPtr<std::map<const System*, SystemNorm > > error_norms =
AutoPtr<std::map<const System*, SystemNorm > >
(new std::map<const System*, SystemNorm>);
if (_es)
{
libmesh_assert(!_system);
libmesh_assert(_es->n_systems());
_system = &(_es->get_system(0));
libmesh_assert(&(_system->get_equation_systems()) == _es);
libmesh_assert(_es->n_systems());
for (unsigned int i=0; i != _es->n_systems(); ++i)
// We have to break the rules here, because we can't refine a const System
system_list.push_back(const_cast<System *>(&(_es->get_system(i))));
// If we're computing one vector, we need to know how to scale
// each variable's contributions to it.
if (_error_norms)
{
libmesh_assert(!errors_per_cell);
}
else
// If we're computing many vectors, we just need to know which
// variables to skip
{
libmesh_assert (errors_per_cell);
_error_norms = error_norms.get();
for (unsigned int i=0; i!= _es->n_systems(); ++i)
{
const System &sys = _es->get_system(i);
unsigned int n_vars = sys.n_vars();
std::vector<Real> weights(n_vars, 0.0);
for (unsigned int v = 0; v != n_vars; ++v)
{
if (errors_per_cell->find(std::make_pair(&sys, v)) ==
errors_per_cell->end())
continue;
weights[v] = 1.0;
}
(*error_norms)[&sys] =
SystemNorm(std::vector<FEMNormType>(n_vars, error_norm.type(0)),
weights);
}
}
}
else
{
libmesh_assert(_system);
// We have to break the rules here, because we can't refine a const System
system_list.push_back(const_cast<System *>(_system));
libmesh_assert(!_error_norms);
(*error_norms)[_system] = error_norm;
_error_norms = error_norms.get();
}
// An EquationSystems reference will be convenient.
// We have to break the rules here, because we can't refine a const System
EquationSystems& es =
const_cast<EquationSystems &>(_system->get_equation_systems());
// The current mesh
MeshBase& mesh = es.get_mesh();
// The dimensionality of the mesh
const unsigned int dim = mesh.mesh_dimension();
// Resize the error_per_cell vectors to be
// the number of elements, initialize them to 0.
if (error_per_cell)
{
error_per_cell->clear();
error_per_cell->resize (mesh.max_elem_id(), 0.);
}
else
{
libmesh_assert(errors_per_cell);
for (ErrorMap::iterator i = errors_per_cell->begin();
i != errors_per_cell->end(); ++i)
{
ErrorVector *e = i->second;
e->clear();
e->resize(mesh.max_elem_id(), 0.);
}
}
// We'll want to back up all coarse grid vectors
std::vector<std::map<std::string, NumericVector<Number> *> >
coarse_vectors(system_list.size());
std::vector<NumericVector<Number> *>
coarse_solutions(system_list.size());
std::vector<NumericVector<Number> *>
coarse_local_solutions(system_list.size());
// And make copies of projected solutions
std::vector<NumericVector<Number> *>
projected_solutions(system_list.size());
// And we'll need to temporarily change solution projection settings
std::vector<bool> old_projection_settings(system_list.size());
// And it'll be best to avoid any repartitioning
AutoPtr<Partitioner> old_partitioner = mesh.partitioner();
mesh.partitioner().reset(NULL);
for (unsigned int i=0; i != system_list.size(); ++i)
{
System &system = *system_list[i];
// Check for valid error_norms
libmesh_assert (_error_norms->find(&system) !=
_error_norms->end());
// Back up the solution vector
coarse_solutions[i] = system.solution->clone().release();
coarse_local_solutions[i] =
system.current_local_solution->clone().release();
// Back up all other coarse grid vectors
for (System::vectors_iterator vec = system.vectors_begin(); vec !=
system.vectors_end(); ++vec)
{
// The (string) name of this vector
const std::string& var_name = vec->first;
coarse_vectors[i][var_name] = vec->second->clone().release();
}
// Use a non-standard solution vector if necessary
if (solution_vectors &&
solution_vectors->find(&system) != solution_vectors->end() &&
solution_vectors->find(&system)->second &&
solution_vectors->find(&system)->second != system.solution.get())
{
NumericVector<Number>* newsol =
const_cast<NumericVector<Number>*>
(solution_vectors->find(&system)->second);
newsol->swap(*system.solution);
system.update();
}
// Make sure the solution is projected when we refine the mesh
old_projection_settings[i] = system.project_solution_on_reinit();
system.project_solution_on_reinit() = true;
}
// Find the number of coarse mesh elements, to make it possible
// to find correct coarse elem ids later
const unsigned int max_coarse_elem_id = mesh.max_elem_id();
#ifndef NDEBUG
// n_coarse_elem is only used in an assertion later so
// avoid declaring it unless asserts are active.
const unsigned int n_coarse_elem = mesh.n_elem();
#endif
// Uniformly refine the mesh
MeshRefinement mesh_refinement(mesh);
libmesh_assert (number_h_refinements > 0 || number_p_refinements > 0);
// FIXME: this may break if there is more than one System
// on this mesh but estimate_error was still called instead of
// estimate_errors
for (unsigned int i = 0; i != number_h_refinements; ++i)
{
mesh_refinement.uniformly_refine(1);
es.reinit();
}
for (unsigned int i = 0; i != number_p_refinements; ++i)
{
mesh_refinement.uniformly_p_refine(1);
es.reinit();
}
for (unsigned int i=0; i != system_list.size(); ++i)
{
System &system = *system_list[i];
// Copy the projected coarse grid solutions, which will be
// overwritten by solve()
// projected_solutions[i] = system.solution->clone().release();
projected_solutions[i] = NumericVector<Number>::build().release();
projected_solutions[i]->init(system.solution->size(), true, SERIAL);
system.solution->localize(*projected_solutions[i],
system.get_dof_map().get_send_list());
}
// Get the uniformly refined solution.
if (_es)
{
// No specified vectors == forward solve
if (!solution_vectors)
es.solve();
else
{
libmesh_assert(solution_vectors->size() == es.n_systems());
libmesh_assert(solution_vectors->find(system_list[0]) !=
solution_vectors->end());
const bool solve_adjoint =
(system_list[0]->have_vector('adjoint_solution') &&
(solution_vectors->find(system_list[0])->second ==
&system_list[0]->get_adjoint_solution()));
libmesh_assert(solve_adjoint ||
(solution_vectors->find(system_list[0])->second ==
system_list[0]->solution.get()) ||
!solution_vectors->find(system_list[0])->second);
#ifdef DEBUG
for (unsigned int i=0; i != system_list.size(); ++i)
{
libmesh_assert(solution_vectors->find(system_list[i]) !=
solution_vectors->end());
libmesh_assert(!solve_adjoint ||
solution_vectors->find(system_list[i])->second ==
&system_list[i]->get_adjoint_solution());
libmesh_assert(solve_adjoint ||
solution_vectors->find(system_list[i])->second ==
system_list[i]->solution.get() ||
!solution_vectors->find(system_list[i])->second);
}
#endif
if (solve_adjoint)
{
// Set up proper initial guesses
for (unsigned int i=0; i != system_list.size(); ++i)
system_list[i]->get_adjoint_solution() = *system_list[i]->solution;
es.adjoint_solve();
// Put the adjoint_solution into solution for
// comparisons
for (unsigned int i=0; i != system_list.size(); ++i)
{
system_list[i]->get_adjoint_solution().swap(*system_list[i]->solution);
system_list[i]->update();
}
}
else
es.solve();
}
}
else
{
// No specified vectors == forward solve
if (!solution_vectors)
system_list[0]->solve();
else
{
libmesh_assert(solution_vectors->find(system_list[0]) !=
solution_vectors->end());
const bool solve_adjoint =
(system_list[0]->have_vector('adjoint_solution') &&
(solution_vectors->find(system_list[0])->second ==
&system_list[0]->get_adjoint_solution()));
libmesh_assert(solve_adjoint ||
(solution_vectors->find(system_list[0])->second ==
system_list[0]->solution.get()) ||
!solution_vectors->find(system_list[0])->second);
if (solve_adjoint)
{
// Set up proper initial guesses
for (unsigned int i=0; i != system_list.size(); ++i)
system_list[0]->get_adjoint_solution() = *system_list[0]->solution;
system_list[0]->adjoint_solve();
// Put the adjoint_solution into solution for
// comparisons
system_list[0]->get_adjoint_solution().swap(*system_list[0]->solution);
system_list[0]->update();
}
else
system_list[0]->solve();
}
}
// Get the error in the uniformly refined solution(s).
for (unsigned int i=0; i != system_list.size(); ++i)
{
System &system = *system_list[i];
unsigned int n_vars = system.n_vars();
DofMap &dof_map = system.get_dof_map();
const SystemNorm &system_i_norm =
_error_norms->find(&system)->second;
NumericVector<Number> *projected_solution = projected_solutions[i];
// Loop over all the variables in the system
for (unsigned int var=0; var<n_vars; var++)
{
// Get the error vector to fill for this system and variable
ErrorVector *err_vec = error_per_cell;
if (!err_vec)
{
libmesh_assert(errors_per_cell);
err_vec =
(*errors_per_cell)[std::make_pair(&system,var)];
}
// The type of finite element to use for this variable
const FEType& fe_type = dof_map.variable_type (var);
// Finite element object for each fine element
AutoPtr<FEBase> fe (FEBase::build (dim, fe_type));
// Build and attach an appropriate quadrature rule
AutoPtr<QBase> qrule = fe_type.default_quadrature_rule(dim);
fe->attach_quadrature_rule (qrule.get());
const std::vector<Real>& JxW = fe->get_JxW();
const std::vector<std::vector<Real> >& phi = fe->get_phi();
const std::vector<std::vector<RealGradient> >& dphi =
fe->get_dphi();
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
const std::vector<std::vector<RealTensor> >& d2phi =
fe->get_d2phi();
#endif
// The global DOF indices for the fine element
std::vector<unsigned int> dof_indices;
// Iterate over all the active elements in the fine mesh
// that live on this processor.
MeshBase::const_element_iterator elem_it = mesh.active_local_elements_begin();
const MeshBase::const_element_iterator elem_end = mesh.active_local_elements_end();
for (; elem_it != elem_end; ++elem_it)
{
// e is necessarily an active element on the local processor
const Elem* elem = *elem_it;
// Find the element id for the corresponding coarse grid element
const Elem* coarse = elem;
unsigned int e_id = coarse->id();
while (e_id >= max_coarse_elem_id)
{
libmesh_assert (coarse->parent());
coarse = coarse->parent();
e_id = coarse->id();
}
double L2normsq = 0., H1seminormsq = 0., H2seminormsq = 0.;
// reinitialize the element-specific data
// for the current element
fe->reinit (elem);
// Get the local to global degree of freedom maps
dof_map.dof_indices (elem, dof_indices, var);
// The number of quadrature points
const unsigned int n_qp = qrule->n_points();
// The number of shape functions
const unsigned int n_sf = dof_indices.size();
//
// Begin the loop over the Quadrature points.
//
for (unsigned int qp=0; qp<n_qp; qp++)
{
Number u_fine = 0., u_coarse = 0.;
Gradient grad_u_fine, grad_u_coarse;
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
Tensor grad2_u_fine, grad2_u_coarse;
#endif
// Compute solution values at the current
// quadrature point. This reqiures a sum
// over all the shape functions evaluated
// at the quadrature point.
for (unsigned int i=0; i<n_sf; i++)
{
u_fine += phi[i][qp]*system.current_solution (dof_indices[i]);
u_coarse += phi[i][qp]*(*projected_solution) (dof_indices[i]);
grad_u_fine += dphi[i][qp]*system.current_solution (dof_indices[i]);
grad_u_coarse += dphi[i][qp]*(*projected_solution) (dof_indices[i]);
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
grad2_u_fine += d2phi[i][qp]*system.current_solution (dof_indices[i]);
grad2_u_coarse += d2phi[i][qp]*(*projected_solution) (dof_indices[i]);
#endif
}
// Compute the value of the error at this quadrature point
const Number val_error = u_fine - u_coarse;
// Add the squares of the error to each contribution
if (system_i_norm.type(var) == L2 ||
system_i_norm.type(var) == H1 ||
system_i_norm.type(var) == H2)
{
L2normsq += JxW[qp] * system_i_norm.weight_sq(var) *
libmesh_norm(val_error);
libmesh_assert (L2normsq >= 0.);
}
// Compute the value of the error in the gradient at this
// quadrature point
if (system_i_norm.type(var) == H1 ||
system_i_norm.type(var) == H2 ||
system_i_norm.type(var) == H1_SEMINORM)
{
Gradient grad_error = grad_u_fine - grad_u_coarse;
H1seminormsq += JxW[qp] * system_i_norm.weight_sq(var) *
grad_error.size_sq();
libmesh_assert (H1seminormsq >= 0.);
}
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
// Compute the value of the error in the hessian at this
// quadrature point
if (system_i_norm.type(var) == H2 ||
system_i_norm.type(var) == H2_SEMINORM)
{
Tensor grad2_error = grad2_u_fine - grad2_u_coarse;
H2seminormsq += JxW[qp] * system_i_norm.weight_sq(var) *
grad2_error.size_sq();
libmesh_assert (H2seminormsq >= 0.);
}
#endif
} // end qp loop
if (system_i_norm.type(var) == L2 ||
system_i_norm.type(var) == H1 ||
system_i_norm.type(var) == H2)
(*err_vec)[e_id] += L2normsq;
if (system_i_norm.type(var) == H1 ||
system_i_norm.type(var) == H2 ||
system_i_norm.type(var) == H1_SEMINORM)
(*err_vec)[e_id] += H1seminormsq;
if (system_i_norm.type(var) == H2 ||
system_i_norm.type(var) == H2_SEMINORM)
(*err_vec)[e_id] += H2seminormsq;
} // End loop over active local elements
} // End loop over variables
// Don't bother projecting the solution; we'll restore from backup
// after coarsening
system.project_solution_on_reinit() = false;
}
// Uniformly coarsen the mesh, without projecting the solution
libmesh_assert (number_h_refinements > 0 || number_p_refinements > 0);
for (unsigned int i = 0; i != number_h_refinements; ++i)
{
mesh_refinement.uniformly_coarsen(1);
// FIXME - should the reinits here be necessary? - RHS
es.reinit();
}
for (unsigned int i = 0; i != number_p_refinements; ++i)
{
mesh_refinement.uniformly_p_coarsen(1);
es.reinit();
}
// We should be back where we started
libmesh_assert(n_coarse_elem == mesh.n_elem());
// Each processor has now computed the error contribuions
// for its local elements. We need to sum the vector
// and then take the square-root of each component. Note
// that we only need to sum if we are running on multiple
// processors, and we only need to take the square-root
// if the value is nonzero. There will in general be many
// zeros for the inactive elements.
if (error_per_cell)
{
// First sum the vector of estimated error values
this->reduce_error(*error_per_cell);
// Compute the square-root of each component.
START_LOG('std::sqrt()', 'UniformRefinementEstimator');
for (unsigned int i=0; i<error_per_cell->size(); i++)
if ((*error_per_cell)[i] != 0.)
(*error_per_cell)[i] = std::sqrt((*error_per_cell)[i]);
STOP_LOG('std::sqrt()', 'UniformRefinementEstimator');
}
else
{
for (ErrorMap::iterator i = errors_per_cell->begin();
i != errors_per_cell->end(); ++i)
{
ErrorVector *e = i->second;
// First sum the vector of estimated error values
this->reduce_error(*e);
// Compute the square-root of each component.
START_LOG('std::sqrt()', 'UniformRefinementEstimator');
for (unsigned int i=0; i<e->size(); i++)
if ((*e)[i] != 0.)
(*e)[i] = std::sqrt((*e)[i]);
STOP_LOG('std::sqrt()', 'UniformRefinementEstimator');
}
}
// Restore old solutions and clean up the heap
for (unsigned int i=0; i != system_list.size(); ++i)
{
System &system = *system_list[i];
system.project_solution_on_reinit() = old_projection_settings[i];
// Restore the coarse solution vectors and delete their copies
*system.solution = *coarse_solutions[i];
delete coarse_solutions[i];
*system.current_local_solution = *coarse_local_solutions[i];
delete coarse_local_solutions[i];
delete projected_solutions[i];
for (System::vectors_iterator vec = system.vectors_begin(); vec !=
system.vectors_end(); ++vec)
{
// The (string) name of this vector
const std::string& var_name = vec->first;
system.get_vector(var_name) = *coarse_vectors[i][var_name];
coarse_vectors[i][var_name]->clear();
delete coarse_vectors[i][var_name];
}
}
// Restore old partitioner settings
mesh.partitioner() = old_partitioner;
}
system.solve() must be called, and so should have no side effects.
Only the provided system is solved on the refined mesh; for problems decoupled into multiple systems, use of estimate_errors() should be more reliable.
The estimated error is output in the vector error_per_cell
Implements ErrorEstimator.
Definition at line 45 of file uniform_refinement_estimator.C.
References _estimate_error().
{
START_LOG('estimate_error()', 'UniformRefinementEstimator');
std::map<const System*, const NumericVector<Number>* > solution_vectors;
solution_vectors[&_system] = solution_vector;
this->_estimate_error (NULL, &_system, &error_per_cell, NULL, NULL,
&solution_vectors, estimate_parent_error);
STOP_LOG('estimate_error()', 'UniformRefinementEstimator');
}
This function is named estimate_errors instead of estimate_error because otherwise C++ can get confused.
Reimplemented from ErrorEstimator.
Definition at line 58 of file uniform_refinement_estimator.C.
References _estimate_error().
{
START_LOG('estimate_errors()', 'UniformRefinementEstimator');
this->_estimate_error (&_es, NULL, &error_per_cell, NULL,
&error_norms, solution_vectors,
estimate_parent_error);
STOP_LOG('estimate_errors()', 'UniformRefinementEstimator');
}
The user selects which errors get computed by filling a map with error vectors: If errors_per_cell[&system][v] exists, it will be filled with the error values in variable v of system
Reimplemented from ErrorEstimator.
Definition at line 71 of file uniform_refinement_estimator.C.
References _estimate_error().
{
START_LOG('estimate_errors()', 'UniformRefinementEstimator');
this->_estimate_error (&_es, NULL, NULL, &errors_per_cell, NULL,
solution_vectors, estimate_parent_error);
STOP_LOG('estimate_errors()', 'UniformRefinementEstimator');
}
Definition at line 32 of file error_estimator.C.
Referenced by _estimate_error(), PatchRecoveryErrorEstimator::estimate_error(), and JumpErrorEstimator::estimate_error().
{
// This function must be run on all processors at once
parallel_only();
// Each processor has now computed the error contribuions
// for its local elements. We may need to sum the vector to
// recover the error for each element.
Parallel::sum(error_per_cell);
}
Part of this functionality was supported via component_scale and sobolev_order in older libMesh versions, and a small part was supported via component_mask in even older versions. Hopefully the encapsulation here will allow us to avoid changing this API again.
Definition at line 137 of file error_estimator.h.
Referenced by _estimate_error(), KellyErrorEstimator::boundary_side_integration(), DiscontinuityMeasure::boundary_side_integration(), DiscontinuityMeasure::DiscontinuityMeasure(), JumpErrorEstimator::estimate_error(), ErrorEstimator::estimate_errors(), ExactErrorEstimator::ExactErrorEstimator(), ExactErrorEstimator::find_squared_element_error(), KellyErrorEstimator::internal_side_integration(), LaplacianErrorEstimator::internal_side_integration(), DiscontinuityMeasure::internal_side_integration(), KellyErrorEstimator::KellyErrorEstimator(), LaplacianErrorEstimator::LaplacianErrorEstimator(), PatchRecoveryErrorEstimator::EstimateError::operator()(), PatchRecoveryErrorEstimator::PatchRecoveryErrorEstimator(), and UniformRefinementEstimator().
Definition at line 109 of file uniform_refinement_estimator.h.
Referenced by _estimate_error().
Definition at line 114 of file uniform_refinement_estimator.h.
Referenced by _estimate_error().
Generated automatically by Doxygen for libMesh from the source code.