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 Systemsystem
 
const PatchRecoveryErrorEstimatorerror_estimator
 
ErrorVectorerror_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 130 of file patch_recovery_error_estimator.h.

Constructor & Destructor Documentation

◆ EstimateError()

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

Member Function Documentation

◆ operator()()

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

Definition at line 211 of file patch_recovery_error_estimator.C.

References libMesh::PatchRecoveryErrorEstimator::_extra_order, 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(), libMesh::Utility::enum_to_string(), error_estimator, libMesh::ErrorEstimator::error_norm, error_per_cell, libMesh::ErrorVectorReal, 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, 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.

212 {
213  // The current mesh
214  const MeshBase & mesh = system.get_mesh();
215 
216  // The dimensionality of the mesh
217  const unsigned int dim = mesh.mesh_dimension();
218 
219  // The number of variables in the system
220  const unsigned int n_vars = system.n_vars();
221 
222  // The DofMap for this system
223  const DofMap & dof_map = system.get_dof_map();
224 
225  //------------------------------------------------------------
226  // Iterate over all the elements in the range.
227  for (const auto & elem : range)
228  {
229  // We'll need an index into the error vector
230  const dof_id_type e_id=elem->id();
231 
232  // We are going to build a patch containing the current element
233  // and its neighbors on the local processor
234  Patch patch(mesh.processor_id());
235 
236  // If we are reusing patches and the current element
237  // already has an estimate associated with it, move on the
238  // next element
239  if (this->error_estimator.patch_reuse && error_per_cell[e_id] != 0)
240  continue;
241 
242  // If we are not reusing patches or haven't built one containing this element, we build one
243 
244  // Use user specified patch size and growth strategy
245  patch.build_around_element (elem, error_estimator.target_patch_size,
247 
248  // Declare a new_error_per_cell vector to hold error estimates
249  // from each element in this patch, or one estimate if we are
250  // not reusing patches since we will only be computing error for
251  // one cell
252  std::vector<Real> new_error_per_cell(1, 0.);
253  if (this->error_estimator.patch_reuse)
254  new_error_per_cell.resize(patch.size(), 0.);
255 
256  //------------------------------------------------------------
257  // Process each variable in the system using the current patch
258  for (unsigned int var=0; var<n_vars; var++)
259  {
260  const auto norm_type = error_estimator.error_norm.type(var);
261 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
262 #ifdef DEBUG
263  bool is_valid_norm_type =
264  norm_type == L2 ||
265  norm_type == H1_SEMINORM ||
266  norm_type == H2_SEMINORM ||
267  norm_type == H1_X_SEMINORM ||
268  norm_type == H1_Y_SEMINORM ||
269  norm_type == H1_Z_SEMINORM ||
270  norm_type == L_INF ||
271  norm_type == W1_INF_SEMINORM ||
272  norm_type == W2_INF_SEMINORM;
273  libmesh_assert (is_valid_norm_type);
274 #endif // DEBUG
275 #else
276  libmesh_assert (norm_type == L2 ||
277  norm_type == L_INF ||
278  norm_type == H1_SEMINORM ||
279  norm_type == H1_X_SEMINORM ||
280  norm_type == H1_Y_SEMINORM ||
281  norm_type == H1_Z_SEMINORM ||
282  norm_type == W1_INF_SEMINORM);
283 #endif
284 
285 
286 #ifdef DEBUG
287  if (var > 0)
288  {
289  // We can't mix L_inf and L_2 norms
290  bool is_valid_norm_combo =
291  ((norm_type == L2 ||
292  norm_type == H1_SEMINORM ||
293  norm_type == H1_X_SEMINORM ||
294  norm_type == H1_Y_SEMINORM ||
295  norm_type == H1_Z_SEMINORM ||
296  norm_type == H2_SEMINORM) &&
297  (error_estimator.error_norm.type(var-1) == L2 ||
303  ((norm_type == L_INF ||
304  norm_type == W1_INF_SEMINORM ||
305  norm_type == W2_INF_SEMINORM) &&
306  (error_estimator.error_norm.type(var-1) == L_INF ||
309  libmesh_assert (is_valid_norm_combo);
310  }
311 #endif // DEBUG
312 
313  // Possibly skip this variable
314  if (error_estimator.error_norm.weight(var) == 0.0) continue;
315 
316  // The type of finite element to use for this variable
317  const FEType & fe_type = dof_map.variable_type (var);
318 
319  const Order element_order = fe_type.order + elem->p_level();
320 
321  // Finite element object for use in this patch
322  std::unique_ptr<FEBase> fe (FEBase::build (dim, fe_type));
323 
324  // Build an appropriate Gaussian quadrature rule
325  std::unique_ptr<QBase> qrule (fe_type.default_quadrature_rule(dim, error_estimator._extra_order));
326 
327  // Tell the finite element about the quadrature rule.
328  fe->attach_quadrature_rule (qrule.get());
329 
330  // Get Jacobian values, etc..
331  const std::vector<Real> & JxW = fe->get_JxW();
332  const std::vector<Point> & q_point = fe->get_xyz();
333 
334  // Get whatever phi/dphi/d2phi values we need. Avoid
335  // getting them unless the requested norm is actually going
336  // to use them.
337 
338  const std::vector<std::vector<Real>> * phi = nullptr;
339  // If we're using phi to assert the correct dof_indices
340  // vector size later, then we'll need to get_phi whether we
341  // plan to use it or not.
342 #ifdef NDEBUG
343  if (norm_type == L2 ||
344  norm_type == L_INF)
345 #endif
346  phi = &(fe->get_phi());
347 
348  const std::vector<std::vector<RealGradient>> * dphi = nullptr;
349  if (norm_type == H1_SEMINORM ||
350  norm_type == H1_X_SEMINORM ||
351  norm_type == H1_Y_SEMINORM ||
352  norm_type == H1_Z_SEMINORM ||
353  norm_type == W1_INF_SEMINORM)
354  dphi = &(fe->get_dphi());
355 
356 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
357  const std::vector<std::vector<RealTensor>> * d2phi = nullptr;
358  if (norm_type == H2_SEMINORM ||
359  norm_type == W2_INF_SEMINORM)
360  d2phi = &(fe->get_d2phi());
361 #endif
362 
363  // global DOF indices
364  std::vector<dof_id_type> dof_indices;
365 
366  // Compute the appropriate size for the patch projection matrices
367  // and vectors;
368  unsigned int matsize = element_order + 1;
369  if (dim > 1)
370  {
371  matsize *= (element_order + 2);
372  matsize /= 2;
373  }
374  if (dim > 2)
375  {
376  matsize *= (element_order + 3);
377  matsize /= 3;
378  }
379 
380  DenseMatrix<Number> Kp(matsize,matsize);
381  DenseVector<Number> F, Fx, Fy, Fz, Fxy, Fxz, Fyz;
382  DenseVector<Number> Pu_h, Pu_x_h, Pu_y_h, Pu_z_h, Pu_xy_h, Pu_xz_h, Pu_yz_h;
383  if (norm_type == L2 ||
384  norm_type == L_INF)
385  {
386  F.resize(matsize); Pu_h.resize(matsize);
387  }
388  else if (norm_type == H1_SEMINORM ||
389  norm_type == W1_INF_SEMINORM ||
390  norm_type == H2_SEMINORM ||
391  norm_type == W2_INF_SEMINORM)
392  {
393  Fx.resize(matsize); Pu_x_h.resize(matsize); // stores xx in W2 cases
394 #if LIBMESH_DIM > 1
395  Fy.resize(matsize); Pu_y_h.resize(matsize); // stores yy in W2 cases
396 #endif
397 #if LIBMESH_DIM > 2
398  Fz.resize(matsize); Pu_z_h.resize(matsize); // stores zz in W2 cases
399 #endif
400  }
401  else if (norm_type == H1_X_SEMINORM)
402  {
403  Fx.resize(matsize); Pu_x_h.resize(matsize); // Only need to compute the x gradient for the x component seminorm
404  }
405  else if (norm_type == H1_Y_SEMINORM)
406  {
407  libmesh_assert_greater (LIBMESH_DIM, 1);
408  Fy.resize(matsize); Pu_y_h.resize(matsize); // Only need to compute the y gradient for the y component seminorm
409  }
410  else if (norm_type == H1_Z_SEMINORM)
411  {
412  libmesh_assert_greater (LIBMESH_DIM, 2);
413  Fz.resize(matsize); Pu_z_h.resize(matsize); // Only need to compute the z gradient for the z component seminorm
414  }
415 
416 #if LIBMESH_DIM > 1
417  if (norm_type == H2_SEMINORM ||
418  norm_type == W2_INF_SEMINORM)
419  {
420  Fxy.resize(matsize); Pu_xy_h.resize(matsize);
421 #if LIBMESH_DIM > 2
422  Fxz.resize(matsize); Pu_xz_h.resize(matsize);
423  Fyz.resize(matsize); Pu_yz_h.resize(matsize);
424 #endif
425  }
426 #endif
427 
428  //------------------------------------------------------
429  // Loop over each element in the patch and compute their
430  // contribution to the patch gradient projection.
431  for (const auto & e_p : patch)
432  {
433  // Reinitialize the finite element data for this element
434  fe->reinit (e_p);
435 
436  // Get the global DOF indices for the current variable
437  // in the current element
438  dof_map.dof_indices (e_p, dof_indices, var);
439  libmesh_assert_equal_to (dof_indices.size(), phi->size());
440 
441  const unsigned int n_dofs =
442  cast_int<unsigned int>(dof_indices.size());
443  const unsigned int n_qp = qrule->n_points();
444 
445  // Compute the projection components from this cell.
446  // \int_{Omega_e} \psi_i \psi_j = \int_{Omega_e} du_h/dx_k \psi_i
447  for (unsigned int qp=0; qp<n_qp; qp++)
448  {
449  // Construct the shape function values for the patch projection
450  std::vector<Real> psi(specpoly(dim, element_order, q_point[qp], matsize));
451 
452  const unsigned int psi_size = cast_int<unsigned int>(psi.size());
453 
454  // Patch matrix contribution
455  const unsigned int m = Kp.m(), n = Kp.n();
456  for (unsigned int i=0; i<m; i++)
457  for (unsigned int j=0; j<n; j++)
458  Kp(i,j) += JxW[qp]*psi[i]*psi[j];
459 
460  if (norm_type == L2 ||
461  norm_type == L_INF)
462  {
463  // Compute the solution on the current patch element
464  // the quadrature point
465  Number u_h = libMesh::zero;
466 
467  for (unsigned int i=0; i<n_dofs; i++)
468  u_h += (*phi)[i][qp]*system.current_solution (dof_indices[i]);
469 
470  // Patch RHS contributions
471  for (unsigned int i=0; i != psi_size; i++)
472  F(i) += JxW[qp]*u_h*psi[i];
473 
474  }
475  else if (norm_type == H1_SEMINORM ||
476  norm_type == W1_INF_SEMINORM)
477  {
478  // Compute the gradient on the current patch element
479  // at the quadrature point
480  Gradient grad_u_h;
481 
482  for (unsigned int i=0; i<n_dofs; i++)
483  grad_u_h.add_scaled ((*dphi)[i][qp],
484  system.current_solution(dof_indices[i]));
485 
486  // Patch RHS contributions
487  for (unsigned int i=0; i != psi_size; i++)
488  {
489  Fx(i) += JxW[qp]*grad_u_h(0)*psi[i];
490 #if LIBMESH_DIM > 1
491  Fy(i) += JxW[qp]*grad_u_h(1)*psi[i];
492 #endif
493 #if LIBMESH_DIM > 2
494  Fz(i) += JxW[qp]*grad_u_h(2)*psi[i];
495 #endif
496  }
497  }
498  else if (norm_type == H1_X_SEMINORM)
499  {
500  // Compute the gradient on the current patch element
501  // at the quadrature point
502  Gradient grad_u_h;
503 
504  for (unsigned int i=0; i<n_dofs; i++)
505  grad_u_h.add_scaled ((*dphi)[i][qp],
506  system.current_solution(dof_indices[i]));
507 
508  // Patch RHS contributions
509  for (unsigned int i=0; i != psi_size; i++)
510  {
511  Fx(i) += JxW[qp]*grad_u_h(0)*psi[i];
512  }
513  }
514 #if LIBMESH_DIM > 1
515  else if (norm_type == H1_Y_SEMINORM)
516  {
517  // Compute the gradient on the current patch element
518  // at the quadrature point
519  Gradient grad_u_h;
520 
521  for (unsigned int i=0; i<n_dofs; i++)
522  grad_u_h.add_scaled ((*dphi)[i][qp],
523  system.current_solution(dof_indices[i]));
524 
525  // Patch RHS contributions
526  for (unsigned int i=0; i != psi_size; i++)
527  {
528  Fy(i) += JxW[qp]*grad_u_h(1)*psi[i];
529  }
530  }
531 #endif // LIBMESH_DIM > 1
532 #if LIBMESH_DIM > 2
533  else if (norm_type == H1_Z_SEMINORM)
534  {
535  // Compute the gradient on the current patch element
536  // at the quadrature point
537  Gradient grad_u_h;
538 
539  for (unsigned int i=0; i<n_dofs; i++)
540  grad_u_h.add_scaled ((*dphi)[i][qp],
541  system.current_solution(dof_indices[i]));
542 
543  // Patch RHS contributions
544  for (unsigned int i=0; i != psi_size; i++)
545  {
546  Fz(i) += JxW[qp]*grad_u_h(2)*psi[i];
547  }
548  }
549 #endif // LIBMESH_DIM > 2
550  else if (norm_type == H2_SEMINORM ||
551  norm_type == W2_INF_SEMINORM)
552  {
553 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
554  // Compute the hessian on the current patch element
555  // at the quadrature point
556  Tensor hess_u_h;
557 
558  for (unsigned int i=0; i<n_dofs; i++)
559  hess_u_h.add_scaled ((*d2phi)[i][qp],
560  system.current_solution(dof_indices[i]));
561 
562  // Patch RHS contributions
563  for (unsigned int i=0; i != psi_size; i++)
564  {
565  Fx(i) += JxW[qp]*hess_u_h(0,0)*psi[i];
566 #if LIBMESH_DIM > 1
567  Fy(i) += JxW[qp]*hess_u_h(1,1)*psi[i];
568  Fxy(i) += JxW[qp]*hess_u_h(0,1)*psi[i];
569 #endif
570 #if LIBMESH_DIM > 2
571  Fz(i) += JxW[qp]*hess_u_h(2,2)*psi[i];
572  Fxz(i) += JxW[qp]*hess_u_h(0,2)*psi[i];
573  Fyz(i) += JxW[qp]*hess_u_h(1,2)*psi[i];
574 #endif
575  }
576 #else
577  libmesh_error_msg("ERROR: --enable-second-derivatives is required \nfor _sobolev_order == 2!");
578 #endif
579  }
580  else
581  libmesh_error_msg("Unsupported error norm type == " << Utility::enum_to_string(norm_type));
582  } // end quadrature loop
583  } // end patch loop
584 
585 
586 
587  //--------------------------------------------------
588  // Now we have fully assembled the projection system
589  // for this patch. Project the gradient components.
590  // MAY NEED TO USE PARTIAL PIVOTING!
591  if (norm_type == L2 ||
592  norm_type == L_INF)
593  {
594  Kp.lu_solve(F, Pu_h);
595  }
596  else if (norm_type == H1_SEMINORM ||
597  norm_type == W1_INF_SEMINORM ||
598  norm_type == H2_SEMINORM ||
599  norm_type == W2_INF_SEMINORM)
600  {
601  Kp.lu_solve (Fx, Pu_x_h);
602 #if LIBMESH_DIM > 1
603  Kp.lu_solve (Fy, Pu_y_h);
604 #endif
605 #if LIBMESH_DIM > 2
606  Kp.lu_solve (Fz, Pu_z_h);
607 #endif
608  }
609  else if (norm_type == H1_X_SEMINORM)
610  {
611  Kp.lu_solve (Fx, Pu_x_h);
612  }
613  else if (norm_type == H1_Y_SEMINORM)
614  {
615  Kp.lu_solve (Fy, Pu_y_h);
616  }
617  else if (norm_type == H1_Z_SEMINORM)
618  {
619  Kp.lu_solve (Fz, Pu_z_h);
620  }
621 
622 #if LIBMESH_DIM > 1
623  if (norm_type == H2_SEMINORM ||
624  norm_type == W2_INF_SEMINORM)
625  {
626  Kp.lu_solve(Fxy, Pu_xy_h);
627 #if LIBMESH_DIM > 2
628  Kp.lu_solve(Fxz, Pu_xz_h);
629  Kp.lu_solve(Fyz, Pu_yz_h);
630 #endif
631  }
632 #endif
633 
634  // If we are reusing patches, reuse the current patch to loop
635  // over all elements in the current patch, otherwise build a new
636  // patch containing just the current element and loop over it
637  // Note that C++ will not allow patch_re_end to be a const here
638  Patch::const_iterator patch_re_it;
639  Patch::const_iterator patch_re_end;
640 
641  // Declare a new patch
642  Patch patch_re(mesh.processor_id());
643 
644  if (this->error_estimator.patch_reuse)
645  {
646  // Just get the iterators from the current patch
647  patch_re_it = patch.begin();
648  patch_re_end = patch.end();
649  }
650  else
651  {
652  // Use a target patch size of just 0, this will contain
653  // just the current element
654  patch_re.build_around_element (elem, 0,
656 
657  // Get the iterators from this newly constructed patch
658  patch_re_it = patch_re.begin();
659  patch_re_end = patch_re.end();
660  }
661 
662  // If we are reusing patches, loop over all the elements
663  // in the current patch and develop an estimate
664  // for all the elements by computing ||P u_h - u_h|| or ||P grad_u_h - grad_u_h||
665  // or ||P hess_u_h - hess_u_h|| according to the requested
666  // seminorm, otherwise just compute it for the current element
667 
668  // Loop over every element in the patch
669  for (unsigned int e = 0 ; patch_re_it != patch_re_end; ++patch_re_it, ++e)
670  {
671  // Build the Finite Element for the current element
672 
673  // The pth element in the patch
674  const Elem * e_p = *patch_re_it;
675 
676  // We'll need an index into the error vector for this element
677  const dof_id_type e_p_id = e_p->id();
678 
679  // We will update the new_error_per_cell vector with element_error if the
680  // error_per_cell[e_p_id] entry is non-zero, otherwise update it
681  // with 0. i.e. leave it unchanged
682 
683  // No need to compute the estimate if we are reusing patches and already have one
684  if (this->error_estimator.patch_reuse && error_per_cell[e_p_id] != 0.)
685  continue;
686 
687  // Reinitialize the finite element data for this element
688  fe->reinit (e_p);
689 
690  // Get the global DOF indices for the current variable
691  // in the current element
692  dof_map.dof_indices (e_p, dof_indices, var);
693  libmesh_assert_equal_to (dof_indices.size(), phi->size());
694 
695  // The number of dofs for this variable on this element
696  const unsigned int n_dofs =
697  cast_int<unsigned int>(dof_indices.size());
698 
699  // Variable to hold the error on the current element
700  Real element_error = 0;
701 
702  const Order qorder = fe_type.order + e_p->p_level();
703 
704  // A quadrature rule for this element
705  QGrid samprule (dim, qorder);
706 
707  if (norm_type == W1_INF_SEMINORM ||
708  norm_type == 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 (norm_type == L2 ||
723  norm_type == 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 (norm_type == H1_SEMINORM ||
741  norm_type == 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 (norm_type == 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 (norm_type == 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 (norm_type == 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 (norm_type == H2_SEMINORM ||
830  norm_type == 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 (norm_type == L_INF)
875  element_error = std::max(element_error, std::abs(temperr[0]));
876  else if (norm_type == 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 (norm_type == 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 (norm_type == L2)
883  element_error += JxW[sp]*TensorTools::norm_sq(temperr[0]);
884  else if (norm_type == 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 (norm_type == H1_X_SEMINORM)
888  element_error += JxW[sp]*TensorTools::norm_sq(temperr[0]);
889  else if (norm_type == H1_Y_SEMINORM)
890  element_error += JxW[sp]*TensorTools::norm_sq(temperr[1]);
891  else if (norm_type == H1_Z_SEMINORM)
892  element_error += JxW[sp]*TensorTools::norm_sq(temperr[2]);
893  else if (norm_type == 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 (norm_type == L_INF ||
905  norm_type == W1_INF_SEMINORM ||
906  norm_type == W2_INF_SEMINORM)
907  new_error_per_cell[e] += error_estimator.error_norm.weight(var) * element_error;
908  else if (norm_type == L2 ||
909  norm_type == H1_SEMINORM ||
910  norm_type == H1_X_SEMINORM ||
911  norm_type == H1_Y_SEMINORM ||
912  norm_type == H1_Z_SEMINORM ||
913  norm_type == 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 == " << Utility::enum_to_string(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,
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 ||
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  {
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
void add_scaled(const TypeVector< T2 > &, const T &)
Add a scaled value to this vector without creating a temporary.
Definition: type_vector.h:629
Order
defines an enum for polynomial orders.
Definition: enum_order.h:40
SystemNorm error_norm
When estimating the error in a single system, the error_norm is used to control the scaling and norm ...
Patch::PMF patch_growth_strategy
The PatchErrorEstimator will use this pointer to a Patch member function when growing patches...
unsigned int dim
MeshBase & mesh
const Number zero
.
Definition: libmesh.h:304
const MeshBase & get_mesh() const
Definition: system.h:2358
DIE A HORRIBLE DEATH HERE typedef float ErrorVectorReal
Number current_solution(const dof_id_type global_dof_number) const
Definition: system.C:165
unsigned int n_vars
FEMNormType type(unsigned int var) const
Definition: system_norm.C:112
NumberVectorValue Gradient
static std::unique_ptr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.
libmesh_assert(ctx)
unsigned int target_patch_size
The PatchErrorEstimator will build patches of at least this many elements to perform estimates...
static std::vector< Real > specpoly(const unsigned int dim, const Order order, const Point p, const unsigned int matsize)
void add_scaled(const TypeTensor< T2 > &, const T &)
Add a scaled tensor to this tensor without creating a temporary.
Definition: type_tensor.h:851
std::string enum_to_string(const T e)
Real weight_sq(unsigned int var) const
Definition: system_norm.C:177
Real weight(unsigned int var) const
Definition: system_norm.C:134
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
NumberTensorValue Tensor
auto norm_sq(const T &a) -> decltype(std::norm(a))
Definition: tensor_tools.h:104
int _extra_order
Extra order to use for quadrature rule.
unsigned int n_vars() const
Definition: system.h:2430
const DofMap & get_dof_map() const
Definition: system.h:2374
uint8_t dof_id_type
Definition: id_types.h:67
spin_mutex spin_mtx
A convenient spin mutex object which can be used for obtaining locks.
Definition: threads.C:30

Member Data Documentation

◆ error_estimator

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

Definition at line 145 of file patch_recovery_error_estimator.h.

Referenced by operator()().

◆ error_per_cell

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

Definition at line 146 of file patch_recovery_error_estimator.h.

Referenced by operator()().

◆ system

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

Definition at line 144 of file patch_recovery_error_estimator.h.

Referenced by operator()().


The documentation for this class was generated from the following files: