exact_solution.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2018 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 // Local includes
19 #include "libmesh/dof_map.h"
20 #include "libmesh/elem.h"
21 #include "libmesh/exact_solution.h"
23 #include "libmesh/fe_base.h"
24 #include "libmesh/function_base.h"
25 #include "libmesh/mesh_base.h"
26 #include "libmesh/mesh_function.h"
27 #include "libmesh/numeric_vector.h"
28 #include "libmesh/parallel.h"
29 #include "libmesh/quadrature.h"
31 #include "libmesh/fe_interface.h"
32 #include "libmesh/raw_accessor.h"
33 #include "libmesh/tensor_tools.h"
34 #include "libmesh/enum_norm_type.h"
35 
36 namespace libMesh
37 {
38 
40  _equation_systems(es),
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 }
67 
68 
71 
72 
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 }
83 
84 
85 void ExactSolution::attach_exact_value (ValueFunctionPointer fptr)
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 }
102 
103 
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 }
113 
114 
115 void ExactSolution::attach_exact_value (unsigned int sys_num,
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 }
124 
125 
126 void ExactSolution::attach_exact_deriv (GradientFunctionPointer gptr)
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 }
143 
144 
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 }
154 
155 
156 void ExactSolution::attach_exact_deriv (unsigned int sys_num,
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 }
165 
166 
167 void ExactSolution::attach_exact_hessian (HessianFunctionPointer hptr)
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 }
184 
185 
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 }
195 
196 
197 void ExactSolution::attach_exact_hessian (unsigned int sys_num,
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 }
206 
207 
208 std::vector<Real> & ExactSolution::_check_inputs(const std::string & sys_name,
209  const std::string & unknown_name)
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 }
231 
232 
233 
234 void ExactSolution::compute_error(const std::string & sys_name,
235  const std::string & unknown_name)
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 }
266 
267 
268 
269 
270 
271 Real ExactSolution::error_norm(const std::string & sys_name,
272  const std::string & unknown_name,
273  const FEMNormType & norm)
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 }
333 
334 
335 
336 
337 
338 
339 
340 Real ExactSolution::l2_error(const std::string & sys_name,
341  const std::string & unknown_name)
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 }
353 
354 
355 
356 
357 
358 
359 
360 Real ExactSolution::l1_error(const std::string & sys_name,
361  const std::string & unknown_name)
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 }
373 
374 
375 
376 
377 
378 
379 
380 Real ExactSolution::l_inf_error(const std::string & sys_name,
381  const std::string & unknown_name)
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 }
393 
394 
395 
396 
397 
398 
399 
400 Real ExactSolution::h1_error(const std::string & sys_name,
401  const std::string & unknown_name)
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 }
415 
416 
417 Real ExactSolution::hcurl_error(const std::string & sys_name,
418  const std::string & unknown_name)
419 {
420  return this->error_norm(sys_name,unknown_name,HCURL);
421 }
422 
423 
424 Real ExactSolution::hdiv_error(const std::string & sys_name,
425  const std::string & unknown_name)
426 {
427  return this->error_norm(sys_name,unknown_name,HDIV);
428 }
429 
430 
431 
432 Real ExactSolution::h2_error(const std::string & sys_name,
433  const std::string & unknown_name)
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 }
447 
448 
449 
450 
451 
452 
453 
454 template<typename OutputShape>
455 void ExactSolution::_compute_error(const std::string & sys_name,
456  const std::string & unknown_name,
457  std::vector<Real> & error_vals)
458 {
459  // Make sure we aren't "overconfigured"
460  libmesh_assert (!(_exact_values.size() && _equation_systems_fine));
461 
462  // We need a communicator.
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>
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 }
820 
821 // Explicit instantiations of templated member functions
822 template void ExactSolution::_compute_error<Real>(const std::string &, const std::string &, std::vector<Real> &);
823 template void ExactSolution::_compute_error<RealGradient>(const std::string &, const std::string &, std::vector<Real> &);
824 
825 } // namespace libMesh
Manages the family, order, etc. parameters for a given FE.
Definition: fe_type.h:179
Wrap a libMesh-style function pointer into a FunctionBase object.
Manages multiples systems of equations.
void _compute_error(const std::string &sys_name, const std::string &unknown_name, std::vector< Real > &error_vals)
unsigned int variable_scalar_number(const std::string &var, unsigned int component) const
Definition: system.h:2164
void update_global_solution(std::vector< Number > &global_soln) const
Definition: system.C:642
unsigned int n_systems() const
const Variable & variable(unsigned int var) const
Definition: system.h:2133
const std::vector< Point > & get_xyz() const
Definition: fe_abstract.h:238
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
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Definition: dof_map.C:1930
ExactSolution(const EquationSystems &es)
void attach_exact_values(const std::vector< FunctionBase< Number > *> &f)
MPI_Comm communicator
Definition: communicator.h:57
const T_sys & get_system(const std::string &name) const
const FEType & variable_type(const unsigned int c) const
Definition: dof_map.h:1792
Real h2_error(const std::string &sys_name, const std::string &unknown_name)
Real hdiv_error(const std::string &sys_name, const std::string &unknown_name)
static FEFieldType field_type(const FEType &fe_type)
const Parallel::Communicator & comm() const
void attach_exact_derivs(const std::vector< FunctionBase< Gradient > *> &g)
const std::vector< std::vector< OutputDivergence > > & get_div_phi() const
Definition: fe_base.h:231
Real error_norm(const std::string &sys_name, const std::string &unknown_name, const FEMNormType &norm)
const MeshBase & get_mesh() const
Definition: system.h:2033
void compute_error(const std::string &sys_name, const std::string &unknown_name)
const EquationSystems * _equation_systems_fine
Base class for Mesh.
Definition: mesh_base.h:77
Number current_solution(const dof_id_type global_dof_number) const
Definition: system.C:194
std::map< std::string, std::vector< Real > > SystemErrorMap
virtual SimpleRange< element_iterator > active_local_element_ptr_range()=0
Manages the degrees of freedom (DOFs) in a simulation.
Definition: dof_map.h:176
std::unique_ptr< QBase > default_quadrature_rule(const unsigned int dim, const int extraorder=0) const
Definition: fe_type.C:31
const std::vector< std::vector< OutputGradient > > & get_dphi() const
Definition: fe_base.h:215
unsigned int number() const
Definition: system.h:2025
void attach_exact_hessians(std::vector< FunctionBase< Tensor > *> h)
Manages consistently variables, degrees of freedom, and coefficient vectors.
Definition: system.h:92
unsigned short int variable_number(const std::string &var) const
Definition: system.C:1243
const std::vector< Real > & get_JxW() const
Definition: fe_abstract.h:245
std::unique_ptr< NumericVector< Number > > solution
Definition: system.h:1523
static std::unique_ptr< FEGenericBase > build(const unsigned int dim, const FEType &type)
std::vector< std::unique_ptr< FunctionBase< Gradient > > > _exact_derivs
virtual void reinit(const Elem *elem, const std::vector< Point > *const pts=nullptr, const std::vector< Real > *const weights=nullptr)=0
const std::string & variable_name(const unsigned int i) const
Definition: system.h:2153
unsigned int n_points() const
Definition: quadrature.h:127
const std::set< unsigned char > & elem_dimensions() const
Definition: mesh_base.h:220
OStreamProxy err(std::cerr)
const std::vector< std::vector< OutputTensor > > & get_d2phi() const
Definition: fe_base.h:289
bool active_on_subdomain(subdomain_id_type sid) const
Definition: variable.h:136
TensorTools::MakeNumber< OutputShape >::type OutputNumber
Definition: fe_base.h:123
void attach_exact_hessian(unsigned int sys_num, FunctionBase< Tensor > *h)
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)
static std::unique_ptr< NumericVector< T > > build(const Parallel::Communicator &comm, const SolverPackage solver_package=libMesh::default_solver_package())
void attach_exact_deriv(unsigned int sys_num, FunctionBase< Gradient > *g)
std::vector< std::unique_ptr< FunctionBase< Tensor > > > _exact_hessians
const std::vector< std::vector< OutputShape > > & get_curl_phi() const
Definition: fe_base.h:223
Real hcurl_error(const std::string &sys_name, const std::string &unknown_name)
Real norm_sq() const
Definition: type_vector.h:943
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
std::map< std::string, SystemErrorMap > _errors
Real h1_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)
bool has_system(const std::string &name) const
virtual std::unique_ptr< FunctionBase< Output > > clone() const =0
void attach_exact_value(unsigned int sys_num, FunctionBase< Number > *f)
const std::string & name() const
Definition: system.h:2017
void attach_reference_solution(const EquationSystems *es_fine)
Real l2_error(const std::string &sys_name, const std::string &unknown_name)
unsigned int n_vars() const
Definition: system.h:2105
const DofMap & get_dof_map() const
Definition: system.h:2049
Base class for all quadrature families and orders.
Definition: quadrature.h:62
std::vector< std::unique_ptr< FunctionBase< Number > > > _exact_values
Number curl_from_grad(const VectorValue< Number > &)
Definition: tensor_tools.C:28
std::vector< Real > & _check_inputs(const std::string &sys_name, const std::string &unknown_name)
const std::vector< std::vector< OutputShape > > & get_phi() const
Definition: fe_base.h:207