libMesh::ExactSolution Class Reference

#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)
 
typedef Gradient(* GradientFunctionPointer) (const Point &p, const Parameters &parameters, const std::string &sys_name, const std::string &unknown_name)
 
typedef Tensor(* HessianFunctionPointer) (const Point &p, const Parameters &parameters, const std::string &sys_name, const std::string &unknown_name)
 

Public Member Functions

 ExactSolution (const EquationSystems &es)
 
 ExactSolution (const ExactSolution &)=delete
 
ExactSolutionoperator= (const ExactSolution &)=delete
 
ExactSolutionoperator= (ExactSolution &&)=delete
 
 ExactSolution (ExactSolution &&)
 
 ~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 (ValueFunctionPointer fptr)
 
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 (GradientFunctionPointer gptr)
 
void attach_exact_hessians (std::vector< FunctionBase< Tensor > *> h)
 
void attach_exact_hessian (unsigned int sys_num, FunctionBase< Tensor > *h)
 
void attach_exact_hessian (HessianFunctionPointer hptr)
 
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< std::unique_ptr< FunctionBase< Number > > > _exact_values
 
std::vector< std::unique_ptr< FunctionBase< Gradient > > > _exact_derivs
 
std::vector< std::unique_ptr< 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 72 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 148 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 171 of file exact_solution.h.

◆ SystemErrorMap

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

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

Definition at line 336 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 125 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 39 of file exact_solution.C.

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

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

◆ ExactSolution() [2/3]

libMesh::ExactSolution::ExactSolution ( const ExactSolution )
delete

The copy constructor and copy/move assignment operators are deleted. This class has containers of unique_ptrs so it can't be default (shallow) copied, and it has a const reference so it can't be assigned to after creation.

◆ ExactSolution() [3/3]

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

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

◆ ~ExactSolution()

libMesh::ExactSolution::~ExactSolution ( )
default

Member Function Documentation

◆ _check_inputs()

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

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

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

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

210 {
211  // If no exact solution function or fine grid solution has been
212  // attached, we now just compute the solution norm (i.e. the
213  // difference from an "exact solution" of zero
214 
215  // Make sure the requested sys_name exists.
216  std::map<std::string, SystemErrorMap>::iterator sys_iter =
217  _errors.find(sys_name);
218 
219  if (sys_iter == _errors.end())
220  libmesh_error_msg("Sorry, couldn't find the requested system '" << sys_name << "'.");
221 
222  // Make sure the requested unknown_name exists.
223  SystemErrorMap::iterator var_iter = (*sys_iter).second.find(unknown_name);
224 
225  if (var_iter == (*sys_iter).second.end())
226  libmesh_error_msg("Sorry, couldn't find the requested variable '" << unknown_name << "'.");
227 
228  // Return a reference to the proper error entry
229  return (*var_iter).second;
230 }
std::map< std::string, SystemErrorMap > _errors

◆ _compute_error()

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

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

Definition at line 455 of file exact_solution.C.

References _equation_systems, _equation_systems_fine, _exact_derivs, _exact_hessians, _exact_values, _extra_order, libMesh::MeshBase::active_local_element_ptr_range(), 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::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::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::TYPE_VECTOR, libMesh::System::update_global_solution(), libMesh::System::variable(), libMesh::System::variable_number(), libMesh::System::variable_scalar_number(), and libMesh::DofMap::variable_type().

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

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

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

158 {
159  if (_exact_derivs.size() <= sys_num)
160  _exact_derivs.resize(sys_num+1);
161 
162  if (g)
163  _exact_derivs[sys_num] = g->clone();
164 }
std::vector< std::unique_ptr< FunctionBase< Gradient > > > _exact_derivs

◆ attach_exact_deriv() [2/2]

void libMesh::ExactSolution::attach_exact_deriv ( GradientFunctionPointer  gptr)

Definition at line 126 of file exact_solution.C.

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

127 {
128  libmesh_assert(gptr);
129 
130  // Clear out any previous _exact_derivs entries, then add a new
131  // entry for each system.
132  _exact_derivs.clear();
133 
134  for (unsigned int sys=0; sys<_equation_systems.n_systems(); ++sys)
135  {
136  const System & system = _equation_systems.get_system(sys);
137  _exact_derivs.emplace_back(libmesh_make_unique<WrappedFunction<Gradient>>(system, gptr, &_equation_systems.parameters));
138  }
139 
140  // If we're using exact values, we're not using a fine grid solution
141  _equation_systems_fine = nullptr;
142 }
unsigned int n_systems() const
const EquationSystems & _equation_systems
const T_sys & get_system(const std::string &name) const
const EquationSystems * _equation_systems_fine
std::vector< std::unique_ptr< FunctionBase< Gradient > > > _exact_derivs

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

