Loading [MathJax]/extensions/TeX/newcommand.js
 deal.II version GIT relicensing-3087-ga35b476a78 2025-04-19 20:30:01+00:00
\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}} \newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=} \newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]} \newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
fe_rt_bubbles.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2018 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15
18
20
21#include <deal.II/fe/fe.h>
23#include <deal.II/fe/fe_tools.h>
25#include <deal.II/fe/mapping.h>
26
27#include <deal.II/grid/tria.h>
29
30#include <memory>
31#include <sstream>
32
33
35
36// TODO: implement the adjust_quad_dof_index_for_face_orientation_table and
37// adjust_line_dof_index_for_line_orientation_table fields, and write tests
38// similar to bits/face_orientation_and_fe_q_*
39
40template <int dim>
41FE_RT_Bubbles<dim>::FE_RT_Bubbles(const unsigned int deg)
42 : FE_PolyTensor<dim>(
43 PolynomialsRT_Bubbles<dim>(deg),
44 FiniteElementData<dim>(get_dpo_vector(deg),
45 dim,
46 deg + 1,
47 FiniteElementData<dim>::Hdiv),
48 get_ria_vector(deg),
49 std::vector<ComponentMask>(PolynomialsRT_Bubbles<dim>::n_polynomials(deg),
50 ComponentMask(std::vector<bool>(dim, true))))
51{
52 Assert(dim >= 2, ExcImpossibleInDim(dim));
53 Assert(
54 deg >= 1,
56 "Lowest order RT_Bubbles element is degree 1, but you requested for degree 0"));
57 const unsigned int n_dofs = this->n_dofs_per_cell();
58
60 // Initialize support points and quadrature weights
62 // Compute the inverse node matrix to get
63 // the correct basis functions
65 this->inverse_node_matrix.reinit(n_dofs, n_dofs);
67
68 // Reinit the vectors of prolongation matrices to the
69 // right sizes. There are no restriction matrices implemented
70 for (const unsigned int ref_case :
73 {
74 const unsigned int nc = this->reference_cell().template n_children<dim>(
75 RefinementCase<dim>(ref_case));
76
77 for (unsigned int i = 0; i < nc; ++i)
78 this->prolongation[ref_case - 1][i].reinit(n_dofs, n_dofs);
79 }
80
81 // TODO: the implementation makes the assumption that all faces have the
82 // same number of dofs
84 const unsigned int face_no = 0;
85
86 // Fill prolongation matrices with embedding operators
87 // set tolerance to 1, as embedding error accumulate quickly
90 for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_face; ++i)
91 face_embeddings[i].reinit(this->n_dofs_per_face(face_no),
92 this->n_dofs_per_face(face_no));
93 FETools::compute_face_embedding_matrices<dim, double>(*this,
94 face_embeddings,
95 0,
96 0);
97 this->interface_constraints.reinit((1 << (dim - 1)) *
98 this->n_dofs_per_face(face_no),
99 this->n_dofs_per_face(face_no));
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)
103 {
104 for (unsigned int j = 0; j < face_embeddings[d].n(); ++j)
105 this->interface_constraints(target_row, j) = face_embeddings[d](i, j);
106 ++target_row;
107 }
108
109 // We need to initialize the dof permutation table and the one for the sign
110 // change.
112}
113
114
115template <int dim>
116void
118{
119 // for 1d and 2d, do nothing
120 if (dim < 3)
121 return;
122
123 // TODO: Implement this for this class
124 return;
125}
126
127
128
129template <int dim>
130std::string
132{
133 // Note: this->degree is the maximal polynomial degree and is thus one higher
134 // than the argument given to the constructor
135 std::ostringstream namebuf;
136 namebuf << "FE_RT_Bubbles<" << dim << ">(" << this->degree << ")";
137
138 return namebuf.str();
139}
140
141
142
143template <int dim>
144std::unique_ptr<FiniteElement<dim, dim>>
146{
147 return std::make_unique<FE_RT_Bubbles<dim>>(*this);
148}
149
150
151//---------------------------------------------------------------------------
152// Auxiliary and internal functions
153//---------------------------------------------------------------------------
154
155
156
157template <int dim>
158void
160{
161 // TODO: the implementation makes the assumption that all faces have the
162 // same number of dofs
163 AssertDimension(this->n_unique_faces(), 1);
164 const unsigned int face_no = 0;
165
166 this->generalized_support_points.resize(this->n_dofs_per_cell());
167 this->generalized_face_support_points[face_no].resize(
168 this->n_dofs_per_face(face_no));
169
170 // Index of the point being entered
171 unsigned int current = 0;
172
173 // On the faces, we choose as many Gauss-Lobatto points
174 // as required to determine the normal component uniquely.
175 // This is the deg of the RT_Bubble element plus one.
176 if (dim > 1)
177 {
178 const QGaussLobatto<dim - 1> face_points(deg + 1);
179 Assert(face_points.size() == this->n_dofs_per_face(face_no),
181 for (unsigned int k = 0; k < this->n_dofs_per_face(face_no); ++k)
182 this->generalized_face_support_points[face_no][k] =
183 face_points.point(k);
184 Quadrature<dim> faces =
185 QProjector<dim>::project_to_all_faces(this->reference_cell(),
186 face_points);
187 for (unsigned int face_no = 0;
188 face_no < GeometryInfo<dim>::faces_per_cell;
189 ++face_no)
190 {
192 this->reference_cell(),
193 face_no,
195 face_points.size());
196 for (unsigned int face_point = 0; face_point < face_points.size();
197 ++face_point)
198 {
199 // Enter the support point into the vector
200 this->generalized_support_points[current] =
201 faces.point(offset + face_point);
202 ++current;
203 }
204 }
205 }
206
207 if (deg == 1)
208 return;
209
210 // In the interior, we need anisotropic Gauss-Lobatto quadratures,
211 // one for each direction
212 const QGaussLobatto<1> high(deg + 1);
213 std::vector<Point<1>> pts = high.get_points();
214 if (pts.size() > 2)
215 {
216 pts.erase(pts.begin());
217 pts.erase(pts.end() - 1);
218 }
219
220 std::vector<double> wts(pts.size(), 1);
221 const Quadrature<1> low(pts, wts);
222
223 for (unsigned int d = 0; d < dim; ++d)
224 {
225 std::unique_ptr<QAnisotropic<dim>> quadrature;
226 switch (dim)
227 {
228 case 1:
229 quadrature = std::make_unique<QAnisotropic<dim>>(high);
230 break;
231 case 2:
232 quadrature =
233 std::make_unique<QAnisotropic<dim>>(((d == 0) ? low : high),
234 ((d == 1) ? low : high));
235 break;
236 case 3:
237 quadrature =
238 std::make_unique<QAnisotropic<dim>>(((d == 0) ? low : high),
239 ((d == 1) ? low : high),
240 ((d == 2) ? low : high));
241 break;
242 default:
244 }
245
246 for (unsigned int k = 0; k < quadrature->size(); ++k)
247 this->generalized_support_points[current++] = quadrature->point(k);
248 }
249 Assert(current == this->n_dofs_per_cell(), ExcInternalError());
250}
251
252
253
254template <int dim>
255std::vector<unsigned int>
257{
258 // We have (deg+1)^(dim-1) DoFs per face...
259 unsigned int dofs_per_face = 1;
260 for (unsigned int d = 1; d < dim; ++d)
261 dofs_per_face *= deg + 1;
262
263 // ...plus the interior DoFs for the total of dim*(deg+1)^dim
264 const unsigned int interior_dofs =
265 dim * (deg - 1) * Utilities::pow(deg + 1, dim - 1);
266
267 std::vector<unsigned int> dpo(dim + 1);
268 dpo[dim - 1] = dofs_per_face;
269 dpo[dim] = interior_dofs;
270
271 return dpo;
272}
273
274
275
276template <>
277std::vector<bool>
279{
280 Assert(false, ExcImpossibleInDim(1));
281 return std::vector<bool>();
282}
283
284
285
286template <int dim>
287std::vector<bool>
289{
290 const unsigned int dofs_per_cell =
292 unsigned int dofs_per_face = deg + 1;
293 for (unsigned int d = 2; d < dim; ++d)
294 dofs_per_face *= deg + 1;
295 // All face dofs need to be non-additive, since they have
296 // continuity requirements. The interior dofs are
297 // made additive.
298 std::vector<bool> ret_val(dofs_per_cell, false);
299 for (unsigned int i = GeometryInfo<dim>::faces_per_cell * dofs_per_face;
300 i < dofs_per_cell;
301 ++i)
302 ret_val[i] = true;
303
304 return ret_val;
305}
306
307
308
309template <int dim>
310void
312 const std::vector<Vector<double>> &support_point_values,
313 std::vector<double> &nodal_values) const
314{
315 Assert(support_point_values.size() == this->generalized_support_points.size(),
316 ExcDimensionMismatch(support_point_values.size(),
317 this->generalized_support_points.size()));
318 Assert(nodal_values.size() == this->n_dofs_per_cell(),
319 ExcDimensionMismatch(nodal_values.size(), this->n_dofs_per_cell()));
320 Assert(support_point_values[0].size() == this->n_components(),
321 ExcDimensionMismatch(support_point_values[0].size(),
322 this->n_components()));
323
324 // First do interpolation on faces. There, the component
325 // evaluated depends on the face direction and orientation.
326 unsigned int fbase = 0;
327 unsigned int f = 0;
328 for (; f < GeometryInfo<dim>::faces_per_cell;
329 ++f, fbase += this->n_dofs_per_face(f))
330 {
331 for (unsigned int i = 0; i < this->n_dofs_per_face(f); ++i)
332 {
333 nodal_values[fbase + i] = support_point_values[fbase + i](
335 }
336 }
337
338 // The remaining points form dim chunks, one for each component.
339 const unsigned int istep = (this->n_dofs_per_cell() - fbase) / dim;
340 Assert((this->n_dofs_per_cell() - fbase) % dim == 0, ExcInternalError());
341
342 f = 0;
343 while (fbase < this->n_dofs_per_cell())
344 {
345 for (unsigned int i = 0; i < istep; ++i)
346 {
347 nodal_values[fbase + i] = support_point_values[fbase + i](f);
348 }
349 fbase += istep;
350 ++f;
351 }
352 Assert(fbase == this->n_dofs_per_cell(), ExcInternalError());
353}
354
355
356
357// explicit instantiations
358#include "fe/fe_rt_bubbles.inst"
359
360
FullMatrix< double > inverse_node_matrix
std::vector< MappingKind > mapping_kind
FE_RT_Bubbles(const unsigned int k)
virtual std::string get_name() const override
virtual std::unique_ptr< FiniteElement< dim, dim > > clone() const override
void initialize_support_points(const unsigned int rt_degree)
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 std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
static std::vector< bool > get_ria_vector(const unsigned int degree)
void initialize_quad_dof_index_permutation_and_sign_change()
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
FullMatrix< double > interface_constraints
Definition fe.h:2549
std::vector< std::vector< FullMatrix< double > > > prolongation
Definition fe.h:2537
size_type n() const
void invert(const FullMatrix< number2 > &M)
size_type m() const
static unsigned int n_polynomials(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)
const Point< dim > & point(const unsigned int i) const
static std::array< RefinementCase< dim >, n_refinement_cases > all_refinement_cases()
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:35
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
#define DEAL_II_NOT_IMPLEMENTED()
#define Assert(cond, exc)
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
@ mapping_raviart_thomas
Definition mapping.h:136
std::size_t size
Definition mpi.cc:745
void compute_embedding_matrices(const FiniteElement< dim, spacedim > &fe, std::vector< std::vector< FullMatrix< number > > > &matrices, const bool isotropic_only=false, const double threshold=1.e-12)
FullMatrix< double > compute_node_matrix(const FiniteElement< dim, spacedim > &fe)
constexpr T pow(const T base, const int iexp)
Definition utilities.h:967
constexpr types::geometric_orientation default_geometric_orientation
Definition types.h:352
STL namespace.