libMesh::ExactSolution Class Reference

#include <exact_solution.h>

List of all members.

Public Member Functions

 ExactSolution (const EquationSystems &es)
 ~ExactSolution ()
void attach_reference_solution (const EquationSystems *es_fine)
void attach_exact_values (std::vector< FunctionBase< Number > * > f)
void attach_exact_value (unsigned int sys_num, FunctionBase< Number > *f)
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_derivs (std::vector< FunctionBase< Gradient > * > g)
void attach_exact_deriv (unsigned int sys_num, FunctionBase< Gradient > *g)
void attach_exact_deriv (Gradient gptr(const Point &p, const Parameters &parameters, const std::string &sys_name, const std::string &unknown_name))
void attach_exact_hessians (std::vector< FunctionBase< Tensor > * > h)
void attach_exact_hessian (unsigned int sys_num, FunctionBase< Tensor > *h)
void attach_exact_hessian (Tensor hptr(const Point &p, const Parameters &parameters, const std::string &sys_name, const std::string &unknown_name))
void extra_quadrature_order (const int extraorder)
void compute_error (const std::string &sys_name, const std::string &unknown_name)
Real l2_error (const std::string &sys_name, const std::string &unknown_name)
Real l1_error (const std::string &sys_name, const std::string &unknown_name)
Real l_inf_error (const std::string &sys_name, const std::string &unknown_name)
Real h1_error (const std::string &sys_name, const std::string &unknown_name)
Real hcurl_error (const std::string &sys_name, const std::string &unknown_name)
Real hdiv_error (const std::string &sys_name, const std::string &unknown_name)
Real h2_error (const std::string &sys_name, const std::string &unknown_name)
Real error_norm (const std::string &sys_name, const std::string &unknown_name, const FEMNormType &norm)

Private Types

typedef std::map< std::string,
std::vector< Real > > 
SystemErrorMap

Private Member Functions

template<typename OutputShape >
void _compute_error (const std::string &sys_name, const std::string &unknown_name, std::vector< Real > &error_vals)
std::vector< Real > & _check_inputs (const std::string &sys_name, const std::string &unknown_name)

Private Attributes

std::vector< FunctionBase
< Number > * > 
_exact_values
std::vector< FunctionBase
< Gradient > * > 
_exact_derivs
std::vector< FunctionBase
< Tensor > * > 
_exact_hessians
std::map< std::string,
SystemErrorMap
_errors
const EquationSystems_equation_systems
const EquationSystems_equation_systems_fine
int _extra_order

Detailed Description

This class handles the computation of the L2 and/or H1 error for the Systems in the EquationSystems object which is passed to it. Note that for it to be useful, the user must attach at least one, and possibly two functions which can compute the exact solution and its derivative for any component of any system. These are the exact_value and exact_deriv functions below.

Author:
Benjamin S. Kirk w/ modifications for libmesh by John W. Peterson

Definition at line 64 of file exact_solution.h.


Member Typedef Documentation

typedef std::map<std::string, std::vector<Real> > libMesh::ExactSolution::SystemErrorMap [private]

Data structure which stores the errors: ||e|| = ||u - u_h|| ||grad(e)|| = ||grad(u) - grad(u_h)|| for each unknown in a single system. The name of the unknown is the key for the map.

Definition at line 293 of file exact_solution.h.


Constructor & Destructor Documentation

libMesh::ExactSolution::ExactSolution ( const EquationSystems es  )  [explicit]

Constructor. The ExactSolution object must be initialized with an EquationSystems object.

Definition at line 41 of file exact_solution.C.

References _equation_systems, _errors, libMesh::EquationSystems::get_system(), libMesh::EquationSystems::n_systems(), libMesh::System::n_vars(), libMesh::System::name(), and libMesh::System::variable_name().

00041                                                       :
00042   _equation_systems(es),
00043   _equation_systems_fine(NULL),
00044   _extra_order(0)
00045 {
00046   // Initialize the _errors data structure which holds all
00047   // the eventual values of the error.
00048   for (unsigned int sys=0; sys<_equation_systems.n_systems(); ++sys)
00049     {
00050       // Reference to the system
00051       const System& system = _equation_systems.get_system(sys);
00052 
00053       // The name of the system
00054       const std::string& sys_name = system.name();
00055 
00056       // The SystemErrorMap to be inserted
00057       ExactSolution::SystemErrorMap sem;
00058 
00059       for (unsigned int var=0; var<system.n_vars(); ++var)
00060         {
00061           // The name of this variable
00062           const std::string& var_name = system.variable_name(var);
00063           sem[var_name] = std::vector<Real>(5, 0.);
00064         }
00065 
00066       _errors[sys_name] = sem;
00067     }
00068 }

libMesh::ExactSolution::~ExactSolution (  ) 

Destructor.

Definition at line 71 of file exact_solution.C.

References _exact_derivs, _exact_hessians, and _exact_values.

00072 {
00073   // delete will clean up any cloned functors and no-op on any NULL
00074   // pointers
00075 
00076   for (unsigned int i=0; i != _exact_values.size(); ++i)
00077     delete (_exact_values[i]);
00078 
00079   for (unsigned int i=0; i != _exact_derivs.size(); ++i)
00080     delete (_exact_derivs[i]);
00081 
00082   for (unsigned int i=0; i != _exact_hessians.size(); ++i)
00083     delete (_exact_hessians[i]);
00084 }


Member Function Documentation

std::vector< Real > & libMesh::ExactSolution::_check_inputs ( const std::string &  sys_name,
const std::string &  unknown_name 
) [private]