References _exact_derivs.

146 {
147  // Automatically delete any previous _exact_derivs entries, then add a new
148  // entry for each system.
149  _exact_derivs.clear();
150 
151  for (auto ptr : g)
152  _exact_derivs.emplace_back(ptr ? ptr->clone() : nullptr);
153 }
std::vector< std::unique_ptr< FunctionBase< Gradient > > > _exact_derivs

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

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

199 {
200  if (_exact_hessians.size() <= sys_num)
201  _exact_hessians.resize(sys_num+1);
202 
203  if (h)
204  _exact_hessians[sys_num] = h->clone();
205 }
std::vector< std::unique_ptr< FunctionBase< Tensor > > > _exact_hessians

◆ attach_exact_hessian() [2/2]

void libMesh::ExactSolution::attach_exact_hessian ( HessianFunctionPointer  hptr)

Definition at line 167 of file exact_solution.C.

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

168 {
169  libmesh_assert(hptr);
170 
171  // Clear out any previous _exact_hessians entries, then add a new
172  // entry for each system.
173  _exact_hessians.clear();
174 
175  for (unsigned int sys=0; sys<_equation_systems.n_systems(); ++sys)
176  {
177  const System & system = _equation_systems.get_system(sys);
178  _exact_hessians.emplace_back(libmesh_make_unique<WrappedFunction<Tensor>>(system, hptr, &_equation_systems.parameters));
179  }
180 
181  // If we're using exact values, we're not using a fine grid solution
182  _equation_systems_fine = nullptr;
183 }
unsigned int n_systems() const
const EquationSystems & _equation_systems
const T_sys & get_system(const std::string &name) const
const EquationSystems * _equation_systems_fine
std::vector< std::unique_ptr< FunctionBase< Tensor > > > _exact_hessians

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

References _exact_hessians.

187 {
188  // Automatically delete any previous _exact_hessians entries, then add a new
189  // entry for each system.
190  _exact_hessians.clear();
191 
192  for (auto ptr : h)
193  _exact_hessians.emplace_back(ptr ? ptr->clone() : nullptr);
194 }
std::vector< std::unique_ptr< FunctionBase< Tensor > > > _exact_hessians

◆ attach_exact_value() [1/2]

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

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

Definition at line 115 of file exact_solution.C.

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

117 {
118  if (_exact_values.size() <= sys_num)
119  _exact_values.resize(sys_num+1);
120 
121  if (f)
122  _exact_values[sys_num] = f->clone();
123 }
virtual std::unique_ptr< FunctionBase< Output > > clone() const =0
std::vector< std::unique_ptr< FunctionBase< Number > > > _exact_values

◆ attach_exact_value() [2/2]

void libMesh::ExactSolution::attach_exact_value ( ValueFunctionPointer  fptr)

Definition at line 85 of file exact_solution.C.

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

86 {
87  libmesh_assert(fptr);
88 
89  // Clear out any previous _exact_values entries, then add a new
90  // entry for each system.
91  _exact_values.clear();
92 
93  for (unsigned int sys=0; sys<_equation_systems.n_systems(); ++sys)
94  {
95  const System & system = _equation_systems.get_system(sys);
96  _exact_values.emplace_back(libmesh_make_unique<WrappedFunction<Number>>(system, fptr, &_equation_systems.parameters));
97  }
98 
99  // If we're using exact values, we're not using a fine grid solution
100  _equation_systems_fine = nullptr;
101 }
unsigned int n_systems() const
const EquationSystems & _equation_systems
const T_sys & get_system(const std::string &name) const
const EquationSystems * _equation_systems_fine
std::vector< std::unique_ptr< FunctionBase< Number > > > _exact_values

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

References _exact_values.

105 {
106  // Automatically delete any previous _exact_values entries, then add a new
107  // entry for each system.
108  _exact_values.clear();
109 
110  for (auto ptr : f)
111  _exact_values.emplace_back(ptr ? ptr->clone() : nullptr);
112 }
std::vector< std::unique_ptr< FunctionBase< Number > > > _exact_values

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

References _equation_systems_fine, _exact_derivs, _exact_hessians, and _exact_values.

74 {
75  libmesh_assert(es_fine);
76  _equation_systems_fine = es_fine;
77 
78  // If we're using a fine grid solution, we're not using exact values
79  _exact_values.clear();
80  _exact_derivs.clear();
81  _exact_hessians.clear();
82 }
const EquationSystems * _equation_systems_fine
std::vector< std::unique_ptr< FunctionBase< Gradient > > > _exact_derivs
std::vector< std::unique_ptr< FunctionBase< Tensor > > > _exact_hessians
std::vector< std::unique_ptr< FunctionBase< Number > > > _exact_values

