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

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.

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(), libMesh::DofMap::dof_indices(), 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::DenseMatrix< T >::lu_solve(), libMesh::DenseMatrixBase< T >::m(), std::max(), 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.

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  for (unsigned int i=0; i<Kp.m(); i++)
455  for (unsigned int j=0; j<Kp.n(); j++)
456  Kp(i,j) += JxW[qp]*psi[i]*psi[j];
457 
458  if (error_estimator.error_norm.type(var) == L2 ||
459  error_estimator.error_norm.type(var) == L_INF)
460  {
461  // Compute the solution on the current patch element
462  // the quadrature point
463  Number u_h = libMesh::zero;
464 
465  for (unsigned int i=0; i<n_dofs; i++)
466  u_h += (*phi)[i][qp]*system.current_solution (dof_indices[i]);
467 
468  // Patch RHS contributions
469  for (unsigned int i=0; i != psi_size; i++)
470  F(i) += JxW[qp]*u_h*psi[i];
471 
472  }
473  else if (error_estimator.error_norm.type(var) == H1_SEMINORM ||
474  error_estimator.error_norm.type(var) == W1_INF_SEMINORM)
475  {
476  // Compute the gradient on the current patch element
477  // at the quadrature point
478  Gradient grad_u_h;
479 
480  for (unsigned int i=0; i<n_dofs; i++)
481  grad_u_h.add_scaled ((*dphi)[i][qp],
482  system.current_solution(dof_indices[i]));
483 
484  // Patch RHS contributions
485  for (unsigned int i=0; i != psi_size; i++)
486  {
487  Fx(i) += JxW[qp]*grad_u_h(0)*psi[i];
488 #if LIBMESH_DIM > 1
489  Fy(i) += JxW[qp]*grad_u_h(1)*psi[i];
490 #endif
491 #if LIBMESH_DIM > 2
492  Fz(i) += JxW[qp]*grad_u_h(2)*psi[i];
493 #endif
494  }
495  }
496  else if (error_estimator.error_norm.type(var) == H1_X_SEMINORM)
497  {
498  // Compute the gradient on the current patch element
499  // at the quadrature point
500  Gradient grad_u_h;
501 
502  for (unsigned int i=0; i<n_dofs; i++)
503  grad_u_h.add_scaled ((*dphi)[i][qp],
504  system.current_solution(dof_indices[i]));
505 
506  // Patch RHS contributions
507  for (unsigned int i=0; i != psi_size; i++)
508  {
509  Fx(i) += JxW[qp]*grad_u_h(0)*psi[i];
510  }
511  }
512  else if (error_estimator.error_norm.type(var) == H1_Y_SEMINORM)
513  {
514  // Compute the gradient on the current patch element
515  // at the quadrature point
516  Gradient grad_u_h;
517 
518  for (unsigned int i=0; i<n_dofs; i++)
519  grad_u_h.add_scaled ((*dphi)[i][qp],
520  system.current_solution(dof_indices[i]));
521 
522  // Patch RHS contributions
523  for (unsigned int i=0; i != psi_size; i++)
524  {
525  Fy(i) += JxW[qp]*grad_u_h(1)*psi[i];
526  }
527  }
528  else if (error_estimator.error_norm.type(var) == H1_Z_SEMINORM)
529  {
530  // Compute the gradient on the current patch element
531  // at the quadrature point
532  Gradient grad_u_h;
533 
534  for (unsigned int i=0; i<n_dofs; i++)
535  grad_u_h.add_scaled ((*dphi)[i][qp],
536  system.current_solution(dof_indices[i]));
537 
538  // Patch RHS contributions
539  for (unsigned int i=0; i != psi_size; i++)
540  {
541  Fz(i) += JxW[qp]*grad_u_h(2)*psi[i];
542  }
543  }
544  else if (error_estimator.error_norm.type(var) == H2_SEMINORM ||
545  error_estimator.error_norm.type(var) == W2_INF_SEMINORM)
546  {
547 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
548  // Compute the hessian on the current patch element
549  // at the quadrature point
550  Tensor hess_u_h;
551 
552  for (unsigned int i=0; i<n_dofs; i++)
553  hess_u_h.add_scaled ((*d2phi)[i][qp],
554  system.current_solution(dof_indices[i]));
555 
556  // Patch RHS contributions
557  for (unsigned int i=0; i != psi_size; i++)
558  {
559  Fx(i) += JxW[qp]*hess_u_h(0,0)*psi[i];
560 #if LIBMESH_DIM > 1
561  Fy(i) += JxW[qp]*hess_u_h(1,1)*psi[i];
562  Fxy(i) += JxW[qp]*hess_u_h(0,1)*psi[i];
563 #endif
564 #if LIBMESH_DIM > 2
565  Fz(i) += JxW[qp]*hess_u_h(2,2)*psi[i];
566  Fxz(i) += JxW[qp]*hess_u_h(0,2)*psi[i];
567  Fyz(i) += JxW[qp]*hess_u_h(1,2)*psi[i];
568 #endif
569  }
570 #else
571  libmesh_error_msg("ERROR: --enable-second-derivatives is required \nfor _sobolev_order == 2!");
572 #endif
573  }
574  else
575  libmesh_error_msg("Unsupported error norm type!");
576  } // end quadrature loop
577  } // end patch loop
578 
579 
580 
581  //--------------------------------------------------
582  // Now we have fully assembled the projection system
583  // for this patch. Project the gradient components.
584  // MAY NEED TO USE PARTIAL PIVOTING!
585  if (error_estimator.error_norm.type(var) == L2 ||
586  error_estimator.error_norm.type(var) == L_INF)
587  {
588  Kp.lu_solve(F, Pu_h);
589  }
590  else if (error_estimator.error_norm.type(var) == H1_SEMINORM ||
591  error_estimator.error_norm.type(var) == W1_INF_SEMINORM ||
592  error_estimator.error_norm.type(var) == H2_SEMINORM ||
593  error_estimator.error_norm.type(var) == W2_INF_SEMINORM)
594  {
595  Kp.lu_solve (Fx, Pu_x_h);
596 #if LIBMESH_DIM > 1
597  Kp.lu_solve (Fy, Pu_y_h);
598 #endif
599 #if LIBMESH_DIM > 2
600  Kp.lu_solve (Fz, Pu_z_h);
601 #endif
602  }
603  else if (error_estimator.error_norm.type(var) == H1_X_SEMINORM)
604  {
605  Kp.lu_solve (Fx, Pu_x_h);
606  }
607  else if (error_estimator.error_norm.type(var) == H1_Y_SEMINORM)
608  {
609  Kp.lu_solve (Fy, Pu_y_h);
610  }
611  else if (error_estimator.error_norm.type(var) == H1_Z_SEMINORM)
612  {
613  Kp.lu_solve (Fz, Pu_z_h);
614  }
615 
616 #if LIBMESH_DIM > 1
617  if (error_estimator.error_norm.type(var) == H2_SEMINORM ||
618  error_estimator.error_norm.type(var) == W2_INF_SEMINORM)
619  {
620  Kp.lu_solve(Fxy, Pu_xy_h);
621 #if LIBMESH_DIM > 2
622  Kp.lu_solve(Fxz, Pu_xz_h);
623  Kp.lu_solve(Fyz, Pu_yz_h);
624 #endif
625  }
626 #endif
627 
628  // If we are reusing patches, reuse the current patch to loop
629  // over all elements in the current patch, otherwise build a new
630  // patch containing just the current element and loop over it
631  // Note that C++ will not allow patch_re_end to be a const here
632  Patch::const_iterator patch_re_it;
633  Patch::const_iterator patch_re_end;
634 
635  // Declare a new patch
636  Patch patch_re(mesh.processor_id());
637 
638  if (this->error_estimator.patch_reuse)
639  {
640  // Just get the iterators from the current patch
641  patch_re_it = patch.begin();
642  patch_re_end = patch.end();
643  }
644  else
645  {
646  // Use a target patch size of just 0, this will contain
647  // just the current element
648  patch_re.build_around_element (elem, 0,
649  error_estimator.patch_growth_strategy);
650 
651  // Get the iterators from this newly constructed patch
652  patch_re_it = patch_re.begin();
653  patch_re_end = patch_re.end();
654  }
655 
656  // If we are reusing patches, loop over all the elements
657  // in the current patch and develop an estimate
658  // for all the elements by computing ||P u_h - u_h|| or ||P grad_u_h - grad_u_h||
659  // or ||P hess_u_h - hess_u_h|| according to the requested
660  // seminorm, otherwise just compute it for the current element
661 
662  // Loop over every element in the patch
663  for (unsigned int e = 0 ; patch_re_it != patch_re_end; ++patch_re_it, ++e)
664  {
665  // Build the Finite Element for the current element
666 
667  // The pth element in the patch
668  const Elem * e_p = *patch_re_it;
669 
670  // We'll need an index into the error vector for this element
671  const dof_id_type e_p_id = e_p->id();
672 
673  // We will update the new_error_per_cell vector with element_error if the
674  // error_per_cell[e_p_id] entry is non-zero, otherwise update it
675  // with 0. i.e. leave it unchanged
676 
677  // No need to compute the estimate if we are reusing patches and already have one
678  if (this->error_estimator.patch_reuse && error_per_cell[e_p_id] != 0.)
679  continue;
680 
681  // Reinitialize the finite element data for this element
682  fe->reinit (e_p);
683 
684  // Get the global DOF indices for the current variable
685  // in the current element
686  dof_map.dof_indices (e_p, dof_indices, var);
687  libmesh_assert_equal_to (dof_indices.size(), phi->size());
688 
689  // The number of dofs for this variable on this element
690  const unsigned int n_dofs =
691  cast_int<unsigned int>(dof_indices.size());
692 
693  // Variable to hold the error on the current element
694  Real element_error = 0;
695 
696  const Order qorder =
697  static_cast<Order>(fe_type.order + e_p->p_level());
698 
699  // A quadrature rule for this element
700  QGrid samprule (dim, qorder);
701 
702  if (error_estimator.error_norm.type(var) == W1_INF_SEMINORM ||
703  error_estimator.error_norm.type(var) == W2_INF_SEMINORM)
704  fe->attach_quadrature_rule (&samprule);
705 
706  // The number of points we will sample over
707  const unsigned int n_sp =
708  cast_int<unsigned int>(JxW.size());
709 
710  // Loop over every sample point for the current element
711  for (unsigned int sp=0; sp<n_sp; sp++)
712  {
713  // Compute the solution at the current sample point
714 
715  std::vector<Number> temperr(6,0.0); // x,y,z or xx,yy,zz,xy,xz,yz
716 
717  if (error_estimator.error_norm.type(var) == L2 ||
718  error_estimator.error_norm.type(var) == L_INF)
719  {
720  // Compute the value at the current sample point
721  Number u_h = libMesh::zero;
722 
723  for (unsigned int i=0; i<n_dofs; i++)
724  u_h += (*phi)[i][sp]*system.current_solution (dof_indices[i]);
725 
726  // Compute the phi values at the current sample point
727  std::vector<Real> psi(specpoly(dim, element_order, q_point[sp], matsize));
728  for (unsigned int i=0; i<matsize; i++)
729  {
730  temperr[0] += psi[i]*Pu_h(i);
731  }
732 
733  temperr[0] -= u_h;
734  }
735  else if (error_estimator.error_norm.type(var) == H1_SEMINORM ||
736  error_estimator.error_norm.type(var) == W1_INF_SEMINORM)
737  {
738  // Compute the gradient at the current sample point
739  Gradient grad_u_h;
740 
741  for (unsigned int i=0; i<n_dofs; i++)
742  grad_u_h.add_scaled ((*dphi)[i][sp],
743  system.current_solution(dof_indices[i]));
744 
745  // Compute the phi values at the current sample point
746  std::vector<Real> psi(specpoly(dim, element_order, q_point[sp], matsize));
747 
748  for (unsigned int i=0; i<matsize; i++)
749  {
750  temperr[0] += psi[i]*Pu_x_h(i);
751 #if LIBMESH_DIM > 1
752  temperr[1] += psi[i]*Pu_y_h(i);
753 #endif
754 #if LIBMESH_DIM > 2
755  temperr[2] += psi[i]*Pu_z_h(i);
756 #endif
757  }
758  temperr[0] -= grad_u_h(0);
759 #if LIBMESH_DIM > 1
760  temperr[1] -= grad_u_h(1);
761 #endif
762 #if LIBMESH_DIM > 2
763  temperr[2] -= grad_u_h(2);
764 #endif
765  }
766  else if (error_estimator.error_norm.type(var) == H1_X_SEMINORM)
767  {
768  // Compute the gradient at the current sample point
769  Gradient grad_u_h;
770 
771  for (unsigned int i=0; i<n_dofs; i++)
772  grad_u_h.add_scaled ((*dphi)[i][sp],
773  system.current_solution(dof_indices[i]));
774 
775  // Compute the phi values at the current sample point
776  std::vector<Real> psi(specpoly(dim, element_order, q_point[sp], matsize));
777  for (unsigned int i=0; i<matsize; i++)
778  {
779  temperr[0] += psi[i]*Pu_x_h(i);
780  }
781 
782  temperr[0] -= grad_u_h(0);
783  }
784  else if (error_estimator.error_norm.type(var) == H1_Y_SEMINORM)
785  {
786  // Compute the gradient at the current sample point
787  Gradient grad_u_h;
788 
789  for (unsigned int i=0; i<n_dofs; i++)
790  grad_u_h.add_scaled ((*dphi)[i][sp],
791  system.current_solution(dof_indices[i]));
792 
793  // Compute the phi values at the current sample point
794  std::vector<Real> psi(specpoly(dim, element_order, q_point[sp], matsize));
795  for (unsigned int i=0; i<matsize; i++)
796  {
797  temperr[1] += psi[i]*Pu_y_h(i);
798  }
799 
800  temperr[1] -= grad_u_h(1);
801  }
802  else if (error_estimator.error_norm.type(var) == H1_Z_SEMINORM)
803  {
804  // Compute the gradient at the current sample point
805  Gradient grad_u_h;
806 
807  for (unsigned int i=0; i<n_dofs; i++)
808  grad_u_h.add_scaled ((*dphi)[i][sp],
809  system.current_solution(dof_indices[i]));
810 
811  // Compute the phi values at the current sample point
812  std::vector<Real> psi(specpoly(dim, element_order, q_point[sp], matsize));
813  for (unsigned int i=0; i<matsize; i++)
814  {
815  temperr[2] += psi[i]*Pu_z_h(i);
816  }
817 
818  temperr[2] -= grad_u_h(2);
819  }
820  else if (error_estimator.error_norm.type(var) == H2_SEMINORM ||
821  error_estimator.error_norm.type(var) == W2_INF_SEMINORM)
822  {
823 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
824  // Compute the Hessian at the current sample point
825  Tensor hess_u_h;
826 
827  for (unsigned int i=0; i<n_dofs; i++)
828  hess_u_h.add_scaled ((*d2phi)[i][sp],
829  system.current_solution(dof_indices[i]));
830 
831  // Compute the phi values at the current sample point
832  std::vector<Real> psi(specpoly(dim, element_order, q_point[sp], matsize));
833  for (unsigned int i=0; i<matsize; i++)
834  {
835  temperr[0] += psi[i]*Pu_x_h(i);
836 #if LIBMESH_DIM > 1
837  temperr[1] += psi[i]*Pu_y_h(i);
838  temperr[3] += psi[i]*Pu_xy_h(i);
839 #endif
840 #if LIBMESH_DIM > 2
841  temperr[2] += psi[i]*Pu_z_h(i);
842  temperr[4] += psi[i]*Pu_xz_h(i);
843  temperr[5] += psi[i]*Pu_yz_h(i);
844 #endif
845  }
846 
847  temperr[0] -= hess_u_h(0,0);
848 #if LIBMESH_DIM > 1
849  temperr[1] -= hess_u_h(1,1);
850  temperr[3] -= hess_u_h(0,1);
851 #endif
852 #if LIBMESH_DIM > 2
853  temperr[2] -= hess_u_h(2,2);
854  temperr[4] -= hess_u_h(0,2);
855  temperr[5] -= hess_u_h(1,2);
856 #endif
857 #else
858  libmesh_error_msg("ERROR: --enable-second-derivatives is required \nfor _sobolev_order == 2!");
859 #endif
860  }
861  // Add up relevant terms. We can easily optimize the
862  // LIBMESH_DIM < 3 cases a little bit with the exception
863  // of the W2 cases
864 
865  if (error_estimator.error_norm.type(var) == L_INF)
866  element_error = std::max(element_error, std::abs(temperr[0]));
867  else if (error_estimator.error_norm.type(var) == W1_INF_SEMINORM)
868  for (unsigned int i=0; i != LIBMESH_DIM; ++i)
869  element_error = std::max(element_error, std::abs(temperr[i]));
870  else if (error_estimator.error_norm.type(var) == W2_INF_SEMINORM)
871  for (unsigned int i=0; i != 6; ++i)
872  element_error = std::max(element_error, std::abs(temperr[i]));
873  else if (error_estimator.error_norm.type(var) == L2)
874  element_error += JxW[sp]*TensorTools::norm_sq(temperr[0]);
875  else if (error_estimator.error_norm.type(var) == H1_SEMINORM)
876  for (unsigned int i=0; i != LIBMESH_DIM; ++i)
877  element_error += JxW[sp]*TensorTools::norm_sq(temperr[i]);
878  else if (error_estimator.error_norm.type(var) == H1_X_SEMINORM)
879  element_error += JxW[sp]*TensorTools::norm_sq(temperr[0]);
880  else if (error_estimator.error_norm.type(var) == H1_Y_SEMINORM)
881  element_error += JxW[sp]*TensorTools::norm_sq(temperr[1]);
882  else if (error_estimator.error_norm.type(var) == H1_Z_SEMINORM)
883  element_error += JxW[sp]*TensorTools::norm_sq(temperr[2]);
884  else if (error_estimator.error_norm.type(var) == H2_SEMINORM)
885  {
886  for (unsigned int i=0; i != LIBMESH_DIM; ++i)
887  element_error += JxW[sp]*TensorTools::norm_sq(temperr[i]);
888  // Off diagonal terms enter into the Hessian norm twice
889  for (unsigned int i=3; i != 6; ++i)
890  element_error += JxW[sp]*2*TensorTools::norm_sq(temperr[i]);
891  }
892 
893  } // End loop over sample points
894 
895  if (error_estimator.error_norm.type(var) == L_INF ||
896  error_estimator.error_norm.type(var) == W1_INF_SEMINORM ||
897  error_estimator.error_norm.type(var) == W2_INF_SEMINORM)
898  new_error_per_cell[e] += error_estimator.error_norm.weight(var) * element_error;
899  else if (error_estimator.error_norm.type(var) == L2 ||
900  error_estimator.error_norm.type(var) == H1_SEMINORM ||
901  error_estimator.error_norm.type(var) == H1_X_SEMINORM ||
902  error_estimator.error_norm.type(var) == H1_Y_SEMINORM ||
903  error_estimator.error_norm.type(var) == H1_Z_SEMINORM ||
904  error_estimator.error_norm.type(var) == H2_SEMINORM)
905  new_error_per_cell[e] += error_estimator.error_norm.weight_sq(var) * element_error;
906  else
907  libmesh_error_msg("Unsupported error norm type!");
908  } // End (re) loop over patch elements
909 
910  } // end variables loop
911 
912  // Now that we have the contributions from each variable,
913  // we have take square roots of the entries we
914  // added to error_per_cell to get an error norm
915  // If we are reusing patches, once again reuse the current patch to loop
916  // over all elements in the current patch, otherwise build a new
917  // patch containing just the current element and loop over it
918  Patch::const_iterator patch_re_it;
919  Patch::const_iterator patch_re_end;
920 
921  // Build a new patch if necessary
922  Patch current_elem_patch(mesh.processor_id());
923 
924  if (this->error_estimator.patch_reuse)
925  {
926  // Just get the iterators from the current patch
927  patch_re_it = patch.begin();
928  patch_re_end = patch.end();
929  }
930  else
931  {
932  // Use a target patch size of just 0, this will contain
933  // just the current element.
934  current_elem_patch.build_around_element (elem, 0,
935  error_estimator.patch_growth_strategy);
936 
937  // Get the iterators from this newly constructed patch
938  patch_re_it = current_elem_patch.begin();
939  patch_re_end = current_elem_patch.end();
940  }
941 
942  // Loop over every element in the patch we just constructed
943  for (unsigned int i = 0 ; patch_re_it != patch_re_end; ++patch_re_it, ++i)
944  {
945  // The pth element in the patch
946  const Elem * e_p = *patch_re_it;
947 
948  // We'll need an index into the error vector
949  const dof_id_type e_p_id = e_p->id();
950 
951  // Update the error_per_cell vector for this element
952  if (error_estimator.error_norm.type(0) == L2 ||
953  error_estimator.error_norm.type(0) == H1_SEMINORM ||
954  error_estimator.error_norm.type(0) == H1_X_SEMINORM ||
955  error_estimator.error_norm.type(0) == H1_Y_SEMINORM ||
956  error_estimator.error_norm.type(0) == H1_Z_SEMINORM ||
957  error_estimator.error_norm.type(0) == H2_SEMINORM)
958  {
959  Threads::spin_mutex::scoped_lock acquire(Threads::spin_mtx);
960  if (!error_per_cell[e_p_id])
961  error_per_cell[e_p_id] =
962  static_cast<ErrorVectorReal>(std::sqrt(new_error_per_cell[i]));
963  }
964  else
965  {
966  libmesh_assert (error_estimator.error_norm.type(0) == L_INF ||
967  error_estimator.error_norm.type(0) == W1_INF_SEMINORM ||
968  error_estimator.error_norm.type(0) == W2_INF_SEMINORM);
969  Threads::spin_mutex::scoped_lock acquire(Threads::spin_mtx);
970  if (!error_per_cell[e_p_id])
971  error_per_cell[e_p_id] =
972  static_cast<ErrorVectorReal>(new_error_per_cell[i]);
973  }
974 
975  } // End loop over every element in patch
976 
977  } // end element loop
978 
979 } // End () operator definition
std::abs
double abs(double a)
Definition: libmesh_augment_std_namespace.h:59
libMesh::TypeTensor::add_scaled
void add_scaled(const TypeTensor< T2 > &, const T)
Definition: type_tensor.h:808
libMesh::TypeVector::add_scaled
void add_scaled(const TypeVector< T2 > &, const T)
Definition: type_vector.h:627
mesh
MeshBase & mesh
Definition: mesh_communication.C:1245
n_vars
const unsigned int n_vars
Definition: tecplot_io.C:69
libMesh::zero
const Number zero
Definition: libmesh.h:239
std::max
long double max(long double a, double b)
Definition: libmesh_augment_std_namespace.h:27
libMesh::System::get_mesh
const MeshBase & get_mesh() const
Definition: system.h:2033
libMesh::ErrorVectorReal
DIE A HORRIBLE DEATH HERE typedef float ErrorVectorReal
Definition: libmesh_common.h:202
libMesh::System::current_solution
Number current_solution(const dof_id_type global_dof_number) const
Definition: system.C:194
libMesh::Threads::spin_mtx
spin_mutex spin_mtx
Definition: threads.C:29
libMesh::SystemNorm::type
FEMNormType type(unsigned int var) const
Definition: system_norm.C:110
libMesh::Gradient
NumberVectorValue Gradient
Definition: exact_solution.h:57
libMesh::FEGenericBase::build
static std::unique_ptr< FEGenericBase > build(const unsigned int dim, const FEType &type)
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::SystemNorm::weight_sq
Real weight_sq(unsigned int var) const
Definition: system_norm.C:175
libMesh::SystemNorm::weight
Real weight(unsigned int var) const
Definition: system_norm.C:132
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:131
libMesh::Tensor
NumberTensorValue Tensor
Definition: exact_solution.h:55
libMesh::System::n_vars
unsigned int n_vars() const
Definition: system.h:2105
libMesh::System::get_dof_map
const DofMap & get_dof_map() const
Definition: system.h:2049
libMesh::dof_id_type
uint8_t dof_id_type
Definition: id_types.h:64

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:
generated by doxygen