This function is responsible for checking the validity of the sys_name and unknown_name inputs, and returning a reference to the proper vector for storing the values.

Definition at line 261 of file exact_solution.C.

References _errors, and libMesh::err.

Referenced by compute_error(), error_norm(), h1_error(), h2_error(), l1_error(), l2_error(), and l_inf_error().

00263 {
00264   // If no exact solution function or fine grid solution has been
00265   // attached, we now just compute the solution norm (i.e. the
00266   // difference from an "exact solution" of zero
00267 
00268   // Make sure the requested sys_name exists.
00269   std::map<std::string, SystemErrorMap>::iterator sys_iter =
00270     _errors.find(sys_name);
00271 
00272   if (sys_iter == _errors.end())
00273     {
00274       libMesh::err << "Sorry, couldn't find the requested system '"
00275                     << sys_name << "'."
00276                     << std::endl;
00277       libmesh_error();
00278     }
00279 
00280   // Make sure the requested unknown_name exists.
00281   SystemErrorMap::iterator var_iter = (*sys_iter).second.find(unknown_name);
00282 
00283   if (var_iter == (*sys_iter).second.end())
00284     {
00285       libMesh::err << "Sorry, couldn't find the requested variable '"
00286                     << unknown_name << "'."
00287                     << std::endl;
00288       libmesh_error();
00289     }
00290 
00291   // Return a reference to the proper error entry
00292   return (*var_iter).second;
00293 }

template<typename OutputShape >
template void libMesh::ExactSolution::_compute_error< RealGradient > ( const std::string &  sys_name,
const std::string &  unknown_name,
std::vector< Real > &  error_vals 
) [inline, private]

This function computes the error (in the solution and its first derivative) for a single unknown in a single system. It is a private function since it is used by the implementation when solving for several unknowns in several systems.

Definition at line 533 of file exact_solution.C.

References _equation_systems, _equation_systems_fine, _exact_derivs, _exact_hessians, _exact_values, _extra_order, libMesh::MeshBase::active_local_elements_begin(), libMesh::MeshBase::active_local_elements_end(), libMesh::CommWorld, libMesh::TensorTools::curl_from_grad(), libMesh::System::current_solution(), libMesh::FEType::default_quadrature_rule(), libMesh::TensorTools::div_from_grad(), libMesh::DofMap::dof_indices(), libMesh::err, libMesh::FEInterface::field_type(), libMesh::AutoPtr< Tp >::get(), libMesh::System::get_dof_map(), libMesh::System::get_mesh(), libMesh::EquationSystems::get_system(), libMesh::Parallel::Communicator::max(), libMesh::MeshBase::mesh_dimension(), libMesh::FEInterface::n_vec_dim(), libMesh::TensorTools::norm_sq(), libMesh::System::number(), libMesh::Real, libMeshEnums::SERIAL, libMesh::TypeVector< T >::size_sq(), libMesh::System::solution, libMesh::MeshBase::spatial_dimension(), libMesh::Parallel::Communicator::sum(), libMeshEnums::TYPE_VECTOR, libMesh::System::update_global_solution(), libMesh::System::variable_number(), libMesh::System::variable_scalar_number(), and libMesh::DofMap::variable_type().

