91 const unsigned int face_no = 0;
96 for (
unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_face; ++i)
99 FETools::compute_face_embedding_matrices<dim, double>(*
this,
106 unsigned int target_row = 0;
107 for (
unsigned int d = 0; d < GeometryInfo<dim>::max_children_per_face; ++d)
108 for (
unsigned int i = 0; i < face_embeddings[d].
m(); ++i)
110 for (
unsigned int j = 0; j < face_embeddings[d].
n(); ++j)
137 std::ostringstream namebuf;
138 namebuf <<
"FE_RaviartThomas<" << dim <<
">(" << this->degree - 1 <<
")";
140 return namebuf.str();
145std::unique_ptr<FiniteElement<dim, dim>>
148 return std::make_unique<FE_RaviartThomas<dim>>(*this);
162 const unsigned int n_interior_points = (deg > 0) ? cell_quadrature.
size() : 0;
167 const unsigned int face_no = 0;
169 unsigned int n_face_points = (dim > 1) ? 1 : 0;
171 for (
unsigned int d = 1; d < dim; ++d)
172 n_face_points *= deg + 1;
175 this->generalized_support_points.resize(
177 this->generalized_face_support_points[face_no].resize(n_face_points);
180 unsigned int current = 0;
184 QGauss<dim - 1> face_points(deg + 1);
188 boundary_weights.reinit(n_face_points, legendre.n());
190 for (
unsigned int k = 0; k < n_face_points; ++k)
192 this->generalized_face_support_points[face_no][k] =
193 face_points.point(k);
197 for (
unsigned int i = 0; i < legendre.n(); ++i)
199 boundary_weights(k, i) =
200 face_points.weight(k) *
201 legendre.compute_value(i, face_points.point(k));
208 for (; current < GeometryInfo<dim>::faces_per_cell * n_face_points;
213 this->generalized_support_points[current] = faces.
point(
216 this->reference_cell(), 0,
true,
false,
false, n_face_points));
224 std::unique_ptr<AnisotropicPolynomials<dim>> polynomials[dim];
225 for (
unsigned int dd = 0; dd < dim; ++dd)
227 std::vector<std::vector<Polynomials::Polynomial<double>>> poly(dim);
228 for (
unsigned int d = 0; d < dim; ++d)
232 polynomials[dd] = std::make_unique<AnisotropicPolynomials<dim>>(poly);
235 interior_weights.reinit(
238 for (
unsigned int k = 0; k < cell_quadrature.
size(); ++k)
240 this->generalized_support_points[current++] = cell_quadrature.
point(k);
241 for (
unsigned int i = 0; i < polynomials[0]->n(); ++i)
242 for (
unsigned int d = 0; d < dim; ++d)
243 interior_weights(k, i, d) =
244 cell_quadrature.
weight(k) *
245 polynomials[d]->compute_value(i, cell_quadrature.
point(k));
248 Assert(current == this->generalized_support_points.size(),
269 const unsigned int face_no = 0;
272 this->adjust_quad_dof_index_for_face_orientation_table[0].n_elements() ==
273 this->reference_cell().n_face_orientations(face_no) *
274 this->n_dofs_per_quad(face_no),
278 const unsigned int n = this->tensor_degree();
332 for (
unsigned int local_face_dof = 0;
333 local_face_dof < this->n_dofs_per_quad(face_no);
337 unsigned int i = local_face_dof % n;
338 unsigned int j = local_face_dof / n;
343 for (
const bool face_orientation : {
false,
true})
344 for (
const bool face_flip : {
false,
true})
345 for (
const bool face_rotation : {
false,
true})
352 if (((!face_orientation) && (!face_rotation)) ||
353 ((face_orientation) && (face_rotation)))
358 ->adjust_quad_dof_index_for_face_orientation_table[face_no](
359 local_face_dof, case_no) = j + i * n - local_face_dof;
365 ->adjust_quad_dof_index_for_face_orientation_table[face_no](
366 local_face_dof, case_no) = 0;
370 const unsigned int new_local_face_dof =
372 this->adjust_quad_dof_index_for_face_orientation_table[face_no](
373 local_face_dof, case_no);
375 i = new_local_face_dof % n;
376 j = new_local_face_dof / n;
383 if (!face_flip && face_rotation)
386 ->adjust_quad_dof_sign_for_face_orientation_table[face_no](
387 local_face_dof, case_no) = ((j % 2) == 1);
390 else if (face_flip && !face_rotation)
395 ->adjust_quad_dof_sign_for_face_orientation_table[face_no](
396 local_face_dof, case_no) =
397 ((j % 2) == 1) != ((i % 2) == 1);
400 else if (face_flip && face_rotation)
403 ->adjust_quad_dof_sign_for_face_orientation_table[face_no](
404 local_face_dof, case_no) = ((i % 2) == 1);
414 if (!face_orientation)
415 this->adjust_quad_dof_sign_for_face_orientation_table[face_no](
416 local_face_dof, case_no) =
418 ->adjust_quad_dof_sign_for_face_orientation_table[face_no](
419 local_face_dof, case_no);
432 for (
unsigned int i = 0; i < GeometryInfo<1>::max_children_per_cell; ++i)
433 this->restriction[0][i].reinit(0, 0);
454 QGauss<dim - 1> q_base(this->degree);
455 const unsigned int n_face_points = q_base.size();
471 for (
unsigned int k = 0; k < q_face.
size(); ++k)
472 for (
unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
473 cached_values_on_face(i, k) = this->shape_value_component(
476 for (
unsigned int sub = 0; sub < GeometryInfo<dim>::max_children_per_face;
484 this->reference_cell(), q_base, face, sub);
500 for (
unsigned int k = 0; k < n_face_points; ++k)
501 for (
unsigned int i_child = 0; i_child < this->n_dofs_per_cell();
503 for (
unsigned int i_face = 0;
504 i_face < this->n_dofs_per_face(face);
512 this->restriction[iso][child](
513 face * this->n_dofs_per_face(face) + i_face, i_child) +=
514 Utilities::fixed_power<dim - 1>(.5) * q_sub.
weight(k) *
515 cached_values_on_face(i_child, k) *
516 this->shape_value_component(
517 face * this->n_dofs_per_face(face) + i_face,
524 if (this->degree == 1)
529 std::unique_ptr<AnisotropicPolynomials<dim>> polynomials[dim];
530 for (
unsigned int dd = 0; dd < dim; ++dd)
532 std::vector<std::vector<Polynomials::Polynomial<double>>> poly(dim);
533 for (
unsigned int d = 0; d < dim; ++d)
539 polynomials[dd] = std::make_unique<AnisotropicPolynomials<dim>>(poly);
545 const unsigned int face_no = 0;
548 const unsigned int start_cell_dofs =
557 for (
unsigned int k = 0; k < q_cell.
size(); ++k)
558 for (
unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
559 for (
unsigned int d = 0; d < dim; ++d)
560 cached_values_on_cell(i, k, d) =
561 this->shape_value_component(i, q_cell.
point(k), d);
563 for (
unsigned int child = 0; child < GeometryInfo<dim>::max_children_per_cell;
571 for (
unsigned int k = 0; k < q_sub.
size(); ++k)
572 for (
unsigned int i_child = 0; i_child < this->n_dofs_per_cell();
574 for (
unsigned int d = 0; d < dim; ++d)
575 for (
unsigned int i_weight = 0; i_weight < polynomials[d]->n();
578 this->restriction[iso][child](start_cell_dofs + i_weight * dim +
581 q_sub.
weight(k) * cached_values_on_cell(i_child, k, d) *
582 polynomials[d]->compute_value(i_weight, q_sub.
point(k));
590std::vector<unsigned int>
595 unsigned int dofs_per_face = 1;
596 for (
unsigned int d = 1; d < dim; ++d)
597 dofs_per_face *= deg + 1;
600 const unsigned int interior_dofs = dim * deg * dofs_per_face;
602 std::vector<unsigned int> dpo(dim + 1);
603 dpo[dim - 1] = dofs_per_face;
604 dpo[dim] = interior_dofs;
612std::pair<Table<2, bool>, std::vector<unsigned int>>
616 for (
unsigned int d = 0; d < dim; ++d)
617 for (
unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
618 constant_modes(d, i) =
true;
619 std::vector<unsigned int> components;
620 for (
unsigned int d = 0; d < dim; ++d)
621 components.push_back(d);
622 return std::pair<Table<2, bool>, std::vector<unsigned int>>(constant_modes,
636 const unsigned int face_index)
const
644 switch (this->degree)
657 return (face_index !=
679 std::vector<double> &nodal_values)
const
681 Assert(support_point_values.size() == this->generalized_support_points.size(),
683 this->generalized_support_points.size()));
684 Assert(nodal_values.size() == this->n_dofs_per_cell(),
686 Assert(support_point_values[0].size() == this->n_components(),
688 this->n_components()));
690 std::fill(nodal_values.begin(), nodal_values.end(), 0.);
692 const unsigned int n_face_points = boundary_weights.size(0);
694 for (
unsigned int k = 0; k < n_face_points; ++k)
695 for (
unsigned int i = 0; i < boundary_weights.size(1); ++i)
697 nodal_values[i + face * this->n_dofs_per_face(face)] +=
698 boundary_weights(k, i) *
699 support_point_values[face * n_face_points + k](
706 const unsigned int face_no = 0;
708 const unsigned int start_cell_dofs =
710 const unsigned int start_cell_points =
713 for (
unsigned int k = 0; k < interior_weights.size(0); ++k)
714 for (
unsigned int i = 0; i < interior_weights.size(1); ++i)
715 for (
unsigned int d = 0; d < dim; ++d)
716 nodal_values[start_cell_dofs + i * dim + d] +=
717 interior_weights(k, i, d) *
718 support_point_values[k + start_cell_points](d);
734#include "fe_raviart_thomas.inst"
FullMatrix< double > inverse_node_matrix
std::vector< MappingKind > mapping_kind
virtual std::size_t memory_consumption() 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
void initialize_restriction()
virtual std::pair< Table< 2, bool >, std::vector< unsigned int > > get_constant_modes() const override
void initialize_quad_dof_index_permutation_and_sign_change()
virtual std::unique_ptr< FiniteElement< dim, dim > > clone() const override
friend class FE_RaviartThomas
void initialize_support_points(const unsigned int rt_degree)
static std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
virtual std::string get_name() const override
unsigned int n_dofs_per_cell() const
unsigned int n_dofs_per_face(unsigned int face_no=0, unsigned int child=0) const
unsigned int n_unique_faces() const
void reinit_restriction_and_prolongation_matrices(const bool isotropic_restriction_only=false, const bool isotropic_prolongation_only=false)
FullMatrix< double > interface_constraints
std::vector< std::vector< FullMatrix< double > > > prolongation
void invert(const FullMatrix< number2 > &M)
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
static Quadrature< dim > project_to_child(const ReferenceCell &reference_cell, const Quadrature< dim > &quadrature, const unsigned int child_no)
static Quadrature< dim > project_to_all_faces(const ReferenceCell &reference_cell, const hp::QCollection< dim - 1 > &quadrature)
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)
const Point< dim > & point(const unsigned int i) const
double weight(const unsigned int i) const
unsigned int size() const
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
#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)
#define DEAL_II_NOT_IMPLEMENTED()
unsigned char combined_face_orientation(const bool face_orientation, const bool face_rotation, const bool face_flip)
static unsigned int child_cell_on_face(const RefinementCase< dim > &ref_case, const unsigned int face, const unsigned int subface, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false, const RefinementCase< dim - 1 > &face_refinement_case=RefinementCase< dim - 1 >::isotropic_refinement)
static std_cxx20::ranges::iota_view< unsigned int, unsigned int > face_indices()