exact_solution.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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"
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"
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(NULL),
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 (unsigned int i=0; i != _exact_values.size(); ++i)
77  delete (_exact_values[i]);
78 
79  for (unsigned int i=0; i != _exact_derivs.size(); ++i)
80  delete (_exact_derivs[i]);
81 
82  for (unsigned int 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 {
104  libmesh_assert(fptr);
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
119  _equation_systems_fine = NULL;
120 }
121 
122 
124 {
125  // Clear out any previous _exact_values entries, then add a new
126  // entry for each system.
127  for (unsigned int i=0; i != _exact_values.size(); ++i)
128  delete (_exact_values[i]);
129 
130  _exact_values.clear();
131  _exact_values.resize(f.size(), NULL);
132 
133  // We use clone() to get non-sliced copies of FunctionBase
134  // subclasses, but we can't put the resulting AutoPtrs into an STL
135  // container.
136  for (unsigned int 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, NULL);
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 {
158  libmesh_assert(gptr);
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
173  _equation_systems_fine = NULL;
174 }
175 
176 
178 {
179  // Clear out any previous _exact_derivs entries, then add a new
180  // entry for each system.
181  for (unsigned int i=0; i != _exact_derivs.size(); ++i)
182  delete (_exact_derivs[i]);
183 
184  _exact_derivs.clear();
185  _exact_derivs.resize(g.size(), NULL);
186 
187  // We use clone() to get non-sliced copies of FunctionBase
188  // subclasses, but we can't put the resulting AutoPtrs into an STL
189  // container.
190  for (unsigned int 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, NULL);
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
227  _equation_systems_fine = NULL;
228 }
229 
230 
232 {
233  // Clear out any previous _exact_hessians entries, then add a new
234  // entry for each system.
235  for (unsigned int i=0; i != _exact_hessians.size(); ++i)
236  delete (_exact_hessians[i]);
237 
238  _exact_hessians.clear();
239  _exact_hessians.resize(h.size(), NULL);
240 
241  // We use clone() to get non-sliced copies of FunctionBase
242  // subclasses, but we can't put the resulting AutoPtrs into an STL
243  // container.
244  for (unsigned int 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, NULL);
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  {
274  libMesh::err << "Sorry, couldn't find the requested system '"
275  << sys_name << "'."
276  << std::endl;
277  libmesh_error();
278  }
279 
280  // Make sure the requested unknown_name exists.
281  SystemErrorMap::iterator var_iter = (*sys_iter).second.find(unknown_name);
282 
283  if (var_iter == (*sys_iter).second.end())
284  {
285  libMesh::err << "Sorry, couldn't find the requested variable '"
286  << unknown_name << "'."
287  << std::endl;
288  libmesh_error();
289  }
290 
291  // Return a reference to the proper error entry
292  return (*var_iter).second;
293 }
294 
295 
296 
297 void ExactSolution::compute_error(const std::string& sys_name,
298  const std::string& unknown_name)
299 {
300  // Check the inputs for validity, and get a reference
301  // to the proper location to store the error
302  std::vector<Real>& error_vals = this->_check_inputs(sys_name,
303  unknown_name);
304 
306  const System& sys = _equation_systems.get_system<System>( sys_name );
307 
308  libmesh_assert( sys.has_variable( unknown_name ) );
309  switch( FEInterface::field_type(sys.variable_type( unknown_name )) )
310  {
311  case TYPE_SCALAR:
312  {
313  this->_compute_error<Real>(sys_name,
314  unknown_name,
315  error_vals);
316  break;
317  }
318  case TYPE_VECTOR:
319  {
320  this->_compute_error<RealGradient>(sys_name,
321  unknown_name,
322  error_vals);
323  break;
324  }
325  default:
326  libmesh_error();
327  }
328 
329  return;
330 }
331 
332 
333 
334 
335 
336 Real ExactSolution::error_norm(const std::string& sys_name,
337  const std::string& unknown_name,
338  const FEMNormType& norm)
339 {
340  // Check the inputs for validity, and get a reference
341  // to the proper location to store the error
342  std::vector<Real>& error_vals = this->_check_inputs(sys_name,
343  unknown_name);
344 
346  libmesh_assert(_equation_systems.get_system(sys_name).has_variable( unknown_name ));
347  const FEType& fe_type = _equation_systems.get_system(sys_name).variable_type(unknown_name);
348 
349  switch (norm)
350  {
351  case L2:
352  return std::sqrt(error_vals[0]);
353  case H1:
354  return std::sqrt(error_vals[0] + error_vals[1]);
355  case H2:
356  return std::sqrt(error_vals[0] + error_vals[1] + error_vals[2]);
357  case HCURL:
358  {
359  if(FEInterface::field_type(fe_type) == TYPE_SCALAR)
360  {
361  libMesh::err << "Cannot compute HCurl error norm of scalar-valued variables!" << std::endl;
362  libmesh_error();
363  }
364  else
365  return std::sqrt(error_vals[0] + error_vals[5]);
366  }
367  case HDIV:
368  {
369  if(FEInterface::field_type(fe_type) == TYPE_SCALAR)
370  {
371  libMesh::err << "Cannot compute HDiv error norm of scalar-valued variables!" << std::endl;
372  libmesh_error();
373  }
374  else
375  return std::sqrt(error_vals[0] + error_vals[6]);
376  }
377  case H1_SEMINORM:
378  return std::sqrt(error_vals[1]);
379  case H2_SEMINORM:
380  return std::sqrt(error_vals[2]);
381  case HCURL_SEMINORM:
382  {
383  if(FEInterface::field_type(fe_type) == TYPE_SCALAR)
384  {
385  libMesh::err << "Cannot compute HCurl error seminorm of scalar-valued variables!" << std::endl;
386  libmesh_error();
387  }
388  else
389  return std::sqrt(error_vals[5]);
390  }
391  case HDIV_SEMINORM:
392  {
393  if(FEInterface::field_type(fe_type) == TYPE_SCALAR)
394  {
395  libMesh::err << "Cannot compute HDiv error seminorm of scalar-valued variables!" << std::endl;
396  libmesh_error();
397  }
398  else
399  return std::sqrt(error_vals[6]);
400  }
401  case L1:
402  return error_vals[3];
403  case L_INF:
404  return error_vals[4];
405 
406  // Currently only Sobolev norms/seminorms are supported
407  default:
408  libmesh_error();
409  }
410 }
411 
412 
413 
414 
415 
416 
417 
418 Real ExactSolution::l2_error(const std::string& sys_name,
419  const std::string& unknown_name)
420 {
421 
422  // Check the inputs for validity, and get a reference
423  // to the proper location to store the error
424  std::vector<Real>& error_vals = this->_check_inputs(sys_name,
425  unknown_name);
426 
427  // Return the square root of the first component of the
428  // computed error.
429  return std::sqrt(error_vals[0]);
430 }
431 
432 
433 
434 
435 
436 
437 
438 Real ExactSolution::l1_error(const std::string& sys_name,
439  const std::string& unknown_name)
440 {
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 first component of the
448  // computed error.
449  return error_vals[3];
450 }
451 
452 
453 
454 
455 
456 
457 
458 Real ExactSolution::l_inf_error(const std::string& sys_name,
459  const std::string& unknown_name)
460 {
461 
462  // Check the inputs for validity, and get a reference
463  // to the proper location to store the error
464  std::vector<Real>& error_vals = this->_check_inputs(sys_name,
465  unknown_name);
466 
467  // Return the square root of the first component of the
468  // computed error.
469  return error_vals[4];
470 }
471 
472 
473 
474 
475 
476 
477 
478 Real ExactSolution::h1_error(const std::string& sys_name,
479  const std::string& unknown_name)
480 {
481  // If the user has supplied no exact derivative function, we
482  // just integrate the H1 norm of the solution; i.e. its
483  // difference from an "exact solution" of zero.
484 
485  // Check the inputs for validity, and get a reference
486  // to the proper location to store the error
487  std::vector<Real>& error_vals = this->_check_inputs(sys_name,
488  unknown_name);
489 
490  // Return the square root of the sum of the computed errors.
491  return std::sqrt(error_vals[0] + error_vals[1]);
492 }
493 
494 
495 Real ExactSolution::hcurl_error(const std::string& sys_name,
496  const std::string& unknown_name)
497 {
498  return this->error_norm(sys_name,unknown_name,HCURL);
499 }
500 
501 
502 Real ExactSolution::hdiv_error(const std::string& sys_name,
503  const std::string& unknown_name)
504 {
505  return this->error_norm(sys_name,unknown_name,HDIV);
506 }
507 
508 
509 
510 Real ExactSolution::h2_error(const std::string& sys_name,
511  const std::string& unknown_name)
512 {
513  // If the user has supplied no exact derivative functions, we
514  // just integrate the H2 norm of the solution; i.e. its
515  // difference from an "exact solution" of zero.
516 
517  // Check the inputs for validity, and get a reference
518  // to the proper location to store the error
519  std::vector<Real>& error_vals = this->_check_inputs(sys_name,
520  unknown_name);
521 
522  // Return the square root of the sum of the computed errors.
523  return std::sqrt(error_vals[0] + error_vals[1] + error_vals[2]);
524 }
525 
526 
527 
528 
529 
530 
531 
532 template< typename OutputShape>
533 void ExactSolution::_compute_error(const std::string& sys_name,
534  const std::string& unknown_name,
535  std::vector<Real>& error_vals)
536 {
537  // Make sure we aren't "overconfigured"
539 
540  // We need a commmunicator.
542 
543  // This function must be run on all processors at once
545 
546  // Get a reference to the system whose error is being computed.
547  // If we have a fine grid, however, we'll integrate on that instead
548  // for more accuracy.
549  const System& computed_system = _equation_systems_fine ?
551  _equation_systems.get_system (sys_name);
552 
553  const Real time = _equation_systems.get_system(sys_name).time;
554 
555  const unsigned int sys_num = computed_system.number();
556  const unsigned int var = computed_system.variable_number(unknown_name);
557  const unsigned int var_component =
558  computed_system.variable_scalar_number(var, 0);
559 
560  // Prepare a global solution and a MeshFunction of the coarse system if we need one
561  AutoPtr<MeshFunction> coarse_values;
564  {
565  const System& comparison_system
566  = _equation_systems.get_system(sys_name);
567 
568  std::vector<Number> global_soln;
569  comparison_system.update_global_solution(global_soln);
570  comparison_soln->init(comparison_system.solution->size(), true, SERIAL);
571  (*comparison_soln) = global_soln;
572 
573  coarse_values = AutoPtr<MeshFunction>
575  *comparison_soln,
576  comparison_system.get_dof_map(),
577  comparison_system.variable_number(unknown_name)));
578  coarse_values->init();
579  }
580 
581  // Initialize any functors we're going to use
582  for (unsigned int i=0; i != _exact_values.size(); ++i)
583  if (_exact_values[i])
584  _exact_values[i]->init();
585 
586  for (unsigned int i=0; i != _exact_derivs.size(); ++i)
587  if (_exact_derivs[i])
588  _exact_derivs[i]->init();
589 
590  for (unsigned int i=0; i != _exact_hessians.size(); ++i)
591  if (_exact_hessians[i])
592  _exact_hessians[i]->init();
593 
594  // Get a reference to the dofmap and mesh for that system
595  const DofMap& computed_dof_map = computed_system.get_dof_map();
596 
597  const MeshBase& _mesh = computed_system.get_mesh();
598 
599  const unsigned int dim = _mesh.mesh_dimension();
600 
601  // Zero the error before summation
602  // 0 - sum of square of function error (L2)
603  // 1 - sum of square of gradient error (H1 semi)
604  // 2 - sum of square of Hessian error (H2 semi)
605  // 3 - sum of sqrt(square of function error) (L1)
606  // 4 - max of sqrt(square of function error) (Linfty)
607  // 5 - sum of square of curl error (HCurl semi)
608  // 6 - sum of square of div error (HDiv semi)
609  error_vals = std::vector<Real>(7, 0.);
610 
611  // Construct Quadrature rule based on default quadrature order
612  const FEType& fe_type = computed_dof_map.variable_type(var);
613 
614  unsigned int n_vec_dim = FEInterface::n_vec_dim( _mesh, fe_type );
615 
616  // FIXME: MeshFunction needs to be updated to support vector-valued
617  // elements before we can use a reference solution.
618  if( (n_vec_dim > 1) && _equation_systems_fine )
619  {
620  libMesh::err << "Error calculation using reference solution not yet\n"
621  << "supported for vector-valued elements."
622  << std::endl;
623  libmesh_not_implemented();
624  }
625 
626  AutoPtr<QBase> qrule =
627  fe_type.default_quadrature_rule (dim,
628  _extra_order);
629 
630  // Construct finite element object
631 
633 
634  // Attach quadrature rule to FE object
635  fe->attach_quadrature_rule (qrule.get());
636 
637  // The Jacobian*weight at the quadrature points.
638  const std::vector<Real>& JxW = fe->get_JxW();
639 
640  // The value of the shape functions at the quadrature points
641  // i.e. phi(i) = phi_values[i][qp]
642  const std::vector<std::vector<OutputShape> >& phi_values = fe->get_phi();
643 
644  // The value of the shape function gradients at the quadrature points
645  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputGradient> >&
646  dphi_values = fe->get_dphi();
647 
648  // The value of the shape function curls at the quadrature points
649  // Only computed for vector-valued elements
650  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputShape> >* curl_values = NULL;
651 
652  // The value of the shape function divergences at the quadrature points
653  // Only computed for vector-valued elements
654  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputDivergence> >* div_values = NULL;
655 
656  if( FEInterface::field_type(fe_type) == TYPE_VECTOR )
657  {
658  curl_values = &fe->get_curl_phi();
659  div_values = &fe->get_div_phi();
660  }
661 
662 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
663  // The value of the shape function second derivatives at the quadrature points
664  const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputTensor> >&
665  d2phi_values = fe->get_d2phi();
666 #endif
667 
668  // The XYZ locations (in physical space) of the quadrature points
669  const std::vector<Point>& q_point = fe->get_xyz();
670 
671  // The global degree of freedom indices associated
672  // with the local degrees of freedom.
673  std::vector<dof_id_type> dof_indices;
674 
675 
676  //
677  // Begin the loop over the elements
678  //
679  // TODO: this ought to be threaded (and using subordinate
680  // MeshFunction objects in each thread rather than a single
681  // master)
684 
685  for ( ; el != end_el; ++el)
686  {
687  // Store a pointer to the element we are currently
688  // working on. This allows for nicer syntax later.
689  const Elem* elem = *el;
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  libmesh_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 reqiures 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  exact_val = (*coarse_values)(q_point[qp]);
755  }
756  const typename FEGenericBase<OutputShape>::OutputNumber val_error = u_h - exact_val;
757 
758  // Add the squares of the error to each contribution
759  Real error_sq = TensorTools::norm_sq(val_error);
760  error_vals[0] += JxW[qp]*error_sq;
761 
762  Real norm = sqrt(error_sq);
763  error_vals[3] += JxW[qp]*norm;
764 
765  if(error_vals[4]<norm) { error_vals[4] = norm; }
766 
767  // Compute the value of the error in the gradient at this
768  // quadrature point
771  if (_exact_derivs.size() > sys_num && _exact_derivs[sys_num])
772  {
773  for( unsigned int c = 0; c < n_vec_dim; c++)
774  for( unsigned int d = 0; d < _mesh.spatial_dimension(); d++ )
775  exact_grad_accessor(d + c*_mesh.spatial_dimension() ) =
776  _exact_derivs[sys_num]->
777  component(var_component+c, q_point[qp], time)(d);
778  }
779  else if (_equation_systems_fine)
780  {
781  // FIXME: Needs to be updated for vector-valued elements
782  exact_grad = coarse_values->gradient(q_point[qp]);
783  }
784 
785  const typename FEGenericBase<OutputShape>::OutputNumberGradient grad_error = grad_u_h - exact_grad;
786 
787  error_vals[1] += JxW[qp]*grad_error.size_sq();
788 
789 
790  if( FEInterface::field_type(fe_type) == TYPE_VECTOR )
791  {
792  // Compute the value of the error in the curl at this
793  // quadrature point
794  typename FEGenericBase<OutputShape>::OutputNumber exact_curl = 0.0;
795  if (_exact_derivs.size() > sys_num && _exact_derivs[sys_num])
796  {
797  exact_curl = TensorTools::curl_from_grad( exact_grad );
798  }
799  else if (_equation_systems_fine)
800  {
801  // FIXME: Need to implement curl for MeshFunction and support reference
802  // solution for vector-valued elements
803  }
804 
805  const typename FEGenericBase<OutputShape>::OutputNumber curl_error = curl_u_h - exact_curl;
806 
807  error_vals[5] += JxW[qp]*TensorTools::norm_sq(curl_error);
808 
809  // Compute the value of the error in the divergence at this
810  // quadrature point
812  if (_exact_derivs.size() > sys_num && _exact_derivs[sys_num])
813  {
814  exact_div = TensorTools::div_from_grad( exact_grad );
815  }
816  else if (_equation_systems_fine)
817  {
818  // FIXME: Need to implement div for MeshFunction and support reference
819  // solution for vector-valued elements
820  }
821 
822  const typename FEGenericBase<OutputShape>::OutputNumberDivergence div_error = div_u_h - exact_div;
823 
824  error_vals[6] += JxW[qp]*TensorTools::norm_sq(div_error);
825  }
826 
827 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
828  // Compute the value of the error in the hessian at this
829  // quadrature point
831  RawAccessor<typename FEGenericBase<OutputShape>::OutputNumberTensor> exact_hess_accessor( exact_hess, dim );
832  if (_exact_hessians.size() > sys_num && _exact_hessians[sys_num])
833  {
834  //FIXME: This needs to be implemented to support rank 3 tensors
835  // which can't happen until type_n_tensor is fully implemented
836  // and a RawAccessor<TypeNTensor> is fully implemented
837  if( FEInterface::field_type(fe_type) == TYPE_VECTOR )
838  libmesh_not_implemented();
839 
840  for( unsigned int c = 0; c < n_vec_dim; c++)
841  for( unsigned int d = 0; d < dim; d++ )
842  for( unsigned int e =0; e < dim; e++ )
843  exact_hess_accessor(d + e*dim + c*dim*dim) =
844  _exact_hessians[sys_num]->
845  component(var_component+c, q_point[qp], time)(d,e);
846  }
847  else if (_equation_systems_fine)
848  {
849  // FIXME: Needs to be updated for vector-valued elements
850  exact_hess = coarse_values->hessian(q_point[qp]);
851  }
852 
853  const typename FEGenericBase<OutputShape>::OutputNumberTensor grad2_error = grad2_u_h - exact_hess;
854 
855  // FIXME: PB: Is this what we want for rank 3 tensors?
856  error_vals[2] += JxW[qp]*grad2_error.size_sq();
857 #endif
858 
859  } // end qp loop
860  } // end element loop
861 
862  // Add up the error values on all processors, except for the L-infty
863  // norm, for which the maximum is computed.
864  Real l_infty_norm = error_vals[4];
865  communicator.max(l_infty_norm);
866  communicator.sum(error_vals);
867  error_vals[4] = l_infty_norm;
868 }
869 
870  // Explicit instantiations of templated member functions
871  template void ExactSolution::_compute_error<Real>(const std::string&, const std::string&, std::vector<Real>&);
872  template void ExactSolution::_compute_error<RealGradient>(const std::string&, const std::string&, std::vector<Real>&);
873 
874 } // namespace libMesh

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:04 UTC

Hosted By:
SourceForge.net Logo