35 std::vector<unsigned int>
36 get_dpo_vector_fe_p(
const unsigned int dim,
const unsigned int degree)
38 std::vector<unsigned int> dpo(dim + 1, 0
U);
65 std::vector<Point<dim>>
66 unit_support_points_fe_poly(
const unsigned int degree)
68 std::vector<Point<dim>> unit_points;
77 1.0 /
double(dim + 1));
78 unit_points.emplace_back(centroid);
89 unit_points.emplace_back(0.0);
90 unit_points.emplace_back(1.0);
93 unit_points.emplace_back(0.5);
101 unit_points.emplace_back(0.0, 0.0);
102 unit_points.emplace_back(1.0, 0.0);
103 unit_points.emplace_back(0.0, 1.0);
107 unit_points.emplace_back(0.5, 0.0);
108 unit_points.emplace_back(0.5, 0.5);
109 unit_points.emplace_back(0.0, 0.5);
118 unit_points.emplace_back(0.0, 0.0, 0.0);
119 unit_points.emplace_back(1.0, 0.0, 0.0);
120 unit_points.emplace_back(0.0, 1.0, 0.0);
121 unit_points.emplace_back(0.0, 0.0, 1.0);
125 unit_points.emplace_back(0.5, 0.0, 0.0);
126 unit_points.emplace_back(0.5, 0.5, 0.0);
127 unit_points.emplace_back(0.0, 0.5, 0.0);
128 unit_points.emplace_back(0.0, 0.0, 0.5);
129 unit_points.emplace_back(0.5, 0.0, 0.5);
130 unit_points.emplace_back(0.0, 0.5, 0.5);
147 std::vector<std::vector<
Point<dim - 1>>>
148 unit_face_support_points_fe_poly(
const unsigned int degree)
154 std::vector<std::vector<
Point<dim - 1>>> unit_face_points;
162 unit_face_points.emplace_back(
163 unit_support_points_fe_poly<dim - 1>(degree));
166 return unit_face_points;
176 constraints_fe_poly(
const unsigned int )
185 constraints_fe_poly<2>(
const unsigned int degree)
187 const unsigned int dim = 2;
197 std::vector<
Point<dim - 1>> constraint_points;
199 constraint_points.emplace_back(0.5);
203 constraint_points.emplace_back(0.25);
205 constraint_points.emplace_back(0.75);
211 const unsigned int n_dofs_constrained = constraint_points.size();
212 unsigned int n_dofs_per_face = degree + 1;
218 for (
unsigned int i = 0; i < n_dofs_constrained; ++i)
219 for (
unsigned int j = 0; j < n_dofs_per_face; ++j)
221 interface_constraints(i, j) =
222 poly.compute_value(j, constraint_points[i]);
228 if (
std::fabs(interface_constraints(i, j)) < 1
e-13)
229 interface_constraints(i, j) = 0;
231 return interface_constraints;
238 std::vector<unsigned int>
239 get_dpo_vector_fe_dgp(
const unsigned int dim,
const unsigned int degree)
241 std::vector<unsigned int> dpo(dim + 1, 0
U);
244 if (dim == 2 && degree == 1)
246 else if (dim == 2 && degree == 2)
248 else if (dim == 3 && degree == 1)
250 else if (dim == 3 && degree == 2)
263template <
int dim,
int spacedim>
265 const unsigned int degree,
266 const std::vector<unsigned int> & dpo_vector,
297 unit_face_support_points_fe_poly<dim>(
degree);
303template <
int dim,
int spacedim>
304std::pair<Table<2, bool>, std::vector<unsigned int>>
308 constant_modes.fill(
true);
309 return std::pair<Table<2, bool>, std::vector<unsigned int>>(
310 constant_modes, std::vector<unsigned int>(1, 0));
315template <
int dim,
int spacedim>
318 const unsigned int child,
326 if (this->prolongation[refinement_case - 1][child].n() == 0)
328 std::lock_guard<std::mutex> lock(this->mutex);
331 if (this->prolongation[refinement_case - 1][child].n() ==
332 this->n_dofs_per_cell())
333 return this->prolongation[refinement_case - 1][child];
339 std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
341 isotropic_matrices.back().resize(
347 this_nonconst.prolongation[refinement_case - 1].swap(
348 isotropic_matrices.back());
352 return this->prolongation[refinement_case - 1][child];
357template <
int dim,
int spacedim>
360 const unsigned int child,
368 if (this->restriction[refinement_case - 1][child].n() == 0)
370 std::lock_guard<std::mutex> lock(this->mutex);
373 if (this->restriction[refinement_case - 1][child].n() ==
374 this->n_dofs_per_cell())
375 return this->restriction[refinement_case - 1][child];
381 std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
383 isotropic_matrices.back().resize(
389 this_nonconst.restriction[refinement_case - 1].swap(
390 isotropic_matrices.back());
394 return this->restriction[refinement_case - 1][child];
399template <
int dim,
int spacedim>
404 const unsigned int face_no)
const
418 const double eps = 2
e-13 * this->degree * (dim - 1);
420 std::vector<Point<dim>> face_quadrature_points(quad_face_support.size());
424 face_quadrature_points);
427 for (
unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
429 double matrix_entry =
430 this->shape_value(this->face_to_cell_index(j, 0),
431 face_quadrature_points[i]);
441 interpolation_matrix(i, j) = matrix_entry;
449 for (
unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
450 sum += interpolation_matrix(j, i);
464 spacedim>::ExcInterpolationNotImplemented()));
469template <
int dim,
int spacedim>
473 const unsigned int subface,
475 const unsigned int face_no)
const
489 const double eps = 2
e-13 * this->degree * (dim - 1);
491 std::vector<Point<dim>> subface_quadrature_points(
492 quad_face_support.size());
497 subface_quadrature_points);
500 for (
unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
502 double matrix_entry =
503 this->shape_value(this->face_to_cell_index(j, 0),
504 subface_quadrature_points[i]);
514 interpolation_matrix(i, j) = matrix_entry;
522 for (
unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
523 sum += interpolation_matrix(j, i);
537 spacedim>::ExcInterpolationNotImplemented()));
542template <
int dim,
int spacedim>
551template <
int dim,
int spacedim>
556 std::vector<double> & nodal_values)
const
559 this->get_unit_support_points().size());
563 for (
unsigned int i = 0; i < this->dofs_per_cell; ++i)
567 nodal_values[i] = support_point_values[i](0);
573template <
int dim,
int spacedim>
576 get_dpo_vector_fe_p(dim, degree),
582template <
int dim,
int spacedim>
583std::unique_ptr<FiniteElement<dim, spacedim>>
586 return std::make_unique<FE_SimplexP<dim, spacedim>>(*this);
591template <
int dim,
int spacedim>
595 std::ostringstream namebuf;
596 namebuf <<
"FE_SimplexP<" << dim <<
">(" << this->degree <<
")";
598 return namebuf.str();
603template <
int dim,
int spacedim>
607 const unsigned int codim)
const
628 if (this->degree < fe_p_other->degree)
630 else if (this->degree == fe_p_other->degree)
638 if (this->degree < fe_q_other->degree)
640 else if (this->degree == fe_q_other->degree)
648 if (fe_nothing->is_dominating())
663template <
int dim,
int spacedim>
664std::vector<std::pair<unsigned int, unsigned int>>
708template <
int dim,
int spacedim>
709std::vector<std::pair<unsigned int, unsigned int>>
728 std::vector<std::pair<unsigned int, unsigned int>> identities;
730 for (
unsigned int i = 0; i < this->degree - 1; ++i)
731 for (
unsigned int j = 0; j < fe_p_other->degree - 1; ++j)
733 fe_p_other->unit_support_points[i + 3][0]) < 1
e-14)
734 identities.emplace_back(i, j);
753 const std::vector<unsigned int> &index_map_inverse_q_other =
754 fe_q_other->get_poly_space_numbering_inverse();
756 std::vector<std::pair<unsigned int, unsigned int>> identities;
758 for (
unsigned int i = 0; i < this->degree - 1; ++i)
759 for (
unsigned int j = 0; j < fe_q_other->degree - 1; ++j)
761 fe_q_other->get_unit_support_points()
762 [index_map_inverse_q_other[j + 1]][0]) < 1
e-14)
763 identities.emplace_back(i, j);
793template <
int dim,
int spacedim>
796 get_dpo_vector_fe_dgp(dim, degree),
802template <
int dim,
int spacedim>
803std::unique_ptr<FiniteElement<dim, spacedim>>
806 return std::make_unique<FE_SimplexDGP<dim, spacedim>>(*this);
811template <
int dim,
int spacedim>
815 std::ostringstream namebuf;
816 namebuf <<
"FE_SimplexDGP<" << dim <<
">(" << this->degree <<
")";
818 return namebuf.str();
822template <
int dim,
int spacedim>
826 const unsigned int codim)
const
843 if (this->degree < fe_dgp_other->degree)
845 else if (this->degree == fe_dgp_other->degree)
853 if (this->degree < fe_dgq_other->degree)
855 else if (this->degree == fe_dgq_other->degree)
863 if (fe_nothing->is_dominating())
878template <
int dim,
int spacedim>
879std::vector<std::pair<unsigned int, unsigned int>>
890template <
int dim,
int spacedim>
891std::vector<std::pair<unsigned int, unsigned int>>
901#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
virtual const FullMatrix< double > & get_prolongation_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const override
FE_SimplexPoly(const unsigned int degree, const std::vector< unsigned int > &dpo_vector, const typename FiniteElementData< dim >::Conformity conformity)
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_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< Point< dim > > unit_support_points
FullMatrix< double > interface_constraints
static void project_to_subface(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 SubQuadrature &quadrature, const unsigned int face_no, std::vector< Point< dim > > &q_points)
SymmetricTensor< rank, dim, Number > sum(const SymmetricTensor< rank, dim, Number > &local, const MPI_Comm &mpi_communicator)
#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)
Expression fabs(const Expression &x)
@ either_element_can_dominate
@ other_element_dominates
@ neither_element_dominates
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
void L2(Vector< number > &result, const FEValuesBase< dim > &fe, const std::vector< double > &input, const double factor=1.)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
constexpr const ReferenceCell Tetrahedron
constexpr const ReferenceCell Triangle
std::vector< Point< dim > > unit_support_points(const std::vector< Point< 1 > > &line_support_points, const std::vector< unsigned int > &renumbering)