libMesh::ExactSolution Class Reference

#include <exact_solution.h>

Public Member Functions

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

Private Types

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

Private Member Functions

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

Private Attributes

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

Detailed Description

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

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

Member Typedef Documentation

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

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

Definition at line 313 of file exact_solution.h.

Constructor & Destructor Documentation

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

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

Definition at line 41 of file exact_solution.C.

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

41  :
44  _extra_order(0)
45 {
46  // Initialize the _errors data structure which holds all
47  // the eventual values of the error.
48  for (unsigned int sys=0; sys<_equation_systems.n_systems(); ++sys)
49  {
50  // Reference to the system
51  const System & system = _equation_systems.get_system(sys);
52 
53  // The name of the system
54  const std::string & sys_name = system.name();
55 
56  // The SystemErrorMap to be inserted
58 
59  for (unsigned int var=0; var<system.n_vars(); ++var)
60  {
61  // The name of this variable
62  const std::string & var_name = system.variable_name(var);
63  sem[var_name] = std::vector<Real>(5, 0.);
64  }
65 
66  _errors[sys_name] = sem;
67  }
68 }
const EquationSystems & _equation_systems
ImplicitSystem & sys
const class libmesh_nullptr_t libmesh_nullptr
const EquationSystems * _equation_systems_fine
std::map< std::string, std::vector< Real > > SystemErrorMap
unsigned int n_systems() const
std::map< std::string, SystemErrorMap > _errors
const T_sys & get_system(const std::string &name) const
libMesh::ExactSolution::~ExactSolution ( )

Destructor.

Definition at line 71 of file exact_solution.C.

References _exact_derivs, _exact_hessians, and _exact_values.

72 {
73  // delete will clean up any cloned functors and no-op on any NULL
74  // pointers
75 
76  for (std::size_t i=0; i != _exact_values.size(); ++i)
77  delete (_exact_values[i]);
78 
79  for (std::size_t i=0; i != _exact_derivs.size(); ++i)
80  delete (_exact_derivs[i]);
81 
82  for (std::size_t i=0; i != _exact_hessians.size(); ++i)
83  delete (_exact_hessians[i]);
84 }
std::vector< FunctionBase< Tensor > * > _exact_hessians
std::vector< FunctionBase< Number > * > _exact_values
std::vector< FunctionBase< Gradient > * > _exact_derivs

Member Function Documentation

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

This function is responsible for checking the validity of the sys_name and unknown_name inputs.

Returns
A reference to the proper vector for storing the values.

Definition at line 261 of file exact_solution.C.

References _errors.

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

263 {
264  // If no exact solution function or fine grid solution has been
265  // attached, we now just compute the solution norm (i.e. the
266  // difference from an "exact solution" of zero
267 
268  // Make sure the requested sys_name exists.
269  std::map<std::string, SystemErrorMap>::iterator sys_iter =
270  _errors.find(sys_name);
271 
272  if (sys_iter == _errors.end())
273  libmesh_error_msg("Sorry, couldn't find the requested system '" << sys_name << "'.");
274 
275  // Make sure the requested unknown_name exists.
276  SystemErrorMap::iterator var_iter = (*sys_iter).second.find(unknown_name);
277 
278  if (var_iter == (*sys_iter).second.end())
279  libmesh_error_msg("Sorry, couldn't find the requested variable '" << unknown_name << "'.");
280 
281  // Return a reference to the proper error entry
282  return (*var_iter).second;
283 }
std::map< std::string, SystemErrorMap > _errors
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 508 of file exact_solution.C.

References _equation_systems, _equation_systems_fine, _exact_derivs, _exact_hessians, _exact_values, _extra_order, libMesh::MeshBase::active_local_elements_begin(), libMesh::MeshBase::active_local_elements_end(), libMesh::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(), libMesh::Elem::dim(), libMesh::TensorTools::div_from_grad(), libMesh::DofMap::dof_indices(), libMesh::MeshBase::elem_dimensions(), libMesh::err, 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::System::get_mesh(), libMesh::FEGenericBase< OutputType >::get_phi(), libMesh::EquationSystems::get_system(), libMesh::FEAbstract::get_xyz(), libMesh::libmesh_assert(), libmesh_nullptr, libMesh::QBase::n_points(), libMesh::FEInterface::n_vec_dim(), libMesh::TensorTools::norm_sq(), libMesh::TypeVector< T >::norm_sq(), libMesh::System::number(), libMesh::Real, libMesh::FEAbstract::reinit(), libMesh::SERIAL, libMesh::System::solution, libMesh::Elem::subdomain_id(), 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().

Referenced by extra_quadrature_order().

