24 #include "libmesh/hp_coarsentest.h" 25 #include "libmesh/dense_matrix.h" 26 #include "libmesh/dense_vector.h" 27 #include "libmesh/dof_map.h" 28 #include "libmesh/fe_base.h" 29 #include "libmesh/fe_interface.h" 30 #include "libmesh/libmesh_logging.h" 31 #include "libmesh/elem.h" 32 #include "libmesh/error_vector.h" 33 #include "libmesh/mesh_base.h" 34 #include "libmesh/mesh_refinement.h" 35 #include "libmesh/quadrature.h" 36 #include "libmesh/system.h" 37 #include "libmesh/tensor_value.h" 38 #include "libmesh/smoothness_estimator.h" 40 #ifdef LIBMESH_ENABLE_AMR 70 const unsigned int n_dofs =
78 const unsigned int n_coarse_dofs =
83 Ke.
resize(n_coarse_dofs, n_coarse_dofs);
100 for (
unsigned int i=0; i != n_dofs; i++)
103 val += (*phi)[i][qp] *
114 Fe(i) += (*JxW)[qp] *
115 (*phi_coarse)[i][qp]*val;
117 Fe(i) += (*JxW)[qp] *
118 (grad*(*dphi_coarse)[i][qp]);
120 Fe(i) += (*JxW)[qp] *
125 Ke(i,j) += (*JxW)[qp] *
126 (*phi_coarse)[i][qp]*(*phi_coarse)[j][qp];
128 Ke(i,j) += (*JxW)[qp] *
129 (*dphi_coarse)[i][qp]*(*dphi_coarse)[j][qp];
131 Ke(i,j) += (*JxW)[qp] *
132 ((*d2phi_coarse)[i][qp].contract((*
d2phi_coarse)[j][qp]));
140 LOG_SCOPE(
"select_refinement()",
"HPCoarsenTest");
146 const unsigned int dim =
mesh.mesh_dimension();
155 const unsigned int sys_num = system.
number();
161 "ERROR: component_scale is the wrong size:\n" 162 <<
" component_scale.size()=" 181 std::vector<ErrorVectorReal> h_error_per_cell(
mesh.max_elem_id(), 0.);
182 std::vector<ErrorVectorReal> p_error_per_cell(
mesh.max_elem_id(), 0.);
185 for (
unsigned int var=0; var<
n_vars; var++)
201 unsigned int cached_coarse_p_level = 0;
212 fe->attach_quadrature_rule (
qrule.get());
216 JxW = &(
fe->get_JxW());
220 phi = &(
fe->get_phi());
233 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 237 libmesh_error_msg(
"Minimization of H2 error without second derivatives is not possible.");
243 for (
auto & elem :
mesh.active_local_element_ptr_range())
254 if (elem->parent() &&
256 cached_coarse_p_level != elem->
p_level()))
261 cached_coarse_p_level = elem->
p_level();
280 const unsigned int n_qp =
qrule->n_points();
283 const unsigned int n_dofs =
287 const unsigned int n_nodes = elem->n_nodes();
297 if (elem->p_level() == 0)
299 unsigned int n_vertices = 0;
300 for (
unsigned int n = 0; n !=
n_nodes; ++n)
301 if (elem->is_vertex(n))
304 const Node & node = elem->node_ref(n);
308 average_val /= n_vertices;
312 unsigned int old_elem_level = elem->p_level();
313 elem->hack_p_level(old_elem_level - 1);
317 const unsigned int n_coarse_dofs =
320 elem->hack_p_level(old_elem_level);
322 Ke.
resize(n_coarse_dofs, n_coarse_dofs);
335 for (
unsigned int i=0; i != n_dofs; i++)
338 val += (*phi)[i][qp] *
349 Fe(i) += (*JxW)[qp] *
350 (*phi_coarse)[i][qp]*val;
352 Fe(i) += (*JxW)[qp] *
353 grad * (*dphi_coarse)[i][qp];
355 Fe(i) += (*JxW)[qp] *
360 Ke(i,j) += (*JxW)[qp] *
361 (*phi_coarse)[i][qp]*(*phi_coarse)[j][qp];
363 Ke(i,j) += (*JxW)[qp] *
364 (*dphi_coarse)[i][qp]*(*dphi_coarse)[j][qp];
366 Ke(i,j) += (*JxW)[qp] *
367 ((*d2phi_coarse)[i][qp].contract((*
d2phi_coarse)[j][qp]));
377 for (
unsigned int qp=0; qp<n_qp; qp++)
382 for (
unsigned int i=0; i<n_dofs; i++)
385 value_error += (*phi)[i][qp] *
392 if (elem->p_level() == 0)
394 value_error -= average_val;
400 value_error -= (*phi_coarse)[i][qp] *
Up(i);
414 (*JxW)[qp] * grad_error.
norm_sq());
418 (*JxW)[qp] * hessian_error.
norm_sq());
427 h_error_per_cell[e_id] =
428 std::numeric_limits<ErrorVectorReal>::max() / 2;
443 unsigned int n_coarse_dofs =
447 for (
unsigned int qp=0; qp<n_qp; qp++)
454 for (
unsigned int i=0; i != n_dofs; ++i)
457 value_error += (*phi)[i][qp] *
465 for (
unsigned int i=0; i != n_coarse_dofs; ++i)
467 value_error -= (*phi_coarse)[i][qp] *
Uc(i);
480 (*JxW)[qp] * grad_error.
norm_sq());
484 (*JxW)[qp] * hessian_error.
norm_sq());
496 for (
auto & elem :
mesh.active_local_element_ptr_range())
505 unsigned int dofs_per_elem = 0, dofs_per_p_elem = 0;
508 for (
unsigned int var=0; var<
n_vars; var++)
525 const unsigned int new_h_dofs = dofs_per_elem *
526 (elem->n_children() - 1);
528 const unsigned int new_p_dofs = dofs_per_p_elem -
540 std::sqrt(p_error_per_cell[e_id]) *
p_weight / new_p_dofs;
542 std::sqrt(h_error_per_cell[e_id]) /
543 static_cast<Real>(new_h_dofs);
544 if (smoothness[elem->id()] > smoothness.
mean() && p_value > h_value)
561 #endif // #ifdef LIBMESH_ENABLE_AMR class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
dof_id_type dof_number(const unsigned int s, const unsigned int var, const unsigned int comp) const
const std::vector< std::vector< Real > > * phi
The shape functions and their derivatives.
const Elem * parent() const
void add_scaled(const TypeVector< T2 > &, const T &)
Add a scaled value to this vector without creating a temporary.
std::vector< float > component_scale
This vector can be used to "scale" certain variables in a system.
void add_projection(const System &, const Elem *, unsigned int var)
The helper function which adds individual fine element data to the coarse element projection...
A Node is like a Point, but with more information.
static unsigned int n_dofs(const unsigned int dim, const FEType &fe_t, const ElemType t)
virtual void zero() override final
Set every element in the vector to 0.
void subtract_scaled(const TypeVector< T2 > &, const T &)
Subtract a scaled value from this vector without creating a temporary.
virtual void select_refinement(System &system) override
This pure virtual function must be redefined in derived classes to take a mesh flagged for h refineme...
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
virtual void zero() override final
Sets all elements of the matrix to 0 and resets any decomposition flag which may have been previously...
static Point inverse_map(const unsigned int dim, const Elem *elem, const Point &p, const Real tolerance=TOLERANCE, const bool secure=true, const bool extra_checks=true)
The ErrorVector is a specialization of the StatisticsVector for error data computed on a finite eleme...
std::unique_ptr< FEBase > fe_coarse
const std::vector< std::vector< Real > > * phi_coarse
int _extra_order
Extra order to use for quadrature rule.
void resize(const unsigned int n)
Resize the vector.
const FEType & variable_type(const unsigned int c) const
This is the base class from which all geometric element types are derived.
auto norm_sq() const -> decltype(std::norm(T()))
const std::vector< std::vector< RealGradient > > * dphi_coarse
unsigned int p_level() const
This class defines a vector in LIBMESH_DIM dimensional Real or Complex space.
The libMesh namespace provides an interface to certain functionality in the library.
const std::vector< Point > * xyz_values
Quadrature locations.
const MeshBase & get_mesh() const
DIE A HORRIBLE DEATH HERE typedef float ErrorVectorReal
DenseVector< Number > Uc
Coefficients for projected coarse and projected p-derefined solutions.
This is the MeshBase class.
SimpleRange< ChildRefIter > child_ref_range()
Returns a range with all children of a parent element, usable in range-based for loops.
Number current_solution(const dof_id_type global_dof_number) const
const std::vector< Real > * JxW
Mapping jacobians.
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
const dof_id_type n_nodes
unsigned int number() const
auto norm_sq() const -> decltype(std::norm(T()))
void subtract_scaled(const TypeTensor< T2 > &, const T &)
Subtract a scaled value from this tensor without creating a temporary.
Implements (adaptive) mesh refinement algorithms for a MeshBase.
This class implements the Smoothness estimate.
Manages consistently variables, degrees of freedom, and coefficient vectors.
static std::unique_ptr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.
bool make_flags_parallel_consistent()
Copy refinement flags on ghost elements from their local processors.
virtual void estimate_smoothness(const System &system, ErrorVector &smoothness_per_cell, const NumericVector< Number > *solution_vector=nullptr)
This function uses the Legendre expansion of solution to estimate coefficient decay to quantify the s...
void add_scaled(const TypeTensor< T2 > &, const T &)
Add a scaled tensor to this tensor without creating a temporary.
CompareTypes< T, T2 >::supertype contract(const TypeTensor< T2 > &) const
Multiply 2 tensors together to return a scalar, i.e.
std::unique_ptr< FEBase > fe
The finite element objects for fine and coarse elements.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
FEContinuity
defines an enum for finite element types to libmesh_assert a certain level (or type? Hcurl?) of continuity.
void resize(const unsigned int new_m, const unsigned int new_n)
Resizes the matrix to the specified size and calls zero().
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
unsigned int mesh_dimension() const
virtual unsigned int size() const override final
const std::vector< std::vector< RealTensor > > * d2phi_coarse
const std::vector< std::vector< RealTensor > > * d2phi
void cholesky_solve(const DenseVector< T2 > &b, DenseVector< T2 > &x)
For symmetric positive definite (SPD) matrices.
const std::vector< std::vector< RealGradient > > * dphi
unsigned int n_vars() const
std::vector< dof_id_type > dof_indices
Global DOF indices for fine elements.
Real p_weight
Because the coarsening test seems to always choose p refinement, we're providing an option to make h ...
std::vector< Point > coarse_qpoints
void hack_p_level(const unsigned int p)
Sets the value of the p-refinement level for the element without altering the p-level of its ancestor...
const DofMap & get_dof_map() const
Elem * coarse
The coarse element on which a solution projection is cached.
std::unique_ptr< QBase > qrule
The quadrature rule for the fine element.
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
This class defines a tensor in LIBMESH_DIM dimensional Real or Complex space.
virtual Real mean() const override
DenseMatrix< Number > Ke
Linear system for projections.