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