511 {
512  // Make sure we aren't "overconfigured"
514 
515  // We need a communicator.
516  const Parallel::Communicator & communicator(_equation_systems.comm());
517 
518  // This function must be run on all processors at once
519  libmesh_parallel_only(communicator);
520 
521  // Get a reference to the system whose error is being computed.
522  // If we have a fine grid, however, we'll integrate on that instead
523  // for more accuracy.
524  const System & computed_system = _equation_systems_fine ?
526  _equation_systems.get_system (sys_name);
527 
528  const Real time = _equation_systems.get_system(sys_name).time;
529 
530  const unsigned int sys_num = computed_system.number();
531  const unsigned int var = computed_system.variable_number(unknown_name);
532  const unsigned int var_component =
533  computed_system.variable_scalar_number(var, 0);
534 
535  // Prepare a global solution and a MeshFunction of the coarse system if we need one
536  UniquePtr<MeshFunction> coarse_values;
537  UniquePtr<NumericVector<Number> > comparison_soln = NumericVector<Number>::build(_equation_systems.comm());
539  {
540  const System & comparison_system
541  = _equation_systems.get_system(sys_name);
542 
543  std::vector<Number> global_soln;
544  comparison_system.update_global_solution(global_soln);
545  comparison_soln->init(comparison_system.solution->size(), true, SERIAL);
546  (*comparison_soln) = global_soln;
547 
548  coarse_values = UniquePtr<MeshFunction>
549  (new MeshFunction(_equation_systems,
550  *comparison_soln,
551  comparison_system.get_dof_map(),
552  comparison_system.variable_number(unknown_name)));
553  coarse_values->init();
554  }
555 
556  // Initialize any functors we're going to use
557  for (std::size_t i=0; i != _exact_values.size(); ++i)
558  if (_exact_values[i])
559  _exact_values[i]->init();
560 
561  for (std::size_t i=0; i != _exact_derivs.size(); ++i)
562  if (_exact_derivs[i])
563  _exact_derivs[i]->init();
564 
565  for (std::size_t i=0; i != _exact_hessians.size(); ++i)
566  if (_exact_hessians[i])
567  _exact_hessians[i]->init();
568 
569  // Get a reference to the dofmap and mesh for that system
570  const DofMap & computed_dof_map = computed_system.get_dof_map();
571 
572  const MeshBase & _mesh = computed_system.get_mesh();
573 
574  // Grab which element dimensions are present in the mesh
575  const std::set<unsigned char> & elem_dims = _mesh.elem_dimensions();
576 
577  // Zero the error before summation
578  // 0 - sum of square of function error (L2)
579  // 1 - sum of square of gradient error (H1 semi)
580  // 2 - sum of square of Hessian error (H2 semi)
581  // 3 - sum of sqrt(square of function error) (L1)
582  // 4 - max of sqrt(square of function error) (Linfty)
583  // 5 - sum of square of curl error (HCurl semi)
584  // 6 - sum of square of div error (HDiv semi)
585  error_vals = std::vector<Real>(7, 0.);
586 
587  // Construct Quadrature rule based on default quadrature order
588  const FEType & fe_type = computed_dof_map.variable_type(var);
589 
590  unsigned int n_vec_dim = FEInterface::n_vec_dim( _mesh, fe_type );
591 
592  // FIXME: MeshFunction needs to be updated to support vector-valued
593  // elements before we can use a reference solution.
594  if ((n_vec_dim > 1) && _equation_systems_fine)
595  {
596  libMesh::err << "Error calculation using reference solution not yet\n"
597  << "supported for vector-valued elements."
598  << std::endl;
599  libmesh_not_implemented();
600  }
601 
602 
603  // Allow space for dims 0-3, even if we don't use them all
604  std::vector<FEGenericBase<OutputShape> *> fe_ptrs(4, libmesh_nullptr);
605  std::vector<QBase *> q_rules(4, libmesh_nullptr);
606 
607  // Prepare finite elements for each dimension present in the mesh
608  for (std::set<unsigned char>::const_iterator d_it = elem_dims.begin();
609  d_it != elem_dims.end(); ++d_it)
610  {
611  q_rules[*d_it] =
612  fe_type.default_quadrature_rule (*d_it, _extra_order).release();
613 
614  // Construct finite element object
615 
616  fe_ptrs[*d_it] = FEGenericBase<OutputShape>::build(*d_it, fe_type).release();
617 
618  // Attach quadrature rule to FE object
619  fe_ptrs[*d_it]->attach_quadrature_rule (q_rules[*d_it]);
620  }
621 
622  // The global degree of freedom indices associated
623  // with the local degrees of freedom.
624  std::vector<dof_id_type> dof_indices;
625 
626 
627  //
628  // Begin the loop over the elements
629  //
630  // TODO: this ought to be threaded (and using subordinate
631  // MeshFunction objects in each thread rather than a single
632  // master)
633  MeshBase::const_element_iterator el = _mesh.active_local_elements_begin();
634  const MeshBase::const_element_iterator end_el = _mesh.active_local_elements_end();
635 
636  for ( ; el != end_el; ++el)
637  {
638  // Store a pointer to the element we are currently
639  // working on. This allows for nicer syntax later.
640  const Elem * elem = *el;
641  const unsigned int dim = elem->dim();
642 
643  const subdomain_id_type elem_subid = elem->subdomain_id();
644 
645  // If the variable is not active on this subdomain, don't bother
646  if (!computed_system.variable(var).active_on_subdomain(elem_subid))
647  continue;
648 
649  /* If the variable is active, then we're going to restrict the
650  MeshFunction evaluations to the current element subdomain.
651  This is for cases such as mixed dimension meshes where we want
652  to restrict the calculation to one particular domain. */
653  std::set<subdomain_id_type> subdomain_id;
654  subdomain_id.insert(elem_subid);
655 
656  FEGenericBase<OutputShape> * fe = fe_ptrs[dim];
657  QBase * qrule = q_rules[dim];
658  libmesh_assert(fe);
659  libmesh_assert(qrule);
660 
661  // The Jacobian*weight at the quadrature points.
662  const std::vector<Real> & JxW = fe->get_JxW();
663 
664  // The value of the shape functions at the quadrature points
665  // i.e. phi(i) = phi_values[i][qp]
666  const std::vector<std::vector<OutputShape> > & phi_values = fe->get_phi();
667 
668  // The value of the shape function gradients at the quadrature points
669  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputGradient> > &
670  dphi_values = fe->get_dphi();
671 
672  // The value of the shape function curls at the quadrature points
673  // Only computed for vector-valued elements
674  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputShape> > * curl_values = libmesh_nullptr;
675 
676  // The value of the shape function divergences at the quadrature points
677  // Only computed for vector-valued elements
678  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputDivergence> > * div_values = libmesh_nullptr;
679 
680  if (FEInterface::field_type(fe_type) == TYPE_VECTOR)
681  {
682  curl_values = &fe->get_curl_phi();
683  div_values = &fe->get_div_phi();
684  }
685 
686 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
687  // The value of the shape function second derivatives at the quadrature points
688  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputTensor> > &
689  d2phi_values = fe->get_d2phi();
690 #endif
691 
692  // The XYZ locations (in physical space) of the quadrature points
693  const std::vector<Point> & q_point = fe->get_xyz();
694 
695  // reinitialize the element-specific data
696  // for the current element
697  fe->reinit (elem);
698 
699  // Get the local to global degree of freedom maps
700  computed_dof_map.dof_indices (elem, dof_indices, var);
701 
702  // The number of quadrature points
703  const unsigned int n_qp = qrule->n_points();
704 
705  // The number of shape functions
706  const unsigned int n_sf =
707  cast_int<unsigned int>(dof_indices.size());
708 
709  //
710  // Begin the loop over the Quadrature points.
711  //
712  for (unsigned int qp=0; qp<n_qp; qp++)
713  {
714  // Real u_h = 0.;
715  // RealGradient grad_u_h;
716 
718 
720 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
722 #endif
723  typename FEGenericBase<OutputShape>::OutputNumber curl_u_h(0.0);
725 
726  // Compute solution values at the current
727  // quadrature point. This requires a sum
728  // over all the shape functions evaluated
729  // at the quadrature point.
730  for (unsigned int i=0; i<n_sf; i++)
731  {
732  // Values from current solution.
733  u_h += phi_values[i][qp]*computed_system.current_solution (dof_indices[i]);
734  grad_u_h += dphi_values[i][qp]*computed_system.current_solution (dof_indices[i]);
735 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
736  grad2_u_h += d2phi_values[i][qp]*computed_system.current_solution (dof_indices[i]);
737 #endif
738  if (FEInterface::field_type(fe_type) == TYPE_VECTOR)
739  {
740  curl_u_h += (*curl_values)[i][qp]*computed_system.current_solution (dof_indices[i]);
741  div_u_h += (*div_values)[i][qp]*computed_system.current_solution (dof_indices[i]);
742  }
743  }
744 
745  // Compute the value of the error at this quadrature point
746  typename FEGenericBase<OutputShape>::OutputNumber exact_val(0);
747  RawAccessor<typename FEGenericBase<OutputShape>::OutputNumber> exact_val_accessor( exact_val, dim );
748  if (_exact_values.size() > sys_num && _exact_values[sys_num])
749  {
750  for (unsigned int c = 0; c < n_vec_dim; c++)
751  exact_val_accessor(c) =
752  _exact_values[sys_num]->
753  component(var_component+c, q_point[qp], time);
754  }
755  else if (_equation_systems_fine)
756  {
757  // FIXME: Needs to be updated for vector-valued elements
758  DenseVector<Number> output(1);
759  (*coarse_values)(q_point[qp],time,output,&subdomain_id);
760  exact_val = output(0);
761  }
762  const typename FEGenericBase<OutputShape>::OutputNumber val_error = u_h - exact_val;
763 
764  // Add the squares of the error to each contribution
765  Real error_sq = TensorTools::norm_sq(val_error);
766  error_vals[0] += JxW[qp]*error_sq;
767 
768  Real norm = sqrt(error_sq);
769  error_vals[3] += JxW[qp]*norm;
770 
771  if (error_vals[4]<norm) { error_vals[4] = norm; }
772 
773  // Compute the value of the error in the gradient at this
774  // quadrature point
776  RawAccessor<typename FEGenericBase<OutputShape>::OutputNumberGradient> exact_grad_accessor( exact_grad, LIBMESH_DIM );
777  if (_exact_derivs.size() > sys_num && _exact_derivs[sys_num])
778  {
779  for (unsigned int c = 0; c < n_vec_dim; c++)
780  for (unsigned int d = 0; d < LIBMESH_DIM; d++)
781  exact_grad_accessor(d + c*LIBMESH_DIM) =
782  _exact_derivs[sys_num]->
783  component(var_component+c, q_point[qp], time)(d);
784  }
785  else if (_equation_systems_fine)
786  {
787  // FIXME: Needs to be updated for vector-valued elements
788  std::vector<Gradient> output(1);
789  coarse_values->gradient(q_point[qp],time,output,&subdomain_id);
790  exact_grad = output[0];
791  }
792 
793  const typename FEGenericBase<OutputShape>::OutputNumberGradient grad_error = grad_u_h - exact_grad;
794 
795  error_vals[1] += JxW[qp]*grad_error.norm_sq();
796 
797 
798  if (FEInterface::field_type(fe_type) == TYPE_VECTOR)
799  {
800  // Compute the value of the error in the curl at this
801  // quadrature point
802  typename FEGenericBase<OutputShape>::OutputNumber exact_curl(0.0);
803  if (_exact_derivs.size() > sys_num && _exact_derivs[sys_num])
804  {
805  exact_curl = TensorTools::curl_from_grad( exact_grad );
806  }
807  else if (_equation_systems_fine)
808  {
809  // FIXME: Need to implement curl for MeshFunction and support reference
810  // solution for vector-valued elements
811  }
812 
813  const typename FEGenericBase<OutputShape>::OutputNumber curl_error = curl_u_h - exact_curl;
814 
815  error_vals[5] += JxW[qp]*TensorTools::norm_sq(curl_error);
816 
817  // Compute the value of the error in the divergence at this
818  // quadrature point
820  if (_exact_derivs.size() > sys_num && _exact_derivs[sys_num])
821  {
822  exact_div = TensorTools::div_from_grad( exact_grad );
823  }
824  else if (_equation_systems_fine)
825  {
826  // FIXME: Need to implement div for MeshFunction and support reference
827  // solution for vector-valued elements
828  }
829 
830  const typename FEGenericBase<OutputShape>::OutputNumberDivergence div_error = div_u_h - exact_div;
831 
832  error_vals[6] += JxW[qp]*TensorTools::norm_sq(div_error);
833  }
834 
835 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
836  // Compute the value of the error in the hessian at this
837  // quadrature point
839  RawAccessor<typename FEGenericBase<OutputShape>::OutputNumberTensor> exact_hess_accessor( exact_hess, dim );
840  if (_exact_hessians.size() > sys_num && _exact_hessians[sys_num])
841  {
842  //FIXME: This needs to be implemented to support rank 3 tensors
843  // which can't happen until type_n_tensor is fully implemented
844  // and a RawAccessor<TypeNTensor> is fully implemented
845  if (FEInterface::field_type(fe_type) == TYPE_VECTOR)
846  libmesh_not_implemented();
847 
848  for (unsigned int c = 0; c < n_vec_dim; c++)
849  for (unsigned int d = 0; d < dim; d++)
850  for (unsigned int e =0; e < dim; e++)
851  exact_hess_accessor(d + e*dim + c*dim*dim) =
852  _exact_hessians[sys_num]->
853  component(var_component+c, q_point[qp], time)(d,e);
854  }
855  else if (_equation_systems_fine)
856  {
857  // FIXME: Needs to be updated for vector-valued elements
858  std::vector<Tensor> output(1);
859  coarse_values->hessian(q_point[qp],time,output,&subdomain_id);
860  exact_hess = output[0];
861  }
862 
863  const typename FEGenericBase<OutputShape>::OutputNumberTensor grad2_error = grad2_u_h - exact_hess;
864 
865  // FIXME: PB: Is this what we want for rank 3 tensors?
866  error_vals[2] += JxW[qp]*grad2_error.norm_sq();
867 #endif
868 
869  } // end qp loop
870  } // end element loop
871 
872  // Clean up the FE and QBase pointers we created
873  for (std::set<unsigned char>::const_iterator d_it = elem_dims.begin();
874  d_it != elem_dims.end(); ++d_it)
875  {
876  delete fe_ptrs[*d_it];
877  delete q_rules[*d_it];
878  }
879 
880  // Add up the error values on all processors, except for the L-infty
881  // norm, for which the maximum is computed.
882  Real l_infty_norm = error_vals[4];
883  communicator.max(l_infty_norm);
884  communicator.sum(error_vals);
885  error_vals[4] = l_infty_norm;
886 }
std::vector< FunctionBase< Tensor > * > _exact_hessians
TensorTools::DecrementRank< OutputNumber >::type OutputNumberDivergence
Definition: fe_base.h:127
Number div_from_grad(const VectorValue< Number > &)
Dummy. Divergence of a scalar not defined, but is needed for ExactSolution to compile.
Definition: tensor_tools.C:54
std::vector< FunctionBase< Number > * > _exact_values
const EquationSystems & _equation_systems
TestClass subdomain_id_type
Definition: id_types.h:43
MPI_Comm communicator
Definition: parallel.h:177
TensorTools::IncrementRank< OutputNumberGradient >::type OutputNumberTensor
Definition: fe_base.h:126
static FEFieldType field_type(const FEType &fe_type)
const class libmesh_nullptr_t libmesh_nullptr
TensorTools::IncrementRank< OutputNumber >::type OutputNumberGradient
Definition: fe_base.h:125
std::vector< FunctionBase< Gradient > * > _exact_derivs
const EquationSystems * _equation_systems_fine
libmesh_assert(j)
static UniquePtr< NumericVector< T > > build(const Parallel::Communicator &comm, const SolverPackage solver_package=libMesh::default_solver_package())
OStreamProxy err(std::cerr)
TensorTools::MakeNumber< OutputShape >::type OutputNumber
Definition: fe_base.h:124
static unsigned int n_vec_dim(const MeshBase &mesh, const FEType &fe_type)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const Parallel::Communicator & comm() const
const T_sys & get_system(const std::string &name) const
Number curl_from_grad(const VectorValue< Number > &)
Definition: tensor_tools.C:28
static UniquePtr< FEGenericBase > build(const unsigned int dim, const FEType &type)
void libMesh::ExactSolution::attach_exact_deriv ( unsigned int  sys_num,
FunctionBase< Gradient > *  g 
)

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

