libMesh::ExactSolution Class Reference

This class handles the computation of the L2 and/or H1 error for the Systems in the EquationSystems object which is passed to it. More...

#include <exact_solution.h>

Public Types
typedef Number(*	ValueFunctionPointer) (const 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. More...
typedef Gradient(*	GradientFunctionPointer) (const 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. More...
typedef Tensor(*	HessianFunctionPointer) (const 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. More...

Public Member Functions
	ExactSolution (const EquationSystems &es)
	Constructor. More...
	ExactSolution (const ExactSolution &)=delete
	The copy constructor and copy/move assignment operators are deleted. More...
ExactSolution &	operator= (const ExactSolution &)=delete
ExactSolution &	operator= (ExactSolution &&)=delete
	ExactSolution (ExactSolution &&)
	Move constructor and destructor are defaulted out-of-line (in the C file) to play nicely with our forward declarations. More...
	~ExactSolution ()
void	set_excluded_subdomains (const std::set< subdomain_id_type > &excluded)
	The user can indicate that elements in certain subdomains should be excluded from the error calculation by passing in a set of subdomain ids to ignore. More...
void	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. More...
void	attach_exact_values (const std::vector< FunctionBase< Number > * > &f)
	Clone and attach arbitrary functors which compute the exact values of the EquationSystems' solutions at any point. More...
void	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. More...
void	attach_exact_value (ValueFunctionPointer fptr)
void	attach_exact_derivs (const std::vector< FunctionBase< Gradient > * > &g)
	Clone and attach arbitrary functors which compute the exact gradients of the EquationSystems' solutions at any point. More...
void	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. More...
void	attach_exact_deriv (GradientFunctionPointer gptr)
void	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. More...
void	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. More...
void	attach_exact_hessian (HessianFunctionPointer hptr)
void	extra_quadrature_order (const int extraorder)
	Increases or decreases the order of the quadrature rule used for numerical integration. More...
void	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)). More...
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
	Data structure which stores the errors: ||e|| = ||u - u_h|| ||grad(e)|| = ||grad(u) - grad(u_h)|| for each unknown in a single system. More...

Private Member Functions
template<typename OutputShape >
void	_compute_error (const std::string &sys_name, const std::string &unknown_name, std::vector< Real > &error_vals)
	This function computes the error (in the solution and its first derivative) for a single unknown in a single system. More...
std::vector< Real > &	_check_inputs (const std::string &sys_name, const std::string &unknown_name)
	This function is responsible for checking the validity of the sys_name and unknown_name inputs. More...

Private Attributes
std::vector< std::unique_ptr< FunctionBase< Number > > >	_exact_values
	User-provided functors which compute the exact value of the solution for each system. More...
std::vector< std::unique_ptr< FunctionBase< Gradient > > >	_exact_derivs
	User-provided functors which compute the exact derivative of the solution for each system. More...
std::vector< std::unique_ptr< FunctionBase< Tensor > > >	_exact_hessians
	User-provided functors which compute the exact hessians of the solution for each system. More...
std::map< std::string, SystemErrorMap >	_errors
	A map of SystemErrorMaps, which contains entries for each system in the EquationSystems object. More...
const EquationSystems &	_equation_systems
	Constant reference to the EquationSystems object used for the simulation. More...
const EquationSystems *	_equation_systems_fine
	Constant pointer to the EquationSystems object containing the fine grid solution. More...
int	_extra_order
	Extra order to use for quadrature rule. More...
std::set< subdomain_id_type >	_excluded_subdomains
	Elements in a subdomain from this set are skipped during the error computation. More...

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

For this 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

John W. Peterson (modifications for libmesh)

Date

2004

Definition at line 73 of file exact_solution.h.

Member Typedef Documentation

GradientFunctionPointer

typedef Gradient(* libMesh::ExactSolution::GradientFunctionPointer) (const 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 157 of file exact_solution.h.

HessianFunctionPointer

typedef Tensor(* libMesh::ExactSolution::HessianFunctionPointer) (const 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 180 of file exact_solution.h.

SystemErrorMap

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 345 of file exact_solution.h.

ValueFunctionPointer

typedef Number(* libMesh::ExactSolution::ValueFunctionPointer) (const 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.

Definition at line 134 of file exact_solution.h.

Constructor & Destructor Documentation

ExactSolution() [1/3]

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.

                                                       :
  _equation_systems(es),
  _equation_systems_fine(nullptr),
  _extra_order(0)
{
  // Initialize the _errors data structure which holds all
  // the eventual values of the error.
  for (auto sys : IntRange<unsigned int>(0, _equation_systems.n_systems()))
    {
      // Reference to the system
      const System & system = _equation_systems.get_system(sys);
 
      // The name of the system
      const std::string & sys_name = system.name();
 
      // The SystemErrorMap to be inserted
      ExactSolution::SystemErrorMap sem;
 
      for (auto var : IntRange<unsigned int>(0, system.n_vars()))
        {
          // The name of this variable
          const std::string & var_name = system.variable_name(var);
          sem[var_name] = std::vector<Real>(5, 0.);
        }
 
      _errors[sys_name] = sem;
    }
}

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

ExactSolution() [2/3]

libMesh::ExactSolution::ExactSolution ( const ExactSolution &  )

delete

The copy constructor and copy/move assignment operators are deleted.

This class has containers of unique_ptrs so it can't be default (shallow) copied, and it has a const reference so it can't be assigned to after creation.

ExactSolution() [3/3]

libMesh::ExactSolution::ExactSolution ( ExactSolution &&  )

default

Move constructor and destructor are defaulted out-of-line (in the C file) to play nicely with our forward declarations.

~ExactSolution()

libMesh::ExactSolution::~ExactSolution ( )

default

Member Function Documentation

_check_inputs()

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.

Returns

A reference to the proper vector for storing the values.

Definition at line 216 of file exact_solution.C.

{
  // Return a reference to the proper error entry, or throw an error
  // if it doesn't exist.
  auto & system_error_map = libmesh_map_find(_errors, sys_name);
  return libmesh_map_find(system_error_map, unknown_name);
}

References _errors.

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

_compute_error()

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  )

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 448 of file exact_solution.C.

{
  // Make sure we aren't "overconfigured"
  libmesh_assert (!(_exact_values.size() && _equation_systems_fine));
 
  // We need a communicator.
  const Parallel::Communicator & communicator(_equation_systems.comm());
 
  // This function must be run on all processors at once
  libmesh_parallel_only(communicator);
 
  // Get a reference to the system whose error is being computed.
  // If we have a fine grid, however, we'll integrate on that instead
  // for more accuracy.
  const System & computed_system = _equation_systems_fine ?
    _equation_systems_fine->get_system(sys_name) :
    _equation_systems.get_system (sys_name);
 
  const Real time = _equation_systems.get_system(sys_name).time;
 
  const unsigned int sys_num = computed_system.number();
  const unsigned int var = computed_system.variable_number(unknown_name);
  const unsigned int var_component =
    computed_system.variable_scalar_number(var, 0);
 
  // Prepare a global solution and a MeshFunction of the coarse system if we need one
  std::unique_ptr<MeshFunction> coarse_values;
  std::unique_ptr<NumericVector<Number>> comparison_soln = NumericVector<Number>::build(_equation_systems.comm());
  if (_equation_systems_fine)
    {
      const System & comparison_system
        = _equation_systems.get_system(sys_name);
 
      std::vector<Number> global_soln;
      comparison_system.update_global_solution(global_soln);
      comparison_soln->init(comparison_system.solution->size(), true, SERIAL);
      (*comparison_soln) = global_soln;
 
      coarse_values = libmesh_make_unique<MeshFunction>
        (_equation_systems,
         *comparison_soln,
         comparison_system.get_dof_map(),
         comparison_system.variable_number(unknown_name));
      coarse_values->init();
    }
 
  // Initialize any functors we're going to use
  for (auto & ev : _exact_values)
    if (ev)
      ev->init();
 
  for (auto & ed : _exact_derivs)
    if (ed)
      ed->init();
 
  for (auto & eh : _exact_hessians)
    if (eh)
      eh->init();
 
  // Get a reference to the dofmap and mesh for that system
  const DofMap & computed_dof_map = computed_system.get_dof_map();
 
  const MeshBase & mesh = computed_system.get_mesh();
 
  // Grab which element dimensions are present in the mesh
  const std::set<unsigned char> & elem_dims = mesh.elem_dimensions();
 
  // Zero the error before summation
  // 0 - sum of square of function error (L2)
  // 1 - sum of square of gradient error (H1 semi)
  // 2 - sum of square of Hessian error (H2 semi)
  // 3 - sum of sqrt(square of function error) (L1)
  // 4 - max of sqrt(square of function error) (Linfty)
  // 5 - sum of square of curl error (HCurl semi)
  // 6 - sum of square of div error (HDiv semi)
  error_vals = std::vector<Real>(7, 0.);
 
  // Construct Quadrature rule based on default quadrature order
  const FEType & fe_type  = computed_dof_map.variable_type(var);
 
  unsigned int n_vec_dim = FEInterface::n_vec_dim( mesh, fe_type );
 
  // FIXME: MeshFunction needs to be updated to support vector-valued
  //        elements before we can use a reference solution.
  if ((n_vec_dim > 1) && _equation_systems_fine)
    {
      libMesh::err << "Error calculation using reference solution not yet\n"
                   << "supported for vector-valued elements."
                   << std::endl;
      libmesh_not_implemented();
    }
 
 
  // Allow space for dims 0-3, even if we don't use them all
  std::vector<std::unique_ptr<FEGenericBase<OutputShape>>> fe_ptrs(4);
  std::vector<std::unique_ptr<QBase>> q_rules(4);
 
  // Prepare finite elements for each dimension present in the mesh
  for (const auto dim : elem_dims)
    {
      // Build a quadrature rule.
      q_rules[dim] = fe_type.default_quadrature_rule (dim, _extra_order);
 
      // Construct finite element object
      fe_ptrs[dim] = FEGenericBase<OutputShape>::build(dim, fe_type);
 
      // Attach quadrature rule to FE object
      fe_ptrs[dim]->attach_quadrature_rule (q_rules[dim].get());
    }
 
  // The global degree of freedom indices associated
  // with the local degrees of freedom.
  std::vector<dof_id_type> dof_indices;
 
 
  //
  // Begin the loop over the elements
  //
  // TODO: this ought to be threaded (and using subordinate
  // MeshFunction objects in each thread rather than a single
  // master)
  for (const auto & elem : mesh.active_local_element_ptr_range())
    {
      // Skip this element if it is in a subdomain excluded by the user.
      const subdomain_id_type elem_subid = elem->subdomain_id();
      if (_excluded_subdomains.count(elem_subid))
        continue;
 
      // The spatial dimension of the current Elem. FEs and other data
      // are indexed on dim.
      const unsigned int dim = elem->dim();
 
      // If the variable is not active on this subdomain, don't bother
      if (!computed_system.variable(var).active_on_subdomain(elem_subid))
        continue;
 
      /* If the variable is active, then we're going to restrict the
         MeshFunction evaluations to the current element subdomain.
         This is for cases such as mixed dimension meshes where we want
         to restrict the calculation to one particular domain. */
      std::set<subdomain_id_type> subdomain_id;
      subdomain_id.insert(elem_subid);
 
      FEGenericBase<OutputShape> * fe = fe_ptrs[dim].get();
      QBase * qrule = q_rules[dim].get();
      libmesh_assert(fe);
      libmesh_assert(qrule);
 
      // The Jacobian*weight at the quadrature points.
      const std::vector<Real> & JxW = fe->get_JxW();
 
      // The value of the shape functions at the quadrature points
      // i.e. phi(i) = phi_values[i][qp]
      const std::vector<std::vector<OutputShape>> &  phi_values = fe->get_phi();
 
      // The value of the shape function gradients at the quadrature points
      const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputGradient>> &
        dphi_values = fe->get_dphi();
 
      // The value of the shape function curls at the quadrature points
      // Only computed for vector-valued elements
      const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputShape>> * curl_values = nullptr;
 
      // The value of the shape function divergences at the quadrature points
      // Only computed for vector-valued elements
      const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputDivergence>> * div_values = nullptr;
 
      if (FEInterface::field_type(fe_type) == TYPE_VECTOR)
        {
          curl_values = &fe->get_curl_phi();
          div_values = &fe->get_div_phi();
        }
 
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
      // The value of the shape function second derivatives at the quadrature points
      const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputTensor>> &
        d2phi_values = fe->get_d2phi();
#endif
 
      // The XYZ locations (in physical space) of the quadrature points
      const std::vector<Point> & q_point = fe->get_xyz();
 
      // reinitialize the element-specific data
      // for the current element
      fe->reinit (elem);
 
      // Get the local to global degree of freedom maps
      computed_dof_map.dof_indices    (elem, dof_indices, var);
 
      // The number of quadrature points
      const unsigned int n_qp = qrule->n_points();
 
      // The number of shape functions
      const unsigned int n_sf =
        cast_int<unsigned int>(dof_indices.size());
 
      //
      // Begin the loop over the Quadrature points.
      //
      for (unsigned int qp=0; qp<n_qp; qp++)
        {
          // Real u_h = 0.;
          // RealGradient grad_u_h;
 
          typename FEGenericBase<OutputShape>::OutputNumber u_h(0.);
 
          typename FEGenericBase<OutputShape>::OutputNumberGradient grad_u_h;
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
          typename FEGenericBase<OutputShape>::OutputNumberTensor grad2_u_h;
#endif
          typename FEGenericBase<OutputShape>::OutputNumber curl_u_h(0.0);
          typename FEGenericBase<OutputShape>::OutputNumberDivergence div_u_h = 0.0;
 
          // Compute solution values at the current
          // quadrature point.  This requires a sum
          // over all the shape functions evaluated
          // at the quadrature point.
          for (unsigned int i=0; i<n_sf; i++)
            {
              // Values from current solution.
              u_h      += phi_values[i][qp]*computed_system.current_solution  (dof_indices[i]);
              grad_u_h += dphi_values[i][qp]*computed_system.current_solution (dof_indices[i]);
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
              grad2_u_h += d2phi_values[i][qp]*computed_system.current_solution (dof_indices[i]);
#endif
              if (FEInterface::field_type(fe_type) == TYPE_VECTOR)
                {
                  curl_u_h += (*curl_values)[i][qp]*computed_system.current_solution (dof_indices[i]);
                  div_u_h += (*div_values)[i][qp]*computed_system.current_solution (dof_indices[i]);
                }
            }
 
          // Compute the value of the error at this quadrature point
          typename FEGenericBase<OutputShape>::OutputNumber exact_val(0);
          RawAccessor<typename FEGenericBase<OutputShape>::OutputNumber> exact_val_accessor( exact_val, dim );
          if (_exact_values.size() > sys_num && _exact_values[sys_num])
            {
              for (unsigned int c = 0; c < n_vec_dim; c++)
                exact_val_accessor(c) =
                  _exact_values[sys_num]->
                  component(var_component+c, q_point[qp], time);
            }
          else if (_equation_systems_fine)
            {
              // FIXME: Needs to be updated for vector-valued elements
              DenseVector<Number> output(1);
              (*coarse_values)(q_point[qp],time,output,&subdomain_id);
              exact_val = output(0);
            }
          const typename FEGenericBase<OutputShape>::OutputNumber val_error = u_h - exact_val;
 
          // Add the squares of the error to each contribution
          Real error_sq = TensorTools::norm_sq(val_error);
          error_vals[0] += JxW[qp]*error_sq;
 
          Real norm = sqrt(error_sq);
          error_vals[3] += JxW[qp]*norm;
 
          if (error_vals[4]<norm) { error_vals[4] = norm; }
 
          // Compute the value of the error in the gradient at this
          // quadrature point
          typename FEGenericBase<OutputShape>::OutputNumberGradient exact_grad;
          RawAccessor<typename FEGenericBase<OutputShape>::OutputNumberGradient> exact_grad_accessor( exact_grad, LIBMESH_DIM );
          if (_exact_derivs.size() > sys_num && _exact_derivs[sys_num])
            {
              for (unsigned int c = 0; c < n_vec_dim; c++)
                for (unsigned int d = 0; d < LIBMESH_DIM; d++)
                  exact_grad_accessor(d + c*LIBMESH_DIM) =
                    _exact_derivs[sys_num]->
                    component(var_component+c, q_point[qp], time)(d);
            }
          else if (_equation_systems_fine)
            {
              // FIXME: Needs to be updated for vector-valued elements
              std::vector<Gradient> output(1);
              coarse_values->gradient(q_point[qp],time,output,&subdomain_id);
              exact_grad = output[0];
            }
 
          const typename FEGenericBase<OutputShape>::OutputNumberGradient grad_error = grad_u_h - exact_grad;
 
          error_vals[1] += JxW[qp]*grad_error.norm_sq();
 
 
          if (FEInterface::field_type(fe_type) == TYPE_VECTOR)
            {
              // Compute the value of the error in the curl at this
              // quadrature point
              typename FEGenericBase<OutputShape>::OutputNumber exact_curl(0.0);
              if (_exact_derivs.size() > sys_num && _exact_derivs[sys_num])
                {
                  exact_curl = TensorTools::curl_from_grad( exact_grad );
                }
              else if (_equation_systems_fine)
                {
                  // FIXME: Need to implement curl for MeshFunction and support reference
                  //        solution for vector-valued elements
                }
 
              const typename FEGenericBase<OutputShape>::OutputNumber curl_error = curl_u_h - exact_curl;
 
              error_vals[5] += JxW[qp]*TensorTools::norm_sq(curl_error);
 
              // Compute the value of the error in the divergence at this
              // quadrature point
              typename FEGenericBase<OutputShape>::OutputNumberDivergence exact_div = 0.0;
              if (_exact_derivs.size() > sys_num && _exact_derivs[sys_num])
                {
                  exact_div = TensorTools::div_from_grad( exact_grad );
                }
              else if (_equation_systems_fine)
                {
                  // FIXME: Need to implement div for MeshFunction and support reference
                  //        solution for vector-valued elements
                }
 
              const typename FEGenericBase<OutputShape>::OutputNumberDivergence div_error = div_u_h - exact_div;
 
              error_vals[6] += JxW[qp]*TensorTools::norm_sq(div_error);
            }
 
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
          // Compute the value of the error in the hessian at this
          // quadrature point
          typename FEGenericBase<OutputShape>::OutputNumberTensor exact_hess;
          RawAccessor<typename FEGenericBase<OutputShape>::OutputNumberTensor> exact_hess_accessor( exact_hess, dim );
          if (_exact_hessians.size() > sys_num && _exact_hessians[sys_num])
            {
              //FIXME: This needs to be implemented to support rank 3 tensors
              //       which can't happen until type_n_tensor is fully implemented
              //       and a RawAccessor<TypeNTensor> is fully implemented
              if (FEInterface::field_type(fe_type) == TYPE_VECTOR)
                libmesh_not_implemented();
 
              for (unsigned int c = 0; c < n_vec_dim; c++)
                for (unsigned int d = 0; d < dim; d++)
                  for (unsigned int e =0; e < dim; e++)
                    exact_hess_accessor(d + e*dim + c*dim*dim) =
                      _exact_hessians[sys_num]->
                      component(var_component+c, q_point[qp], time)(d,e);
            }
          else if (_equation_systems_fine)
            {
              // FIXME: Needs to be updated for vector-valued elements
              std::vector<Tensor> output(1);
              coarse_values->hessian(q_point[qp],time,output,&subdomain_id);
              exact_hess = output[0];
            }
 
          const typename FEGenericBase<OutputShape>::OutputNumberTensor grad2_error = grad2_u_h - exact_hess;
 
          // FIXME: PB: Is this what we want for rank 3 tensors?
          error_vals[2] += JxW[qp]*grad2_error.norm_sq();
#endif
 
        } // end qp loop
    } // end element loop
 
  // Add up the error values on all processors, except for the L-infty
  // norm, for which the maximum is computed.
  Real l_infty_norm = error_vals[4];
  communicator.max(l_infty_norm);
  communicator.sum(error_vals);
  error_vals[4] = l_infty_norm;
}

References _equation_systems, _equation_systems_fine, _exact_derivs, _exact_hessians, _exact_values, _excluded_subdomains, _extra_order, libMesh::Variable::active_on_subdomain(), libMesh::FEGenericBase< OutputType >::build(), libMesh::NumericVector< T >::build(), libMesh::ParallelObject::comm(), libMesh::TensorTools::curl_from_grad(), libMesh::System::current_solution(), libMesh::FEType::default_quadrature_rule(), dim, libMesh::TensorTools::div_from_grad(), libMesh::DofMap::dof_indices(), libMesh::err, exact_grad(), libMesh::FEInterface::field_type(), libMesh::ReferenceElem::get(), libMesh::FEGenericBase< OutputType >::get_curl_phi(), libMesh::FEGenericBase< OutputType >::get_d2phi(), libMesh::FEGenericBase< OutputType >::get_div_phi(), libMesh::System::get_dof_map(), libMesh::FEGenericBase< OutputType >::get_dphi(), libMesh::FEAbstract::get_JxW(), libMesh::System::get_mesh(), libMesh::FEGenericBase< OutputType >::get_phi(), libMesh::EquationSystems::get_system(), libMesh::FEAbstract::get_xyz(), libMesh::libmesh_assert(), mesh, libMesh::QBase::n_points(), libMesh::FEInterface::n_vec_dim(), std::norm(), libMesh::TensorTools::norm_sq(), libMesh::TypeVector< T >::norm_sq(), libMesh::System::number(), libMesh::Real, libMesh::FEAbstract::reinit(), libMesh::SERIAL, libMesh::System::solution, std::sqrt(), libMesh::TYPE_VECTOR, libMesh::System::update_global_solution(), libMesh::System::variable(), libMesh::System::variable_number(), libMesh::System::variable_scalar_number(), and libMesh::DofMap::variable_type().

attach_exact_deriv() [1/2]

void libMesh::ExactSolution::attach_exact_deriv

(

GradientFunctionPointer

gptr

)

Definition at line 134 of file exact_solution.C.

{
  libmesh_assert(gptr);
 
  // Clear out any previous _exact_derivs entries, then add a new
  // entry for each system.
  _exact_derivs.clear();
 
  for (auto sys : IntRange<unsigned int>(0, _equation_systems.n_systems()))
    {
      const System & system = _equation_systems.get_system(sys);
      _exact_derivs.emplace_back(libmesh_make_unique<WrappedFunction<Gradient>>(system, gptr, &_equation_systems.parameters));
    }
 
  // If we're using exact values, we're not using a fine grid solution
  _equation_systems_fine = nullptr;
}

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

attach_exact_deriv() [2/2]

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 164 of file exact_solution.C.

{
  if (_exact_derivs.size() <= sys_num)
    _exact_derivs.resize(sys_num+1);
 
  if (g)
    _exact_derivs[sys_num] = g->clone();
}

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

Referenced by main().

attach_exact_derivs()

void libMesh::ExactSolution::attach_exact_derivs

(

const 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 153 of file exact_solution.C.

{
  // Automatically delete any previous _exact_derivs entries, then add a new
  // entry for each system.
  _exact_derivs.clear();
 
  for (auto ptr : g)
    _exact_derivs.emplace_back(ptr ? ptr->clone() : nullptr);
}

References _exact_derivs.

Referenced by main().

attach_exact_hessian() [1/2]

void libMesh::ExactSolution::attach_exact_hessian

(

HessianFunctionPointer

hptr

)

Definition at line 175 of file exact_solution.C.

{
  libmesh_assert(hptr);
 
  // Clear out any previous _exact_hessians entries, then add a new
  // entry for each system.
  _exact_hessians.clear();
 
  for (auto sys : IntRange<unsigned int>(0, _equation_systems.n_systems()))
    {
      const System & system = _equation_systems.get_system(sys);
      _exact_hessians.emplace_back(libmesh_make_unique<WrappedFunction<Tensor>>(system, hptr, &_equation_systems.parameters));
    }
 
  // If we're using exact values, we're not using a fine grid solution
  _equation_systems_fine = nullptr;
}

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

attach_exact_hessian() [2/2]

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 205 of file exact_solution.C.

{
  if (_exact_hessians.size() <= sys_num)
    _exact_hessians.resize(sys_num+1);
 
  if (h)
    _exact_hessians[sys_num] = h->clone();
}

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

Referenced by main().

attach_exact_hessians()

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 194 of file exact_solution.C.

{
  // Automatically delete any previous _exact_hessians entries, then add a new
  // entry for each system.
  _exact_hessians.clear();
 
  for (auto ptr : h)
    _exact_hessians.emplace_back(ptr ? ptr->clone() : nullptr);
}

References _exact_hessians.

attach_exact_value() [1/2]

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 123 of file exact_solution.C.

{
  if (_exact_values.size() <= sys_num)
    _exact_values.resize(sys_num+1);
 
  if (f)
    _exact_values[sys_num] = f->clone();
}

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

Referenced by main().

attach_exact_value() [2/2]

void libMesh::ExactSolution::attach_exact_value

(

ValueFunctionPointer

fptr

)

Definition at line 93 of file exact_solution.C.

{
  libmesh_assert(fptr);
 
  // Clear out any previous _exact_values entries, then add a new
  // entry for each system.
  _exact_values.clear();
 
  for (auto sys : IntRange<unsigned int>(0, _equation_systems.n_systems()))
    {
      const System & system = _equation_systems.get_system(sys);
      _exact_values.emplace_back(libmesh_make_unique<WrappedFunction<Number>>(system, fptr, &_equation_systems.parameters));
    }
 
  // If we're using exact values, we're not using a fine grid solution
  _equation_systems_fine = nullptr;
}

References _equation_systems, _equation_systems_fine, _exact_values, fptr(), libMesh::EquationSystems::get_system(), libMesh::libmesh_assert(), libMesh::EquationSystems::n_systems(), and libMesh::EquationSystems::parameters.

attach_exact_values()

void libMesh::ExactSolution::attach_exact_values

(

const 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 112 of file exact_solution.C.

{
  // Automatically delete any previous _exact_values entries, then add a new
  // entry for each system.
  _exact_values.clear();
 
  for (auto ptr : f)
    _exact_values.emplace_back(ptr ? ptr->clone() : nullptr);
}

References _exact_values.

Referenced by main().

attach_reference_solution()

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 81 of file exact_solution.C.

{
  libmesh_assert(es_fine);
  _equation_systems_fine = es_fine;
 
  // If we're using a fine grid solution, we're not using exact values
  _exact_values.clear();
  _exact_derivs.clear();
  _exact_hessians.clear();
}

References _equation_systems_fine, _exact_derivs, _exact_hessians, _exact_values, and libMesh::libmesh_assert().

Referenced by main().

compute_error()

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 227 of file exact_solution.C.

{
  // Check the inputs for validity, and get a reference
  // to the proper location to store the error
  std::vector<Real> & error_vals = this->_check_inputs(sys_name,
                                                       unknown_name);
 
  libmesh_assert( _equation_systems.has_system(sys_name) );
  const System & sys = _equation_systems.get_system<System>( sys_name );
 
  libmesh_assert( sys.has_variable( unknown_name ) );
  switch( FEInterface::field_type(sys.variable_type( unknown_name )) )
    {
    case TYPE_SCALAR:
      {
        this->_compute_error<Real>(sys_name,
                                   unknown_name,
                                   error_vals);
        break;
      }
    case TYPE_VECTOR:
      {
        this->_compute_error<RealGradient>(sys_name,
                                           unknown_name,
                                           error_vals);
        break;
      }
    default:
      libmesh_error_msg("Invalid variable type!");
    }
}

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

Referenced by main().

error_norm()

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

Returns

The error in the requested norm for the system sys_name for the unknown unknown_name.

Note

No error computations are actually performed, you must call compute_error() for that.

The result is not exact, but an approximation based on the chosen quadrature rule.

Definition at line 264 of file exact_solution.C.

{
  // Check the inputs for validity, and get a reference
  // to the proper location to store the error
  std::vector<Real> & error_vals = this->_check_inputs(sys_name,
                                                       unknown_name);
 
  libmesh_assert(_equation_systems.has_system(sys_name));
  libmesh_assert(_equation_systems.get_system(sys_name).has_variable( unknown_name ));
  const FEType & fe_type = _equation_systems.get_system(sys_name).variable_type(unknown_name);
 
  switch (norm)
    {
    case L2:
      return std::sqrt(error_vals[0]);
    case H1:
      return std::sqrt(error_vals[0] + error_vals[1]);
    case H2:
      return std::sqrt(error_vals[0] + error_vals[1] + error_vals[2]);
    case HCURL:
      {
        if (FEInterface::field_type(fe_type) == TYPE_SCALAR)
          libmesh_error_msg("Cannot compute HCurl error norm of scalar-valued variables!");
        else
          return std::sqrt(error_vals[0] + error_vals[5]);
      }
    case HDIV:
      {
        if (FEInterface::field_type(fe_type) == TYPE_SCALAR)
          libmesh_error_msg("Cannot compute HDiv error norm of scalar-valued variables!");
        else
          return std::sqrt(error_vals[0] + error_vals[6]);
      }
    case H1_SEMINORM:
      return std::sqrt(error_vals[1]);
    case H2_SEMINORM:
      return std::sqrt(error_vals[2]);
    case HCURL_SEMINORM:
      {
        if (FEInterface::field_type(fe_type) == TYPE_SCALAR)
          libmesh_error_msg("Cannot compute HCurl error seminorm of scalar-valued variables!");
        else
          return std::sqrt(error_vals[5]);
      }
    case HDIV_SEMINORM:
      {
        if (FEInterface::field_type(fe_type) == TYPE_SCALAR)
          libmesh_error_msg("Cannot compute HDiv error seminorm of scalar-valued variables!");
        else
          return std::sqrt(error_vals[6]);
      }
    case L1:
      return error_vals[3];
    case L_INF:
      return error_vals[4];
 
    default:
      libmesh_error_msg("Currently only Sobolev norms/seminorms are supported!");
    }
}

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

Referenced by hcurl_error(), hdiv_error(), and main().

extra_quadrature_order()

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 190 of file exact_solution.h.

  { _extra_order = extraorder; }

References _extra_order.

Referenced by main().

h1_error()

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

Returns

The H1 error for the system sys_name for the unknown unknown_name.

Note

No error computations are actually performed, you must call compute_error() for that.

Definition at line 393 of file exact_solution.C.

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

References _check_inputs(), and std::sqrt().

Referenced by main().

h2_error()

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

Returns

The H2 error for the system sys_name for the unknown unknown_name.

Note

No error computations are actually performed, you must call compute_error() for that.

Definition at line 425 of file exact_solution.C.

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

References _check_inputs(), and std::sqrt().

Referenced by main().

hcurl_error()

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

Returns

The H(curl) error for the system sys_name for the unknown unknown_name.

Note

No error computations are actually performed, you must call compute_error() for that.

This is only valid for vector-valued elements. An error is thrown if requested for scalar-valued elements.

Definition at line 410 of file exact_solution.C.

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

References error_norm(), and libMesh::HCURL.

Referenced by main().

hdiv_error()

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

Returns

The H(div) error for the system sys_name for the unknown unknown_name.

Note

No error computations are actually performed, you must call compute_error() for that.

This is only valid for vector-valued elements. An error is thrown if requested for scalar-valued elements.

Definition at line 417 of file exact_solution.C.

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

References error_norm(), and libMesh::HDIV.

l1_error()

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

Returns

The integrated L1 error for the system sys_name for the unknown unknown_name.

Note

No error computations are actually performed, you must call compute_error() for that.

Definition at line 353 of file exact_solution.C.

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

References _check_inputs().

l2_error()

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

Returns

The integrated L2 error for the system sys_name for the unknown unknown_name.

Note

No error computations are actually performed, you must call compute_error() for that.

Definition at line 333 of file exact_solution.C.

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

References _check_inputs(), and std::sqrt().

Referenced by main().

l_inf_error()

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

Returns

The L_INF error for the system sys_name for the unknown unknown_name.

Note

No error computations are actually performed, you must call compute_error() for 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 373 of file exact_solution.C.

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

References _check_inputs().

operator=() [1/2]

ExactSolution& libMesh::ExactSolution::operator= ( const ExactSolution &  )

delete

operator=() [2/2]

ExactSolution& libMesh::ExactSolution::operator= ( ExactSolution &&  )

delete

set_excluded_subdomains()

void libMesh::ExactSolution::set_excluded_subdomains

(

const std::set< subdomain_id_type > &

excluded

)

The user can indicate that elements in certain subdomains should be excluded from the error calculation by passing in a set of subdomain ids to ignore.

By default, all subdomains are considered in the error calculation.

Definition at line 76 of file exact_solution.C.

{
  _excluded_subdomains = excluded;
}

References _excluded_subdomains.

Member Data Documentation

_equation_systems

const EquationSystems& libMesh::ExactSolution::_equation_systems

private

Constant reference to the EquationSystems object used for the simulation.

Definition at line 359 of file exact_solution.h.

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

_equation_systems_fine

const EquationSystems* libMesh::ExactSolution::_equation_systems_fine

private

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

Definition at line 365 of file exact_solution.h.

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

_errors

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 353 of file exact_solution.h.

Referenced by _check_inputs(), and ExactSolution().

_exact_derivs

std::vector<std::unique_ptr<FunctionBase<Gradient> > > libMesh::ExactSolution::_exact_derivs

private

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

Definition at line 329 of file exact_solution.h.

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

_exact_hessians

std::vector<std::unique_ptr<FunctionBase<Tensor> > > libMesh::ExactSolution::_exact_hessians

private

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

Definition at line 335 of file exact_solution.h.

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

_exact_values

std::vector<std::unique_ptr<FunctionBase<Number> > > libMesh::ExactSolution::_exact_values

private

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

Definition at line 323 of file exact_solution.h.

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

_excluded_subdomains

std::set<subdomain_id_type> libMesh::ExactSolution::_excluded_subdomains

private

Elements in a subdomain from this set are skipped during the error computation.

Definition at line 376 of file exact_solution.h.

Referenced by _compute_error(), and set_excluded_subdomains().

_extra_order

int libMesh::ExactSolution::_extra_order

private

Extra order to use for quadrature rule.

Definition at line 370 of file exact_solution.h.

Referenced by _compute_error(), and extra_quadrature_order().

---

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

• include/error_estimation/exact_solution.h
• src/error_estimation/exact_solution.C