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(
51 midpoint += reference_cell.template vertex<dim>(vertex_no);
54 midpoint /= reference_cell.face_reference_cell(0).n_vertices();
62 std::vector<unsigned int>
63 get_dpo_vector_fe_p(
const unsigned int dim,
const unsigned int degree)
116 std::vector<Point<dim>>
117 unit_support_points_fe_p(
const unsigned int degree)
121 std::vector<Point<dim>> unit_points;
122 const auto reference_cell = ReferenceCells::get_simplex<dim>();
128 unit_points.emplace_back(reference_cell.template barycenter<dim>());
133 const auto dpo = get_dpo_vector_fe_p(dim, degree);
146 for (
unsigned int p = 0; p < dpo[1]; ++p)
147 unit_points.push_back((
double(dpo[1] - p) / (dpo[1] + 1)) * p0 +
148 (
double(p + 1) / (dpo[1] + 1)) * p1);
151 if (dim > 1 && dpo[2] > 0)
158 unit_points.push_back(face_midpoint<dim>(f));
165 std::vector<Point<0>>
166 unit_support_points_fe_p(
const unsigned int )
176 std::vector<std::vector<
Point<dim - 1>>>
177 unit_face_support_points_fe_p(
178 const unsigned int degree,
189 std::vector<std::vector<
Point<dim - 1>>> unit_face_points;
195 unit_face_points.emplace_back(
196 unit_support_points_fe_p<dim - 1>(degree));
199 return unit_face_points;
209 constraints_fe_p(
const unsigned int )
218 constraints_fe_p<2>(
const unsigned int degree)
220 constexpr int dim = 2;
229 std::vector<
Point<dim - 1>> constraint_points;
231 constraint_points.emplace_back(0.5);
235 constraint_points.emplace_back(0.25);
237 constraint_points.emplace_back(0.75);
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]);
523 for (
unsigned int i = 0; i < this->n_dofs_per_cell(); i++)
527 for (
unsigned int j = 0; j < this->n_dofs_per_cell(); j++)
528 sum += restriction_mat[i][j];
530 Assert(std::fabs(sum - 1) < eps || std::fabs(sum) < eps,
532 "The entries in a row of the local "
533 "restriction matrix do not add to zero or one. "
534 "This typically indicates that the "
535 "polynomial interpolation is "
536 "ill-conditioned such that round-off "
537 "prevents the sum to be one."));
542 for (
unsigned int i = 0; i < restriction_mat.
m(); ++i)
543 for (
unsigned int j = 0; j < restriction_mat.
n(); ++j)
545 if (std::fabs(restriction_mat(i, j)) < eps)
546 restriction_mat(i, j) = 0.;
547 if (std::fabs(restriction_mat(i, j) - 1) < eps)
548 restriction_mat(i, j) = 1.;
552 this->restriction[refinement_case - 1][child]) =
553 std::move(restriction_mat);
557 return this->restriction[refinement_case - 1][child];
562template <
int dim,
int spacedim>
567 const unsigned int face_no)
const
581 const double eps = 2e-13 * this->degree * (dim - 1);
583 std::vector<Point<dim>> face_quadrature_points(quad_face_support.size());
587 face_quadrature_points);
590 for (
unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
592 double matrix_entry =
593 this->shape_value(this->face_to_cell_index(j, 0),
594 face_quadrature_points[i]);
599 if (std::fabs(matrix_entry - 1.0) < eps)
601 if (std::fabs(matrix_entry) < eps)
604 interpolation_matrix(i, j) = matrix_entry;
612 for (
unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
613 sum += interpolation_matrix(j, i);
627 spacedim>::ExcInterpolationNotImplemented()));
632template <
int dim,
int spacedim>
636 const unsigned int subface,
638 const unsigned int face_no)
const
652 const double eps = 2e-13 * this->degree * (dim - 1);
654 std::vector<Point<dim>> subface_quadrature_points(
655 quad_face_support.size());
660 subface_quadrature_points);
663 for (
unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
665 double matrix_entry =
666 this->shape_value(this->face_to_cell_index(j, 0),
667 subface_quadrature_points[i]);
672 if (std::fabs(matrix_entry - 1.0) < eps)
674 if (std::fabs(matrix_entry) < eps)
677 interpolation_matrix(i, j) = matrix_entry;
685 for (
unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
686 sum += interpolation_matrix(j, i);
700 spacedim>::ExcInterpolationNotImplemented()));
705template <
int dim,
int spacedim>
714template <
int dim,
int spacedim>
719 std::vector<double> &nodal_values)
const
722 this->get_unit_support_points().size());
726 for (
unsigned int i = 0; i < this->dofs_per_cell; ++i)
730 nodal_values[i] = support_point_values[i](0);
736template <
int dim,
int spacedim>
746 unit_support_points_fe_p<dim>(degree),
748 constraints_fe_p<dim>(degree))
760template <
int dim,
int spacedim>
761std::unique_ptr<FiniteElement<dim, spacedim>>
764 return std::make_unique<FE_SimplexP<dim, spacedim>>(*this);
769template <
int dim,
int spacedim>
773 std::ostringstream namebuf;
775 << this->degree <<
")";
777 return namebuf.str();
782template <
int dim,
int spacedim>
786 const unsigned int codim)
const
807 if (this->degree < fe_p_other->degree)
809 else if (this->degree == fe_p_other->degree)
817 if (this->degree < fe_q_other->degree)
819 else if (this->degree == fe_q_other->degree)
827 if (fe_nothing->is_dominating())
842template <
int dim,
int spacedim>
843std::vector<std::pair<unsigned int, unsigned int>>
887template <
int dim,
int spacedim>
888std::vector<std::pair<unsigned int, unsigned int>>
904 std::vector<std::pair<unsigned int, unsigned int>> identities;
906 for (
unsigned int i = 0; i < this->degree - 1; ++i)
907 for (
unsigned int j = 0; j < fe_p_other->degree - 1; ++j)
908 if (std::fabs(this->unit_support_points[i + 3][0] -
909 fe_p_other->unit_support_points[i + 3][0]) < 1e-14)
910 identities.emplace_back(i, j);
938 const std::vector<unsigned int> &index_map_inverse_q_other =
939 fe_q_other->get_poly_space_numbering_inverse();
941 std::vector<std::pair<unsigned int, unsigned int>> identities;
943 for (
unsigned int i = 0; i < this->degree - 1; ++i)
944 for (
unsigned int j = 0; j < fe_q_other->degree - 1; ++j)
945 if (std::fabs(this->unit_support_points[i + 3][0] -
946 fe_q_other->get_unit_support_points()
947 [index_map_inverse_q_other[j + 1]][0]) < 1e-14)
948 identities.emplace_back(i, j);
991template <
int dim,
int spacedim>
1001 unit_support_points_fe_p<dim>(degree),
1003 constraints_fe_p<dim>(degree))
1008template <
int dim,
int spacedim>
1009std::unique_ptr<FiniteElement<dim, spacedim>>
1012 return std::make_unique<FE_SimplexDGP<dim, spacedim>>(*this);
1017template <
int dim,
int spacedim>
1021 std::ostringstream namebuf;
1023 << this->degree <<
")";
1025 return namebuf.str();
1029template <
int dim,
int spacedim>
1033 const unsigned int codim)
const
1050 if (this->degree < fe_dgp_other->degree)
1052 else if (this->degree == fe_dgp_other->degree)
1060 if (this->degree < fe_dgq_other->degree)
1062 else if (this->degree == fe_dgq_other->degree)
1070 if (fe_nothing->is_dominating())
1085template <
int dim,
int spacedim>
1086std::vector<std::pair<unsigned int, unsigned int>>
1097template <
int dim,
int spacedim>
1098std::vector<std::pair<unsigned int, unsigned int>>
1109template <
int dim,
int spacedim>
1112 const unsigned int child,
1132 if (this->restriction[refinement_case - 1][child].n() == 0)
1134 std::lock_guard<std::mutex> lock(this->restriction_matrix_mutex);
1137 if (this->restriction[refinement_case - 1][child].n() ==
1138 this->n_dofs_per_cell())
1139 return this->restriction[refinement_case - 1][child];
1147 std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
1149 isotropic_matrices.back().resize(
1150 this->reference_cell().n_children(
1153 this->n_dofs_per_cell()));
1157 this_nonconst.restriction[refinement_case - 1] =
1158 std::move(isotropic_matrices.back());
1162 std::vector<std::vector<FullMatrix<double>>> matrices(
1165 this->reference_cell().n_children(
1168 this->n_dofs_per_cell())));
1170 for (
unsigned int refinement_direction =
static_cast<unsigned int>(
1172 refinement_direction <=
1174 refinement_direction++)
1175 this_nonconst.restriction[refinement_direction - 1] =
1176 std::move(matrices[refinement_direction - 1]);
1183 return this->restriction[refinement_case - 1][child];
1187#include "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
static constexpr unsigned char default_combined_face_orientation()
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
#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_ASSERT_UNREACHABLE()
#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)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
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)