#include <exact_error_estimator.h>
typedef std::map< std::pair< const System *, unsigned int >, ErrorVector * > ErrorMap
ExactErrorEstimator ()
~ExactErrorEstimator ()
void attach_exact_value (Number fptr(const Point &p, const Parameters &Parameters, const std::string &sys_name, const std::string &unknown_name))
void attach_exact_deriv (Gradient fptr(const Point &p, const Parameters ¶meters, const std::string &sys_name, const std::string &unknown_name))
void attach_exact_hessian (Tensor fptr(const Point &p, const Parameters ¶meters, const std::string &sys_name, const std::string &unknown_name))
void attach_reference_solution (EquationSystems *es_fine)
void extra_quadrature_order (const int extraorder)
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)
void reduce_error (std::vector< float > &error_per_cell) const
Real find_squared_element_error (const System &system, const std::string &var_name, const Elem *elem, const DenseVector< Number > &Uelem, FEBase *fe, MeshFunction *fine_values) const
Number(* _exact_value )(const Point &p, const Parameters ¶meters, const std::string &sys_name, const std::string &unknown_name)
Gradient(* _exact_deriv )(const Point &p, const Parameters ¶meters, const std::string &sys_name, const std::string &unknown_name)
Tensor(* _exact_hessian )(const Point &p, const Parameters ¶meters, const std::string &sys_name, const std::string &unknown_name)
EquationSystems * _equation_systems_fine
int _extra_order
This class implements an 'error estimator' based on the difference between the approximate and exact solution. In theory the quadrature error in this estimate should be much lower than the approximation error in other estimates, so this estimator can be used to calculate effectivity.
Author:
Definition at line 62 of file exact_error_estimator.h.
Definition at line 105 of file error_estimator.h.
Definition at line 70 of file exact_error_estimator.h.
References ErrorEstimator::error_norm, and libMeshEnums::H1.
: _exact_value(NULL),
_exact_deriv(NULL),
_exact_hessian(NULL),
_extra_order(0)
{ error_norm = H1; }
Definition at line 79 of file exact_error_estimator.h.
{}
Definition at line 58 of file exact_error_estimator.C.
References _equation_systems_fine, and _exact_deriv.
{
libmesh_assert (fptr != NULL);
_exact_deriv = fptr;
// If we're not using a fine grid solution
_equation_systems_fine = NULL;
}
Definition at line 72 of file exact_error_estimator.C.
References _equation_systems_fine, and _exact_hessian.
{
libmesh_assert (fptr != NULL);
_exact_hessian = fptr;
// If we're not using a fine grid solution
_equation_systems_fine = NULL;
}
Definition at line 84 of file exact_error_estimator.C.
References _equation_systems_fine, _exact_deriv, _exact_hessian, and _exact_value.
{
libmesh_assert (es_fine != NULL);
_equation_systems_fine = es_fine;
// If we're using a fine grid solution, we're not using exact values
_exact_value = NULL;
_exact_deriv = NULL;
_exact_hessian = NULL;
}
Implements ErrorEstimator.
Definition at line 96 of file exact_error_estimator.C.
References Elem::active(), MeshBase::active_local_elements_begin(), MeshBase::active_local_elements_end(), NumericVector< T >::build(), FEBase::build(), Elem::child(), FEBase::coarsened_dof_values(), System::current_local_solution, System::current_solution(), FEType::default_quadrature_rule(), DofMap::dof_indices(), AutoPtr< Tp >::get(), System::get_dof_map(), System::get_mesh(), DofObject::id(), MeshBase::max_elem_id(), MeshBase::mesh_dimension(), Elem::n_children(), System::n_vars(), System::name(), Elem::parent(), libMeshEnums::SERIAL, System::solution, NumericVector< T >::swap(), System::update_global_solution(), System::variable_name(), System::variable_number(), and DofMap::variable_type().
{
// The current mesh
const MeshBase& mesh = system.get_mesh();
// The dimensionality of the mesh
const unsigned int dim = mesh.mesh_dimension();
// The number of variables in the system
const unsigned int n_vars = system.n_vars();
// The DofMap for this system
const DofMap& dof_map = system.get_dof_map();
// Resize the error_per_cell vector to be
// the number of elements, initialize it to 0.
error_per_cell.resize (mesh.max_elem_id());
std::fill (error_per_cell.begin(), error_per_cell.end(), 0.);
// Prepare current_local_solution to localize a non-standard
// solution vector if necessary
if (solution_vector && solution_vector != system.solution.get())
{
NumericVector<Number>* newsol =
const_cast<NumericVector<Number>*>(solution_vector);
System &sys = const_cast<System&>(system);
newsol->swap(*sys.solution);
sys.update();
}
// Loop over all the variables in the system
for (unsigned int var=0; var<n_vars; var++)
{
// Possibly skip this variable
if (error_norm.weight(var) == 0.0) continue;
// The (string) name of this variable
const std::string& var_name = system.variable_name(var);
// The type of finite element to use for this variable
const FEType& fe_type = dof_map.variable_type (var);
AutoPtr<FEBase> fe (FEBase::build (dim, fe_type));
// Build an appropriate Gaussian quadrature rule
AutoPtr<QBase> qrule =
fe_type.default_quadrature_rule (dim,
_extra_order);
fe->attach_quadrature_rule (qrule.get());
// Prepare a global solution and a MeshFunction of the fine system if we need one
AutoPtr<MeshFunction> fine_values;
AutoPtr<NumericVector<Number> > fine_soln = NumericVector<Number>::build();
if (_equation_systems_fine)
{
const System& fine_system = _equation_systems_fine->get_system(system.name());
std::vector<Number> global_soln;
// FIXME - we're assuming that the fine system solution gets
// used even when a different vector is used for the coarse
// system
fine_system.update_global_solution(global_soln);
fine_soln->init (global_soln.size(), true, SERIAL);
(*fine_soln) = global_soln;
fine_values = AutoPtr<MeshFunction>
(new MeshFunction(*_equation_systems_fine,
*fine_soln,
fine_system.get_dof_map(),
fine_system.variable_number(var_name)));
fine_values->init();
}
// Request the data we'll need to compute with
fe->get_JxW();
fe->get_phi();
fe->get_dphi();
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
fe->get_d2phi();
#endif
fe->get_xyz();
// Iterate over all the active elements in the 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;
const unsigned int e_id = elem->id();
#ifdef LIBMESH_ENABLE_AMR
// See if the parent of element e has been examined yet;
// if not, we may want to compute the estimator on it
const Elem* parent = elem->parent();
// We only can compute and only need to compute on
// parents with all active children
bool compute_on_parent = true;
if (!parent || !estimate_parent_error)
compute_on_parent = false;
else
for (unsigned int c=0; c != parent->n_children(); ++c)
if (!parent->child(c)->active())
compute_on_parent = false;
if (compute_on_parent &&
!error_per_cell[parent->id()])
{
// Compute a projection onto the parent
DenseVector<Number> Uparent;
FEBase::coarsened_dof_values(*(system.current_local_solution),
dof_map, parent, Uparent,
var, false);
error_per_cell[parent->id()] =
find_squared_element_error(system, var_name, parent, Uparent,
fe.get(), fine_values.get());
}
#endif
// Get the local to global degree of freedom maps
std::vector<unsigned int> dof_indices;
dof_map.dof_indices (elem, dof_indices, var);
DenseVector<Number> Uelem(dof_indices.size());
for (unsigned int i=0; i != dof_indices.size(); ++i)
Uelem(i) = system.current_solution(dof_indices[i]);
error_per_cell[e_id] =
find_squared_element_error(system, var_name, elem, Uelem,
fe.get(), fine_values.get());
} // End loop over active local elements
} // End loop over variables
// 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.
// 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()', 'ExactErrorEstimator');
for (unsigned int i=0; i<error_per_cell.size(); i++)
{
if (error_per_cell[i] != 0.)
{
libmesh_assert (error_per_cell[i] > 0.);
error_per_cell[i] = std::sqrt(error_per_cell[i]);
}
}
STOP_LOG('std::sqrt()', 'ExactErrorEstimator');
// If we used a non-standard solution before, now is the time to fix
// the current_local_solution
if (solution_vector && solution_vector != system.solution.get())
{
NumericVector<Number>* newsol =
const_cast<NumericVector<Number>*>(solution_vector);
System &sys = const_cast<System&>(system);
newsol->swap(*sys.solution);
sys.update();
}
}
Currently this function ignores the error_norm member variable, and uses the function argument error_norms instead.
This function is named estimate_errors instead of estimate_error because otherwise C++ can get confused.
Reimplemented in UniformRefinementEstimator.
Definition at line 46 of file error_estimator.C.
References ErrorEstimator::error_norm, ErrorEstimator::estimate_error(), EquationSystems::get_system(), and EquationSystems::n_systems().
{
SystemNorm old_error_norm = this->error_norm;
// Sum the error values from each system
for (unsigned int s = 0; s != equation_systems.n_systems(); ++s)
{
ErrorVector system_error_per_cell;
const System &sys = equation_systems.get_system(s);
if (error_norms.find(&sys) == error_norms.end())
this->error_norm = old_error_norm;
else
this->error_norm = error_norms.find(&sys)->second;
const NumericVector<Number>* solution_vector = NULL;
if (solution_vectors &&
solution_vectors->find(&sys) != solution_vectors->end())
solution_vector = solution_vectors->find(&sys)->second;
this->estimate_error(sys, system_error_per_cell,
solution_vector, estimate_parent_error);
if (s)
{
libmesh_assert(error_per_cell.size() == system_error_per_cell.size());
for (unsigned int i=0; i != error_per_cell.size(); ++i)
error_per_cell[i] += system_error_per_cell[i];
}
else
error_per_cell = system_error_per_cell;
}
// Restore our old state before returning
this->error_norm = old_error_norm;
}
Currently this function ignores the error_norm.weight() values because it calculates each variable's error individually, unscaled.
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
FIXME: This is a default implementation - derived classes should reimplement it for efficiency.
Reimplemented in UniformRefinementEstimator.
Definition at line 92 of file error_estimator.C.
References ErrorEstimator::error_norm, ErrorEstimator::estimate_error(), EquationSystems::get_system(), EquationSystems::n_systems(), System::n_vars(), and SystemNorm::type().
{
SystemNorm old_error_norm = this->error_norm;
// Find the requested error values from each system
for (unsigned int s = 0; s != equation_systems.n_systems(); ++s)
{
const System &sys = equation_systems.get_system(s);
unsigned int n_vars = sys.n_vars();
for (unsigned int v = 0; v != n_vars; ++v)
{
// Only fill in ErrorVectors the user asks for
if (errors_per_cell.find(std::make_pair(&sys, v)) ==
errors_per_cell.end())
continue;
// Calculate error in only one variable
std::vector<Real> weights(n_vars, 0.0);
weights[v] = 1.0;
this->error_norm =
SystemNorm(std::vector<FEMNormType>(n_vars, old_error_norm.type(0)),
weights);
const NumericVector<Number>* solution_vector = NULL;
if (solution_vectors &&
solution_vectors->find(&sys) != solution_vectors->end())
solution_vector = solution_vectors->find(&sys)->second;
this->estimate_error
(sys, *errors_per_cell[std::make_pair(&sys, v)],
solution_vector, estimate_parent_error);
}
}
// Restore our old state before returning
this->error_norm = old_error_norm;
}
Definition at line 127 of file exact_error_estimator.h.
References _extra_order.
{ _extra_order = extraorder; }
Definition at line 287 of file exact_error_estimator.C.
References _equation_systems_fine, _exact_deriv, _exact_hessian, _exact_value, ErrorEstimator::error_norm, FEBase::get_d2phi(), FEBase::get_dphi(), System::get_equation_systems(), FEBase::get_JxW(), FEBase::get_phi(), FEBase::get_xyz(), MeshFunction::gradient(), libMeshEnums::H1, libMeshEnums::H1_SEMINORM, libMeshEnums::H2, libMeshEnums::H2_SEMINORM, MeshFunction::hessian(), libMeshEnums::L2, libmesh_norm(), System::name(), EquationSystems::parameters, FEBase::reinit(), DenseVector< T >::size(), TypeTensor< T >::size_sq(), TypeVector< T >::size_sq(), SystemNorm::type(), and System::variable_number().
{
// The (string) name of this system
const std::string& sys_name = system.name();
unsigned int var = system.variable_number(var_name);
const Parameters& parameters = system.get_equation_systems().parameters;
// reinitialize the element-specific data
// for the current element
fe->reinit (elem);
// Get the data we need to compute with
const std::vector<Real> & JxW = fe->get_JxW();
const std::vector<std::vector<Real> >& phi_values = fe->get_phi();
const std::vector<std::vector<RealGradient> >& dphi_values = fe->get_dphi();
const std::vector<Point>& q_point = fe->get_xyz();
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
const std::vector<std::vector<RealTensor> >& d2phi_values = fe->get_d2phi();
#endif
// The number of shape functions
const unsigned int n_sf = Uelem.size();
// The number of quadrature points
const unsigned int n_qp = JxW.size();
Real error_val = 0;
// Begin the loop over the Quadrature points.
//
for (unsigned int qp=0; qp<n_qp; qp++)
{
// Real u_h = 0.;
// RealGradient grad_u_h;
Number u_h = 0.;
Gradient grad_u_h;
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
Tensor grad2_u_h;
#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++)
{
// Values from current solution.
u_h += phi_values[i][qp]*Uelem(i);
grad_u_h += dphi_values[i][qp]*Uelem(i);
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
grad2_u_h += d2phi_values[i][qp]*Uelem(i);
#endif
}
// Compute the value of the error at this quadrature point
if (error_norm.type(var) == L2 ||
error_norm.type(var) == H1 ||
error_norm.type(var) == H2)
{
Number val_error = 0;
if(_exact_value)
val_error = u_h - _exact_value(q_point[qp],parameters,sys_name,var_name);
else if(_equation_systems_fine)
val_error = u_h - (*fine_values)(q_point[qp]);
// Add the squares of the error to each contribution
error_val += JxW[qp]*libmesh_norm(val_error);
}
// Compute the value of the error in the gradient at this
// quadrature point
if ((_exact_deriv || _equation_systems_fine) &&
(error_norm.type(var) == H1 ||
error_norm.type(var) == H1_SEMINORM ||
error_norm.type(var) == H2))
{
Gradient grad_error;
if(_exact_deriv)
grad_error = grad_u_h - _exact_deriv(q_point[qp],parameters,sys_name,var_name);
else if(_equation_systems_fine)
grad_error = grad_u_h - fine_values->gradient(q_point[qp]);
error_val += JxW[qp]*grad_error.size_sq();
}
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
// Compute the value of the error in the hessian at this
// quadrature point
if ((_exact_hessian || _equation_systems_fine) &&
(error_norm.type(var) == H2_SEMINORM ||
error_norm.type(var) == H2))
{
Tensor grad2_error;
if(_exact_hessian)
grad2_error = grad2_u_h - _exact_hessian(q_point[qp],parameters,sys_name,var_name);
else if (_equation_systems_fine)
grad2_error = grad2_u_h - fine_values->hessian(q_point[qp]);
error_val += JxW[qp]*grad2_error.size_sq();
}
#endif
} // end qp loop
libmesh_assert (error_val >= 0.);
return error_val;
}
Definition at line 32 of file error_estimator.C.
Referenced by UniformRefinementEstimator::_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);
}
Definition at line 182 of file exact_error_estimator.h.
Referenced by attach_exact_deriv(), attach_exact_hessian(), ExactSolution::attach_exact_value(), attach_reference_solution(), and find_squared_element_error().
Definition at line 164 of file exact_error_estimator.h.
Referenced by attach_exact_deriv(), attach_reference_solution(), and find_squared_element_error().
Definition at line 173 of file exact_error_estimator.h.
Referenced by attach_exact_hessian(), attach_reference_solution(), and find_squared_element_error().
Definition at line 155 of file exact_error_estimator.h.
Referenced by ExactSolution::attach_exact_value(), attach_reference_solution(), and find_squared_element_error().
Definition at line 197 of file exact_error_estimator.h.
Referenced by extra_quadrature_order().
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 UniformRefinementEstimator::_estimate_error(), KellyErrorEstimator::boundary_side_integration(), DiscontinuityMeasure::boundary_side_integration(), DiscontinuityMeasure::DiscontinuityMeasure(), JumpErrorEstimator::estimate_error(), ErrorEstimator::estimate_errors(), 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::UniformRefinementEstimator().
Generated automatically by Doxygen for libMesh from the source code.