00536 {
00537   // This function must be run on all processors at once
00538   parallel_only();
00539 
00540   // Make sure we aren't "overconfigured"
00541   libmesh_assert (!(_exact_values.size() && _equation_systems_fine));
00542 
00543   // Get a reference to the system whose error is being computed.
00544   // If we have a fine grid, however, we'll integrate on that instead
00545   // for more accuracy.
00546   const System& computed_system = _equation_systems_fine ?
00547     _equation_systems_fine->get_system(sys_name) :
00548     _equation_systems.get_system (sys_name);
00549 
00550   const Real time = _equation_systems.get_system(sys_name).time;
00551 
00552   const unsigned int sys_num = computed_system.number();
00553   const unsigned int var = computed_system.variable_number(unknown_name);
00554   const unsigned int var_component =
00555     computed_system.variable_scalar_number(var, 0);
00556 
00557   // Prepare a global solution and a MeshFunction of the coarse system if we need one
00558   AutoPtr<MeshFunction> coarse_values;
00559   AutoPtr<NumericVector<Number> > comparison_soln = NumericVector<Number>::build();
00560   if (_equation_systems_fine)
00561     {
00562       const System& comparison_system
00563         = _equation_systems.get_system(sys_name);
00564 
00565       std::vector<Number> global_soln;
00566       comparison_system.update_global_solution(global_soln);
00567       comparison_soln->init(comparison_system.solution->size(), true, SERIAL);
00568       (*comparison_soln) = global_soln;
00569 
00570       coarse_values = AutoPtr<MeshFunction>
00571         (new MeshFunction(_equation_systems,
00572                           *comparison_soln,
00573                           comparison_system.get_dof_map(),
00574                           comparison_system.variable_number(unknown_name)));
00575       coarse_values->init();
00576     }
00577 
00578   // Initialize any functors we're going to use
00579   for (unsigned int i=0; i != _exact_values.size(); ++i)
00580     if (_exact_values[i])
00581       _exact_values[i]->init();
00582 
00583   for (unsigned int i=0; i != _exact_derivs.size(); ++i)
00584     if (_exact_derivs[i])
00585       _exact_derivs[i]->init();
00586 
00587   for (unsigned int i=0; i != _exact_hessians.size(); ++i)
00588     if (_exact_hessians[i])
00589       _exact_hessians[i]->init();
00590 
00591   // Get a reference to the dofmap and mesh for that system
00592   const DofMap& computed_dof_map = computed_system.get_dof_map();
00593 
00594   const MeshBase& _mesh = computed_system.get_mesh();
00595 
00596   const unsigned int dim = _mesh.mesh_dimension();
00597 
00598   // Zero the error before summation
00599   // 0 - sum of square of function error (L2)
00600   // 1 - sum of square of gradient error (H1 semi)
00601   // 2 - sum of square of Hessian error (H2 semi)
00602   // 3 - sum of sqrt(square of function error) (L1)
00603   // 4 - max of sqrt(square of function error) (Linfty)
00604   // 5 - sum of square of curl error (HCurl semi)
00605   // 6 - sum of square of div error (HDiv semi)
00606   error_vals = std::vector<Real>(7, 0.);
00607 
00608   // Construct Quadrature rule based on default quadrature order
00609   const FEType& fe_type  = computed_dof_map.variable_type(var);
00610 
00611   unsigned int n_vec_dim = FEInterface::n_vec_dim( _mesh, fe_type );
00612 
00613   // FIXME: MeshFunction needs to be updated to support vector-valued
00614   //        elements before we can use a reference solution.
00615   if( (n_vec_dim > 1) && _equation_systems_fine )
00616     {
00617       libMesh::err << "Error calculation using reference solution not yet\n"
00618                    << "supported for vector-valued elements."
00619                    << std::endl;
00620       libmesh_not_implemented();
00621     }
00622 
00623   AutoPtr<QBase> qrule =
00624     fe_type.default_quadrature_rule (dim,
00625                                      _extra_order);
00626 
00627   // Construct finite element object
00628 
00629   AutoPtr<FEGenericBase<OutputShape> > fe(FEGenericBase<OutputShape>::build(dim, fe_type));
00630 
00631   // Attach quadrature rule to FE object
00632   fe->attach_quadrature_rule (qrule.get());
00633 
00634   // The Jacobian*weight at the quadrature points.
00635   const std::vector<Real>& JxW                               = fe->get_JxW();
00636 
00637   // The value of the shape functions at the quadrature points
00638   // i.e. phi(i) = phi_values[i][qp]
00639   const std::vector<std::vector<OutputShape> >&  phi_values         = fe->get_phi();
00640 
00641   // The value of the shape function gradients at the quadrature points
00642   const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputGradient> >& 
00643                     dphi_values = fe->get_dphi();
00644 
00645   // The value of the shape function curls at the quadrature points
00646   // Only computed for vector-valued elements
00647   const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputShape> >* curl_values = NULL;
00648 
00649   // The value of the shape function divergences at the quadrature points
00650   // Only computed for vector-valued elements
00651   const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputDivergence> >* div_values = NULL;
00652 
00653   if( FEInterface::field_type(fe_type) == TYPE_VECTOR )
00654     {
00655       curl_values = &fe->get_curl_phi();
00656       div_values = &fe->get_div_phi();
00657     }
00658 
00659 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
00660   // The value of the shape function second derivatives at the quadrature points
00661   const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputTensor> >& 
00662                     d2phi_values = fe->get_d2phi(); 
00663 #endif
00664 
00665   // The XYZ locations (in physical space) of the quadrature points
00666   const std::vector<Point>& q_point                          = fe->get_xyz();
00667 
00668   // The global degree of freedom indices associated
00669   // with the local degrees of freedom.
00670   std::vector<dof_id_type> dof_indices;
00671 
00672 
00673   //
00674   // Begin the loop over the elements
00675   //
00676   // TODO: this ought to be threaded (and using subordinate
00677   // MeshFunction objects in each thread rather than a single
00678   // master)
00679   MeshBase::const_element_iterator       el     = _mesh.active_local_elements_begin();
00680   const MeshBase::const_element_iterator end_el = _mesh.active_local_elements_end();
00681 
00682   for ( ; el != end_el; ++el)
00683     {
00684       // Store a pointer to the element we are currently
00685       // working on.  This allows for nicer syntax later.
00686       const Elem* elem = *el;
00687 
00688       // reinitialize the element-specific data
00689       // for the current element
00690       fe->reinit (elem);
00691 
00692       // Get the local to global degree of freedom maps
00693       computed_dof_map.dof_indices    (elem, dof_indices, var);
00694 
00695       // The number of quadrature points
00696       const unsigned int n_qp = qrule->n_points();
00697 
00698       // The number of shape functions
00699       const unsigned int n_sf = 
00700         libmesh_cast_int<unsigned int>(dof_indices.size());
00701 
00702       //
00703       // Begin the loop over the Quadrature points.
00704       //
00705       for (unsigned int qp=0; qp<n_qp; qp++)
00706         {
00707           // Real u_h = 0.;
00708           // RealGradient grad_u_h;
00709 
00710           typename FEGenericBase<OutputShape>::OutputNumber u_h = 0.;
00711 
00712           typename FEGenericBase<OutputShape>::OutputNumberGradient grad_u_h;
00713 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
00714           typename FEGenericBase<OutputShape>::OutputNumberTensor grad2_u_h;
00715 #endif
00716           typename FEGenericBase<OutputShape>::OutputNumber curl_u_h = 0.0;
00717           typename FEGenericBase<OutputShape>::OutputNumberDivergence div_u_h = 0.0;
00718 
00719           // Compute solution values at the current
00720           // quadrature point.  This reqiures a sum
00721           // over all the shape functions evaluated
00722           // at the quadrature point.
00723           for (unsigned int i=0; i<n_sf; i++)
00724             {
00725               // Values from current solution.
00726               u_h      += phi_values[i][qp]*computed_system.current_solution  (dof_indices[i]);
00727               grad_u_h += dphi_values[i][qp]*computed_system.current_solution (dof_indices[i]);
00728 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
00729               grad2_u_h += d2phi_values[i][qp]*computed_system.current_solution (dof_indices[i]);
00730 #endif
00731               if( FEInterface::field_type(fe_type) == TYPE_VECTOR )
00732                 {
00733                   curl_u_h += (*curl_values)[i][qp]*computed_system.current_solution (dof_indices[i]);
00734                   div_u_h += (*div_values)[i][qp]*computed_system.current_solution (dof_indices[i]);
00735                 }
00736             }
00737 
00738           // Compute the value of the error at this quadrature point
00739           typename FEGenericBase<OutputShape>::OutputNumber exact_val = 0;
00740           RawAccessor<typename FEGenericBase<OutputShape>::OutputNumber> exact_val_accessor( exact_val, dim );
00741           if (_exact_values.size() > sys_num && _exact_values[sys_num])
00742             {
00743               for( unsigned int c = 0; c < n_vec_dim; c++)
00744                 exact_val_accessor(c) =
00745                   _exact_values[sys_num]->
00746                   component(var_component+c, q_point[qp], time);
00747             }
00748           else if (_equation_systems_fine)
00749             {
00750               // FIXME: Needs to be updated for vector-valued elements
00751               exact_val = (*coarse_values)(q_point[qp]);
00752             }
00753           const typename FEGenericBase<OutputShape>::OutputNumber val_error = u_h - exact_val;
00754 
00755           // Add the squares of the error to each contribution
00756           Real error_sq = TensorTools::norm_sq(val_error);
00757           error_vals[0] += JxW[qp]*error_sq;
00758 
00759           Real norm = sqrt(error_sq);
00760           error_vals[3] += JxW[qp]*norm;
00761 
00762           if(error_vals[4]<norm) { error_vals[4] = norm; }
00763 
00764           // Compute the value of the error in the gradient at this
00765           // quadrature point
00766           typename FEGenericBase<OutputShape>::OutputNumberGradient exact_grad;
00767           RawAccessor<typename FEGenericBase<OutputShape>::OutputNumberGradient> exact_grad_accessor( exact_grad, _mesh.spatial_dimension() );
00768           if (_exact_derivs.size() > sys_num && _exact_derivs[sys_num])
00769             {
00770               for( unsigned int c = 0; c < n_vec_dim; c++)
00771                 for( unsigned int d = 0; d < _mesh.spatial_dimension(); d++ )
00772                   exact_grad_accessor(d + c*_mesh.spatial_dimension() ) =
00773                     _exact_derivs[sys_num]->
00774                     component(var_component+c, q_point[qp], time)(d);
00775             }
00776           else if (_equation_systems_fine)
00777             {
00778               // FIXME: Needs to be updated for vector-valued elements
00779               exact_grad = coarse_values->gradient(q_point[qp]);
00780             }
00781 
00782           const typename FEGenericBase<OutputShape>::OutputNumberGradient grad_error = grad_u_h - exact_grad;
00783           
00784           error_vals[1] += JxW[qp]*grad_error.size_sq();
00785 
00786 
00787           if( FEInterface::field_type(fe_type) == TYPE_VECTOR )
00788             {
00789               // Compute the value of the error in the curl at this
00790               // quadrature point
00791               typename FEGenericBase<OutputShape>::OutputNumber exact_curl = 0.0;
00792               if (_exact_derivs.size() > sys_num && _exact_derivs[sys_num])
00793                 {
00794                   exact_curl = TensorTools::curl_from_grad( exact_grad );
00795                 }
00796               else if (_equation_systems_fine)
00797                 {
00798                   // FIXME: Need to implement curl for MeshFunction and support reference
00799                   //        solution for vector-valued elements
00800                 }
00801 
00802               const typename FEGenericBase<OutputShape>::OutputNumber curl_error = curl_u_h - exact_curl;
00803               
00804               error_vals[5] += JxW[qp]*TensorTools::norm_sq(curl_error);
00805 
00806               // Compute the value of the error in the divergence at this
00807               // quadrature point
00808               typename FEGenericBase<OutputShape>::OutputNumberDivergence exact_div = 0.0;
00809               if (_exact_derivs.size() > sys_num && _exact_derivs[sys_num])
00810                 {
00811                   exact_div = TensorTools::div_from_grad( exact_grad );
00812                 }
00813               else if (_equation_systems_fine)
00814                 {
00815                   // FIXME: Need to implement div for MeshFunction and support reference
00816                   //        solution for vector-valued elements
00817                 }
00818               
00819               const typename FEGenericBase<OutputShape>::OutputNumberDivergence div_error = div_u_h - exact_div;
00820               
00821               error_vals[6] += JxW[qp]*TensorTools::norm_sq(div_error);
00822             }
00823 
00824 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
00825           // Compute the value of the error in the hessian at this
00826           // quadrature point
00827           typename FEGenericBase<OutputShape>::OutputNumberTensor exact_hess;
00828           RawAccessor<typename FEGenericBase<OutputShape>::OutputNumberTensor> exact_hess_accessor( exact_hess, dim );
00829           if (_exact_hessians.size() > sys_num && _exact_hessians[sys_num])
00830             {
00831               //FIXME: This needs to be implemented to support rank 3 tensors 
00832               //       which can't happen until type_n_tensor is fully implemented
00833               //       and a RawAccessor<TypeNTensor> is fully implemented
00834               if( FEInterface::field_type(fe_type) == TYPE_VECTOR )
00835                 libmesh_not_implemented();
00836 
00837               for( unsigned int c = 0; c < n_vec_dim; c++)
00838                 for( unsigned int d = 0; d < dim; d++ )
00839                   for( unsigned int e =0; e < dim; e++ )
00840                     exact_hess_accessor(d + e*dim + c*dim*dim) =
00841                       _exact_hessians[sys_num]->
00842                       component(var_component+c, q_point[qp], time)(d,e);
00843             }
00844           else if (_equation_systems_fine)
00845             {
00846               // FIXME: Needs to be updated for vector-valued elements
00847               exact_hess = coarse_values->hessian(q_point[qp]);
00848             }
00849           
00850           const typename FEGenericBase<OutputShape>::OutputNumberTensor grad2_error = grad2_u_h - exact_hess;
00851           
00852           // FIXME: PB: Is this what we want for rank 3 tensors?
00853           error_vals[2] += JxW[qp]*grad2_error.size_sq();
00854 #endif
00855 
00856         } // end qp loop
00857     } // end element loop
00858 
00859   // Add up the error values on all processors, except for the L-infty
00860   // norm, for which the maximum is computed.
00861   Real l_infty_norm = error_vals[4];
00862   CommWorld.max(l_infty_norm);
00863   CommWorld.sum(error_vals);
00864   error_vals[4] = l_infty_norm;
00865 }