Definition at line 196 of file exact_solution.C.

References _equation_systems, _equation_systems_fine, _exact_derivs, _exact_hessians, attach_exact_hessian(), libMesh::FunctionBase< Output >::clone(), libMesh::EquationSystems::get_system(), libMesh::libmesh_assert(), libmesh_nullptr, libMesh::EquationSystems::n_systems(), libMesh::EquationSystems::parameters, and libMesh::sys.

Referenced by attach_exact_value().

198 {
199  if (_exact_derivs.size() <= sys_num)
200  _exact_derivs.resize(sys_num+1, libmesh_nullptr);
201 
202  if (g)
203  _exact_derivs[sys_num] = g->clone().release();
204 }
const class libmesh_nullptr_t libmesh_nullptr
std::vector< FunctionBase< Gradient > * > _exact_derivs
void libMesh::ExactSolution::attach_exact_deriv ( Gradient   gptrconst Point &p, const Parameters &parameters, const std::string &sys_name, const std::string &unknown_name)

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

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

References _exact_derivs, and libmesh_nullptr.

178 {
179  // Clear out any previous _exact_derivs entries, then add a new
180  // entry for each system.
181  for (std::size_t i=0; i != _exact_derivs.size(); ++i)
182  delete (_exact_derivs[i]);
183 
184  _exact_derivs.clear();
185  _exact_derivs.resize(g.size(), libmesh_nullptr);
186 
187  // We use clone() to get non-sliced copies of FunctionBase
188  // subclasses, but we don't currently put the resulting UniquePtrs
189  // into an STL container.
190  for (std::size_t i=0; i != g.size(); ++i)
191  if (g[i])
192  _exact_derivs[i] = g[i]->clone().release();
193 }
const class libmesh_nullptr_t libmesh_nullptr
std::vector< FunctionBase< Gradient > * > _exact_derivs
void libMesh::ExactSolution::attach_exact_hessian ( unsigned int  sys_num,
FunctionBase< Tensor > *  h 
)

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

