libMesh
patch_recovery_error_estimator.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
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"
27 #include "libmesh/patch_recovery_error_estimator.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"
33 #include "libmesh/libmesh_logging.h"
34 #include "libmesh/elem.h"
35 #include "libmesh/patch.h"
36 #include "libmesh/quadrature_grid.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  unsigned int npows = order+1;
62  std::vector<Real> xpow(npows,1.), ypow, zpow;
63  {
64  Real x = p(0);
65  for (unsigned int i=1; i != npows; ++i)
66  xpow[i] = xpow[i-1] * x;
67  }
68  if (dim > 1)
69  {
70  Real y = p(1);
71  ypow.resize(npows,1.);
72  for (unsigned int i=1; i != npows; ++i)
73  ypow[i] = ypow[i-1] * y;
74  }
75  if (dim > 2)
76  {
77  Real z = p(2);
78  zpow.resize(npows,1.);
79  for (unsigned int i=1; i != npows; ++i)
80  zpow[i] = zpow[i-1] * z;
81  }
82 
83  // builds psi vector of form 1 x y z x^2 xy xz y^2 yz z^2 etc..
84  // I haven't added 1D support here
85  for (unsigned int poly_deg=0; poly_deg <= static_cast<unsigned int>(order) ; poly_deg++)
86  { // loop over all polynomials of total degree = poly_deg
87 
88  switch (dim)
89  {
90  // 3D spectral polynomial basis functions
91  case 3:
92  {
93  for (int xexp=poly_deg; xexp >= 0; xexp--) // use an int for xexp since we -- it
94  for (int yexp=poly_deg-xexp; yexp >= 0; yexp--)
95  {
96  int zexp = poly_deg - xexp - yexp;
97  psi.push_back(xpow[xexp]*ypow[yexp]*zpow[zexp]);
98  }
99  break;
100  }
101 
102  // 2D spectral polynomial basis functions
103  case 2:
104  {
105  for (int xexp=poly_deg; xexp >= 0; xexp--) // use an int for xexp since we -- it
106  {
107  int yexp = poly_deg - xexp;
108  psi.push_back(xpow[xexp]*ypow[yexp]);
109  }
110  break;
111  }
112 
113  // 1D spectral polynomial basis functions
114  case 1:
115  {
116  int xexp = poly_deg;
117  psi.push_back(xpow[xexp]);
118  break;
119  }
120 
121  default:
122  libmesh_error_msg("Invalid dimension dim " << dim);
123  }
124  }
125 
126  return psi;
127 }
128 
129 
130 
132  ErrorVector & error_per_cell,
133  const NumericVector<Number> * solution_vector,
134  bool)
135 {
136  LOG_SCOPE("estimate_error()", "PatchRecoveryErrorEstimator");
137 
138  // The current mesh
139  const MeshBase & mesh = system.get_mesh();
140 
141  // Resize the error_per_cell vector to be
142  // the number of elements, initialize it to 0.
143  error_per_cell.resize (mesh.max_elem_id());
144  std::fill (error_per_cell.begin(), error_per_cell.end(), 0.);
145 
146  // Prepare current_local_solution to localize a non-standard
147  // solution vector if necessary
148  if (solution_vector && solution_vector != system.solution.get())
149  {
150  NumericVector<Number> * newsol =
151  const_cast<NumericVector<Number> *>(solution_vector);
152  System & sys = const_cast<System &>(system);
153  newsol->swap(*sys.solution);
154  sys.update();
155  }
156 
157  //------------------------------------------------------------
158  // Iterate over all the active elements in the mesh
159  // that live on this processor.
162  200),
163  EstimateError(system,
164  *this,
165  error_per_cell)
166  );
167 
168  // Each processor has now computed the error contributions
169  // for its local elements, and error_per_cell contains 0 for all the
170  // non-local elements. Summing the vector will provide the true
171  // value for each element, local or remote
172  this->reduce_error(error_per_cell, system.comm());
173 
174  // If we used a non-standard solution before, now is the time to fix
175  // the current_local_solution
176  if (solution_vector && solution_vector != system.solution.get())
177  {
178  NumericVector<Number> * newsol =
179  const_cast<NumericVector<Number> *>(solution_vector);
180  System & sys = const_cast<System &>(system);
181  newsol->swap(*sys.solution);
182  sys.update();
183  }
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
225  patch.build_around_element (elem, error_estimator.target_patch_size,
226  error_estimator.patch_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 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
241 #ifdef DEBUG
242  bool is_valid_norm_type =
243  error_estimator.error_norm.type(var) == L2 ||
244  error_estimator.error_norm.type(var) == H1_SEMINORM ||
245  error_estimator.error_norm.type(var) == H2_SEMINORM ||
246  error_estimator.error_norm.type(var) == H1_X_SEMINORM ||
247  error_estimator.error_norm.type(var) == H1_Y_SEMINORM ||
248  error_estimator.error_norm.type(var) == H1_Z_SEMINORM ||
249  error_estimator.error_norm.type(var) == L_INF ||
250  error_estimator.error_norm.type(var) == W1_INF_SEMINORM ||
251  error_estimator.error_norm.type(var) == W2_INF_SEMINORM;
252  libmesh_assert (is_valid_norm_type);
253 #endif // DEBUG
254 #else
255  libmesh_assert (error_estimator.error_norm.type(var) == L2 ||
256  error_estimator.error_norm.type(var) == L_INF ||
257  error_estimator.error_norm.type(var) == H1_SEMINORM ||
258  error_estimator.error_norm.type(var) == H1_X_SEMINORM ||
259  error_estimator.error_norm.type(var) == H1_Y_SEMINORM ||
260  error_estimator.error_norm.type(var) == H1_Z_SEMINORM ||
261  error_estimator.error_norm.type(var) == W1_INF_SEMINORM);
262 #endif
263 
264 
265 #ifdef DEBUG
266  if (var > 0)
267  {
268  // We can't mix L_inf and L_2 norms
269  bool is_valid_norm_type =
270  ((error_estimator.error_norm.type(var) == L2 ||
271  error_estimator.error_norm.type(var) == H1_SEMINORM ||
272  error_estimator.error_norm.type(var) == H1_X_SEMINORM ||
273  error_estimator.error_norm.type(var) == H1_Y_SEMINORM ||
274  error_estimator.error_norm.type(var) == H1_Z_SEMINORM ||
275  error_estimator.error_norm.type(var) == H2_SEMINORM) &&
276  (error_estimator.error_norm.type(var-1) == L2 ||
277  error_estimator.error_norm.type(var-1) == H1_SEMINORM ||
278  error_estimator.error_norm.type(var-1) == H1_X_SEMINORM ||
279  error_estimator.error_norm.type(var-1) == H1_Y_SEMINORM ||
280  error_estimator.error_norm.type(var-1) == H1_Z_SEMINORM ||
281  error_estimator.error_norm.type(var-1) == H2_SEMINORM)) ||
282  ((error_estimator.error_norm.type(var) == L_INF ||
283  error_estimator.error_norm.type(var) == W1_INF_SEMINORM ||
284  error_estimator.error_norm.type(var) == W2_INF_SEMINORM) &&
285  (error_estimator.error_norm.type(var-1) == L_INF ||
286  error_estimator.error_norm.type(var-1) == W1_INF_SEMINORM ||
287  error_estimator.error_norm.type(var-1) == W2_INF_SEMINORM));
288  libmesh_assert (is_valid_norm_type);
289  }
290 #endif // DEBUG
291 
292  // Possibly skip this variable
293  if (error_estimator.error_norm.weight(var) == 0.0) continue;
294 
295  // The type of finite element to use for this variable
296  const FEType & fe_type = dof_map.variable_type (var);
297 
298  const Order element_order = static_cast<Order>
299  (fe_type.order + elem->p_level());
300 
301  // Finite element object for use in this patch
302  UniquePtr<FEBase> fe (FEBase::build (dim, fe_type));
303 
304  // Build an appropriate Gaussian quadrature rule
305  UniquePtr<QBase> qrule (fe_type.default_quadrature_rule(dim));
306 
307  // Tell the finite element about the quadrature rule.
308  fe->attach_quadrature_rule (qrule.get());
309 
310  // Get Jacobian values, etc..
311  const std::vector<Real> & JxW = fe->get_JxW();
312  const std::vector<Point> & q_point = fe->get_xyz();
313 
314  // Get whatever phi/dphi/d2phi values we need. Avoid
315  // getting them unless the requested norm is actually going
316  // to use them.
317 
318  const std::vector<std::vector<Real>> * phi = libmesh_nullptr;
319  // If we're using phi to assert the correct dof_indices
320  // vector size later, then we'll need to get_phi whether we
321  // plan to use it or not.
322 #ifdef NDEBUG
323  if (error_estimator.error_norm.type(var) == L2 ||
324  error_estimator.error_norm.type(var) == L_INF)
325 #endif
326  phi = &(fe->get_phi());
327 
328  const std::vector<std::vector<RealGradient>> * dphi = libmesh_nullptr;
329  if (error_estimator.error_norm.type(var) == H1_SEMINORM ||
330  error_estimator.error_norm.type(var) == H1_X_SEMINORM ||
331  error_estimator.error_norm.type(var) == H1_Y_SEMINORM ||
332  error_estimator.error_norm.type(var) == H1_Z_SEMINORM ||
333  error_estimator.error_norm.type(var) == W1_INF_SEMINORM)
334  dphi = &(fe->get_dphi());
335 
336 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
337  const std::vector<std::vector<RealTensor>> * d2phi = libmesh_nullptr;
338  if (error_estimator.error_norm.type(var) == H2_SEMINORM ||
339  error_estimator.error_norm.type(var) == W2_INF_SEMINORM)
340  d2phi = &(fe->get_d2phi());
341 #endif
342 
343  // global DOF indices
344  std::vector<dof_id_type> dof_indices;
345 
346  // Compute the appropriate size for the patch projection matrices
347  // and vectors;
348  unsigned int matsize = element_order + 1;
349  if (dim > 1)
350  {
351  matsize *= (element_order + 2);
352  matsize /= 2;
353  }
354  if (dim > 2)
355  {
356  matsize *= (element_order + 3);
357  matsize /= 3;
358  }
359 
360  DenseMatrix<Number> Kp(matsize,matsize);
361  DenseVector<Number> F, Fx, Fy, Fz, Fxy, Fxz, Fyz;
362  DenseVector<Number> Pu_h, Pu_x_h, Pu_y_h, Pu_z_h, Pu_xy_h, Pu_xz_h, Pu_yz_h;
363  if (error_estimator.error_norm.type(var) == L2 ||
364  error_estimator.error_norm.type(var) == L_INF)
365  {
366  F.resize(matsize); Pu_h.resize(matsize);
367  }
368  else if (error_estimator.error_norm.type(var) == H1_SEMINORM ||
369  error_estimator.error_norm.type(var) == W1_INF_SEMINORM ||
370  error_estimator.error_norm.type(var) == H2_SEMINORM ||
371  error_estimator.error_norm.type(var) == W2_INF_SEMINORM)
372  {
373  Fx.resize(matsize); Pu_x_h.resize(matsize); // stores xx in W2 cases
374 #if LIBMESH_DIM > 1
375  Fy.resize(matsize); Pu_y_h.resize(matsize); // stores yy in W2 cases
376 #endif
377 #if LIBMESH_DIM > 2
378  Fz.resize(matsize); Pu_z_h.resize(matsize); // stores zz in W2 cases
379 #endif
380  }
381  else if (error_estimator.error_norm.type(var) == H1_X_SEMINORM)
382  {
383  Fx.resize(matsize); Pu_x_h.resize(matsize); // Only need to compute the x gradient for the x component seminorm
384  }
385  else if (error_estimator.error_norm.type(var) == H1_Y_SEMINORM)
386  {
387  libmesh_assert_greater (LIBMESH_DIM, 1);
388  Fy.resize(matsize); Pu_y_h.resize(matsize); // Only need to compute the y gradient for the y component seminorm
389  }
390  else if (error_estimator.error_norm.type(var) == H1_Z_SEMINORM)
391  {
392  libmesh_assert_greater (LIBMESH_DIM, 2);
393  Fz.resize(matsize); Pu_z_h.resize(matsize); // Only need to compute the z gradient for the z component seminorm
394  }
395 
396 #if LIBMESH_DIM > 1
397  if (error_estimator.error_norm.type(var) == H2_SEMINORM ||
398  error_estimator.error_norm.type(var) == W2_INF_SEMINORM)
399  {
400  Fxy.resize(matsize); Pu_xy_h.resize(matsize);
401 #if LIBMESH_DIM > 2
402  Fxz.resize(matsize); Pu_xz_h.resize(matsize);
403  Fyz.resize(matsize); Pu_yz_h.resize(matsize);
404 #endif
405  }
406 #endif
407 
408  //------------------------------------------------------
409  // Loop over each element in the patch and compute their
410  // contribution to the patch gradient projection.
411  Patch::const_iterator patch_it = patch.begin();
412  const Patch::const_iterator patch_end = patch.end();
413 
414  for (; patch_it != patch_end; ++patch_it)
415  {
416  // The pth element in the patch
417  const Elem * e_p = *patch_it;
418 
419  // Reinitialize the finite element data for this element
420  fe->reinit (e_p);
421 
422  // Get the global DOF indices for the current variable
423  // in the current element
424  dof_map.dof_indices (e_p, dof_indices, var);
425  libmesh_assert_equal_to (dof_indices.size(), phi->size());
426 
427  const unsigned int n_dofs =
428  cast_int<unsigned int>(dof_indices.size());
429  const unsigned int n_qp = qrule->n_points();
430 
431  // Compute the projection components from this cell.
432  // \int_{Omega_e} \psi_i \psi_j = \int_{Omega_e} du_h/dx_k \psi_i
433  for (unsigned int qp=0; qp<n_qp; qp++)
434  {
435  // Construct the shape function values for the patch projection
436  std::vector<Real> psi(specpoly(dim, element_order, q_point[qp], matsize));
437 
438  // Patch matrix contribution
439  for (unsigned int i=0; i<Kp.m(); i++)
440  for (unsigned int j=0; j<Kp.n(); j++)
441  Kp(i,j) += JxW[qp]*psi[i]*psi[j];
442 
443  if (error_estimator.error_norm.type(var) == L2 ||
444  error_estimator.error_norm.type(var) == L_INF)
445  {
446  // Compute the solution on the current patch element
447  // the quadrature point
448  Number u_h = libMesh::zero;
449 
450  for (unsigned int i=0; i<n_dofs; i++)
451  u_h += (*phi)[i][qp]*system.current_solution (dof_indices[i]);
452 
453  // Patch RHS contributions
454  for (std::size_t i=0; i<psi.size(); i++)
455  F(i) += JxW[qp]*u_h*psi[i];
456 
457  }
458  else if (error_estimator.error_norm.type(var) == H1_SEMINORM ||
459  error_estimator.error_norm.type(var) == W1_INF_SEMINORM)
460  {
461  // Compute the gradient on the current patch element
462  // at the quadrature point
463  Gradient grad_u_h;
464 
465  for (unsigned int i=0; i<n_dofs; i++)
466  grad_u_h.add_scaled ((*dphi)[i][qp],
467  system.current_solution(dof_indices[i]));
468 
469  // Patch RHS contributions
470  for (std::size_t i=0; i<psi.size(); i++)
471  {
472  Fx(i) += JxW[qp]*grad_u_h(0)*psi[i];
473 #if LIBMESH_DIM > 1
474  Fy(i) += JxW[qp]*grad_u_h(1)*psi[i];
475 #endif
476 #if LIBMESH_DIM > 2
477  Fz(i) += JxW[qp]*grad_u_h(2)*psi[i];
478 #endif
479  }
480  }
481  else if (error_estimator.error_norm.type(var) == H1_X_SEMINORM)
482  {
483  // Compute the gradient on the current patch element
484  // at the quadrature point
485  Gradient grad_u_h;
486 
487  for (unsigned int i=0; i<n_dofs; i++)
488  grad_u_h.add_scaled ((*dphi)[i][qp],
489  system.current_solution(dof_indices[i]));
490 
491  // Patch RHS contributions
492  for (std::size_t i=0; i<psi.size(); i++)
493  {
494  Fx(i) += JxW[qp]*grad_u_h(0)*psi[i];
495  }
496  }
497  else if (error_estimator.error_norm.type(var) == H1_Y_SEMINORM)
498  {
499  // Compute the gradient on the current patch element
500  // at the quadrature point
501  Gradient grad_u_h;
502 
503  for (unsigned int i=0; i<n_dofs; i++)
504  grad_u_h.add_scaled ((*dphi)[i][qp],
505  system.current_solution(dof_indices[i]));
506 
507  // Patch RHS contributions
508  for (std::size_t i=0; i<psi.size(); i++)
509  {
510  Fy(i) += JxW[qp]*grad_u_h(1)*psi[i];
511  }
512  }
513  else if (error_estimator.error_norm.type(var) == H1_Z_SEMINORM)
514  {
515  // Compute the gradient on the current patch element
516  // at the quadrature point
517  Gradient grad_u_h;
518 
519  for (unsigned int i=0; i<n_dofs; i++)
520  grad_u_h.add_scaled ((*dphi)[i][qp],
521  system.current_solution(dof_indices[i]));
522 
523  // Patch RHS contributions
524  for (std::size_t i=0; i<psi.size(); i++)
525  {
526  Fz(i) += JxW[qp]*grad_u_h(2)*psi[i];
527  }
528  }
529  else if (error_estimator.error_norm.type(var) == H2_SEMINORM ||
530  error_estimator.error_norm.type(var) == W2_INF_SEMINORM)
531  {
532 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
533  // Compute the hessian on the current patch element
534  // at the quadrature point
535  Tensor hess_u_h;
536 
537  for (unsigned int i=0; i<n_dofs; i++)
538  hess_u_h.add_scaled ((*d2phi)[i][qp],
539  system.current_solution(dof_indices[i]));
540 
541  // Patch RHS contributions
542  for (std::size_t i=0; i<psi.size(); i++)
543  {
544  Fx(i) += JxW[qp]*hess_u_h(0,0)*psi[i];
545 #if LIBMESH_DIM > 1
546  Fy(i) += JxW[qp]*hess_u_h(1,1)*psi[i];
547  Fxy(i) += JxW[qp]*hess_u_h(0,1)*psi[i];
548 #endif
549 #if LIBMESH_DIM > 2
550  Fz(i) += JxW[qp]*hess_u_h(2,2)*psi[i];
551  Fxz(i) += JxW[qp]*hess_u_h(0,2)*psi[i];
552  Fyz(i) += JxW[qp]*hess_u_h(1,2)*psi[i];
553 #endif
554  }
555 #else
556  libmesh_error_msg("ERROR: --enable-second-derivatives is required \nfor _sobolev_order == 2!");
557 #endif
558  }
559  else
560  libmesh_error_msg("Unsupported error norm type!");
561  } // end quadrature loop
562  } // end patch loop
563 
564 
565 
566  //--------------------------------------------------
567  // Now we have fully assembled the projection system
568  // for this patch. Project the gradient components.
569  // MAY NEED TO USE PARTIAL PIVOTING!
570  if (error_estimator.error_norm.type(var) == L2 ||
571  error_estimator.error_norm.type(var) == L_INF)
572  {
573  Kp.lu_solve(F, Pu_h);
574  }
575  else if (error_estimator.error_norm.type(var) == H1_SEMINORM ||
576  error_estimator.error_norm.type(var) == W1_INF_SEMINORM ||
577  error_estimator.error_norm.type(var) == H2_SEMINORM ||
578  error_estimator.error_norm.type(var) == W2_INF_SEMINORM)
579  {
580  Kp.lu_solve (Fx, Pu_x_h);
581 #if LIBMESH_DIM > 1
582  Kp.lu_solve (Fy, Pu_y_h);
583 #endif
584 #if LIBMESH_DIM > 2
585  Kp.lu_solve (Fz, Pu_z_h);
586 #endif
587  }
588  else if (error_estimator.error_norm.type(var) == H1_X_SEMINORM)
589  {
590  Kp.lu_solve (Fx, Pu_x_h);
591  }
592  else if (error_estimator.error_norm.type(var) == H1_Y_SEMINORM)
593  {
594  Kp.lu_solve (Fy, Pu_y_h);
595  }
596  else if (error_estimator.error_norm.type(var) == H1_Z_SEMINORM)
597  {
598  Kp.lu_solve (Fz, Pu_z_h);
599  }
600 
601 #if LIBMESH_DIM > 1
602  if (error_estimator.error_norm.type(var) == H2_SEMINORM ||
603  error_estimator.error_norm.type(var) == W2_INF_SEMINORM)
604  {
605  Kp.lu_solve(Fxy, Pu_xy_h);
606 #if LIBMESH_DIM > 2
607  Kp.lu_solve(Fxz, Pu_xz_h);
608  Kp.lu_solve(Fyz, Pu_yz_h);
609 #endif
610  }
611 #endif
612 
613  // If we are reusing patches, reuse the current patch to loop
614  // over all elements in the current patch, otherwise build a new
615  // patch containing just the current element and loop over it
616  // Note that C++ will not allow patch_re_end to be a const here
617  Patch::const_iterator patch_re_it;
618  Patch::const_iterator patch_re_end;
619 
620  // Declare a new patch
621  Patch patch_re(mesh.processor_id());
622 
623  if (this->error_estimator.patch_reuse)
624  {
625  // Just get the iterators from the current patch
626  patch_re_it = patch.begin();
627  patch_re_end = patch.end();
628  }
629  else
630  {
631  // Use a target patch size of just 0, this will contain
632  // just the current element
633  patch_re.build_around_element (elem, 0,
634  error_estimator.patch_growth_strategy);
635 
636  // Get the iterators from this newly constructed patch
637  patch_re_it = patch_re.begin();
638  patch_re_end = patch_re.end();
639  }
640 
641  // If we are reusing patches, loop over all the elements
642  // in the current patch and develop an estimate
643  // for all the elements by computing ||P u_h - u_h|| or ||P grad_u_h - grad_u_h||
644  // or ||P hess_u_h - hess_u_h|| according to the requested
645  // seminorm, otherwise just compute it for the current element
646 
647  // Loop over every element in the patch
648  for (unsigned int e = 0 ; patch_re_it != patch_re_end; ++patch_re_it, ++e)
649  {
650  // Build the Finite Element for the current element
651 
652  // The pth element in the patch
653  const Elem * e_p = *patch_re_it;
654 
655  // We'll need an index into the error vector for this element
656  const dof_id_type e_p_id = e_p->id();
657 
658  // We will update the new_error_per_cell vector with element_error if the
659  // error_per_cell[e_p_id] entry is non-zero, otherwise update it
660  // with 0. i.e. leave it unchanged
661 
662  // No need to compute the estimate if we are reusing patches and already have one
663  if (this->error_estimator.patch_reuse && error_per_cell[e_p_id] != 0.)
664  continue;
665 
666  // Reinitialize the finite element data for this element
667  fe->reinit (e_p);
668 
669  // Get the global DOF indices for the current variable
670  // in the current element
671  dof_map.dof_indices (e_p, dof_indices, var);
672  libmesh_assert_equal_to (dof_indices.size(), phi->size());
673 
674  // The number of dofs for this variable on this element
675  const unsigned int n_dofs =
676  cast_int<unsigned int>(dof_indices.size());
677 
678  // Variable to hold the error on the current element
679  Real element_error = 0;
680 
681  const Order qorder =
682  static_cast<Order>(fe_type.order + e_p->p_level());
683 
684  // A quadrature rule for this element
685  QGrid samprule (dim, qorder);
686 
687  if (error_estimator.error_norm.type(var) == W1_INF_SEMINORM ||
688  error_estimator.error_norm.type(var) == W2_INF_SEMINORM)
689  fe->attach_quadrature_rule (&samprule);
690 
691  // The number of points we will sample over
692  const unsigned int n_sp =
693  cast_int<unsigned int>(JxW.size());
694 
695  // Loop over every sample point for the current element
696  for (unsigned int sp=0; sp<n_sp; sp++)
697  {
698  // Compute the solution at the current sample point
699 
700  std::vector<Number> temperr(6,0.0); // x,y,z or xx,yy,zz,xy,xz,yz
701 
702  if (error_estimator.error_norm.type(var) == L2 ||
703  error_estimator.error_norm.type(var) == L_INF)
704  {
705  // Compute the value at the current sample point
706  Number u_h = libMesh::zero;
707 
708  for (unsigned int i=0; i<n_dofs; i++)
709  u_h += (*phi)[i][sp]*system.current_solution (dof_indices[i]);
710 
711  // Compute the phi values at the current sample point
712  std::vector<Real> psi(specpoly(dim, element_order, q_point[sp], matsize));
713  for (unsigned int i=0; i<matsize; i++)
714  {
715  temperr[0] += psi[i]*Pu_h(i);
716  }
717 
718  temperr[0] -= u_h;
719  }
720  else if (error_estimator.error_norm.type(var) == H1_SEMINORM ||
721  error_estimator.error_norm.type(var) == W1_INF_SEMINORM)
722  {
723  // Compute the gradient at the current sample point
724  Gradient grad_u_h;
725 
726  for (unsigned int i=0; i<n_dofs; i++)
727  grad_u_h.add_scaled ((*dphi)[i][sp],
728  system.current_solution(dof_indices[i]));
729 
730  // Compute the phi values at the current sample point
731  std::vector<Real> psi(specpoly(dim, element_order, q_point[sp], matsize));
732 
733  for (unsigned int i=0; i<matsize; i++)
734  {
735  temperr[0] += psi[i]*Pu_x_h(i);
736 #if LIBMESH_DIM > 1
737  temperr[1] += psi[i]*Pu_y_h(i);
738 #endif
739 #if LIBMESH_DIM > 2
740  temperr[2] += psi[i]*Pu_z_h(i);
741 #endif
742  }
743  temperr[0] -= grad_u_h(0);
744 #if LIBMESH_DIM > 1
745  temperr[1] -= grad_u_h(1);
746 #endif
747 #if LIBMESH_DIM > 2
748  temperr[2] -= grad_u_h(2);
749 #endif
750  }
751  else if (error_estimator.error_norm.type(var) == H1_X_SEMINORM)
752  {
753  // Compute the gradient at the current sample point
754  Gradient grad_u_h;
755 
756  for (unsigned int i=0; i<n_dofs; i++)
757  grad_u_h.add_scaled ((*dphi)[i][sp],
758  system.current_solution(dof_indices[i]));
759 
760  // Compute the phi values at the current sample point
761  std::vector<Real> psi(specpoly(dim, element_order, q_point[sp], matsize));
762  for (unsigned int i=0; i<matsize; i++)
763  {
764  temperr[0] += psi[i]*Pu_x_h(i);
765  }
766 
767  temperr[0] -= grad_u_h(0);
768  }
769  else if (error_estimator.error_norm.type(var) == H1_Y_SEMINORM)
770  {
771  // Compute the gradient at the current sample point
772  Gradient grad_u_h;
773 
774  for (unsigned int i=0; i<n_dofs; i++)
775  grad_u_h.add_scaled ((*dphi)[i][sp],
776  system.current_solution(dof_indices[i]));
777 
778  // Compute the phi values at the current sample point
779  std::vector<Real> psi(specpoly(dim, element_order, q_point[sp], matsize));
780  for (unsigned int i=0; i<matsize; i++)
781  {
782  temperr[1] += psi[i]*Pu_y_h(i);
783  }
784 
785  temperr[1] -= grad_u_h(1);
786  }
787  else if (error_estimator.error_norm.type(var) == H1_Z_SEMINORM)
788  {
789  // Compute the gradient at the current sample point
790  Gradient grad_u_h;
791 
792  for (unsigned int i=0; i<n_dofs; i++)
793  grad_u_h.add_scaled ((*dphi)[i][sp],
794  system.current_solution(dof_indices[i]));
795 
796  // Compute the phi values at the current sample point
797  std::vector<Real> psi(specpoly(dim, element_order, q_point[sp], matsize));
798  for (unsigned int i=0; i<matsize; i++)
799  {
800  temperr[2] += psi[i]*Pu_z_h(i);
801  }
802 
803  temperr[2] -= grad_u_h(2);
804  }
805  else if (error_estimator.error_norm.type(var) == H2_SEMINORM ||
806  error_estimator.error_norm.type(var) == W2_INF_SEMINORM)
807  {
808 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
809  // Compute the Hessian at the current sample point
810  Tensor hess_u_h;
811 
812  for (unsigned int i=0; i<n_dofs; i++)
813  hess_u_h.add_scaled ((*d2phi)[i][sp],
814  system.current_solution(dof_indices[i]));
815 
816  // Compute the phi values at the current sample point
817  std::vector<Real> psi(specpoly(dim, element_order, q_point[sp], matsize));
818  for (unsigned int i=0; i<matsize; i++)
819  {
820  temperr[0] += psi[i]*Pu_x_h(i);
821 #if LIBMESH_DIM > 1
822  temperr[1] += psi[i]*Pu_y_h(i);
823  temperr[3] += psi[i]*Pu_xy_h(i);
824 #endif
825 #if LIBMESH_DIM > 2
826  temperr[2] += psi[i]*Pu_z_h(i);
827  temperr[4] += psi[i]*Pu_xz_h(i);
828  temperr[5] += psi[i]*Pu_yz_h(i);
829 #endif
830  }
831 
832  temperr[0] -= hess_u_h(0,0);
833 #if LIBMESH_DIM > 1
834  temperr[1] -= hess_u_h(1,1);
835  temperr[3] -= hess_u_h(0,1);
836 #endif
837 #if LIBMESH_DIM > 2
838  temperr[2] -= hess_u_h(2,2);
839  temperr[4] -= hess_u_h(0,2);
840  temperr[5] -= hess_u_h(1,2);
841 #endif
842 #else
843  libmesh_error_msg("ERROR: --enable-second-derivatives is required \nfor _sobolev_order == 2!");
844 #endif
845  }
846  // Add up relevant terms. We can easily optimize the
847  // LIBMESH_DIM < 3 cases a little bit with the exception
848  // of the W2 cases
849 
850  if (error_estimator.error_norm.type(var) == L_INF)
851  element_error = std::max(element_error, std::abs(temperr[0]));
852  else if (error_estimator.error_norm.type(var) == W1_INF_SEMINORM)
853  for (unsigned int i=0; i != LIBMESH_DIM; ++i)
854  element_error = std::max(element_error, std::abs(temperr[i]));
855  else if (error_estimator.error_norm.type(var) == W2_INF_SEMINORM)
856  for (unsigned int i=0; i != 6; ++i)
857  element_error = std::max(element_error, std::abs(temperr[i]));
858  else if (error_estimator.error_norm.type(var) == L2)
859  element_error += JxW[sp]*TensorTools::norm_sq(temperr[0]);
860  else if (error_estimator.error_norm.type(var) == H1_SEMINORM)
861  for (unsigned int i=0; i != LIBMESH_DIM; ++i)
862  element_error += JxW[sp]*TensorTools::norm_sq(temperr[i]);
863  else if (error_estimator.error_norm.type(var) == H1_X_SEMINORM)
864  element_error += JxW[sp]*TensorTools::norm_sq(temperr[0]);
865  else if (error_estimator.error_norm.type(var) == H1_Y_SEMINORM)
866  element_error += JxW[sp]*TensorTools::norm_sq(temperr[1]);
867  else if (error_estimator.error_norm.type(var) == H1_Z_SEMINORM)
868  element_error += JxW[sp]*TensorTools::norm_sq(temperr[2]);
869  else if (error_estimator.error_norm.type(var) == H2_SEMINORM)
870  {
871  for (unsigned int i=0; i != LIBMESH_DIM; ++i)
872  element_error += JxW[sp]*TensorTools::norm_sq(temperr[i]);
873  // Off diagonal terms enter into the Hessian norm twice
874  for (unsigned int i=3; i != 6; ++i)
875  element_error += JxW[sp]*2*TensorTools::norm_sq(temperr[i]);
876  }
877 
878  } // End loop over sample points
879 
880  if (error_estimator.error_norm.type(var) == L_INF ||
881  error_estimator.error_norm.type(var) == W1_INF_SEMINORM ||
882  error_estimator.error_norm.type(var) == W2_INF_SEMINORM)
883  new_error_per_cell[e] += error_estimator.error_norm.weight(var) * element_error;
884  else if (error_estimator.error_norm.type(var) == L2 ||
885  error_estimator.error_norm.type(var) == H1_SEMINORM ||
886  error_estimator.error_norm.type(var) == H1_X_SEMINORM ||
887  error_estimator.error_norm.type(var) == H1_Y_SEMINORM ||
888  error_estimator.error_norm.type(var) == H1_Z_SEMINORM ||
889  error_estimator.error_norm.type(var) == H2_SEMINORM)
890  new_error_per_cell[e] += error_estimator.error_norm.weight_sq(var) * element_error;
891  else
892  libmesh_error_msg("Unsupported error norm type!");
893  } // End (re) loop over patch elements
894 
895  } // end variables loop
896 
897  // Now that we have the contributions from each variable,
898  // we have take square roots of the entries we
899  // added to error_per_cell to get an error norm
900  // If we are reusing patches, once again reuse the current patch to loop
901  // over all elements in the current patch, otherwise build a new
902  // patch containing just the current element and loop over it
903  Patch::const_iterator patch_re_it;
904  Patch::const_iterator patch_re_end;
905 
906  // Build a new patch if necessary
907  Patch current_elem_patch(mesh.processor_id());
908 
909  if (this->error_estimator.patch_reuse)
910  {
911  // Just get the iterators from the current patch
912  patch_re_it = patch.begin();
913  patch_re_end = patch.end();
914  }
915  else
916  {
917  // Use a target patch size of just 0, this will contain
918  // just the current element.
919  current_elem_patch.build_around_element (elem, 0,
920  error_estimator.patch_growth_strategy);
921 
922  // Get the iterators from this newly constructed patch
923  patch_re_it = current_elem_patch.begin();
924  patch_re_end = current_elem_patch.end();
925  }
926 
927  // Loop over every element in the patch we just constructed
928  for (unsigned int i = 0 ; patch_re_it != patch_re_end; ++patch_re_it, ++i)
929  {
930  // The pth element in the patch
931  const Elem * e_p = *patch_re_it;
932 
933  // We'll need an index into the error vector
934  const dof_id_type e_p_id = e_p->id();
935 
936  // Update the error_per_cell vector for this element
937  if (error_estimator.error_norm.type(0) == L2 ||
938  error_estimator.error_norm.type(0) == H1_SEMINORM ||
939  error_estimator.error_norm.type(0) == H1_X_SEMINORM ||
940  error_estimator.error_norm.type(0) == H1_Y_SEMINORM ||
941  error_estimator.error_norm.type(0) == H1_Z_SEMINORM ||
942  error_estimator.error_norm.type(0) == H2_SEMINORM)
943  {
944  Threads::spin_mutex::scoped_lock acquire(Threads::spin_mtx);
945  if (!error_per_cell[e_p_id])
946  error_per_cell[e_p_id] =
947  static_cast<ErrorVectorReal>(std::sqrt(new_error_per_cell[i]));
948  }
949  else
950  {
951  libmesh_assert (error_estimator.error_norm.type(0) == L_INF ||
952  error_estimator.error_norm.type(0) == W1_INF_SEMINORM ||
953  error_estimator.error_norm.type(0) == W2_INF_SEMINORM);
954  Threads::spin_mutex::scoped_lock acquire(Threads::spin_mtx);
955  if (!error_per_cell[e_p_id])
956  error_per_cell[e_p_id] =
957  static_cast<ErrorVectorReal>(new_error_per_cell[i]);
958  }
959 
960  } // End loop over every element in patch
961 
962  } // end element loop
963 
964 } // End () operator definition
965 
966 } // namespace libMesh
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
double abs(double a)
unsigned int n() const
This class implements useful utility functions for a patch of elements.
Definition: patch.h:47
unsigned int p_level() const
Definition: elem.h:2422
const FEType & variable_type(const unsigned int c) const
Definition: dof_map.h:1697
void parallel_for(const Range &range, const Body &body)
Execute the provided function object in parallel on the specified range.
Definition: threads_none.h:73
void add_scaled(const TypeTensor< T2 > &, const T)
Add a scaled tensor to this tensor without creating a temporary.
Definition: type_tensor.h:777
void add_scaled(const TypeVector< T2 > &, const T)
Add a scaled value to this vector without creating a temporary.
Definition: type_vector.h:624
unsigned int dim
The ErrorVector is a specialization of the StatisticsVector for error data computed on a finite eleme...
Definition: error_vector.h:50
ImplicitSystem & sys
void resize(const unsigned int n)
Resize the vector.
Definition: dense_vector.h:350
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
MeshBase & mesh
const class libmesh_nullptr_t libmesh_nullptr
OrderWrapper order
The approximation order of the element.
Definition: fe_type.h:197
The StoredRange class defines a contiguous, divisible set of objects.
Definition: stored_range.h:52
const unsigned int n_vars
Definition: tecplot_io.C:68
The libMesh namespace provides an interface to certain functionality in the library.
const_iterator begin() const
Beginning of the range.
Definition: stored_range.h:261
const Number zero
.
Definition: libmesh.h:178
long double max(long double a, double b)
DIE A HORRIBLE DEATH HERE typedef float ErrorVectorReal
void build_around_element(const Elem *elem, const unsigned int target_patch_size=10, PMF patchtype=&Patch::add_local_face_neighbors)
Erases any elements in the current patch, then builds a new patch containing element elem by repeated...
Definition: patch.C:194
This is the MeshBase class.
Definition: mesh_base.h:68
StoredRange< MeshBase::const_element_iterator, const Elem * > ConstElemRange
Definition: elem_range.h:34
libmesh_assert(j)
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:167
virtual dof_id_type max_elem_id() const =0
const MeshBase & get_mesh() const
Definition: system.h:2014
spin_mutex spin_mtx
A convenient spin mutex object which can be used for obtaining locks.
Definition: threads.C:29
PetscErrorCode Vec x
unsigned int m() const
virtual element_iterator active_local_elements_begin()=0
This is the base class for classes which contain information related to any physical process that mig...
Definition: system.h:76
UniquePtr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1523
virtual void estimate_error(const System &system, ErrorVector &error_per_cell, const NumericVector< Number > *solution_vector=libmesh_nullptr, bool estimate_parent_error=false) libmesh_override
This function uses the Patch Recovery error estimate to estimate the error on each cell...
const_iterator end() const
End of the range.
Definition: stored_range.h:266
static std::vector< Real > specpoly(const unsigned int dim, const Order order, const Point p, const unsigned int matsize)
void reduce_error(std::vector< ErrorVectorReal > &error_per_cell, const Parallel::Communicator &comm LIBMESH_CAN_DEFAULT_TO_COMMWORLD) const
This method takes the local error contributions in error_per_cell from each processor and combines th...
std::vector< object_type >::const_iterator const_iterator
Allows an StoredRange to behave like an STL container.
Definition: stored_range.h:58
void lu_solve(const DenseVector< T > &b, DenseVector< T > &x)
Solve the system Ax=b given the input vector b.
Definition: dense_matrix.C:589
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual void swap(NumericVector< T > &v)
Swaps the contents of this with v.
UniquePtr< QBase > default_quadrature_rule(const unsigned int dim, const int extraorder=0) const
Definition: fe_type.C:30
const Parallel::Communicator & comm() const
unsigned int mesh_dimension() const
Definition: mesh_base.C:148
Order
defines an enum for polynomial orders.
Definition: enum_order.h:32
dof_id_type id() const
Definition: dof_object.h:632
virtual element_iterator active_local_elements_end()=0
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
This class defines a tensor in LIBMESH_DIM dimensional Real or Complex space.
processor_id_type processor_id() const
uint8_t dof_id_type
Definition: id_types.h:64
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Fills the vector di with the global degree of freedom indices for the element.
Definition: dof_map.C:1917
static UniquePtr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.