Reference documentation for deal.II version 9.2.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\}}\)
fe_rt_bubbles.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2018 - 2019 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 
20 
22 
23 #include <deal.II/fe/fe.h>
25 #include <deal.II/fe/fe_tools.h>
26 #include <deal.II/fe/fe_values.h>
27 #include <deal.II/fe/mapping.h>
28 
29 #include <deal.II/grid/tria.h>
31 
32 #include <sstream>
33 
34 
36 
37 
38 template <int dim>
39 FE_RT_Bubbles<dim>::FE_RT_Bubbles(const unsigned int deg)
40  : FE_PolyTensor<dim>(
41  PolynomialsRT_Bubbles<dim>(deg),
42  FiniteElementData<dim>(get_dpo_vector(deg),
43  dim,
44  deg + 1,
45  FiniteElementData<dim>::Hdiv),
46  get_ria_vector(deg),
47  std::vector<ComponentMask>(PolynomialsRT_Bubbles<dim>::n_polynomials(deg),
48  std::vector<bool>(dim, true)))
49 {
50  Assert(dim >= 2, ExcImpossibleInDim(dim));
51  Assert(
52  deg >= 1,
53  ExcMessage(
54  "Lowest order RT_Bubbles element is degree 1, but you requested for degree 0"));
55  const unsigned int n_dofs = this->dofs_per_cell;
56 
58  // Initialize support points and quadrature weights
60  // Compute the inverse node matrix to get
61  // the correct basis functions
63  this->inverse_node_matrix.reinit(n_dofs, n_dofs);
64  this->inverse_node_matrix.invert(M);
65 
66  // Reinit the vectors of prolongation matrices to the
67  // right sizes. There are no restriction matrices implemented
68  for (unsigned int ref_case = RefinementCase<dim>::cut_x;
69  ref_case < RefinementCase<dim>::isotropic_refinement + 1;
70  ++ref_case)
71  {
72  const unsigned int nc =
74 
75  for (unsigned int i = 0; i < nc; ++i)
76  this->prolongation[ref_case - 1][i].reinit(n_dofs, n_dofs);
77  }
78  // Fill prolongation matrices with embedding operators
79  // set tolerance to 1, as embedding error accumulate quickly
80  FETools::compute_embedding_matrices(*this, this->prolongation, true, 1.0);
82  for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_face; ++i)
83  face_embeddings[i].reinit(this->dofs_per_face, this->dofs_per_face);
84  FETools::compute_face_embedding_matrices<dim, double>(*this,
85  face_embeddings,
86  0,
87  0);
88  this->interface_constraints.reinit((1 << (dim - 1)) * this->dofs_per_face,
89  this->dofs_per_face);
90  unsigned int target_row = 0;
91  for (unsigned int d = 0; d < GeometryInfo<dim>::max_children_per_face; ++d)
92  for (unsigned int i = 0; i < face_embeddings[d].m(); ++i)
93  {
94  for (unsigned int j = 0; j < face_embeddings[d].n(); ++j)
95  this->interface_constraints(target_row, j) = face_embeddings[d](i, j);
96  ++target_row;
97  }
98 }
99 
100 
101 
102 template <int dim>
103 std::string
105 {
106  // Note: this->degree is the maximal polynomial degree and is thus one higher
107  // than the argument given to the constructor
108  std::ostringstream namebuf;
109  namebuf << "FE_RT_Bubbles<" << dim << ">(" << this->degree << ")";
110 
111  return namebuf.str();
112 }
113 
114 
115 
116 template <int dim>
117 std::unique_ptr<FiniteElement<dim, dim>>
119 {
120  return std_cxx14::make_unique<FE_RT_Bubbles<dim>>(*this);
121 }
122 
123 
124 //---------------------------------------------------------------------------
125 // Auxiliary and internal functions
126 //---------------------------------------------------------------------------
127 
128 
129 
130 template <int dim>
131 void
133 {
134  this->generalized_support_points.resize(this->dofs_per_cell);
135  this->generalized_face_support_points.resize(this->dofs_per_face);
136 
137  // Index of the point being entered
138  unsigned int current = 0;
139 
140  // On the faces, we choose as many Gauss-Lobatto points
141  // as required to determine the normal component uniquely.
142  // This is the deg of the RT_Bubble element plus one.
143  if (dim > 1)
144  {
145  QGaussLobatto<dim - 1> face_points(deg + 1);
146  Assert(face_points.size() == this->dofs_per_face, ExcInternalError());
147  for (unsigned int k = 0; k < this->dofs_per_face; ++k)
148  this->generalized_face_support_points[k] = face_points.point(k);
149  Quadrature<dim> faces =
151  for (unsigned int k = 0;
152  k < this->dofs_per_face * GeometryInfo<dim>::faces_per_cell;
153  ++k)
154  this->generalized_support_points[k] =
156  0, true, false, false, this->dofs_per_face));
157 
158  current = this->dofs_per_face * GeometryInfo<dim>::faces_per_cell;
159  }
160 
161  if (deg == 1)
162  return;
163 
164  // In the interior, we need anisotropic Gauss-Lobatto quadratures,
165  // one for each direction
166  QGaussLobatto<1> high(deg + 1);
167  std::vector<Point<1>> pts = high.get_points();
168  pts.erase(pts.begin());
169  pts.erase(pts.end() - 1);
170 
171  std::vector<double> wts(pts.size(), 1);
172  Quadrature<1> low(pts, wts);
173 
174  for (unsigned int d = 0; d < dim; ++d)
175  {
176  std::unique_ptr<QAnisotropic<dim>> quadrature;
177  switch (dim)
178  {
179  case 1:
180  quadrature = std_cxx14::make_unique<QAnisotropic<dim>>(high);
181  break;
182  case 2:
183  quadrature = std_cxx14::make_unique<QAnisotropic<dim>>(
184  ((d == 0) ? low : high), ((d == 1) ? low : high));
185  break;
186  case 3:
187  quadrature =
188  std_cxx14::make_unique<QAnisotropic<dim>>(((d == 0) ? low : high),
189  ((d == 1) ? low : high),
190  ((d == 2) ? low :
191  high));
192  break;
193  default:
194  Assert(false, ExcNotImplemented());
195  }
196 
197  for (unsigned int k = 0; k < quadrature->size(); ++k)
198  this->generalized_support_points[current++] = quadrature->point(k);
199  }
200  Assert(current == this->dofs_per_cell, ExcInternalError());
201 }
202 
203 
204 
205 template <int dim>
206 std::vector<unsigned int>
207 FE_RT_Bubbles<dim>::get_dpo_vector(const unsigned int deg)
208 {
209  // We have (deg+1)^(dim-1) DoFs per face...
210  unsigned int dofs_per_face = 1;
211  for (unsigned int d = 1; d < dim; ++d)
212  dofs_per_face *= deg + 1;
213 
214  // ...plus the interior DoFs for the total of dim*(deg+1)^dim
215  const unsigned int interior_dofs =
216  dim * (deg - 1) * Utilities::pow(deg + 1, dim - 1);
217 
218  std::vector<unsigned int> dpo(dim + 1);
219  dpo[dim - 1] = dofs_per_face;
220  dpo[dim] = interior_dofs;
221 
222  return dpo;
223 }
224 
225 
226 
227 template <>
228 std::vector<bool>
230 {
231  Assert(false, ExcImpossibleInDim(1));
232  return std::vector<bool>();
233 }
234 
235 
236 
237 template <int dim>
238 std::vector<bool>
239 FE_RT_Bubbles<dim>::get_ria_vector(const unsigned int deg)
240 {
241  const unsigned int dofs_per_cell =
243  unsigned int dofs_per_face = deg + 1;
244  for (unsigned int d = 2; d < dim; ++d)
245  dofs_per_face *= deg + 1;
246  // All face dofs need to be non-additive, since they have
247  // continuity requirements. The interior dofs are
248  // made additive.
249  std::vector<bool> ret_val(dofs_per_cell, false);
250  for (unsigned int i = GeometryInfo<dim>::faces_per_cell * dofs_per_face;
251  i < dofs_per_cell;
252  ++i)
253  ret_val[i] = true;
254 
255  return ret_val;
256 }
257 
258 
259 
260 template <int dim>
261 void
263  const std::vector<Vector<double>> &support_point_values,
264  std::vector<double> & nodal_values) const
265 {
266  Assert(support_point_values.size() == this->generalized_support_points.size(),
267  ExcDimensionMismatch(support_point_values.size(),
268  this->generalized_support_points.size()));
269  Assert(nodal_values.size() == this->dofs_per_cell,
270  ExcDimensionMismatch(nodal_values.size(), this->dofs_per_cell));
271  Assert(support_point_values[0].size() == this->n_components(),
272  ExcDimensionMismatch(support_point_values[0].size(),
273  this->n_components()));
274 
275  // First do interpolation on faces. There, the component
276  // evaluated depends on the face direction and orientation.
277  unsigned int fbase = 0;
278  unsigned int f = 0;
279  for (; f < GeometryInfo<dim>::faces_per_cell;
280  ++f, fbase += this->dofs_per_face)
281  {
282  for (unsigned int i = 0; i < this->dofs_per_face; ++i)
283  {
284  nodal_values[fbase + i] = support_point_values[fbase + i](
286  }
287  }
288 
289  // The remaining points form dim chunks, one for each component.
290  const unsigned int istep = (this->dofs_per_cell - fbase) / dim;
291  Assert((this->dofs_per_cell - fbase) % dim == 0, ExcInternalError());
292 
293  f = 0;
294  while (fbase < this->dofs_per_cell)
295  {
296  for (unsigned int i = 0; i < istep; ++i)
297  {
298  nodal_values[fbase + i] = support_point_values[fbase + i](f);
299  }
300  fbase += istep;
301  ++f;
302  }
303  Assert(fbase == this->dofs_per_cell, ExcInternalError());
304 }
305 
306 
307 
308 // explicit instantiations
309 #include "fe_rt_bubbles.inst"
310 
311 
Utilities::pow
constexpr T pow(const T base, const int iexp)
Definition: utilities.h:476
fe_values.h
fe_rt_bubbles.h
GeometryInfo::n_children
static unsigned int n_children(const RefinementCase< dim > &refinement_case)
PolynomialsRT_Bubbles
Definition: polynomials_rt_bubbles.h:88
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
tria.h
FiniteElementData
Definition: fe_base.h:147
FiniteElement< dim, dim >::interface_constraints
FullMatrix< double > interface_constraints
Definition: fe.h:2440
FullMatrix::m
size_type m() const
tria_iterator.h
GeometryInfo
Definition: geometry_info.h:1224
FE_RT_Bubbles::clone
virtual std::unique_ptr< FiniteElement< dim, dim > > clone() const override
Definition: fe_rt_bubbles.cc:118
TableBase::reinit
void reinit(const TableIndices< N > &new_size, const bool omit_default_initialization=false)
FE_PolyTensor< dim >::mapping_kind
std::vector< MappingKind > mapping_kind
Definition: fe_poly_tensor.h:220
DoFHandler::n_components
unsigned int n_components(const DoFHandler< dim, spacedim > &dh)
ComponentMask
Definition: component_mask.h:83
FullMatrix::n
size_type n() const
FETools::compute_embedding_matrices
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)
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
FE_RT_Bubbles::convert_generalized_support_point_values_to_dof_values
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
Definition: fe_rt_bubbles.cc:262
quadrature_lib.h
QProjector::project_to_all_faces
static Quadrature< dim > project_to_all_faces(const SubQuadrature &quadrature)
QProjector
Definition: qprojector.h:78
bool
mapping_raviart_thomas
@ mapping_raviart_thomas
Definition: mapping.h:118
Quadrature::point
const Point< dim > & point(const unsigned int i) const
FETools::compute_node_matrix
FullMatrix< double > compute_node_matrix(const FiniteElement< dim, spacedim > &fe)
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
internal::reinit
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:621
FiniteElementData::dofs_per_cell
const unsigned int dofs_per_cell
Definition: fe_base.h:282
PolynomialsRT_Bubbles::n_polynomials
static unsigned int n_polynomials(const unsigned int degree)
Definition: polynomials_rt_bubbles.cc:838
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
FiniteElement< dim, dim >::prolongation
std::vector< std::vector< FullMatrix< double > > > prolongation
Definition: fe.h:2428
RefinementCase
Definition: geometry_info.h:795
FE_PolyTensor< dim >::inverse_node_matrix
FullMatrix< double > inverse_node_matrix
Definition: fe_poly_tensor.h:490
fe.h
FE_RT_Bubbles::get_dpo_vector
static std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
Definition: fe_rt_bubbles.cc:207
QGaussLobatto
Definition: quadrature_lib.h:76
fe_tools.h
FiniteElementData::dofs_per_face
const unsigned int dofs_per_face
Definition: fe_base.h:275
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
mapping.h
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1419
StandardExceptions::ExcImpossibleInDim
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
dof_accessor.h
FE_RT_Bubbles::initialize_support_points
void initialize_support_points(const unsigned int rt_degree)
Definition: fe_rt_bubbles.cc:132
FullMatrix::invert
void invert(const FullMatrix< number2 > &M)
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
FE_PolyTensor
Definition: fe_poly_tensor.h:143
FullMatrix< double >
FE_RT_Bubbles::FE_RT_Bubbles
FE_RT_Bubbles(const unsigned int k)
Definition: fe_rt_bubbles.cc:39
memory.h
Quadrature::get_points
const std::vector< Point< dim > > & get_points() const
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
FE_RT_Bubbles::get_name
virtual std::string get_name() const override
Definition: fe_rt_bubbles.cc:104
Quadrature
Definition: quadrature.h:85
Vector< double >
qprojector.h
FE_RT_Bubbles::get_ria_vector
static std::vector< bool > get_ria_vector(const unsigned int degree)
Definition: fe_rt_bubbles.cc:239