libMesh
Public Member Functions | Private Attributes | List of all members
libMesh::WeightedPatchRecoveryErrorEstimator::EstimateError Class Reference

Class to compute the error contribution for a range of elements. More...

Public Member Functions

 EstimateError (const System &sys, const WeightedPatchRecoveryErrorEstimator &ee, ErrorVector &epc)
 
void operator() (const ConstElemRange &range) const
 

Private Attributes

const System & system
 Function to set the boolean patch_reuse in case the user wants to change the default behaviour of patch_recovery_error_estimator. More...
 
const WeightedPatchRecoveryErrorEstimator & error_estimator
 
ErrorVector & error_per_cell
 

Detailed Description

Class to compute the error contribution for a range of elements.

May be executed in parallel on separate threads.

Definition at line 93 of file weighted_patch_recovery_error_estimator.h.

Constructor & Destructor Documentation

◆ EstimateError()

libMesh::WeightedPatchRecoveryErrorEstimator::EstimateError::EstimateError ( const System &  sys,
const WeightedPatchRecoveryErrorEstimator &  ee,
ErrorVector &  epc 
)
inline

Definition at line 96 of file weighted_patch_recovery_error_estimator.h.

98  :
99  system(sys),
100  error_estimator(ee),
101  error_per_cell(epc)
102  {}

Member Function Documentation

◆ operator()()

void libMesh::WeightedPatchRecoveryErrorEstimator::EstimateError::operator() ( const ConstElemRange &  range) const

Definition at line 112 of file weighted_patch_recovery_error_estimator.C.

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  const unsigned int m = Kp.m(), n = Kp.n();
357  for (unsigned int i=0; i<m; i++)
358  for (unsigned int j=0; j<n; j++)
359  Kp(i,j) += JxW[qp]*psi[i]*psi[j];
360 
361  if (error_estimator.error_norm.type(var) == L2 ||
362  error_estimator.error_norm.type(var) == L_INF)
363  {
364  // Compute the solution on the current patch element
365  // the quadrature point
366  Number u_h = libMesh::zero;
367 
368  for (unsigned int i=0; i<n_dofs; i++)
369  u_h += (*phi)[i][qp]*system.current_solution (dof_indices[i]);
370 
371  // Patch RHS contributions
372  for (unsigned int i=0; i != psi_size; i++)
373  F(i) = JxW[qp]*u_h*psi[i];
374 
375  }
376  else if (error_estimator.error_norm.type(var) == H1_SEMINORM ||
377  error_estimator.error_norm.type(var) == W1_INF_SEMINORM)
378  {
379  // Compute the gradient on the current patch element
380  // at the quadrature point
381  Gradient grad_u_h;
382 
383  for (std::size_t i=0; i<n_dofs; i++)
384  grad_u_h.add_scaled ((*dphi)[i][qp],
385  system.current_solution(dof_indices[i]));
386 
387 
388 
389  // Patch RHS contributions
390  for (unsigned int i=0; i != psi_size; i++)
391  {
392  Fx(i) += JxW[qp]*grad_u_h(0)*psi[i];
393 #if LIBMESH_DIM > 1
394  Fy(i) += JxW[qp]*grad_u_h(1)*psi[i];
395 #endif
396 #if LIBMESH_DIM > 2
397  Fz(i) += JxW[qp]*grad_u_h(2)*psi[i];
398 #endif
399  }
400  }
401  else if (error_estimator.error_norm.type(var) == H1_X_SEMINORM)
402  {
403  // Compute the gradient on the current patch element
404  // at the quadrature point
405  Gradient grad_u_h;
406 
407  for (unsigned int i=0; i<n_dofs; i++)
408  grad_u_h.add_scaled ((*dphi)[i][qp],
409  system.current_solution(dof_indices[i]));
410 
411 
412 
413  // Patch RHS contributions
414  for (unsigned int i=0; i != psi_size; i++)
415  {
416  Fx(i) += JxW[qp]*grad_u_h(0)*psi[i];
417  }
418  }
419 #if LIBMESH_DIM > 1
420  else if (error_estimator.error_norm.type(var) == H1_Y_SEMINORM)
421  {
422  // Compute the gradient on the current patch element
423  // at the quadrature point
424  Gradient grad_u_h;
425 
426  for (unsigned int i=0; i<n_dofs; i++)
427  grad_u_h.add_scaled ((*dphi)[i][qp],
428  system.current_solution(dof_indices[i]));
429 
430 
431 
432  // Patch RHS contributions
433  for (unsigned int i=0; i != psi_size; i++)
434  {
435  Fy(i) += JxW[qp]*grad_u_h(1)*psi[i];
436  }
437  }
438 #endif // LIBMESH_DIM > 1
439 #if LIBMESH_DIM > 2
440  else if (error_estimator.error_norm.type(var) == H1_Z_SEMINORM)
441  {
442  // Compute the gradient on the current patch element
443  // at the quadrature point
444  Gradient grad_u_h;
445 
446  for (unsigned int i=0; i<n_dofs; i++)
447  grad_u_h.add_scaled ((*dphi)[i][qp],
448  system.current_solution(dof_indices[i]));
449 
450 
451 
452  // Patch RHS contributions
453  for (unsigned int i=0; i != psi_size; i++)
454  {
455  Fz(i) += JxW[qp]*grad_u_h(2)*psi[i];
456  }
457  }
458 #endif // LIBMESH_DIM > 2
459  else if (error_estimator.error_norm.type(var) == H2_SEMINORM ||
460  error_estimator.error_norm.type(var) == W2_INF_SEMINORM)
461  {
462 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
463  // Compute the hessian on the current patch element
464  // at the quadrature point
465  Tensor hess_u_h;
466 
467  for (unsigned int i=0; i<n_dofs; i++)
468  hess_u_h.add_scaled ((*d2phi)[i][qp],
469  system.current_solution(dof_indices[i]));
470 
471 
472 
473  // Patch RHS contributions
474  for (unsigned int i=0; i != psi_size; i++)
475  {
476  Fx(i) += JxW[qp]*hess_u_h(0,0)*psi[i];
477 #if LIBMESH_DIM > 1
478  Fy(i) += JxW[qp]*hess_u_h(1,1)*psi[i];
479  Fxy(i) += JxW[qp]*hess_u_h(0,1)*psi[i];
480 #endif
481 #if LIBMESH_DIM > 2
482  Fz(i) += JxW[qp]*hess_u_h(2,2)*psi[i];
483  Fxz(i) += JxW[qp]*hess_u_h(0,2)*psi[i];
484  Fyz(i) += JxW[qp]*hess_u_h(1,2)*psi[i];
485 #endif
486  }
487 #else
488  libmesh_error_msg("ERROR: --enable-second-derivatives is required \nfor _sobolev_order == 2!");
489 #endif
490  }
491  else
492  libmesh_error_msg("Unsupported error norm type!");
493  } // end quadrature loop
494  } // end patch loop
495 
496 
497 
498  //--------------------------------------------------
499  // Now we have fully assembled the projection system
500  // for this patch. Project the gradient components.
501  // MAY NEED TO USE PARTIAL PIVOTING!
502  if (error_estimator.error_norm.type(var) == L2 ||
503  error_estimator.error_norm.type(var) == L_INF)
504  {
505  Kp.lu_solve(F, Pu_h);
506  }
507  else if (error_estimator.error_norm.type(var) == H1_SEMINORM ||
508  error_estimator.error_norm.type(var) == W1_INF_SEMINORM ||
509  error_estimator.error_norm.type(var) == H2_SEMINORM ||
510  error_estimator.error_norm.type(var) == W2_INF_SEMINORM)
511  {
512  Kp.lu_solve (Fx, Pu_x_h);
513 #if LIBMESH_DIM > 1
514  Kp.lu_solve (Fy, Pu_y_h);
515 #endif
516 #if LIBMESH_DIM > 2
517  Kp.lu_solve (Fz, Pu_z_h);
518 #endif
519  }
520  else if (error_estimator.error_norm.type(var) == H1_X_SEMINORM)
521  {
522  Kp.lu_solve (Fx, Pu_x_h);
523  }
524  else if (error_estimator.error_norm.type(var) == H1_Y_SEMINORM)
525  {
526  Kp.lu_solve (Fy, Pu_y_h);
527  }
528  else if (error_estimator.error_norm.type(var) == H1_Z_SEMINORM)
529  {
530  Kp.lu_solve (Fz, Pu_z_h);
531  }
532 
533 #if LIBMESH_DIM > 1
534  if (error_estimator.error_norm.type(var) == H2_SEMINORM ||
535  error_estimator.error_norm.type(var) == W2_INF_SEMINORM)
536  {
537  Kp.lu_solve(Fxy, Pu_xy_h);
538 #if LIBMESH_DIM > 2
539  Kp.lu_solve(Fxz, Pu_xz_h);
540  Kp.lu_solve(Fyz, Pu_yz_h);
541 #endif
542  }
543 #endif
544 
545  // If we are reusing patches, reuse the current patch to loop
546  // over all elements in the current patch, otherwise build a new
547  // patch containing just the current element and loop over it
548  // Note that C++ will not allow patch_re_end to be a const here
549  Patch::const_iterator patch_re_it;
550  Patch::const_iterator patch_re_end;
551 
552  // Declare a new patch
553  Patch patch_re(mesh.processor_id());
554 
555  if (this->error_estimator.patch_reuse)
556  {
557  // Just get the iterators from the current patch
558  patch_re_it = patch.begin();
559  patch_re_end = patch.end();
560  }
561  else
562  {
563  // Use a target patch size of just 0, this will contain
564  // just the current element
565  patch_re.build_around_element (elem, 0,
566  error_estimator.patch_growth_strategy);
567 
568  // Get the iterators from this newly constructed patch
569  patch_re_it = patch_re.begin();
570  patch_re_end = patch_re.end();
571  }
572 
573  // If we are reusing patches, loop over all the elements
574  // in the current patch and develop an estimate
575  // for all the elements by computing || w * (P u_h - u_h)|| or ||w *(P grad_u_h - grad_u_h)||
576  // or ||w * (P hess_u_h - hess_u_h)|| according to the requested
577  // seminorm, otherwise just compute it for the current element
578 
579  // Get an FEMContext for this system, this will help us in
580  // obtaining the weights from the user code
581  FEMContext femcontext(system);
582  error_estimator.weight_functions[var]->init_context(femcontext);
583 
584  // Loop over every element in the patch
585  for (unsigned int e = 0 ; patch_re_it != patch_re_end; ++patch_re_it, ++e)
586  {
587  // Build the Finite Element for the current element
588 
589  // The pth element in the patch
590  const Elem * e_p = *patch_re_it;
591 
592  // We'll need an index into the error vector for this element
593  const dof_id_type e_p_id = e_p->id();
594 
595  // Initialize the FEMContext
596  femcontext.pre_fe_reinit(system, e_p);
597 
598  // We will update the new_error_per_cell vector with element_error if the
599  // error_per_cell[e_p_id] entry is non-zero, otherwise update it
600  // with 0. i.e. leave it unchanged
601 
602  // No need to compute the estimate if we are reusing patches and already have one
603  if (this->error_estimator.patch_reuse && error_per_cell[e_p_id] != 0.)
604  continue;
605 
606  // Reinitialize the finite element data for this element
607  fe->reinit (e_p);
608 
609  // Get the global DOF indices for the current variable
610  // in the current element
611  dof_map.dof_indices (e_p, dof_indices, var);
612  libmesh_assert (dof_indices.size() == phi->size());
613 
614  // The number of dofs for this variable on this element
615  const unsigned int n_dofs =
616  cast_int<unsigned int>(dof_indices.size());
617 
618  // Variable to hold the error on the current element
619  Real element_error = 0;
620 
621  const Order qorder =
622  static_cast<Order>(fe_type.order + e_p->p_level());
623 
624  // A quadrature rule for this element
625  QGrid samprule (dim, qorder);
626 
627  if (error_estimator.error_norm.type(var) == W1_INF_SEMINORM ||
628  error_estimator.error_norm.type(var) == W2_INF_SEMINORM)
629  fe->attach_quadrature_rule (&samprule);
630 
631  // The number of points we will sample over
632  const unsigned int n_sp =
633  cast_int<unsigned int>(JxW.size());
634 
635  // Loop over every sample point for the current element
636  for (unsigned int sp=0; sp<n_sp; sp++)
637  {
638  // Compute the solution at the current sample point
639 
640  std::vector<Number> temperr(6,0.0); // x,y,z or xx,yy,zz,xy,xz,yz
641 
642  if (error_estimator.error_norm.type(var) == L2 ||
643  error_estimator.error_norm.type(var) == L_INF)
644  {
645  // Compute the value at the current sample point
646  Number u_h = libMesh::zero;
647 
648  for (unsigned int i=0; i<n_dofs; i++)
649  u_h += (*phi)[i][sp]*system.current_solution (dof_indices[i]);
650 
651  // Compute the phi values at the current sample point
652  std::vector<Real> psi(specpoly(dim, element_order, q_point[sp], matsize));
653  for (unsigned int i=0; i<matsize; i++)
654  {
655  temperr[0] += psi[i]*Pu_h(i);
656  }
657 
658  temperr[0] -= u_h;
659  }
660  else if (error_estimator.error_norm.type(var) == H1_SEMINORM ||
661  error_estimator.error_norm.type(var) == W1_INF_SEMINORM)
662  {
663  // Compute the gradient at the current sample point
664  Gradient grad_u_h;
665 
666  for (unsigned int i=0; i<n_dofs; i++)
667  grad_u_h.add_scaled ((*dphi)[i][sp],
668  system.current_solution(dof_indices[i]));
669 
670  // Compute the phi values at the current sample point
671  std::vector<Real> psi(specpoly(dim, element_order, q_point[sp], matsize));
672 
673  for (unsigned int i=0; i<matsize; i++)
674  {
675  temperr[0] += psi[i]*Pu_x_h(i);
676 #if LIBMESH_DIM > 1
677  temperr[1] += psi[i]*Pu_y_h(i);
678 #endif
679 #if LIBMESH_DIM > 2
680  temperr[2] += psi[i]*Pu_z_h(i);
681 #endif
682  }
683  temperr[0] -= grad_u_h(0);
684 #if LIBMESH_DIM > 1
685  temperr[1] -= grad_u_h(1);
686 #endif
687 #if LIBMESH_DIM > 2
688  temperr[2] -= grad_u_h(2);
689 #endif
690  }
691  else if (error_estimator.error_norm.type(var) == H1_X_SEMINORM)
692  {
693  // Compute the gradient at the current sample point
694  Gradient grad_u_h;
695 
696  for (unsigned int i=0; i<n_dofs; i++)
697  grad_u_h.add_scaled ((*dphi)[i][sp],
698  system.current_solution(dof_indices[i]));
699 
700  // Compute the phi values at the current sample point
701  std::vector<Real> psi(specpoly(dim, element_order, q_point[sp], matsize));
702  for (unsigned int i=0; i<matsize; i++)
703  {
704  temperr[0] += psi[i]*Pu_x_h(i);
705  }
706 
707  temperr[0] -= grad_u_h(0);
708  }
709 #if LIBMESH_DIM > 1
710  else if (error_estimator.error_norm.type(var) == H1_Y_SEMINORM)
711  {
712  // Compute the gradient at the current sample point
713  Gradient grad_u_h;
714 
715  for (unsigned int i=0; i<n_dofs; i++)
716  grad_u_h.add_scaled ((*dphi)[i][sp],
717  system.current_solution(dof_indices[i]));
718 
719  // Compute the phi values at the current sample point
720  std::vector<Real> psi(specpoly(dim, element_order, q_point[sp], matsize));
721  for (unsigned int i=0; i<matsize; i++)
722  {
723  temperr[1] += psi[i]*Pu_y_h(i);
724  }
725 
726  temperr[1] -= grad_u_h(1);
727  }
728 #endif // LIBMESH_DIM > 1
729 #if LIBMESH_DIM > 2
730  else if (error_estimator.error_norm.type(var) == H1_Z_SEMINORM)
731  {
732  // Compute the gradient at the current sample point
733  Gradient grad_u_h;
734 
735  for (unsigned int i=0; i<n_dofs; i++)
736  grad_u_h.add_scaled ((*dphi)[i][sp],
737  system.current_solution(dof_indices[i]));
738 
739  // Compute the phi values at the current sample point
740  std::vector<Real> psi(specpoly(dim, element_order, q_point[sp], matsize));
741  for (unsigned int i=0; i<matsize; i++)
742  {
743  temperr[2] += psi[i]*Pu_z_h(i);
744  }
745 
746  temperr[2] -= grad_u_h(2);
747  }
748 #endif // LIBMESH_DIM > 2
749  else if (error_estimator.error_norm.type(var) == H2_SEMINORM ||
750  error_estimator.error_norm.type(var) == W2_INF_SEMINORM)
751  {
752 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
753  // Compute the Hessian at the current sample point
754  Tensor hess_u_h;
755 
756  for (unsigned int i=0; i<n_dofs; i++)
757  hess_u_h.add_scaled ((*d2phi)[i][sp],
758  system.current_solution(dof_indices[i]));
759 
760  // Compute the phi values at the current sample point
761  std::vector<Real> psi(specpoly(dim, element_order, q_point[sp], matsize));
762  for (unsigned int i=0; i<matsize; i++)
763  {
764  temperr[0] += psi[i]*Pu_x_h(i);
765 #if LIBMESH_DIM > 1
766  temperr[1] += psi[i]*Pu_y_h(i);
767  temperr[3] += psi[i]*Pu_xy_h(i);
768 #endif
769 #if LIBMESH_DIM > 2
770  temperr[2] += psi[i]*Pu_z_h(i);
771  temperr[4] += psi[i]*Pu_xz_h(i);
772  temperr[5] += psi[i]*Pu_yz_h(i);
773 #endif
774  }
775 
776  temperr[0] -= hess_u_h(0,0);
777 #if LIBMESH_DIM > 1
778  temperr[1] -= hess_u_h(1,1);
779  temperr[3] -= hess_u_h(0,1);
780 #endif
781 #if LIBMESH_DIM > 2
782  temperr[2] -= hess_u_h(2,2);
783  temperr[4] -= hess_u_h(0,2);
784  temperr[5] -= hess_u_h(1,2);
785 #endif
786 #else
787  libmesh_error_msg("ERROR: --enable-second-derivatives is required \nfor _sobolev_order == 2!");
788 #endif
789  }
790 
791  // Get the weight from the user code
792  Number weight = (*error_estimator.weight_functions[var])(femcontext, q_point[sp], system.time);
793 
794  // Add up relevant terms. We can easily optimize the
795  // LIBMESH_DIM < 3 cases a little bit with the exception
796  // of the W2 cases
797 
798  if (error_estimator.error_norm.type(var) == L_INF)
799  element_error = std::max(element_error, std::abs(weight*temperr[0]));
800  else if (error_estimator.error_norm.type(var) == W1_INF_SEMINORM)
801  for (unsigned int i=0; i != LIBMESH_DIM; ++i)
802  element_error = std::max(element_error, std::abs(weight*temperr[i]));
803  else if (error_estimator.error_norm.type(var) == W2_INF_SEMINORM)
804  for (unsigned int i=0; i != 6; ++i)
805  element_error = std::max(element_error, std::abs(weight*temperr[i]));
806  else if (error_estimator.error_norm.type(var) == L2)
807  element_error += JxW[sp]*TensorTools::norm_sq(weight*temperr[0]);
808  else if (error_estimator.error_norm.type(var) == H1_SEMINORM)
809  for (unsigned int i=0; i != LIBMESH_DIM; ++i)
810  element_error += JxW[sp]*TensorTools::norm_sq(weight*temperr[i]);
811  else if (error_estimator.error_norm.type(var) == H1_X_SEMINORM)
812  element_error += JxW[sp]*TensorTools::norm_sq(weight*temperr[0]);
813  else if (error_estimator.error_norm.type(var) == H1_Y_SEMINORM)
814  element_error += JxW[sp]*TensorTools::norm_sq(weight*temperr[1]);
815  else if (error_estimator.error_norm.type(var) == H1_Z_SEMINORM)
816  element_error += JxW[sp]*TensorTools::norm_sq(weight*temperr[2]);
817  else if (error_estimator.error_norm.type(var) == H2_SEMINORM)
818  {
819  for (unsigned int i=0; i != LIBMESH_DIM; ++i)
820  element_error += JxW[sp]*TensorTools::norm_sq(weight*temperr[i]);
821  // Off diagonal terms enter into the Hessian norm twice
822  for (unsigned int i=3; i != 6; ++i)
823  element_error += JxW[sp]*2*TensorTools::norm_sq(weight*temperr[i]);
824  }
825 
826  } // End loop over sample points
827 
828  if (error_estimator.error_norm.type(var) == L_INF ||
829  error_estimator.error_norm.type(var) == W1_INF_SEMINORM ||
830  error_estimator.error_norm.type(var) == W2_INF_SEMINORM)
831  new_error_per_cell[e] += error_estimator.error_norm.weight(var) * element_error;
832  else if (error_estimator.error_norm.type(var) == L2 ||
833  error_estimator.error_norm.type(var) == H1_SEMINORM ||
834  error_estimator.error_norm.type(var) == H1_X_SEMINORM ||
835  error_estimator.error_norm.type(var) == H1_Y_SEMINORM ||
836  error_estimator.error_norm.type(var) == H1_Z_SEMINORM ||
837  error_estimator.error_norm.type(var) == H2_SEMINORM)
838  new_error_per_cell[e] += error_estimator.error_norm.weight_sq(var) * element_error;
839  else
840  libmesh_error_msg("Unsupported error norm type!");
841  } // End (re) loop over patch elements
842 
843  } // end variables loop
844 
845  // Now that we have the contributions from each variable,
846  // we have take square roots of the entries we
847  // added to error_per_cell to get an error norm
848  // If we are reusing patches, once again reuse the current patch to loop
849  // over all elements in the current patch, otherwise build a new
850  // patch containing just the current element and loop over it
851  Patch::const_iterator patch_re_it;
852  Patch::const_iterator patch_re_end;
853 
854  // Build a new patch if necessary
855  Patch current_elem_patch(mesh.processor_id());
856 
857  if (this->error_estimator.patch_reuse)
858  {
859  // Just get the iterators from the current patch
860  patch_re_it = patch.begin();
861  patch_re_end = patch.end();
862  }
863  else
864  {
865  // Use a target patch size of just 0, this will contain
866  // just the current element.
867  current_elem_patch.build_around_element (elem, 0,
868  error_estimator.patch_growth_strategy);
869 
870  // Get the iterators from this newly constructed patch
871  patch_re_it = current_elem_patch.begin();
872  patch_re_end = current_elem_patch.end();
873  }
874 
875  // Loop over every element in the patch we just constructed
876  for (unsigned int i = 0 ; patch_re_it != patch_re_end; ++patch_re_it, ++i)
877  {
878  // The pth element in the patch
879  const Elem * e_p = *patch_re_it;
880 
881  // We'll need an index into the error vector
882  const dof_id_type e_p_id = e_p->id();
883 
884  // Update the error_per_cell vector for this element
885  if (error_estimator.error_norm.type(0) == L2 ||
886  error_estimator.error_norm.type(0) == H1_SEMINORM ||
887  error_estimator.error_norm.type(0) == H1_X_SEMINORM ||
888  error_estimator.error_norm.type(0) == H1_Y_SEMINORM ||
889  error_estimator.error_norm.type(0) == H1_Z_SEMINORM ||
890  error_estimator.error_norm.type(0) == H2_SEMINORM)
891  {
892  Threads::spin_mutex::scoped_lock acquire(Threads::spin_mtx);
893  if (!error_per_cell[e_p_id])
894  error_per_cell[e_p_id] = static_cast<ErrorVectorReal>
895  (std::sqrt(new_error_per_cell[i]));
896  }
897  else
898  {
899  libmesh_assert (error_estimator.error_norm.type(0) == L_INF ||
900  error_estimator.error_norm.type(0) == W1_INF_SEMINORM ||
901  error_estimator.error_norm.type(0) == W2_INF_SEMINORM);
902  Threads::spin_mutex::scoped_lock acquire(Threads::spin_mtx);
903  if (!error_per_cell[e_p_id])
904  error_per_cell[e_p_id] = static_cast<ErrorVectorReal>
905  (new_error_per_cell[i]);
906  }
907 
908  } // End loop over every element in patch
909 
910 
911  } // end element loop
912 
913 } // End () operator definition

