ExactSolution Class Reference

#include <exact_solution.h>

List of all members.

Public Member Functions
	ExactSolution (const EquationSystems &es)
	~ExactSolution ()
void	attach_reference_solution (const EquationSystems *es_fine)
void	attach_exact_value (Number fptr(const Point &p, const Parameters &Parameters, const std::string &sys_name, const std::string &unknown_name))
void	attach_exact_deriv (Gradient fptr(const Point &p, const Parameters &parameters, const std::string &sys_name, const std::string &unknown_name))
void	attach_exact_hessian (Tensor fptr(const Point &p, const Parameters &parameters, const std::string &sys_name, const std::string &unknown_name))
void	extra_quadrature_order (const int extraorder)
void	compute_error (const std::string &sys_name, const std::string &unknown_name)
Real	l2_error (const std::string &sys_name, const std::string &unknown_name)
Real	l1_error (const std::string &sys_name, const std::string &unknown_name)
Real	l_inf_error (const std::string &sys_name, const std::string &unknown_name)
Real	h1_error (const std::string &sys_name, const std::string &unknown_name)
Real	h2_error (const std::string &sys_name, const std::string &unknown_name)
Real	error_norm (const std::string &sys_name, const std::string &unknown_name, const FEMNormType &norm)
Private Types
typedef std::map< std::string, std::vector< Real > >	SystemErrorMap
Private Member Functions
void	_compute_error (const std::string &sys_name, const std::string &unknown_name, std::vector< Real > &error_vals)
std::vector< Real > &	_check_inputs (const std::string &sys_name, const std::string &unknown_name)
Private Attributes
Number(*	_exact_value )(const Point &p, const Parameters &parameters, const std::string &sys_name, const std::string &unknown_name)
Gradient(*	_exact_deriv )(const Point &p, const Parameters &parameters, const std::string &sys_name, const std::string &unknown_name)
Tensor(*	_exact_hessian )(const Point &p, const Parameters &parameters, const std::string &sys_name, const std::string &unknown_name)
std::map< std::string, SystemErrorMap >	_errors
const EquationSystems &	_equation_systems
const EquationSystems *	_equation_systems_fine
int	_extra_order

---

Detailed Description

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

Author:

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

Definition at line 62 of file exact_solution.h.

---

Member Typedef Documentation

typedef std::map<std::string, std::vector<Real> > 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 245 of file exact_solution.h.

---

Constructor & Destructor Documentation

ExactSolution::ExactSolution

(

const EquationSystems &

es

)

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

