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

Member Typedef Documentation

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

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

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

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 105 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(), 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
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.

72 {
73 }

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 211 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().

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

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

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

161 {
162  if (_exact_derivs.size() <= sys_num)
163  _exact_derivs.resize(sys_num+1);
164 
165  if (g)
166  _exact_derivs[sys_num] = g->clone();
167 }
std::vector< std::unique_ptr< FunctionBase< Gradient > > > _exact_derivs
void libMesh::ExactSolution::attach_exact_deriv ( GradientFunctionPointer  gptr)

Definition at line 129 of file exact_solution.C.

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

130 {
131  libmesh_assert(gptr);
132 
133  // Clear out any previous _exact_derivs entries, then add a new
134  // entry for each system.
135  _exact_derivs.clear();
136 
137  for (unsigned int sys=0; sys<_equation_systems.n_systems(); ++sys)
138  {
139  const System & system = _equation_systems.get_system(sys);
140  _exact_derivs.emplace_back(libmesh_make_unique<WrappedFunction<Gradient>>(system, gptr, &_equation_systems.parameters));
141  }
142 
143  // If we're using exact values, we're not using a fine grid solution
145 }
const EquationSystems & _equation_systems
const class libmesh_nullptr_t libmesh_nullptr
const EquationSystems * _equation_systems_fine
std::vector< std::unique_ptr< FunctionBase< Gradient > > > _exact_derivs
unsigned int n_systems() const
const T_sys & get_system(const std::string &name) const
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 148 of file exact_solution.C.

References _exact_derivs, and libmesh_nullptr.

149 {
150  // Automatically delete any previous _exact_derivs entries, then add a new
151  // entry for each system.
152  _exact_derivs.clear();
153 
154  for (auto ptr : g)
155  _exact_derivs.emplace_back(ptr ? ptr->clone() : libmesh_nullptr);
156 }
const class libmesh_nullptr_t libmesh_nullptr
std::vector< std::unique_ptr< 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 200 of file exact_solution.C.

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

202 {
203  if (_exact_hessians.size() <= sys_num)
204  _exact_hessians.resize(sys_num+1);
205 
206  if (h)
207  _exact_hessians[sys_num] = h->clone();
208 }
std::vector< std::unique_ptr< FunctionBase< Tensor > > > _exact_hessians
void libMesh::ExactSolution::attach_exact_hessian ( HessianFunctionPointer  hptr)

Definition at line 170 of file exact_solution.C.

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

171 {
172  libmesh_assert(hptr);
173 
174  // Clear out any previous _exact_hessians entries, then add a new
175  // entry for each system.
176  _exact_hessians.clear();
177 
178  for (unsigned int sys=0; sys<_equation_systems.n_systems(); ++sys)
179  {
180  const System & system = _equation_systems.get_system(sys);
181  _exact_hessians.emplace_back(libmesh_make_unique<WrappedFunction<Tensor>>(system, hptr, &_equation_systems.parameters));
182  }
183 
184  // If we're using exact values, we're not using a fine grid solution
186 }
const EquationSystems & _equation_systems
const class libmesh_nullptr_t libmesh_nullptr
const EquationSystems * _equation_systems_fine
unsigned int n_systems() const
std::vector< std::unique_ptr< FunctionBase< Tensor > > > _exact_hessians
const T_sys & get_system(const std::string &name) const
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 189 of file exact_solution.C.

References _exact_hessians, and libmesh_nullptr.

190 {
191  // Automatically delete any previous _exact_hessians entries, then add a new
192  // entry for each system.
193  _exact_hessians.clear();
194 
195  for (auto ptr : h)
196  _exact_hessians.emplace_back(ptr ? ptr->clone() : libmesh_nullptr);
197 }
const class libmesh_nullptr_t libmesh_nullptr
std::vector< std::unique_ptr< FunctionBase< Tensor > > > _exact_hessians
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 118 of file exact_solution.C.

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

120 {
121  if (_exact_values.size() <= sys_num)
122  _exact_values.resize(sys_num+1);
123 
124  if (f)
125  _exact_values[sys_num] = f->clone();
126 }
virtual std::unique_ptr< FunctionBase< Output > > clone() const =0
std::vector< std::unique_ptr< FunctionBase< Number > > > _exact_values
void libMesh::ExactSolution::attach_exact_value ( ValueFunctionPointer  fptr)

