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