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

◆ EstimateError()

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

Definition at line 96 of file weighted_patch_recovery_error_estimator.h.

98  :
99  system(sys),
100  error_estimator(ee),
101  error_per_cell(epc)
102  {}
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

◆ operator()()

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

Definition at line 112 of file weighted_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::FEMContext::pre_fe_reinit(), libMesh::Real, libMesh::DenseVector< T >::resize(), libMesh::PatchRecoveryErrorEstimator::specpoly(), libMesh::Threads::spin_mtx, system, libMesh::PatchRecoveryErrorEstimator::target_patch_size, libMesh::System::time, libMesh::SystemNorm::type(), libMesh::DofMap::variable_type(), libMesh::W1_INF_SEMINORM, libMesh::W2_INF_SEMINORM, libMesh::MeshTools::weight(), libMesh::SystemNorm::weight(), libMesh::WeightedPatchRecoveryErrorEstimator::weight_functions, libMesh::SystemNorm::weight_sq(), and libMesh::zero.

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

Definition at line 114 of file weighted_patch_recovery_error_estimator.h.

Referenced by operator()().

◆ error_per_cell

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

Definition at line 115 of file weighted_patch_recovery_error_estimator.h.

Referenced by operator()().

◆ system

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

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

Definition at line 113 of file weighted_patch_recovery_error_estimator.h.

Referenced by operator()().


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