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

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

Public Member Functions

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

Private Attributes

const Systemsystem
 
const PatchRecoveryErrorEstimatorerror_estimator
 
ErrorVectorerror_per_cell
 

Detailed Description

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

May be executed in parallel on separate threads.

Definition at line 118 of file patch_recovery_error_estimator.h.

Constructor & Destructor Documentation

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

Member Function Documentation

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

Definition at line 188 of file patch_recovery_error_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::ParallelObject::processor_id(), libMesh::Real, libMesh::DenseVector< T >::resize(), libMesh::PatchRecoveryErrorEstimator::specpoly(), libMesh::Threads::spin_mtx, libMesh::DofMap::variable_type(), libMesh::W1_INF_SEMINORM, libMesh::W2_INF_SEMINORM, and libMesh::zero.

Referenced by EstimateError().

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

Definition at line 133 of file patch_recovery_error_estimator.h.

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

Definition at line 134 of file patch_recovery_error_estimator.h.

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

Definition at line 132 of file patch_recovery_error_estimator.h.


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