void libMesh::ExactSolution::attach_exact_deriv ( Gradient   gptrconst Point &p,const Parameters &parameters,const std::string &sys_name,const std::string &unknown_name  ) 

Attach an arbitrary function which computes the exact gradient of the solution at any point.

Definition at line 153 of file exact_solution.C.

References _equation_systems, _equation_systems_fine, _exact_derivs, libMesh::EquationSystems::get_system(), libMesh::EquationSystems::n_systems(), and libMesh::EquationSystems::parameters.

00157 {
00158   libmesh_assert(gptr);
00159 
00160   // Clear out any previous _exact_derivs entries, then add a new
00161   // entry for each system.
00162   _exact_derivs.clear();
00163 
00164   for (unsigned int sys=0; sys<_equation_systems.n_systems(); ++sys)
00165     {
00166       const System& system = _equation_systems.get_system(sys);
00167       _exact_derivs.push_back
00168         (new WrappedFunction<Gradient>
00169           (system, gptr, &_equation_systems.parameters));
00170     }
00171 
00172   // If we're using exact values, we're not using a fine grid solution
00173   _equation_systems_fine = NULL;
00174 }

void libMesh::ExactSolution::attach_exact_deriv ( unsigned int  sys_num,
FunctionBase< Gradient > *  g 
)

