118 const unsigned int dim =
mesh.mesh_dimension();
128 for (
const auto & elem : range)
135 Patch patch(
mesh.processor_id());
153 std::vector<Real> new_error_per_cell(1, 0.);
155 new_error_per_cell.resize(patch.size(), 0.);
159 for (
unsigned int var=0; var<
n_vars; var++)
162 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 164 bool is_valid_norm_type =
171 norm_type ==
L_INF ||
178 norm_type ==
L_INF ||
190 bool is_valid_norm_combo =
203 ((norm_type ==
L_INF ||
217 const FEType & fe_type = dof_map.variable_type (var);
219 const Order element_order =
static_cast<Order> 220 (fe_type.order + elem->p_level());
226 std::unique_ptr<QBase> qrule =
230 fe->attach_quadrature_rule (qrule.get());
233 const std::vector<Real> & JxW = fe->get_JxW();
234 const std::vector<Point> & q_point = fe->get_xyz();
240 const std::vector<std::vector<Real>> * phi =
nullptr;
245 if (norm_type ==
L2 ||
248 phi = &(fe->get_phi());
250 const std::vector<std::vector<RealGradient>> * dphi =
nullptr;
256 dphi = &(fe->get_dphi());
258 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 259 const std::vector<std::vector<RealTensor>> * d2phi =
nullptr;
262 d2phi = &(fe->get_d2phi());
266 std::vector<dof_id_type> dof_indices;
270 unsigned int matsize = element_order + 1;
273 matsize *= (element_order + 2);
278 matsize *= (element_order + 3);
282 DenseMatrix<Number> Kp(matsize,matsize);
283 DenseVector<Number> F, Fx, Fy, Fz, Fxy, Fxz, Fyz;
284 DenseVector<Number> Pu_h, Pu_x_h, Pu_y_h, Pu_z_h, Pu_xy_h, Pu_xz_h, Pu_yz_h;
285 if (norm_type ==
L2 ||
288 F.resize(matsize); Pu_h.resize(matsize);
295 Fx.resize(matsize); Pu_x_h.resize(matsize);
297 Fy.resize(matsize); Pu_y_h.resize(matsize);
300 Fz.resize(matsize); Pu_z_h.resize(matsize);
305 Fx.resize(matsize); Pu_x_h.resize(matsize);
310 Fy.resize(matsize); Pu_y_h.resize(matsize);
315 Fz.resize(matsize); Pu_z_h.resize(matsize);
322 Fxy.resize(matsize); Pu_xy_h.resize(matsize);
324 Fxz.resize(matsize); Pu_xz_h.resize(matsize);
325 Fyz.resize(matsize); Pu_yz_h.resize(matsize);
334 for (
const auto & e_p : patch)
341 dof_map.dof_indices (e_p, dof_indices, var);
344 const unsigned int n_dofs =
345 cast_int<unsigned int>(dof_indices.size());
346 const unsigned int n_qp = qrule->n_points();
350 for (
unsigned int qp=0; qp<n_qp; qp++)
353 std::vector<Real> psi(
specpoly(
dim, element_order, q_point[qp], matsize));
355 const unsigned int psi_size = cast_int<unsigned int>(psi.size());
358 const unsigned int m = Kp.m(), n = Kp.n();
359 for (
unsigned int i=0; i<m; i++)
360 for (
unsigned int j=0; j<n; j++)
361 Kp(i,j) += JxW[qp]*psi[i]*psi[j];
363 if (norm_type ==
L2 ||
370 for (
unsigned int i=0; i<n_dofs; i++)
374 for (
unsigned int i=0; i != psi_size; i++)
375 F(i) = JxW[qp]*u_h*psi[i];
385 for (std::size_t i=0; i<n_dofs; i++)
392 for (
unsigned int i=0; i != psi_size; i++)
394 Fx(i) += JxW[qp]*grad_u_h(0)*psi[i];
396 Fy(i) += JxW[qp]*grad_u_h(1)*psi[i];
399 Fz(i) += JxW[qp]*grad_u_h(2)*psi[i];
409 for (
unsigned int i=0; i<n_dofs; i++)
416 for (
unsigned int i=0; i != psi_size; i++)
418 Fx(i) += JxW[qp]*grad_u_h(0)*psi[i];
428 for (
unsigned int i=0; i<n_dofs; i++)
435 for (
unsigned int i=0; i != psi_size; i++)
437 Fy(i) += JxW[qp]*grad_u_h(1)*psi[i];
440 #endif // LIBMESH_DIM > 1 448 for (
unsigned int i=0; i<n_dofs; i++)
455 for (
unsigned int i=0; i != psi_size; i++)
457 Fz(i) += JxW[qp]*grad_u_h(2)*psi[i];
460 #endif // LIBMESH_DIM > 2 464 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 469 for (
unsigned int i=0; i<n_dofs; i++)
476 for (
unsigned int i=0; i != psi_size; i++)
478 Fx(i) += JxW[qp]*hess_u_h(0,0)*psi[i];
480 Fy(i) += JxW[qp]*hess_u_h(1,1)*psi[i];
481 Fxy(i) += JxW[qp]*hess_u_h(0,1)*psi[i];
484 Fz(i) += JxW[qp]*hess_u_h(2,2)*psi[i];
485 Fxz(i) += JxW[qp]*hess_u_h(0,2)*psi[i];
486 Fyz(i) += JxW[qp]*hess_u_h(1,2)*psi[i];
490 libmesh_error_msg(
"ERROR: --enable-second-derivatives is required \nfor _sobolev_order == 2!");
504 if (norm_type ==
L2 ||
507 Kp.lu_solve(F, Pu_h);
514 Kp.lu_solve (Fx, Pu_x_h);
516 Kp.lu_solve (Fy, Pu_y_h);
519 Kp.lu_solve (Fz, Pu_z_h);
524 Kp.lu_solve (Fx, Pu_x_h);
528 Kp.lu_solve (Fy, Pu_y_h);
532 Kp.lu_solve (Fz, Pu_z_h);
539 Kp.lu_solve(Fxy, Pu_xy_h);
541 Kp.lu_solve(Fxz, Pu_xz_h);
542 Kp.lu_solve(Fyz, Pu_yz_h);
551 Patch::const_iterator patch_re_it;
552 Patch::const_iterator patch_re_end;
555 Patch patch_re(
mesh.processor_id());
560 patch_re_it = patch.begin();
561 patch_re_end = patch.end();
567 patch_re.build_around_element (elem, 0,
571 patch_re_it = patch_re.begin();
572 patch_re_end = patch_re.end();
584 FEMContext femcontext(
system,
nullptr,
589 for (
unsigned int e = 0 ; patch_re_it != patch_re_end; ++patch_re_it, ++e)
594 const Elem * e_p = *patch_re_it;
600 femcontext.pre_fe_reinit(
system, e_p);
615 dof_map.dof_indices (e_p, dof_indices, var);
619 const unsigned int n_dofs =
620 cast_int<unsigned int>(dof_indices.size());
623 Real element_error = 0;
626 static_cast<Order>(fe_type.order + e_p->p_level());
629 QGrid samprule (
dim, qorder);
633 fe->attach_quadrature_rule (&samprule);
636 const unsigned int n_sp =
637 cast_int<unsigned int>(JxW.size());
640 for (
unsigned int sp=0; sp<n_sp; sp++)
644 std::vector<Number> temperr(6,0.0);
646 if (norm_type ==
L2 ||
652 for (
unsigned int i=0; i<n_dofs; i++)
656 std::vector<Real> psi(
specpoly(
dim, element_order, q_point[sp], matsize));
657 for (
unsigned int i=0; i<matsize; i++)
659 temperr[0] += psi[i]*Pu_h(i);
670 for (
unsigned int i=0; i<n_dofs; i++)
675 std::vector<Real> psi(
specpoly(
dim, element_order, q_point[sp], matsize));
677 for (
unsigned int i=0; i<matsize; i++)
679 temperr[0] += psi[i]*Pu_x_h(i);
681 temperr[1] += psi[i]*Pu_y_h(i);
684 temperr[2] += psi[i]*Pu_z_h(i);
687 temperr[0] -= grad_u_h(0);
689 temperr[1] -= grad_u_h(1);
692 temperr[2] -= grad_u_h(2);
700 for (
unsigned int i=0; i<n_dofs; i++)
705 std::vector<Real> psi(
specpoly(
dim, element_order, q_point[sp], matsize));
706 for (
unsigned int i=0; i<matsize; i++)
708 temperr[0] += psi[i]*Pu_x_h(i);
711 temperr[0] -= grad_u_h(0);
719 for (
unsigned int i=0; i<n_dofs; i++)
724 std::vector<Real> psi(
specpoly(
dim, element_order, q_point[sp], matsize));
725 for (
unsigned int i=0; i<matsize; i++)
727 temperr[1] += psi[i]*Pu_y_h(i);
730 temperr[1] -= grad_u_h(1);
732 #endif // LIBMESH_DIM > 1 739 for (
unsigned int i=0; i<n_dofs; i++)
744 std::vector<Real> psi(
specpoly(
dim, element_order, q_point[sp], matsize));
745 for (
unsigned int i=0; i<matsize; i++)
747 temperr[2] += psi[i]*Pu_z_h(i);
750 temperr[2] -= grad_u_h(2);
752 #endif // LIBMESH_DIM > 2 756 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 760 for (
unsigned int i=0; i<n_dofs; i++)
765 std::vector<Real> psi(
specpoly(
dim, element_order, q_point[sp], matsize));
766 for (
unsigned int i=0; i<matsize; i++)
768 temperr[0] += psi[i]*Pu_x_h(i);
770 temperr[1] += psi[i]*Pu_y_h(i);
771 temperr[3] += psi[i]*Pu_xy_h(i);
774 temperr[2] += psi[i]*Pu_z_h(i);
775 temperr[4] += psi[i]*Pu_xz_h(i);
776 temperr[5] += psi[i]*Pu_yz_h(i);
780 temperr[0] -= hess_u_h(0,0);
782 temperr[1] -= hess_u_h(1,1);
783 temperr[3] -= hess_u_h(0,1);
786 temperr[2] -= hess_u_h(2,2);
787 temperr[4] -= hess_u_h(0,2);
788 temperr[5] -= hess_u_h(1,2);
791 libmesh_error_msg(
"ERROR: --enable-second-derivatives is required \nfor _sobolev_order == 2!");
802 if (norm_type ==
L_INF)
803 element_error = std::max(element_error,
std::abs(
weight*temperr[0]));
805 for (
unsigned int i=0; i != LIBMESH_DIM; ++i)
806 element_error = std::max(element_error,
std::abs(
weight*temperr[i]));
808 for (
unsigned int i=0; i != 6; ++i)
809 element_error = std::max(element_error,
std::abs(
weight*temperr[i]));
810 else if (norm_type ==
L2)
813 for (
unsigned int i=0; i != LIBMESH_DIM; ++i)
823 for (
unsigned int i=0; i != LIBMESH_DIM; ++i)
826 for (
unsigned int i=3; i != 6; ++i)
832 if (norm_type ==
L_INF ||
836 else if (norm_type ==
L2 ||
855 Patch::const_iterator patch_re_it;
856 Patch::const_iterator patch_re_end;
859 Patch current_elem_patch(
mesh.processor_id());
864 patch_re_it = patch.begin();
865 patch_re_end = patch.end();
871 current_elem_patch.build_around_element (elem, 0,
875 patch_re_it = current_elem_patch.begin();
876 patch_re_end = current_elem_patch.end();
880 for (
unsigned int i = 0 ; patch_re_it != patch_re_end; ++patch_re_it, ++i)
883 const Elem * e_p = *patch_re_it;
909 (new_error_per_cell[i]);
Real time
For time-dependent problems, this is the time t at the beginning of the current timestep.
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 System & system
Function to set the boolean patch_reuse in case the user wants to change the default behaviour of pat...
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...
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
Real weight(unsigned int var) const
ErrorVector & error_per_cell
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const WeightedPatchRecoveryErrorEstimator & error_estimator
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
const DofMap & get_dof_map() const
spin_mutex spin_mtx
A convenient spin mutex object which can be used for obtaining locks.