EstimateError (const System &sys, const PatchRecoveryErrorEstimator &ee, ErrorVector &epc)
void operator() (const ConstElemRange &range) const
const System & system
const PatchRecoveryErrorEstimator & error_estimator
ErrorVector & error_per_cell
Class to compute the error contribution for a range of elements. May be executed in parallel on separate threads.
Definition at line 105 of file patch_recovery_error_estimator.h.
Definition at line 108 of file patch_recovery_error_estimator.h.
:
system(sys),
error_estimator(ee),
error_per_cell(epc)
{}
Definition at line 170 of file patch_recovery_error_estimator.C.
References TypeTensor< T >::add_scaled(), TypeVector< T >::add_scaled(), StoredRange< iterator_type, object_type >::begin(), FEBase::build(), Patch::build_around_element(), System::current_solution(), FEType::default_quadrature_rule(), DofMap::dof_indices(), StoredRange< iterator_type, object_type >::end(), error_estimator, ErrorEstimator::error_norm, error_per_cell, System::get_dof_map(), System::get_mesh(), libMeshEnums::H1_SEMINORM, libMeshEnums::H2_SEMINORM, DofObject::id(), libmesh_norm(), DenseMatrix< T >::lu_solve(), DenseMatrixBase< T >::m(), std::max(), MeshBase::mesh_dimension(), DenseMatrixBase< T >::n(), System::n_vars(), FEType::order, Elem::p_level(), PatchRecoveryErrorEstimator::patch_growth_strategy, DenseVector< T >::resize(), PatchRecoveryErrorEstimator::specpoly(), system, PatchRecoveryErrorEstimator::target_patch_size, SystemNorm::type(), DofMap::variable_type(), libMeshEnums::W1_INF_SEMINORM, libMeshEnums::W2_INF_SEMINORM, SystemNorm::weight(), and SystemNorm::weight_sq().
{
// 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();
//------------------------------------------------------------
// Iterate over all the elements in the range.
for (ConstElemRange::const_iterator elem_it=range.begin(); elem_it!=range.end(); ++elem_it)
{
// elem is necessarily an active element on the local processor
const Elem* elem = *elem_it;
// We'll need an index into the error vector
const int e_id=elem->id();
// Build a patch containing the current element
// and its neighbors on the local processor
Patch patch;
// Use user specified patch size and growth strategy
patch.build_around_element (elem, error_estimator.target_patch_size,
error_estimator.patch_growth_strategy);
//------------------------------------------------------------
// Process each variable in the system using the current patch
for (unsigned int var=0; var<n_vars; var++)
{
#ifndef DEBUG
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
libmesh_assert (error_estimator.error_norm.type(var) == H1_SEMINORM ||
error_estimator.error_norm.type(var) == H2_SEMINORM ||
error_estimator.error_norm.type(var) == W1_INF_SEMINORM ||
error_estimator.error_norm.type(var) == W2_INF_SEMINORM);
#else
libmesh_assert (error_estimator.error_norm.type(var) == H1_SEMINORM ||
error_estimator.error_norm.type(var) == W1_INF_SEMINORM);
#endif
if (var > 0)
// We can't mix L_inf and L_2 norms
libmesh_assert (((error_estimator.error_norm.type(var) == H1_SEMINORM ||
error_estimator.error_norm.type(var) == H2_SEMINORM) &&
(error_estimator.error_norm.type(var-1) == H1_SEMINORM ||
error_estimator.error_norm.type(var-1) == H2_SEMINORM)) ||
((error_estimator.error_norm.type(var) == W1_INF_SEMINORM ||
error_estimator.error_norm.type(var) == W2_INF_SEMINORM) &&
(error_estimator.error_norm.type(var-1) == W1_INF_SEMINORM ||
error_estimator.error_norm.type(var-1) == W2_INF_SEMINORM)));
#endif
// Possibly skip this variable
if (error_estimator.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);
const Order element_order = static_cast<Order>
(fe_type.order + elem->p_level());
// Finite element object for use in this patch
AutoPtr<FEBase> fe (FEBase::build (dim, fe_type));
// Build an appropriate Gaussian quadrature rule
AutoPtr<QBase> qrule (fe_type.default_quadrature_rule(dim));
// Tell the finite element about the quadrature rule.
fe->attach_quadrature_rule (qrule.get());
// Get Jacobian values, etc..
const std::vector<Real>& JxW = fe->get_JxW();
const std::vector<Point>& q_point = fe->get_xyz();
#ifndef NDEBUG
// We only use phi below to assert that the dof_indices
// vector has the correct size. So we avoid declaring it
// here unless asserts are active.
const std::vector<std::vector<Real> >& phi = fe->get_phi();
#endif
const std::vector<std::vector<RealGradient> > *dphi = NULL;
if (error_estimator.error_norm.type(var) == H1_SEMINORM ||
error_estimator.error_norm.type(var) == W1_INF_SEMINORM)
dphi = &(fe->get_dphi());
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
const std::vector<std::vector<RealTensor> > *d2phi = NULL;
if (error_estimator.error_norm.type(var) == H2_SEMINORM ||
error_estimator.error_norm.type(var) == W1_INF_SEMINORM)
d2phi = &(fe->get_d2phi());
#endif
// global DOF indices
std::vector<unsigned int> dof_indices;
// Compute the approprite size for the patch projection matrices
// and vectors;
unsigned int matsize = element_order + 1;
if (dim > 1)
{
matsize *= (element_order + 2);
matsize /= 2;
}
if (dim > 2)
{
matsize *= (element_order + 3);
matsize /= 3;
}
DenseMatrix<Number> Kp(matsize,matsize);
DenseVector<Number>
Fx(matsize), Pu_x_h(matsize), // Also xx
Fy(matsize), Pu_y_h(matsize), // Also yy
Fz(matsize), Pu_z_h(matsize); // Also zz
DenseVector<Number> Fxy, Pu_xy_h, Fxz, Pu_xz_h, Fyz, Pu_yz_h;
if (error_estimator.error_norm.type(var) == H2_SEMINORM ||
error_estimator.error_norm.type(var) == W2_INF_SEMINORM)
{
Fxy.resize(matsize); Pu_xy_h.resize(matsize);
Fxz.resize(matsize); Pu_xz_h.resize(matsize);
Fyz.resize(matsize); Pu_yz_h.resize(matsize);
}
//------------------------------------------------------
// Loop over each element in the patch and compute their
// contribution to the patch gradient projection.
Patch::const_iterator patch_it = patch.begin();
const Patch::const_iterator patch_end = patch.end();
for (; patch_it != patch_end; ++patch_it)
{
// The pth element in the patch
const Elem* e_p = *patch_it;
// Reinitialize the finite element data for this element
fe->reinit (e_p);
// Get the global DOF indices for the current variable
// in the current element
dof_map.dof_indices (e_p, dof_indices, var);
libmesh_assert (dof_indices.size() == phi.size());
const unsigned int n_dofs = dof_indices.size();
const unsigned int n_qp = qrule->n_points();
// Compute the projection components from this cell.
// int_{Omega_e}
si_i
si_j = int_{Omega_e} du_h/dx_k
si_i
for (unsigned int qp=0; qp<n_qp; qp++)
{
// Construct the shape function values for the patch projection
std::vector<Real> psi(specpoly(dim, element_order, q_point[qp], matsize));
// Patch matrix contribution
for (unsigned int i=0; i<Kp.m(); i++)
for (unsigned int j=0; j<Kp.n(); j++)
Kp(i,j) += JxW[qp]*psi[i]*psi[j];
if (error_estimator.error_norm.type(var) == H1_SEMINORM ||
error_estimator.error_norm.type(var) == W1_INF_SEMINORM)
{
// Compute the gradient on the current patch element
// at the quadrature point
Gradient grad_u_h;
for (unsigned int i=0; i<n_dofs; i++)
grad_u_h.add_scaled ((*dphi)[i][qp],
system.current_solution(dof_indices[i]));
// Patch RHS contributions
for (unsigned int i=0; i<psi.size(); i++)
{
Fx(i) += JxW[qp]*grad_u_h(0)*psi[i];
Fy(i) += JxW[qp]*grad_u_h(1)*psi[i];
Fz(i) += JxW[qp]*grad_u_h(2)*psi[i];
}
}
else if (error_estimator.error_norm.type(var) == H2_SEMINORM ||
error_estimator.error_norm.type(var) == W2_INF_SEMINORM)
{
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
// Compute the hessian on the current patch element
// at the quadrature point
Tensor hess_u_h;
for (unsigned int i=0; i<n_dofs; i++)
hess_u_h.add_scaled ((*d2phi)[i][qp],
system.current_solution(dof_indices[i]));
// Patch RHS contributions
for (unsigned int i=0; i<psi.size(); i++)
{
Fx(i) += JxW[qp]*hess_u_h(0,0)*psi[i];
Fy(i) += JxW[qp]*hess_u_h(1,1)*psi[i];
Fz(i) += JxW[qp]*hess_u_h(2,2)*psi[i];
Fxy(i) += JxW[qp]*hess_u_h(0,1)*psi[i];
Fxz(i) += JxW[qp]*hess_u_h(0,2)*psi[i];
Fyz(i) += JxW[qp]*hess_u_h(1,2)*psi[i];
}
#else
std::cerr << 'ERROR: --enable-second-derivatives is required
<< ' for _sobolev_order == 2!;
libmesh_error();
#endif
}
else
libmesh_error();
} // end quadrature loop
} // end patch loop
//--------------------------------------------------
// Now we have fully assembled the projection system
// for this patch. Project the gradient components.
// MAY NEED TO USE PARTIAL PIVOTING!
Kp.lu_solve (Fx, Pu_x_h);
Kp.lu_solve (Fy, Pu_y_h);
Kp.lu_solve (Fz, Pu_z_h);
if (error_estimator.error_norm.type(var) == H2_SEMINORM ||
error_estimator.error_norm.type(var) == W2_INF_SEMINORM)
{
Kp.lu_solve(Fxy, Pu_xy_h);
Kp.lu_solve(Fxz, Pu_xz_h);
Kp.lu_solve(Fyz, Pu_yz_h);
}
//--------------------------------------------------
// Finally, estimate the error in the current variable
// for the current element by computing ||P grad_u_h - grad_u_h||
// or ||P hess_u_h - hess_u_h|| according to the requested
// seminorm
dof_map.dof_indices (elem, dof_indices, var);
const unsigned int n_dofs = dof_indices.size();
Real element_error = 0;
// we approximate the max norm by sampling over a set of points
// in future we may add specialized routines for specific cases
// or use some optimization package
const Order qorder = element_order;
// build a 'fake' quadrature rule for the element
//
// For linear elements, grad is a constant, so we need
// to compute grad u_h once on the element. Also as G_H
// u_h - gradu_h is linear on an element, it assumes its
// maximum at a vertex of the element
QGrid samprule (dim, qorder);
if (error_estimator.error_norm.type(var) == W1_INF_SEMINORM ||
error_estimator.error_norm.type(var) == W2_INF_SEMINORM)
fe->attach_quadrature_rule (&samprule);
// reinitialize element for integration or sampling
fe->reinit(elem);
const unsigned int n_sp = JxW.size();
for (unsigned int sp=0; sp< n_sp; sp++)
{
std::vector<Number> temperr(6,0.0); // x,y,z or xx,yy,zz,xy,xz,yz
if (error_estimator.error_norm.type(var) == H1_SEMINORM ||
error_estimator.error_norm.type(var) == W1_INF_SEMINORM)
{
// Compute the gradient at the current sample point
Gradient grad_u_h;
for (unsigned int i=0; i<n_dofs; i++)
grad_u_h.add_scaled ((*dphi)[i][sp],
system.current_solution(dof_indices[i]));
// Compute the phi values at the current sample point
std::vector<Real> psi(specpoly(dim, element_order, q_point[sp], matsize));
for (unsigned int i=0; i<matsize; i++)
{
temperr[0] += psi[i]*Pu_x_h(i);
temperr[1] += psi[i]*Pu_y_h(i);
temperr[2] += psi[i]*Pu_z_h(i);
}
temperr[0] -= grad_u_h(0);
temperr[1] -= grad_u_h(1);
temperr[2] -= grad_u_h(2);
}
else if (error_estimator.error_norm.type(var) == H2_SEMINORM ||
error_estimator.error_norm.type(var) == W2_INF_SEMINORM)
{
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
// Compute the Hessian at the current sample point
Tensor hess_u_h;
for (unsigned int i=0; i<n_dofs; i++)
hess_u_h.add_scaled ((*d2phi)[i][sp],
system.current_solution(dof_indices[i]));
// Compute the phi values at the current sample point
std::vector<Real> psi(specpoly(dim, element_order, q_point[sp], matsize));
for (unsigned int i=0; i<matsize; i++)
{
temperr[0] += psi[i]*Pu_x_h(i);
temperr[1] += psi[i]*Pu_y_h(i);
temperr[2] += psi[i]*Pu_z_h(i);
temperr[3] += psi[i]*Pu_xy_h(i);
temperr[4] += psi[i]*Pu_xz_h(i);
temperr[5] += psi[i]*Pu_yz_h(i);
}
temperr[0] -= hess_u_h(0,0);
temperr[1] -= hess_u_h(1,1);
temperr[2] -= hess_u_h(2,2);
temperr[3] -= hess_u_h(0,1);
temperr[4] -= hess_u_h(0,2);
temperr[5] -= hess_u_h(1,2);
#else
std::cerr << 'ERROR: --enable-second-derivatives is required
<< ' for _sobolev_order == 2!;
libmesh_error();
#endif
}
if (error_estimator.error_norm.type(var) == W1_INF_SEMINORM)
for (unsigned int i=0; i != 3; ++i)
element_error = std::max(element_error, std::abs(temperr[i]));
else if (error_estimator.error_norm.type(var) == W2_INF_SEMINORM)
for (unsigned int i=0; i != 6; ++i)
element_error = std::max(element_error, std::abs(temperr[i]));
else if (error_estimator.error_norm.type(var) == H1_SEMINORM)
for (unsigned int i=0; i != 3; ++i)
element_error += JxW[sp]*libmesh_norm(temperr[i]);
else if (error_estimator.error_norm.type(var) == H2_SEMINORM)
{
for (unsigned int i=0; i != 3; ++i)
element_error += JxW[sp]*libmesh_norm(temperr[i]);
// Off diagonal terms enter into the Hessian norm twice
for (unsigned int i=3; i != 6; ++i)
element_error += JxW[sp]*2*libmesh_norm(temperr[i]);
}
} // end sample_point_loop
// The patch error estimator works element-by-element --
// there is no need to get a mutex on error_per_cell!
if (error_estimator.error_norm.type(var) == W1_INF_SEMINORM ||
error_estimator.error_norm.type(var) == W2_INF_SEMINORM)
error_per_cell[e_id] += error_estimator.error_norm.weight(var) * element_error;
else if (error_estimator.error_norm.type(var) == H1_SEMINORM ||
error_estimator.error_norm.type(var) == H2_SEMINORM)
error_per_cell[e_id] += error_estimator.error_norm.weight_sq(var) * element_error;
else
libmesh_error();
} // end variable loop
if (error_estimator.error_norm.type(0) == H1_SEMINORM ||
error_estimator.error_norm.type(0) == H2_SEMINORM)
error_per_cell[e_id] = std::sqrt(error_per_cell[e_id]);
} // end element loop
}
Definition at line 121 of file patch_recovery_error_estimator.h.
Definition at line 122 of file patch_recovery_error_estimator.h.
Definition at line 120 of file patch_recovery_error_estimator.h.
Referenced by operator()().
Generated automatically by Doxygen for libMesh from the source code.