Clone and attach an arbitrary functor which computes the exact gradient of the system sys_num solution at any point.

Definition at line 196 of file exact_solution.C.

References _exact_derivs, libMesh::FunctionBase< Output >::clone(), and libMesh::AutoPtr< Tp >::release().

00198 {
00199   if (_exact_derivs.size() <= sys_num)
00200   _exact_derivs.resize(sys_num+1, NULL);
00201 
00202   if (g)
00203     _exact_derivs[sys_num] = g->clone().release();
00204 }

void libMesh::ExactSolution::attach_exact_derivs ( std::vector< FunctionBase< Gradient > * >  g  ) 

Clone and attach arbitrary functors which compute the exact gradients of the EquationSystems' solutions at any point.

Definition at line 177 of file exact_solution.C.

References _exact_derivs.

00178 {
00179   // Clear out any previous _exact_derivs entries, then add a new
00180   // entry for each system.
00181   for (unsigned int i=0; i != _exact_derivs.size(); ++i)
00182     delete (_exact_derivs[i]);
00183 
00184   _exact_derivs.clear();
00185   _exact_derivs.resize(g.size(), NULL);
00186 
00187   // We use clone() to get non-sliced copies of FunctionBase
00188   // subclasses, but we can't put the resulting AutoPtrs into an STL
00189   // container.
00190   for (unsigned int i=0; i != g.size(); ++i)
00191     if (g[i])
00192       _exact_derivs[i] = g[i]->clone().release();
00193 }

void libMesh::ExactSolution::attach_exact_hessian ( Tensor   hptrconst Point &p,const Parameters &parameters,const std::string &sys_name,const std::string &unknown_name  ) 

Attach an arbitrary function which computes the exact second derivatives of the solution at any point.

Definition at line 207 of file exact_solution.C.

References _equation_systems, _equation_systems_fine, _exact_hessians, libMesh::EquationSystems::get_system(), libMesh::EquationSystems::n_systems(), and libMesh::EquationSystems::parameters.

00211 {
00212   libmesh_assert(hptr);
00213 
00214   // Clear out any previous _exact_hessians entries, then add a new
00215   // entry for each system.
00216   _exact_hessians.clear();
00217 
00218   for (unsigned int sys=0; sys<_equation_systems.n_systems(); ++sys)
00219     {
00220       const System& system = _equation_systems.get_system(sys);
00221       _exact_hessians.push_back
00222         (new WrappedFunction<Tensor>
00223           (system, hptr, &_equation_systems.parameters));
00224     }
00225 
00226   // If we're using exact values, we're not using a fine grid solution
00227   _equation_systems_fine = NULL;
00228 }

void libMesh::ExactSolution::attach_exact_hessian ( unsigned int  sys_num,
FunctionBase< Tensor > *  h 
)

