40 face_midpoint(
const unsigned int face_no)
43 const auto face_reference_cell =
44 reference_cell.face_reference_cell(face_no);
47 for (
const auto face_vertex_no : face_reference_cell.
vertex_indices())
49 const auto vertex_no = reference_cell.face_to_cell_vertices(
52 midpoint += reference_cell.template vertex<dim>(vertex_no);
55 midpoint /= reference_cell.face_reference_cell(0).n_vertices();
63 std::vector<unsigned int>
64 get_dpo_vector_fe_p(
const unsigned int dim,
const unsigned int degree)
117 std::vector<Point<dim>>
118 unit_support_points_fe_p(
const unsigned int degree)
122 std::vector<Point<dim>> unit_points;
123 const auto reference_cell = ReferenceCells::get_simplex<dim>();
129 unit_points.emplace_back(reference_cell.template barycenter<dim>());
134 const auto dpo = get_dpo_vector_fe_p(dim, degree);
139 unit_points.push_back(reference_cell.template vertex<dim>(d));
141 for (
const unsigned int l : reference_cell.line_indices())
144 unit_points[reference_cell.line_to_cell_vertices(l, 0)];
146 unit_points[reference_cell.line_to_cell_vertices(l, 1)];
147 for (
unsigned int p = 0; p < dpo[1]; ++p)
148 unit_points.push_back((
double(dpo[1] - p) / (dpo[1] + 1)) * p0 +
149 (
double(p + 1) / (dpo[1] + 1)) * p1);
152 if (dim > 1 && dpo[2] > 0)
159 unit_points.push_back(face_midpoint<dim>(f));
166 std::vector<Point<0>>
167 unit_support_points_fe_p(
const unsigned int )
177 std::vector<std::vector<
Point<dim - 1>>>
178 unit_face_support_points_fe_p(
179 const unsigned int degree,
190 std::vector<std::vector<
Point<dim - 1>>> unit_face_points;
196 unit_face_points.emplace_back(
197 unit_support_points_fe_p<dim - 1>(degree));
200 return unit_face_points;
210 constraints_fe_p(
const unsigned int )
219 constraints_fe_p<2>(
const unsigned int degree)
221 constexpr int dim = 2;
229 std::vector<
Point<dim - 1>> constraint_points;
231 constraint_points.emplace_back(0.5);
233 for (
unsigned int i = 1; i < degree; ++i)
234 constraint_points.push_back(
238 for (
unsigned int i = 1; i < degree; ++i)
239 constraint_points.push_back(
246 const unsigned int n_dofs_constrained = constraint_points.size();
247 unsigned int n_dofs_per_face = degree + 1;
253 for (
unsigned int i = 0; i < n_dofs_constrained; ++i)
254 for (
unsigned int j = 0; j < n_dofs_per_face; ++j)
256 interface_constraints(i, j) =
257 poly.compute_value(j, constraint_points[i]);
263 if (std::fabs(interface_constraints(i, j)) < 1e-13)
264 interface_constraints(i, j) = 0;
266 return interface_constraints;
275 std::vector<unsigned int>
276 get_dpo_vector_fe_dgp(
const unsigned int dim,
const unsigned int degree)
281 std::vector<unsigned int> dpo(dim + 1, 0U);
292 const auto continuous_dpo = get_dpo_vector_fe_p(dim, degree);
299 continuous_dpo[dim]};
307 continuous_dpo[dim]};
317 continuous_dpo[dim]};
328template <
int dim,
int spacedim>
332 const bool prolongation_is_additive,
333 const std::vector<
Point<dim>> &unit_support_points,
334 const std::vector<std::vector<
Point<dim - 1>>> unit_face_support_points,
339 std::vector<
bool>(fe_data.dofs_per_cell, prolongation_is_additive),
350template <
int dim,
int spacedim>
351std::pair<Table<2, bool>, std::vector<unsigned int>>
355 constant_modes.fill(
true);
356 return std::pair<Table<2, bool>, std::vector<unsigned int>>(
357 constant_modes, std::vector<unsigned int>(1, 0));
362template <
int dim,
int spacedim>
365 const unsigned int child,
386 if (this->prolongation[refinement_case - 1][child].n() == 0)
388 std::lock_guard<std::mutex> lock(prolongation_matrix_mutex);
391 if (this->prolongation[refinement_case - 1][child].n() ==
392 this->n_dofs_per_cell())
393 return this->prolongation[refinement_case - 1][child];
401 std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
403 isotropic_matrices.back().resize(
404 this->reference_cell().n_children(
407 this->n_dofs_per_cell()));
411 this_nonconst.prolongation[refinement_case - 1] =
412 std::move(isotropic_matrices.back());
416 std::vector<std::vector<FullMatrix<double>>> matrices(
419 this->reference_cell().n_children(
422 this->n_dofs_per_cell())));
424 for (
unsigned int refinement_direction =
static_cast<unsigned int>(
426 refinement_direction <=
428 refinement_direction++)
429 this_nonconst.prolongation[refinement_direction - 1] =
430 std::move(matrices[refinement_direction - 1]);
437 return this->prolongation[refinement_case - 1][child];
442template <
int dim,
int spacedim>
445 const unsigned int face_dof_index,
446 const unsigned int face,
457 ExcMessage(
"In 1D, all faces must have the default orientation."));
460 this->reference_cell().n_face_orientations(face));
464 if (face_dof_index < this->get_first_face_line_index(face))
469 const unsigned int face_vertex =
470 face_dof_index / this->n_dofs_per_vertex();
471 const unsigned int dof_index_on_vertex =
472 face_dof_index % this->n_dofs_per_vertex();
476 return this->reference_cell().face_to_cell_vertices(
477 face, face_vertex, combined_orientation) *
478 this->n_dofs_per_vertex() +
481 else if (face_dof_index < this->get_first_face_quad_index(face))
486 const unsigned int index =
487 face_dof_index - this->get_first_face_line_index(face);
489 const unsigned int face_line = index / this->n_dofs_per_line();
490 const unsigned int dof_index_on_line = index % this->n_dofs_per_line();
494 unsigned int adjusted_dof_index_on_line = 0;
503 adjusted_dof_index_on_line = dof_index_on_line;
505 adjusted_dof_index_on_line =
506 this->n_dofs_per_line() - dof_index_on_line - 1;
526 Assert((this->n_dofs_per_line() <= 1) ||
527 combined_orientation ==
529 combined_orientation ==
535 adjusted_dof_index_on_line = dof_index_on_line;
539 adjusted_dof_index_on_line =
540 this->n_dofs_per_line() - dof_index_on_line - 1;
548 return (this->get_first_line_index() +
549 this->reference_cell().face_to_cell_lines(face,
551 combined_orientation) *
552 this->n_dofs_per_line() +
553 adjusted_dof_index_on_line);
561 const unsigned int index =
562 face_dof_index - this->get_first_face_quad_index(face);
567 Assert((this->n_dofs_per_quad(face) <= 1) ||
570 return (this->get_first_quad_index(face) + index);
576template <
int dim,
int spacedim>
579 const unsigned int child,
599 if (this->restriction[refinement_case - 1][child].n() == 0)
601 std::lock_guard<std::mutex> lock(restriction_matrix_mutex);
604 if (this->restriction[refinement_case - 1][child].n() ==
605 this->n_dofs_per_cell())
606 return this->restriction[refinement_case - 1][child];
616 const double eps = 1e-12;
618 this->n_dofs_per_cell());
621 const std::vector<Point<dim>> unit_support_points =
622 this->get_unit_support_points();
630 tria.
begin_active()->set_refine_choice(refinement_case);
633 const auto &child_cell = tria.
begin(0)->child(child);
637 for (
unsigned int i = 0; i < unit_support_points.size(); i++)
639 std::vector<Point<dim>> transformed_point(1);
640 const std::vector<Point<spacedim>> unit_support_point = {
642 unit_support_points[i][1]) :
644 unit_support_points[i][1],
645 unit_support_points[i][2])};
646 this->reference_cell()
647 .template get_default_linear_mapping<dim, spacedim>()
648 .transform_points_real_to_unit_cell(
654 if (this->reference_cell().contains_point(transformed_point[0], eps))
655 for (
unsigned int j = 0; j < this->n_dofs_per_cell(); j++)
656 restriction_mat[i][j] =
657 this->shape_value(j, transformed_point[0]);
661 for (
unsigned int i = 0; i < this->n_dofs_per_cell(); i++)
665 for (
unsigned int j = 0; j < this->n_dofs_per_cell(); j++)
666 sum += restriction_mat[i][j];
668 Assert(std::fabs(sum - 1) < eps || std::fabs(sum) < eps,
670 "The entries in a row of the local "
671 "restriction matrix do not add to zero or one. "
672 "This typically indicates that the "
673 "polynomial interpolation is "
674 "ill-conditioned such that round-off "
675 "prevents the sum to be one."));
680 for (
unsigned int i = 0; i < restriction_mat.
m(); ++i)
681 for (
unsigned int j = 0; j < restriction_mat.
n(); ++j)
683 if (std::fabs(restriction_mat(i, j)) < eps)
684 restriction_mat(i, j) = 0.;
685 if (std::fabs(restriction_mat(i, j) - 1) < eps)
686 restriction_mat(i, j) = 1.;
690 this->restriction[refinement_case - 1][child]) =
691 std::move(restriction_mat);
695 return this->restriction[refinement_case - 1][child];
700template <
int dim,
int spacedim>
705 const unsigned int face_no)
const
719 const double eps = 2e-13 * this->degree * (dim - 1);
721 const std::vector<Point<dim>> face_quadrature_points =
729 for (
unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
731 double matrix_entry =
732 this->shape_value(this->face_to_cell_index(j, 0),
733 face_quadrature_points[i]);
738 if (std::fabs(matrix_entry - 1.0) < eps)
740 if (std::fabs(matrix_entry) < eps)
743 interpolation_matrix(i, j) = matrix_entry;
752 for (
unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
753 sum += interpolation_matrix(j, i);
767 spacedim>::ExcInterpolationNotImplemented()));
772template <
int dim,
int spacedim>
776 const unsigned int subface,
778 const unsigned int face_no)
const
792 const double eps = 2e-13 * this->degree * (dim - 1);
796 this->reference_cell(),
804 for (
unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
806 double matrix_entry =
807 this->shape_value(this->face_to_cell_index(j, 0),
808 subface_quadrature.point(i));
813 if (std::fabs(matrix_entry - 1.0) < eps)
815 if (std::fabs(matrix_entry) < eps)
818 interpolation_matrix(i, j) = matrix_entry;
827 for (
unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
828 sum += interpolation_matrix(j, i);
842 spacedim>::ExcInterpolationNotImplemented()));
847template <
int dim,
int spacedim>
856template <
int dim,
int spacedim>
861 std::vector<double> &nodal_values)
const
864 this->get_unit_support_points().size());
868 for (
unsigned int i = 0; i < this->dofs_per_cell; ++i)
872 nodal_values[i] = support_point_values[i](0);
878template <
int dim,
int spacedim>
888 unit_support_points_fe_p<dim>(degree),
890 constraints_fe_p<dim>(degree))
902template <
int dim,
int spacedim>
903std::unique_ptr<FiniteElement<dim, spacedim>>
906 return std::make_unique<FE_SimplexP<dim, spacedim>>(*this);
911template <
int dim,
int spacedim>
915 std::ostringstream namebuf;
917 << this->degree <<
")";
919 return namebuf.str();
924template <
int dim,
int spacedim>
928 const unsigned int codim)
const
949 if (this->degree < fe_p_other->degree)
951 else if (this->degree == fe_p_other->degree)
959 if (this->degree < fe_q_other->degree)
961 else if (this->degree == fe_q_other->degree)
969 if (fe_nothing->is_dominating())
984template <
int dim,
int spacedim>
985std::vector<std::pair<unsigned int, unsigned int>>
1029template <
int dim,
int spacedim>
1030std::vector<std::pair<unsigned int, unsigned int>>
1046 std::vector<std::pair<unsigned int, unsigned int>> identities;
1048 for (
unsigned int i = 0; i < this->degree - 1; ++i)
1049 for (
unsigned int j = 0; j < fe_p_other->degree - 1; ++j)
1050 if (std::fabs(this->unit_support_points[i + 3][0] -
1051 fe_p_other->unit_support_points[i + 3][0]) < 1e-14)
1052 identities.emplace_back(i, j);
1080 const std::vector<unsigned int> &index_map_inverse_q_other =
1081 fe_q_other->get_poly_space_numbering_inverse();
1083 std::vector<std::pair<unsigned int, unsigned int>> identities;
1085 for (
unsigned int i = 0; i < this->degree - 1; ++i)
1086 for (
unsigned int j = 0; j < fe_q_other->degree - 1; ++j)
1087 if (std::fabs(this->unit_support_points[i + 3][0] -
1088 fe_q_other->get_unit_support_points()
1089 [index_map_inverse_q_other[j + 1]][0]) < 1e-14)
1090 identities.emplace_back(i, j);
1133template <
int dim,
int spacedim>
1143 unit_support_points_fe_p<dim>(degree),
1145 constraints_fe_p<dim>(degree))
1150template <
int dim,
int spacedim>
1151std::unique_ptr<FiniteElement<dim, spacedim>>
1154 return std::make_unique<FE_SimplexDGP<dim, spacedim>>(*this);
1159template <
int dim,
int spacedim>
1163 std::ostringstream namebuf;
1165 << this->degree <<
")";
1167 return namebuf.str();
1171template <
int dim,
int spacedim>
1175 const unsigned int codim)
const
1192 if (this->degree < fe_dgp_other->degree)
1194 else if (this->degree == fe_dgp_other->degree)
1202 if (this->degree < fe_dgq_other->degree)
1204 else if (this->degree == fe_dgq_other->degree)
1212 if (fe_nothing->is_dominating())
1227template <
int dim,
int spacedim>
1228std::vector<std::pair<unsigned int, unsigned int>>
1239template <
int dim,
int spacedim>
1240std::vector<std::pair<unsigned int, unsigned int>>
1251template <
int dim,
int spacedim>
1254 const unsigned int child,
1274 if (this->restriction[refinement_case - 1][child].n() == 0)
1276 std::lock_guard<std::mutex> lock(this->restriction_matrix_mutex);
1279 if (this->restriction[refinement_case - 1][child].n() ==
1280 this->n_dofs_per_cell())
1281 return this->restriction[refinement_case - 1][child];
1289 std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
1291 isotropic_matrices.back().resize(
1292 this->reference_cell().n_children(
1295 this->n_dofs_per_cell()));
1299 this_nonconst.restriction[refinement_case - 1] =
1300 std::move(isotropic_matrices.back());
1304 std::vector<std::vector<FullMatrix<double>>> matrices(
1307 this->reference_cell().n_children(
1310 this->n_dofs_per_cell())));
1312 for (
unsigned int refinement_direction =
static_cast<unsigned int>(
1314 refinement_direction <=
1316 refinement_direction++)
1317 this_nonconst.restriction[refinement_direction - 1] =
1318 std::move(matrices[refinement_direction - 1]);
1325 return this->restriction[refinement_case - 1][child];
1329#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
virtual unsigned int face_to_cell_index(const unsigned int face_dof_index, const unsigned int face, const types::geometric_orientation combined_orientation=numbers::default_geometric_orientation) 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)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
#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 ReferenceCell Triangle
constexpr ReferenceCell Tetrahedron
constexpr const ReferenceCell & get_simplex()
constexpr ReferenceCell Line
std::string dim_string(const int dim, const int spacedim)
types::geometric_orientation combined_face_orientation(const bool face_orientation, const bool face_rotation, const bool face_flip)
constexpr types::geometric_orientation default_geometric_orientation
unsigned char geometric_orientation