References std::abs(), libMesh::TypeVector< T >::add_scaled(), libMesh::TypeTensor< T >::add_scaled(), libMesh::FEGenericBase< OutputType >::build(), libMesh::Patch::build_around_element(), libMesh::System::current_solution(), libMesh::FEType::default_quadrature_rule(), dim, libMesh::DofMap::dof_indices(), error_estimator, libMesh::ErrorEstimator::error_norm, error_per_cell, libMesh::System::get_dof_map(), libMesh::System::get_mesh(), libMesh::H1_SEMINORM, libMesh::H1_X_SEMINORM, libMesh::H1_Y_SEMINORM, libMesh::H1_Z_SEMINORM, libMesh::H2_SEMINORM, libMesh::DofObject::id(), libMesh::L2, libMesh::L_INF, libMesh::libmesh_assert(), libMesh::DenseMatrix< T >::lu_solve(), libMesh::DenseMatrixBase< T >::m(), mesh, libMesh::DenseMatrixBase< T >::n(), n_vars, libMesh::System::n_vars(), libMesh::TensorTools::norm_sq(), libMesh::FEType::order, libMesh::Elem::p_level(), libMesh::PatchRecoveryErrorEstimator::patch_growth_strategy, libMesh::PatchRecoveryErrorEstimator::patch_reuse, libMesh::FEMContext::pre_fe_reinit(), libMesh::Real, libMesh::DenseVector< T >::resize(), libMesh::PatchRecoveryErrorEstimator::specpoly(), libMesh::Threads::spin_mtx, std::sqrt(), system, libMesh::PatchRecoveryErrorEstimator::target_patch_size, libMesh::System::time, libMesh::SystemNorm::type(), libMesh::DofMap::variable_type(), libMesh::W1_INF_SEMINORM, libMesh::W2_INF_SEMINORM, libMesh::MeshTools::weight(), libMesh::SystemNorm::weight(), libMesh::WeightedPatchRecoveryErrorEstimator::weight_functions, libMesh::SystemNorm::weight_sq(), and libMesh::zero.

