91 const unsigned int face_no = 0;
95 std::vector<FullMatrix<double>> face_embeddings(
96 this->
reference_cell().face_reference_cell(0).
template n_children<dim - 1>(
98 for (
auto &face_embedding : face_embeddings)
99 face_embedding.reinit(this->n_dofs_per_face(face_no),
108 unsigned int target_row = 0;
109 for (
const auto &face_embedding : face_embeddings)
110 for (
unsigned int i = 0; i < face_embedding.m(); ++i)
112 for (
unsigned int j = 0; j < face_embedding.n(); ++j)
139 std::ostringstream namebuf;
140 namebuf <<
"FE_RaviartThomas<" << dim <<
">(" << this->degree - 1 <<
")";
142 return namebuf.str();
147std::unique_ptr<FiniteElement<dim, dim>>
150 return std::make_unique<FE_RaviartThomas<dim>>(*this);
164 const unsigned int n_interior_points = (deg > 0) ? cell_quadrature.size() : 0;
169 const unsigned int face_no = 0;
171 unsigned int n_face_points = (dim > 1) ? 1 : 0;
173 for (
unsigned int d = 1; d < dim; ++d)
174 n_face_points *= deg + 1;
177 this->generalized_support_points.resize(
178 this->reference_cell().n_faces() * n_face_points + n_interior_points);
179 this->generalized_face_support_points[face_no].resize(n_face_points);
182 unsigned int current = 0;
186 const QGauss<dim - 1> face_points(deg + 1);
190 boundary_weights.reinit(n_face_points, legendre.n());
192 for (
unsigned int k = 0; k < n_face_points; ++k)
194 this->generalized_face_support_points[face_no][k] =
195 face_points.point(k);
199 for (
unsigned int i = 0; i < legendre.n(); ++i)
201 boundary_weights(k, i) =
202 face_points.weight(k) *
203 legendre.compute_value(i, face_points.point(k));
211 for (
unsigned int face_no = 0;
212 face_no < GeometryInfo<dim>::faces_per_cell;
216 this->reference_cell(),
220 for (
unsigned int face_point = 0; face_point < n_face_points;
224 this->generalized_support_points[current] =
225 faces.
point(offset + face_point);
235 std::unique_ptr<AnisotropicPolynomials<dim>> polynomials[dim];
236 for (
unsigned int dd = 0; dd < dim; ++dd)
238 std::vector<std::vector<Polynomials::Polynomial<double>>> poly(dim);
239 for (
unsigned int d = 0; d < dim; ++d)
243 polynomials[dd] = std::make_unique<AnisotropicPolynomials<dim>>(poly);
246 interior_weights.reinit(
249 for (
unsigned int k = 0; k < cell_quadrature.size(); ++k)
251 this->generalized_support_points[current++] = cell_quadrature.point(k);
252 for (
unsigned int i = 0; i < polynomials[0]->n(); ++i)
253 for (
unsigned int d = 0; d < dim; ++d)
254 interior_weights(k, i, d) =
255 cell_quadrature.weight(k) *
256 polynomials[d]->compute_value(i, cell_quadrature.point(k));
259 Assert(current == this->generalized_support_points.size(),
280 const unsigned int face_no = 0;
283 this->adjust_quad_dof_index_for_face_orientation_table[0].n_elements() ==
284 this->reference_cell().n_face_orientations(face_no) *
285 this->n_dofs_per_quad(face_no),
289 const unsigned int n = this->tensor_degree();
343 for (
unsigned int local_face_dof = 0;
344 local_face_dof < this->n_dofs_per_quad(face_no);
348 unsigned int i = local_face_dof % n;
349 unsigned int j = local_face_dof / n;
354 for (
const bool face_orientation : {
false,
true})
355 for (
const bool face_flip : {
false,
true})
356 for (
const bool face_rotation : {
false,
true})
363 if (((!face_orientation) && (!face_rotation)) ||
364 ((face_orientation) && (face_rotation)))
369 ->adjust_quad_dof_index_for_face_orientation_table[face_no](
370 local_face_dof, case_no) = j + i * n - local_face_dof;
376 ->adjust_quad_dof_index_for_face_orientation_table[face_no](
377 local_face_dof, case_no) = 0;
381 const unsigned int new_local_face_dof =
383 this->adjust_quad_dof_index_for_face_orientation_table[face_no](
384 local_face_dof, case_no);
386 i = new_local_face_dof % n;
387 j = new_local_face_dof / n;
394 if (!face_flip && face_rotation)
397 ->adjust_quad_dof_sign_for_face_orientation_table[face_no](
398 local_face_dof, case_no) = ((j % 2) == 1);
401 else if (face_flip && !face_rotation)
406 ->adjust_quad_dof_sign_for_face_orientation_table[face_no](
407 local_face_dof, case_no) =
408 ((j % 2) == 1) != ((i % 2) == 1);
411 else if (face_flip && face_rotation)
414 ->adjust_quad_dof_sign_for_face_orientation_table[face_no](
415 local_face_dof, case_no) = ((i % 2) == 1);
425 if (!face_orientation)
426 this->adjust_quad_dof_sign_for_face_orientation_table[face_no](
427 local_face_dof, case_no) =
429 ->adjust_quad_dof_sign_for_face_orientation_table[face_no](
430 local_face_dof, case_no);
443 for (
auto &restriction_matrix : this->restriction[0])
444 restriction_matrix.reinit(0, 0);
465 const QGauss<dim - 1> q_base(this->degree);
466 const unsigned int n_face_points = q_base.size();
469 for (
const unsigned int face : this->reference_cell().face_indices())
476 this->reference_cell(),
485 for (
unsigned int k = 0; k < q_face.
size(); ++k)
486 for (
unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
487 cached_values_on_face(i, k) = this->shape_value_component(
490 for (
unsigned int sub = 0; sub < this->reference_cell()
491 .face_reference_cell(face)
492 .template n_children<dim - 1>();
500 this->reference_cell(),
521 for (
unsigned int k = 0; k < n_face_points; ++k)
522 for (
unsigned int i_child = 0; i_child < this->n_dofs_per_cell();
524 for (
unsigned int i_face = 0;
525 i_face < this->n_dofs_per_face(face);
533 this->restriction[iso][child](
534 face * this->n_dofs_per_face(face) + i_face, i_child) +=
535 Utilities::fixed_power<dim - 1>(.5) * q_sub.
weight(k) *
536 cached_values_on_face(i_child, k) *
537 this->shape_value_component(
538 face * this->n_dofs_per_face(face) + i_face,
545 if (this->degree == 1)
550 std::unique_ptr<AnisotropicPolynomials<dim>> polynomials[dim];
551 for (
unsigned int dd = 0; dd < dim; ++dd)
553 std::vector<std::vector<Polynomials::Polynomial<double>>> poly(dim);
554 for (
unsigned int d = 0; d < dim; ++d)
560 polynomials[dd] = std::make_unique<AnisotropicPolynomials<dim>>(poly);
566 const unsigned int face_no = 0;
569 const unsigned int start_cell_dofs =
570 this->reference_cell().n_faces() * this->n_dofs_per_face(face_no);
578 for (
unsigned int k = 0; k < q_cell.size(); ++k)
579 for (
unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
580 for (
unsigned int d = 0; d < dim; ++d)
581 cached_values_on_cell(i, k, d) =
582 this->shape_value_component(i, q_cell.point(k), d);
584 for (
unsigned int child = 0;
585 child < this->reference_cell().template n_children<dim>();
593 for (
unsigned int k = 0; k < q_sub.
size(); ++k)
594 for (
unsigned int i_child = 0; i_child < this->n_dofs_per_cell();
596 for (
unsigned int d = 0; d < dim; ++d)
597 for (
unsigned int i_weight = 0; i_weight < polynomials[d]->n();
600 this->restriction[iso][child](start_cell_dofs + i_weight * dim +
603 q_sub.
weight(k) * cached_values_on_cell(i_child, k, d) *
604 polynomials[d]->compute_value(i_weight, q_sub.
point(k));
612std::vector<unsigned int>
617 unsigned int dofs_per_face = 1;
618 for (
unsigned int d = 1; d < dim; ++d)
619 dofs_per_face *= deg + 1;
622 const unsigned int interior_dofs = dim * deg * dofs_per_face;
624 std::vector<unsigned int> dpo(dim + 1);
625 dpo[dim - 1] = dofs_per_face;
626 dpo[dim] = interior_dofs;
634std::pair<Table<2, bool>, std::vector<unsigned int>>
638 for (
unsigned int d = 0; d < dim; ++d)
639 for (
unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
640 constant_modes(d, i) =
true;
641 std::vector<unsigned int> components;
642 components.reserve(dim);
643 for (
unsigned int d = 0; d < dim; ++d)
644 components.push_back(d);
645 return std::pair<Table<2, bool>, std::vector<unsigned int>>(constant_modes,
659 const unsigned int face_index)
const
667 switch (this->degree)
680 return (face_index !=
702 std::vector<double> &nodal_values)
const
704 Assert(support_point_values.size() == this->generalized_support_points.size(),
706 this->generalized_support_points.size()));
707 Assert(nodal_values.size() == this->n_dofs_per_cell(),
709 Assert(support_point_values[0].
size() == this->n_components(),
711 this->n_components()));
713 std::fill(nodal_values.begin(), nodal_values.end(), 0.);
715 const unsigned int n_face_points = boundary_weights.size(0);
716 for (
const unsigned int face : this->reference_cell().face_indices())
717 for (
unsigned int k = 0; k < n_face_points; ++k)
718 for (
unsigned int i = 0; i < boundary_weights.size(1); ++i)
720 nodal_values[i + face * this->n_dofs_per_face(face)] +=
721 boundary_weights(k, i) *
722 support_point_values[face * n_face_points + k](
729 const unsigned int face_no = 0;
731 const unsigned int start_cell_dofs =
732 this->reference_cell().n_faces() * this->n_dofs_per_face(face_no);
733 const unsigned int start_cell_points =
734 this->reference_cell().n_faces() * n_face_points;
736 for (
unsigned int k = 0; k < interior_weights.size(0); ++k)
737 for (
unsigned int i = 0; i < interior_weights.size(1); ++i)
738 for (
unsigned int d = 0; d < dim; ++d)
739 nodal_values[start_cell_dofs + i * dim + d] +=
740 interior_weights(k, i, d) *
741 support_point_values[k + start_cell_points](d);
757#include "fe/fe_raviart_thomas.inst"
ArrayView< std::remove_reference_t< typename std::iterator_traits< Iterator >::reference >, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
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
ReferenceCell reference_cell() 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 DataSetDescriptor face(const ReferenceCell &reference_cell, const unsigned int face_no, const bool face_orientation, const bool face_flip, const bool face_rotation, const unsigned int n_quadrature_points)
static Quadrature< dim > project_to_all_faces(const ReferenceCell &reference_cell, const hp::QCollection< dim - 1 > &quadrature)
static Quadrature< dim > project_to_child(const ReferenceCell &reference_cell, const Quadrature< dim > &quadrature, const unsigned int child_no)
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 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)
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
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)