Clone and attach an arbitrary functor which computes the exact second derivatives of the system sys_num solution at any point.

Definition at line 250 of file exact_solution.C.

References _exact_hessians, and libMesh::FunctionBase< Output >::clone().

00252 {
00253   if (_exact_hessians.size() <= sys_num)
00254   _exact_hessians.resize(sys_num+1, NULL);
00255 
00256   if (h)
00257     _exact_hessians[sys_num] = h->clone().release();
00258 }

void libMesh::ExactSolution::attach_exact_hessians ( std::vector< FunctionBase< Tensor > * >  h  ) 

Clone and attach arbitrary functors which compute the exact second derivatives of the EquationSystems' solutions at any point.

Definition at line 231 of file exact_solution.C.

References _exact_hessians.

00232 {
00233   // Clear out any previous _exact_hessians entries, then add a new
00234   // entry for each system.
00235   for (unsigned int i=0; i != _exact_hessians.size(); ++i)
00236     delete (_exact_hessians[i]);
00237 
00238   _exact_hessians.clear();
00239   _exact_hessians.resize(h.size(), NULL);
00240 
00241   // We use clone() to get non-sliced copies of FunctionBase
00242   // subclasses, but we can't put the resulting AutoPtrs into an STL
00243   // container.
00244   for (unsigned int i=0; i != h.size(); ++i)
00245     if (h[i])
00246       _exact_hessians[i] = h[i]->clone().release();
00247 }

void libMesh::ExactSolution::attach_exact_value ( Number   fptrconst Point &p, const Parameters &Parameters, const std::string &sys_name, const std::string &unknown_name  ) 

Attach an arbitrary function which computes the exact value of the solution at any point.

void libMesh::ExactSolution::attach_exact_value ( unsigned int  sys_num,
FunctionBase< Number > *  f 
)

Clone and attach an arbitrary functor which computes the exact value of the system sys_num solution at any point.

Definition at line 142 of file exact_solution.C.

References _exact_values, libMesh::FunctionBase< Output >::clone(), and libMesh::AutoPtr< Tp >::release().

00144 {
00145   if (_exact_values.size() <= sys_num)
00146   _exact_values.resize(sys_num+1, NULL);
00147 
00148   if (f)
00149     _exact_values[sys_num] = f->clone().release();
00150 }

void libMesh::ExactSolution::attach_exact_values ( std::vector< FunctionBase< Number > * >  f  ) 

Clone and attach arbitrary functors which compute the exact values of the EquationSystems' solutions at any point.

Definition at line 123 of file exact_solution.C.

References _exact_values.

00124 {
00125   // Clear out any previous _exact_values entries, then add a new
00126   // entry for each system.
00127   for (unsigned int i=0; i != _exact_values.size(); ++i)
00128     delete (_exact_values[i]);
00129 
00130   _exact_values.clear();
00131   _exact_values.resize(f.size(), NULL);
00132 
00133   // We use clone() to get non-sliced copies of FunctionBase
00134   // subclasses, but we can't put the resulting AutoPtrs into an STL
00135   // container.
00136   for (unsigned int i=0; i != f.size(); ++i)
00137     if (f[i])
00138       _exact_values[i] = f[i]->clone().release();
00139 }

void libMesh::ExactSolution::attach_reference_solution ( const EquationSystems es_fine  ) 

Attach function similar to system.h which allows the user to attach a second EquationSystems object with a reference fine grid solution.

Definition at line 87 of file exact_solution.C.

References _equation_systems_fine, _exact_derivs, _exact_hessians, and _exact_values.

00088 {
00089   libmesh_assert(es_fine);
00090   _equation_systems_fine = es_fine;
00091 
00092   // If we're using a fine grid solution, we're not using exact values
00093   _exact_values.clear();
00094   _exact_derivs.clear();
00095   _exact_hessians.clear();
00096 }

void libMesh::ExactSolution::compute_error ( const std::string &  sys_name,
const std::string &  unknown_name 
)

Computes and stores the error in the solution value e = u-u_h, the gradient grad(e) = grad(u) - grad(u_h), and possibly the hessian grad(grad(e)) = grad(grad(u)) - grad(grad(u_h)). Does not return any value. For that you need to call the l2_error(), h1_error() or h2_error() functions respectively.

Definition at line 297 of file exact_solution.C.

References _check_inputs(), _equation_systems, libMesh::FEInterface::field_type(), libMesh::EquationSystems::get_system(), libMesh::EquationSystems::has_system(), libMeshEnums::TYPE_SCALAR, and libMeshEnums::TYPE_VECTOR.

00299 {
00300   // Check the inputs for validity, and get a reference
00301   // to the proper location to store the error
00302   std::vector<Real>& error_vals = this->_check_inputs(sys_name,
00303                                                       unknown_name);
00304 
00305   libmesh_assert( _equation_systems.has_system(sys_name) );
00306   const System& sys = _equation_systems.get_system<System>( sys_name );
00307   
00308   libmesh_assert( sys.has_variable( unknown_name ) );
00309   switch( FEInterface::field_type(sys.variable_type( unknown_name )) )
00310     {
00311     case TYPE_SCALAR:
00312       {
00313         this->_compute_error<Real>(sys_name,
00314                                    unknown_name,
00315                                    error_vals);
00316         break;
00317       }
00318     case TYPE_VECTOR:
00319       {
00320         this->_compute_error<RealGradient>(sys_name,
00321                                            unknown_name,
00322                                            error_vals);
00323         break;
00324       }
00325     default:
00326       libmesh_error();
00327     }
00328 
00329   return;
00330 }