Definition at line 38 of file exact_solution.C.

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

                                                      :
  _exact_value (NULL),
  _exact_deriv (NULL),
  _exact_hessian (NULL),
  _equation_systems(es),
  _equation_systems_fine(NULL),
  _extra_order(0)
{
  // Initialize the _errors data structure which holds all
  // the eventual values of the error.
  for (unsigned int sys=0; sys<_equation_systems.n_systems(); ++sys)
    {
      // 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 (unsigned int var=0; var<system.n_vars(); ++var)
        {
          // 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::~ExactSolution

(

)

[inline]

Destructor.

Definition at line 76 of file exact_solution.h.

{}

---

Member Function Documentation

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

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

Definition at line 124 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().

{
  // If no exact solution function or fine grid solution has been
  // attached, we now just compute the solution norm (i.e. the
  // difference from an "exact solution" of zero
/*
  if (_exact_value == NULL)
    {
      std::cerr << "Cannot compute error, you must provide a "
                << "function which computes the exact solution."
                << std::endl;
      libmesh_error();
    }
*/
  
  // Make sure the requested sys_name exists.
  std::map<std::string, SystemErrorMap>::iterator sys_iter =
    _errors.find(sys_name);
  
  if (sys_iter == _errors.end())
    {
      std::cerr << "Sorry, couldn't find the requested system '"
                << sys_name << "'."
                << std::endl;
      libmesh_error();
    }
  
  // Make sure the requested unknown_name exists.
  SystemErrorMap::iterator var_iter = (*sys_iter).second.find(unknown_name);

  if (var_iter == (*sys_iter).second.end())
    {
      std::cerr << "Sorry, couldn't find the requested variable '"
                << unknown_name << "'."
                << std::endl;
      libmesh_error();
    }

  // Return a reference to the proper error entry
  return (*var_iter).second;
}

void ExactSolution::_compute_error	(	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 349 of file exact_solution.C.

References _equation_systems, _equation_systems_fine, _exact_deriv, _exact_hessian, _exact_value, _extra_order, _mesh, MeshBase::active_local_elements_begin(), MeshBase::active_local_elements_end(), FEBase::build(), System::current_solution(), FEType::default_quadrature_rule(), DofMap::dof_indices(), AutoPtr< Tp >::get(), System::get_dof_map(), System::get_mesh(), EquationSystems::get_system(), libmesh_norm(), std::max(), MeshBase::mesh_dimension(), EquationSystems::parameters, libMeshEnums::SERIAL, TypeTensor< T >::size_sq(), TypeVector< T >::size_sq(), System::solution, System::update_global_solution(), System::variable_number(), and DofMap::variable_type().

Referenced by compute_error().

{
  // This function must be run on all processors at once
  parallel_only();

  // Make sure we aren't "overconfigured"
  libmesh_assert (!(_exact_value && _equation_systems_fine));

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

  // Prepare a global solution and a MeshFunction of the coarse system if we need one
  AutoPtr<MeshFunction> coarse_values;
  AutoPtr<NumericVector<Number> > comparison_soln = NumericVector<Number>::build();
  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 = AutoPtr<MeshFunction>
        (new MeshFunction(_equation_systems,
                          *comparison_soln,
                          comparison_system.get_dof_map(),
                          comparison_system.variable_number(unknown_name)));
      coarse_values->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();

  // Zero the error before summation
  error_vals = std::vector<Real>(5, 0.);

  // get the EquationSystems parameters for use with _exact_value
  const Parameters& parameters = this->_equation_systems.parameters;

  // Get the current time, in case the exact solution depends on it.
  // Steady systems of equations do not have a time parameter, so this
  // routine needs to take that into account.
  // FIXME!!!
  // const Real time = 0.;//_equation_systems.parameter("time");  

  // Construct Quadrature rule based on default quadrature order
  const unsigned int var = computed_system.variable_number(unknown_name);
  const FEType& fe_type  = computed_dof_map.variable_type(var);

  AutoPtr<QBase> qrule =
    fe_type.default_quadrature_rule (_mesh.mesh_dimension(),
                                     _extra_order);

  // Construct finite element object
  
  AutoPtr<FEBase> fe(FEBase::build(_mesh.mesh_dimension(), fe_type));

  // Attach quadrature rule to FE object
  fe->attach_quadrature_rule (qrule.get());
  
  // 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<Real> >&  phi_values         = fe->get_phi();
  
  // The value of the shape function gradients at the quadrature points
  const std::vector<std::vector<RealGradient> >& dphi_values = fe->get_dphi();

#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
  // The value of the shape function second derivatives at the quadrature points
  const std::vector<std::vector<RealTensor> >& 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();
            
  // The global degree of freedom indices associated
  // with the local degrees of freedom.
  std::vector<unsigned int> dof_indices;


  //
  // Begin the loop over the elements
  //
  MeshBase::const_element_iterator       el     = _mesh.active_local_elements_begin();
  const MeshBase::const_element_iterator end_el = _mesh.active_local_elements_end(); 

  for ( ; el != end_el; ++el)
    {
      // Store a pointer to the element we are currently
      // working on.  This allows for nicer syntax later.
      const Elem* elem = *el;

      // 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 = 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;

          Number u_h = 0.;

          Gradient grad_u_h;
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
          Tensor grad2_u_h;
#endif
  

          // Compute solution values at the current
          // quadrature point.  This reqiures 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
            }

          // Compute the value of the error at this quadrature point
          Number exact_val = 0.0;
          if (_exact_value)
            exact_val = _exact_value(q_point[qp], parameters,
                                     sys_name, unknown_name);
          else if (_equation_systems_fine)
            exact_val = (*coarse_values)(q_point[qp]);

          const Number val_error = u_h - exact_val;

          // Add the squares of the error to each contribution
          error_vals[0] += JxW[qp]*libmesh_norm(val_error);
          Real norm = sqrt(libmesh_norm(val_error));
          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
          Gradient exact_grad;
          if (_exact_deriv)
            exact_grad = _exact_deriv(q_point[qp], parameters,
                                      sys_name, unknown_name);
          else if (_equation_systems_fine)
            exact_grad = coarse_values->gradient(q_point[qp]);

          const Gradient grad_error = grad_u_h - exact_grad;

          error_vals[1] += JxW[qp]*grad_error.size_sq();


#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
          // Compute the value of the error in the hessian at this
          // quadrature point
          Tensor exact_hess;
          if (_exact_hessian)
            exact_hess = _exact_hessian(q_point[qp], parameters,
                                        sys_name, unknown_name);
          else if (_equation_systems_fine)
            exact_hess = coarse_values->hessian(q_point[qp]);

          const Tensor grad2_error = grad2_u_h - exact_hess;

          error_vals[2] += JxW[qp]*grad2_error.size_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];
  Parallel::max(l_infty_norm);
  Parallel::sum(error_vals);
  error_vals[4] = l_infty_norm;
}

void ExactSolution::attach_exact_deriv

(

Gradient

fptrconst Point &p,const Parameters &parameters,const std::string &sys_name,const std::string &unknown_name

)

Attach function similar to system.h which allows the user to attach an arbitrary function which computes the exact derivative of the solution at any point.

Definition at line 96 of file exact_solution.C.

References _equation_systems_fine, and _exact_deriv.

{
  libmesh_assert (fptr != NULL);
  _exact_deriv = fptr;

  // If we're using exact values, we're not using a fine grid solution
  _equation_systems_fine = NULL;
}

void ExactSolution::attach_exact_hessian

(

Tensor

fptrconst Point &p,const Parameters &parameters,const std::string &sys_name,const std::string &unknown_name

)

Attach function similar to system.h which allows the user to attach an arbitrary function which computes the exact second derivatives of the solution at any point.

Definition at line 109 of file exact_solution.C.

References _equation_systems_fine, and _exact_hessian.

{
  libmesh_assert (fptr != NULL);
  _exact_hessian = fptr;

  // If we're using exact values, we're not using a fine grid solution
  _equation_systems_fine = NULL;
}

void ExactSolution::attach_exact_value

(

Number

fptrconst Point &p, const Parameters &Parameters, const std::string &sys_name, const std::string &unknown_name

)

Attach function similar to system.h which allows the user to attach an arbitrary function which computes the exact value of the solution at any point.

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

References _equation_systems_fine, _exact_deriv, _exact_hessian, and _exact_value.

{
  libmesh_assert (es_fine != NULL);
  _equation_systems_fine = es_fine;

  // If we're using a fine grid solution, we're not using exact values
  _exact_value = NULL;
  _exact_deriv = NULL;
  _exact_hessian = NULL;
}

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

References _check_inputs(), and _compute_error().

{
  // 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);
  this->_compute_error(sys_name,
                       unknown_name,
                       error_vals);
}

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

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

Definition at line 188 of file exact_solution.C.

References _check_inputs(), libMeshEnums::H1, libMeshEnums::H1_SEMINORM, libMeshEnums::H2, libMeshEnums::H2_SEMINORM, libMeshEnums::L1, libMeshEnums::L2, and libMeshEnums::L_INF.

{
  // 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);
  
  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 H1_SEMINORM:
      return std::sqrt(error_vals[1]);
    case H2_SEMINORM:
      return std::sqrt(error_vals[2]);
    case L1:
      return error_vals[3];
    case L_INF:
      return error_vals[4];
    // Currently only Sobolev norms/seminorms are supported
    default:
      libmesh_error();
    }
}

void ExactSolution::extra_quadrature_order

(

const int

extraorder

)

[inline]

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

Definition at line 123 of file exact_solution.h.

References _extra_order.

    { _extra_order = extraorder; }

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

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

Definition at line 285 of file exact_solution.C.

References _check_inputs().

{
  // 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.
/*
  if (_exact_deriv == NULL)
    {
      std::cerr << "Cannot compute H1 error, you must provide a "
                << "function which computes the gradient of the exact solution."
                << std::endl;
      libmesh_error();
    }
*/
  
  // 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]);
}

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

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

Definition at line 316 of file exact_solution.C.

References _check_inputs().

{
  // If the user has supplied no exact derivative functions, we
  // just integrate the H1 norm of the solution; i.e. its
  // difference from an "exact solution" of zero.
/*
  if (_exact_deriv == NULL || _exact_hessian == NULL)
    {
      std::cerr << "Cannot compute H2 error, you must provide functions "
                << "which computes the gradient and hessian of the "
                << "exact solution."
                << std::endl;
      libmesh_error();
    }
*/
  
  // 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]);
}

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

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

Definition at line 245 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];
}

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

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

