patch_recovery_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 std::sqrt std::pow std::abs
23 
24 
25 // Local Includes
26 #include "libmesh/libmesh_common.h"
28 #include "libmesh/dof_map.h"
29 #include "libmesh/fe_base.h"
30 #include "libmesh/dense_matrix.h"
31 #include "libmesh/dense_vector.h"
32 #include "libmesh/error_vector.h"
34 #include "libmesh/elem.h"
35 #include "libmesh/patch.h"
37 #include "libmesh/system.h"
38 #include "libmesh/mesh_base.h"
39 #include "libmesh/numeric_vector.h"
40 #include "libmesh/tensor_value.h"
41 #include "libmesh/threads.h"
42 #include "libmesh/tensor_tools.h"
43 
44 namespace libMesh
45 {
46  // Setter function for the patch_reuse flag
48  {
49  patch_reuse = patch_reuse_flag;
50  }
51 
52 //-----------------------------------------------------------------
53 // PatchRecoveryErrorEstimator implementations
54 std::vector<Real> PatchRecoveryErrorEstimator::specpoly(const unsigned int dim,
55  const Order order,
56  const Point p,
57  const unsigned int matsize)
58 {
59  std::vector<Real> psi;
60  psi.reserve(matsize);
61  Real x = p(0), y=0., z=0.;
62  unsigned int npows = order+1;
63  std::vector<Real> xpow(npows,1.), ypow, zpow;
64  for (unsigned int i=1; i != npows; ++i)
65  xpow[i] = xpow[i-1] * x;
66  if (dim > 1)
67  {
68  y = p(1);
69  ypow.resize(npows,1.);
70  for (unsigned int i=1; i != npows; ++i)
71  ypow[i] = ypow[i-1] * y;
72  }
73  if (dim > 2)
74  {
75  z = p(2);
76  zpow.resize(npows,1.);
77  for (unsigned int i=1; i != npows; ++i)
78  zpow[i] = zpow[i-1] * z;
79  }
80 
81  // builds psi vector of form 1 x y z x^2 xy xz y^2 yz z^2 etc..
82  // I haven't added 1D support here
83  for (unsigned int poly_deg=0; poly_deg <= static_cast<unsigned int>(order) ; poly_deg++)
84  { // loop over all polynomials of total degreee = poly_deg
85 
86  switch (dim)
87  {
88  // 3D spectral polynomial basis functions
89  case 3:
90  {
91  for (int xexp=poly_deg; xexp >= 0; xexp--) // use an int for xexp since we -- it
92  for (int yexp=poly_deg-xexp; yexp >= 0; yexp--)
93  {
94  int zexp = poly_deg - xexp - yexp;
95  psi.push_back(xpow[xexp]*ypow[yexp]*zpow[zexp]);
96  }
97  break;
98  }
99 
100  // 2D spectral polynomial basis functions
101  case 2:
102  {
103  for (int xexp=poly_deg; xexp >= 0; xexp--) // use an int for xexp since we -- it
104  {
105  int yexp = poly_deg - xexp;
106  psi.push_back(xpow[xexp]*ypow[yexp]);
107  }
108  break;
109  }
110 
111  // 1D spectral polynomial basis functions
112  case 1:
113  {
114  int xexp = poly_deg;
115  psi.push_back(xpow[xexp]);
116  break;
117  }
118 
119  default:
120  libmesh_error();
121  }
122  }
123 
124  return psi;
125 }
126 
127 
128 
130  ErrorVector& error_per_cell,
131  const NumericVector<Number>* solution_vector,
132  bool)
133 {
134  START_LOG("estimate_error()", "PatchRecoveryErrorEstimator");
135 
136  // The current mesh
137  const MeshBase& mesh = system.get_mesh();
138 
139  // Resize the error_per_cell vector to be
140  // the number of elements, initialize it to 0.
141  error_per_cell.resize (mesh.max_elem_id());
142  std::fill (error_per_cell.begin(), error_per_cell.end(), 0.);
143 
144  // Prepare current_local_solution to localize a non-standard
145  // solution vector if necessary
146  if (solution_vector && solution_vector != system.solution.get())
147  {
148  NumericVector<Number>* newsol =
149  const_cast<NumericVector<Number>*>(solution_vector);
150  System &sys = const_cast<System&>(system);
151  newsol->swap(*sys.solution);
152  sys.update();
153  }
154 
155  //------------------------------------------------------------
156  // Iterate over all the active elements in the mesh
157  // that live on this processor.
160  200),
161  EstimateError(system,
162  *this,
163  error_per_cell)
164  );
165 
166  // Each processor has now computed the error contributions
167  // for its local elements, and error_per_cell contains 0 for all the
168  // non-local elements. Summing the vector will provide the true
169  // value for each element, local or remote
170  this->reduce_error(error_per_cell, system.comm());
171 
172  // If we used a non-standard solution before, now is the time to fix
173  // the current_local_solution
174  if (solution_vector && solution_vector != system.solution.get())
175  {
176  NumericVector<Number>* newsol =
177  const_cast<NumericVector<Number>*>(solution_vector);
178  System &sys = const_cast<System&>(system);
179  newsol->swap(*sys.solution);
180  sys.update();
181  }
182 
183  STOP_LOG("estimate_error()", "PatchRecoveryErrorEstimator");
184 }
185 
186 
187 
189 {
190  // The current mesh
191  const MeshBase& mesh = system.get_mesh();
192 
193  // The dimensionality of the mesh
194  const unsigned int dim = mesh.mesh_dimension();
195 
196  // The number of variables in the system
197  const unsigned int n_vars = system.n_vars();
198 
199  // The DofMap for this system
200  const DofMap& dof_map = system.get_dof_map();
201 
202  //------------------------------------------------------------
203  // Iterate over all the elements in the range.
204  for (ConstElemRange::const_iterator elem_it=range.begin(); elem_it!=range.end(); ++elem_it)
205  {
206  // elem is necessarily an active element on the local processor
207  const Elem* elem = *elem_it;
208 
209  // We'll need an index into the error vector
210  const dof_id_type e_id=elem->id();
211 
212  // We are going to build a patch containing the current element
213  // and its neighbors on the local processor
214  Patch patch(mesh.processor_id());
215 
216  // If we are reusing patches and the current element
217  // already has an estimate associated with it, move on the
218  // next element
219  if(this->error_estimator.patch_reuse && error_per_cell[e_id] != 0)
220  continue;
221 
222  // If we are not reusing patches or havent built one containing this element, we build one
223 
224  // Use user specified patch size and growth strategy
227 
228  // Declare a new_error_per_cell vector to hold error estimates
229  // from each element in this patch, or one estimate if we are
230  // not reusing patches since we will only be computing error for
231  // one cell
232  std::vector<Real> new_error_per_cell(1, 0.);
233  if(this->error_estimator.patch_reuse)
234  new_error_per_cell.resize(patch.size(), 0.);
235 
236  //------------------------------------------------------------
237  // Process each variable in the system using the current patch
238  for (unsigned int var=0; var<n_vars; var++)
239  {
240 #ifndef DEBUG
241  #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
251  #else
259  #endif
260  if (var > 0)
261  // We can't mix L_inf and L_2 norms
268  (error_estimator.error_norm.type(var-1) == L2 ||
274  ((error_estimator.error_norm.type(var) == L_INF ||
277  (error_estimator.error_norm.type(var-1) == L_INF ||
280 #endif
281 
282  // Possibly skip this variable
283  if (error_estimator.error_norm.weight(var) == 0.0) continue;
284 
285  // The type of finite element to use for this variable
286  const FEType& fe_type = dof_map.variable_type (var);
287 
288  const Order element_order = static_cast<Order>
289  (fe_type.order + elem->p_level());
290 
291  // Finite element object for use in this patch
292  AutoPtr<FEBase> fe (FEBase::build (dim, fe_type));
293 
294  // Build an appropriate Gaussian quadrature rule
295  AutoPtr<QBase> qrule (fe_type.default_quadrature_rule(dim));
296 
297  // Tell the finite element about the quadrature rule.
298  fe->attach_quadrature_rule (qrule.get());
299 
300  // Get Jacobian values, etc..
301  const std::vector<Real>& JxW = fe->get_JxW();
302  const std::vector<Point>& q_point = fe->get_xyz();
303 
304  // Get whatever phi/dphi/d2phi values we need. Avoid
305  // getting them unless the requested norm is actually going
306  // to use them.
307 
308  const std::vector<std::vector<Real> > *phi = NULL;
309  // If we're using phi to assert the correct dof_indices
310  // vector size later, then we'll need to get_phi whether we
311  // plan to use it or not.
312 #ifdef NDEBUG
313  if (error_estimator.error_norm.type(var) == L2 ||
315 #endif
316  phi = &(fe->get_phi());
317 
318  const std::vector<std::vector<RealGradient> > *dphi = NULL;
324  dphi = &(fe->get_dphi());
325 
326 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
327  const std::vector<std::vector<RealTensor> > *d2phi = NULL;
330  d2phi = &(fe->get_d2phi());
331 #endif
332 
333  // global DOF indices
334  std::vector<dof_id_type> dof_indices;
335 
336  // Compute the approprite size for the patch projection matrices
337  // and vectors;
338  unsigned int matsize = element_order + 1;
339  if (dim > 1)
340  {
341  matsize *= (element_order + 2);
342  matsize /= 2;
343  }
344  if (dim > 2)
345  {
346  matsize *= (element_order + 3);
347  matsize /= 3;
348  }
349 
350  DenseMatrix<Number> Kp(matsize,matsize);
351  DenseVector<Number> F, Fx, Fy, Fz, Fxy, Fxz, Fyz;
352  DenseVector<Number> Pu_h, Pu_x_h, Pu_y_h, Pu_z_h, Pu_xy_h, Pu_xz_h, Pu_yz_h;
353  if (error_estimator.error_norm.type(var) == L2 ||
355  {
356  F.resize(matsize); Pu_h.resize(matsize);
357  }
358  else if (error_estimator.error_norm.type(var) == H1_SEMINORM ||
362  {
363  Fx.resize(matsize); Pu_x_h.resize(matsize); // stores xx in W2 cases
364 #if LIBMESH_DIM > 1
365  Fy.resize(matsize); Pu_y_h.resize(matsize); // stores yy in W2 cases
366 #endif
367 #if LIBMESH_DIM > 2
368  Fz.resize(matsize); Pu_z_h.resize(matsize); // stores zz in W2 cases
369 #endif
370  }
371  else if (error_estimator.error_norm.type(var) == H1_X_SEMINORM)
372  {
373  Fx.resize(matsize); Pu_x_h.resize(matsize); // Only need to compute the x gradient for the x component seminorm
374  }
375  else if (error_estimator.error_norm.type(var) == H1_Y_SEMINORM)
376  {
377  libmesh_assert_greater (LIBMESH_DIM, 1);
378  Fy.resize(matsize); Pu_y_h.resize(matsize); // Only need to compute the y gradient for the y component seminorm
379  }
380  else if (error_estimator.error_norm.type(var) == H1_Z_SEMINORM)
381  {
382  libmesh_assert_greater (LIBMESH_DIM, 2);
383  Fz.resize(matsize); Pu_z_h.resize(matsize); // Only need to compute the z gradient for the z component seminorm
384  }
385 
386 #if LIBMESH_DIM > 1
389  {
390  Fxy.resize(matsize); Pu_xy_h.resize(matsize);
391 #if LIBMESH_DIM > 2
392  Fxz.resize(matsize); Pu_xz_h.resize(matsize);
393  Fyz.resize(matsize); Pu_yz_h.resize(matsize);
394 #endif
395  }
396 #endif
397 
398  //------------------------------------------------------
399  // Loop over each element in the patch and compute their
400  // contribution to the patch gradient projection.
401  Patch::const_iterator patch_it = patch.begin();
402  const Patch::const_iterator patch_end = patch.end();
403 
404  for (; patch_it != patch_end; ++patch_it)
405  {
406  // The pth element in the patch
407  const Elem* e_p = *patch_it;
408 
409  // Reinitialize the finite element data for this element
410  fe->reinit (e_p);
411 
412  // Get the global DOF indices for the current variable
413  // in the current element
414  dof_map.dof_indices (e_p, dof_indices, var);
415  libmesh_assert_equal_to (dof_indices.size(), phi->size());
416 
417  const unsigned int n_dofs =
418  libmesh_cast_int<unsigned int>(dof_indices.size());
419  const unsigned int n_qp = qrule->n_points();
420 
421  // Compute the projection components from this cell.
422  // \int_{Omega_e} \psi_i \psi_j = \int_{Omega_e} du_h/dx_k \psi_i
423  for (unsigned int qp=0; qp<n_qp; qp++)
424  {
425  // Construct the shape function values for the patch projection
426  std::vector<Real> psi(specpoly(dim, element_order, q_point[qp], matsize));
427 
428  // Patch matrix contribution
429  for (unsigned int i=0; i<Kp.m(); i++)
430  for (unsigned int j=0; j<Kp.n(); j++)
431  Kp(i,j) += JxW[qp]*psi[i]*psi[j];
432 
433  if (error_estimator.error_norm.type(var) == L2 ||
435  {
436  // Compute the solution on the current patch element
437  // the quadrature point
438  Number u_h = libMesh::zero;
439 
440  for (unsigned int i=0; i<n_dofs; i++)
441  u_h += (*phi)[i][qp]*system.current_solution (dof_indices[i]);
442 
443  // Patch RHS contributions
444  for (unsigned int i=0; i<psi.size(); i++)
445  F(i) = JxW[qp]*u_h*psi[i];
446 
447  }
448  else if (error_estimator.error_norm.type(var) == H1_SEMINORM ||
450  {
451  // Compute the gradient on the current patch element
452  // at the quadrature point
453  Gradient grad_u_h;
454 
455  for (unsigned int i=0; i<n_dofs; i++)
456  grad_u_h.add_scaled ((*dphi)[i][qp],
457  system.current_solution(dof_indices[i]));
458 
459  // Patch RHS contributions
460  for (unsigned int i=0; i<psi.size(); i++)
461  {
462  Fx(i) += JxW[qp]*grad_u_h(0)*psi[i];
463 #if LIBMESH_DIM > 1
464  Fy(i) += JxW[qp]*grad_u_h(1)*psi[i];
465 #endif
466 #if LIBMESH_DIM > 2
467  Fz(i) += JxW[qp]*grad_u_h(2)*psi[i];
468 #endif
469  }
470  }
471  else if (error_estimator.error_norm.type(var) == H1_X_SEMINORM)
472  {
473  // Compute the gradient on the current patch element
474  // at the quadrature point
475  Gradient grad_u_h;
476 
477  for (unsigned int i=0; i<n_dofs; i++)
478  grad_u_h.add_scaled ((*dphi)[i][qp],
479  system.current_solution(dof_indices[i]));
480 
481  // Patch RHS contributions
482  for (unsigned int i=0; i<psi.size(); i++)
483  {
484  Fx(i) += JxW[qp]*grad_u_h(0)*psi[i];
485  }
486  }
487  else if (error_estimator.error_norm.type(var) == H1_Y_SEMINORM)
488  {
489  // Compute the gradient on the current patch element
490  // at the quadrature point
491  Gradient grad_u_h;
492 
493  for (unsigned int i=0; i<n_dofs; i++)
494  grad_u_h.add_scaled ((*dphi)[i][qp],
495  system.current_solution(dof_indices[i]));
496 
497  // Patch RHS contributions
498  for (unsigned int i=0; i<psi.size(); i++)
499  {
500  Fy(i) += JxW[qp]*grad_u_h(1)*psi[i];
501  }
502  }
503  else if (error_estimator.error_norm.type(var) == H1_Z_SEMINORM)
504  {
505  // Compute the gradient on the current patch element
506  // at the quadrature point
507  Gradient grad_u_h;
508 
509  for (unsigned int i=0; i<n_dofs; i++)
510  grad_u_h.add_scaled ((*dphi)[i][qp],
511  system.current_solution(dof_indices[i]));
512 
513  // Patch RHS contributions
514  for (unsigned int i=0; i<psi.size(); i++)
515  {
516  Fz(i) += JxW[qp]*grad_u_h(2)*psi[i];
517  }
518  }
519  else if (error_estimator.error_norm.type(var) == H2_SEMINORM ||
521  {
522 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
523  // Compute the hessian on the current patch element
524  // at the quadrature point
525  Tensor hess_u_h;
526 
527  for (unsigned int i=0; i<n_dofs; i++)
528  hess_u_h.add_scaled ((*d2phi)[i][qp],
529  system.current_solution(dof_indices[i]));
530 
531  // Patch RHS contributions
532  for (unsigned int i=0; i<psi.size(); i++)
533  {
534  Fx(i) += JxW[qp]*hess_u_h(0,0)*psi[i];
535 #if LIBMESH_DIM > 1
536  Fy(i) += JxW[qp]*hess_u_h(1,1)*psi[i];
537  Fxy(i) += JxW[qp]*hess_u_h(0,1)*psi[i];
538 #endif
539 #if LIBMESH_DIM > 2
540  Fz(i) += JxW[qp]*hess_u_h(2,2)*psi[i];
541  Fxz(i) += JxW[qp]*hess_u_h(0,2)*psi[i];
542  Fyz(i) += JxW[qp]*hess_u_h(1,2)*psi[i];
543 #endif
544  }
545 #else
546  libMesh::err << "ERROR: --enable-second-derivatives is required\n"
547  << " for _sobolev_order == 2!\n";
548  libmesh_error();
549 #endif
550  }
551  else
552  libmesh_error();
553  } // end quadrature loop
554  } // end patch loop
555 
556 
557 
558  //--------------------------------------------------
559  // Now we have fully assembled the projection system
560  // for this patch. Project the gradient components.
561  // MAY NEED TO USE PARTIAL PIVOTING!
562  if (error_estimator.error_norm.type(var) == L2 ||
564  {
565  Kp.lu_solve(F, Pu_h);
566  }
567  else if (error_estimator.error_norm.type(var) == H1_SEMINORM ||
571  {
572  Kp.lu_solve (Fx, Pu_x_h);
573 #if LIBMESH_DIM > 1
574  Kp.lu_solve (Fy, Pu_y_h);
575 #endif
576 #if LIBMESH_DIM > 2
577  Kp.lu_solve (Fz, Pu_z_h);
578 #endif
579  }
580  else if (error_estimator.error_norm.type(var) == H1_X_SEMINORM)
581  {
582  Kp.lu_solve (Fx, Pu_x_h);
583  }
584  else if (error_estimator.error_norm.type(var) == H1_Y_SEMINORM)
585  {
586  Kp.lu_solve (Fy, Pu_y_h);
587  }
588  else if (error_estimator.error_norm.type(var) == H1_Z_SEMINORM)
589  {
590  Kp.lu_solve (Fz, Pu_z_h);
591  }
592 
593 #if LIBMESH_DIM > 1
596  {
597  Kp.lu_solve(Fxy, Pu_xy_h);
598 #if LIBMESH_DIM > 2
599  Kp.lu_solve(Fxz, Pu_xz_h);
600  Kp.lu_solve(Fyz, Pu_yz_h);
601 #endif
602  }
603 #endif
604 
605  // If we are reusing patches, reuse the current patch to loop
606  // over all elements in the current patch, otherwise build a new
607  // patch containing just the current element and loop over it
608  // Note that C++ will not allow patch_re_end to be a const here
609  Patch::const_iterator patch_re_it;
610  Patch::const_iterator patch_re_end;
611 
612  // Declare a new patch
613  Patch patch_re(mesh.processor_id());
614 
615  if(this->error_estimator.patch_reuse)
616  {
617  // Just get the iterators from the current patch
618  patch_re_it = patch.begin();
619  patch_re_end = patch.end();
620  }
621  else
622  {
623  // Use a target patch size of just 0, this will contain
624  // just the current element
625  patch_re.build_around_element (elem, 0,
627 
628  // Get the iterators from this newly constructed patch
629  patch_re_it = patch_re.begin();
630  patch_re_end = patch_re.end();
631  }
632 
633  // If we are reusing patches, loop over all the elements
634  // in the current patch and develop an estimate
635  // for all the elements by computing ||P u_h - u_h|| or ||P grad_u_h - grad_u_h||
636  // or ||P hess_u_h - hess_u_h|| according to the requested
637  // seminorm, otherwise just compute it for the current element
638 
639  // Loop over every element in the patch
640  for (unsigned int e = 0 ; patch_re_it != patch_re_end; patch_re_it++, ++e)
641  {
642  // Build the Finite Element for the current element
643 
644  // The pth element in the patch
645  const Elem* e_p = *patch_re_it;
646 
647  // We'll need an index into the error vector for this element
648  const dof_id_type e_p_id = e_p->id();
649 
650  // We will update the new_error_per_cell vector with element_error if the
651  // error_per_cell[e_p_id] entry is non-zero, otherwise update it
652  // with 0. i.e. leave it unchanged
653 
654  // No need to compute the estimate if we are reusing patches and already have one
655  if (this->error_estimator.patch_reuse && error_per_cell[e_p_id] != 0.)
656  continue;
657 
658  // Reinitialize the finite element data for this element
659  fe->reinit (e_p);
660 
661  // Get the global DOF indices for the current variable
662  // in the current element
663  dof_map.dof_indices (e_p, dof_indices, var);
664  libmesh_assert_equal_to (dof_indices.size(), phi->size());
665 
666  // The number of dofs for this variable on this element
667  const unsigned int n_dofs =
668  libmesh_cast_int<unsigned int>(dof_indices.size());
669 
670  // Variable to hold the error on the current element
671  Real element_error = 0;
672 
673  const Order qorder =
674  static_cast<Order>(fe_type.order + e_p->p_level());
675 
676  // A quadrature rule for this element
677  QGrid samprule (dim, qorder);
678 
681  fe->attach_quadrature_rule (&samprule);
682 
683  // The number of points we will sample over
684  const unsigned int n_sp =
685  libmesh_cast_int<unsigned int>(JxW.size());
686 
687  // Loop over every sample point for the current element
688  for (unsigned int sp=0; sp<n_sp; sp++)
689  {
690  // Compute the solution at the current sample point
691 
692  std::vector<Number> temperr(6,0.0); // x,y,z or xx,yy,zz,xy,xz,yz
693 
694  if (error_estimator.error_norm.type(var) == L2 ||
696  {
697  // Compute the value at the current sample point
698  Number u_h = libMesh::zero;
699 
700  for (unsigned int i=0; i<n_dofs; i++)
701  u_h += (*phi)[i][sp]*system.current_solution (dof_indices[i]);
702 
703  // Compute the phi values at the current sample point
704  std::vector<Real> psi(specpoly(dim, element_order, q_point[sp], matsize));
705  for (unsigned int i=0; i<matsize; i++)
706  {
707  temperr[0] += psi[i]*Pu_h(i);
708  }
709 
710  temperr[0] -= u_h;
711  }
712  else if (error_estimator.error_norm.type(var) == H1_SEMINORM ||
714  {
715  // Compute the gradient at the current sample point
716  Gradient grad_u_h;
717 
718  for (unsigned int i=0; i<n_dofs; i++)
719  grad_u_h.add_scaled ((*dphi)[i][sp],
720  system.current_solution(dof_indices[i]));
721 
722  // Compute the phi values at the current sample point
723  std::vector<Real> psi(specpoly(dim, element_order, q_point[sp], matsize));
724 
725  for (unsigned int i=0; i<matsize; i++)
726  {
727  temperr[0] += psi[i]*Pu_x_h(i);
728 #if LIBMESH_DIM > 1
729  temperr[1] += psi[i]*Pu_y_h(i);
730 #endif
731 #if LIBMESH_DIM > 2
732  temperr[2] += psi[i]*Pu_z_h(i);
733 #endif
734  }
735  temperr[0] -= grad_u_h(0);
736 #if LIBMESH_DIM > 1
737  temperr[1] -= grad_u_h(1);
738 #endif
739 #if LIBMESH_DIM > 2
740  temperr[2] -= grad_u_h(2);
741 #endif
742  }
743  else if (error_estimator.error_norm.type(var) == H1_X_SEMINORM)
744  {
745  // Compute the gradient at the current sample point
746  Gradient grad_u_h;
747 
748  for (unsigned int i=0; i<n_dofs; i++)
749  grad_u_h.add_scaled ((*dphi)[i][sp],
750  system.current_solution(dof_indices[i]));
751 
752  // Compute the phi values at the current sample point
753  std::vector<Real> psi(specpoly(dim, element_order, q_point[sp], matsize));
754  for (unsigned int i=0; i<matsize; i++)
755  {
756  temperr[0] += psi[i]*Pu_x_h(i);
757  }
758 
759  temperr[0] -= grad_u_h(0);
760  }
761  else if (error_estimator.error_norm.type(var) == H1_Y_SEMINORM)
762  {
763  // Compute the gradient at the current sample point
764  Gradient grad_u_h;
765 
766  for (unsigned int i=0; i<n_dofs; i++)
767  grad_u_h.add_scaled ((*dphi)[i][sp],
768  system.current_solution(dof_indices[i]));
769 
770  // Compute the phi values at the current sample point
771  std::vector<Real> psi(specpoly(dim, element_order, q_point[sp], matsize));
772  for (unsigned int i=0; i<matsize; i++)
773  {
774  temperr[1] += psi[i]*Pu_y_h(i);
775  }
776 
777  temperr[1] -= grad_u_h(1);
778  }
779  else if (error_estimator.error_norm.type(var) == H1_Z_SEMINORM)
780  {
781  // Compute the gradient at the current sample point
782  Gradient grad_u_h;
783 
784  for (unsigned int i=0; i<n_dofs; i++)
785  grad_u_h.add_scaled ((*dphi)[i][sp],
786  system.current_solution(dof_indices[i]));
787 
788  // Compute the phi values at the current sample point
789  std::vector<Real> psi(specpoly(dim, element_order, q_point[sp], matsize));
790  for (unsigned int i=0; i<matsize; i++)
791  {
792  temperr[2] += psi[i]*Pu_z_h(i);
793  }
794 
795  temperr[2] -= grad_u_h(2);
796  }
797  else if (error_estimator.error_norm.type(var) == H2_SEMINORM ||
799  {
800 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
801  // Compute the Hessian at the current sample point
802  Tensor hess_u_h;
803 
804  for (unsigned int i=0; i<n_dofs; i++)
805  hess_u_h.add_scaled ((*d2phi)[i][sp],
806  system.current_solution(dof_indices[i]));
807 
808  // Compute the phi values at the current sample point
809  std::vector<Real> psi(specpoly(dim, element_order, q_point[sp], matsize));
810  for (unsigned int i=0; i<matsize; i++)
811  {
812  temperr[0] += psi[i]*Pu_x_h(i);
813 #if LIBMESH_DIM > 1
814  temperr[1] += psi[i]*Pu_y_h(i);
815  temperr[3] += psi[i]*Pu_xy_h(i);
816 #endif
817 #if LIBMESH_DIM > 2
818  temperr[2] += psi[i]*Pu_z_h(i);
819  temperr[4] += psi[i]*Pu_xz_h(i);
820  temperr[5] += psi[i]*Pu_yz_h(i);
821 #endif
822  }
823 
824  temperr[0] -= hess_u_h(0,0);
825 #if LIBMESH_DIM > 1
826  temperr[1] -= hess_u_h(1,1);
827  temperr[3] -= hess_u_h(0,1);
828 #endif
829 #if LIBMESH_DIM > 2
830  temperr[2] -= hess_u_h(2,2);
831  temperr[4] -= hess_u_h(0,2);
832  temperr[5] -= hess_u_h(1,2);
833 #endif
834 #else
835  libMesh::err << "ERROR: --enable-second-derivatives is required\n"
836  << " for _sobolev_order == 2!\n";
837  libmesh_error();
838 #endif
839  }
840  // Add up relevant terms. We can easily optimize the
841  // LIBMESH_DIM < 3 cases a little bit with the exception
842  // of the W2 cases
843 
844  if (error_estimator.error_norm.type(var) == L_INF)
845  element_error = std::max(element_error, std::abs(temperr[0]));
847  for (unsigned int i=0; i != LIBMESH_DIM; ++i)
848  element_error = std::max(element_error, std::abs(temperr[i]));
850  for (unsigned int i=0; i != 6; ++i)
851  element_error = std::max(element_error, std::abs(temperr[i]));
852  else if (error_estimator.error_norm.type(var) == L2)
853  element_error += JxW[sp]*TensorTools::norm_sq(temperr[0]);
854  else if (error_estimator.error_norm.type(var) == H1_SEMINORM)
855  for (unsigned int i=0; i != LIBMESH_DIM; ++i)
856  element_error += JxW[sp]*TensorTools::norm_sq(temperr[i]);
857  else if (error_estimator.error_norm.type(var) == H1_X_SEMINORM)
858  element_error += JxW[sp]*TensorTools::norm_sq(temperr[0]);
859  else if (error_estimator.error_norm.type(var) == H1_Y_SEMINORM)
860  element_error += JxW[sp]*TensorTools::norm_sq(temperr[1]);
861  else if (error_estimator.error_norm.type(var) == H1_Z_SEMINORM)
862  element_error += JxW[sp]*TensorTools::norm_sq(temperr[2]);
863  else if (error_estimator.error_norm.type(var) == H2_SEMINORM)
864  {
865  for (unsigned int i=0; i != LIBMESH_DIM; ++i)
866  element_error += JxW[sp]*TensorTools::norm_sq(temperr[i]);
867  // Off diagonal terms enter into the Hessian norm twice
868  for (unsigned int i=3; i != 6; ++i)
869  element_error += JxW[sp]*2*TensorTools::norm_sq(temperr[i]);
870  }
871 
872  } // End loop over sample points
873 
874  if (error_estimator.error_norm.type(var) == L_INF ||
877  new_error_per_cell[e] += error_estimator.error_norm.weight(var) * element_error;
878  else if (error_estimator.error_norm.type(var) == L2 ||
884  new_error_per_cell[e] += error_estimator.error_norm.weight_sq(var) * element_error;
885  else
886  libmesh_error();
887  } // End (re) loop over patch elements
888 
889  } // end variables loop
890 
891  // Now that we have the contributions from each variable,
892  // we have take square roots of the entries we
893  // added to error_per_cell to get an error norm
894  // If we are reusing patches, once again reuse the current patch to loop
895  // over all elements in the current patch, otherwise build a new
896  // patch containing just the current element and loop over it
897  Patch::const_iterator patch_re_it;
898  Patch::const_iterator patch_re_end;
899 
900  // Build a new patch if necessary
901  Patch current_elem_patch(mesh.processor_id());
902 
903  if(this->error_estimator.patch_reuse)
904  {
905  // Just get the iterators from the current patch
906  patch_re_it = patch.begin();
907  patch_re_end = patch.end();
908  }
909  else
910  {
911  // Use a target patch size of just 0, this will contain
912  // just the current element.
913  current_elem_patch.build_around_element (elem, 0,
915 
916  // Get the iterators from this newly constructed patch
917  patch_re_it = current_elem_patch.begin();
918  patch_re_end = current_elem_patch.end();
919  }
920 
921  // Loop over every element in the patch we just constructed
922  for (unsigned int i = 0 ; patch_re_it != patch_re_end; ++patch_re_it, ++i)
923  {
924  // The pth element in the patch
925  const Elem* e_p = *patch_re_it;
926 
927  // We'll need an index into the error vector
928  const dof_id_type e_p_id = e_p->id();
929 
930  // Update the error_per_cell vector for this element
931  if (error_estimator.error_norm.type(0) == L2 ||
937  {
938  Threads::spin_mutex::scoped_lock acquire(Threads::spin_mtx);
939  if (!error_per_cell[e_p_id])
940  error_per_cell[e_p_id] =
941  static_cast<ErrorVectorReal>(std::sqrt(new_error_per_cell[i]));
942  }
943  else
944  {
948  Threads::spin_mutex::scoped_lock acquire(Threads::spin_mtx);
949  if (!error_per_cell[e_p_id])
950  error_per_cell[e_p_id] =
951  static_cast<ErrorVectorReal>(new_error_per_cell[i]);
952  }
953 
954  } // End loop over every element in patch
955 
956  } // end element loop
957 
958 } // End () operator definition
959 
960 } // namespace libMesh

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

Hosted By:
SourceForge.net Logo