Definition at line 250 of file exact_solution.C.

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

Referenced by attach_exact_deriv().

252 {
253  if (_exact_hessians.size() <= sys_num)
254  _exact_hessians.resize(sys_num+1, libmesh_nullptr);
255 
256  if (h)
257  _exact_hessians[sys_num] = h->clone().release();
258 }
std::vector< FunctionBase< Tensor > * > _exact_hessians
const class libmesh_nullptr_t libmesh_nullptr
void libMesh::ExactSolution::attach_exact_hessian ( Tensor   hptrconst Point &p, const Parameters &parameters, const std::string &sys_name, const std::string &unknown_name)

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

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

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

Definition at line 231 of file exact_solution.C.

References _exact_hessians, and libmesh_nullptr.

232 {
233  // Clear out any previous _exact_hessians entries, then add a new
234  // entry for each system.
235  for (std::size_t i=0; i != _exact_hessians.size(); ++i)
236  delete (_exact_hessians[i]);
237 
238  _exact_hessians.clear();
239  _exact_hessians.resize(h.size(), libmesh_nullptr);
240 
241  // We use clone() to get non-sliced copies of FunctionBase
242  // subclasses, but we don't currently put the resulting UniquePtrs
243  // into an STL container.
244  for (std::size_t i=0; i != h.size(); ++i)
245  if (h[i])
246  _exact_hessians[i] = h[i]->clone().release();
247 }
std::vector< FunctionBase< Tensor > * > _exact_hessians
const class libmesh_nullptr_t libmesh_nullptr
void libMesh::ExactSolution::attach_exact_value ( unsigned int  sys_num,
FunctionBase< Number > *  f 
)

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