Definition at line 88 of file exact_solution.C.

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

89 {
90  libmesh_assert(fptr);
91 
92  // Clear out any previous _exact_values entries, then add a new
93  // entry for each system.
94  _exact_values.clear();
95 
96  for (unsigned int sys=0; sys<_equation_systems.n_systems(); ++sys)
97  {
98  const System & system = _equation_systems.get_system(sys);
99  _exact_values.emplace_back(libmesh_make_unique<WrappedFunction<Number>>(system, fptr, &_equation_systems.parameters));
100  }
101 
102  // If we're using exact values, we're not using a fine grid solution
104 }
const EquationSystems & _equation_systems
const class libmesh_nullptr_t libmesh_nullptr
const EquationSystems * _equation_systems_fine
unsigned int n_systems() const
const T_sys & get_system(const std::string &name) const
std::vector< std::unique_ptr< FunctionBase< Number > > > _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 107 of file exact_solution.C.

References _exact_values, and libmesh_nullptr.

108 {
109  // Automatically delete any previous _exact_values entries, then add a new
110  // entry for each system.
111  _exact_values.clear();
112 
113  for (auto ptr : f)
114  _exact_values.emplace_back(ptr ? ptr->clone() : libmesh_nullptr);
115 }
const class libmesh_nullptr_t libmesh_nullptr
std::vector< std::unique_ptr< FunctionBase< Number > > > _exact_values
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 76 of file exact_solution.C.

References _equation_systems_fine, _exact_derivs, _exact_hessians, and _exact_values.

77 {
78  libmesh_assert(es_fine);
79  _equation_systems_fine = es_fine;
80 
81  // If we're using a fine grid solution, we're not using exact values
82  _exact_values.clear();
83  _exact_derivs.clear();
84  _exact_hessians.clear();
85 }
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
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 237 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.

Referenced by extra_quadrature_order().

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

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

162  { _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 403 of file exact_solution.C.

References _check_inputs().

Referenced by extra_quadrature_order().

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

References _check_inputs().

Referenced by extra_quadrature_order().

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

References error_norm(), and libMesh::HCURL.

Referenced by extra_quadrature_order().

422 {
423  return this->error_norm(sys_name,unknown_name,HCURL);
424 }
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 427 of file exact_solution.C.

References error_norm(), and libMesh::HDIV.

Referenced by extra_quadrature_order().

429 {
430  return this->error_norm(sys_name,unknown_name,HDIV);
431 }
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 363 of file exact_solution.C.

References _check_inputs().

Referenced by extra_quadrature_order().

365 {
366 
367  // Check the inputs for validity, and get a reference
368  // to the proper location to store the error
369  std::vector<Real> & error_vals = this->_check_inputs(sys_name,
370  unknown_name);
371 
372  // Return the square root of the first component of the
373  // computed error.
374  return error_vals[3];
375 }
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 343 of file exact_solution.C.

References _check_inputs().

Referenced by extra_quadrature_order().

345 {
346 
347  // Check the inputs for validity, and get a reference
348  // to the proper location to store the error
349  std::vector<Real> & error_vals = this->_check_inputs(sys_name,
350  unknown_name);
351 
352  // Return the square root of the first component of the
353  // computed error.
354  return std::sqrt(error_vals[0]);
355 }
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 383 of file exact_solution.C.

References _check_inputs().

Referenced by extra_quadrature_order().

385 {
386 
387  // Check the inputs for validity, and get a reference
388  // to the proper location to store the error
389  std::vector<Real> & error_vals = this->_check_inputs(sys_name,
390  unknown_name);
391 
392  // Return the square root of the first component of the
393  // computed error.
394  return error_vals[4];
395 }
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 330 of file exact_solution.h.

Referenced by _compute_error(), attach_exact_deriv(), attach_exact_hessian(), attach_exact_value(), 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 336 of file exact_solution.h.

Referenced by _compute_error(), attach_exact_deriv(), attach_exact_hessian(), 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 324 of file exact_solution.h.

Referenced by _check_inputs(), and ExactSolution().

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

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

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

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

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

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

int libMesh::ExactSolution::_extra_order
private

Extra order to use for quadrature rule

Definition at line 341 of file exact_solution.h.

Referenced by _compute_error(), and extra_quadrature_order().


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