20 #include "libmesh/dense_matrix.h" 21 #include "libmesh/dense_vector.h" 22 #include "libmesh/dof_map.h" 23 #include "libmesh/elem.h" 24 #include "libmesh/error_vector.h" 25 #include "libmesh/fe_base.h" 26 #include "libmesh/fem_context.h" 27 #include "libmesh/libmesh_logging.h" 28 #include "libmesh/mesh_base.h" 29 #include "libmesh/numeric_vector.h" 30 #include "libmesh/quadrature_grid.h" 31 #include "libmesh/system.h" 32 #include "libmesh/tensor_value.h" 33 #include "libmesh/threads.h" 34 #include "libmesh/weighted_patch_recovery_error_estimator.h" 35 #include "libmesh/enum_error_estimator_type.h" 36 #include "libmesh/enum_order.h" 37 #include "libmesh/enum_norm_type.h" 38 #include "libmesh/enum_to_string.h" 60 LOG_SCOPE(
"estimate_error()",
"WeightedPatchRecoveryErrorEstimator");
67 error_per_cell.resize (
mesh.max_elem_id());
68 std::fill (error_per_cell.begin(), error_per_cell.end(), 0.);
72 if (solution_vector && solution_vector != system.
solution.get())
77 newsol->
swap(*sys.solution);
85 mesh.active_local_elements_end(),
100 if (solution_vector && solution_vector != system.
solution.get())
105 newsol->
swap(*sys.solution);
118 const unsigned int dim =
mesh.mesh_dimension();
128 for (
const auto & elem : range)
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 ||
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);
284 if (norm_type ==
L2 ||
333 for (
const auto & e_p : patch)
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 ||
550 Patch::const_iterator patch_re_it;
551 Patch::const_iterator patch_re_end;
559 patch_re_it = patch.begin();
560 patch_re_end = patch.end();
570 patch_re_it = patch_re.begin();
571 patch_re_end = patch_re.end();
588 for (
unsigned int e = 0 ; patch_re_it != patch_re_end; ++patch_re_it, ++e)
593 const Elem * e_p = *patch_re_it;
618 const unsigned int n_dofs =
619 cast_int<unsigned int>(dof_indices.size());
622 Real element_error = 0;
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();
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]);
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
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.
virtual void pre_fe_reinit(const System &, const Elem *e)
Reinitializes local data vectors/matrices on the current geometric element.
void operator()(const ConstElemRange &range) const
virtual void estimate_error(const System &system, ErrorVector &error_per_cell, const NumericVector< Number > *solution_vector=nullptr, bool estimate_parent_error=false) override
This function uses the Patch Recovery error estimate to estimate the error on each cell...
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...
This class implements useful utility functions for a patch of elements.
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
void parallel_for(const Range &range, const Body &body)
Execute the provided function object in parallel on the specified range.
The ErrorVector is a specialization of the StatisticsVector for error data computed on a finite eleme...
ErrorEstimatorType
Defines an enum for the different types of error estimators which are available.
void resize(const unsigned int n)
Resize the vector.
const FEType & variable_type(const unsigned int c) const
friend class EstimateError
const System & system
Function to set the boolean patch_reuse in case the user wants to change the default behaviour of pat...
This is the base class from which all geometric element types are derived.
const Parallel::Communicator & comm() const
unsigned int p_level() const
OrderWrapper order
The approximation order of the element.
This class defines a vector in LIBMESH_DIM dimensional Real or Complex space.
The StoredRange class defines a contiguous, divisible set of objects.
virtual ErrorEstimatorType type() const override
The libMesh namespace provides an interface to certain functionality in the library.
const MeshBase & get_mesh() const
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)
Erases any elements in the current patch, then builds a new patch containing element elem by repeated...
This is the MeshBase class.
StoredRange< MeshBase::const_element_iterator, const Elem * > ConstElemRange
Number current_solution(const dof_id_type global_dof_number) const
This class handles the numbering of degrees of freedom on a mesh.
std::unique_ptr< QBase > default_quadrature_rule(const unsigned int dim, const int extraorder=0) const
Manages consistently variables, degrees of freedom, and coefficient vectors.
FEMNormType type(unsigned int var) const
std::unique_ptr< NumericVector< Number > > solution
Data structure to hold solution values.
static std::unique_ptr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.
This class provides all data required for a physics package (e.g.
void reduce_error(std::vector< ErrorVectorReal > &error_per_cell, const Parallel::Communicator &comm) const
This method takes the local error contributions in error_per_cell from each processor and combines th...
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.
This class creates quadrature points on a uniform grid, with order+1 points on an edge...
std::string enum_to_string(const T e)
Real weight_sq(unsigned int var) const
Real weight(unsigned int var) const
void lu_solve(const DenseVector< T > &b, DenseVector< T > &x)
Solve the system Ax=b given the input vector b.
ErrorVector & error_per_cell
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual void swap(NumericVector< T > &v)
Swaps the contents of this with v.
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
This class defines a tensor in LIBMESH_DIM dimensional Real or Complex space.
spin_mutex spin_mtx
A convenient spin mutex object which can be used for obtaining locks.