#include <jump_error_estimator.h>
Inherits ErrorEstimator.
Inherited by DiscontinuityMeasure, KellyErrorEstimator, and LaplacianErrorEstimator.
typedef std::map< std::pair< const System *, unsigned int >, ErrorVector * > ErrorMap
JumpErrorEstimator ()
virtual ~JumpErrorEstimator ()
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)
bool scale_by_n_flux_faces
SystemNorm error_norm
void reinit_sides ()
float coarse_n_flux_faces_increment ()
virtual void initialize (const System &system, ErrorVector &error_per_cell, bool estimate_parent_error)
virtual void internal_side_integration ()=0
virtual bool boundary_side_integration ()
void reduce_error (std::vector< float > &error_per_cell) const
bool integrate_boundary_sides
const Elem * fine_elem
const Elem * coarse_elem
Real fine_error
Real coarse_error
unsigned int fine_side
unsigned int var
DenseVector< Number > Ufine
DenseVector< Number > Ucoarse
AutoPtr< FEBase > fe_fine
AutoPtr< FEBase > fe_coarse
This abstract base class implements utility functions for error estimators which are based on integrated jumps between elements.
Author:
Definition at line 49 of file jump_error_estimator.h.
Definition at line 105 of file error_estimator.h.
Definition at line 56 of file jump_error_estimator.h.
: scale_by_n_flux_faces(false),
integrate_boundary_sides(false),
fe_fine(NULL), fe_coarse(NULL) {}
Definition at line 64 of file jump_error_estimator.h.
{}
Reimplemented in DiscontinuityMeasure, and KellyErrorEstimator.
Definition at line 119 of file jump_error_estimator.h.
Referenced by estimate_error().
{ return false; }
Definition at line 433 of file jump_error_estimator.C.
References coarse_elem, Elem::dim(), FEBase::dim, fine_elem, and Elem::level().
Referenced by estimate_error().
{
// Keep track of the number of internal flux sides found on each
// element
unsigned int dim = coarse_elem->dim();
const unsigned int divisor =
1 << (dim-1)*(fine_elem->level() - coarse_elem->level());
// With a difference of n levels between fine and coarse elements,
// we compute a fractional flux face for the coarse element by adding:
// 1/2^n in 2D
// 1/4^n in 3D
// each time. This code will get hit 2^n times in 2D and 4^n
// times in 3D so that the final flux face count for the coarse
// element will be an integer value.
return 1.0 / static_cast<Real>(divisor);
}
Implements ErrorEstimator.
Definition at line 52 of file jump_error_estimator.C.
References Elem::active(), MeshBase::active_local_elements_begin(), MeshBase::active_local_elements_end(), boundary_side_integration(), FEBase::build(), Elem::child(), coarse_elem, coarse_error, coarse_n_flux_faces_increment(), FEBase::coarsened_dof_values(), System::current_solution(), FEType::default_quadrature_order(), FEBase::dim, DofMap::dof_indices(), ErrorEstimator::error_norm, fe_coarse, fe_fine, FEBase::fe_type, fine_elem, fine_error, fine_side, AutoPtr< Tp >::get(), System::get_dof_map(), System::get_mesh(), DofObject::id(), initialize(), integrate_boundary_sides, internal_side_integration(), Elem::level(), MeshBase::max_elem_id(), MeshBase::mesh_dimension(), Elem::n_children(), Elem::n_neighbors(), System::n_vars(), Elem::neighbor(), Elem::parent(), FEBase::qrule, ErrorEstimator::reduce_error(), reinit_sides(), DenseVector< T >::resize(), scale_by_n_flux_faces, System::solution, NumericVector< T >::swap(), Ucoarse, Ufine, var, DofMap::variable_type(), and SystemNorm::weight().
Referenced by Adaptive< T >::solve().
{
START_LOG('estimate_error()', 'JumpErrorEstimator');
/*
Conventions for assigning the direction of the normal:
- e & f are global element ids
Case (1.) Elements are at the same level, e<f
Compute the flux jump on the face and
add it as a contribution to error_per_cell[e]
and error_per_cell[f]
----------------------
| | |
| | f |
| | |
| e |---> n |
| | |
| | |
----------------------
Case (2.) The neighbor is at a higher level.
Compute the flux jump on e's face and
add it as a contribution to error_per_cell[e]
and error_per_cell[f]
----------------------
| | | |
| | e |---> n |
| | | |
|-----------| f |
| | | |
| | | |
| | | |
----------------------
*/
// 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.);
// Declare a vector of floats which is as long as
// error_per_cell above, and fill with zeros. This vector will be
// used to keep track of the number of edges (faces) on each active
// element which are either:
// 1) an internal edge
// 2) an edge on a Neumann boundary for which a boundary condition
// function has been specified.
// The error estimator can be scaled by the number of flux edges (faces)
// which the element actually has to obtain a more uniform measure
// of the error. Use floats instead of ints since in case 2 (above)
// f gets 1/2 of a flux face contribution from each of his
// neighbors
std::vector<float> n_flux_faces (error_per_cell.size());
// 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 (var=0; var<n_vars; var++)
{
// Possibly skip this variable
if (error_norm.weight(var) == 0.0) continue;
// The type of finite element to use for this variable
const FEType& fe_type = dof_map.variable_type (var);
// Finite element objects for the same face from
// different sides
fe_fine = FEBase::build (dim, fe_type);
fe_coarse = FEBase::build (dim, fe_type);
// Build an appropriate Gaussian quadrature rule
QGauss qrule (dim-1, fe_type.default_quadrature_order());
// Tell the finite element for the fine element about the quadrature
// rule. The finite element for the coarse element need not know about it
fe_fine->attach_quadrature_rule (&qrule);
// By convention we will always do the integration
// on the face of element e. We'll need its Jacobian values and
// physical point locations, at least
fe_fine->get_JxW();
fe_fine->get_xyz();
// Our derived classes may want to do some initialization here
this->initialize(system, error_per_cell, estimate_parent_error);
// The global DOF indices for elements e & f
std::vector<unsigned int> dof_indices_fine;
std::vector<unsigned int> dof_indices_coarse;
// 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* e = *elem_it;
const unsigned int e_id = e->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 = e->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.solution),
dof_map, parent, Uparent,
var, false);
// Loop over the neighbors of the parent
for (unsigned int n_p=0; n_p<parent->n_neighbors(); n_p++)
{
if (parent->neighbor(n_p) != NULL) // parent has a neighbor here
{
// Find the active neighbors in this direction
std::vector<const Elem*> active_neighbors;
parent->neighbor(n_p)->
active_family_tree_by_neighbor(active_neighbors,
parent);
// Compute the flux to each active neighbor
for (unsigned int a=0;
a != active_neighbors.size(); ++a)
{
const Elem *f = active_neighbors[a];
// FIXME - what about when f->level <
// parent->level()??
if (f->level() >= parent->level())
{
fine_elem = f;
coarse_elem = parent;
Ucoarse = Uparent;
dof_map.dof_indices (fine_elem, dof_indices_fine, var);
const unsigned int n_dofs_fine = dof_indices_fine.size();
Ufine.resize(n_dofs_fine);
for (unsigned int i=0; i<n_dofs_fine; i++)
Ufine(i) = system.current_solution(dof_indices_fine[i]);
this->reinit_sides();
this->internal_side_integration();
error_per_cell[fine_elem->id()] += fine_error;
error_per_cell[coarse_elem->id()] += coarse_error;
// Keep track of the number of internal flux
// sides found on each element
n_flux_faces[fine_elem->id()]++;
n_flux_faces[coarse_elem->id()] += this->coarse_n_flux_faces_increment();
}
}
}
else if (integrate_boundary_sides)
{
fine_elem = parent;
Ufine = Uparent;
// Reinitialize shape functions on the fine element side
fe_fine->reinit (fine_elem, fine_side);
if (this->boundary_side_integration())
{
error_per_cell[fine_elem->id()] += fine_error;
n_flux_faces[fine_elem->id()]++;
}
}
}
}
#endif // #ifdef LIBMESH_ENABLE_AMR
// If we do any more flux integration, e will be the fine element
fine_elem = e;
// Loop over the neighbors of element e
for (unsigned int n_e=0; n_e<e->n_neighbors(); n_e++)
{
fine_side = n_e;
if (e->neighbor(n_e) != NULL) // e is not on the boundary
{
const Elem* f = e->neighbor(n_e);
const unsigned int f_id = f->id();
// Compute flux jumps if we are in case 1 or case 2.
if ((f->active() && (f->level() == e->level()) && (e_id < f_id))
|| (f->level() < e->level()))
{
// f is now the coarse element
coarse_elem = f;
// Get the DOF indices for the two elements
dof_map.dof_indices (fine_elem, dof_indices_fine, var);
dof_map.dof_indices (coarse_elem, dof_indices_coarse, var);
// The number of DOFS on each element
const unsigned int n_dofs_fine = dof_indices_fine.size();
const unsigned int n_dofs_coarse = dof_indices_coarse.size();
Ufine.resize(n_dofs_fine);
Ucoarse.resize(n_dofs_coarse);
// The local solutions on each element
for (unsigned int i=0; i<n_dofs_fine; i++)
Ufine(i) = system.current_solution(dof_indices_fine[i]);
for (unsigned int i=0; i<n_dofs_coarse; i++)
Ucoarse(i) = system.current_solution(dof_indices_coarse[i]);
this->reinit_sides();
this->internal_side_integration();
error_per_cell[fine_elem->id()] += fine_error;
error_per_cell[coarse_elem->id()] += coarse_error;
// Keep track of the number of internal flux
// sides found on each element
n_flux_faces[fine_elem->id()]++;
n_flux_faces[coarse_elem->id()] += this->coarse_n_flux_faces_increment();
} // end if (case1 || case2)
} // if (e->neigbor(n_e) != NULL)
// Otherwise, e is on the boundary. If it happens to
// be on a Dirichlet boundary, we need not do anything.
// On the other hand, if e is on a Neumann (flux) boundary
// with grad(u).n = g, we need to compute the additional residual
// (h * int |g - grad(u_h).n|^2 dS)^(1/2).
// We can only do this with some knowledge of the boundary
// conditions, i.e. the user must have attached an appropriate
// BC function.
else
{
if (integrate_boundary_sides)
{
// Reinitialize shape functions on the fine element side
fe_fine->reinit (fine_elem, fine_side);
// Get the DOF indices
dof_map.dof_indices (fine_elem, dof_indices_fine, var);
// The number of DOFS on each element
const unsigned int n_dofs_fine = dof_indices_fine.size();
Ufine.resize(n_dofs_fine);
for (unsigned int i=0; i<n_dofs_fine; i++)
Ufine(i) = system.current_solution(dof_indices_fine[i]);
if (this->boundary_side_integration())
{
error_per_cell[fine_elem->id()] += fine_error;
n_flux_faces[fine_elem->id()]++;
}
} // end if _bc_function != NULL
} // end if (e->neighbor(n_e) == NULL)
} // end loop over neighbors
} // 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.
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]);
if (this->scale_by_n_flux_faces)
{
// Sum the vector of flux face counts
this->reduce_error(n_flux_faces);
// Sanity check: Make sure the number of flux faces is
// always an integer value
#ifdef DEBUG
for (unsigned int i=0; i<n_flux_faces.size(); ++i)
libmesh_assert (n_flux_faces[i] == static_cast<float>(static_cast<unsigned int>(n_flux_faces[i])) );
#endif
// Scale the error by the number of flux faces for each element
for (unsigned int i=0; i<n_flux_faces.size(); ++i)
{
if (n_flux_faces[i] == 0.0) // inactive or non-local element
continue;
//std::cout << 'Element ' << i << ' has ' << n_flux_faces[i] << ' flux faces.' << std::endl;
error_per_cell[i] /= static_cast<Real>(n_flux_faces[i]);
}
}
// 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();
}
STOP_LOG('estimate_error()', 'JumpErrorEstimator');
}
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;
}
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;
}
Reimplemented in DiscontinuityMeasure, LaplacianErrorEstimator, and KellyErrorEstimator.
Definition at line 44 of file jump_error_estimator.C.
Referenced by estimate_error().
{
}
Implemented in DiscontinuityMeasure, LaplacianErrorEstimator, and KellyErrorEstimator.
Referenced by estimate_error().
Definition at line 32 of file error_estimator.C.
Referenced by UniformRefinementEstimator::_estimate_error(), PatchRecoveryErrorEstimator::estimate_error(), and 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 412 of file jump_error_estimator.C.
References coarse_elem, Elem::dim(), fe_coarse, fe_fine, fine_elem, fine_side, and InfFE< Dim, T_radial, T_map >::inverse_map().
Referenced by estimate_error().
{
// The master quadrature point locations on the coarse element
std::vector<Point> qp_coarse;
// Reinitialize shape functions on the fine element side
fe_fine->reinit (fine_elem, fine_side);
// Get the physical locations of the fine element quadrature points
std::vector<Point> qface_point = fe_fine->get_xyz();
// Find their locations on the coarse element
FEInterface::inverse_map (coarse_elem->dim(), fe_coarse->get_fe_type(),
coarse_elem, qface_point, qp_coarse);
// Calculate the coarse element shape functions at those locations
fe_coarse->reinit (coarse_elem, &qp_coarse);
}
Definition at line 130 of file jump_error_estimator.h.
Referenced by coarse_n_flux_faces_increment(), estimate_error(), KellyErrorEstimator::internal_side_integration(), LaplacianErrorEstimator::internal_side_integration(), DiscontinuityMeasure::internal_side_integration(), and reinit_sides().
Definition at line 135 of file jump_error_estimator.h.
Referenced by estimate_error(), KellyErrorEstimator::internal_side_integration(), LaplacianErrorEstimator::internal_side_integration(), and DiscontinuityMeasure::internal_side_integration().
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(), 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::UniformRefinementEstimator().
Definition at line 155 of file jump_error_estimator.h.
Referenced by estimate_error(), KellyErrorEstimator::initialize(), LaplacianErrorEstimator::initialize(), DiscontinuityMeasure::initialize(), KellyErrorEstimator::internal_side_integration(), LaplacianErrorEstimator::internal_side_integration(), DiscontinuityMeasure::internal_side_integration(), and reinit_sides().
Definition at line 155 of file jump_error_estimator.h.
Referenced by KellyErrorEstimator::boundary_side_integration(), DiscontinuityMeasure::boundary_side_integration(), estimate_error(), KellyErrorEstimator::initialize(), LaplacianErrorEstimator::initialize(), DiscontinuityMeasure::initialize(), KellyErrorEstimator::internal_side_integration(), LaplacianErrorEstimator::internal_side_integration(), DiscontinuityMeasure::internal_side_integration(), and reinit_sides().
Definition at line 130 of file jump_error_estimator.h.
Referenced by KellyErrorEstimator::boundary_side_integration(), DiscontinuityMeasure::boundary_side_integration(), coarse_n_flux_faces_increment(), estimate_error(), KellyErrorEstimator::internal_side_integration(), LaplacianErrorEstimator::internal_side_integration(), DiscontinuityMeasure::internal_side_integration(), and reinit_sides().
Definition at line 135 of file jump_error_estimator.h.
Referenced by KellyErrorEstimator::boundary_side_integration(), DiscontinuityMeasure::boundary_side_integration(), estimate_error(), KellyErrorEstimator::internal_side_integration(), LaplacianErrorEstimator::internal_side_integration(), and DiscontinuityMeasure::internal_side_integration().
Definition at line 140 of file jump_error_estimator.h.
Referenced by estimate_error(), and reinit_sides().
Definition at line 125 of file jump_error_estimator.h.
Referenced by KellyErrorEstimator::attach_flux_bc_function(), and estimate_error().
Definition at line 86 of file jump_error_estimator.h.
Referenced by estimate_error().
Definition at line 150 of file jump_error_estimator.h.
Referenced by estimate_error(), KellyErrorEstimator::internal_side_integration(), LaplacianErrorEstimator::internal_side_integration(), and DiscontinuityMeasure::internal_side_integration().
Definition at line 150 of file jump_error_estimator.h.
Referenced by KellyErrorEstimator::boundary_side_integration(), DiscontinuityMeasure::boundary_side_integration(), estimate_error(), KellyErrorEstimator::internal_side_integration(), LaplacianErrorEstimator::internal_side_integration(), and DiscontinuityMeasure::internal_side_integration().
Definition at line 145 of file jump_error_estimator.h.
Referenced by KellyErrorEstimator::boundary_side_integration(), DiscontinuityMeasure::boundary_side_integration(), estimate_error(), KellyErrorEstimator::internal_side_integration(), LaplacianErrorEstimator::internal_side_integration(), and DiscontinuityMeasure::internal_side_integration().
Generated automatically by Doxygen for libMesh from the source code.