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