17 #include <deal.II/base/qprojector.h> 18 #include <deal.II/base/quadrature.h> 19 #include <deal.II/base/quadrature_lib.h> 20 #include <deal.II/base/std_cxx14/memory.h> 21 #include <deal.II/base/table.h> 22 #include <deal.II/base/utilities.h> 24 #include <deal.II/dofs/dof_accessor.h> 26 #include <deal.II/fe/fe.h> 27 #include <deal.II/fe/fe_raviart_thomas.h> 28 #include <deal.II/fe/fe_tools.h> 29 #include <deal.II/fe/fe_values.h> 30 #include <deal.II/fe/mapping.h> 32 #include <deal.II/grid/tria.h> 33 #include <deal.II/grid/tria_iterator.h> 43 DEAL_II_NAMESPACE_OPEN
58 std::vector<bool>(dim, true)))
93 for (
unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_face; ++i)
95 FETools::compute_face_embedding_matrices<dim, double>(*
this,
101 unsigned int target_row = 0;
102 for (
unsigned int d = 0; d < GeometryInfo<dim>::max_children_per_face; ++d)
103 for (
unsigned int i = 0; i < face_embeddings[d].m(); ++i)
105 for (
unsigned int j = 0; j < face_embeddings[d].n(); ++j)
128 std::ostringstream namebuf;
129 namebuf <<
"FE_RaviartThomas<" << dim <<
">(" << this->
degree - 1 <<
")";
131 return namebuf.str();
136 std::unique_ptr<FiniteElement<dim, dim>>
139 return std_cxx14::make_unique<FE_RaviartThomas<dim>>(*this);
153 const unsigned int n_interior_points = (deg > 0) ? cell_quadrature.
size() : 0;
155 unsigned int n_face_points = (dim > 1) ? 1 : 0;
157 for (
unsigned int d = 1; d < dim; ++d)
158 n_face_points *= deg + 1;
166 unsigned int current = 0;
170 QGauss<dim - 1> face_points(deg + 1);
179 for (
unsigned int k = 0; k < n_face_points; ++k)
185 for (
unsigned int i = 0; i < legendre.n(); ++i)
188 face_points.weight(k) *
189 legendre.compute_value(i, face_points.point(k));
195 for (; current < GeometryInfo<dim>::faces_per_cell * n_face_points;
202 0,
true,
false,
false, n_face_points));
210 std::unique_ptr<AnisotropicPolynomials<dim>> polynomials[dim];
211 for (
unsigned int dd = 0; dd < dim; ++dd)
213 std::vector<std::vector<Polynomials::Polynomial<double>>> poly(dim);
214 for (
unsigned int d = 0; d < dim; ++d)
219 std_cxx14::make_unique<AnisotropicPolynomials<dim>>(poly);
225 for (
unsigned int k = 0; k < cell_quadrature.
size(); ++k)
228 for (
unsigned int i = 0; i < polynomials[0]->n(); ++i)
229 for (
unsigned int d = 0; d < dim; ++d)
231 cell_quadrature.
weight(k) *
232 polynomials[d]->compute_value(i, cell_quadrature.
point(k));
248 for (
unsigned int i = 0; i < GeometryInfo<1>::max_children_per_cell; ++i)
271 const unsigned int n_face_points = q_base.size();
274 for (
unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
286 for (
unsigned int k = 0; k < q_face.
size(); ++k)
291 for (
unsigned int sub = 0; sub < GeometryInfo<dim>::max_children_per_face;
315 for (
unsigned int k = 0; k < n_face_points; ++k)
316 for (
unsigned int i_child = 0; i_child < this->
dofs_per_cell;
326 this->
restriction[iso][child](face * this->dofs_per_face +
329 Utilities::fixed_power<dim - 1>(.5) * q_sub.
weight(k) *
330 cached_values_on_face(i_child, k) *
332 face * this->dofs_per_face + i_face,
344 std::unique_ptr<AnisotropicPolynomials<dim>> polynomials[dim];
345 for (
unsigned int dd = 0; dd < dim; ++dd)
347 std::vector<std::vector<Polynomials::Polynomial<double>>> poly(dim);
348 for (
unsigned int d = 0; d < dim; ++d)
355 std_cxx14::make_unique<AnisotropicPolynomials<dim>>(poly);
359 const unsigned int start_cell_dofs =
368 for (
unsigned int k = 0; k < q_cell.size(); ++k)
370 for (
unsigned int d = 0; d < dim; ++d)
371 cached_values_on_cell(i, k, d) =
374 for (
unsigned int child = 0; child < GeometryInfo<dim>::max_children_per_cell;
379 for (
unsigned int k = 0; k < q_sub.
size(); ++k)
380 for (
unsigned int i_child = 0; i_child < this->
dofs_per_cell; ++i_child)
381 for (
unsigned int d = 0; d < dim; ++d)
382 for (
unsigned int i_weight = 0; i_weight < polynomials[d]->n();
385 this->
restriction[iso][child](start_cell_dofs + i_weight * dim +
388 q_sub.
weight(k) * cached_values_on_cell(i_child, k, d) *
389 polynomials[d]->compute_value(i_weight, q_sub.
point(k));
397 std::vector<unsigned int>
403 for (
unsigned int d = 1; d < dim; ++d)
404 dofs_per_face *= deg + 1;
407 const unsigned int interior_dofs = dim * deg *
dofs_per_face;
409 std::vector<unsigned int> dpo(dim + 1);
411 dpo[dim] = interior_dofs;
419 std::pair<Table<2, bool>, std::vector<unsigned int>>
423 for (
unsigned int d = 0; d < dim; ++d)
425 constant_modes(d, i) =
true;
427 for (
unsigned int d = 0; d < dim; ++d)
428 components.push_back(d);
429 return std::pair<Table<2, bool>, std::vector<unsigned int>>(constant_modes,
443 const unsigned int face_index)
const 466 return (face_index !=
488 std::vector<double> & nodal_values)
const 499 std::fill(nodal_values.begin(), nodal_values.end(), 0.);
502 for (
unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
503 for (
unsigned int k = 0; k < n_face_points; ++k)
508 support_point_values[face * n_face_points + k](
512 const unsigned int start_cell_dofs =
514 const unsigned int start_cell_points =
519 for (
unsigned int d = 0; d < dim; ++d)
520 nodal_values[start_cell_dofs + i * dim + d] +=
522 support_point_values[k + start_cell_points](d);
538 #include "fe_raviart_thomas.inst" 541 DEAL_II_NAMESPACE_CLOSE
static Quadrature< dim > project_to_child(const Quadrature< dim > &quadrature, const unsigned int child_no)
std::vector< std::vector< FullMatrix< double > > > restriction
const unsigned int components
std::vector< Point< dim > > generalized_support_points
FullMatrix< double > interface_constraints
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
const unsigned int degree
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
virtual std::size_t memory_consumption() const override
virtual std::string get_name() const override
unsigned int size() const
static::ExceptionBase & ExcImpossibleInDim(int arg1)
friend class FE_RaviartThomas
std::vector< std::vector< FullMatrix< double > > > prolongation
void invert(const FullMatrix< number2 > &M)
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)
void reinit(const TableIndices< N > &new_size, const bool omit_default_initialization=false)
size_type size(const unsigned int i) const
#define Assert(cond, exc)
void initialize_support_points(const unsigned int rt_degree)
static::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
double weight(const unsigned int i) const
virtual double shape_value_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
unsigned int n_components() const
void initialize_restriction()
const unsigned int dofs_per_cell
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)
virtual std::unique_ptr< FiniteElement< dim, dim > > clone() const override
static void project_to_face(const SubQuadrature &quadrature, const unsigned int face_no, std::vector< Point< dim >> &q_points)
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
const unsigned int dofs_per_face
void reinit_restriction_and_prolongation_matrices(const bool isotropic_restriction_only=false, const bool isotropic_prolongation_only=false)
static::ExceptionBase & ExcNotImplemented()
static std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
virtual std::pair< Table< 2, bool >, std::vector< unsigned int > > get_constant_modes() const override
FullMatrix< double > inverse_node_matrix
Table< 3, double > interior_weights
Table< 2, double > boundary_weights
static::ExceptionBase & ExcInternalError()