46 std::vector<unsigned int>
47 get_rt_dpo_vector(
const unsigned int dim,
const unsigned int degree)
49 std::vector<unsigned int> dpo(dim + 1);
52 unsigned int dofs_per_face = 1;
53 for (
unsigned int d = 1;
d < dim; ++
d)
54 dofs_per_face *= (degree + 1);
56 dpo[dim - 1] = dofs_per_face;
57 dpo[dim] = dim * degree * dofs_per_face;
87 const std::vector<unsigned int> numbering =
96 this->n_dofs_per_cell());
98 const unsigned int face_no = 0;
105 for (
unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_face; ++i)
108 FETools::compute_face_embedding_matrices<dim, double>(*
this,
115 unsigned int target_row = 0;
116 for (
unsigned int d = 0; d < GeometryInfo<dim>::max_children_per_face; ++d)
117 for (
unsigned int i = 0; i < face_embeddings[d].
m(); ++i)
119 for (
unsigned int j = 0; j < face_embeddings[d].
n(); ++j)
141 return "FE_RaviartThomasNodal<" + std::to_string(dim) +
">(" +
142 std::to_string(this->degree - 1) +
")";
147std::unique_ptr<FiniteElement<dim, dim>>
150 return std::make_unique<FE_RaviartThomasNodal<dim>>(*this);
163 dim>::initialize_quad_dof_index_permutation_and_sign_change()
169 const unsigned int n = this->degree;
170 const unsigned int face_no = 0;
172 for (
unsigned int local = 0; local < this->n_dofs_per_quad(face_no); ++local)
176 unsigned int i = local % n, j = local / n;
179 this->adjust_quad_dof_index_for_face_orientation_table[face_no](
183 this->adjust_quad_dof_index_for_face_orientation_table[face_no](
185 i + (n - 1 - j) * n - local;
187 this->adjust_quad_dof_index_for_face_orientation_table[face_no](
189 (n - 1 - j) + (n - 1 - i) * n - local;
191 this->adjust_quad_dof_index_for_face_orientation_table[face_no](
193 (n - 1 - i) + j * n - local;
195 this->adjust_quad_dof_index_for_face_orientation_table[face_no](
198 this->adjust_quad_dof_index_for_face_orientation_table[face_no](
200 j + (n - 1 - i) * n - local;
202 this->adjust_quad_dof_index_for_face_orientation_table[face_no](
204 (n - 1 - i) + (n - 1 - j) * n - local;
206 this->adjust_quad_dof_index_for_face_orientation_table[face_no](
208 (n - 1 - j) + i * n - local;
211 for (
const bool rotation : {
false,
true})
212 for (
const bool flip : {
false,
true})
213 this->adjust_quad_dof_sign_for_face_orientation_table[face_no](
224 const unsigned int shape_index,
225 const unsigned int face_index)
const
232 const unsigned int support_face = shape_index / this->n_dofs_per_face();
250 std::vector<double> &nodal_values)
const
252 Assert(support_point_values.size() == this->generalized_support_points.size(),
254 this->generalized_support_points.size()));
255 Assert(nodal_values.size() == this->n_dofs_per_cell(),
257 Assert(support_point_values[0].
size() == this->n_components(),
259 this->n_components()));
263 unsigned int fbase = 0;
265 for (; f < GeometryInfo<dim>::faces_per_cell;
266 ++f, fbase += this->n_dofs_per_face(f))
268 for (
unsigned int i = 0; i < this->n_dofs_per_face(f); ++i)
270 nodal_values[fbase + i] = support_point_values[fbase + i](
276 const unsigned int istep = (this->n_dofs_per_cell() - fbase) / dim;
280 while (fbase < this->n_dofs_per_cell())
282 for (
unsigned int i = 0; i < istep; ++i)
284 nodal_values[fbase + i] = support_point_values[fbase + i](f);
309std::vector<std::pair<unsigned int, unsigned int>>
317 return std::vector<std::pair<unsigned int, unsigned int>>();
319 return std::vector<std::pair<unsigned int, unsigned int>>();
323 return std::vector<std::pair<unsigned int, unsigned int>>();
330std::vector<std::pair<unsigned int, unsigned int>>
341 return std::vector<std::pair<unsigned int, unsigned int>>();
356 const unsigned int p = this->degree - 1;
357 const unsigned int q = fe_q_other->degree - 1;
359 std::vector<std::pair<unsigned int, unsigned int>> identities;
362 for (
unsigned int i = 0; i < p + 1; ++i)
363 identities.emplace_back(i, i);
365 else if (p % 2 == 0 && q % 2 == 0)
366 identities.emplace_back(p / 2, q / 2);
374 return std::vector<std::pair<unsigned int, unsigned int>>();
379 return std::vector<std::pair<unsigned int, unsigned int>>();
385std::vector<std::pair<unsigned int, unsigned int>>
388 const unsigned int face_no)
const
397 return std::vector<std::pair<unsigned int, unsigned int>>();
400 const unsigned int p = this->n_dofs_per_quad(face_no);
403 const unsigned int q = fe_q_other->n_dofs_per_quad(0);
405 std::vector<std::pair<unsigned int, unsigned int>> identities;
408 for (
unsigned int i = 0; i < p; ++i)
409 identities.emplace_back(i, i);
411 else if (p % 2 != 0 && q % 2 != 0)
412 identities.emplace_back(p / 2, q / 2);
420 return std::vector<std::pair<unsigned int, unsigned int>>();
425 return std::vector<std::pair<unsigned int, unsigned int>>();
434 const unsigned int codim)
const
444 if (this->degree < fe_rt_nodal_other->degree)
446 else if (this->degree == fe_rt_nodal_other->degree)
454 if (fe_nothing->is_dominating())
474 const unsigned int)
const
486 const unsigned int)
const
498 const unsigned int face_no)
const
505 &x_source_fe) !=
nullptr),
508 Assert(interpolation_matrix.
n() == this->n_dofs_per_face(face_no),
510 this->n_dofs_per_face(face_no)));
536 double eps = 2e-13 * this->degree * (dim - 1);
548 const Point<dim> &p = face_projection.point(i);
550 for (
unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
552 double matrix_entry =
553 this->shape_value_component(this->face_to_cell_index(j, 0), p, 0);
558 if (std::fabs(matrix_entry - 1.0) < eps)
560 if (std::fabs(matrix_entry) < eps)
563 interpolation_matrix(i, j) = matrix_entry;
575 for (
unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
576 sum += interpolation_matrix(j, i);
578 Assert(std::fabs(sum - 1) < 2e-13 * this->degree * (dim - 1),
589 const unsigned int subface,
591 const unsigned int face_no)
const
597 &x_source_fe) !=
nullptr),
600 Assert(interpolation_matrix.
n() == this->n_dofs_per_face(face_no),
602 this->n_dofs_per_face(face_no)));
628 double eps = 2e-13 * this->degree * (dim - 1);
634 this->reference_cell(),
643 const Point<dim> &p = subface_projection.point(i);
645 for (
unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
647 double matrix_entry =
648 this->shape_value_component(this->face_to_cell_index(j, 0), p, 0);
653 if (std::fabs(matrix_entry - 1.0) < eps)
655 if (std::fabs(matrix_entry) < eps)
658 interpolation_matrix(i, j) = matrix_entry;
670 for (
unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
671 sum += interpolation_matrix(j, i);
673 Assert(std::fabs(sum - 1) < 2e-13 * this->degree * (dim - 1),
684 const unsigned int child,
691 "Prolongation matrices are only available for refined cells!"));
693 child, this->reference_cell().
template n_children<dim>(refinement_case));
696 if (this->prolongation[refinement_case - 1][child].n() == 0)
698 std::lock_guard<std::mutex> lock(prolongation_matrix_mutex);
701 if (this->prolongation[refinement_case - 1][child].n() ==
702 this->n_dofs_per_cell())
703 return this->prolongation[refinement_case - 1][child];
711 std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
713 isotropic_matrices.back().resize(
714 this->reference_cell().
template n_children<dim>(
717 this->n_dofs_per_cell()));
720 std::move(isotropic_matrices.back());
736 return this->prolongation[refinement_case - 1][child];
744 const unsigned int child,
751 "Restriction matrices are only available for refined cells!"));
753 child, this->reference_cell().
template n_children<dim>(refinement_case));
756 if (this->restriction[refinement_case - 1][child].n() == 0)
758 std::lock_guard<std::mutex> lock(restriction_matrix_mutex);
761 if (this->restriction[refinement_case - 1][child].n() ==
762 this->n_dofs_per_cell())
763 return this->restriction[refinement_case - 1][child];
771 std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
773 isotropic_matrices.back().resize(
774 this->reference_cell().
template n_children<dim>(
777 this->n_dofs_per_cell()));
780 std::move(isotropic_matrices.back());
796 return this->restriction[refinement_case - 1][child];
802#include "fe/fe_raviart_thomas_nodal.inst"
std::vector< MappingKind > mapping_kind
virtual std::unique_ptr< FiniteElement< dim, dim > > clone() 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
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const override
virtual const FullMatrix< double > & get_restriction_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const override
void initialize_quad_dof_index_permutation_and_sign_change()
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_line_dof_identities(const FiniteElement< dim > &fe_other) const override
virtual void get_face_interpolation_matrix(const FiniteElement< dim > &source, FullMatrix< double > &matrix, const unsigned int face_no=0) const override
virtual std::string get_name() const override
virtual bool hp_constraints_are_implemented() const override
virtual void get_subface_interpolation_matrix(const FiniteElement< dim > &source, const unsigned int subface, FullMatrix< double > &matrix, const unsigned int face_no=0) const override
virtual FiniteElementDomination::Domination compare_for_domination(const FiniteElement< dim > &fe_other, const unsigned int codim=0) const override final
virtual const FullMatrix< double > & get_prolongation_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const override
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_vertex_dof_identities(const FiniteElement< dim > &fe_other) const override
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_quad_dof_identities(const FiniteElement< dim > &fe_other, const unsigned int face_no=0) const override
FE_RaviartThomasNodal(const unsigned int p)
static std::vector< unsigned int > get_lexicographic_numbering(const unsigned int degree)
const unsigned int degree
unsigned int n_dofs_per_face(unsigned int face_no=0, unsigned int child=0) const
virtual std::string get_name() const =0
std::vector< std::vector< FullMatrix< double > > > restriction
void reinit_restriction_and_prolongation_matrices(const bool isotropic_restriction_only=false, const bool isotropic_prolongation_only=false)
std::vector< std::vector< Point< dim - 1 > > > generalized_face_support_points
FullMatrix< double > interface_constraints
std::vector< Point< dim > > generalized_support_points
std::vector< std::vector< FullMatrix< double > > > prolongation
std::vector< Point< dim > > get_polynomial_support_points() const
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)
#define DEAL_II_NAMESPACE_OPEN
constexpr bool running_in_debug_mode()
#define DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NOT_IMPLEMENTED()
#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
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
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