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 (std::string_view sys_name, std::string_view 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 (std::string_view sys_name, std::string_view unknown_name)
Real	l1_error (std::string_view sys_name, std::string_view unknown_name)
Real	l_inf_error (std::string_view sys_name, std::string_view unknown_name)
Real	h1_error (std::string_view sys_name, std::string_view unknown_name)
Real	hcurl_error (std::string_view sys_name, std::string_view unknown_name)
Real	hdiv_error (std::string_view sys_name, std::string_view unknown_name)
Real	h2_error (std::string_view sys_name, std::string_view unknown_name)
Real	error_norm (std::string_view sys_name, std::string_view unknown_name, const FEMNormType &norm)

Private Types
typedef std::map< std::string, std::vector< Real >, std::less<> >	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 (std::string_view sys_name, std::string_view 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 (std::string_view sys_name, std::string_view 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, std::less<> >	_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 65 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 149 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 172 of file exact_solution.h.

SystemErrorMap

typedef std::map<std::string, std::vector<Real>, std::less<> > 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 341 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 126 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 45 of file exact_solution.C.

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

                                                       :
  _equation_systems(es),
  _equation_systems_fine(nullptr),
  _extra_order(1)
{
  // Initialize the _errors data structure which holds all
  // the eventual values of the error.
  for (auto sys : make_range(_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 : make_range(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;
    }
}

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 ( std::string_view  sys_name, std::string_view  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 220 of file exact_solution.C.

References _errors.

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

{
  // 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);
}

_compute_error()

template<typename OutputShape >

template LIBMESH_EXPORT void libMesh::ExactSolution::_compute_error< RealGradient > ( std::string_view  sys_name, std::string_view  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 452 of file exact_solution.C.

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(), communicator, 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::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::EquationSystems::get_mesh(), libMesh::System::get_mesh(), libMesh::FEGenericBase< OutputType >::get_phi(), libMesh::EquationSystems::get_system(), libMesh::FEAbstract::get_xyz(), libMesh::libmesh_assert(), libMesh::make_range(), mesh, libMesh::QBase::n_points(), libMesh::FEInterface::n_vec_dim(), libMesh::TensorTools::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::DofMap::variable_type(), and libMesh::System::variable_type().

{
  // 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 MeshBase & mesh = computed_system.get_mesh();

  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);
  unsigned int var_component = 0;
  for (const auto var_num : make_range(var))
  {
    const auto & var_fe_type = computed_system.variable_type(var_num);
    const auto var_vec_dim = FEInterface::n_vec_dim(mesh, var_fe_type);
    var_component += var_vec_dim;
  }

  // Prepare a global solution, a serialized mesh, and a MeshFunction
  // of the coarse system, if we need them for fine-system integration
  std::unique_ptr<MeshFunction> coarse_values;
  std::unique_ptr<NumericVector<Number>> comparison_soln =
    NumericVector<Number>::build(_equation_systems.comm());
  MeshSerializer
    serial(const_cast<MeshBase&>(_equation_systems.get_mesh()),
           _equation_systems_fine);
  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 = std::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();

  // 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);
  const auto field_type = FEInterface::field_type(fe_type);

  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 (field_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
      // Not computed for vector-valued elements
      const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputTensor>> *
        d2phi_values = nullptr;

      if (field_type != TYPE_VECTOR)
        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]);
              if (field_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]);
                }
              else
                {
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
                  grad2_u_h += (*d2phi_values)[i][qp]*computed_system.current_solution (dof_indices[i]);
#endif
                }
            }

          // 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 (field_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 (field_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);

              // FIXME: operator- is not currently implemented for TypeNTensor
              const typename FEGenericBase<OutputShape>::OutputNumberTensor grad2_error = grad2_u_h - exact_hess;
              error_vals[2] += JxW[qp]*grad2_error.norm_sq();
            }
          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];

              // FIXME: operator- is not currently implemented for TypeNTensor
              const typename FEGenericBase<OutputShape>::OutputNumberTensor grad2_error = grad2_u_h - exact_hess;
              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;
}