Definition at line 142 of file exact_solution.C.

References _equation_systems, _equation_systems_fine, _exact_derivs, _exact_values, attach_exact_deriv(), libMesh::FunctionBase< Output >::clone(), libMesh::EquationSystems::get_system(), libMesh::libmesh_assert(), libmesh_nullptr, libMesh::EquationSystems::n_systems(), libMesh::EquationSystems::parameters, and libMesh::sys.

Referenced by attach_reference_solution().

144 {
145  if (_exact_values.size() <= sys_num)
146  _exact_values.resize(sys_num+1, libmesh_nullptr);
147 
148  if (f)
149  _exact_values[sys_num] = f->clone().release();
150 }
std::vector< FunctionBase< Number > * > _exact_values
const class libmesh_nullptr_t libmesh_nullptr
virtual UniquePtr< FunctionBase< Output > > clone() const =0
void libMesh::ExactSolution::attach_exact_value ( Number   fptrconst Point &p, const Parameters &Parameters, const std::string &sys_name, const std::string &unknown_name)

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

void libMesh::ExactSolution::attach_exact_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 123 of file exact_solution.C.

References _exact_values, and libmesh_nullptr.

124 {
125  // Clear out any previous _exact_values entries, then add a new
126  // entry for each system.
127  for (std::size_t i=0; i != _exact_values.size(); ++i)
128  delete (_exact_values[i]);
129 
130  _exact_values.clear();
131  _exact_values.resize(f.size(), libmesh_nullptr);
132 
133  // We use clone() to get non-sliced copies of FunctionBase
134  // subclasses, but we don't currently put the resulting UniquePtrs
135  // into an STL container.
136  for (std::size_t i=0; i != f.size(); ++i)
137  if (f[i])
138  _exact_values[i] = f[i]->clone().release();
139 }
std::vector< FunctionBase< Number > * > _exact_values
const class libmesh_nullptr_t libmesh_nullptr
void libMesh::ExactSolution::attach_reference_solution ( const EquationSystems es_fine)

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

