Reference documentation for deal.II version GIT f0f8c7fe18 2023-03-21 21:25: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\}}\)
fe_rt_bubbles.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2018 - 2022 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>
25 #include <deal.II/fe/fe_values.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 
41 template <int dim>
42 FE_RT_Bubbles<dim>::FE_RT_Bubbles(const unsigned int deg)
43  : FE_PolyTensor<dim>(
44  PolynomialsRT_Bubbles<dim>(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,
56  ExcMessage(
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);
67  this->inverse_node_matrix.invert(M);
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
84  AssertDimension(this->n_unique_faces(), 1);
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
89  FETools::compute_embedding_matrices(*this, this->prolongation, true, 1.0);
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 
116 template <int dim>
117 void
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 
130 template <int dim>
131 std::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 
144 template <int dim>
145 std::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 
158 template <int dim>
159 void
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),
181  ExcInternalError());
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 =
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 
251 template <int dim>
252 std::vector<unsigned int>
253 FE_RT_Bubbles<dim>::get_dpo_vector(const unsigned int deg)
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 
273 template <>
274 std::vector<bool>
276 {
277  Assert(false, ExcImpossibleInDim(1));
278  return std::vector<bool>();
279 }
280 
281 
282 
283 template <int dim>
284 std::vector<bool>
285 FE_RT_Bubbles<dim>::get_ria_vector(const unsigned int deg)
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 
306 template <int dim>
307 void
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
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
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)
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 std::vector< Point< dim > > & get_points() const
const Point< dim > & point(const unsigned int i) const
Definition: vector.h:109
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1586
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1759
static ::ExceptionBase & ExcMessage(std::string arg1)
@ mapping_raviart_thomas
Definition: mapping.h:130
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)
std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
constexpr T pow(const T base, const int iexp)
Definition: utilities.h:447
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618
static unsigned int n_children(const RefinementCase< dim > &refinement_case)