attach_exact_deriv() [1/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 168 of file exact_solution.C.

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

Referenced by main().

{
  if (_exact_derivs.size() <= sys_num)
    _exact_derivs.resize(sys_num+1);

  if (g)
    _exact_derivs[sys_num] = g->clone();
}

attach_exact_deriv() [2/2]

void libMesh::ExactSolution::attach_exact_deriv

(

GradientFunctionPointer

gptr

)

Definition at line 138 of file exact_solution.C.

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

{
  libmesh_assert(gptr);

  // Clear out any previous _exact_derivs entries, then add a new
  // entry for each system.
  _exact_derivs.clear();

  for (auto sys : make_range(_equation_systems.n_systems()))
    {
      const System & system = _equation_systems.get_system(sys);
      _exact_derivs.emplace_back(std::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;
}

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

References _exact_derivs.

Referenced by main().

{
  // 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);
}

attach_exact_hessian() [1/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 209 of file exact_solution.C.

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

Referenced by main().

{
  if (_exact_hessians.size() <= sys_num)
    _exact_hessians.resize(sys_num+1);

  if (h)
    _exact_hessians[sys_num] = h->clone();
}

attach_exact_hessian() [2/2]

void libMesh::ExactSolution::attach_exact_hessian

(

HessianFunctionPointer

hptr

)

Definition at line 179 of file exact_solution.C.

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

{
  libmesh_assert(hptr);

  // Clear out any previous _exact_hessians entries, then add a new
  // entry for each system.
  _exact_hessians.clear();

  for (auto sys : make_range(_equation_systems.n_systems()))
    {
      const System & system = _equation_systems.get_system(sys);
      _exact_hessians.emplace_back(std::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;
}

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

References _exact_hessians.

{
  // 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);
}

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

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

Referenced by main(), and ConstraintOperatorTest::test1DCoarseningOperator().

{
  if (_exact_values.size() <= sys_num)
    _exact_values.resize(sys_num+1);

  if (f)
    _exact_values[sys_num] = f->clone();
}

attach_exact_value() [2/2]

void libMesh::ExactSolution::attach_exact_value

(

ValueFunctionPointer

fptr

)

Definition at line 97 of file exact_solution.C.

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

{
  libmesh_assert(fptr);

  // Clear out any previous _exact_values entries, then add a new
  // entry for each system.
  _exact_values.clear();

  for (auto sys : make_range(_equation_systems.n_systems()))
    {
      const System & system = _equation_systems.get_system(sys);
      _exact_values.emplace_back(std::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;
}

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

References _exact_values.

Referenced by main().

{
  // 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);
}

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

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

Referenced by main().

{
  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();
}

compute_error()

void libMesh::ExactSolution::compute_error	(	std::string_view	sys_name,
		std::string_view	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 231 of file exact_solution.C.

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(), and ConstraintOperatorTest::test1DCoarseningOperator().

{
  // 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!");
    }
}

error_norm()

Real libMesh::ExactSolution::error_norm	(	std::string_view	sys_name,
		std::string_view	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 268 of file exact_solution.C.

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(), libMesh::TensorTools::norm(), std::sqrt(), and libMesh::TYPE_SCALAR.

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

{
  // 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:
      {
        libmesh_error_msg_if(FEInterface::field_type(fe_type) == TYPE_SCALAR,
                             "Cannot compute HCurl error norm of scalar-valued variables!");

        return std::sqrt(error_vals[0] + error_vals[5]);
      }
    case HDIV:
      {
        libmesh_error_msg_if(FEInterface::field_type(fe_type) == TYPE_SCALAR,
                             "Cannot compute HDiv error norm of scalar-valued variables!");

        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:
      {
        libmesh_error_msg_if(FEInterface::field_type(fe_type) == TYPE_SCALAR,
                             "Cannot compute HCurl error seminorm of scalar-valued variables!");

        return std::sqrt(error_vals[5]);
      }
    case HDIV_SEMINORM:
      {
        libmesh_error_msg_if(FEInterface::field_type(fe_type) == TYPE_SCALAR,
                             "Cannot compute HDiv error seminorm of scalar-valued variables!");

        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!");
    }
}

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.

The default extraorder is 1, because properly integrating L2 error requires integrating the squares of terms with order p+1, and 2p+2 is 1 higher than what we default to using for reasonable mass matrix integration.

Definition at line 185 of file exact_solution.h.

References _extra_order.

Referenced by main().

  { _extra_order = extraorder; }

h1_error()

Real libMesh::ExactSolution::h1_error	(	std::string_view	sys_name,
		std::string_view	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 397 of file exact_solution.C.

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

Referenced by main().

{
  // 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]);
}

h2_error()

Real libMesh::ExactSolution::h2_error	(	std::string_view	sys_name,
		std::string_view	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 429 of file exact_solution.C.

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

Referenced by main().

{
  // 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]);
}

hcurl_error()

Real libMesh::ExactSolution::hcurl_error	(	std::string_view	sys_name,
		std::string_view	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 414 of file exact_solution.C.

References error_norm(), and libMesh::HCURL.

Referenced by main().

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

hdiv_error()

Real libMesh::ExactSolution::hdiv_error	(	std::string_view	sys_name,
		std::string_view	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 421 of file exact_solution.C.

References error_norm(), and libMesh::HDIV.

Referenced by main().

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

l1_error()

Real libMesh::ExactSolution::l1_error	(	std::string_view	sys_name,
		std::string_view	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 357 of file exact_solution.C.

References _check_inputs().

{

  // 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];
}

l2_error()

Real libMesh::ExactSolution::l2_error	(	std::string_view	sys_name,
		std::string_view	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 337 of file exact_solution.C.

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

Referenced by main(), and ConstraintOperatorTest::test1DCoarseningOperator().

{

  // 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]);
}

l_inf_error()

Real libMesh::ExactSolution::l_inf_error	(	std::string_view	sys_name,
		std::string_view	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 377 of file exact_solution.C.

References _check_inputs().

{

  // 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];
}

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

References _excluded_subdomains.

{
  _excluded_subdomains = excluded;
}

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 355 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 361 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, std::less<> > 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 349 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 324 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 330 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 318 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 372 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 366 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