Real libMesh::ExactSolution::error_norm ( const std::string &  sys_name,
const std::string &  unknown_name,
const FEMNormType &  norm 
)

This function returns the error in the requested norm for the system sys_name for the unknown unknown_name. Note that no error computations are actually performed, you must call compute_error() for that. Note also that the result is not exact, but an approximation based on the chosen quadrature rule.

Definition at line 336 of file exact_solution.C.

References _check_inputs(), _equation_systems, libMesh::err, libMesh::FEInterface::field_type(), libMesh::EquationSystems::get_system(), libMeshEnums::H1, libMeshEnums::H1_SEMINORM, libMeshEnums::H2, libMeshEnums::H2_SEMINORM, libMesh::EquationSystems::has_system(), libMeshEnums::HCURL, libMeshEnums::HCURL_SEMINORM, libMeshEnums::HDIV, libMeshEnums::HDIV_SEMINORM, libMeshEnums::L1, libMeshEnums::L2, libMeshEnums::L_INF, and libMeshEnums::TYPE_SCALAR.

Referenced by hcurl_error(), and hdiv_error().

00339 {
00340   // Check the inputs for validity, and get a reference
00341   // to the proper location to store the error
00342   std::vector<Real>& error_vals = this->_check_inputs(sys_name,
00343                                                       unknown_name);
00344 
00345   libmesh_assert(_equation_systems.has_system(sys_name));
00346   libmesh_assert(_equation_systems.get_system(sys_name).has_variable( unknown_name ));
00347   const FEType& fe_type = _equation_systems.get_system(sys_name).variable_type(unknown_name);
00348 
00349   switch (norm)
00350     {
00351     case L2:
00352       return std::sqrt(error_vals[0]);
00353     case H1:
00354       return std::sqrt(error_vals[0] + error_vals[1]);
00355     case H2:
00356       return std::sqrt(error_vals[0] + error_vals[1] + error_vals[2]);
00357     case HCURL:
00358       {
00359         if(FEInterface::field_type(fe_type) == TYPE_SCALAR)
00360           {
00361             libMesh::err << "Cannot compute HCurl error norm of scalar-valued variables!" << std::endl;
00362             libmesh_error();
00363           }
00364         else
00365           return std::sqrt(error_vals[0] + error_vals[5]);
00366       }
00367     case HDIV:
00368       {
00369         if(FEInterface::field_type(fe_type) == TYPE_SCALAR)
00370           {
00371             libMesh::err << "Cannot compute HDiv error norm of scalar-valued variables!" << std::endl;
00372             libmesh_error();
00373           }
00374         else
00375           return std::sqrt(error_vals[0] + error_vals[6]);
00376       }
00377     case H1_SEMINORM:
00378       return std::sqrt(error_vals[1]);
00379     case H2_SEMINORM:
00380       return std::sqrt(error_vals[2]);
00381     case HCURL_SEMINORM:
00382       {
00383         if(FEInterface::field_type(fe_type) == TYPE_SCALAR)
00384           {
00385             libMesh::err << "Cannot compute HCurl error seminorm of scalar-valued variables!" << std::endl;
00386             libmesh_error();
00387           }
00388         else
00389           return std::sqrt(error_vals[5]);
00390       }
00391     case HDIV_SEMINORM:
00392       {
00393         if(FEInterface::field_type(fe_type) == TYPE_SCALAR)
00394           {
00395             libMesh::err << "Cannot compute HDiv error seminorm of scalar-valued variables!" << std::endl;
00396             libmesh_error();
00397           }
00398         else
00399           return std::sqrt(error_vals[6]);
00400       }
00401     case L1:
00402       return error_vals[3];
00403     case L_INF:
00404       return error_vals[4];
00405     
00406     // Currently only Sobolev norms/seminorms are supported
00407     default:
00408       libmesh_error();
00409     }
00410 }

void libMesh::ExactSolution::extra_quadrature_order ( const int  extraorder  )  [inline]

Increases or decreases the order of the quadrature rule used for numerical integration.

Definition at line 159 of file exact_solution.h.

References _extra_order.

00160     { _extra_order = extraorder; }

Real libMesh::ExactSolution::h1_error ( const std::string &  sys_name,
const std::string &  unknown_name 
)

This function computes and returns the H1 error for the system sys_name for the unknown unknown_name. Note that no error computations are actually performed, you must call compute_error() for that.

Definition at line 478 of file exact_solution.C.

References _check_inputs().

00480 {
00481   // If the user has supplied no exact derivative function, we
00482   // just integrate the H1 norm of the solution; i.e. its
00483   // difference from an "exact solution" of zero.
00484 
00485   // Check the inputs for validity, and get a reference
00486   // to the proper location to store the error
00487   std::vector<Real>& error_vals = this->_check_inputs(sys_name,
00488                                                       unknown_name);
00489 
00490   // Return the square root of the sum of the computed errors.
00491   return std::sqrt(error_vals[0] + error_vals[1]);
00492 }

Real libMesh::ExactSolution::h2_error ( const std::string &  sys_name,
const std::string &  unknown_name 
)

This function computes and returns the H2 error for the system sys_name for the unknown unknown_name. Note that no error computations are actually performed, you must call compute_error() for that.

Definition at line 510 of file exact_solution.C.

References _check_inputs().

