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

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

Public Member Functions

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

Private Attributes

const Systemsystem
 Function to set the boolean patch_reuse in case the user wants to change the default behaviour of patch_recovery_error_estimator. More...
 
const WeightedPatchRecoveryErrorEstimatorerror_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 93 of file weighted_patch_recovery_error_estimator.h.

Constructor & Destructor Documentation

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

Definition at line 96 of file weighted_patch_recovery_error_estimator.h.

References operator()().

98  :
99  system(sys),
100  error_estimator(ee),
101  error_per_cell(epc)
102  {}
ImplicitSystem & sys
const System & system
Function to set the boolean patch_reuse in case the user wants to change the default behaviour of pat...

Member Function Documentation

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

Definition at line 101 of file weighted_patch_recovery_estimator.C.

References std::abs(), libMesh::TypeVector< T >::add_scaled(), libMesh::TypeTensor< T >::add_scaled(), libMesh::StoredRange< iterator_type, object_type >::begin(), libMesh::FEGenericBase< OutputType >::build(), libMesh::Patch::build_around_element(), libMesh::FEType::default_quadrature_rule(), dim, libMesh::DofMap::dof_indices(), libMesh::StoredRange< iterator_type, object_type >::end(), libMesh::ErrorVectorReal, 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_nullptr, libMesh::DenseMatrix< T >::lu_solve(), libMesh::DenseMatrixBase< T >::m(), std::max(), mesh, libMesh::MeshBase::mesh_dimension(), libMesh::DenseMatrixBase< T >::n(), n_vars, libMesh::TensorTools::norm_sq(), libMesh::FEType::order, libMesh::Elem::p_level(), libMesh::FEMContext::pre_fe_reinit(), libMesh::ParallelObject::processor_id(), libMesh::Real, libMesh::DenseVector< T >::resize(), libMesh::PatchRecoveryErrorEstimator::specpoly(), libMesh::Threads::spin_mtx, libMesh::DiffContext::time, libMesh::DofMap::variable_type(), libMesh::W1_INF_SEMINORM, libMesh::W2_INF_SEMINORM, libMesh::MeshTools::weight(), and libMesh::zero.

Referenced by EstimateError().

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

Member Data Documentation

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

Definition at line 114 of file weighted_patch_recovery_error_estimator.h.

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

Definition at line 115 of file weighted_patch_recovery_error_estimator.h.

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

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

Definition at line 113 of file weighted_patch_recovery_error_estimator.h.


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