Definition at line 87 of file exact_solution.C.

References _equation_systems, _equation_systems_fine, _exact_derivs, _exact_hessians, _exact_values, attach_exact_value(), libMesh::EquationSystems::get_system(), libMesh::libmesh_assert(), libmesh_nullptr, libMesh::EquationSystems::n_systems(), libMesh::EquationSystems::parameters, and libMesh::sys.

88 {
89  libmesh_assert(es_fine);
90  _equation_systems_fine = es_fine;
91 
92  // If we're using a fine grid solution, we're not using exact values
93  _exact_values.clear();
94  _exact_derivs.clear();
95  _exact_hessians.clear();
96 }
std::vector< FunctionBase< Tensor > * > _exact_hessians
std::vector< FunctionBase< Number > * > _exact_values
std::vector< FunctionBase< Gradient > * > _exact_derivs
const EquationSystems * _equation_systems_fine
libmesh_assert(j)
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 287 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::sys, libMesh::TYPE_SCALAR, and libMesh::TYPE_VECTOR.

Referenced by extra_quadrature_order().

289 {
290  // Check the inputs for validity, and get a reference
291  // to the proper location to store the error
292  std::vector<Real> & error_vals = this->_check_inputs(sys_name,
293  unknown_name);
294 
296  const System & sys = _equation_systems.get_system<System>( sys_name );
297 
298  libmesh_assert( sys.has_variable( unknown_name ) );
299  switch( FEInterface::field_type(sys.variable_type( unknown_name )) )
300  {
301  case TYPE_SCALAR:
302  {
303  this->_compute_error<Real>(sys_name,
304  unknown_name,
305  error_vals);
306  break;
307  }
308  case TYPE_VECTOR:
309  {
310  this->_compute_error<RealGradient>(sys_name,
311  unknown_name,
312  error_vals);
313  break;
314  }
315  default:
316  libmesh_error_msg("Invalid variable type!");
317  }
318 }
bool has_system(const std::string &name) const
const EquationSystems & _equation_systems
ImplicitSystem & sys
static FEFieldType field_type(const FEType &fe_type)
libmesh_assert(j)
const T_sys & get_system(const std::string &name) const
std::vector< Real > & _check_inputs(const std::string &sys_name, const std::string &unknown_name)
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 324 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(), and libMesh::TYPE_SCALAR.

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

