Poster of Linux kernelThe best gift for a Linux geek
PatchRecoveryErrorEstimator::EstimateError

PatchRecoveryErrorEstimator::EstimateError

Section: C Library Functions (3) Updated: Thu Apr 7 2011
Local index Up
 

NAME

PatchRecoveryErrorEstimator::EstimateError -  

SYNOPSIS


 

Public Member Functions


EstimateError (const System &sys, const PatchRecoveryErrorEstimator &ee, ErrorVector &epc)

void operator() (const ConstElemRange &range) const
 

Private Attributes


const System & system

const PatchRecoveryErrorEstimator & error_estimator

ErrorVector & error_per_cell
 

Detailed Description

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.  

Constructor & Destructor Documentation

 

PatchRecoveryErrorEstimator::EstimateError::EstimateError (const System &sys, const PatchRecoveryErrorEstimator &ee, ErrorVector &epc) [inline]

Definition at line 108 of file patch_recovery_error_estimator.h.

                                     :
      system(sys),
      error_estimator(ee),
      error_per_cell(epc)
    {}
 

Member Function Documentation

 

void PatchRecoveryErrorEstimator::EstimateError::operator() (const ConstElemRange &range) const

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 }
 

Member Data Documentation

 

const PatchRecoveryErrorEstimator& PatchRecoveryErrorEstimator::EstimateError::error_estimator [private]

Definition at line 121 of file patch_recovery_error_estimator.h.

Referenced by operator()().  

ErrorVector& PatchRecoveryErrorEstimator::EstimateError::error_per_cell [private]

Definition at line 122 of file patch_recovery_error_estimator.h.

Referenced by operator()().  

const System& PatchRecoveryErrorEstimator::EstimateError::system [private]

Definition at line 120 of file patch_recovery_error_estimator.h.

Referenced by operator()().

 

Author

Generated automatically by Doxygen for libMesh from the source code.


 

Index

NAME
SYNOPSIS
Public Member Functions
Private Attributes
Detailed Description
Constructor & Destructor Documentation
PatchRecoveryErrorEstimator::EstimateError::EstimateError (const System &sys, const PatchRecoveryErrorEstimator &ee, ErrorVector &epc) [inline]
Member Function Documentation
void PatchRecoveryErrorEstimator::EstimateError::operator() (const ConstElemRange &range) const
Member Data Documentation
const PatchRecoveryErrorEstimator& PatchRecoveryErrorEstimator::EstimateError::error_estimator [private]
ErrorVector& PatchRecoveryErrorEstimator::EstimateError::error_per_cell [private]
const System& PatchRecoveryErrorEstimator::EstimateError::system [private]
Author

This document was created by man2html, using the manual pages.
Time: 21:52:04 GMT, April 16, 2011