Reference documentation for deal.II version GIT relicensing-384-g4769914377 2024-04-13 09:00:02+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\}}\)
Loading...
Searching...
No Matches
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 =
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 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 k = 0; k < this->n_dofs_per_face(face_no) *
189 ++k)
190 this->generalized_support_points[k] =
192 this->reference_cell(),
193 0,
195 this->n_dofs_per_face(face_no)));
196
197 current =
198 this->n_dofs_per_face(face_no) * GeometryInfo<dim>::faces_per_cell;
199 }
200
201 if (deg == 1)
202 return;
203
204 // In the interior, we need anisotropic Gauss-Lobatto quadratures,
205 // one for each direction
206 QGaussLobatto<1> high(deg + 1);
207 std::vector<Point<1>> pts = high.get_points();
208 if (pts.size() > 2)
209 {
210 pts.erase(pts.begin());
211 pts.erase(pts.end() - 1);
212 }
213
214 std::vector<double> wts(pts.size(), 1);
215 Quadrature<1> low(pts, wts);
216
217 for (unsigned int d = 0; d < dim; ++d)
218 {
219 std::unique_ptr<QAnisotropic<dim>> quadrature;
220 switch (dim)
221 {
222 case 1:
223 quadrature = std::make_unique<QAnisotropic<dim>>(high);
224 break;
225 case 2:
226 quadrature =
227 std::make_unique<QAnisotropic<dim>>(((d == 0) ? low : high),
228 ((d == 1) ? low : high));
229 break;
230 case 3:
231 quadrature =
232 std::make_unique<QAnisotropic<dim>>(((d == 0) ? low : high),
233 ((d == 1) ? low : high),
234 ((d == 2) ? low : high));
235 break;
236 default:
238 }
239
240 for (unsigned int k = 0; k < quadrature->size(); ++k)
241 this->generalized_support_points[current++] = quadrature->point(k);
242 }
243 Assert(current == this->n_dofs_per_cell(), ExcInternalError());
244}
245
246
247
248template <int dim>
249std::vector<unsigned int>
251{
252 // We have (deg+1)^(dim-1) DoFs per face...
253 unsigned int dofs_per_face = 1;
254 for (unsigned int d = 1; d < dim; ++d)
255 dofs_per_face *= deg + 1;
256
257 // ...plus the interior DoFs for the total of dim*(deg+1)^dim
258 const unsigned int interior_dofs =
259 dim * (deg - 1) * Utilities::pow(deg + 1, dim - 1);
260
261 std::vector<unsigned int> dpo(dim + 1);
262 dpo[dim - 1] = dofs_per_face;
263 dpo[dim] = interior_dofs;
264
265 return dpo;
266}
267
268
269
270template <>
271std::vector<bool>
273{
274 Assert(false, ExcImpossibleInDim(1));
275 return std::vector<bool>();
276}
277
278
279
280template <int dim>
281std::vector<bool>
283{
284 const unsigned int dofs_per_cell =
286 unsigned int dofs_per_face = deg + 1;
287 for (unsigned int d = 2; d < dim; ++d)
288 dofs_per_face *= deg + 1;
289 // All face dofs need to be non-additive, since they have
290 // continuity requirements. The interior dofs are
291 // made additive.
292 std::vector<bool> ret_val(dofs_per_cell, false);
293 for (unsigned int i = GeometryInfo<dim>::faces_per_cell * dofs_per_face;
294 i < dofs_per_cell;
295 ++i)
296 ret_val[i] = true;
297
298 return ret_val;
299}
300
301
302
303template <int dim>
304void
306 const std::vector<Vector<double>> &support_point_values,
307 std::vector<double> &nodal_values) const
308{
309 Assert(support_point_values.size() == this->generalized_support_points.size(),
310 ExcDimensionMismatch(support_point_values.size(),
311 this->generalized_support_points.size()));
312 Assert(nodal_values.size() == this->n_dofs_per_cell(),
313 ExcDimensionMismatch(nodal_values.size(), this->n_dofs_per_cell()));
314 Assert(support_point_values[0].size() == this->n_components(),
315 ExcDimensionMismatch(support_point_values[0].size(),
316 this->n_components()));
317
318 // First do interpolation on faces. There, the component
319 // evaluated depends on the face direction and orientation.
320 unsigned int fbase = 0;
321 unsigned int f = 0;
322 for (; f < GeometryInfo<dim>::faces_per_cell;
323 ++f, fbase += this->n_dofs_per_face(f))
324 {
325 for (unsigned int i = 0; i < this->n_dofs_per_face(f); ++i)
326 {
327 nodal_values[fbase + i] = support_point_values[fbase + i](
329 }
330 }
331
332 // The remaining points form dim chunks, one for each component.
333 const unsigned int istep = (this->n_dofs_per_cell() - fbase) / dim;
334 Assert((this->n_dofs_per_cell() - fbase) % dim == 0, ExcInternalError());
335
336 f = 0;
337 while (fbase < this->n_dofs_per_cell())
338 {
339 for (unsigned int i = 0; i < istep; ++i)
340 {
341 nodal_values[fbase + i] = support_point_values[fbase + i](f);
342 }
343 fbase += istep;
344 ++f;
345 }
346 Assert(fbase == this->n_dofs_per_cell(), ExcInternalError());
347}
348
349
350
351// explicit instantiations
352#include "fe_rt_bubbles.inst"
353
354
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
unsigned int n_unique_faces() const
FullMatrix< double > interface_constraints
Definition fe.h:2423
std::vector< std::vector< FullMatrix< double > > > prolongation
Definition fe.h:2411
size_type n() const
void invert(const FullMatrix< number2 > &M)
size_type m() const
static unsigned int n_polynomials(const unsigned int degree)
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 constexpr unsigned char default_combined_face_orientation()
static std::array< RefinementCase< dim >, n_refinement_cases > all_refinement_cases()
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
#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:134
#define DEAL_II_NOT_IMPLEMENTED()
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:966
STL namespace.
static unsigned int n_children(const RefinementCase< dim > &refinement_case)