49 LOG_SCOPE(
"estimate_error()",
"WeightedPatchRecoveryErrorEstimator");
57 std::fill (error_per_cell.begin(), error_per_cell.end(), 0.);
61 if (solution_vector && solution_vector != system.
solution.get())
66 newsol->
swap(*sys.solution);
89 if (solution_vector && solution_vector != system.
solution.get())
94 newsol->
swap(*sys.solution);
110 const unsigned int n_vars = system.n_vars();
113 const DofMap & dof_map = system.get_dof_map();
120 const Elem * elem = *elem_it;
132 if (this->error_estimator.patch_reuse && error_per_cell[e_id] != 0)
139 error_estimator.patch_growth_strategy);
145 std::vector<Real> new_error_per_cell(1, 0.);
146 if (this->error_estimator.patch_reuse)
147 new_error_per_cell.resize(patch.size(), 0.);
151 for (
unsigned int var=0; var<
n_vars; var++)
153 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 155 bool is_valid_norm_type =
156 error_estimator.error_norm.type(var) ==
L2 ||
157 error_estimator.error_norm.type(var) ==
H1_SEMINORM ||
158 error_estimator.error_norm.type(var) ==
H2_SEMINORM ||
162 error_estimator.error_norm.type(var) ==
L_INF ||
165 libmesh_assert (is_valid_norm_type);
168 libmesh_assert (error_estimator.error_norm.type(var) ==
L2 ||
169 error_estimator.error_norm.type(var) ==
L_INF ||
170 error_estimator.error_norm.type(var) ==
H1_SEMINORM ||
181 bool is_valid_norm_type =
182 ((error_estimator.error_norm.type(var) ==
L2 ||
183 error_estimator.error_norm.type(var) ==
H1_SEMINORM ||
187 error_estimator.error_norm.type(var) ==
H2_SEMINORM) &&
188 (error_estimator.error_norm.type(var-1) ==
L2 ||
189 error_estimator.error_norm.type(var-1) ==
H1_SEMINORM ||
193 error_estimator.error_norm.type(var-1) ==
H2_SEMINORM)) ||
194 ((error_estimator.error_norm.type(var) ==
L_INF ||
197 (error_estimator.error_norm.type(var-1) ==
L_INF ||
200 libmesh_assert (is_valid_norm_type);
205 if (error_estimator.error_norm.weight(var) == 0.0)
continue;
210 const Order element_order =
static_cast<Order> 220 fe->attach_quadrature_rule (qrule.get());
223 const std::vector<Real> & JxW = fe->get_JxW();
224 const std::vector<Point> & q_point = fe->get_xyz();
235 if (error_estimator.error_norm.type(var) ==
L2 ||
236 error_estimator.error_norm.type(var) ==
L_INF)
238 phi = &(fe->get_phi());
241 if (error_estimator.error_norm.type(var) ==
H1_SEMINORM ||
246 dphi = &(fe->get_dphi());
248 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 250 if (error_estimator.error_norm.type(var) ==
H2_SEMINORM ||
252 d2phi = &(fe->get_d2phi());
256 std::vector<dof_id_type> dof_indices;
260 unsigned int matsize = element_order + 1;
263 matsize *= (element_order + 2);
268 matsize *= (element_order + 3);
275 if (error_estimator.error_norm.type(var) ==
L2 ||
276 error_estimator.error_norm.type(var) ==
L_INF)
280 else if (error_estimator.error_norm.type(var) ==
H1_SEMINORM ||
282 error_estimator.error_norm.type(var) ==
H2_SEMINORM ||
293 else if (error_estimator.error_norm.type(var) ==
H1_X_SEMINORM)
297 else if (error_estimator.error_norm.type(var) ==
H1_Y_SEMINORM)
299 libmesh_assert(LIBMESH_DIM > 1);
302 else if (error_estimator.error_norm.type(var) ==
H1_Z_SEMINORM)
304 libmesh_assert(LIBMESH_DIM > 2);
309 if (error_estimator.error_norm.type(var) ==
H2_SEMINORM ||
324 Patch::const_iterator patch_it = patch.begin();
325 const Patch::const_iterator patch_end = patch.end();
327 for (; patch_it != patch_end; ++patch_it)
330 const Elem * e_p = *patch_it;
338 libmesh_assert (dof_indices.size() == phi->size());
340 const unsigned int n_dofs =
341 cast_int<unsigned int>(dof_indices.size());
342 const unsigned int n_qp = qrule->n_points();
346 for (
unsigned int qp=0; qp<n_qp; qp++)
349 std::vector<Real> psi(
specpoly(dim, element_order, q_point[qp], matsize));
352 for (
unsigned int i=0; i<Kp.
m(); i++)
353 for (
unsigned int j=0; j<Kp.
n(); j++)
354 Kp(i,j) += JxW[qp]*psi[i]*psi[j];
356 if (error_estimator.error_norm.type(var) ==
L2 ||
357 error_estimator.error_norm.type(var) ==
L_INF)
363 for (
unsigned int i=0; i<n_dofs; i++)
364 u_h += (*phi)[i][qp]*system.current_solution (dof_indices[i]);
367 for (std::size_t i=0; i<psi.size(); i++)
368 F(i) = JxW[qp]*u_h*psi[i];
371 else if (error_estimator.error_norm.type(var) ==
H1_SEMINORM ||
378 for (std::size_t i=0; i<n_dofs; i++)
380 system.current_solution(dof_indices[i]));
385 for (std::size_t i=0; i<psi.size(); i++)
387 Fx(i) += JxW[qp]*grad_u_h(0)*psi[i];
389 Fy(i) += JxW[qp]*grad_u_h(1)*psi[i];
392 Fz(i) += JxW[qp]*grad_u_h(2)*psi[i];
396 else if (error_estimator.error_norm.type(var) ==
H1_X_SEMINORM)
402 for (
unsigned int i=0; i<n_dofs; i++)
404 system.current_solution(dof_indices[i]));
409 for (std::size_t i=0; i<psi.size(); i++)
411 Fx(i) += JxW[qp]*grad_u_h(0)*psi[i];
414 else if (error_estimator.error_norm.type(var) ==
H1_Y_SEMINORM)
420 for (
unsigned int i=0; i<n_dofs; i++)
422 system.current_solution(dof_indices[i]));
427 for (std::size_t i=0; i<psi.size(); i++)
429 Fy(i) += JxW[qp]*grad_u_h(1)*psi[i];
432 else if (error_estimator.error_norm.type(var) ==
H1_Z_SEMINORM)
438 for (
unsigned int i=0; i<n_dofs; i++)
440 system.current_solution(dof_indices[i]));
445 for (std::size_t i=0; i<psi.size(); i++)
447 Fz(i) += JxW[qp]*grad_u_h(2)*psi[i];
450 else if (error_estimator.error_norm.type(var) ==
H2_SEMINORM ||
453 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 458 for (
unsigned int i=0; i<n_dofs; i++)
460 system.current_solution(dof_indices[i]));
465 for (std::size_t i=0; i<psi.size(); i++)
467 Fx(i) += JxW[qp]*hess_u_h(0,0)*psi[i];
469 Fy(i) += JxW[qp]*hess_u_h(1,1)*psi[i];
470 Fxy(i) += JxW[qp]*hess_u_h(0,1)*psi[i];
473 Fz(i) += JxW[qp]*hess_u_h(2,2)*psi[i];
474 Fxz(i) += JxW[qp]*hess_u_h(0,2)*psi[i];
475 Fyz(i) += JxW[qp]*hess_u_h(1,2)*psi[i];
479 libmesh_error_msg(
"ERROR: --enable-second-derivatives is required \nfor _sobolev_order == 2!");
483 libmesh_error_msg(
"Unsupported error norm type!");
493 if (error_estimator.error_norm.type(var) ==
L2 ||
494 error_estimator.error_norm.type(var) ==
L_INF)
498 else if (error_estimator.error_norm.type(var) ==
H1_SEMINORM ||
500 error_estimator.error_norm.type(var) ==
H2_SEMINORM ||
511 else if (error_estimator.error_norm.type(var) ==
H1_X_SEMINORM)
515 else if (error_estimator.error_norm.type(var) ==
H1_Y_SEMINORM)
519 else if (error_estimator.error_norm.type(var) ==
H1_Z_SEMINORM)
525 if (error_estimator.error_norm.type(var) ==
H2_SEMINORM ||
540 Patch::const_iterator patch_re_it;
541 Patch::const_iterator patch_re_end;
546 if (this->error_estimator.patch_reuse)
549 patch_re_it = patch.begin();
550 patch_re_end = patch.end();
557 error_estimator.patch_growth_strategy);
560 patch_re_it = patch_re.begin();
561 patch_re_end = patch_re.end();
573 error_estimator.weight_functions[var]->init_context(femcontext);
576 for (
unsigned int e = 0 ; patch_re_it != patch_re_end; ++patch_re_it, ++e)
581 const Elem * e_p = *patch_re_it;
588 Elem * e_p_cast =
const_cast<Elem *
>(*patch_re_it);
598 if (this->error_estimator.patch_reuse && error_per_cell[e_p_id] != 0.)
607 libmesh_assert (dof_indices.size() == phi->size());
610 const unsigned int n_dofs =
611 cast_int<unsigned int>(dof_indices.size());
614 Real element_error = 0;
620 QGrid samprule (dim, qorder);
624 fe->attach_quadrature_rule (&samprule);
627 const unsigned int n_sp =
628 cast_int<unsigned int>(JxW.size());
631 for (
unsigned int sp=0; sp<n_sp; sp++)
635 std::vector<Number> temperr(6,0.0);
637 if (error_estimator.error_norm.type(var) ==
L2 ||
638 error_estimator.error_norm.type(var) ==
L_INF)
643 for (
unsigned int i=0; i<n_dofs; i++)
644 u_h += (*phi)[i][sp]*system.current_solution (dof_indices[i]);
647 std::vector<Real> psi(
specpoly(dim, element_order, q_point[sp], matsize));
648 for (
unsigned int i=0; i<matsize; i++)
650 temperr[0] += psi[i]*Pu_h(i);
655 else if (error_estimator.error_norm.type(var) ==
H1_SEMINORM ||
661 for (
unsigned int i=0; i<n_dofs; i++)
663 system.current_solution(dof_indices[i]));
666 std::vector<Real> psi(
specpoly(dim, element_order, q_point[sp], matsize));
668 for (
unsigned int i=0; i<matsize; i++)
670 temperr[0] += psi[i]*Pu_x_h(i);
672 temperr[1] += psi[i]*Pu_y_h(i);
675 temperr[2] += psi[i]*Pu_z_h(i);
678 temperr[0] -= grad_u_h(0);
680 temperr[1] -= grad_u_h(1);
683 temperr[2] -= grad_u_h(2);
686 else if (error_estimator.error_norm.type(var) ==
H1_X_SEMINORM)
691 for (
unsigned int i=0; i<n_dofs; i++)
693 system.current_solution(dof_indices[i]));
696 std::vector<Real> psi(
specpoly(dim, element_order, q_point[sp], matsize));
697 for (
unsigned int i=0; i<matsize; i++)
699 temperr[0] += psi[i]*Pu_x_h(i);
702 temperr[0] -= grad_u_h(0);
704 else if (error_estimator.error_norm.type(var) ==
H1_Y_SEMINORM)
709 for (
unsigned int i=0; i<n_dofs; i++)
711 system.current_solution(dof_indices[i]));
714 std::vector<Real> psi(
specpoly(dim, element_order, q_point[sp], matsize));
715 for (
unsigned int i=0; i<matsize; i++)
717 temperr[1] += psi[i]*Pu_y_h(i);
720 temperr[1] -= grad_u_h(1);
722 else if (error_estimator.error_norm.type(var) ==
H1_Z_SEMINORM)
727 for (
unsigned int i=0; i<n_dofs; i++)
729 system.current_solution(dof_indices[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[2] += psi[i]*Pu_z_h(i);
738 temperr[2] -= grad_u_h(2);
740 else if (error_estimator.error_norm.type(var) ==
H2_SEMINORM ||
743 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 747 for (
unsigned int i=0; i<n_dofs; i++)
749 system.current_solution(dof_indices[i]));
752 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);
758 temperr[3] += psi[i]*Pu_xy_h(i);
761 temperr[2] += psi[i]*Pu_z_h(i);
762 temperr[4] += psi[i]*Pu_xz_h(i);
763 temperr[5] += psi[i]*Pu_yz_h(i);
767 temperr[0] -= hess_u_h(0,0);
769 temperr[1] -= hess_u_h(1,1);
770 temperr[3] -= hess_u_h(0,1);
773 temperr[2] -= hess_u_h(2,2);
774 temperr[4] -= hess_u_h(0,2);
775 temperr[5] -= hess_u_h(1,2);
778 libmesh_error_msg(
"ERROR: --enable-second-derivatives is required \nfor _sobolev_order == 2!");
783 Number weight = (*error_estimator.weight_functions[var])(femcontext, q_point[sp], system.
time);
789 if (error_estimator.error_norm.type(var) ==
L_INF)
792 for (
unsigned int i=0; i != LIBMESH_DIM; ++i)
795 for (
unsigned int i=0; i != 6; ++i)
797 else if (error_estimator.error_norm.type(var) ==
L2)
799 else if (error_estimator.error_norm.type(var) ==
H1_SEMINORM)
800 for (
unsigned int i=0; i != LIBMESH_DIM; ++i)
802 else if (error_estimator.error_norm.type(var) ==
H1_X_SEMINORM)
804 else if (error_estimator.error_norm.type(var) ==
H1_Y_SEMINORM)
806 else if (error_estimator.error_norm.type(var) ==
H1_Z_SEMINORM)
808 else if (error_estimator.error_norm.type(var) ==
H2_SEMINORM)
810 for (
unsigned int i=0; i != LIBMESH_DIM; ++i)
813 for (
unsigned int i=3; i != 6; ++i)
819 if (error_estimator.error_norm.type(var) ==
L_INF ||
822 new_error_per_cell[e] += error_estimator.error_norm.weight(var) * element_error;
823 else if (error_estimator.error_norm.type(var) ==
L2 ||
824 error_estimator.error_norm.type(var) ==
H1_SEMINORM ||
828 error_estimator.error_norm.type(var) ==
H2_SEMINORM)
829 new_error_per_cell[e] += error_estimator.error_norm.weight_sq(var) * element_error;
831 libmesh_error_msg(
"Unsupported error norm type!");
842 Patch::const_iterator patch_re_it;
843 Patch::const_iterator patch_re_end;
848 if (this->error_estimator.patch_reuse)
851 patch_re_it = patch.begin();
852 patch_re_end = patch.end();
859 error_estimator.patch_growth_strategy);
862 patch_re_it = current_elem_patch.begin();
863 patch_re_end = current_elem_patch.end();
867 for (
unsigned int i = 0 ; patch_re_it != patch_re_end; ++patch_re_it, ++i)
870 const Elem * e_p = *patch_re_it;
876 if (error_estimator.error_norm.type(0) ==
L2 ||
877 error_estimator.error_norm.type(0) ==
H1_SEMINORM ||
884 if (!error_per_cell[e_p_id])
886 (std::sqrt(new_error_per_cell[i]));
890 libmesh_assert (error_estimator.error_norm.type(0) ==
L_INF ||
894 if (!error_per_cell[e_p_id])
896 (new_error_per_cell[i]);
Manages the family, order, etc. parameters for a given FE.
virtual void pre_fe_reinit(const System &, const Elem *e)
unsigned int p_level() const
const FEType & variable_type(const unsigned int c) const
void parallel_for(const Range &range, const Body &body)
void add_scaled(const TypeTensor< T2 > &, const T)
void add_scaled(const TypeVector< T2 > &, const T)
void resize(const unsigned int n)
friend class EstimateError
The base class for all geometric element types.
const class libmesh_nullptr_t libmesh_nullptr
Utility class for defining generic ranges for threading.
const unsigned int n_vars
const_iterator begin() const
long double max(long double a, double b)
DIE A HORRIBLE DEATH HERE typedef float ErrorVectorReal
void build_around_element(const Elem *elem, const unsigned int target_patch_size=10, PMF patchtype=&Patch::add_local_face_neighbors)
StoredRange< MeshBase::const_element_iterator, const Elem * > ConstElemRange
Manages the degrees of freedom (DOFs) in a simulation.
virtual dof_id_type max_elem_id() const =0
const MeshBase & get_mesh() const
virtual element_iterator active_local_elements_begin()=0
std::unique_ptr< NumericVector< Number > > solution
static std::unique_ptr< FEGenericBase > build(const unsigned int dim, const FEType &type)
const_iterator end() const
static std::vector< Real > specpoly(const unsigned int dim, const Order order, const Point p, const unsigned int matsize)
void reduce_error(std::vector< ErrorVectorReal > &error_per_cell, const Parallel::Communicator &comm LIBMESH_CAN_DEFAULT_TO_COMMWORLD) const
Implements grid-based quadrature rules suitable for non-smooth functions.
std::vector< object_type >::const_iterator const_iterator
std::unique_ptr< QBase > default_quadrature_rule(const unsigned int dim, const int extraorder=0) const
void lu_solve(const DenseVector< T > &b, DenseVector< T > &x)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual void swap(NumericVector< T > &v)
const Parallel::Communicator & comm() const
unsigned int mesh_dimension() const
virtual void estimate_error(const System &system, ErrorVector &error_per_cell, const NumericVector< Number > *solution_vector=libmesh_nullptr, bool estimate_parent_error=false) libmesh_override
void operator()(const ConstElemRange &range) const
virtual element_iterator active_local_elements_end()=0
processor_id_type processor_id() const
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const