libMesh
exact_solution.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 // C++ includes
19 
20 
21 // Local includes
22 #include "libmesh/dof_map.h"
23 #include "libmesh/elem.h"
24 #include "libmesh/exact_solution.h"
25 #include "libmesh/equation_systems.h"
26 #include "libmesh/fe_base.h"
27 #include "libmesh/function_base.h"
28 #include "libmesh/mesh_base.h"
29 #include "libmesh/mesh_function.h"
30 #include "libmesh/numeric_vector.h"
31 #include "libmesh/parallel.h"
32 #include "libmesh/quadrature.h"
33 #include "libmesh/wrapped_function.h"
34 #include "libmesh/fe_interface.h"
35 #include "libmesh/raw_accessor.h"
36 #include "libmesh/tensor_tools.h"
37 
38 namespace libMesh
39 {
40 
42  _equation_systems(es),
43  _equation_systems_fine(libmesh_nullptr),
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 }
69 
70 
72 {
73  // delete will clean up any cloned functors and no-op on any NULL
74  // pointers
75 
76  for (std::size_t i=0; i != _exact_values.size(); ++i)
77  delete (_exact_values[i]);
78 
79  for (std::size_t i=0; i != _exact_derivs.size(); ++i)
80  delete (_exact_derivs[i]);
81 
82  for (std::size_t i=0; i != _exact_hessians.size(); ++i)
83  delete (_exact_hessians[i]);
84 }
85 
86 
88 {
89  libmesh_assert(es_fine);
90  _equation_systems_fine = es_fine;
91 
92  // If we're using a fine grid solution, we're not using exact values
93  _exact_values.clear();
94  _exact_derivs.clear();
95  _exact_hessians.clear();
96 }
97 
98 
100  const Parameters & parameters,
101  const std::string & sys_name,
102  const std::string & unknown_name))
103 {
105 
106  // Clear out any previous _exact_values entries, then add a new
107  // entry for each system.
108  _exact_values.clear();
109 
110  for (unsigned int sys=0; sys<_equation_systems.n_systems(); ++sys)
111  {
112  const System & system = _equation_systems.get_system(sys);
113  _exact_values.push_back
115  (system, fptr, &_equation_systems.parameters));
116  }
117 
118  // If we're using exact values, we're not using a fine grid solution
120 }
121 
122 
124 {
125  // Clear out any previous _exact_values entries, then add a new
126  // entry for each system.
127  for (std::size_t i=0; i != _exact_values.size(); ++i)
128  delete (_exact_values[i]);
129 
130  _exact_values.clear();
131  _exact_values.resize(f.size(), libmesh_nullptr);
132 
133  // We use clone() to get non-sliced copies of FunctionBase
134  // subclasses, but we don't currently put the resulting UniquePtrs
135  // into an STL container.
136  for (std::size_t i=0; i != f.size(); ++i)
137  if (f[i])
138  _exact_values[i] = f[i]->clone().release();
139 }
140 
141 
142 void ExactSolution::attach_exact_value (unsigned int sys_num,
144 {
145  if (_exact_values.size() <= sys_num)
146  _exact_values.resize(sys_num+1, libmesh_nullptr);
147 
148  if (f)
149  _exact_values[sys_num] = f->clone().release();
150 }
151 
152 
154  const Parameters & parameters,
155  const std::string & sys_name,
156  const std::string & unknown_name))
157 {
159 
160  // Clear out any previous _exact_derivs entries, then add a new
161  // entry for each system.
162  _exact_derivs.clear();
163 
164  for (unsigned int sys=0; sys<_equation_systems.n_systems(); ++sys)
165  {
166  const System & system = _equation_systems.get_system(sys);
167  _exact_derivs.push_back
169  (system, gptr, &_equation_systems.parameters));
170  }
171 
172  // If we're using exact values, we're not using a fine grid solution
174 }
175 
176 
178 {
179  // Clear out any previous _exact_derivs entries, then add a new
180  // entry for each system.
181  for (std::size_t i=0; i != _exact_derivs.size(); ++i)
182  delete (_exact_derivs[i]);
183 
184  _exact_derivs.clear();
185  _exact_derivs.resize(g.size(), libmesh_nullptr);
186 
187  // We use clone() to get non-sliced copies of FunctionBase
188  // subclasses, but we don't currently put the resulting UniquePtrs
189  // into an STL container.
190  for (std::size_t i=0; i != g.size(); ++i)
191  if (g[i])
192  _exact_derivs[i] = g[i]->clone().release();
193 }
194 
195 
196 void ExactSolution::attach_exact_deriv (unsigned int sys_num,
198 {
199  if (_exact_derivs.size() <= sys_num)
200  _exact_derivs.resize(sys_num+1, libmesh_nullptr);
201 
202  if (g)
203  _exact_derivs[sys_num] = g->clone().release();
204 }
205 
206 
208  const Parameters & parameters,
209  const std::string & sys_name,
210  const std::string & unknown_name))
211 {
212  libmesh_assert(hptr);
213 
214  // Clear out any previous _exact_hessians entries, then add a new
215  // entry for each system.
216  _exact_hessians.clear();
217 
218  for (unsigned int sys=0; sys<_equation_systems.n_systems(); ++sys)
219  {
220  const System & system = _equation_systems.get_system(sys);
221  _exact_hessians.push_back
223  (system, hptr, &_equation_systems.parameters));
224  }
225 
226  // If we're using exact values, we're not using a fine grid solution
228 }
229 
230 
232 {
233  // Clear out any previous _exact_hessians entries, then add a new
234  // entry for each system.
235  for (std::size_t i=0; i != _exact_hessians.size(); ++i)
236  delete (_exact_hessians[i]);
237 
238  _exact_hessians.clear();
239  _exact_hessians.resize(h.size(), libmesh_nullptr);
240 
241  // We use clone() to get non-sliced copies of FunctionBase
242  // subclasses, but we don't currently put the resulting UniquePtrs
243  // into an STL container.
244  for (std::size_t i=0; i != h.size(); ++i)
245  if (h[i])
246  _exact_hessians[i] = h[i]->clone().release();
247 }
248 
249 
250 void ExactSolution::attach_exact_hessian (unsigned int sys_num,
252 {
253  if (_exact_hessians.size() <= sys_num)
254  _exact_hessians.resize(sys_num+1, libmesh_nullptr);
255 
256  if (h)
257  _exact_hessians[sys_num] = h->clone().release();
258 }
259 
260 
261 std::vector<Real> & ExactSolution::_check_inputs(const std::string & sys_name,
262  const std::string & unknown_name)
263 {
264  // If no exact solution function or fine grid solution has been
265  // attached, we now just compute the solution norm (i.e. the
266  // difference from an "exact solution" of zero
267 
268  // Make sure the requested sys_name exists.
269  std::map<std::string, SystemErrorMap>::iterator sys_iter =
270  _errors.find(sys_name);
271 
272  if (sys_iter == _errors.end())
273  libmesh_error_msg("Sorry, couldn't find the requested system '" << sys_name << "'.");
274 
275  // Make sure the requested unknown_name exists.
276  SystemErrorMap::iterator var_iter = (*sys_iter).second.find(unknown_name);
277 
278  if (var_iter == (*sys_iter).second.end())
279  libmesh_error_msg("Sorry, couldn't find the requested variable '" << unknown_name << "'.");
280 
281  // Return a reference to the proper error entry
282  return (*var_iter).second;
283 }
284 
285 
286 
287 void ExactSolution::compute_error(const std::string & sys_name,
288  const std::string & unknown_name)
289 {
290  // Check the inputs for validity, and get a reference
291  // to the proper location to store the error
292  std::vector<Real> & error_vals = this->_check_inputs(sys_name,
293  unknown_name);
294 
296  const System & sys = _equation_systems.get_system<System>( sys_name );
297 
298  libmesh_assert( sys.has_variable( unknown_name ) );
299  switch( FEInterface::field_type(sys.variable_type( unknown_name )) )
300  {
301  case TYPE_SCALAR:
302  {
303  this->_compute_error<Real>(sys_name,
304  unknown_name,
305  error_vals);
306  break;
307  }
308  case TYPE_VECTOR:
309  {
310  this->_compute_error<RealGradient>(sys_name,
311  unknown_name,
312  error_vals);
313  break;
314  }
315  default:
316  libmesh_error_msg("Invalid variable type!");
317  }
318 }
319 
320 
321 
322 
323 
324 Real ExactSolution::error_norm(const std::string & sys_name,
325  const std::string & unknown_name,
326  const FEMNormType & norm)
327 {
328  // Check the inputs for validity, and get a reference
329  // to the proper location to store the error
330  std::vector<Real> & error_vals = this->_check_inputs(sys_name,
331  unknown_name);
332 
334  libmesh_assert(_equation_systems.get_system(sys_name).has_variable( unknown_name ));
335  const FEType & fe_type = _equation_systems.get_system(sys_name).variable_type(unknown_name);
336 
337  switch (norm)
338  {
339  case L2:
340  return std::sqrt(error_vals[0]);
341  case H1:
342  return std::sqrt(error_vals[0] + error_vals[1]);
343  case H2:
344  return std::sqrt(error_vals[0] + error_vals[1] + error_vals[2]);
345  case HCURL:
346  {
347  if (FEInterface::field_type(fe_type) == TYPE_SCALAR)
348  libmesh_error_msg("Cannot compute HCurl error norm of scalar-valued variables!");
349  else
350  return std::sqrt(error_vals[0] + error_vals[5]);
351  }
352  case HDIV:
353  {
354  if (FEInterface::field_type(fe_type) == TYPE_SCALAR)
355  libmesh_error_msg("Cannot compute HDiv error norm of scalar-valued variables!");
356  else
357  return std::sqrt(error_vals[0] + error_vals[6]);
358  }
359  case H1_SEMINORM:
360  return std::sqrt(error_vals[1]);
361  case H2_SEMINORM:
362  return std::sqrt(error_vals[2]);
363  case HCURL_SEMINORM:
364  {
365  if (FEInterface::field_type(fe_type) == TYPE_SCALAR)
366  libmesh_error_msg("Cannot compute HCurl error seminorm of scalar-valued variables!");
367  else
368  return std::sqrt(error_vals[5]);
369  }
370  case HDIV_SEMINORM:
371  {
372  if (FEInterface::field_type(fe_type) == TYPE_SCALAR)
373  libmesh_error_msg("Cannot compute HDiv error seminorm of scalar-valued variables!");
374  else
375  return std::sqrt(error_vals[6]);
376  }
377  case L1:
378  return error_vals[3];
379  case L_INF:
380  return error_vals[4];
381 
382  default:
383  libmesh_error_msg("Currently only Sobolev norms/seminorms are supported!");
384  }
385 }
386 
387 
388 
389 
390 
391 
392 
393 Real ExactSolution::l2_error(const std::string & sys_name,
394  const std::string & unknown_name)
395 {
396 
397  // Check the inputs for validity, and get a reference
398  // to the proper location to store the error
399  std::vector<Real> & error_vals = this->_check_inputs(sys_name,
400  unknown_name);
401 
402  // Return the square root of the first component of the
403  // computed error.
404  return std::sqrt(error_vals[0]);
405 }
406 
407 
408 
409 
410 
411 
412 
413 Real ExactSolution::l1_error(const std::string & sys_name,
414  const std::string & unknown_name)
415 {
416 
417  // Check the inputs for validity, and get a reference
418  // to the proper location to store the error
419  std::vector<Real> & error_vals = this->_check_inputs(sys_name,
420  unknown_name);
421 
422  // Return the square root of the first component of the
423  // computed error.
424  return error_vals[3];
425 }
426 
427 
428 
429 
430 
431 
432 
433 Real ExactSolution::l_inf_error(const std::string & sys_name,
434  const std::string & unknown_name)
435 {
436 
437  // Check the inputs for validity, and get a reference
438  // to the proper location to store the error
439  std::vector<Real> & error_vals = this->_check_inputs(sys_name,
440  unknown_name);
441 
442  // Return the square root of the first component of the
443  // computed error.
444  return error_vals[4];
445 }
446 
447 
448 
449 
450 
451 
452 
453 Real ExactSolution::h1_error(const std::string & sys_name,
454  const std::string & unknown_name)
455 {
456  // If the user has supplied no exact derivative function, we
457  // just integrate the H1 norm of the solution; i.e. its
458  // difference from an "exact solution" of zero.
459 
460  // Check the inputs for validity, and get a reference
461  // to the proper location to store the error
462  std::vector<Real> & error_vals = this->_check_inputs(sys_name,
463  unknown_name);
464 
465  // Return the square root of the sum of the computed errors.
466  return std::sqrt(error_vals[0] + error_vals[1]);
467 }
468 
469 
470 Real ExactSolution::hcurl_error(const std::string & sys_name,
471  const std::string & unknown_name)
472 {
473  return this->error_norm(sys_name,unknown_name,HCURL);
474 }
475 
476 
477 Real ExactSolution::hdiv_error(const std::string & sys_name,
478  const std::string & unknown_name)
479 {
480  return this->error_norm(sys_name,unknown_name,HDIV);
481 }
482 
483 
484 
485 Real ExactSolution::h2_error(const std::string & sys_name,
486  const std::string & unknown_name)
487 {
488  // If the user has supplied no exact derivative functions, we
489  // just integrate the H2 norm of the solution; i.e. its
490  // difference from an "exact solution" of zero.
491 
492  // Check the inputs for validity, and get a reference
493  // to the proper location to store the error
494  std::vector<Real> & error_vals = this->_check_inputs(sys_name,
495  unknown_name);
496 
497  // Return the square root of the sum of the computed errors.
498  return std::sqrt(error_vals[0] + error_vals[1] + error_vals[2]);
499 }
500 
501 
502 
503 
504 
505 
506 
507 template<typename OutputShape>
508 void ExactSolution::_compute_error(const std::string & sys_name,
509  const std::string & unknown_name,
510  std::vector<Real> & error_vals)
511 {
512  // Make sure we aren't "overconfigured"
514 
515  // We need a communicator.
517 
518  // This function must be run on all processors at once
519  libmesh_parallel_only(communicator);
520 
521  // Get a reference to the system whose error is being computed.
522  // If we have a fine grid, however, we'll integrate on that instead
523  // for more accuracy.
524  const System & computed_system = _equation_systems_fine ?
526  _equation_systems.get_system (sys_name);
527 
528  const Real time = _equation_systems.get_system(sys_name).time;
529 
530  const unsigned int sys_num = computed_system.number();
531  const unsigned int var = computed_system.variable_number(unknown_name);
532  const unsigned int var_component =
533  computed_system.variable_scalar_number(var, 0);
534 
535  // Prepare a global solution and a MeshFunction of the coarse system if we need one
536  UniquePtr<MeshFunction> coarse_values;
539  {
540  const System & comparison_system
541  = _equation_systems.get_system(sys_name);
542 
543  std::vector<Number> global_soln;
544  comparison_system.update_global_solution(global_soln);
545  comparison_soln->init(comparison_system.solution->size(), true, SERIAL);
546  (*comparison_soln) = global_soln;
547 
548  coarse_values = UniquePtr<MeshFunction>
550  *comparison_soln,
551  comparison_system.get_dof_map(),
552  comparison_system.variable_number(unknown_name)));
553  coarse_values->init();
554  }
555 
556  // Initialize any functors we're going to use
557  for (std::size_t i=0; i != _exact_values.size(); ++i)
558  if (_exact_values[i])
559  _exact_values[i]->init();
560 
561  for (std::size_t i=0; i != _exact_derivs.size(); ++i)
562  if (_exact_derivs[i])
563  _exact_derivs[i]->init();
564 
565  for (std::size_t i=0; i != _exact_hessians.size(); ++i)
566  if (_exact_hessians[i])
567  _exact_hessians[i]->init();
568 
569  // Get a reference to the dofmap and mesh for that system
570  const DofMap & computed_dof_map = computed_system.get_dof_map();
571 
572  const MeshBase & _mesh = computed_system.get_mesh();
573 
574  // Grab which element dimensions are present in the mesh
575  const std::set<unsigned char> & elem_dims = _mesh.elem_dimensions();
576 
577  // Zero the error before summation
578  // 0 - sum of square of function error (L2)
579  // 1 - sum of square of gradient error (H1 semi)
580  // 2 - sum of square of Hessian error (H2 semi)
581  // 3 - sum of sqrt(square of function error) (L1)
582  // 4 - max of sqrt(square of function error) (Linfty)
583  // 5 - sum of square of curl error (HCurl semi)
584  // 6 - sum of square of div error (HDiv semi)
585  error_vals = std::vector<Real>(7, 0.);
586 
587  // Construct Quadrature rule based on default quadrature order
588  const FEType & fe_type = computed_dof_map.variable_type(var);
589 
590  unsigned int n_vec_dim = FEInterface::n_vec_dim( _mesh, fe_type );
591 
592  // FIXME: MeshFunction needs to be updated to support vector-valued
593  // elements before we can use a reference solution.
594  if ((n_vec_dim > 1) && _equation_systems_fine)
595  {
596  libMesh::err << "Error calculation using reference solution not yet\n"
597  << "supported for vector-valued elements."
598  << std::endl;
599  libmesh_not_implemented();
600  }
601 
602 
603  // Allow space for dims 0-3, even if we don't use them all
604  std::vector<FEGenericBase<OutputShape> *> fe_ptrs(4, libmesh_nullptr);
605  std::vector<QBase *> q_rules(4, libmesh_nullptr);
606 
607  // Prepare finite elements for each dimension present in the mesh
608  for (std::set<unsigned char>::const_iterator d_it = elem_dims.begin();
609  d_it != elem_dims.end(); ++d_it)
610  {
611  q_rules[*d_it] =
612  fe_type.default_quadrature_rule (*d_it, _extra_order).release();
613 
614  // Construct finite element object
615 
616  fe_ptrs[*d_it] = FEGenericBase<OutputShape>::build(*d_it, fe_type).release();
617 
618  // Attach quadrature rule to FE object
619  fe_ptrs[*d_it]->attach_quadrature_rule (q_rules[*d_it]);
620  }
621 
622  // The global degree of freedom indices associated
623  // with the local degrees of freedom.
624  std::vector<dof_id_type> dof_indices;
625 
626 
627  //
628  // Begin the loop over the elements
629  //
630  // TODO: this ought to be threaded (and using subordinate
631  // MeshFunction objects in each thread rather than a single
632  // master)
633  for (const auto & elem : _mesh.active_local_element_ptr_range())
634  {
635  // Store a pointer to the element we are currently
636  // working on. This allows for nicer syntax later.
637  const unsigned int dim = elem->dim();
638 
639  const subdomain_id_type elem_subid = elem->subdomain_id();
640 
641  // If the variable is not active on this subdomain, don't bother
642  if (!computed_system.variable(var).active_on_subdomain(elem_subid))
643  continue;
644 
645  /* If the variable is active, then we're going to restrict the
646  MeshFunction evaluations to the current element subdomain.
647  This is for cases such as mixed dimension meshes where we want
648  to restrict the calculation to one particular domain. */
649  std::set<subdomain_id_type> subdomain_id;
650  subdomain_id.insert(elem_subid);
651 
652  FEGenericBase<OutputShape> * fe = fe_ptrs[dim];
653  QBase * qrule = q_rules[dim];
654  libmesh_assert(fe);
655  libmesh_assert(qrule);
656 
657  // The Jacobian*weight at the quadrature points.
658  const std::vector<Real> & JxW = fe->get_JxW();
659 
660  // The value of the shape functions at the quadrature points
661  // i.e. phi(i) = phi_values[i][qp]
662  const std::vector<std::vector<OutputShape>> & phi_values = fe->get_phi();
663 
664  // The value of the shape function gradients at the quadrature points
665  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputGradient>> &
666  dphi_values = fe->get_dphi();
667 
668  // The value of the shape function curls at the quadrature points
669  // Only computed for vector-valued elements
670  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputShape>> * curl_values = libmesh_nullptr;
671 
672  // The value of the shape function divergences at the quadrature points
673  // Only computed for vector-valued elements
674  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputDivergence>> * div_values = libmesh_nullptr;
675 
676  if (FEInterface::field_type(fe_type) == TYPE_VECTOR)
677  {
678  curl_values = &fe->get_curl_phi();
679  div_values = &fe->get_div_phi();
680  }
681 
682 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
683  // The value of the shape function second derivatives at the quadrature points
684  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputTensor>> &
685  d2phi_values = fe->get_d2phi();
686 #endif
687 
688  // The XYZ locations (in physical space) of the quadrature points
689  const std::vector<Point> & q_point = fe->get_xyz();
690 
691  // reinitialize the element-specific data
692  // for the current element
693  fe->reinit (elem);
694 
695  // Get the local to global degree of freedom maps
696  computed_dof_map.dof_indices (elem, dof_indices, var);
697 
698  // The number of quadrature points
699  const unsigned int n_qp = qrule->n_points();
700 
701  // The number of shape functions
702  const unsigned int n_sf =
703  cast_int<unsigned int>(dof_indices.size());
704 
705  //
706  // Begin the loop over the Quadrature points.
707  //
708  for (unsigned int qp=0; qp<n_qp; qp++)
709  {
710  // Real u_h = 0.;
711  // RealGradient grad_u_h;
712 
714 
716 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
718 #endif
719  typename FEGenericBase<OutputShape>::OutputNumber curl_u_h(0.0);
721 
722  // Compute solution values at the current
723  // quadrature point. This requires a sum
724  // over all the shape functions evaluated
725  // at the quadrature point.
726  for (unsigned int i=0; i<n_sf; i++)
727  {
728  // Values from current solution.
729  u_h += phi_values[i][qp]*computed_system.current_solution (dof_indices[i]);
730  grad_u_h += dphi_values[i][qp]*computed_system.current_solution (dof_indices[i]);
731 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
732  grad2_u_h += d2phi_values[i][qp]*computed_system.current_solution (dof_indices[i]);
733 #endif
734  if (FEInterface::field_type(fe_type) == TYPE_VECTOR)
735  {
736  curl_u_h += (*curl_values)[i][qp]*computed_system.current_solution (dof_indices[i]);
737  div_u_h += (*div_values)[i][qp]*computed_system.current_solution (dof_indices[i]);
738  }
739  }
740 
741  // Compute the value of the error at this quadrature point
742  typename FEGenericBase<OutputShape>::OutputNumber exact_val(0);
743  RawAccessor<typename FEGenericBase<OutputShape>::OutputNumber> exact_val_accessor( exact_val, dim );
744  if (_exact_values.size() > sys_num && _exact_values[sys_num])
745  {
746  for (unsigned int c = 0; c < n_vec_dim; c++)
747  exact_val_accessor(c) =
748  _exact_values[sys_num]->
749  component(var_component+c, q_point[qp], time);
750  }
751  else if (_equation_systems_fine)
752  {
753  // FIXME: Needs to be updated for vector-valued elements
754  DenseVector<Number> output(1);
755  (*coarse_values)(q_point[qp],time,output,&subdomain_id);
756  exact_val = output(0);
757  }
758  const typename FEGenericBase<OutputShape>::OutputNumber val_error = u_h - exact_val;
759 
760  // Add the squares of the error to each contribution
761  Real error_sq = TensorTools::norm_sq(val_error);
762  error_vals[0] += JxW[qp]*error_sq;
763 
764  Real norm = sqrt(error_sq);
765  error_vals[3] += JxW[qp]*norm;
766 
767  if (error_vals[4]<norm) { error_vals[4] = norm; }
768 
769  // Compute the value of the error in the gradient at this
770  // quadrature point
772  RawAccessor<typename FEGenericBase<OutputShape>::OutputNumberGradient> exact_grad_accessor( exact_grad, LIBMESH_DIM );
773  if (_exact_derivs.size() > sys_num && _exact_derivs[sys_num])
774  {
775  for (unsigned int c = 0; c < n_vec_dim; c++)
776  for (unsigned int d = 0; d < LIBMESH_DIM; d++)
777  exact_grad_accessor(d + c*LIBMESH_DIM) =
778  _exact_derivs[sys_num]->
779  component(var_component+c, q_point[qp], time)(d);
780  }
781  else if (_equation_systems_fine)
782  {
783  // FIXME: Needs to be updated for vector-valued elements
784  std::vector<Gradient> output(1);
785  coarse_values->gradient(q_point[qp],time,output,&subdomain_id);
786  exact_grad = output[0];
787  }
788 
789  const typename FEGenericBase<OutputShape>::OutputNumberGradient grad_error = grad_u_h - exact_grad;
790 
791  error_vals[1] += JxW[qp]*grad_error.norm_sq();
792 
793 
794  if (FEInterface::field_type(fe_type) == TYPE_VECTOR)
795  {
796  // Compute the value of the error in the curl at this
797  // quadrature point
798  typename FEGenericBase<OutputShape>::OutputNumber exact_curl(0.0);
799  if (_exact_derivs.size() > sys_num && _exact_derivs[sys_num])
800  {
801  exact_curl = TensorTools::curl_from_grad( exact_grad );
802  }
803  else if (_equation_systems_fine)
804  {
805  // FIXME: Need to implement curl for MeshFunction and support reference
806  // solution for vector-valued elements
807  }
808 
809  const typename FEGenericBase<OutputShape>::OutputNumber curl_error = curl_u_h - exact_curl;
810 
811  error_vals[5] += JxW[qp]*TensorTools::norm_sq(curl_error);
812 
813  // Compute the value of the error in the divergence at this
814  // quadrature point
816  if (_exact_derivs.size() > sys_num && _exact_derivs[sys_num])
817  {
818  exact_div = TensorTools::div_from_grad( exact_grad );
819  }
820  else if (_equation_systems_fine)
821  {
822  // FIXME: Need to implement div for MeshFunction and support reference
823  // solution for vector-valued elements
824  }
825 
826  const typename FEGenericBase<OutputShape>::OutputNumberDivergence div_error = div_u_h - exact_div;
827 
828  error_vals[6] += JxW[qp]*TensorTools::norm_sq(div_error);
829  }
830 
831 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
832  // Compute the value of the error in the hessian at this
833  // quadrature point
835  RawAccessor<typename FEGenericBase<OutputShape>::OutputNumberTensor> exact_hess_accessor( exact_hess, dim );
836  if (_exact_hessians.size() > sys_num && _exact_hessians[sys_num])
837  {
838  //FIXME: This needs to be implemented to support rank 3 tensors
839  // which can't happen until type_n_tensor is fully implemented
840  // and a RawAccessor<TypeNTensor> is fully implemented
841  if (FEInterface::field_type(fe_type) == TYPE_VECTOR)
842  libmesh_not_implemented();
843 
844  for (unsigned int c = 0; c < n_vec_dim; c++)
845  for (unsigned int d = 0; d < dim; d++)
846  for (unsigned int e =0; e < dim; e++)
847  exact_hess_accessor(d + e*dim + c*dim*dim) =
848  _exact_hessians[sys_num]->
849  component(var_component+c, q_point[qp], time)(d,e);
850  }
851  else if (_equation_systems_fine)
852  {
853  // FIXME: Needs to be updated for vector-valued elements
854  std::vector<Tensor> output(1);
855  coarse_values->hessian(q_point[qp],time,output,&subdomain_id);
856  exact_hess = output[0];
857  }
858 
859  const typename FEGenericBase<OutputShape>::OutputNumberTensor grad2_error = grad2_u_h - exact_hess;
860 
861  // FIXME: PB: Is this what we want for rank 3 tensors?
862  error_vals[2] += JxW[qp]*grad2_error.norm_sq();
863 #endif
864 
865  } // end qp loop
866  } // end element loop
867 
868  // Clean up the FE and QBase pointers we created
869  for (std::set<unsigned char>::const_iterator d_it = elem_dims.begin();
870  d_it != elem_dims.end(); ++d_it)
871  {
872  delete fe_ptrs[*d_it];
873  delete q_rules[*d_it];
874  }
875 
876  // Add up the error values on all processors, except for the L-infty
877  // norm, for which the maximum is computed.
878  Real l_infty_norm = error_vals[4];
879  communicator.max(l_infty_norm);
880  communicator.sum(error_vals);
881  error_vals[4] = l_infty_norm;
882 }
883 
884 // Explicit instantiations of templated member functions
885 template void ExactSolution::_compute_error<Real>(const std::string &, const std::string &, std::vector<Real> &);
886 template void ExactSolution::_compute_error<RealGradient>(const std::string &, const std::string &, std::vector<Real> &);
887 
888 } // namespace libMesh
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
OStreamProxy err
std::vector< FunctionBase< Tensor > * > _exact_hessians
User-provided functors which compute the exact hessians of the solution for each system.
This class provides a wrapper with which to evaluate a (libMesh-style) function pointer in a Function...
This is the EquationSystems class.
void _compute_error(const std::string &sys_name, const std::string &unknown_name, std::vector< Real > &error_vals)
This function computes the error (in the solution and its first derivative) for a single unknown in a...
const std::vector< std::vector< OutputDivergence > > & get_div_phi() const
Definition: fe_base.h:232
Encapsulates the MPI_Comm object.
Definition: parallel.h:657
bool has_system(const std::string &name) const
TensorTools::DecrementRank< OutputNumber >::type OutputNumberDivergence
Definition: fe_base.h:127
Number div_from_grad(const VectorValue< Number > &grad)
Dummy. Divergence of a scalar not defined, but is needed for ExactSolution to compile.
Definition: tensor_tools.C:54
const std::string & variable_name(const unsigned int i) const
Definition: system.h:2134
This class provides the ability to map between arbitrary, user-defined strings and several data types...
Definition: parameters.h:63
std::vector< FunctionBase< Number > * > _exact_values
User-provided functors which compute the exact value of the solution for each system.
const EquationSystems & _equation_systems
Constant reference to the EquationSystems object used for the simulation.
const FEType & variable_type(const unsigned int c) const
Definition: dof_map.h:1697
ExactSolution(const EquationSystems &es)
Constructor.
unsigned int dim
ImplicitSystem & sys
MPI_Comm communicator
Communicator object for talking with subsets of processors.
Definition: parallel.h:181
int _extra_order
Extra order to use for quadrature rule.
Real h2_error(const std::string &sys_name, const std::string &unknown_name)
void attach_exact_values(const std::vector< FunctionBase< Number > * > &f)
Clone and attach arbitrary functors which compute the exact values of the EquationSystems&#39; solutions ...
Real hdiv_error(const std::string &sys_name, const std::string &unknown_name)
static FEFieldType field_type(const FEType &fe_type)
const class libmesh_nullptr_t libmesh_nullptr
void attach_exact_hessians(std::vector< FunctionBase< Tensor > * > h)
Clone and attach arbitrary functors which compute the exact second derivatives of the EquationSystems...
The libMesh namespace provides an interface to certain functionality in the library.
std::vector< FunctionBase< Gradient > * > _exact_derivs
User-provided functors which compute the exact derivative of the solution for each system...
Real error_norm(const std::string &sys_name, const std::string &unknown_name, const FEMNormType &norm)
void 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(...
const EquationSystems * _equation_systems_fine
Constant pointer to the EquationSystems object containing the fine grid solution. ...
~ExactSolution()
Destructor.
const std::string & name() const
Definition: system.h:1998
This is the MeshBase class.
Definition: mesh_base.h:68
Gradient exact_grad(const Point &, const Parameters &, const std::string &, const std::string &)
Definition: mysystems.h:25
Number fptr(const Point &p, const Parameters &, const std::string &libmesh_dbg_var(sys_name), const std::string &unknown_name)
Definition: projection.C:78
libmesh_assert(j)
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
static UniquePtr< NumericVector< T > > build(const Parallel::Communicator &comm, const SolverPackage solver_package=libMesh::default_solver_package())
Builds a NumericVector on the processors in communicator comm using the linear solver package specifi...
std::map< std::string, std::vector< Real > > SystemErrorMap
Data structure which stores the errors: ||e|| = ||u - u_h|| ||grad(e)|| = ||grad(u) - grad(u_h)|| for...
virtual SimpleRange< element_iterator > active_local_element_ptr_range()=0
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:167
unsigned short int variable_number(const std::string &var) const
Definition: system.C:1263
void update_global_solution(std::vector< Number > &global_soln) const
Fill the input vector global_soln so that it contains the global solution on all processors.
Definition: system.C:662
const MeshBase & get_mesh() const
Definition: system.h:2014
const DofMap & get_dof_map() const
Definition: system.h:2030
const std::vector< Real > & get_JxW() const
Definition: fe_abstract.h:237
const Variable & variable(unsigned int var) const
Return a constant reference to Variable var.
Definition: system.h:2114
This is the base class for classes which contain information related to any physical process that mig...
Definition: system.h:76
const std::vector< std::vector< OutputShape > > & get_phi() const
Definition: fe_base.h:208
bool active_on_subdomain(subdomain_id_type sid) const
Definition: variable.h:136
const std::vector< std::vector< OutputShape > > & get_curl_phi() const
Definition: fe_base.h:224
UniquePtr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1523
Real norm_sq() const
Definition: type_vector.h:940
const std::set< unsigned char > & elem_dimensions() const
Definition: mesh_base.h:206
TensorTools::MakeNumber< OutputShape >::type OutputNumber
Definition: fe_base.h:124
void attach_exact_hessian(unsigned int sys_num, FunctionBase< Tensor > *h)
Clone and attach an arbitrary functor which computes the exact second derivatives of the system sys_n...
Real l1_error(const std::string &sys_name, const std::string &unknown_name)
static unsigned int n_vec_dim(const MeshBase &mesh, const FEType &fe_type)
void attach_exact_deriv(unsigned int sys_num, FunctionBase< Gradient > *g)
Clone and attach an arbitrary functor which computes the exact gradient of the system sys_num solutio...
void attach_exact_derivs(const std::vector< FunctionBase< Gradient > * > &g)
Clone and attach arbitrary functors which compute the exact gradients of the EquationSystems&#39; solutio...
unsigned int n_systems() const
Real hcurl_error(const std::string &sys_name, const std::string &unknown_name)
const std::vector< std::vector< OutputGradient > > & get_dphi() const
Definition: fe_base.h:216
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Number current_solution(const dof_id_type global_dof_number) const
Definition: system.C:192
std::map< std::string, SystemErrorMap > _errors
A map of SystemErrorMaps, which contains entries for each system in the EquationSystems object...
unsigned int number() const
Definition: system.h:2006
Real h1_error(const std::string &sys_name, const std::string &unknown_name)
const std::vector< std::vector< OutputTensor > > & get_d2phi() const
Definition: fe_base.h:290
Real l_inf_error(const std::string &sys_name, const std::string &unknown_name)
UniquePtr< QBase > default_quadrature_rule(const unsigned int dim, const int extraorder=0) const
Definition: fe_type.C:30
const Parallel::Communicator & comm() const
Gradient gptr(const Point &p, const Parameters &, const std::string &libmesh_dbg_var(sys_name), const std::string &unknown_name)
Definition: projection.C:93
virtual UniquePtr< FunctionBase< Output > > clone() const =0
Parameters parameters
Data structure holding arbitrary parameters.
unsigned int variable_scalar_number(const std::string &var, unsigned int component) const
Definition: system.h:2145
const T_sys & get_system(const std::string &name) const
void attach_exact_value(unsigned int sys_num, FunctionBase< Number > *f)
Clone and attach an arbitrary functor which computes the exact value of the system sys_num solution a...
unsigned int n_vars() const
Definition: system.h:2086
This class provides function-like objects for data distributed over a mesh.
Definition: mesh_function.h:53
void attach_reference_solution(const EquationSystems *es_fine)
Attach function similar to system.h which allows the user to attach a second EquationSystems object w...
Real l2_error(const std::string &sys_name, const std::string &unknown_name)
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
unsigned int n_points() const
Definition: quadrature.h:113
The QBase class provides the basic functionality from which various quadrature rules can be derived...
Definition: quadrature.h:53
const std::vector< Point > & get_xyz() const
Definition: fe_abstract.h:230
This class forms the foundation from which generic finite elements may be derived.
This class provides single index access to FieldType (i.e.
Definition: raw_accessor.h:93
Number curl_from_grad(const VectorValue< Number > &)
Definition: tensor_tools.C:28
This class defines a tensor in LIBMESH_DIM dimensional Real or Complex space.
virtual void reinit(const Elem *elem, const std::vector< Point > *const pts=libmesh_nullptr, const std::vector< Real > *const weights=libmesh_nullptr)=0
This is at the core of this class.
FEMNormType
defines an enum for norms defined on vectors of finite element coefficients
std::vector< Real > & _check_inputs(const std::string &sys_name, const std::string &unknown_name)
This function is responsible for checking the validity of the sys_name and unknown_name inputs...
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Fills the vector di with the global degree of freedom indices for the element.
Definition: dof_map.C:1917
static UniquePtr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.