00512 {
00513   // If the user has supplied no exact derivative functions, we
00514   // just integrate the H2 norm of the solution; i.e. its
00515   // difference from an "exact solution" of zero.
00516 
00517   // Check the inputs for validity, and get a reference
00518   // to the proper location to store the error
00519   std::vector<Real>& error_vals = this->_check_inputs(sys_name,
00520                                                       unknown_name);
00521 
00522   // Return the square root of the sum of the computed errors.
00523   return std::sqrt(error_vals[0] + error_vals[1] + error_vals[2]);
00524 }

Real libMesh::ExactSolution::hcurl_error ( const std::string &  sys_name,
const std::string &  unknown_name 
)

This function computes and returns the HCurl error for the system sys_name for the unknown unknown_name. Note that no error computations are actually performed, you must call compute_error() for that. This is only valid for vector-valued element. An error is thrown if requested for scalar-valued elements.

Definition at line 495 of file exact_solution.C.

References error_norm(), and libMeshEnums::HCURL.

00497 {
00498   return this->error_norm(sys_name,unknown_name,HCURL);
00499 }

Real libMesh::ExactSolution::hdiv_error ( const std::string &  sys_name,
const std::string &  unknown_name 
)

This function computes and returns the HDiv error for the system sys_name for the unknown unknown_name. Note that no error computations are actually performed, you must call compute_error() for that. This is only valid for vector-valued element. An error is thrown if requested for scalar-valued elements.

Definition at line 502 of file exact_solution.C.

References error_norm(), and libMeshEnums::HDIV.

00504 {
00505   return this->error_norm(sys_name,unknown_name,HDIV);
00506 }

Real libMesh::ExactSolution::l1_error ( const std::string &  sys_name,
const std::string &  unknown_name 
)

This function returns the integrated L1 error for the system sys_name for the unknown unknown_name. Note that no error computations are actually performed, you must call compute_error() for that.

Definition at line 438 of file exact_solution.C.

References _check_inputs().

00440 {
00441 
00442   // Check the inputs for validity, and get a reference
00443   // to the proper location to store the error
00444   std::vector<Real>& error_vals = this->_check_inputs(sys_name,
00445                                                       unknown_name);
00446 
00447   // Return the square root of the first component of the
00448   // computed error.
00449   return error_vals[3];
00450 }

Real libMesh::ExactSolution::l2_error ( const std::string &  sys_name,
const std::string &  unknown_name 
)

This function returns the integrated L2 error for the system sys_name for the unknown unknown_name. Note that no error computations are actually performed, you must call compute_error() for that.

Definition at line 418 of file exact_solution.C.

References _check_inputs().

00420 {
00421 
00422   // Check the inputs for validity, and get a reference
00423   // to the proper location to store the error
00424   std::vector<Real>& error_vals = this->_check_inputs(sys_name,
00425                                                       unknown_name);
00426 
00427   // Return the square root of the first component of the
00428   // computed error.
00429   return std::sqrt(error_vals[0]);
00430 }

Real libMesh::ExactSolution::l_inf_error ( const std::string &  sys_name,
const std::string &  unknown_name 
)

This function returns the L_INF error for the system sys_name for the unknown unknown_name. Note that no error computations are actually performed, you must call compute_error() for that. Note also that the result (as for the other norms as well) is not exact, but an approximation based on the chosen quadrature rule: to compute it, we take the max of the absolute value of the error over all the quadrature points.

Definition at line 458 of file exact_solution.C.

References _check_inputs().

00460 {
00461 
00462   // Check the inputs for validity, and get a reference
00463   // to the proper location to store the error
00464   std::vector<Real>& error_vals = this->_check_inputs(sys_name,
00465                                                       unknown_name);
00466 
00467   // Return the square root of the first component of the
00468   // computed error.
00469   return error_vals[4];
00470 }


Member Data Documentation

Constant reference to the EquationSystems object used for the simulation.

Definition at line 307 of file exact_solution.h.

Referenced by _compute_error(), attach_exact_deriv(), attach_exact_hessian(), compute_error(), error_norm(), and ExactSolution().

Constant pointer to the EquationSystems object containing the fine grid solution.

Definition at line 313 of file exact_solution.h.

Referenced by _compute_error(), attach_exact_deriv(), attach_exact_hessian(), and attach_reference_solution().

std::map<std::string, SystemErrorMap> libMesh::ExactSolution::_errors [private]

A map of SystemErrorMaps, which contains entries for each system in the EquationSystems object. This is required, since it is possible for two systems to have unknowns with the *same name*.

Definition at line 301 of file exact_solution.h.

Referenced by _check_inputs(), and ExactSolution().

User-provided functors which compute the exact derivative of the solution for each system.

Definition at line 277 of file exact_solution.h.

Referenced by _compute_error(), attach_exact_deriv(), attach_exact_derivs(), attach_reference_solution(), and ~ExactSolution().

User-provided functors which compute the exact hessians of the solution for each system.

Definition at line 283 of file exact_solution.h.

Referenced by _compute_error(), attach_exact_hessian(), attach_exact_hessians(), attach_reference_solution(), and ~ExactSolution().

User-provided functors which compute the exact value of the solution for each system.

Definition at line 271 of file exact_solution.h.

Referenced by _compute_error(), attach_exact_value(), attach_exact_values(), attach_reference_solution(), and ~ExactSolution().

Extra order to use for quadrature rule

Definition at line 318 of file exact_solution.h.

Referenced by _compute_error(), and extra_quadrature_order().


The documentation for this class was generated from the following files:

Site Created By: libMesh Developers
Last modified: February 05 2013 19:55:18 UTC

Hosted By:
SourceForge.net Logo