327 {
328  // Check the inputs for validity, and get a reference
329  // to the proper location to store the error
330  std::vector<Real> & error_vals = this->_check_inputs(sys_name,
331  unknown_name);
332 
334  libmesh_assert(_equation_systems.get_system(sys_name).has_variable( unknown_name ));
335  const FEType & fe_type = _equation_systems.get_system(sys_name).variable_type(unknown_name);
336 
337  switch (norm)
338  {
339  case L2:
340  return std::sqrt(error_vals[0]);
341  case H1:
342  return std::sqrt(error_vals[0] + error_vals[1]);
343  case H2:
344  return std::sqrt(error_vals[0] + error_vals[1] + error_vals[2]);
345  case HCURL:
346  {
347  if (FEInterface::field_type(fe_type) == TYPE_SCALAR)
348  libmesh_error_msg("Cannot compute HCurl error norm of scalar-valued variables!");
349  else
350  return std::sqrt(error_vals[0] + error_vals[5]);
351  }
352  case HDIV:
353  {
354  if (FEInterface::field_type(fe_type) == TYPE_SCALAR)
355  libmesh_error_msg("Cannot compute HDiv error norm of scalar-valued variables!");
356  else
357  return std::sqrt(error_vals[0] + error_vals[6]);
358  }
359  case H1_SEMINORM:
360  return std::sqrt(error_vals[1]);
361  case H2_SEMINORM:
362  return std::sqrt(error_vals[2]);
363  case HCURL_SEMINORM:
364  {
365  if (FEInterface::field_type(fe_type) == TYPE_SCALAR)
366  libmesh_error_msg("Cannot compute HCurl error seminorm of scalar-valued variables!");
367  else
368  return std::sqrt(error_vals[5]);
369  }
370  case HDIV_SEMINORM:
371  {
372  if (FEInterface::field_type(fe_type) == TYPE_SCALAR)
373  libmesh_error_msg("Cannot compute HDiv error seminorm of scalar-valued variables!");
374  else
375  return std::sqrt(error_vals[6]);
376  }
377  case L1:
378  return error_vals[3];
379  case L_INF:
380  return error_vals[4];
381 
382  default:
383  libmesh_error_msg("Currently only Sobolev norms/seminorms are supported!");
384  }
385 }
bool has_system(const std::string &name) const
const EquationSystems & _equation_systems
static FEFieldType field_type(const FEType &fe_type)
libmesh_assert(j)
const T_sys & get_system(const std::string &name) const
std::vector< Real > & _check_inputs(const std::string &sys_name, const std::string &unknown_name)
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 158 of file exact_solution.h.

References _check_inputs(), _compute_error(), _extra_order, compute_error(), error_norm(), h1_error(), h2_error(), hcurl_error(), hdiv_error(), l1_error(), l2_error(), l_inf_error(), and libMesh::Real.

159  { _extra_order = extraorder; }
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 453 of file exact_solution.C.

