17 #include <deal.II/base/qprojector.h> 18 #include <deal.II/base/quadrature_lib.h> 19 #include <deal.II/base/std_cxx14/memory.h> 20 #include <deal.II/base/table.h> 22 #include <deal.II/dofs/dof_accessor.h> 24 #include <deal.II/fe/fe.h> 25 #include <deal.II/fe/fe_nothing.h> 26 #include <deal.II/fe/fe_raviart_thomas.h> 27 #include <deal.II/fe/fe_tools.h> 28 #include <deal.II/fe/fe_values.h> 29 #include <deal.II/fe/mapping.h> 31 #include <deal.II/grid/tria.h> 32 #include <deal.II/grid/tria_iterator.h> 37 DEAL_II_NAMESPACE_OPEN
51 std::vector<bool>(dim, true)))
78 ref_case < RefinementCase<dim>::isotropic_refinement + 1;
81 const unsigned int nc =
84 for (
unsigned int i = 0; i < nc; ++i)
85 this->
prolongation[ref_case - 1][i].reinit(n_dofs, n_dofs);
92 for (
unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_face; ++i)
94 FETools::compute_face_embedding_matrices<dim, double>(*
this,
100 unsigned int target_row = 0;
101 for (
unsigned int d = 0; d < GeometryInfo<dim>::max_children_per_face; ++d)
102 for (
unsigned int i = 0; i < face_embeddings[d].m(); ++i)
104 for (
unsigned int j = 0; j < face_embeddings[d].n(); ++j)
127 std::ostringstream namebuf;
128 namebuf <<
"FE_RaviartThomasNodal<" << dim <<
">(" << this->
degree - 1 <<
")";
130 return namebuf.str();
135 std::unique_ptr<FiniteElement<dim, dim>>
138 return std_cxx14::make_unique<FE_RaviartThomasNodal<dim>>(*this);
156 unsigned int current = 0;
166 QGauss<dim - 1> face_points(deg + 1);
172 for (
unsigned int k = 0;
177 0,
true,
false,
false, this->dofs_per_face));
179 current = this->dofs_per_face * GeometryInfo<dim>::faces_per_cell;
190 for (
unsigned int d = 0; d < dim; ++d)
192 std::unique_ptr<QAnisotropic<dim>> quadrature;
196 quadrature = std_cxx14::make_unique<QAnisotropic<dim>>(high);
199 quadrature = std_cxx14::make_unique<QAnisotropic<dim>>(
200 ((d == 0) ? low : high), ((d == 1) ? low : high));
204 std_cxx14::make_unique<QAnisotropic<dim>>(((d == 0) ? low : high),
205 ((d == 1) ? low : high),
213 for (
unsigned int k = 0; k < quadrature->size(); ++k)
222 std::vector<unsigned int>
228 for (
unsigned int d = 1; d < dim; ++d)
229 dofs_per_face *= deg + 1;
232 const unsigned int interior_dofs = dim * deg *
dofs_per_face;
234 std::vector<unsigned int> dpo(dim + 1);
236 dpo[dim] = interior_dofs;
248 return std::vector<bool>();
260 for (
unsigned int d = 2; d < dim; ++d)
261 dofs_per_face *= deg + 1;
267 std::vector<bool> ret_val(dofs_per_cell,
false);
280 const unsigned int shape_index,
281 const unsigned int face_index)
const 291 const unsigned int support_face = shape_index / this->
degree;
312 std::vector<double> & nodal_values)
const 327 unsigned int fbase = 0;
329 for (; f < GeometryInfo<dim>::faces_per_cell;
334 nodal_values[fbase + i] = support_point_values[fbase + i](
341 const unsigned int istep = (this->
dofs_per_cell - fbase) / dim;
347 for (
unsigned int i = 0; i < istep; ++i)
349 nodal_values[fbase + i] = support_point_values[fbase + i](f);
374 std::vector<std::pair<unsigned int, unsigned int>>
384 return std::vector<std::pair<unsigned int, unsigned int>>();
386 return std::vector<std::pair<unsigned int, unsigned int>>();
390 return std::vector<std::pair<unsigned int, unsigned int>>();
397 std::vector<std::pair<unsigned int, unsigned int>>
411 return std::vector<std::pair<unsigned int, unsigned int>>();
432 const unsigned int p = this->
degree - 1;
433 const unsigned int q = fe_q_other->degree - 1;
435 std::vector<std::pair<unsigned int, unsigned int>> identities;
438 for (
unsigned int i = 0; i < p + 1; ++i)
439 identities.emplace_back(i, i);
441 else if (p % 2 == 0 && q % 2 == 0)
442 identities.emplace_back(p / 2, q / 2);
450 return std::vector<std::pair<unsigned int, unsigned int>>();
455 return std::vector<std::pair<unsigned int, unsigned int>>();
461 std::vector<std::pair<unsigned int, unsigned int>>
475 return std::vector<std::pair<unsigned int, unsigned int>>();
480 const unsigned int q = fe_q_other->dofs_per_quad;
482 std::vector<std::pair<unsigned int, unsigned int>> identities;
485 for (
unsigned int i = 0; i < p; ++i)
486 identities.emplace_back(i, i);
488 else if (p % 2 != 0 && q % 2 != 0)
489 identities.emplace_back(p / 2, q / 2);
497 return std::vector<std::pair<unsigned int, unsigned int>>();
502 return std::vector<std::pair<unsigned int, unsigned int>>();
515 if (this->degree < fe_q_other->
degree)
517 else if (this->degree == fe_q_other->degree)
525 if (fe_q_other->is_dominating())
577 &x_source_fe) !=
nullptr),
619 double eps = 2e-13 * this->
degree * (dim - 1);
629 const Point<dim> &p = face_projection.point(i);
633 double matrix_entry =
642 if (std::fabs(matrix_entry - 1.0) < eps)
644 if (std::fabs(matrix_entry) < eps)
647 interpolation_matrix(i, j) = matrix_entry;
661 sum += interpolation_matrix(j, i);
663 Assert(std::fabs(sum - 1) < 2e-13 * this->
degree * (dim - 1),
673 const unsigned int subface,
681 &x_source_fe) !=
nullptr),
723 double eps = 2e-13 * this->
degree * (dim - 1);
734 const Point<dim> &p = subface_projection.point(i);
738 double matrix_entry =
747 if (std::fabs(matrix_entry - 1.0) < eps)
749 if (std::fabs(matrix_entry) < eps)
752 interpolation_matrix(i, j) = matrix_entry;
766 sum += interpolation_matrix(j, i);
768 Assert(std::fabs(sum - 1) < 2e-13 * this->
degree * (dim - 1),
776 #include "fe_raviart_thomas_nodal.inst" 779 DEAL_II_NAMESPACE_CLOSE
std::vector< Point< dim > > generalized_support_points
FullMatrix< double > interface_constraints
const std::vector< Point< dim-1 > > & get_generalized_face_support_points() const
FE_RaviartThomasNodal(const unsigned int p)
const unsigned int dofs_per_quad
const unsigned int degree
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
static unsigned int compute_n_pols(unsigned int degree)
#define AssertThrow(cond, exc)
std::vector< Point< dim-1 > > generalized_face_support_points
static::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
const Point< dim > & point(const unsigned int i) const
static std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
virtual bool hp_constraints_are_implemented() const override
static::ExceptionBase & ExcImpossibleInDim(int arg1)
std::vector< std::vector< FullMatrix< double > > > prolongation
void invert(const FullMatrix< number2 > &M)
void reinit(const TableIndices< N > &new_size, const bool omit_default_initialization=false)
#define Assert(cond, exc)
static::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
void initialize_support_points(const unsigned int rt_degree)
virtual double shape_value_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
virtual std::string get_name() const =0
unsigned int n_components() const
const unsigned int dofs_per_cell
virtual std::string get_name() const override
static unsigned int n_children(const RefinementCase< dim > &refinement_case)
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 Quadrature< dim > project_to_all_faces(const SubQuadrature &quadrature)
static void project_to_face(const SubQuadrature &quadrature, const unsigned int face_no, std::vector< Point< dim >> &q_points)
virtual unsigned int face_to_cell_index(const unsigned int face_dof_index, const unsigned int face, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false) const
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const override
const unsigned int dofs_per_face
static::ExceptionBase & ExcNotImplemented()
FullMatrix< double > inverse_node_matrix
static::ExceptionBase & ExcInternalError()
static std::vector< bool > get_ria_vector(const unsigned int degree)