libMesh::PatchRecoveryErrorEstimator::EstimateError Class Reference

Public Member Functions

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

Private Attributes

const 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 115 of file patch_recovery_error_estimator.h.

Constructor & Destructor Documentation

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

Definition at line 118 of file patch_recovery_error_estimator.h.

120  :
121  system(sys),
122  error_estimator(ee),
123  error_per_cell(epc)
124  {}

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< T >::build(), libMesh::Patch::build_around_element(), libMesh::System::current_solution(), libMesh::FEType::default_quadrature_rule(), libMesh::dim, libMesh::DofMap::dof_indices(), libMesh::dof_map, libMesh::StoredRange< iterator_type, object_type >::end(), libMesh::err, error_estimator, libMesh::ErrorEstimator::error_norm, error_per_cell, libMesh::ErrorVectorReal, libMesh::System::get_dof_map(), libMesh::System::get_mesh(), libMeshEnums::H1_SEMINORM, libMeshEnums::H1_X_SEMINORM, libMeshEnums::H1_Y_SEMINORM, libMeshEnums::H1_Z_SEMINORM, libMeshEnums::H2_SEMINORM, libMesh::DofObject::id(), libMeshEnums::L2, libMeshEnums::L_INF, libMesh::libmesh_assert(), libMesh::libmesh_assert_greater(), libMesh::DenseMatrix< T >::lu_solve(), libMesh::DenseMatrixBase< T >::m(), std::max(), mesh, libMesh::MeshBase::mesh_dimension(), libMesh::DenseMatrixBase< T >::n(), libMesh::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::ParallelObject::processor_id(), libMesh::Real, libMesh::DenseVector< T >::resize(), libMesh::PatchRecoveryErrorEstimator::specpoly(), libMesh::Threads::spin_mtx, system, libMesh::PatchRecoveryErrorEstimator::target_patch_size, libMesh::SystemNorm::type(), libMesh::DofMap::variable_type(), libMeshEnums::W1_INF_SEMINORM, libMeshEnums::W2_INF_SEMINORM, libMesh::SystemNorm::weight(), libMesh::SystemNorm::weight_sq(), and libMesh::zero.

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

Member Data Documentation

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

Definition at line 136 of file patch_recovery_error_estimator.h.

Referenced by operator()().

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

Definition at line 137 of file patch_recovery_error_estimator.h.

Referenced by operator()().

const System& libMesh::PatchRecoveryErrorEstimator::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 135 of file patch_recovery_error_estimator.h.

Referenced by operator()().


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

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:29 UTC

Hosted By:
SourceForge.net Logo