◆ compute_error()

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

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

Definition at line 234 of file exact_solution.C.

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

236 {
237  // Check the inputs for validity, and get a reference
238  // to the proper location to store the error
239  std::vector<Real> & error_vals = this->_check_inputs(sys_name,
240  unknown_name);
241 
242  libmesh_assert( _equation_systems.has_system(sys_name) );
243  const System & sys = _equation_systems.get_system<System>( sys_name );
244 
245  libmesh_assert( sys.has_variable( unknown_name ) );
246  switch( FEInterface::field_type(sys.variable_type( unknown_name )) )
247  {
248  case TYPE_SCALAR:
249  {
250  this->_compute_error<Real>(sys_name,
251  unknown_name,
252  error_vals);
253  break;
254  }
255  case TYPE_VECTOR:
256  {
257  this->_compute_error<RealGradient>(sys_name,
258  unknown_name,
259  error_vals);
260  break;
261  }
262  default:
263  libmesh_error_msg("Invalid variable type!");
264  }
265 }
const EquationSystems & _equation_systems
const T_sys & get_system(const std::string &name) const
static FEFieldType field_type(const FEType &fe_type)
bool has_system(const std::string &name) const
std::vector< Real > & _check_inputs(const std::string &sys_name, const std::string &unknown_name)

◆ error_norm()

Real libMesh::ExactSolution::error_norm ( const std::string &  sys_name,
const std::string &  unknown_name,
const FEMNormType norm 
)
Returns
The error in the requested norm for the system sys_name for the unknown unknown_name.
Note
No error computations are actually performed, you must call compute_error() for that.
The result is not exact, but an approximation based on the chosen quadrature rule.

Definition at line 271 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, and libMesh::TYPE_SCALAR.

Referenced by hcurl_error(), and hdiv_error().

274 {
275  // Check the inputs for validity, and get a reference
276  // to the proper location to store the error
277  std::vector<Real> & error_vals = this->_check_inputs(sys_name,
278  unknown_name);
279 
280  libmesh_assert(_equation_systems.has_system(sys_name));
281  libmesh_assert(_equation_systems.get_system(sys_name).has_variable( unknown_name ));
282  const FEType & fe_type = _equation_systems.get_system(sys_name).variable_type(unknown_name);
283 
284  switch (norm)
285  {
286  case L2:
287  return std::sqrt(error_vals[0]);
288  case H1:
289  return std::sqrt(error_vals[0] + error_vals[1]);
290  case H2:
291  return std::sqrt(error_vals[0] + error_vals[1] + error_vals[2]);
292  case HCURL:
293  {
294  if (FEInterface::field_type(fe_type) == TYPE_SCALAR)
295  libmesh_error_msg("Cannot compute HCurl error norm of scalar-valued variables!");
296  else
297  return std::sqrt(error_vals[0] + error_vals[5]);
298  }
299  case HDIV:
300  {
301  if (FEInterface::field_type(fe_type) == TYPE_SCALAR)
302  libmesh_error_msg("Cannot compute HDiv error norm of scalar-valued variables!");
303  else
304  return std::sqrt(error_vals[0] + error_vals[6]);
305  }
306  case H1_SEMINORM:
307  return std::sqrt(error_vals[1]);
308  case H2_SEMINORM:
309  return std::sqrt(error_vals[2]);
310  case HCURL_SEMINORM:
311  {
312  if (FEInterface::field_type(fe_type) == TYPE_SCALAR)
313  libmesh_error_msg("Cannot compute HCurl error seminorm of scalar-valued variables!");
314  else
315  return std::sqrt(error_vals[5]);
316  }
317  case HDIV_SEMINORM:
318  {
319  if (FEInterface::field_type(fe_type) == TYPE_SCALAR)
320  libmesh_error_msg("Cannot compute HDiv error seminorm of scalar-valued variables!");
321  else
322  return std::sqrt(error_vals[6]);
323  }
324  case L1:
325  return error_vals[3];
326  case L_INF:
327  return error_vals[4];
328 
329  default:
330  libmesh_error_msg("Currently only Sobolev norms/seminorms are supported!");
331  }
332 }
const EquationSystems & _equation_systems
const T_sys & get_system(const std::string &name) const
static FEFieldType field_type(const FEType &fe_type)
bool has_system(const std::string &name) const
std::vector< Real > & _check_inputs(const std::string &sys_name, const std::string &unknown_name)

◆ extra_quadrature_order()

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

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

Definition at line 181 of file exact_solution.h.

References _extra_order.

182  { _extra_order = extraorder; }

◆ h1_error()

