34 face_midpoint(
const unsigned int face_no)
37 const auto face_reference_cell =
41 for (
const auto face_vertex_no : face_reference_cell.
vertex_indices())
43 const auto vertex_no = reference_cell.face_to_cell_vertices(
48 midpoint += reference_cell.template vertex<dim>(vertex_no);
51 midpoint /= reference_cell.face_reference_cell(0).n_vertices();
59 std::vector<unsigned int>
60 get_dpo_vector_fe_p(
const unsigned int dim,
const unsigned int degree)
113 std::vector<Point<dim>>
114 unit_support_points_fe_p(
const unsigned int degree)
118 std::vector<Point<dim>> unit_points;
119 const auto reference_cell = ReferenceCells::get_simplex<dim>();
125 unit_points.emplace_back(reference_cell.template barycenter<dim>());
130 const auto dpo = get_dpo_vector_fe_p(dim, degree);
134 for (
const unsigned int d : reference_cell.vertex_indices())
135 unit_points.push_back(reference_cell.template vertex<dim>(d));
137 for (
const unsigned int l : reference_cell.line_indices())
140 unit_points[reference_cell.line_to_cell_vertices(l, 0)];
142 unit_points[reference_cell.line_to_cell_vertices(l, 1)];
143 for (
unsigned int p = 0; p < dpo[1]; ++p)
144 unit_points.push_back((
double(dpo[1] - p) / (dpo[1] + 1)) * p0 +
145 (double(p + 1) / (dpo[1] + 1)) * p1);
148 if (dim > 1 && dpo[2] > 0)
152 unit_points.push_back(reference_cell.template barycenter<dim>());
154 for (
const unsigned int f : reference_cell.face_indices())
155 unit_points.push_back(face_midpoint<dim>(f));
162 std::vector<Point<0>>
163 unit_support_points_fe_p(
const unsigned int )
173 std::vector<std::vector<
Point<dim - 1>>>
174 unit_face_support_points_fe_p(
175 const unsigned int degree,
186 std::vector<std::vector<
Point<dim - 1>>> unit_face_points;
192 unit_face_points.emplace_back(
193 unit_support_points_fe_p<dim - 1>(degree));
196 return unit_face_points;
206 constraints_fe_p(
const unsigned int )
215 constraints_fe_p<2>(
const unsigned int degree)
217 constexpr int dim = 2;
226 std::vector<
Point<dim - 1>> constraint_points;
228 constraint_points.emplace_back(0.5);
232 constraint_points.emplace_back(0.25);
234 constraint_points.emplace_back(0.75);
240 const unsigned int n_dofs_constrained = constraint_points.size();
241 unsigned int n_dofs_per_face = degree + 1;
247 for (
unsigned int i = 0; i < n_dofs_constrained; ++i)
248 for (
unsigned int j = 0; j < n_dofs_per_face; ++j)
250 interface_constraints(i, j) =
251 poly.compute_value(j, constraint_points[i]);
257 if (std::fabs(interface_constraints(i, j)) < 1e-13)
258 interface_constraints(i, j) = 0;
260 return interface_constraints;
269 std::vector<unsigned int>
270 get_dpo_vector_fe_dgp(
const unsigned int dim,
const unsigned int degree)
275 std::vector<unsigned int> dpo(dim + 1, 0U);
286 const auto continuous_dpo = get_dpo_vector_fe_p(dim, degree);
293 continuous_dpo[dim]};
301 continuous_dpo[dim]};
311 continuous_dpo[dim]};
322template <
int dim,
int spacedim>
326 const std::vector<
Point<dim>> &unit_support_points,
327 const std::vector<std::vector<
Point<dim - 1>>> unit_face_support_points,
332 std::vector<
bool>(fe_data.dofs_per_cell),
343template <
int dim,
int spacedim>
344std::pair<Table<2, bool>, std::vector<unsigned int>>
348 constant_modes.fill(
true);
349 return std::pair<Table<2, bool>, std::vector<unsigned int>>(
350 constant_modes, std::vector<unsigned int>(1, 0));
355template <
int dim,
int spacedim>
358 const unsigned int child,
366 if (this->prolongation[refinement_case - 1][child].n() == 0)
368 std::lock_guard<std::mutex> lock(prolongation_matrix_mutex);
371 if (this->prolongation[refinement_case - 1][child].n() ==
372 this->n_dofs_per_cell())
373 return this->prolongation[refinement_case - 1][child];
379 std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
381 isotropic_matrices.back().resize(
387 this_nonconst.prolongation[refinement_case - 1] =
388 std::move(isotropic_matrices.back());
392 return this->prolongation[refinement_case - 1][child];
397template <
int dim,
int spacedim>
400 const unsigned int child,
408 if (this->restriction[refinement_case - 1][child].n() == 0)
410 std::lock_guard<std::mutex> lock(restriction_matrix_mutex);
413 if (this->restriction[refinement_case - 1][child].n() ==
414 this->n_dofs_per_cell())
415 return this->restriction[refinement_case - 1][child];
421 std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
423 isotropic_matrices.back().resize(
429 this_nonconst.restriction[refinement_case - 1] =
430 std::move(isotropic_matrices.back());
434 return this->restriction[refinement_case - 1][child];
439template <
int dim,
int spacedim>
444 const unsigned int face_no)
const
458 const double eps = 2e-13 * this->degree * (dim - 1);
460 std::vector<Point<dim>> face_quadrature_points(quad_face_support.size());
464 face_quadrature_points);
467 for (
unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
469 double matrix_entry =
470 this->shape_value(this->face_to_cell_index(j, 0),
471 face_quadrature_points[i]);
476 if (std::fabs(matrix_entry - 1.0) < eps)
478 if (std::fabs(matrix_entry) < eps)
481 interpolation_matrix(i, j) = matrix_entry;
489 for (
unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
490 sum += interpolation_matrix(j, i);
504 spacedim>::ExcInterpolationNotImplemented()));
509template <
int dim,
int spacedim>
513 const unsigned int subface,
515 const unsigned int face_no)
const
529 const double eps = 2e-13 * this->degree * (dim - 1);
531 std::vector<Point<dim>> subface_quadrature_points(
532 quad_face_support.size());
537 subface_quadrature_points);
540 for (
unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
542 double matrix_entry =
543 this->shape_value(this->face_to_cell_index(j, 0),
544 subface_quadrature_points[i]);
549 if (std::fabs(matrix_entry - 1.0) < eps)
551 if (std::fabs(matrix_entry) < eps)
554 interpolation_matrix(i, j) = matrix_entry;
562 for (
unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
563 sum += interpolation_matrix(j, i);
577 spacedim>::ExcInterpolationNotImplemented()));
582template <
int dim,
int spacedim>
591template <
int dim,
int spacedim>
596 std::vector<double> &nodal_values)
const
599 this->get_unit_support_points().size());
603 for (
unsigned int i = 0; i < this->dofs_per_cell; ++i)
607 nodal_values[i] = support_point_values[i](0);
613template <
int dim,
int spacedim>
622 unit_support_points_fe_p<dim>(degree),
624 constraints_fe_p<dim>(degree))
636template <
int dim,
int spacedim>
637std::unique_ptr<FiniteElement<dim, spacedim>>
640 return std::make_unique<FE_SimplexP<dim, spacedim>>(*this);
645template <
int dim,
int spacedim>
649 std::ostringstream namebuf;
650 namebuf <<
"FE_SimplexP<" << dim <<
">(" << this->degree <<
")";
652 return namebuf.str();
657template <
int dim,
int spacedim>
661 const unsigned int codim)
const
682 if (this->degree < fe_p_other->degree)
684 else if (this->degree == fe_p_other->degree)
692 if (this->degree < fe_q_other->degree)
694 else if (this->degree == fe_q_other->degree)
702 if (fe_nothing->is_dominating())
717template <
int dim,
int spacedim>
718std::vector<std::pair<unsigned int, unsigned int>>
762template <
int dim,
int spacedim>
763std::vector<std::pair<unsigned int, unsigned int>>
779 std::vector<std::pair<unsigned int, unsigned int>> identities;
781 for (
unsigned int i = 0; i < this->degree - 1; ++i)
782 for (
unsigned int j = 0; j < fe_p_other->degree - 1; ++j)
783 if (std::fabs(this->unit_support_points[i + 3][0] -
784 fe_p_other->unit_support_points[i + 3][0]) < 1e-14)
785 identities.emplace_back(i, j);
813 const std::vector<unsigned int> &index_map_inverse_q_other =
814 fe_q_other->get_poly_space_numbering_inverse();
816 std::vector<std::pair<unsigned int, unsigned int>> identities;
818 for (
unsigned int i = 0; i < this->degree - 1; ++i)
819 for (
unsigned int j = 0; j < fe_q_other->degree - 1; ++j)
820 if (std::fabs(this->unit_support_points[i + 3][0] -
821 fe_q_other->get_unit_support_points()
822 [index_map_inverse_q_other[j + 1]][0]) < 1e-14)
823 identities.emplace_back(i, j);
866template <
int dim,
int spacedim>
875 unit_support_points_fe_p<dim>(degree),
877 constraints_fe_p<dim>(degree))
882template <
int dim,
int spacedim>
883std::unique_ptr<FiniteElement<dim, spacedim>>
886 return std::make_unique<FE_SimplexDGP<dim, spacedim>>(*this);
891template <
int dim,
int spacedim>
895 std::ostringstream namebuf;
896 namebuf <<
"FE_SimplexDGP<" << dim <<
">(" << this->degree <<
")";
898 return namebuf.str();
902template <
int dim,
int spacedim>
906 const unsigned int codim)
const
923 if (this->degree < fe_dgp_other->degree)
925 else if (this->degree == fe_dgp_other->degree)
933 if (this->degree < fe_dgq_other->degree)
935 else if (this->degree == fe_dgq_other->degree)
943 if (fe_nothing->is_dominating())
958template <
int dim,
int spacedim>
959std::vector<std::pair<unsigned int, unsigned int>>
970template <
int dim,
int spacedim>
971std::vector<std::pair<unsigned int, unsigned int>>
981#include "fe_simplex_p.inst"
static BarycentricPolynomials< dim > get_fe_p_basis(const unsigned int degree)
FiniteElementDomination::Domination compare_for_domination(const FiniteElement< dim, spacedim > &fe_other, const unsigned int codim) const override
std::string get_name() const override
std::vector< std::pair< unsigned int, unsigned int > > hp_line_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
std::vector< std::pair< unsigned int, unsigned int > > hp_vertex_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
FE_SimplexDGP(const unsigned int degree)
std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
FE_SimplexP(const unsigned int degree)
std::vector< std::pair< unsigned int, unsigned int > > hp_line_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
FiniteElementDomination::Domination compare_for_domination(const FiniteElement< dim, spacedim > &fe_other, const unsigned int codim) const override
std::vector< std::pair< unsigned int, unsigned int > > hp_vertex_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
std::string get_name() const override
std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
virtual std::pair< Table< 2, bool >, std::vector< unsigned int > > get_constant_modes() const override
void get_subface_interpolation_matrix(const FiniteElement< dim, spacedim > &x_source_fe, const unsigned int subface, FullMatrix< double > &interpolation_matrix, const unsigned int face_no) const override
virtual const FullMatrix< double > & get_restriction_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const override
virtual void convert_generalized_support_point_values_to_dof_values(const std::vector< Vector< double > > &support_point_values, std::vector< double > &nodal_values) const override
bool hp_constraints_are_implemented() const override
FE_SimplexPoly(const BarycentricPolynomials< dim > polynomials, const FiniteElementData< dim > &fe_data, const std::vector< Point< dim > > &unit_support_points, const std::vector< std::vector< Point< dim - 1 > > > unit_face_support_points, const FullMatrix< double > &interface_constraints)
virtual const FullMatrix< double > & get_prolongation_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const override
void get_face_interpolation_matrix(const FiniteElement< dim, spacedim > &source_fe, FullMatrix< double > &interpolation_matrix, const unsigned int face_no) const override
const unsigned int degree
unsigned int n_dofs_per_line() const
unsigned int n_dofs_per_face(unsigned int face_no=0, unsigned int child=0) const
unsigned int n_unique_faces() const
std::vector< std::vector< Point< dim - 1 > > > unit_face_support_points
const std::vector< Point< dim - 1 > > & get_unit_face_support_points(const unsigned int face_no=0) const
std::vector< int > adjust_line_dof_index_for_line_orientation_table
std::vector< Point< dim > > unit_support_points
FullMatrix< double > interface_constraints
static void project_to_subface(const ReferenceCell &reference_cell, const SubQuadrature &quadrature, const unsigned int face_no, const unsigned int subface_no, std::vector< Point< dim > > &q_points, const RefinementCase< dim - 1 > &ref_case=RefinementCase< dim - 1 >::isotropic_refinement)
static void project_to_face(const ReferenceCell &reference_cell, const SubQuadrature &quadrature, const unsigned int face_no, std::vector< Point< dim > > &q_points)
unsigned int n_vertices() const
static constexpr unsigned char default_combined_face_orientation()
unsigned int n_faces() const
unsigned int n_lines() const
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define AssertThrow(cond, exc)
#define DEAL_II_NOT_IMPLEMENTED()
@ either_element_can_dominate
@ other_element_dominates
@ neither_element_dominates
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
constexpr const ReferenceCell Tetrahedron
constexpr const ReferenceCell Triangle
constexpr const ReferenceCell Line
constexpr const ReferenceCell & get_simplex()