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 =
static_cast<Order> 320 (fe_type.order + elem->p_level());
329 fe->attach_quadrature_rule (qrule.get());
332 const std::vector<Real> & JxW = fe->get_JxW();
333 const std::vector<Point> & q_point = fe->get_xyz();
339 const std::vector<std::vector<Real>> * phi =
nullptr;
344 if (norm_type ==
L2 ||
347 phi = &(fe->get_phi());
349 const std::vector<std::vector<RealGradient>> * dphi =
nullptr;
355 dphi = &(fe->get_dphi());
357 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 358 const std::vector<std::vector<RealTensor>> * d2phi =
nullptr;
361 d2phi = &(fe->get_d2phi());
365 std::vector<dof_id_type> dof_indices;
369 unsigned int matsize = element_order + 1;
372 matsize *= (element_order + 2);
377 matsize *= (element_order + 3);
381 DenseMatrix<Number> Kp(matsize,matsize);
382 DenseVector<Number> F, Fx, Fy, Fz, Fxy, Fxz, Fyz;
383 DenseVector<Number> Pu_h, Pu_x_h, Pu_y_h, Pu_z_h, Pu_xy_h, Pu_xz_h, Pu_yz_h;
384 if (norm_type ==
L2 ||
387 F.resize(matsize); Pu_h.resize(matsize);
394 Fx.resize(matsize); Pu_x_h.resize(matsize);
396 Fy.resize(matsize); Pu_y_h.resize(matsize);
399 Fz.resize(matsize); Pu_z_h.resize(matsize);
404 Fx.resize(matsize); Pu_x_h.resize(matsize);
408 libmesh_assert_greater (LIBMESH_DIM, 1);
409 Fy.resize(matsize); Pu_y_h.resize(matsize);
413 libmesh_assert_greater (LIBMESH_DIM, 2);
414 Fz.resize(matsize); Pu_z_h.resize(matsize);
421 Fxy.resize(matsize); Pu_xy_h.resize(matsize);
423 Fxz.resize(matsize); Pu_xz_h.resize(matsize);
424 Fyz.resize(matsize); Pu_yz_h.resize(matsize);
432 for (
const auto & e_p : patch)
439 dof_map.dof_indices (e_p, dof_indices, var);
440 libmesh_assert_equal_to (dof_indices.size(), phi->size());
442 const unsigned int n_dofs =
443 cast_int<unsigned int>(dof_indices.size());
444 const unsigned int n_qp = qrule->n_points();
448 for (
unsigned int qp=0; qp<n_qp; qp++)
451 std::vector<Real> psi(
specpoly(
dim, element_order, q_point[qp], matsize));
453 const unsigned int psi_size = cast_int<unsigned int>(psi.size());
456 const unsigned int m = Kp.m(), n = Kp.n();
457 for (
unsigned int i=0; i<m; i++)
458 for (
unsigned int j=0; j<n; j++)
459 Kp(i,j) += JxW[qp]*psi[i]*psi[j];
461 if (norm_type ==
L2 ||
468 for (
unsigned int i=0; i<n_dofs; i++)
472 for (
unsigned int i=0; i != psi_size; i++)
473 F(i) += JxW[qp]*u_h*psi[i];
483 for (
unsigned int i=0; i<n_dofs; i++)
488 for (
unsigned int i=0; i != psi_size; i++)
490 Fx(i) += JxW[qp]*grad_u_h(0)*psi[i];
492 Fy(i) += JxW[qp]*grad_u_h(1)*psi[i];
495 Fz(i) += JxW[qp]*grad_u_h(2)*psi[i];
505 for (
unsigned int i=0; i<n_dofs; i++)
510 for (
unsigned int i=0; i != psi_size; i++)
512 Fx(i) += JxW[qp]*grad_u_h(0)*psi[i];
522 for (
unsigned int i=0; i<n_dofs; i++)
527 for (
unsigned int i=0; i != psi_size; i++)
529 Fy(i) += JxW[qp]*grad_u_h(1)*psi[i];
532 #endif // LIBMESH_DIM > 1 540 for (
unsigned int i=0; i<n_dofs; i++)
545 for (
unsigned int i=0; i != psi_size; i++)
547 Fz(i) += JxW[qp]*grad_u_h(2)*psi[i];
550 #endif // LIBMESH_DIM > 2 554 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 559 for (
unsigned int i=0; i<n_dofs; i++)
564 for (
unsigned int i=0; i != psi_size; i++)
566 Fx(i) += JxW[qp]*hess_u_h(0,0)*psi[i];
568 Fy(i) += JxW[qp]*hess_u_h(1,1)*psi[i];
569 Fxy(i) += JxW[qp]*hess_u_h(0,1)*psi[i];
572 Fz(i) += JxW[qp]*hess_u_h(2,2)*psi[i];
573 Fxz(i) += JxW[qp]*hess_u_h(0,2)*psi[i];
574 Fyz(i) += JxW[qp]*hess_u_h(1,2)*psi[i];
578 libmesh_error_msg(
"ERROR: --enable-second-derivatives is required \nfor _sobolev_order == 2!");
592 if (norm_type ==
L2 ||
595 Kp.lu_solve(F, Pu_h);
602 Kp.lu_solve (Fx, Pu_x_h);
604 Kp.lu_solve (Fy, Pu_y_h);
607 Kp.lu_solve (Fz, Pu_z_h);
612 Kp.lu_solve (Fx, Pu_x_h);
616 Kp.lu_solve (Fy, Pu_y_h);
620 Kp.lu_solve (Fz, Pu_z_h);
627 Kp.lu_solve(Fxy, Pu_xy_h);
629 Kp.lu_solve(Fxz, Pu_xz_h);
630 Kp.lu_solve(Fyz, Pu_yz_h);
639 Patch::const_iterator patch_re_it;
640 Patch::const_iterator patch_re_end;
643 Patch patch_re(
mesh.processor_id());
648 patch_re_it = patch.begin();
649 patch_re_end = patch.end();
655 patch_re.build_around_element (elem, 0,
659 patch_re_it = patch_re.begin();
660 patch_re_end = patch_re.end();
670 for (
unsigned int e = 0 ; patch_re_it != patch_re_end; ++patch_re_it, ++e)
675 const Elem * e_p = *patch_re_it;
693 dof_map.dof_indices (e_p, dof_indices, var);
694 libmesh_assert_equal_to (dof_indices.size(), phi->size());
697 const unsigned int n_dofs =
698 cast_int<unsigned int>(dof_indices.size());
701 Real element_error = 0;
704 static_cast<Order>(fe_type.order + e_p->p_level());
707 QGrid samprule (
dim, qorder);
711 fe->attach_quadrature_rule (&samprule);
714 const unsigned int n_sp =
715 cast_int<unsigned int>(JxW.size());
718 for (
unsigned int sp=0; sp<n_sp; sp++)
722 std::vector<Number> temperr(6,0.0);
724 if (norm_type ==
L2 ||
730 for (
unsigned int i=0; i<n_dofs; i++)
734 std::vector<Real> psi(
specpoly(
dim, element_order, q_point[sp], matsize));
735 for (
unsigned int i=0; i<matsize; i++)
737 temperr[0] += psi[i]*Pu_h(i);
748 for (
unsigned int i=0; i<n_dofs; i++)
753 std::vector<Real> psi(
specpoly(
dim, element_order, q_point[sp], matsize));
755 for (
unsigned int i=0; i<matsize; i++)
757 temperr[0] += psi[i]*Pu_x_h(i);
759 temperr[1] += psi[i]*Pu_y_h(i);
762 temperr[2] += psi[i]*Pu_z_h(i);
765 temperr[0] -= grad_u_h(0);
767 temperr[1] -= grad_u_h(1);
770 temperr[2] -= grad_u_h(2);
778 for (
unsigned int i=0; i<n_dofs; i++)
783 std::vector<Real> psi(
specpoly(
dim, element_order, q_point[sp], matsize));
784 for (
unsigned int i=0; i<matsize; i++)
786 temperr[0] += psi[i]*Pu_x_h(i);
789 temperr[0] -= grad_u_h(0);
797 for (
unsigned int i=0; i<n_dofs; i++)
802 std::vector<Real> psi(
specpoly(
dim, element_order, q_point[sp], matsize));
803 for (
unsigned int i=0; i<matsize; i++)
805 temperr[1] += psi[i]*Pu_y_h(i);
808 temperr[1] -= grad_u_h(1);
810 #endif // LIBMESH_DIM > 1 817 for (
unsigned int i=0; i<n_dofs; i++)
822 std::vector<Real> psi(
specpoly(
dim, element_order, q_point[sp], matsize));
823 for (
unsigned int i=0; i<matsize; i++)
825 temperr[2] += psi[i]*Pu_z_h(i);
828 temperr[2] -= grad_u_h(2);
830 #endif // LIBMESH_DIM > 2 834 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 838 for (
unsigned int i=0; i<n_dofs; i++)
843 std::vector<Real> psi(
specpoly(
dim, element_order, q_point[sp], matsize));
844 for (
unsigned int i=0; i<matsize; i++)
846 temperr[0] += psi[i]*Pu_x_h(i);
848 temperr[1] += psi[i]*Pu_y_h(i);
849 temperr[3] += psi[i]*Pu_xy_h(i);
852 temperr[2] += psi[i]*Pu_z_h(i);
853 temperr[4] += psi[i]*Pu_xz_h(i);
854 temperr[5] += psi[i]*Pu_yz_h(i);
858 temperr[0] -= hess_u_h(0,0);
860 temperr[1] -= hess_u_h(1,1);
861 temperr[3] -= hess_u_h(0,1);
864 temperr[2] -= hess_u_h(2,2);
865 temperr[4] -= hess_u_h(0,2);
866 temperr[5] -= hess_u_h(1,2);
869 libmesh_error_msg(
"ERROR: --enable-second-derivatives is required \nfor _sobolev_order == 2!");
876 if (norm_type ==
L_INF)
877 element_error = std::max(element_error,
std::abs(temperr[0]));
879 for (
unsigned int i=0; i != LIBMESH_DIM; ++i)
880 element_error = std::max(element_error,
std::abs(temperr[i]));
882 for (
unsigned int i=0; i != 6; ++i)
883 element_error = std::max(element_error,
std::abs(temperr[i]));
884 else if (norm_type ==
L2)
887 for (
unsigned int i=0; i != LIBMESH_DIM; ++i)
897 for (
unsigned int i=0; i != LIBMESH_DIM; ++i)
900 for (
unsigned int i=3; i != 6; ++i)
906 if (norm_type ==
L_INF ||
910 else if (norm_type ==
L2 ||
929 Patch::const_iterator patch_re_it;
930 Patch::const_iterator patch_re_end;
933 Patch current_elem_patch(
mesh.processor_id());
938 patch_re_it = patch.begin();
939 patch_re_end = patch.end();
945 current_elem_patch.build_around_element (elem, 0,
949 patch_re_it = current_elem_patch.begin();
950 patch_re_end = current_elem_patch.end();
954 for (
unsigned int i = 0 ; patch_re_it != patch_re_end; ++patch_re_it, ++i)
957 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...
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
const MeshBase & get_mesh() const
DIE A HORRIBLE DEATH HERE typedef float ErrorVectorReal
Number current_solution(const dof_id_type global_dof_number) const
ADRealEigenVector< T, D, asd > abs(const ADRealEigenVector< T, D, asd > &)
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.