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 (; current < this->reference_cell().n_faces() * n_face_points;
215 this->generalized_support_points[current] =
217 this->reference_cell(),
228 std::unique_ptr<AnisotropicPolynomials<dim>> polynomials[dim];
229 for (
unsigned int dd = 0; dd < dim; ++dd)
231 std::vector<std::vector<Polynomials::Polynomial<double>>> poly(dim);
232 for (
unsigned int d = 0; d < dim; ++d)
236 polynomials[dd] = std::make_unique<AnisotropicPolynomials<dim>>(poly);
239 interior_weights.reinit(
242 for (
unsigned int k = 0; k < cell_quadrature.size(); ++k)
244 this->generalized_support_points[current++] = cell_quadrature.point(k);
245 for (
unsigned int i = 0; i < polynomials[0]->n(); ++i)
246 for (
unsigned int d = 0; d < dim; ++d)
247 interior_weights(k, i, d) =
248 cell_quadrature.weight(k) *
249 polynomials[d]->compute_value(i, cell_quadrature.point(k));
252 Assert(current == this->generalized_support_points.size(),
273 const unsigned int face_no = 0;
276 this->adjust_quad_dof_index_for_face_orientation_table[0].n_elements() ==
277 this->reference_cell().n_face_orientations(face_no) *
278 this->n_dofs_per_quad(face_no),
282 const unsigned int n = this->tensor_degree();
336 for (
unsigned int local_face_dof = 0;
337 local_face_dof < this->n_dofs_per_quad(face_no);
341 unsigned int i = local_face_dof % n;
342 unsigned int j = local_face_dof / n;
347 for (
const bool face_orientation : {
false,
true})
348 for (
const bool face_flip : {
false,
true})
349 for (
const bool face_rotation : {
false,
true})
356 if (((!face_orientation) && (!face_rotation)) ||
357 ((face_orientation) && (face_rotation)))
362 ->adjust_quad_dof_index_for_face_orientation_table[face_no](
363 local_face_dof, case_no) = j + i * n - local_face_dof;
369 ->adjust_quad_dof_index_for_face_orientation_table[face_no](
370 local_face_dof, case_no) = 0;
374 const unsigned int new_local_face_dof =
376 this->adjust_quad_dof_index_for_face_orientation_table[face_no](
377 local_face_dof, case_no);
379 i = new_local_face_dof % n;
380 j = new_local_face_dof / n;
387 if (!face_flip && face_rotation)
390 ->adjust_quad_dof_sign_for_face_orientation_table[face_no](
391 local_face_dof, case_no) = ((j % 2) == 1);
394 else if (face_flip && !face_rotation)
399 ->adjust_quad_dof_sign_for_face_orientation_table[face_no](
400 local_face_dof, case_no) =
401 ((j % 2) == 1) != ((i % 2) == 1);
404 else if (face_flip && face_rotation)
407 ->adjust_quad_dof_sign_for_face_orientation_table[face_no](
408 local_face_dof, case_no) = ((i % 2) == 1);
418 if (!face_orientation)
419 this->adjust_quad_dof_sign_for_face_orientation_table[face_no](
420 local_face_dof, case_no) =
422 ->adjust_quad_dof_sign_for_face_orientation_table[face_no](
423 local_face_dof, case_no);
436 for (
auto &restriction_matrix : this->restriction[0])
437 restriction_matrix.reinit(0, 0);
458 const QGauss<dim - 1> q_base(this->degree);
459 const unsigned int n_face_points = q_base.size();
462 for (
const unsigned int face : this->reference_cell().face_indices())
469 this->reference_cell(),
478 for (
unsigned int k = 0; k < q_face.
size(); ++k)
479 for (
unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
480 cached_values_on_face(i, k) = this->shape_value_component(
483 for (
unsigned int sub = 0; sub < this->reference_cell()
484 .face_reference_cell(face)
485 .template n_children<dim - 1>();
493 this->reference_cell(),
514 for (
unsigned int k = 0; k < n_face_points; ++k)
515 for (
unsigned int i_child = 0; i_child < this->n_dofs_per_cell();
517 for (
unsigned int i_face = 0;
518 i_face < this->n_dofs_per_face(face);
526 this->restriction[iso][child](
527 face * this->n_dofs_per_face(face) + i_face, i_child) +=
528 Utilities::fixed_power<dim - 1>(.5) * q_sub.
weight(k) *
529 cached_values_on_face(i_child, k) *
530 this->shape_value_component(
531 face * this->n_dofs_per_face(face) + i_face,
538 if (this->degree == 1)
543 std::unique_ptr<AnisotropicPolynomials<dim>> polynomials[dim];
544 for (
unsigned int dd = 0; dd < dim; ++dd)
546 std::vector<std::vector<Polynomials::Polynomial<double>>> poly(dim);
547 for (
unsigned int d = 0; d < dim; ++d)
553 polynomials[dd] = std::make_unique<AnisotropicPolynomials<dim>>(poly);
559 const unsigned int face_no = 0;
562 const unsigned int start_cell_dofs =
563 this->reference_cell().n_faces() * this->n_dofs_per_face(face_no);
571 for (
unsigned int k = 0; k < q_cell.size(); ++k)
572 for (
unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
573 for (
unsigned int d = 0; d < dim; ++d)
574 cached_values_on_cell(i, k, d) =
575 this->shape_value_component(i, q_cell.point(k), d);
577 for (
unsigned int child = 0;
578 child < this->reference_cell().template n_children<dim>();
586 for (
unsigned int k = 0; k < q_sub.
size(); ++k)
587 for (
unsigned int i_child = 0; i_child < this->n_dofs_per_cell();
589 for (
unsigned int d = 0; d < dim; ++d)
590 for (
unsigned int i_weight = 0; i_weight < polynomials[d]->n();
593 this->restriction[iso][child](start_cell_dofs + i_weight * dim +
596 q_sub.
weight(k) * cached_values_on_cell(i_child, k, d) *
597 polynomials[d]->compute_value(i_weight, q_sub.
point(k));
605std::vector<unsigned int>
610 unsigned int dofs_per_face = 1;
611 for (
unsigned int d = 1; d < dim; ++d)
612 dofs_per_face *= deg + 1;
615 const unsigned int interior_dofs = dim * deg * dofs_per_face;
617 std::vector<unsigned int> dpo(dim + 1);
618 dpo[dim - 1] = dofs_per_face;
619 dpo[dim] = interior_dofs;
627std::pair<Table<2, bool>, std::vector<unsigned int>>
631 for (
unsigned int d = 0; d < dim; ++d)
632 for (
unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
633 constant_modes(d, i) =
true;
634 std::vector<unsigned int> components;
635 components.reserve(dim);
636 for (
unsigned int d = 0; d < dim; ++d)
637 components.push_back(d);
638 return std::pair<Table<2, bool>, std::vector<unsigned int>>(constant_modes,
652 const unsigned int face_index)
const
660 switch (this->degree)
673 return (face_index !=
695 std::vector<double> &nodal_values)
const
697 Assert(support_point_values.size() == this->generalized_support_points.size(),
699 this->generalized_support_points.size()));
700 Assert(nodal_values.size() == this->n_dofs_per_cell(),
702 Assert(support_point_values[0].
size() == this->n_components(),
704 this->n_components()));
706 std::fill(nodal_values.begin(), nodal_values.end(), 0.);
708 const unsigned int n_face_points = boundary_weights.size(0);
709 for (
const unsigned int face : this->reference_cell().face_indices())
710 for (
unsigned int k = 0; k < n_face_points; ++k)
711 for (
unsigned int i = 0; i < boundary_weights.size(1); ++i)
713 nodal_values[i + face * this->n_dofs_per_face(face)] +=
714 boundary_weights(k, i) *
715 support_point_values[face * n_face_points + k](
722 const unsigned int face_no = 0;
724 const unsigned int start_cell_dofs =
725 this->reference_cell().n_faces() * this->n_dofs_per_face(face_no);
726 const unsigned int start_cell_points =
727 this->reference_cell().n_faces() * n_face_points;
729 for (
unsigned int k = 0; k < interior_weights.size(0); ++k)
730 for (
unsigned int i = 0; i < interior_weights.size(1); ++i)
731 for (
unsigned int d = 0; d < dim; ++d)
732 nodal_values[start_cell_dofs + i * dim + d] +=
733 interior_weights(k, i, d) *
734 support_point_values[k + start_cell_points](d);
750#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)
Class which transforms dim - 1-dimensional quadrature rules to dim-dimensional face quadratures.
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 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)