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