exact_error_estimator.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 
19 // C++ includes
20 #include <algorithm> // for std::fill
21 #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
22 #include <cmath> // for sqrt
23 
24 
25 // Local Includes
26 #include "libmesh/libmesh_common.h"
28 #include "libmesh/dof_map.h"
30 #include "libmesh/error_vector.h"
31 #include "libmesh/fe_base.h"
33 #include "libmesh/elem.h"
34 #include "libmesh/mesh_base.h"
35 #include "libmesh/mesh_function.h"
36 #include "libmesh/numeric_vector.h"
37 #include "libmesh/quadrature.h"
38 #include "libmesh/system.h"
39 #include "libmesh/tensor_tools.h"
40 
41 namespace libMesh
42 {
43 
44 //-----------------------------------------------------------------
45 // ErrorEstimator implementations
46 void ExactErrorEstimator::attach_exact_value (Number fptr(const Point& p,
47  const Parameters& parameters,
48  const std::string& sys_name,
49  const std::string& unknown_name))
50 {
51  libmesh_assert(fptr);
52  _exact_value = fptr;
53 
54  // We're not using a fine grid solution
56 
57  // We're not using user-provided functors
58  this->clear_functors();
59 }
60 
61 
63 {
64  // Clear out any previous _exact_values entries, then add a new
65  // entry for each system.
66  for (unsigned int i=0; i != _exact_values.size(); ++i)
67  delete (_exact_values[i]);
68 
69  _exact_values.clear();
70  _exact_values.resize(f.size(), NULL);
71 
72  // We use clone() to get non-sliced copies of FunctionBase
73  // subclasses, but we can't put the resulting AutoPtrs into an STL
74  // container.
75  for (unsigned int i=0; i != f.size(); ++i)
76  if (f[i])
77  _exact_values[i] = f[i]->clone().release();
78 }
79 
80 
81 void ExactErrorEstimator::attach_exact_value (unsigned int sys_num,
83 {
84  if (_exact_values.size() <= sys_num)
85  _exact_values.resize(sys_num+1, NULL);
86 
87  if (f)
88  _exact_values[sys_num] = f->clone().release();
89 }
90 
91 
93  const Parameters& parameters,
94  const std::string& sys_name,
95  const std::string& unknown_name))
96 {
97  libmesh_assert(gptr);
98  _exact_deriv = gptr;
99 
100  // We're not using a fine grid solution
101  _equation_systems_fine = NULL;
102 
103  // We're not using user-provided functors
104  this->clear_functors();
105 }
106 
107 
109 {
110  // Clear out any previous _exact_derivs entries, then add a new
111  // entry for each system.
112  for (unsigned int i=0; i != _exact_derivs.size(); ++i)
113  delete (_exact_derivs[i]);
114 
115  _exact_derivs.clear();
116  _exact_derivs.resize(g.size(), NULL);
117 
118  // We use clone() to get non-sliced copies of FunctionBase
119  // subclasses, but we can't put the resulting AutoPtrs into an STL
120  // container.
121  for (unsigned int i=0; i != g.size(); ++i)
122  if (g[i])
123  _exact_derivs[i] = g[i]->clone().release();
124 }
125 
126 
127 void ExactErrorEstimator::attach_exact_deriv (unsigned int sys_num,
129 {
130  if (_exact_derivs.size() <= sys_num)
131  _exact_derivs.resize(sys_num+1, NULL);
132 
133  if (g)
134  _exact_derivs[sys_num] = g->clone().release();
135 }
136 
137 
138 
139 
141  const Parameters& parameters,
142  const std::string& sys_name,
143  const std::string& unknown_name))
144 {
145  libmesh_assert(hptr);
146  _exact_hessian = hptr;
147 
148  // We're not using a fine grid solution
149  _equation_systems_fine = NULL;
150 
151  // We're not using user-provided functors
152  this->clear_functors();
153 }
154 
155 
157 {
158  // Clear out any previous _exact_hessians entries, then add a new
159  // entry for each system.
160  for (unsigned int i=0; i != _exact_hessians.size(); ++i)
161  delete (_exact_hessians[i]);
162 
163  _exact_hessians.clear();
164  _exact_hessians.resize(h.size(), NULL);
165 
166  // We use clone() to get non-sliced copies of FunctionBase
167  // subclasses, but we can't put the resulting AutoPtrs into an STL
168  // container.
169  for (unsigned int i=0; i != h.size(); ++i)
170  if (h[i])
171  _exact_hessians[i] = h[i]->clone().release();
172 }
173 
174 
175 void ExactErrorEstimator::attach_exact_hessian (unsigned int sys_num,
177 {
178  if (_exact_hessians.size() <= sys_num)
179  _exact_hessians.resize(sys_num+1, NULL);
180 
181  if (h)
182  _exact_hessians[sys_num] = h->clone().release();
183 }
184 
185 
187 {
188  libmesh_assert(es_fine);
189  _equation_systems_fine = es_fine;
190 
191  // If we're using a fine grid solution, we're not using exact value
192  // function pointers or functors.
193  _exact_value = NULL;
194  _exact_deriv = NULL;
195  _exact_hessian = NULL;
196 
197  this->clear_functors();
198 }
199 
200 #ifdef LIBMESH_ENABLE_AMR
201 void ExactErrorEstimator::estimate_error (const System& system,
202  ErrorVector& error_per_cell,
203  const NumericVector<Number>* solution_vector,
204  bool estimate_parent_error)
205 #else
206 void ExactErrorEstimator::estimate_error (const System& system,
207  ErrorVector& error_per_cell,
208  const NumericVector<Number>* solution_vector,
209  bool /* estimate_parent_error */ )
210 #endif
211 {
212  // The current mesh
213  const MeshBase& mesh = system.get_mesh();
214 
215  // The dimensionality of the mesh
216  const unsigned int dim = mesh.mesh_dimension();
217 
218  // The number of variables in the system
219  const unsigned int n_vars = system.n_vars();
220 
221  // The DofMap for this system
222  const DofMap& dof_map = system.get_dof_map();
223 
224  // Resize the error_per_cell vector to be
225  // the number of elements, initialize it to 0.
226  error_per_cell.resize (mesh.max_elem_id());
227  std::fill (error_per_cell.begin(), error_per_cell.end(), 0.);
228 
229  // Prepare current_local_solution to localize a non-standard
230  // solution vector if necessary
231  if (solution_vector && solution_vector != system.solution.get())
232  {
233  NumericVector<Number>* newsol =
234  const_cast<NumericVector<Number>*>(solution_vector);
235  System &sys = const_cast<System&>(system);
236  newsol->swap(*sys.solution);
237  sys.update();
238  }
239 
240  // Loop over all the variables in the system
241  for (unsigned int var=0; var<n_vars; var++)
242  {
243  // Possibly skip this variable
244  if (error_norm.weight(var) == 0.0) continue;
245 
246  // The (string) name of this variable
247  const std::string& var_name = system.variable_name(var);
248 
249  // The type of finite element to use for this variable
250  const FEType& fe_type = dof_map.variable_type (var);
251 
252  AutoPtr<FEBase> fe (FEBase::build (dim, fe_type));
253 
254  // Build an appropriate Gaussian quadrature rule
255  AutoPtr<QBase> qrule =
256  fe_type.default_quadrature_rule (dim,
257  _extra_order);
258 
259  fe->attach_quadrature_rule (qrule.get());
260 
261  // Prepare a global solution and a MeshFunction of the fine system if we need one
262  AutoPtr<MeshFunction> fine_values;
264  if (_equation_systems_fine)
265  {
266  const System& fine_system = _equation_systems_fine->get_system(system.name());
267 
268  std::vector<Number> global_soln;
269  // FIXME - we're assuming that the fine system solution gets
270  // used even when a different vector is used for the coarse
271  // system
272  fine_system.update_global_solution(global_soln);
273  fine_soln->init (global_soln.size(), true, SERIAL);
274  (*fine_soln) = global_soln;
275 
276  fine_values = AutoPtr<MeshFunction>
277  (new MeshFunction(*_equation_systems_fine,
278  *fine_soln,
279  fine_system.get_dof_map(),
280  fine_system.variable_number(var_name)));
281  fine_values->init();
282  } else {
283  // Initialize functors if we're using them
284  for (unsigned int i=0; i != _exact_values.size(); ++i)
285  if (_exact_values[i])
286  _exact_values[i]->init();
287 
288  for (unsigned int i=0; i != _exact_derivs.size(); ++i)
289  if (_exact_derivs[i])
290  _exact_derivs[i]->init();
291 
292  for (unsigned int i=0; i != _exact_hessians.size(); ++i)
293  if (_exact_hessians[i])
294  _exact_hessians[i]->init();
295  }
296 
297  // Request the data we'll need to compute with
298  fe->get_JxW();
299  fe->get_phi();
300  fe->get_dphi();
301 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
302  fe->get_d2phi();
303 #endif
304  fe->get_xyz();
305 
306  // If we compute on parent elements, we'll want to do so only
307  // once on each, so we need to keep track of which we've done.
308  std::vector<bool> computed_var_on_parent;
309 
310 #ifdef LIBMESH_ENABLE_AMR
311  if (estimate_parent_error)
312  computed_var_on_parent.resize(error_per_cell.size(), false);
313 #endif
314 
315  // TODO: this ought to be threaded (and using subordinate
316  // MeshFunction objects in each thread rather than a single
317  // master)
318 
319  // Iterate over all the active elements in the mesh
320  // that live on this processor.
322  elem_it = mesh.active_local_elements_begin();
324  elem_end = mesh.active_local_elements_end();
325 
326  for (;elem_it != elem_end; ++elem_it)
327  {
328  // e is necessarily an active element on the local processor
329  const Elem* elem = *elem_it;
330  const dof_id_type e_id = elem->id();
331 
332 #ifdef LIBMESH_ENABLE_AMR
333  // See if the parent of element e has been examined yet;
334  // if not, we may want to compute the estimator on it
335  const Elem* parent = elem->parent();
336 
337  // We only can compute and only need to compute on
338  // parents with all active children
339  bool compute_on_parent = true;
340  if (!parent || !estimate_parent_error)
341  compute_on_parent = false;
342  else
343  for (unsigned int c=0; c != parent->n_children(); ++c)
344  if (!parent->child(c)->active())
345  compute_on_parent = false;
346 
347  if (compute_on_parent &&
348  !computed_var_on_parent[parent->id()])
349  {
350  computed_var_on_parent[parent->id()] = true;
351 
352  // Compute a projection onto the parent
353  DenseVector<Number> Uparent;
355  dof_map, parent, Uparent,
356  var, false);
357 
358  error_per_cell[parent->id()] +=
359  static_cast<ErrorVectorReal>
360  (find_squared_element_error(system, var_name,
361  parent, Uparent,
362  fe.get(),
363  fine_values.get()));
364  }
365 #endif
366 
367  // Get the local to global degree of freedom maps
368  std::vector<dof_id_type> dof_indices;
369  dof_map.dof_indices (elem, dof_indices, var);
370  const unsigned int n_dofs =
371  libmesh_cast_int<unsigned int>(dof_indices.size());
372  DenseVector<Number> Uelem(n_dofs);
373  for (unsigned int i=0; i != n_dofs; ++i)
374  Uelem(i) = system.current_solution(dof_indices[i]);
375 
376  error_per_cell[e_id] +=
377  static_cast<ErrorVectorReal>
378  (find_squared_element_error(system, var_name, elem,
379  Uelem, fe.get(),
380  fine_values.get()));
381 
382  } // End loop over active local elements
383  } // End loop over variables
384 
385 
386 
387  // Each processor has now computed the error contribuions
388  // for its local elements. We need to sum the vector
389  // and then take the square-root of each component. Note
390  // that we only need to sum if we are running on multiple
391  // processors, and we only need to take the square-root
392  // if the value is nonzero. There will in general be many
393  // zeros for the inactive elements.
394 
395  // First sum the vector of estimated error values
396  this->reduce_error(error_per_cell, system.comm());
397 
398  // Compute the square-root of each component.
399  START_LOG("std::sqrt()", "ExactErrorEstimator");
400  for (dof_id_type i=0; i<error_per_cell.size(); i++)
401  {
402 
403  if (error_per_cell[i] != 0.)
404  {
405  libmesh_assert_greater (error_per_cell[i], 0.);
406  error_per_cell[i] = std::sqrt(error_per_cell[i]);
407  }
408 
409 
410  }
411  STOP_LOG("std::sqrt()", "ExactErrorEstimator");
412 
413  // If we used a non-standard solution before, now is the time to fix
414  // the current_local_solution
415  if (solution_vector && solution_vector != system.solution.get())
416  {
417  NumericVector<Number>* newsol =
418  const_cast<NumericVector<Number>*>(solution_vector);
419  System &sys = const_cast<System&>(system);
420  newsol->swap(*sys.solution);
421  sys.update();
422  }
423 }
424 
425 
426 
428  const std::string& var_name,
429  const Elem *elem,
430  const DenseVector<Number> &Uelem,
431  FEBase *fe,
432  MeshFunction *fine_values) const
433 {
434  // The (string) name of this system
435  const std::string& sys_name = system.name();
436  const unsigned int sys_num = system.number();
437 
438  const unsigned int var = system.variable_number(var_name);
439  const unsigned int var_component =
440  system.variable_scalar_number(var, 0);
441 
442  const Parameters& parameters = system.get_equation_systems().parameters;
443 
444  // reinitialize the element-specific data
445  // for the current element
446  fe->reinit (elem);
447 
448  // Get the data we need to compute with
449  const std::vector<Real> & JxW = fe->get_JxW();
450  const std::vector<std::vector<Real> >& phi_values = fe->get_phi();
451  const std::vector<std::vector<RealGradient> >& dphi_values = fe->get_dphi();
452  const std::vector<Point>& q_point = fe->get_xyz();
453 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
454  const std::vector<std::vector<RealTensor> >& d2phi_values = fe->get_d2phi();
455 #endif
456 
457  // The number of shape functions
458  const unsigned int n_sf =
459  libmesh_cast_int<unsigned int>(Uelem.size());
460 
461  // The number of quadrature points
462  const unsigned int n_qp =
463  libmesh_cast_int<unsigned int>(JxW.size());
464 
465  Real error_val = 0;
466 
467  // Begin the loop over the Quadrature points.
468  //
469  for (unsigned int qp=0; qp<n_qp; qp++)
470  {
471  // Real u_h = 0.;
472  // RealGradient grad_u_h;
473 
474  Number u_h = 0.;
475 
476  Gradient grad_u_h;
477 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
478  Tensor grad2_u_h;
479 #endif
480 
481  // Compute solution values at the current
482  // quadrature point. This reqiures a sum
483  // over all the shape functions evaluated
484  // at the quadrature point.
485  for (unsigned int i=0; i<n_sf; i++)
486  {
487  // Values from current solution.
488  u_h += phi_values[i][qp]*Uelem(i);
489  grad_u_h += dphi_values[i][qp]*Uelem(i);
490 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
491  grad2_u_h += d2phi_values[i][qp]*Uelem(i);
492 #endif
493  }
494 
495  // Compute the value of the error at this quadrature point
496  if (error_norm.type(var) == L2 ||
497  error_norm.type(var) == H1 ||
498  error_norm.type(var) == H2)
499  {
500  Number val_error = u_h;
501  if (_exact_value)
502  val_error -= _exact_value(q_point[qp],parameters,sys_name,var_name);
503  else if (_exact_values.size() > sys_num && _exact_values[sys_num])
504  val_error -= _exact_values[sys_num]->
505  component(var_component, q_point[qp], system.time);
506  else if (_equation_systems_fine)
507  val_error -= (*fine_values)(q_point[qp]);
508 
509  // Add the squares of the error to each contribution
510  error_val += JxW[qp]*TensorTools::norm_sq(val_error);
511  }
512 
513  // Compute the value of the error in the gradient at this
514  // quadrature point
515  if (error_norm.type(var) == H1 ||
516  error_norm.type(var) == H1_SEMINORM ||
517  error_norm.type(var) == H2)
518  {
519  Gradient grad_error = grad_u_h;
520  if(_exact_deriv)
521  grad_error -= _exact_deriv(q_point[qp],parameters,sys_name,var_name);
522  else if (_exact_derivs.size() > sys_num && _exact_derivs[sys_num])
523  grad_error -= _exact_derivs[sys_num]->
524  component(var_component, q_point[qp], system.time);
525  else if(_equation_systems_fine)
526  grad_error -= fine_values->gradient(q_point[qp]);
527 
528  error_val += JxW[qp]*grad_error.size_sq();
529  }
530 
531 
532 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
533  // Compute the value of the error in the hessian at this
534  // quadrature point
535  if ((error_norm.type(var) == H2_SEMINORM ||
536  error_norm.type(var) == H2))
537  {
538  Tensor grad2_error = grad2_u_h;
539  if(_exact_hessian)
540  grad2_error -= _exact_hessian(q_point[qp],parameters,sys_name,var_name);
541  else if (_exact_hessians.size() > sys_num && _exact_hessians[sys_num])
542  grad2_error -= _exact_hessians[sys_num]->
543  component(var_component, q_point[qp], system.time);
544  else if (_equation_systems_fine)
545  grad2_error -= fine_values->hessian(q_point[qp]);
546 
547  error_val += JxW[qp]*grad2_error.size_sq();
548  }
549 #endif
550 
551  } // end qp loop
552 
553  libmesh_assert_greater_equal (error_val, 0.);
554 
555  return error_val;
556 }
557 
558 
559 
561 {
562  // delete will clean up any cloned functors and no-op on any NULL
563  // pointers
564 
565  for (unsigned int i=0; i != _exact_values.size(); ++i)
566  delete (_exact_values[i]);
567 
568  for (unsigned int i=0; i != _exact_derivs.size(); ++i)
569  delete (_exact_derivs[i]);
570 
571  for (unsigned int i=0; i != _exact_hessians.size(); ++i)
572  delete (_exact_hessians[i]);
573 }
574 
575 
576 
577 } // namespace libMesh

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

Hosted By:
SourceForge.net Logo