217 const unsigned int dim =
mesh.mesh_dimension();
227 for (
const auto & elem : range)
234 Patch patch(
mesh.processor_id());
252 std::vector<Real> new_error_per_cell(1, 0.);
254 new_error_per_cell.resize(patch.size(), 0.);
258 for (
unsigned int var=0; var<
n_vars; var++)
261 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 263 bool is_valid_norm_type =
270 norm_type ==
L_INF ||
277 norm_type ==
L_INF ||
290 bool is_valid_norm_combo =
303 ((norm_type ==
L_INF ||
317 const FEType & fe_type = dof_map.variable_type (var);
319 const Order element_order = fe_type.order + elem->p_level();
328 fe->attach_quadrature_rule (qrule.get());
331 const std::vector<Real> & JxW = fe->get_JxW();
332 const std::vector<Point> & q_point = fe->get_xyz();
338 const std::vector<std::vector<Real>> * phi =
nullptr;
343 if (norm_type ==
L2 ||
346 phi = &(fe->get_phi());
348 const std::vector<std::vector<RealGradient>> * dphi =
nullptr;
354 dphi = &(fe->get_dphi());
356 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 357 const std::vector<std::vector<RealTensor>> * d2phi =
nullptr;
360 d2phi = &(fe->get_d2phi());
364 std::vector<dof_id_type> dof_indices;
368 unsigned int matsize = element_order + 1;
371 matsize *= (element_order + 2);
376 matsize *= (element_order + 3);
380 DenseMatrix<Number> Kp(matsize,matsize);
381 DenseVector<Number> F, Fx, Fy, Fz, Fxy, Fxz, Fyz;
382 DenseVector<Number> Pu_h, Pu_x_h, Pu_y_h, Pu_z_h, Pu_xy_h, Pu_xz_h, Pu_yz_h;
383 if (norm_type ==
L2 ||
386 F.resize(matsize); Pu_h.resize(matsize);
393 Fx.resize(matsize); Pu_x_h.resize(matsize);
395 Fy.resize(matsize); Pu_y_h.resize(matsize);
398 Fz.resize(matsize); Pu_z_h.resize(matsize);
403 Fx.resize(matsize); Pu_x_h.resize(matsize);
407 libmesh_assert_greater (LIBMESH_DIM, 1);
408 Fy.resize(matsize); Pu_y_h.resize(matsize);
412 libmesh_assert_greater (LIBMESH_DIM, 2);
413 Fz.resize(matsize); Pu_z_h.resize(matsize);
420 Fxy.resize(matsize); Pu_xy_h.resize(matsize);
422 Fxz.resize(matsize); Pu_xz_h.resize(matsize);
423 Fyz.resize(matsize); Pu_yz_h.resize(matsize);
431 for (
const auto & e_p : patch)
438 dof_map.dof_indices (e_p, dof_indices, var);
439 libmesh_assert_equal_to (dof_indices.size(), phi->size());
441 const unsigned int n_dofs =
442 cast_int<unsigned int>(dof_indices.size());
443 const unsigned int n_qp = qrule->n_points();
447 for (
unsigned int qp=0; qp<n_qp; qp++)
450 std::vector<Real> psi(
specpoly(
dim, element_order, q_point[qp], matsize));
452 const unsigned int psi_size = cast_int<unsigned int>(psi.size());
455 const unsigned int m = Kp.m(), n = Kp.n();
456 for (
unsigned int i=0; i<m; i++)
457 for (
unsigned int j=0; j<n; j++)
458 Kp(i,j) += JxW[qp]*psi[i]*psi[j];
460 if (norm_type ==
L2 ||
467 for (
unsigned int i=0; i<n_dofs; i++)
471 for (
unsigned int i=0; i != psi_size; i++)
472 F(i) += JxW[qp]*u_h*psi[i];
482 for (
unsigned int i=0; i<n_dofs; i++)
487 for (
unsigned int i=0; i != psi_size; i++)
489 Fx(i) += JxW[qp]*grad_u_h(0)*psi[i];
491 Fy(i) += JxW[qp]*grad_u_h(1)*psi[i];
494 Fz(i) += JxW[qp]*grad_u_h(2)*psi[i];
504 for (
unsigned int i=0; i<n_dofs; i++)
509 for (
unsigned int i=0; i != psi_size; i++)
511 Fx(i) += JxW[qp]*grad_u_h(0)*psi[i];
521 for (
unsigned int i=0; i<n_dofs; i++)
526 for (
unsigned int i=0; i != psi_size; i++)
528 Fy(i) += JxW[qp]*grad_u_h(1)*psi[i];
531 #endif // LIBMESH_DIM > 1 539 for (
unsigned int i=0; i<n_dofs; i++)
544 for (
unsigned int i=0; i != psi_size; i++)
546 Fz(i) += JxW[qp]*grad_u_h(2)*psi[i];
549 #endif // LIBMESH_DIM > 2 553 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 558 for (
unsigned int i=0; i<n_dofs; i++)
563 for (
unsigned int i=0; i != psi_size; i++)
565 Fx(i) += JxW[qp]*hess_u_h(0,0)*psi[i];
567 Fy(i) += JxW[qp]*hess_u_h(1,1)*psi[i];
568 Fxy(i) += JxW[qp]*hess_u_h(0,1)*psi[i];
571 Fz(i) += JxW[qp]*hess_u_h(2,2)*psi[i];
572 Fxz(i) += JxW[qp]*hess_u_h(0,2)*psi[i];
573 Fyz(i) += JxW[qp]*hess_u_h(1,2)*psi[i];
577 libmesh_error_msg(
"ERROR: --enable-second-derivatives is required \nfor _sobolev_order == 2!");
591 if (norm_type ==
L2 ||
594 Kp.lu_solve(F, Pu_h);
601 Kp.lu_solve (Fx, Pu_x_h);
603 Kp.lu_solve (Fy, Pu_y_h);
606 Kp.lu_solve (Fz, Pu_z_h);
611 Kp.lu_solve (Fx, Pu_x_h);
615 Kp.lu_solve (Fy, Pu_y_h);
619 Kp.lu_solve (Fz, Pu_z_h);
626 Kp.lu_solve(Fxy, Pu_xy_h);
628 Kp.lu_solve(Fxz, Pu_xz_h);
629 Kp.lu_solve(Fyz, Pu_yz_h);
638 Patch::const_iterator patch_re_it;
639 Patch::const_iterator patch_re_end;
642 Patch patch_re(
mesh.processor_id());
647 patch_re_it = patch.begin();
648 patch_re_end = patch.end();
654 patch_re.build_around_element (elem, 0,
658 patch_re_it = patch_re.begin();
659 patch_re_end = patch_re.end();
669 for (
unsigned int e = 0 ; patch_re_it != patch_re_end; ++patch_re_it, ++e)
674 const Elem * e_p = *patch_re_it;
692 dof_map.dof_indices (e_p, dof_indices, var);
693 libmesh_assert_equal_to (dof_indices.size(), phi->size());
696 const unsigned int n_dofs =
697 cast_int<unsigned int>(dof_indices.size());
700 Real element_error = 0;
702 const Order qorder = fe_type.order + e_p->p_level();
705 QGrid samprule (
dim, qorder);
709 fe->attach_quadrature_rule (&samprule);
712 const unsigned int n_sp =
713 cast_int<unsigned int>(JxW.size());
716 for (
unsigned int sp=0; sp<n_sp; sp++)
720 std::vector<Number> temperr(6,0.0);
722 if (norm_type ==
L2 ||
728 for (
unsigned int i=0; i<n_dofs; i++)
732 std::vector<Real> psi(
specpoly(
dim, element_order, q_point[sp], matsize));
733 for (
unsigned int i=0; i<matsize; i++)
735 temperr[0] += psi[i]*Pu_h(i);
746 for (
unsigned int i=0; i<n_dofs; i++)
751 std::vector<Real> psi(
specpoly(
dim, element_order, q_point[sp], matsize));
753 for (
unsigned int i=0; i<matsize; i++)
755 temperr[0] += psi[i]*Pu_x_h(i);
757 temperr[1] += psi[i]*Pu_y_h(i);
760 temperr[2] += psi[i]*Pu_z_h(i);
763 temperr[0] -= grad_u_h(0);
765 temperr[1] -= grad_u_h(1);
768 temperr[2] -= grad_u_h(2);
776 for (
unsigned int i=0; i<n_dofs; i++)
781 std::vector<Real> psi(
specpoly(
dim, element_order, q_point[sp], matsize));
782 for (
unsigned int i=0; i<matsize; i++)
784 temperr[0] += psi[i]*Pu_x_h(i);
787 temperr[0] -= grad_u_h(0);
795 for (
unsigned int i=0; i<n_dofs; i++)
800 std::vector<Real> psi(
specpoly(
dim, element_order, q_point[sp], matsize));
801 for (
unsigned int i=0; i<matsize; i++)
803 temperr[1] += psi[i]*Pu_y_h(i);
806 temperr[1] -= grad_u_h(1);
808 #endif // LIBMESH_DIM > 1 815 for (
unsigned int i=0; i<n_dofs; i++)
820 std::vector<Real> psi(
specpoly(
dim, element_order, q_point[sp], matsize));
821 for (
unsigned int i=0; i<matsize; i++)
823 temperr[2] += psi[i]*Pu_z_h(i);
826 temperr[2] -= grad_u_h(2);
828 #endif // LIBMESH_DIM > 2 832 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 836 for (
unsigned int i=0; i<n_dofs; i++)
841 std::vector<Real> psi(
specpoly(
dim, element_order, q_point[sp], matsize));
842 for (
unsigned int i=0; i<matsize; i++)
844 temperr[0] += psi[i]*Pu_x_h(i);
846 temperr[1] += psi[i]*Pu_y_h(i);
847 temperr[3] += psi[i]*Pu_xy_h(i);
850 temperr[2] += psi[i]*Pu_z_h(i);
851 temperr[4] += psi[i]*Pu_xz_h(i);
852 temperr[5] += psi[i]*Pu_yz_h(i);
856 temperr[0] -= hess_u_h(0,0);
858 temperr[1] -= hess_u_h(1,1);
859 temperr[3] -= hess_u_h(0,1);
862 temperr[2] -= hess_u_h(2,2);
863 temperr[4] -= hess_u_h(0,2);
864 temperr[5] -= hess_u_h(1,2);
867 libmesh_error_msg(
"ERROR: --enable-second-derivatives is required \nfor _sobolev_order == 2!");
874 if (norm_type ==
L_INF)
875 element_error = std::max(element_error, std::abs(temperr[0]));
877 for (
unsigned int i=0; i != LIBMESH_DIM; ++i)
878 element_error = std::max(element_error, std::abs(temperr[i]));
880 for (
unsigned int i=0; i != 6; ++i)
881 element_error = std::max(element_error, std::abs(temperr[i]));
882 else if (norm_type ==
L2)
885 for (
unsigned int i=0; i != LIBMESH_DIM; ++i)
895 for (
unsigned int i=0; i != LIBMESH_DIM; ++i)
898 for (
unsigned int i=3; i != 6; ++i)
904 if (norm_type ==
L_INF ||
908 else if (norm_type ==
L2 ||
927 Patch::const_iterator patch_re_it;
928 Patch::const_iterator patch_re_end;
931 Patch current_elem_patch(
mesh.processor_id());
936 patch_re_it = patch.begin();
937 patch_re_end = patch.end();
943 current_elem_patch.build_around_element (elem, 0,
947 patch_re_it = current_elem_patch.begin();
948 patch_re_end = current_elem_patch.end();
952 for (
unsigned int i = 0 ; patch_re_it != patch_re_end; ++patch_re_it, ++i)
955 const Elem * e_p = *patch_re_it;
void add_scaled(const TypeVector< T2 > &, const T &)
Add a scaled value to this vector without creating a temporary.
Order
defines an enum for polynomial orders.
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...
const MeshBase & get_mesh() const
DIE A HORRIBLE DEATH HERE typedef float ErrorVectorReal
Number current_solution(const dof_id_type global_dof_number) const
FEMNormType type(unsigned int var) const
NumberVectorValue Gradient
static std::unique_ptr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.
unsigned int target_patch_size
The PatchErrorEstimator will build patches of at least this many elements to perform estimates...
const PatchRecoveryErrorEstimator & error_estimator
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.
std::string enum_to_string(const T e)
Real weight_sq(unsigned int var) const
ErrorVector & error_per_cell
Real weight(unsigned int var) const
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
int _extra_order
Extra order to use for quadrature rule.
unsigned int n_vars() const
const DofMap & get_dof_map() const
spin_mutex spin_mtx
A convenient spin mutex object which can be used for obtaining locks.