References _check_inputs().

Referenced by extra_quadrature_order().

455 {
456  // If the user has supplied no exact derivative function, we
457  // just integrate the H1 norm of the solution; i.e. its
458  // difference from an "exact solution" of zero.
459 
460  // Check the inputs for validity, and get a reference
461  // to the proper location to store the error
462  std::vector<Real> & error_vals = this->_check_inputs(sys_name,
463  unknown_name);
464 
465  // Return the square root of the sum of the computed errors.
466  return std::sqrt(error_vals[0] + error_vals[1]);
467 }
std::vector< Real > & _check_inputs(const std::string &sys_name, const std::string &unknown_name)
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 485 of file exact_solution.C.

References _check_inputs().

Referenced by extra_quadrature_order().

487 {
488  // If the user has supplied no exact derivative functions, we
489  // just integrate the H2 norm of the solution; i.e. its
490  // difference from an "exact solution" of zero.
491 
492  // Check the inputs for validity, and get a reference
493  // to the proper location to store the error
494  std::vector<Real> & error_vals = this->_check_inputs(sys_name,
495  unknown_name);
496 
497  // Return the square root of the sum of the computed errors.
498  return std::sqrt(error_vals[0] + error_vals[1] + error_vals[2]);
499 }
std::vector< Real > & _check_inputs(const std::string &sys_name, const std::string &unknown_name)
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 470 of file exact_solution.C.

References error_norm(), and libMesh::HCURL.

Referenced by extra_quadrature_order().

472 {
473  return this->error_norm(sys_name,unknown_name,HCURL);
474 }
Real error_norm(const std::string &sys_name, const std::string &unknown_name, const FEMNormType &norm)
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 477 of file exact_solution.C.

References error_norm(), and libMesh::HDIV.

Referenced by extra_quadrature_order().

479 {
480  return this->error_norm(sys_name,unknown_name,HDIV);
481 }
Real error_norm(const std::string &sys_name, const std::string &unknown_name, const FEMNormType &norm)
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 413 of file exact_solution.C.

References _check_inputs().

Referenced by extra_quadrature_order().

415 {
416 
417  // Check the inputs for validity, and get a reference
418  // to the proper location to store the error
419  std::vector<Real> & error_vals = this->_check_inputs(sys_name,
420  unknown_name);
421 
422  // Return the square root of the first component of the
423  // computed error.
424  return error_vals[3];
425 }
std::vector< Real > & _check_inputs(const std::string &sys_name, const std::string &unknown_name)
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 393 of file exact_solution.C.

References _check_inputs().

Referenced by extra_quadrature_order().

395 {
396 
397  // Check the inputs for validity, and get a reference
398  // to the proper location to store the error
399  std::vector<Real> & error_vals = this->_check_inputs(sys_name,
400  unknown_name);
401 
402  // Return the square root of the first component of the
403  // computed error.
404  return std::sqrt(error_vals[0]);
405 }
std::vector< Real > & _check_inputs(const std::string &sys_name, const std::string &unknown_name)
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 433 of file exact_solution.C.

References _check_inputs().

Referenced by extra_quadrature_order().

435 {
436 
437  // Check the inputs for validity, and get a reference
438  // to the proper location to store the error
439  std::vector<Real> & error_vals = this->_check_inputs(sys_name,
440  unknown_name);
441 
442  // Return the square root of the first component of the
443  // computed error.
444  return error_vals[4];
445 }
std::vector< Real > & _check_inputs(const std::string &sys_name, const std::string &unknown_name)

Member Data Documentation

const EquationSystems& libMesh::ExactSolution::_equation_systems
private

Constant reference to the EquationSystems object used for the simulation.

Definition at line 327 of file exact_solution.h.

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

const EquationSystems* libMesh::ExactSolution::_equation_systems_fine
private

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

Definition at line 333 of file exact_solution.h.

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

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

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

Definition at line 321 of file exact_solution.h.

Referenced by _check_inputs(), and ExactSolution().

std::vector<FunctionBase<Gradient> *> libMesh::ExactSolution::_exact_derivs
private

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

Definition at line 297 of file exact_solution.h.

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

std::vector<FunctionBase<Tensor> *> libMesh::ExactSolution::_exact_hessians
private

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

Definition at line 303 of file exact_solution.h.

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

std::vector<FunctionBase<Number> *> libMesh::ExactSolution::_exact_values
private

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

Definition at line 291 of file exact_solution.h.

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

int libMesh::ExactSolution::_extra_order
private

Extra order to use for quadrature rule

Definition at line 338 of file exact_solution.h.

Referenced by _compute_error(), and extra_quadrature_order().


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