#include <patch_recovery_error_estimator.h>
typedef std::map< std::pair< const System *, unsigned int >, ErrorVector * > ErrorMap
PatchRecoveryErrorEstimator ()
~PatchRecoveryErrorEstimator ()
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 int target_patch_size
Patch::PMF patch_growth_strategy
SystemNorm error_norm
void reduce_error (std::vector< float > &error_per_cell) const
static std::vector< Real > specpoly (const unsigned int dim, const Order order, const Point p, const unsigned int matsize)
This class implements the Patch Recovery error indicator.
Author:
Definition at line 45 of file patch_recovery_error_estimator.h.
Definition at line 105 of file error_estimator.h.
Definition at line 54 of file patch_recovery_error_estimator.h.
References ErrorEstimator::error_norm, and libMeshEnums::H1_SEMINORM.
:
target_patch_size(20),
patch_growth_strategy(&Patch::add_local_face_neighbors)
{ error_norm = H1_SEMINORM; }
Definition at line 62 of file patch_recovery_error_estimator.h.
{}
Implements ErrorEstimator.
Definition at line 111 of file patch_recovery_error_estimator.C.
References MeshBase::active_local_elements_begin(), MeshBase::active_local_elements_end(), AutoPtr< Tp >::get(), System::get_mesh(), MeshBase::max_elem_id(), Threads::parallel_for(), ErrorEstimator::reduce_error(), System::solution, and NumericVector< T >::swap().
{
START_LOG('estimate_error()', 'PatchRecoveryErrorEstimator');
// The current mesh
const MeshBase& mesh = system.get_mesh();
// 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();
}
//------------------------------------------------------------
// Iterate over all the active elements in the mesh
// that live on this processor.
Threads::parallel_for (ConstElemRange(mesh.active_local_elements_begin(),
mesh.active_local_elements_end(),
200),
EstimateError(system,
*this,
error_per_cell)
);
// Each processor has now computed the error contribuions
// for its local elements, and error_per_cell contains 0 for all the
// non-local elements. Summing the vector will provide the L_oo value
// for each element, local or remote
this->reduce_error(error_per_cell);
// 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()', 'PatchRecoveryErrorEstimator');
}
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;
}
Definition at line 32 of file error_estimator.C.
Referenced by UniformRefinementEstimator::_estimate_error(), 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 47 of file patch_recovery_error_estimator.C.
References Utility::pow().
Referenced by PatchRecoveryErrorEstimator::EstimateError::operator()().
{
std::vector<Real> psi;
psi.reserve(matsize);
Real x = p(0), y=0., z=0.;
if (dim > 1)
y = p(1);
if (dim > 2)
z = p(2);
// builds psi vector of form 1 x y z x^2 xy xz y^2 yz z^2 etc..
// I haven't added 1D support here
for (unsigned int poly_deg=0; poly_deg <= static_cast<unsigned int>(order) ; poly_deg++)
{ // loop over all polynomials of total degreee = poly_deg
switch (dim)
{
// 3D spectral polynomial basis functions
case 3:
{
for (int xexp=poly_deg; xexp >= 0; xexp--) // use an int for xexp since we -- it
for (int yexp=poly_deg-xexp; yexp >= 0; yexp--)
{
int zexp = poly_deg - xexp - yexp;
psi.push_back(std::pow(x,static_cast<Real>(xexp))*
std::pow(y,static_cast<Real>(yexp))*
std::pow(z,static_cast<Real>(zexp)));
}
break;
}
// 2D spectral polynomial basis functions
case 2:
{
for (int xexp=poly_deg; xexp >= 0; xexp--) // use an int for xexp since we -- it
{
int yexp = poly_deg - xexp;
psi.push_back(std::pow(x,static_cast<Real>(xexp))*
std::pow(y,static_cast<Real>(yexp)));
}
break;
}
// 1D spectral polynomial basis functions
case 1:
{
int xexp = poly_deg;
psi.push_back(std::pow(x,static_cast<Real>(xexp)));
break;
}
default:
libmesh_error();
}
}
return psi;
}
Definition at line 125 of file patch_recovery_error_estimator.h.
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::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(), and UniformRefinementEstimator::UniformRefinementEstimator().
Definition at line 88 of file patch_recovery_error_estimator.h.
Referenced by PatchRecoveryErrorEstimator::EstimateError::operator()().
Definition at line 80 of file patch_recovery_error_estimator.h.
Referenced by PatchRecoveryErrorEstimator::EstimateError::operator()().
Generated automatically by Doxygen for libMesh from the source code.