Definition at line 225 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 std::sqrt(error_vals[0]);
}

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

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

Definition at line 265 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];
}

---

Member Data Documentation

const EquationSystems& ExactSolution::_equation_systems [private]

Constant reference to the EquationSystems object used for the simulation.

Definition at line 259 of file exact_solution.h.

Referenced by _compute_error(), and ExactSolution().

const EquationSystems* ExactSolution::_equation_systems_fine [private]

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

Definition at line 265 of file exact_solution.h.

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

std::map<std::string, SystemErrorMap> 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 253 of file exact_solution.h.

Referenced by _check_inputs(), and ExactSolution().

Gradient(* ExactSolution::_exact_deriv)(const Point &p, const Parameters &parameters, const std::string &sys_name, const std::string &unknown_name) [private]

Function pointer to user-provided function which computes the exact derivative of the solution.

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

Tensor(* ExactSolution::_exact_hessian)(const Point &p, const Parameters &parameters, const std::string &sys_name, const std::string &unknown_name) [private]

Function pointer to user-provided function which computes the exact hessian of the solution.

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

Number(* ExactSolution::_exact_value)(const Point &p, const Parameters &parameters, const std::string &sys_name, const std::string &unknown_name) [private]

Function pointer to user-provided function which computes the exact value of the solution.

Referenced by _compute_error(), and attach_reference_solution().

int ExactSolution::_extra_order [private]

Extra order to use for quadrature rule

Definition at line 270 of file exact_solution.h.

Referenced by _compute_error(), and extra_quadrature_order().

---

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

• exact_solution.h
• exact_solution.C