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 = fe_type.order + elem->p_level();
225 std::unique_ptr<QBase> qrule =
229 fe->attach_quadrature_rule (qrule.get());
232 const std::vector<Real> & JxW = fe->get_JxW();
233 const std::vector<Point> & q_point = fe->get_xyz();
239 const std::vector<std::vector<Real>> * phi =
nullptr;
244 if (norm_type ==
L2 ||
247 phi = &(fe->get_phi());
249 const std::vector<std::vector<RealGradient>> * dphi =
nullptr;
255 dphi = &(fe->get_dphi());
257 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 258 const std::vector<std::vector<RealTensor>> * d2phi =
nullptr;
261 d2phi = &(fe->get_d2phi());
265 std::vector<dof_id_type> dof_indices;
269 unsigned int matsize = element_order + 1;
272 matsize *= (element_order + 2);
277 matsize *= (element_order + 3);
281 DenseMatrix<Number> Kp(matsize,matsize);
282 DenseVector<Number> F, Fx, Fy, Fz, Fxy, Fxz, Fyz;
283 DenseVector<Number> Pu_h, Pu_x_h, Pu_y_h, Pu_z_h, Pu_xy_h, Pu_xz_h, Pu_yz_h;
284 if (norm_type ==
L2 ||
287 F.resize(matsize); Pu_h.resize(matsize);
294 Fx.resize(matsize); Pu_x_h.resize(matsize);
296 Fy.resize(matsize); Pu_y_h.resize(matsize);
299 Fz.resize(matsize); Pu_z_h.resize(matsize);
304 Fx.resize(matsize); Pu_x_h.resize(matsize);
309 Fy.resize(matsize); Pu_y_h.resize(matsize);
314 Fz.resize(matsize); Pu_z_h.resize(matsize);
321 Fxy.resize(matsize); Pu_xy_h.resize(matsize);
323 Fxz.resize(matsize); Pu_xz_h.resize(matsize);
324 Fyz.resize(matsize); Pu_yz_h.resize(matsize);
333 for (
const auto & e_p : patch)
340 dof_map.dof_indices (e_p, dof_indices, var);
343 const unsigned int n_dofs =
344 cast_int<unsigned int>(dof_indices.size());
345 const unsigned int n_qp = qrule->n_points();
349 for (
unsigned int qp=0; qp<n_qp; qp++)
352 std::vector<Real> psi(
specpoly(
dim, element_order, q_point[qp], matsize));
354 const unsigned int psi_size = cast_int<unsigned int>(psi.size());
357 const unsigned int m = Kp.m(), n = Kp.n();
358 for (
unsigned int i=0; i<m; i++)
359 for (
unsigned int j=0; j<n; j++)
360 Kp(i,j) += JxW[qp]*psi[i]*psi[j];
362 if (norm_type ==
L2 ||
369 for (
unsigned int i=0; i<n_dofs; i++)
373 for (
unsigned int i=0; i != psi_size; i++)
374 F(i) = JxW[qp]*u_h*psi[i];
384 for (std::size_t i=0; i<n_dofs; i++)
391 for (
unsigned int i=0; i != psi_size; i++)
393 Fx(i) += JxW[qp]*grad_u_h(0)*psi[i];
395 Fy(i) += JxW[qp]*grad_u_h(1)*psi[i];
398 Fz(i) += JxW[qp]*grad_u_h(2)*psi[i];
408 for (
unsigned int i=0; i<n_dofs; i++)
415 for (
unsigned int i=0; i != psi_size; i++)
417 Fx(i) += JxW[qp]*grad_u_h(0)*psi[i];
427 for (
unsigned int i=0; i<n_dofs; i++)
434 for (
unsigned int i=0; i != psi_size; i++)
436 Fy(i) += JxW[qp]*grad_u_h(1)*psi[i];
439 #endif // LIBMESH_DIM > 1 447 for (
unsigned int i=0; i<n_dofs; i++)
454 for (
unsigned int i=0; i != psi_size; i++)
456 Fz(i) += JxW[qp]*grad_u_h(2)*psi[i];
459 #endif // LIBMESH_DIM > 2 463 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 468 for (
unsigned int i=0; i<n_dofs; i++)
475 for (
unsigned int i=0; i != psi_size; i++)
477 Fx(i) += JxW[qp]*hess_u_h(0,0)*psi[i];
479 Fy(i) += JxW[qp]*hess_u_h(1,1)*psi[i];
480 Fxy(i) += JxW[qp]*hess_u_h(0,1)*psi[i];
483 Fz(i) += JxW[qp]*hess_u_h(2,2)*psi[i];
484 Fxz(i) += JxW[qp]*hess_u_h(0,2)*psi[i];
485 Fyz(i) += JxW[qp]*hess_u_h(1,2)*psi[i];
489 libmesh_error_msg(
"ERROR: --enable-second-derivatives is required \nfor _sobolev_order == 2!");
503 if (norm_type ==
L2 ||
506 Kp.lu_solve(F, Pu_h);
513 Kp.lu_solve (Fx, Pu_x_h);
515 Kp.lu_solve (Fy, Pu_y_h);
518 Kp.lu_solve (Fz, Pu_z_h);
523 Kp.lu_solve (Fx, Pu_x_h);
527 Kp.lu_solve (Fy, Pu_y_h);
531 Kp.lu_solve (Fz, Pu_z_h);
538 Kp.lu_solve(Fxy, Pu_xy_h);
540 Kp.lu_solve(Fxz, Pu_xz_h);
541 Kp.lu_solve(Fyz, Pu_yz_h);
550 Patch::const_iterator patch_re_it;
551 Patch::const_iterator patch_re_end;
554 Patch patch_re(
mesh.processor_id());
559 patch_re_it = patch.begin();
560 patch_re_end = patch.end();
566 patch_re.build_around_element (elem, 0,
570 patch_re_it = patch_re.begin();
571 patch_re_end = patch_re.end();
583 FEMContext femcontext(
system,
nullptr,
588 for (
unsigned int e = 0 ; patch_re_it != patch_re_end; ++patch_re_it, ++e)
593 const Elem * e_p = *patch_re_it;
599 femcontext.pre_fe_reinit(
system, e_p);
614 dof_map.dof_indices (e_p, dof_indices, var);
618 const unsigned int n_dofs =
619 cast_int<unsigned int>(dof_indices.size());
622 Real element_error = 0;
624 const Order qorder = fe_type.order + e_p->p_level();
627 QGrid samprule (
dim, qorder);
631 fe->attach_quadrature_rule (&samprule);
634 const unsigned int n_sp =
635 cast_int<unsigned int>(JxW.size());
638 for (
unsigned int sp=0; sp<n_sp; sp++)
642 std::vector<Number> temperr(6,0.0);
644 if (norm_type ==
L2 ||
650 for (
unsigned int i=0; i<n_dofs; i++)
654 std::vector<Real> psi(
specpoly(
dim, element_order, q_point[sp], matsize));
655 for (
unsigned int i=0; i<matsize; i++)
657 temperr[0] += psi[i]*Pu_h(i);
668 for (
unsigned int i=0; i<n_dofs; i++)
673 std::vector<Real> psi(
specpoly(
dim, element_order, q_point[sp], matsize));
675 for (
unsigned int i=0; i<matsize; i++)
677 temperr[0] += psi[i]*Pu_x_h(i);
679 temperr[1] += psi[i]*Pu_y_h(i);
682 temperr[2] += psi[i]*Pu_z_h(i);
685 temperr[0] -= grad_u_h(0);
687 temperr[1] -= grad_u_h(1);
690 temperr[2] -= grad_u_h(2);
698 for (
unsigned int i=0; i<n_dofs; i++)
703 std::vector<Real> psi(
specpoly(
dim, element_order, q_point[sp], matsize));
704 for (
unsigned int i=0; i<matsize; i++)
706 temperr[0] += psi[i]*Pu_x_h(i);
709 temperr[0] -= grad_u_h(0);
717 for (
unsigned int i=0; i<n_dofs; i++)
722 std::vector<Real> psi(
specpoly(
dim, element_order, q_point[sp], matsize));
723 for (
unsigned int i=0; i<matsize; i++)
725 temperr[1] += psi[i]*Pu_y_h(i);
728 temperr[1] -= grad_u_h(1);
730 #endif // LIBMESH_DIM > 1 737 for (
unsigned int i=0; i<n_dofs; i++)
742 std::vector<Real> psi(
specpoly(
dim, element_order, q_point[sp], matsize));
743 for (
unsigned int i=0; i<matsize; i++)
745 temperr[2] += psi[i]*Pu_z_h(i);
748 temperr[2] -= grad_u_h(2);
750 #endif // LIBMESH_DIM > 2 754 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 758 for (
unsigned int i=0; i<n_dofs; i++)
763 std::vector<Real> psi(
specpoly(
dim, element_order, q_point[sp], matsize));
764 for (
unsigned int i=0; i<matsize; i++)
766 temperr[0] += psi[i]*Pu_x_h(i);
768 temperr[1] += psi[i]*Pu_y_h(i);
769 temperr[3] += psi[i]*Pu_xy_h(i);
772 temperr[2] += psi[i]*Pu_z_h(i);
773 temperr[4] += psi[i]*Pu_xz_h(i);
774 temperr[5] += psi[i]*Pu_yz_h(i);
778 temperr[0] -= hess_u_h(0,0);
780 temperr[1] -= hess_u_h(1,1);
781 temperr[3] -= hess_u_h(0,1);
784 temperr[2] -= hess_u_h(2,2);
785 temperr[4] -= hess_u_h(0,2);
786 temperr[5] -= hess_u_h(1,2);
789 libmesh_error_msg(
"ERROR: --enable-second-derivatives is required \nfor _sobolev_order == 2!");
800 if (norm_type ==
L_INF)
801 element_error = std::max(element_error, std::abs(
weight*temperr[0]));
803 for (
unsigned int i=0; i != LIBMESH_DIM; ++i)
804 element_error = std::max(element_error, std::abs(
weight*temperr[i]));
806 for (
unsigned int i=0; i != 6; ++i)
807 element_error = std::max(element_error, std::abs(
weight*temperr[i]));
808 else if (norm_type ==
L2)
811 for (
unsigned int i=0; i != LIBMESH_DIM; ++i)
821 for (
unsigned int i=0; i != LIBMESH_DIM; ++i)
824 for (
unsigned int i=3; i != 6; ++i)
830 if (norm_type ==
L_INF ||
834 else if (norm_type ==
L2 ||
853 Patch::const_iterator patch_re_it;
854 Patch::const_iterator patch_re_end;
857 Patch current_elem_patch(
mesh.processor_id());
862 patch_re_it = patch.begin();
863 patch_re_end = patch.end();
869 current_elem_patch.build_around_element (elem, 0,
873 patch_re_it = current_elem_patch.begin();
874 patch_re_end = current_elem_patch.end();
878 for (
unsigned int i = 0 ; patch_re_it != patch_re_end; ++patch_re_it, ++i)
881 const Elem * e_p = *patch_re_it;
897 (std::sqrt(new_error_per_cell[i]));
907 (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...
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...
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.