Real libMesh::ExactSolution::h1_error ( const std::string &  sys_name,
const std::string &  unknown_name 
)
Returns
The H1 error for the system sys_name for the unknown unknown_name.
Note
No error computations are actually performed, you must call compute_error() for that.

Definition at line 400 of file exact_solution.C.

References _check_inputs().

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

◆ h2_error()

Real libMesh::ExactSolution::h2_error ( const std::string &  sys_name,
const std::string &  unknown_name 
)
Returns
The H2 error for the system sys_name for the unknown unknown_name.
Note
No error computations are actually performed, you must call compute_error() for that.

Definition at line 432 of file exact_solution.C.

References _check_inputs().

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

◆ hcurl_error()

Real libMesh::ExactSolution::hcurl_error ( const std::string &  sys_name,
const std::string &  unknown_name 
)
Returns
The H(curl) error for the system sys_name for the unknown unknown_name.
Note
No error computations are actually performed, you must call compute_error() for that.
This is only valid for vector-valued elements. An error is thrown if requested for scalar-valued elements.

Definition at line 417 of file exact_solution.C.

References error_norm(), and libMesh::HCURL.

419 {
420  return this->error_norm(sys_name,unknown_name,HCURL);
421 }
Real error_norm(const std::string &sys_name, const std::string &unknown_name, const FEMNormType &norm)

◆ hdiv_error()

Real libMesh::ExactSolution::hdiv_error ( const std::string &  sys_name,
const std::string &  unknown_name 
)
Returns
The H(div) error for the system sys_name for the unknown unknown_name.
Note
No error computations are actually performed, you must call compute_error() for that.
This is only valid for vector-valued elements. An error is thrown if requested for scalar-valued elements.

Definition at line 424 of file exact_solution.C.

References error_norm(), and libMesh::HDIV.

426 {
427  return this->error_norm(sys_name,unknown_name,HDIV);
428 }
Real error_norm(const std::string &sys_name, const std::string &unknown_name, const FEMNormType &norm)

◆ l1_error()

Real libMesh::ExactSolution::l1_error ( const std::string &  sys_name,
const std::string &  unknown_name 
)
Returns
The integrated L1 error for the system sys_name for the unknown unknown_name.
Note
No error computations are actually performed, you must call compute_error() for that.

Definition at line 360 of file exact_solution.C.

References _check_inputs().

362 {
363 
364  // Check the inputs for validity, and get a reference
365  // to the proper location to store the error
366  std::vector<Real> & error_vals = this->_check_inputs(sys_name,
367  unknown_name);
368 
369  // Return the square root of the first component of the
370  // computed error.
371  return error_vals[3];
372 }
std::vector< Real > & _check_inputs(const std::string &sys_name, const std::string &unknown_name)

◆ l2_error()

Real libMesh::ExactSolution::l2_error ( const std::string &  sys_name,
const std::string &  unknown_name 
)
Returns
The integrated L2 error for the system sys_name for the unknown unknown_name.
Note
No error computations are actually performed, you must call compute_error() for that.

Definition at line 340 of file exact_solution.C.

References _check_inputs().

342 {
343 
344  // Check the inputs for validity, and get a reference
345  // to the proper location to store the error
346  std::vector<Real> & error_vals = this->_check_inputs(sys_name,
347  unknown_name);
348 
349  // Return the square root of the first component of the
350  // computed error.
351  return std::sqrt(error_vals[0]);
352 }
std::vector< Real > & _check_inputs(const std::string &sys_name, const std::string &unknown_name)

◆ l_inf_error()

Real libMesh::ExactSolution::l_inf_error ( const std::string &  sys_name,
const std::string &  unknown_name 
)
Returns
The L_INF error for the system sys_name for the unknown unknown_name.
Note
No error computations are actually performed, you must call compute_error() for that.
The result (as for the other norms as well) is not exact, but an approximation based on the chosen quadrature rule: to compute it, we take the max of the absolute value of the error over all the quadrature points.

Definition at line 380 of file exact_solution.C.

References _check_inputs().

382 {
383 
384  // Check the inputs for validity, and get a reference
385  // to the proper location to store the error
386  std::vector<Real> & error_vals = this->_check_inputs(sys_name,
387  unknown_name);
388 
389  // Return the square root of the first component of the
390  // computed error.
391  return error_vals[4];
392 }
std::vector< Real > & _check_inputs(const std::string &sys_name, const std::string &unknown_name)

◆ operator=() [1/2]

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

◆ operator=() [2/2]

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

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

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

◆ _errors

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

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

Definition at line 344 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 320 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 326 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 314 of file exact_solution.h.

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

◆ _extra_order

int libMesh::ExactSolution::_extra_order
private

Extra order to use for quadrature rule

Definition at line 361 of file exact_solution.h.

Referenced by _compute_error(), and extra_quadrature_order().


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