37 face_midpoint(
const unsigned int face_no)
40 const auto face_reference_cell =
44 for (
const auto face_vertex_no : face_reference_cell.
vertex_indices())
46 const auto vertex_no = reference_cell.face_to_cell_vertices(
49 midpoint += reference_cell.template vertex<dim>(vertex_no);
52 midpoint /= reference_cell.face_reference_cell(0).n_vertices();
60 std::vector<unsigned int>
61 get_dpo_vector_fe_p(
const unsigned int dim,
const unsigned int degree)
114 std::vector<Point<dim>>
115 unit_support_points_fe_p(
const unsigned int degree)
119 std::vector<Point<dim>> unit_points;
120 const auto reference_cell = ReferenceCells::get_simplex<dim>();
126 unit_points.emplace_back(reference_cell.template barycenter<dim>());
131 const auto dpo = get_dpo_vector_fe_p(dim, degree);
135 for (
const unsigned int d : reference_cell.vertex_indices())
136 unit_points.push_back(reference_cell.template vertex<dim>(d));
138 for (
const unsigned int l : reference_cell.line_indices())
141 unit_points[reference_cell.line_to_cell_vertices(l, 0)];
143 unit_points[reference_cell.line_to_cell_vertices(l, 1)];
144 for (
unsigned int p = 0; p < dpo[1]; ++p)
145 unit_points.push_back((
double(dpo[1] - p) / (dpo[1] + 1)) * p0 +
146 (double(p + 1) / (dpo[1] + 1)) * p1);
149 if (dim > 1 && dpo[2] > 0)
153 unit_points.push_back(reference_cell.template barycenter<dim>());
155 for (
const unsigned int f : reference_cell.face_indices())
156 unit_points.push_back(face_midpoint<dim>(f));
163 std::vector<Point<0>>
164 unit_support_points_fe_p(
const unsigned int )
174 std::vector<std::vector<
Point<dim - 1>>>
175 unit_face_support_points_fe_p(
176 const unsigned int degree,
187 std::vector<std::vector<
Point<dim - 1>>> unit_face_points;
193 unit_face_points.emplace_back(
194 unit_support_points_fe_p<dim - 1>(degree));
197 return unit_face_points;
207 constraints_fe_p(
const unsigned int )
216 constraints_fe_p<2>(
const unsigned int degree)
218 constexpr int dim = 2;
226 std::vector<
Point<dim - 1>> constraint_points;
228 constraint_points.emplace_back(0.5);
230 for (
unsigned int i = 1; i < degree; ++i)
231 constraint_points.push_back(
235 for (
unsigned int i = 1; i < degree; ++i)
236 constraint_points.push_back(
243 const unsigned int n_dofs_constrained = constraint_points.size();
244 unsigned int n_dofs_per_face = degree + 1;
250 for (
unsigned int i = 0; i < n_dofs_constrained; ++i)
251 for (
unsigned int j = 0; j < n_dofs_per_face; ++j)
253 interface_constraints(i, j) =
254 poly.compute_value(j, constraint_points[i]);
260 if (std::fabs(interface_constraints(i, j)) < 1e-13)
261 interface_constraints(i, j) = 0;
263 return interface_constraints;
272 std::vector<unsigned int>
273 get_dpo_vector_fe_dgp(
const unsigned int dim,
const unsigned int degree)
278 std::vector<unsigned int> dpo(dim + 1, 0U);
289 const auto continuous_dpo = get_dpo_vector_fe_p(dim, degree);
296 continuous_dpo[dim]};
304 continuous_dpo[dim]};
314 continuous_dpo[dim]};
325template <
int dim,
int spacedim>
329 const bool prolongation_is_additive,
330 const std::vector<
Point<dim>> &unit_support_points,
331 const std::vector<std::vector<
Point<dim - 1>>> unit_face_support_points,
336 std::vector<
bool>(fe_data.dofs_per_cell, prolongation_is_additive),
347template <
int dim,
int spacedim>
348std::pair<Table<2, bool>, std::vector<unsigned int>>
352 constant_modes.fill(
true);
353 return std::pair<Table<2, bool>, std::vector<unsigned int>>(
354 constant_modes, std::vector<unsigned int>(1, 0));
359template <
int dim,
int spacedim>
362 const unsigned int child,
383 if (this->prolongation[refinement_case - 1][child].n() == 0)
385 std::lock_guard<std::mutex> lock(prolongation_matrix_mutex);
388 if (this->prolongation[refinement_case - 1][child].n() ==
389 this->n_dofs_per_cell())
390 return this->prolongation[refinement_case - 1][child];
398 std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
400 isotropic_matrices.back().resize(
401 this->reference_cell().n_children(
404 this->n_dofs_per_cell()));
408 this_nonconst.prolongation[refinement_case - 1] =
409 std::move(isotropic_matrices.back());
413 std::vector<std::vector<FullMatrix<double>>> matrices(
416 this->reference_cell().n_children(
419 this->n_dofs_per_cell())));
421 for (
unsigned int refinement_direction =
static_cast<unsigned int>(
423 refinement_direction <=
425 refinement_direction++)
426 this_nonconst.prolongation[refinement_direction - 1] =
427 std::move(matrices[refinement_direction - 1]);
434 return this->prolongation[refinement_case - 1][child];
439template <
int dim,
int spacedim>
442 const unsigned int child,
462 if (this->restriction[refinement_case - 1][child].n() == 0)
464 std::lock_guard<std::mutex> lock(restriction_matrix_mutex);
467 if (this->restriction[refinement_case - 1][child].n() ==
468 this->n_dofs_per_cell())
469 return this->restriction[refinement_case - 1][child];
479 const double eps = 1e-12;
481 this->n_dofs_per_cell());
484 const std::vector<Point<dim>> unit_support_points =
485 this->get_unit_support_points();
493 tria.
begin_active()->set_refine_choice(refinement_case);
496 const auto &child_cell = tria.
begin(0)->child(child);
500 for (
unsigned int i = 0; i < unit_support_points.size(); i++)
502 std::vector<Point<dim>> transformed_point(1);
503 const std::vector<Point<spacedim>> unit_support_point = {
505 unit_support_points[i][1]) :
507 unit_support_points[i][1],
508 unit_support_points[i][2])};
509 this->reference_cell()
510 .template get_default_linear_mapping<dim, spacedim>()
511 .transform_points_real_to_unit_cell(
517 if (this->reference_cell().contains_point(transformed_point[0], eps))
518 for (
unsigned int j = 0; j < this->n_dofs_per_cell(); j++)
519 restriction_mat[i][j] =
520 this->shape_value(j, transformed_point[0]);
524 for (
unsigned int i = 0; i < this->n_dofs_per_cell(); i++)
528 for (
unsigned int j = 0; j < this->n_dofs_per_cell(); j++)
529 sum += restriction_mat[i][j];
531 Assert(std::fabs(sum - 1) < eps || std::fabs(sum) < eps,
533 "The entries in a row of the local "
534 "restriction matrix do not add to zero or one. "
535 "This typically indicates that the "
536 "polynomial interpolation is "
537 "ill-conditioned such that round-off "
538 "prevents the sum to be one."));
543 for (
unsigned int i = 0; i < restriction_mat.
m(); ++i)
544 for (
unsigned int j = 0; j < restriction_mat.
n(); ++j)
546 if (std::fabs(restriction_mat(i, j)) < eps)
547 restriction_mat(i, j) = 0.;
548 if (std::fabs(restriction_mat(i, j) - 1) < eps)
549 restriction_mat(i, j) = 1.;
553 this->restriction[refinement_case - 1][child]) =
554 std::move(restriction_mat);
558 return this->restriction[refinement_case - 1][child];
563template <
int dim,
int spacedim>
568 const unsigned int face_no)
const
582 const double eps = 2e-13 * this->degree * (dim - 1);
584 const std::vector<Point<dim>> face_quadrature_points =
592 for (
unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
594 double matrix_entry =
595 this->shape_value(this->face_to_cell_index(j, 0),
596 face_quadrature_points[i]);
601 if (std::fabs(matrix_entry - 1.0) < eps)
603 if (std::fabs(matrix_entry) < eps)
606 interpolation_matrix(i, j) = matrix_entry;
615 for (
unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
616 sum += interpolation_matrix(j, i);
630 spacedim>::ExcInterpolationNotImplemented()));
635template <
int dim,
int spacedim>
639 const unsigned int subface,
641 const unsigned int face_no)
const
655 const double eps = 2e-13 * this->degree * (dim - 1);
659 this->reference_cell(),
667 for (
unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
669 double matrix_entry =
670 this->shape_value(this->face_to_cell_index(j, 0),
671 subface_quadrature.point(i));
676 if (std::fabs(matrix_entry - 1.0) < eps)
678 if (std::fabs(matrix_entry) < eps)
681 interpolation_matrix(i, j) = matrix_entry;
690 for (
unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
691 sum += interpolation_matrix(j, i);
705 spacedim>::ExcInterpolationNotImplemented()));
710template <
int dim,
int spacedim>
719template <
int dim,
int spacedim>
724 std::vector<double> &nodal_values)
const
727 this->get_unit_support_points().size());
731 for (
unsigned int i = 0; i < this->dofs_per_cell; ++i)
735 nodal_values[i] = support_point_values[i](0);
741template <
int dim,
int spacedim>
751 unit_support_points_fe_p<dim>(degree),
753 constraints_fe_p<dim>(degree))
765template <
int dim,
int spacedim>
766std::unique_ptr<FiniteElement<dim, spacedim>>
769 return std::make_unique<FE_SimplexP<dim, spacedim>>(*this);
774template <
int dim,
int spacedim>
778 std::ostringstream namebuf;
780 << this->degree <<
")";
782 return namebuf.str();
787template <
int dim,
int spacedim>
791 const unsigned int codim)
const
812 if (this->degree < fe_p_other->degree)
814 else if (this->degree == fe_p_other->degree)
822 if (this->degree < fe_q_other->degree)
824 else if (this->degree == fe_q_other->degree)
832 if (fe_nothing->is_dominating())
847template <
int dim,
int spacedim>
848std::vector<std::pair<unsigned int, unsigned int>>
892template <
int dim,
int spacedim>
893std::vector<std::pair<unsigned int, unsigned int>>
909 std::vector<std::pair<unsigned int, unsigned int>> identities;
911 for (
unsigned int i = 0; i < this->degree - 1; ++i)
912 for (
unsigned int j = 0; j < fe_p_other->degree - 1; ++j)
913 if (std::fabs(this->unit_support_points[i + 3][0] -
914 fe_p_other->unit_support_points[i + 3][0]) < 1e-14)
915 identities.emplace_back(i, j);
943 const std::vector<unsigned int> &index_map_inverse_q_other =
944 fe_q_other->get_poly_space_numbering_inverse();
946 std::vector<std::pair<unsigned int, unsigned int>> identities;
948 for (
unsigned int i = 0; i < this->degree - 1; ++i)
949 for (
unsigned int j = 0; j < fe_q_other->degree - 1; ++j)
950 if (std::fabs(this->unit_support_points[i + 3][0] -
951 fe_q_other->get_unit_support_points()
952 [index_map_inverse_q_other[j + 1]][0]) < 1e-14)
953 identities.emplace_back(i, j);
996template <
int dim,
int spacedim>
1006 unit_support_points_fe_p<dim>(degree),
1008 constraints_fe_p<dim>(degree))
1013template <
int dim,
int spacedim>
1014std::unique_ptr<FiniteElement<dim, spacedim>>
1017 return std::make_unique<FE_SimplexDGP<dim, spacedim>>(*this);
1022template <
int dim,
int spacedim>
1026 std::ostringstream namebuf;
1028 << this->degree <<
")";
1030 return namebuf.str();
1034template <
int dim,
int spacedim>
1038 const unsigned int codim)
const
1055 if (this->degree < fe_dgp_other->degree)
1057 else if (this->degree == fe_dgp_other->degree)
1065 if (this->degree < fe_dgq_other->degree)
1067 else if (this->degree == fe_dgq_other->degree)
1075 if (fe_nothing->is_dominating())
1090template <
int dim,
int spacedim>
1091std::vector<std::pair<unsigned int, unsigned int>>
1102template <
int dim,
int spacedim>
1103std::vector<std::pair<unsigned int, unsigned int>>
1114template <
int dim,
int spacedim>
1117 const unsigned int child,
1137 if (this->restriction[refinement_case - 1][child].n() == 0)
1139 std::lock_guard<std::mutex> lock(this->restriction_matrix_mutex);
1142 if (this->restriction[refinement_case - 1][child].n() ==
1143 this->n_dofs_per_cell())
1144 return this->restriction[refinement_case - 1][child];
1152 std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
1154 isotropic_matrices.back().resize(
1155 this->reference_cell().n_children(
1158 this->n_dofs_per_cell()));
1162 this_nonconst.restriction[refinement_case - 1] =
1163 std::move(isotropic_matrices.back());
1167 std::vector<std::vector<FullMatrix<double>>> matrices(
1170 this->reference_cell().n_children(
1173 this->n_dofs_per_cell())));
1175 for (
unsigned int refinement_direction =
static_cast<unsigned int>(
1177 refinement_direction <=
1179 refinement_direction++)
1180 this_nonconst.restriction[refinement_direction - 1] =
1181 std::move(matrices[refinement_direction - 1]);
1188 return this->restriction[refinement_case - 1][child];
1192#include "fe/fe_simplex_p.inst"
ArrayView< std::remove_reference_t< typename std::iterator_traits< Iterator >::reference >, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
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)
virtual const FullMatrix< double > & get_restriction_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const override
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
FE_SimplexPoly(const BarycentricPolynomials< dim > polynomials, const FiniteElementData< dim > &fe_data, const bool prolongation_is_additive, 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 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
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
unsigned int n_faces() const
unsigned int n_lines() const
cell_iterator begin(const unsigned int level=0) const
virtual void execute_coarsening_and_refinement()
active_cell_iterator begin_active(const unsigned int level=0) const
#define DEAL_II_NAMESPACE_OPEN
constexpr bool running_in_debug_mode()
#define DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
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)
@ 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()
std::string dim_string(const int dim, const int spacedim)
constexpr types::geometric_orientation default_geometric_orientation