Member Data Documentation

◆ error_estimator

const WeightedPatchRecoveryErrorEstimator& libMesh::WeightedPatchRecoveryErrorEstimator::EstimateError::error_estimator
private

Definition at line 114 of file weighted_patch_recovery_error_estimator.h.

Referenced by operator()().

◆ error_per_cell

ErrorVector& libMesh::WeightedPatchRecoveryErrorEstimator::EstimateError::error_per_cell
private

Definition at line 115 of file weighted_patch_recovery_error_estimator.h.

Referenced by operator()().

◆ system

const System& libMesh::WeightedPatchRecoveryErrorEstimator::EstimateError::system
private

Function to set the boolean patch_reuse in case the user wants to change the default behaviour of patch_recovery_error_estimator.

Definition at line 113 of file weighted_patch_recovery_error_estimator.h.

Referenced by operator()().

The documentation for this class was generated from the following files:
libMesh::PatchRecoveryErrorEstimator::patch_growth_strategy
Patch::PMF patch_growth_strategy
The PatchErrorEstimator will use this pointer to a Patch member function when growing patches.
Definition: patch_recovery_error_estimator.h:99
libMesh::Number
Real Number
Definition: libmesh_common.h:195
libMesh::dof_id_type
uint8_t dof_id_type
Definition: id_types.h:67
libMesh::System::n_vars
unsigned int n_vars() const
Definition: system.h:2155
libMesh::ErrorEstimator::error_norm
SystemNorm error_norm
When estimating the error in a single system, the error_norm is used to control the scaling and norm ...
Definition: error_estimator.h:164
n_vars
unsigned int n_vars
Definition: adaptivity_ex3.C:116
libMesh::L_INF
Definition: enum_norm_type.h:42
libMesh::H1_Y_SEMINORM
Definition: enum_norm_type.h:58
libMesh::H1_SEMINORM
Definition: enum_norm_type.h:43
libMesh::Order
Order
Definition: enum_order.h:40
libMesh::SystemNorm::weight_sq
Real weight_sq(unsigned int var) const
Definition: system_norm.C:176
std::sqrt
MetaPhysicL::DualNumber< T, D > sqrt(const MetaPhysicL::DualNumber< T, D > &in)
libMesh::W2_INF_SEMINORM
Definition: enum_norm_type.h:50
mesh
MeshBase & mesh
Definition: mesh_communication.C:1257
libMesh::SystemNorm::weight
Real weight(unsigned int var) const
Definition: system_norm.C:133
libMesh::W1_INF_SEMINORM
Definition: enum_norm_type.h:49
libMesh::WeightedPatchRecoveryErrorEstimator::weight_functions
std::vector< FEMFunctionBase< Number > * > weight_functions
Vector of fem function base pointers, the user will fill this in with pointers to the appropriate wei...
Definition: weighted_patch_recovery_error_estimator.h:83
dim
unsigned int dim
Definition: adaptivity_ex3.C:113
libMesh::WeightedPatchRecoveryErrorEstimator::EstimateError::error_estimator
const WeightedPatchRecoveryErrorEstimator & error_estimator
Definition: weighted_patch_recovery_error_estimator.h:114
libMesh::SystemNorm::type
FEMNormType type(unsigned int var) const
Definition: system_norm.C:111
libMesh::zero
const Number zero
.
Definition: libmesh.h:243
libMesh::System::current_solution
Number current_solution(const dof_id_type global_dof_number) const
Definition: system.C:194
libMesh::libmesh_assert
libmesh_assert(ctx)
libMesh::WeightedPatchRecoveryErrorEstimator::EstimateError::error_per_cell
ErrorVector & error_per_cell
Definition: weighted_patch_recovery_error_estimator.h:115
std::abs
MetaPhysicL::DualNumber< T, D > abs(const MetaPhysicL::DualNumber< T, D > &in)
libMesh::Threads::spin_mtx
spin_mutex spin_mtx
A convenient spin mutex object which can be used for obtaining locks.
Definition: threads.C:29
libMesh::System::get_mesh
const MeshBase & get_mesh() const
Definition: system.h:2083
libMesh::FEGenericBase::build
static std::unique_ptr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.
libMesh::TypeVector::add_scaled
void add_scaled(const TypeVector< T2 > &, const T &)
Add a scaled value to this vector without creating a temporary.
Definition: type_vector.h:665
libMesh::TypeTensor::add_scaled
void add_scaled(const TypeTensor< T2 > &, const T &)
Add a scaled tensor to this tensor without creating a temporary.
Definition: type_tensor.h:869
libMesh::H1_X_SEMINORM
Definition: enum_norm_type.h:57
libMesh::PatchRecoveryErrorEstimator::target_patch_size
unsigned int target_patch_size
The PatchErrorEstimator will build patches of at least this many elements to perform estimates.
Definition: patch_recovery_error_estimator.h:91
libMesh::H2_SEMINORM
Definition: enum_norm_type.h:44
libMesh::System::time
Real time
For time-dependent problems, this is the time t at the beginning of the current timestep.
Definition: system.h:1561
libMesh::Gradient
NumberVectorValue Gradient
Definition: exact_solution.h:58
libMesh::TensorTools::norm_sq
T norm_sq(std::complex< T > a)
Definition: tensor_tools.h:85
libMesh::PatchRecoveryErrorEstimator::patch_reuse
bool patch_reuse
Definition: patch_recovery_error_estimator.h:115
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::H1_Z_SEMINORM
Definition: enum_norm_type.h:59
libMesh::L2
Definition: enum_norm_type.h:36
libMesh::System::get_dof_map
const DofMap & get_dof_map() const
Definition: system.h:2099
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:121
libMesh::Tensor
NumberTensorValue Tensor
Definition: exact_solution.h:56
libMesh::MeshTools::weight
dof_id_type weight(const MeshBase &mesh, const processor_id_type pid)
Definition: mesh_tools.C:236
libMesh::WeightedPatchRecoveryErrorEstimator::EstimateError::system
const System & system
Function to set the boolean patch_reuse in case the user wants to change the default behaviour of pat...
Definition: weighted_patch_recovery_error_estimator.h:113
Generated on Sat Jan 25 2020 12:08:57 for libMesh by   doxygen 1.8.16