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

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

Public Member Functions

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

Private Attributes

const System & system
 
const PatchRecoveryErrorEstimator & 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 123 of file patch_recovery_error_estimator.h.

Constructor & Destructor Documentation

◆ EstimateError()

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

Definition at line 126 of file patch_recovery_error_estimator.h.

128  :
129  system(sys),
130  error_estimator(ee),
131  error_per_cell(epc)
132  {}

Member Function Documentation

◆ operator()()

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

Definition at line 210 of file patch_recovery_error_estimator.C.

211 {
212  // The current mesh
213  const MeshBase & mesh = system.get_mesh();
214 
215  // The dimensionality of the mesh
216  const unsigned int dim = mesh.mesh_dimension();
217 
218  // The number of variables in the system
219  const unsigned int n_vars = system.n_vars();
220 
221  // The DofMap for this system
222  const DofMap & dof_map = system.get_dof_map();
223 
224  //------------------------------------------------------------
225  // Iterate over all the elements in the range.
226  for (const auto & elem : range)
227  {
228  // We'll need an index into the error vector
229  const dof_id_type e_id=elem->id();
230 
231  // We are going to build a patch containing the current element
232  // and its neighbors on the local processor
233  Patch patch(mesh.processor_id());
234 
235  // If we are reusing patches and the current element
236  // already has an estimate associated with it, move on the
237  // next element
238  if (this->error_estimator.patch_reuse && error_per_cell[e_id] != 0)
239  continue;
240 
241  // If we are not reusing patches or havent built one containing this element, we build one
242 
243  // Use user specified patch size and growth strategy
244  patch.build_around_element (elem, error_estimator.target_patch_size,
245  error_estimator.patch_growth_strategy);
246 
247  // Declare a new_error_per_cell vector to hold error estimates
248  // from each element in this patch, or one estimate if we are
249  // not reusing patches since we will only be computing error for
250  // one cell
251  std::vector<Real> new_error_per_cell(1, 0.);
252  if (this->error_estimator.patch_reuse)
253  new_error_per_cell.resize(patch.size(), 0.);
254 
255  //------------------------------------------------------------
256  // Process each variable in the system using the current patch
257  for (unsigned int var=0; var<n_vars; var++)
258  {
259 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
260 #ifdef DEBUG
261  bool is_valid_norm_type =
262  error_estimator.error_norm.type(var) == L2 ||
263  error_estimator.error_norm.type(var) == H1_SEMINORM ||
264  error_estimator.error_norm.type(var) == H2_SEMINORM ||
265  error_estimator.error_norm.type(var) == H1_X_SEMINORM ||
266  error_estimator.error_norm.type(var) == H1_Y_SEMINORM ||
267  error_estimator.error_norm.type(var) == H1_Z_SEMINORM ||
268  error_estimator.error_norm.type(var) == L_INF ||
269  error_estimator.error_norm.type(var) == W1_INF_SEMINORM ||
270  error_estimator.error_norm.type(var) == W2_INF_SEMINORM;
271  libmesh_assert (is_valid_norm_type);
272 #endif // DEBUG
273 #else
274  libmesh_assert (error_estimator.error_norm.type(var) == L2 ||
275  error_estimator.error_norm.type(var) == L_INF ||
276  error_estimator.error_norm.type(var) == H1_SEMINORM ||
277  error_estimator.error_norm.type(var) == H1_X_SEMINORM ||
278  error_estimator.error_norm.type(var) == H1_Y_SEMINORM ||
279  error_estimator.error_norm.type(var) == H1_Z_SEMINORM ||
280  error_estimator.error_norm.type(var) == W1_INF_SEMINORM);
281 #endif
282 
283 
284 #ifdef DEBUG
285  if (var > 0)
286  {
287  // We can't mix L_inf and L_2 norms
288  bool is_valid_norm_combo =
289  ((error_estimator.error_norm.type(var) == L2 ||
290  error_estimator.error_norm.type(var) == H1_SEMINORM ||
291  error_estimator.error_norm.type(var) == H1_X_SEMINORM ||
292  error_estimator.error_norm.type(var) == H1_Y_SEMINORM ||
293  error_estimator.error_norm.type(var) == H1_Z_SEMINORM ||
294  error_estimator.error_norm.type(var) == H2_SEMINORM) &&
295  (error_estimator.error_norm.type(var-1) == L2 ||
296  error_estimator.error_norm.type(var-1) == H1_SEMINORM ||
297  error_estimator.error_norm.type(var-1) == H1_X_SEMINORM ||
298  error_estimator.error_norm.type(var-1) == H1_Y_SEMINORM ||
299  error_estimator.error_norm.type(var-1) == H1_Z_SEMINORM ||
300  error_estimator.error_norm.type(var-1) == H2_SEMINORM)) ||
301  ((error_estimator.error_norm.type(var) == L_INF ||
302  error_estimator.error_norm.type(var) == W1_INF_SEMINORM ||
303  error_estimator.error_norm.type(var) == W2_INF_SEMINORM) &&
304  (error_estimator.error_norm.type(var-1) == L_INF ||
305  error_estimator.error_norm.type(var-1) == W1_INF_SEMINORM ||
306  error_estimator.error_norm.type(var-1) == W2_INF_SEMINORM));
307  libmesh_assert (is_valid_norm_combo);
308  }
309 #endif // DEBUG
310 
311  // Possibly skip this variable
312  if (error_estimator.error_norm.weight(var) == 0.0) continue;
313 
314  // The type of finite element to use for this variable
315  const FEType & fe_type = dof_map.variable_type (var);
316 
317  const Order element_order = static_cast<Order>
318  (fe_type.order + elem->p_level());
319 
320  // Finite element object for use in this patch
321  std::unique_ptr<FEBase> fe (FEBase::build (dim, fe_type));
322 
323  // Build an appropriate Gaussian quadrature rule
324  std::unique_ptr<QBase> qrule (fe_type.default_quadrature_rule(dim));
325 
326  // Tell the finite element about the quadrature rule.
327  fe->attach_quadrature_rule (qrule.get());
328 
329  // Get Jacobian values, etc..
330  const std::vector<Real> & JxW = fe->get_JxW();
331  const std::vector<Point> & q_point = fe->get_xyz();
332 
333  // Get whatever phi/dphi/d2phi values we need. Avoid
334  // getting them unless the requested norm is actually going
335  // to use them.
336 
337  const std::vector<std::vector<Real>> * phi = nullptr;
338  // If we're using phi to assert the correct dof_indices
339  // vector size later, then we'll need to get_phi whether we
340  // plan to use it or not.
341 #ifdef NDEBUG
342  if (error_estimator.error_norm.type(var) == L2 ||
343  error_estimator.error_norm.type(var) == L_INF)
344 #endif
345  phi = &(fe->get_phi());
346 
347  const std::vector<std::vector<RealGradient>> * dphi = nullptr;
348  if (error_estimator.error_norm.type(var) == H1_SEMINORM ||
349  error_estimator.error_norm.type(var) == H1_X_SEMINORM ||
350  error_estimator.error_norm.type(var) == H1_Y_SEMINORM ||
351  error_estimator.error_norm.type(var) == H1_Z_SEMINORM ||
352  error_estimator.error_norm.type(var) == W1_INF_SEMINORM)
353  dphi = &(fe->get_dphi());
354 
355 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
356  const std::vector<std::vector<RealTensor>> * d2phi = nullptr;
357  if (error_estimator.error_norm.type(var) == H2_SEMINORM ||
358  error_estimator.error_norm.type(var) == W2_INF_SEMINORM)
359  d2phi = &(fe->get_d2phi());
360 #endif
361 
362  // global DOF indices
363  std::vector<dof_id_type> dof_indices;
364 
365  // Compute the appropriate size for the patch projection matrices
366  // and vectors;
367  unsigned int matsize = element_order + 1;
368  if (dim > 1)
369  {
370  matsize *= (element_order + 2);
371  matsize /= 2;
372  }
373  if (dim > 2)
374  {
375  matsize *= (element_order + 3);
376  matsize /= 3;
377  }
378 
379  DenseMatrix<Number> Kp(matsize,matsize);
380  DenseVector<Number> F, Fx, Fy, Fz, Fxy, Fxz, Fyz;
381  DenseVector<Number> Pu_h, Pu_x_h, Pu_y_h, Pu_z_h, Pu_xy_h, Pu_xz_h, Pu_yz_h;
382  if (error_estimator.error_norm.type(var) == L2 ||
383  error_estimator.error_norm.type(var) == L_INF)
384  {
385  F.resize(matsize); Pu_h.resize(matsize);
386  }
387  else if (error_estimator.error_norm.type(var) == H1_SEMINORM ||
388  error_estimator.error_norm.type(var) == W1_INF_SEMINORM ||
389  error_estimator.error_norm.type(var) == H2_SEMINORM ||
390  error_estimator.error_norm.type(var) == W2_INF_SEMINORM)
391  {
392  Fx.resize(matsize); Pu_x_h.resize(matsize); // stores xx in W2 cases
393 #if LIBMESH_DIM > 1
394  Fy.resize(matsize); Pu_y_h.resize(matsize); // stores yy in W2 cases
395 #endif
396 #if LIBMESH_DIM > 2
397  Fz.resize(matsize); Pu_z_h.resize(matsize); // stores zz in W2 cases
398 #endif
399  }
400  else if (error_estimator.error_norm.type(var) == H1_X_SEMINORM)
401  {
402  Fx.resize(matsize); Pu_x_h.resize(matsize); // Only need to compute the x gradient for the x component seminorm
403  }
404  else if (error_estimator.error_norm.type(var) == H1_Y_SEMINORM)
405  {
406  libmesh_assert_greater (LIBMESH_DIM, 1);
407  Fy.resize(matsize); Pu_y_h.resize(matsize); // Only need to compute the y gradient for the y component seminorm
408  }
409  else if (error_estimator.error_norm.type(var) == H1_Z_SEMINORM)
410  {
411  libmesh_assert_greater (LIBMESH_DIM, 2);
412  Fz.resize(matsize); Pu_z_h.resize(matsize); // Only need to compute the z gradient for the z component seminorm
413  }
414 
415 #if LIBMESH_DIM > 1
416  if (error_estimator.error_norm.type(var) == H2_SEMINORM ||
417  error_estimator.error_norm.type(var) == W2_INF_SEMINORM)
418  {
419  Fxy.resize(matsize); Pu_xy_h.resize(matsize);
420 #if LIBMESH_DIM > 2
421  Fxz.resize(matsize); Pu_xz_h.resize(matsize);
422  Fyz.resize(matsize); Pu_yz_h.resize(matsize);
423 #endif
424  }
425 #endif
426 
427  //------------------------------------------------------
428  // Loop over each element in the patch and compute their
429  // contribution to the patch gradient projection.
430  for (const auto & e_p : patch)
431  {
432  // Reinitialize the finite element data for this element
433  fe->reinit (e_p);
434 
435  // Get the global DOF indices for the current variable
436  // in the current element
437  dof_map.dof_indices (e_p, dof_indices, var);
438  libmesh_assert_equal_to (dof_indices.size(), phi->size());
439 
440  const unsigned int n_dofs =
441  cast_int<unsigned int>(dof_indices.size());
442  const unsigned int n_qp = qrule->n_points();
443 
444  // Compute the projection components from this cell.
445  // \int_{Omega_e} \psi_i \psi_j = \int_{Omega_e} du_h/dx_k \psi_i
446  for (unsigned int qp=0; qp<n_qp; qp++)
447  {
448  // Construct the shape function values for the patch projection
449  std::vector<Real> psi(specpoly(dim, element_order, q_point[qp], matsize));
450 
451  const unsigned int psi_size = cast_int<unsigned int>(psi.size());
452 
453  // Patch matrix contribution
454  const unsigned int m = Kp.m(), n = Kp.n();
455  for (unsigned int i=0; i<m; i++)
456  for (unsigned int j=0; j<n; j++)
457  Kp(i,j) += JxW[qp]*psi[i]*psi[j];
458 
459  if (error_estimator.error_norm.type(var) == L2 ||
460  error_estimator.error_norm.type(var) == L_INF)
461  {
462  // Compute the solution on the current patch element
463  // the quadrature point
464  Number u_h = libMesh::zero;
465 
466  for (unsigned int i=0; i<n_dofs; i++)
467  u_h += (*phi)[i][qp]*system.current_solution (dof_indices[i]);
468 
469  // Patch RHS contributions
470  for (unsigned int i=0; i != psi_size; i++)
471  F(i) += JxW[qp]*u_h*psi[i];
472 
473  }
474  else if (error_estimator.error_norm.type(var) == H1_SEMINORM ||
475  error_estimator.error_norm.type(var) == W1_INF_SEMINORM)
476  {
477  // Compute the gradient on the current patch element
478  // at the quadrature point
479  Gradient grad_u_h;
480 
481  for (unsigned int i=0; i<n_dofs; i++)
482  grad_u_h.add_scaled ((*dphi)[i][qp],
483  system.current_solution(dof_indices[i]));
484 
485  // Patch RHS contributions
486  for (unsigned int i=0; i != psi_size; i++)
487  {
488  Fx(i) += JxW[qp]*grad_u_h(0)*psi[i];
489 #if LIBMESH_DIM > 1
490  Fy(i) += JxW[qp]*grad_u_h(1)*psi[i];
491 #endif
492 #if LIBMESH_DIM > 2
493  Fz(i) += JxW[qp]*grad_u_h(2)*psi[i];
494 #endif
495  }
496  }
497  else if (error_estimator.error_norm.type(var) == H1_X_SEMINORM)
498  {
499  // Compute the gradient on the current patch element
500  // at the quadrature point
501  Gradient grad_u_h;
502 
503  for (unsigned int i=0; i<n_dofs; i++)
504  grad_u_h.add_scaled ((*dphi)[i][qp],
505  system.current_solution(dof_indices[i]));
506 
507  // Patch RHS contributions
508  for (unsigned int i=0; i != psi_size; i++)
509  {
510  Fx(i) += JxW[qp]*grad_u_h(0)*psi[i];
511  }
512  }
513 #if LIBMESH_DIM > 1
514  else if (error_estimator.error_norm.type(var) == H1_Y_SEMINORM)
515  {
516  // Compute the gradient on the current patch element
517  // at the quadrature point
518  Gradient grad_u_h;
519 
520  for (unsigned int i=0; i<n_dofs; i++)
521  grad_u_h.add_scaled ((*dphi)[i][qp],
522  system.current_solution(dof_indices[i]));
523 
524  // Patch RHS contributions
525  for (unsigned int i=0; i != psi_size; i++)
526  {
527  Fy(i) += JxW[qp]*grad_u_h(1)*psi[i];
528  }
529  }
530 #endif // LIBMESH_DIM > 1
531 #if LIBMESH_DIM > 2
532  else if (error_estimator.error_norm.type(var) == H1_Z_SEMINORM)
533  {
534  // Compute the gradient on the current patch element
535  // at the quadrature point
536  Gradient grad_u_h;
537 
538  for (unsigned int i=0; i<n_dofs; i++)
539  grad_u_h.add_scaled ((*dphi)[i][qp],
540  system.current_solution(dof_indices[i]));
541 
542  // Patch RHS contributions
543  for (unsigned int i=0; i != psi_size; i++)
544  {
545  Fz(i) += JxW[qp]*grad_u_h(2)*psi[i];
546  }
547  }
548 #endif // LIBMESH_DIM > 2
549  else if (error_estimator.error_norm.type(var) == H2_SEMINORM ||
550  error_estimator.error_norm.type(var) == W2_INF_SEMINORM)
551  {
552 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
553  // Compute the hessian on the current patch element
554  // at the quadrature point
555  Tensor hess_u_h;
556 
557  for (unsigned int i=0; i<n_dofs; i++)
558  hess_u_h.add_scaled ((*d2phi)[i][qp],
559  system.current_solution(dof_indices[i]));
560 
561  // Patch RHS contributions
562  for (unsigned int i=0; i != psi_size; i++)
563  {
564  Fx(i) += JxW[qp]*hess_u_h(0,0)*psi[i];
565 #if LIBMESH_DIM > 1
566  Fy(i) += JxW[qp]*hess_u_h(1,1)*psi[i];
567  Fxy(i) += JxW[qp]*hess_u_h(0,1)*psi[i];
568 #endif
569 #if LIBMESH_DIM > 2
570  Fz(i) += JxW[qp]*hess_u_h(2,2)*psi[i];
571  Fxz(i) += JxW[qp]*hess_u_h(0,2)*psi[i];
572  Fyz(i) += JxW[qp]*hess_u_h(1,2)*psi[i];
573 #endif
574  }
575 #else
576  libmesh_error_msg("ERROR: --enable-second-derivatives is required \nfor _sobolev_order == 2!");
577 #endif
578  }
579  else
580  libmesh_error_msg("Unsupported error norm type!");
581  } // end quadrature loop
582  } // end patch loop
583 
584 
585 
586  //--------------------------------------------------
587  // Now we have fully assembled the projection system
588  // for this patch. Project the gradient components.
589  // MAY NEED TO USE PARTIAL PIVOTING!
590  if (error_estimator.error_norm.type(var) == L2 ||
591  error_estimator.error_norm.type(var) == L_INF)
592  {
593  Kp.lu_solve(F, Pu_h);
594  }
595  else if (error_estimator.error_norm.type(var) == H1_SEMINORM ||
596  error_estimator.error_norm.type(var) == W1_INF_SEMINORM ||
597  error_estimator.error_norm.type(var) == H2_SEMINORM ||
598  error_estimator.error_norm.type(var) == W2_INF_SEMINORM)
599  {
600  Kp.lu_solve (Fx, Pu_x_h);
601 #if LIBMESH_DIM > 1
602  Kp.lu_solve (Fy, Pu_y_h);
603 #endif
604 #if LIBMESH_DIM > 2
605  Kp.lu_solve (Fz, Pu_z_h);
606 #endif
607  }
608  else if (error_estimator.error_norm.type(var) == H1_X_SEMINORM)
609  {
610  Kp.lu_solve (Fx, Pu_x_h);
611  }
612  else if (error_estimator.error_norm.type(var) == H1_Y_SEMINORM)
613  {
614  Kp.lu_solve (Fy, Pu_y_h);
615  }
616  else if (error_estimator.error_norm.type(var) == H1_Z_SEMINORM)
617  {
618  Kp.lu_solve (Fz, Pu_z_h);
619  }
620 
621 #if LIBMESH_DIM > 1
622  if (error_estimator.error_norm.type(var) == H2_SEMINORM ||
623  error_estimator.error_norm.type(var) == W2_INF_SEMINORM)
624  {
625  Kp.lu_solve(Fxy, Pu_xy_h);
626 #if LIBMESH_DIM > 2
627  Kp.lu_solve(Fxz, Pu_xz_h);
628  Kp.lu_solve(Fyz, Pu_yz_h);
629 #endif
630  }
631 #endif
632 
633  // If we are reusing patches, reuse the current patch to loop
634  // over all elements in the current patch, otherwise build a new
635  // patch containing just the current element and loop over it
636  // Note that C++ will not allow patch_re_end to be a const here
637  Patch::const_iterator patch_re_it;
638  Patch::const_iterator patch_re_end;
639 
640  // Declare a new patch
641  Patch patch_re(mesh.processor_id());
642 
643  if (this->error_estimator.patch_reuse)
644  {
645  // Just get the iterators from the current patch
646  patch_re_it = patch.begin();
647  patch_re_end = patch.end();
648  }
649  else
650  {
651  // Use a target patch size of just 0, this will contain
652  // just the current element
653  patch_re.build_around_element (elem, 0,
654  error_estimator.patch_growth_strategy);
655 
656  // Get the iterators from this newly constructed patch
657  patch_re_it = patch_re.begin();
658  patch_re_end = patch_re.end();
659  }
660 
661  // If we are reusing patches, loop over all the elements
662  // in the current patch and develop an estimate
663  // for all the elements by computing ||P u_h - u_h|| or ||P grad_u_h - grad_u_h||
664  // or ||P hess_u_h - hess_u_h|| according to the requested
665  // seminorm, otherwise just compute it for the current element
666 
667  // Loop over every element in the patch
668  for (unsigned int e = 0 ; patch_re_it != patch_re_end; ++patch_re_it, ++e)
669  {
670  // Build the Finite Element for the current element
671 
672  // The pth element in the patch
673  const Elem * e_p = *patch_re_it;
674 
675  // We'll need an index into the error vector for this element
676  const dof_id_type e_p_id = e_p->id();
677 
678  // We will update the new_error_per_cell vector with element_error if the
679  // error_per_cell[e_p_id] entry is non-zero, otherwise update it
680  // with 0. i.e. leave it unchanged
681 
682  // No need to compute the estimate if we are reusing patches and already have one
683  if (this->error_estimator.patch_reuse && error_per_cell[e_p_id] != 0.)
684  continue;
685 
686  // Reinitialize the finite element data for this element
687  fe->reinit (e_p);
688 
689  // Get the global DOF indices for the current variable
690  // in the current element
691  dof_map.dof_indices (e_p, dof_indices, var);
692  libmesh_assert_equal_to (dof_indices.size(), phi->size());
693 
694  // The number of dofs for this variable on this element
695  const unsigned int n_dofs =
696  cast_int<unsigned int>(dof_indices.size());
697 
698  // Variable to hold the error on the current element
699  Real element_error = 0;
700 
701  const Order qorder =
702  static_cast<Order>(fe_type.order + e_p->p_level());
703 
704  // A quadrature rule for this element
705  QGrid samprule (dim, qorder);
706 
707  if (error_estimator.error_norm.type(var) == W1_INF_SEMINORM ||
708  error_estimator.error_norm.type(var) == W2_INF_SEMINORM)
709  fe->attach_quadrature_rule (&samprule);
710 
711  // The number of points we will sample over
712  const unsigned int n_sp =
713  cast_int<unsigned int>(JxW.size());
714 
715  // Loop over every sample point for the current element
716  for (unsigned int sp=0; sp<n_sp; sp++)
717  {
718  // Compute the solution at the current sample point
719 
720  std::vector<Number> temperr(6,0.0); // x,y,z or xx,yy,zz,xy,xz,yz
721 
722  if (error_estimator.error_norm.type(var) == L2 ||
723  error_estimator.error_norm.type(var) == L_INF)
724  {
725  // Compute the value at the current sample point
726  Number u_h = libMesh::zero;
727 
728  for (unsigned int i=0; i<n_dofs; i++)
729  u_h += (*phi)[i][sp]*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[0] += psi[i]*Pu_h(i);
736  }
737 
738  temperr[0] -= u_h;
739  }
740  else if (error_estimator.error_norm.type(var) == H1_SEMINORM ||
741  error_estimator.error_norm.type(var) == W1_INF_SEMINORM)
742  {
743  // Compute the gradient at the current sample point
744  Gradient grad_u_h;
745 
746  for (unsigned int i=0; i<n_dofs; i++)
747  grad_u_h.add_scaled ((*dphi)[i][sp],
748  system.current_solution(dof_indices[i]));
749 
750  // Compute the phi values at the current sample point
751  std::vector<Real> psi(specpoly(dim, element_order, q_point[sp], matsize));
752 
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 #endif
759 #if LIBMESH_DIM > 2
760  temperr[2] += psi[i]*Pu_z_h(i);
761 #endif
762  }
763  temperr[0] -= grad_u_h(0);
764 #if LIBMESH_DIM > 1
765  temperr[1] -= grad_u_h(1);
766 #endif
767 #if LIBMESH_DIM > 2
768  temperr[2] -= grad_u_h(2);
769 #endif
770  }
771  else if (error_estimator.error_norm.type(var) == H1_X_SEMINORM)
772  {
773  // Compute the gradient at the current sample point
774  Gradient grad_u_h;
775 
776  for (unsigned int i=0; i<n_dofs; i++)
777  grad_u_h.add_scaled ((*dphi)[i][sp],
778  system.current_solution(dof_indices[i]));
779 
780  // Compute the phi values at the current sample point
781  std::vector<Real> psi(specpoly(dim, element_order, q_point[sp], matsize));
782  for (unsigned int i=0; i<matsize; i++)
783  {
784  temperr[0] += psi[i]*Pu_x_h(i);
785  }
786 
787  temperr[0] -= grad_u_h(0);
788  }
789 #if LIBMESH_DIM > 1
790  else if (error_estimator.error_norm.type(var) == H1_Y_SEMINORM)
791  {
792  // Compute the gradient at the current sample point
793  Gradient grad_u_h;
794 
795  for (unsigned int i=0; i<n_dofs; i++)
796  grad_u_h.add_scaled ((*dphi)[i][sp],
797  system.current_solution(dof_indices[i]));
798 
799  // Compute the phi values at the current sample point
800  std::vector<Real> psi(specpoly(dim, element_order, q_point[sp], matsize));
801  for (unsigned int i=0; i<matsize; i++)
802  {
803  temperr[1] += psi[i]*Pu_y_h(i);
804  }
805 
806  temperr[1] -= grad_u_h(1);
807  }
808 #endif // LIBMESH_DIM > 1
809 #if LIBMESH_DIM > 2
810  else if (error_estimator.error_norm.type(var) == H1_Z_SEMINORM)
811  {
812  // Compute the gradient at the current sample point
813  Gradient grad_u_h;
814 
815  for (unsigned int i=0; i<n_dofs; i++)
816  grad_u_h.add_scaled ((*dphi)[i][sp],
817  system.current_solution(dof_indices[i]));
818 
819  // Compute the phi values at the current sample point
820  std::vector<Real> psi(specpoly(dim, element_order, q_point[sp], matsize));
821  for (unsigned int i=0; i<matsize; i++)
822  {
823  temperr[2] += psi[i]*Pu_z_h(i);
824  }
825 
826  temperr[2] -= grad_u_h(2);
827  }
828 #endif // LIBMESH_DIM > 2
829  else if (error_estimator.error_norm.type(var) == H2_SEMINORM ||
830  error_estimator.error_norm.type(var) == W2_INF_SEMINORM)
831  {
832 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
833  // Compute the Hessian at the current sample point
834  Tensor hess_u_h;
835 
836  for (unsigned int i=0; i<n_dofs; i++)
837  hess_u_h.add_scaled ((*d2phi)[i][sp],
838  system.current_solution(dof_indices[i]));
839 
840  // Compute the phi values at the current sample point
841  std::vector<Real> psi(specpoly(dim, element_order, q_point[sp], matsize));
842  for (unsigned int i=0; i<matsize; i++)
843  {
844  temperr[0] += psi[i]*Pu_x_h(i);
845 #if LIBMESH_DIM > 1
846  temperr[1] += psi[i]*Pu_y_h(i);
847  temperr[3] += psi[i]*Pu_xy_h(i);
848 #endif
849 #if LIBMESH_DIM > 2
850  temperr[2] += psi[i]*Pu_z_h(i);
851  temperr[4] += psi[i]*Pu_xz_h(i);
852  temperr[5] += psi[i]*Pu_yz_h(i);
853 #endif
854  }
855 
856  temperr[0] -= hess_u_h(0,0);
857 #if LIBMESH_DIM > 1
858  temperr[1] -= hess_u_h(1,1);
859  temperr[3] -= hess_u_h(0,1);
860 #endif
861 #if LIBMESH_DIM > 2
862  temperr[2] -= hess_u_h(2,2);
863  temperr[4] -= hess_u_h(0,2);
864  temperr[5] -= hess_u_h(1,2);
865 #endif
866 #else
867  libmesh_error_msg("ERROR: --enable-second-derivatives is required \nfor _sobolev_order == 2!");
868 #endif
869  }
870  // Add up relevant terms. We can easily optimize the
871  // LIBMESH_DIM < 3 cases a little bit with the exception
872  // of the W2 cases
873 
874  if (error_estimator.error_norm.type(var) == L_INF)
875  element_error = std::max(element_error, std::abs(temperr[0]));
876  else if (error_estimator.error_norm.type(var) == W1_INF_SEMINORM)
877  for (unsigned int i=0; i != LIBMESH_DIM; ++i)
878  element_error = std::max(element_error, std::abs(temperr[i]));
879  else if (error_estimator.error_norm.type(var) == W2_INF_SEMINORM)
880  for (unsigned int i=0; i != 6; ++i)
881  element_error = std::max(element_error, std::abs(temperr[i]));
882  else if (error_estimator.error_norm.type(var) == L2)
883  element_error += JxW[sp]*TensorTools::norm_sq(temperr[0]);
884  else if (error_estimator.error_norm.type(var) == H1_SEMINORM)
885  for (unsigned int i=0; i != LIBMESH_DIM; ++i)
886  element_error += JxW[sp]*TensorTools::norm_sq(temperr[i]);
887  else if (error_estimator.error_norm.type(var) == H1_X_SEMINORM)
888  element_error += JxW[sp]*TensorTools::norm_sq(temperr[0]);
889  else if (error_estimator.error_norm.type(var) == H1_Y_SEMINORM)
890  element_error += JxW[sp]*TensorTools::norm_sq(temperr[1]);
891  else if (error_estimator.error_norm.type(var) == H1_Z_SEMINORM)
892  element_error += JxW[sp]*TensorTools::norm_sq(temperr[2]);
893  else if (error_estimator.error_norm.type(var) == H2_SEMINORM)
894  {
895  for (unsigned int i=0; i != LIBMESH_DIM; ++i)
896  element_error += JxW[sp]*TensorTools::norm_sq(temperr[i]);
897  // Off diagonal terms enter into the Hessian norm twice
898  for (unsigned int i=3; i != 6; ++i)
899  element_error += JxW[sp]*2*TensorTools::norm_sq(temperr[i]);
900  }
901 
902  } // End loop over sample points
903 
904  if (error_estimator.error_norm.type(var) == L_INF ||
905  error_estimator.error_norm.type(var) == W1_INF_SEMINORM ||
906  error_estimator.error_norm.type(var) == W2_INF_SEMINORM)
907  new_error_per_cell[e] += error_estimator.error_norm.weight(var) * element_error;
908  else if (error_estimator.error_norm.type(var) == L2 ||
909  error_estimator.error_norm.type(var) == H1_SEMINORM ||
910  error_estimator.error_norm.type(var) == H1_X_SEMINORM ||
911  error_estimator.error_norm.type(var) == H1_Y_SEMINORM ||
912  error_estimator.error_norm.type(var) == H1_Z_SEMINORM ||
913  error_estimator.error_norm.type(var) == H2_SEMINORM)
914  new_error_per_cell[e] += error_estimator.error_norm.weight_sq(var) * element_error;
915  else
916  libmesh_error_msg("Unsupported error norm type!");
917  } // End (re) loop over patch elements
918 
919  } // end variables loop
920 
921  // Now that we have the contributions from each variable,
922  // we have take square roots of the entries we
923  // added to error_per_cell to get an error norm
924  // If we are reusing patches, once again reuse the current patch to loop
925  // over all elements in the current patch, otherwise build a new
926  // patch containing just the current element and loop over it
927  Patch::const_iterator patch_re_it;
928  Patch::const_iterator patch_re_end;
929 
930  // Build a new patch if necessary
931  Patch current_elem_patch(mesh.processor_id());
932 
933  if (this->error_estimator.patch_reuse)
934  {
935  // Just get the iterators from the current patch
936  patch_re_it = patch.begin();
937  patch_re_end = patch.end();
938  }
939  else
940  {
941  // Use a target patch size of just 0, this will contain
942  // just the current element.
943  current_elem_patch.build_around_element (elem, 0,
944  error_estimator.patch_growth_strategy);
945 
946  // Get the iterators from this newly constructed patch
947  patch_re_it = current_elem_patch.begin();
948  patch_re_end = current_elem_patch.end();
949  }
950 
951  // Loop over every element in the patch we just constructed
952  for (unsigned int i = 0 ; patch_re_it != patch_re_end; ++patch_re_it, ++i)
953  {
954  // The pth element in the patch
955  const Elem * e_p = *patch_re_it;
956 
957  // We'll need an index into the error vector
958  const dof_id_type e_p_id = e_p->id();
959 
960  // Update the error_per_cell vector for this element
961  if (error_estimator.error_norm.type(0) == L2 ||
962  error_estimator.error_norm.type(0) == H1_SEMINORM ||
963  error_estimator.error_norm.type(0) == H1_X_SEMINORM ||
964  error_estimator.error_norm.type(0) == H1_Y_SEMINORM ||
965  error_estimator.error_norm.type(0) == H1_Z_SEMINORM ||
966  error_estimator.error_norm.type(0) == H2_SEMINORM)
967  {
968  Threads::spin_mutex::scoped_lock acquire(Threads::spin_mtx);
969  if (!error_per_cell[e_p_id])
970  error_per_cell[e_p_id] =
971  static_cast<ErrorVectorReal>(std::sqrt(new_error_per_cell[i]));
972  }
973  else
974  {
975  libmesh_assert (error_estimator.error_norm.type(0) == L_INF ||
976  error_estimator.error_norm.type(0) == W1_INF_SEMINORM ||
977  error_estimator.error_norm.type(0) == W2_INF_SEMINORM);
978  Threads::spin_mutex::scoped_lock acquire(Threads::spin_mtx);
979  if (!error_per_cell[e_p_id])
980  error_per_cell[e_p_id] =
981  static_cast<ErrorVectorReal>(new_error_per_cell[i]);
982  }
983 
984  } // End loop over every element in patch
985 
986  } // end element loop
987 
988 } // 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::Real, libMesh::DenseVector< T >::resize(), libMesh::PatchRecoveryErrorEstimator::specpoly(), libMesh::Threads::spin_mtx, std::sqrt(), system, libMesh::PatchRecoveryErrorEstimator::target_patch_size, libMesh::SystemNorm::type(), libMesh::DofMap::variable_type(), libMesh::W1_INF_SEMINORM, libMesh::W2_INF_SEMINORM, libMesh::SystemNorm::weight(), libMesh::SystemNorm::weight_sq(), and libMesh::zero.

Member Data Documentation

◆ error_estimator

const PatchRecoveryErrorEstimator& libMesh::PatchRecoveryErrorEstimator::EstimateError::error_estimator
private

Definition at line 138 of file patch_recovery_error_estimator.h.

Referenced by operator()().

◆ error_per_cell

ErrorVector& libMesh::PatchRecoveryErrorEstimator::EstimateError::error_per_cell
private

Definition at line 139 of file patch_recovery_error_estimator.h.

Referenced by operator()().

◆ system

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

Definition at line 137 of file 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::PatchRecoveryErrorEstimator::EstimateError::error_estimator
const PatchRecoveryErrorEstimator & error_estimator
Definition: patch_recovery_error_estimator.h:138
dim
unsigned int dim
Definition: adaptivity_ex3.C:113
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::PatchRecoveryErrorEstimator::EstimateError::error_per_cell
ErrorVector & error_per_cell
Definition: patch_recovery_error_estimator.h:139
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::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::PatchRecoveryErrorEstimator::EstimateError::system
const System & system
Definition: patch_recovery_error_estimator.h:137
Generated on Sat Jan 25 2020 12:08:34 for libMesh by   doxygen 1.8.16