Reference documentation for deal.II version 9.5.0
\(\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// Copyright (C) 2018 - 2023 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16
19
21
22#include <deal.II/fe/fe.h>
24#include <deal.II/fe/fe_tools.h>
26#include <deal.II/fe/mapping.h>
27
28#include <deal.II/grid/tria.h>
30
31#include <memory>
32#include <sstream>
33
34
36
37// TODO: implement the adjust_quad_dof_index_for_face_orientation_table and
38// adjust_line_dof_index_for_line_orientation_table fields, and write tests
39// similar to bits/face_orientation_and_fe_q_*
40
41template <int dim>
42FE_RT_Bubbles<dim>::FE_RT_Bubbles(const unsigned int deg)
43 : FE_PolyTensor<dim>(
44 PolynomialsRT_Bubbles<dim>(deg),
45 FiniteElementData<dim>(get_dpo_vector(deg),
46 dim,
47 deg + 1,
48 FiniteElementData<dim>::Hdiv),
49 get_ria_vector(deg),
50 std::vector<ComponentMask>(PolynomialsRT_Bubbles<dim>::n_polynomials(deg),
51 std::vector<bool>(dim, true)))
52{
53 Assert(dim >= 2, ExcImpossibleInDim(dim));
54 Assert(
55 deg >= 1,
57 "Lowest order RT_Bubbles element is degree 1, but you requested for degree 0"));
58 const unsigned int n_dofs = this->n_dofs_per_cell();
59
61 // Initialize support points and quadrature weights
63 // Compute the inverse node matrix to get
64 // the correct basis functions
66 this->inverse_node_matrix.reinit(n_dofs, n_dofs);
68
69 // Reinit the vectors of prolongation matrices to the
70 // right sizes. There are no restriction matrices implemented
71 for (unsigned int ref_case = RefinementCase<dim>::cut_x;
72 ref_case < RefinementCase<dim>::isotropic_refinement + 1;
73 ++ref_case)
74 {
75 const unsigned int nc =
77
78 for (unsigned int i = 0; i < nc; ++i)
79 this->prolongation[ref_case - 1][i].reinit(n_dofs, n_dofs);
80 }
81
82 // TODO: the implementation makes the assumption that all faces have the
83 // same number of dofs
85 const unsigned int face_no = 0;
86
87 // Fill prolongation matrices with embedding operators
88 // set tolerance to 1, as embedding error accumulate quickly
91 for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_face; ++i)
92 face_embeddings[i].reinit(this->n_dofs_per_face(face_no),
93 this->n_dofs_per_face(face_no));
94 FETools::compute_face_embedding_matrices<dim, double>(*this,
95 face_embeddings,
96 0,
97 0);
98 this->interface_constraints.reinit((1 << (dim - 1)) *
99 this->n_dofs_per_face(face_no),
100 this->n_dofs_per_face(face_no));
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)
104 {
105 for (unsigned int j = 0; j < face_embeddings[d].n(); ++j)
106 this->interface_constraints(target_row, j) = face_embeddings[d](i, j);
107 ++target_row;
108 }
109
110 // We need to initialize the dof permutation table and the one for the sign
111 // change.
113}
114
115
116template <int dim>
117void
119{
120 // for 1d and 2d, do nothing
121 if (dim < 3)
122 return;
123
124 // TODO: Implement this for this class
125 return;
126}
127
128
129
130template <int dim>
131std::string
133{
134 // Note: this->degree is the maximal polynomial degree and is thus one higher
135 // than the argument given to the constructor
136 std::ostringstream namebuf;
137 namebuf << "FE_RT_Bubbles<" << dim << ">(" << this->degree << ")";
138
139 return namebuf.str();
140}
141
142
143
144template <int dim>
145std::unique_ptr<FiniteElement<dim, dim>>
147{
148 return std::make_unique<FE_RT_Bubbles<dim>>(*this);
149}
150
151
152//---------------------------------------------------------------------------
153// Auxiliary and internal functions
154//---------------------------------------------------------------------------
155
156
157
158template <int dim>
159void
161{
162 // TODO: the implementation makes the assumption that all faces have the
163 // same number of dofs
164 AssertDimension(this->n_unique_faces(), 1);
165 const unsigned int face_no = 0;
166
167 this->generalized_support_points.resize(this->n_dofs_per_cell());
168 this->generalized_face_support_points[face_no].resize(
169 this->n_dofs_per_face(face_no));
170
171 // Index of the point being entered
172 unsigned int current = 0;
173
174 // On the faces, we choose as many Gauss-Lobatto points
175 // as required to determine the normal component uniquely.
176 // This is the deg of the RT_Bubble element plus one.
177 if (dim > 1)
178 {
179 QGaussLobatto<dim - 1> face_points(deg + 1);
180 Assert(face_points.size() == this->n_dofs_per_face(face_no),
182 for (unsigned int k = 0; k < this->n_dofs_per_face(face_no); ++k)
183 this->generalized_face_support_points[face_no][k] =
184 face_points.point(k);
185 Quadrature<dim> faces =
186 QProjector<dim>::project_to_all_faces(this->reference_cell(),
187 face_points);
188 for (unsigned int k = 0; k < this->n_dofs_per_face(face_no) *
190 ++k)
191 this->generalized_support_points[k] = faces.point(
192 k + QProjector<dim>::DataSetDescriptor::face(this->reference_cell(),
193 0,
194 true,
195 false,
196 false,
197 this->n_dofs_per_face(
198 face_no)));
199
200 current =
201 this->n_dofs_per_face(face_no) * GeometryInfo<dim>::faces_per_cell;
202 }
203
204 if (deg == 1)
205 return;
206
207 // In the interior, we need anisotropic Gauss-Lobatto quadratures,
208 // one for each direction
209 QGaussLobatto<1> high(deg + 1);
210 std::vector<Point<1>> pts = high.get_points();
211 if (pts.size() > 2)
212 {
213 pts.erase(pts.begin());
214 pts.erase(pts.end() - 1);
215 }
216
217 std::vector<double> wts(pts.size(), 1);
218 Quadrature<1> low(pts, wts);
219
220 for (unsigned int d = 0; d < dim; ++d)
221 {
222 std::unique_ptr<QAnisotropic<dim>> quadrature;
223 switch (dim)
224 {
225 case 1:
226 quadrature = std::make_unique<QAnisotropic<dim>>(high);
227 break;
228 case 2:
229 quadrature =
230 std::make_unique<QAnisotropic<dim>>(((d == 0) ? low : high),
231 ((d == 1) ? low : high));
232 break;
233 case 3:
234 quadrature =
235 std::make_unique<QAnisotropic<dim>>(((d == 0) ? low : high),
236 ((d == 1) ? low : high),
237 ((d == 2) ? low : high));
238 break;
239 default:
240 Assert(false, ExcNotImplemented());
241 }
242
243 for (unsigned int k = 0; k < quadrature->size(); ++k)
244 this->generalized_support_points[current++] = quadrature->point(k);
245 }
246 Assert(current == this->n_dofs_per_cell(), ExcInternalError());
247}
248
249
250
251template <int dim>
252std::vector<unsigned int>
254{
255 // We have (deg+1)^(dim-1) DoFs per face...
256 unsigned int dofs_per_face = 1;
257 for (unsigned int d = 1; d < dim; ++d)
258 dofs_per_face *= deg + 1;
259
260 // ...plus the interior DoFs for the total of dim*(deg+1)^dim
261 const unsigned int interior_dofs =
262 dim * (deg - 1) * Utilities::pow(deg + 1, dim - 1);
263
264 std::vector<unsigned int> dpo(dim + 1);
265 dpo[dim - 1] = dofs_per_face;
266 dpo[dim] = interior_dofs;
267
268 return dpo;
269}
270
271
272
273template <>
274std::vector<bool>
276{
277 Assert(false, ExcImpossibleInDim(1));
278 return std::vector<bool>();
279}
280
281
282
283template <int dim>
284std::vector<bool>
286{
287 const unsigned int dofs_per_cell =
289 unsigned int dofs_per_face = deg + 1;
290 for (unsigned int d = 2; d < dim; ++d)
291 dofs_per_face *= deg + 1;
292 // All face dofs need to be non-additive, since they have
293 // continuity requirements. The interior dofs are
294 // made additive.
295 std::vector<bool> ret_val(dofs_per_cell, false);
296 for (unsigned int i = GeometryInfo<dim>::faces_per_cell * dofs_per_face;
297 i < dofs_per_cell;
298 ++i)
299 ret_val[i] = true;
300
301 return ret_val;
302}
303
304
305
306template <int dim>
307void
309 const std::vector<Vector<double>> &support_point_values,
310 std::vector<double> & nodal_values) const
311{
312 Assert(support_point_values.size() == this->generalized_support_points.size(),
313 ExcDimensionMismatch(support_point_values.size(),
314 this->generalized_support_points.size()));
315 Assert(nodal_values.size() == this->n_dofs_per_cell(),
316 ExcDimensionMismatch(nodal_values.size(), this->n_dofs_per_cell()));
317 Assert(support_point_values[0].size() == this->n_components(),
318 ExcDimensionMismatch(support_point_values[0].size(),
319 this->n_components()));
320
321 // First do interpolation on faces. There, the component
322 // evaluated depends on the face direction and orientation.
323 unsigned int fbase = 0;
324 unsigned int f = 0;
325 for (; f < GeometryInfo<dim>::faces_per_cell;
326 ++f, fbase += this->n_dofs_per_face(f))
327 {
328 for (unsigned int i = 0; i < this->n_dofs_per_face(f); ++i)
329 {
330 nodal_values[fbase + i] = support_point_values[fbase + i](
332 }
333 }
334
335 // The remaining points form dim chunks, one for each component.
336 const unsigned int istep = (this->n_dofs_per_cell() - fbase) / dim;
337 Assert((this->n_dofs_per_cell() - fbase) % dim == 0, ExcInternalError());
338
339 f = 0;
340 while (fbase < this->n_dofs_per_cell())
341 {
342 for (unsigned int i = 0; i < istep; ++i)
343 {
344 nodal_values[fbase + i] = support_point_values[fbase + i](f);
345 }
346 fbase += istep;
347 ++f;
348 }
349 Assert(fbase == this->n_dofs_per_cell(), ExcInternalError());
350}
351
352
353
354// explicit instantiations
355#include "fe_rt_bubbles.inst"
356
357
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:2428
std::vector< std::vector< FullMatrix< double > > > prolongation
Definition fe.h:2416
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
const std::vector< Point< dim > > & get_points() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
static ::ExceptionBase & ExcNotImplemented()
#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:133
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:447
STL namespace.
static unsigned int n_children(const RefinementCase< dim > &refinement_case)