Reference documentation for deal.II version Git a73d35cdfe 2020-05-28 10:34:58 +0200
\(\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_raviart_thomas.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2003 - 2020 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 #include <deal.II/base/table.h>
21 #include <deal.II/base/utilities.h>
22 
24 
25 #include <deal.II/fe/fe.h>
27 #include <deal.II/fe/fe_tools.h>
28 #include <deal.II/fe/fe_values.h>
29 #include <deal.II/fe/mapping.h>
30 
31 #include <deal.II/grid/tria.h>
33 
34 #include <iostream>
35 #include <memory>
36 #include <sstream>
37 
38 // TODO: implement the adjust_quad_dof_index_for_face_orientation_table and
39 // adjust_line_dof_index_for_line_orientation_table fields, and write tests
40 // similar to bits/face_orientation_and_fe_q_*
41 
42 
44 
45 
46 template <int dim>
48  : FE_PolyTensor<dim>(
49  PolynomialsRaviartThomas<dim>(deg),
50  FiniteElementData<dim>(get_dpo_vector(deg),
51  dim,
52  deg + 1,
53  FiniteElementData<dim>::Hdiv),
54  std::vector<bool>(PolynomialsRaviartThomas<dim>::n_polynomials(deg),
55  true),
56  std::vector<ComponentMask>(PolynomialsRaviartThomas<dim>::n_polynomials(
57  deg),
58  std::vector<bool>(dim, true)))
59 {
60  Assert(dim >= 2, ExcImpossibleInDim(dim));
61  const unsigned int n_dofs = this->dofs_per_cell;
62 
64  // First, initialize the
65  // generalized support points and
66  // quadrature weights, since they
67  // are required for interpolation.
69 
70  // Now compute the inverse node matrix, generating the correct
71  // basis functions from the raw ones. For a discussion of what
72  // exactly happens here, see FETools::compute_node_matrix.
74  this->inverse_node_matrix.reinit(n_dofs, n_dofs);
75  this->inverse_node_matrix.invert(M);
76  // From now on, the shape functions provided by FiniteElement::shape_value
77  // and similar functions will be the correct ones, not
78  // the raw shape functions from the polynomial space anymore.
79 
80  // Reinit the vectors of
81  // restriction and prolongation
82  // matrices to the right sizes.
83  // Restriction only for isotropic
84  // refinement
86  // Fill prolongation matrices with embedding operators
89 
90  // TODO[TL]: for anisotropic refinement we will probably need a table of
91  // submatrices with an array for each refine case
93  for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_face; ++i)
94  face_embeddings[i].reinit(this->dofs_per_face, this->dofs_per_face);
95  FETools::compute_face_embedding_matrices<dim, double>(*this,
96  face_embeddings,
97  0,
98  0);
99  this->interface_constraints.reinit((1 << (dim - 1)) * this->dofs_per_face,
100  this->dofs_per_face);
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 
111 
112 
113 template <int dim>
114 std::string
116 {
117  // note that the
118  // FETools::get_fe_by_name
119  // function depends on the
120  // particular format of the string
121  // this function returns, so they
122  // have to be kept in synch
123 
124  // note that this->degree is the maximal
125  // polynomial degree and is thus one higher
126  // than the argument given to the
127  // constructor
128  std::ostringstream namebuf;
129  namebuf << "FE_RaviartThomas<" << dim << ">(" << this->degree - 1 << ")";
130 
131  return namebuf.str();
132 }
133 
134 
135 template <int dim>
136 std::unique_ptr<FiniteElement<dim, dim>>
138 {
139  return std::make_unique<FE_RaviartThomas<dim>>(*this);
140 }
141 
142 
143 //---------------------------------------------------------------------------
144 // Auxiliary and internal functions
145 //---------------------------------------------------------------------------
146 
147 
148 template <int dim>
149 void
151 {
152  QGauss<dim> cell_quadrature(deg + 1);
153  const unsigned int n_interior_points = (deg > 0) ? cell_quadrature.size() : 0;
154 
155  unsigned int n_face_points = (dim > 1) ? 1 : 0;
156  // compute (deg+1)^(dim-1)
157  for (unsigned int d = 1; d < dim; ++d)
158  n_face_points *= deg + 1;
159 
160 
161  this->generalized_support_points.resize(
162  GeometryInfo<dim>::faces_per_cell * n_face_points + n_interior_points);
163  this->generalized_face_support_points.resize(n_face_points);
164 
165  // Number of the point being entered
166  unsigned int current = 0;
167 
168  if (dim > 1)
169  {
170  QGauss<dim - 1> face_points(deg + 1);
171  TensorProductPolynomials<dim - 1> legendre =
173 
174  boundary_weights.reinit(n_face_points, legendre.n());
175 
176  // Assert (face_points.size() == this->dofs_per_face,
177  // ExcInternalError());
178 
179  for (unsigned int k = 0; k < n_face_points; ++k)
180  {
181  this->generalized_face_support_points[k] = face_points.point(k);
182  // Compute its quadrature
183  // contribution for each
184  // moment.
185  for (unsigned int i = 0; i < legendre.n(); ++i)
186  {
187  boundary_weights(k, i) =
188  face_points.weight(k) *
189  legendre.compute_value(i, face_points.point(k));
190  }
191  }
192 
193  Quadrature<dim> faces =
195  for (; current < GeometryInfo<dim>::faces_per_cell * n_face_points;
196  ++current)
197  {
198  // Enter the support point
199  // into the vector
200  this->generalized_support_points[current] =
202  0, true, false, false, n_face_points));
203  }
204  }
205 
206  if (deg == 0)
207  return;
208 
209  // Create Legendre basis for the space D_xi Q_k
210  std::unique_ptr<AnisotropicPolynomials<dim>> polynomials[dim];
211  for (unsigned int dd = 0; dd < dim; ++dd)
212  {
213  std::vector<std::vector<Polynomials::Polynomial<double>>> poly(dim);
214  for (unsigned int d = 0; d < dim; ++d)
217 
218  polynomials[dd] = std::make_unique<AnisotropicPolynomials<dim>>(poly);
219  }
220 
222  TableIndices<3>(n_interior_points, polynomials[0]->n(), dim));
223 
224  for (unsigned int k = 0; k < cell_quadrature.size(); ++k)
225  {
226  this->generalized_support_points[current++] = cell_quadrature.point(k);
227  for (unsigned int i = 0; i < polynomials[0]->n(); ++i)
228  for (unsigned int d = 0; d < dim; ++d)
229  interior_weights(k, i, d) =
230  cell_quadrature.weight(k) *
231  polynomials[d]->compute_value(i, cell_quadrature.point(k));
232  }
233 
234  Assert(current == this->generalized_support_points.size(),
235  ExcInternalError());
236 }
237 
238 
239 
240 template <>
241 void
243 {
244  // there is only one refinement case in 1d,
245  // which is the isotropic one (first index of
246  // the matrix array has to be 0)
247  for (unsigned int i = 0; i < GeometryInfo<1>::max_children_per_cell; ++i)
248  this->restriction[0][i].reinit(0, 0);
249 }
250 
251 
252 
253 // This function is the same Raviart-Thomas interpolation performed by
254 // interpolate. Still, we cannot use interpolate, since it was written
255 // for smooth functions. The functions interpolated here are not
256 // smooth, maybe even not continuous. Therefore, we must double the
257 // number of quadrature points in each direction in order to integrate
258 // only smooth functions.
259 
260 // Then again, the interpolated function is chosen such that the
261 // moments coincide with the function to be interpolated.
262 
263 template <int dim>
264 void
266 {
267  const unsigned int iso = RefinementCase<dim>::isotropic_refinement - 1;
268 
269  QGauss<dim - 1> q_base(this->degree);
270  const unsigned int n_face_points = q_base.size();
271  // First, compute interpolation on
272  // subfaces
273  for (unsigned int face : GeometryInfo<dim>::face_indices())
274  {
275  // The shape functions of the
276  // child cell are evaluated
277  // in the quadrature points
278  // of a full face.
279  Quadrature<dim> q_face = QProjector<dim>::project_to_face(q_base, face);
280  // Store shape values, since the
281  // evaluation suffers if not
282  // ordered by point
283  Table<2, double> cached_values_on_face(this->dofs_per_cell,
284  q_face.size());
285  for (unsigned int k = 0; k < q_face.size(); ++k)
286  for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
287  cached_values_on_face(i, k) = this->shape_value_component(
288  i, q_face.point(k), GeometryInfo<dim>::unit_normal_direction[face]);
289 
290  for (unsigned int sub = 0; sub < GeometryInfo<dim>::max_children_per_face;
291  ++sub)
292  {
293  // The weight functions for
294  // the coarse face are
295  // evaluated on the subface
296  // only.
297  Quadrature<dim> q_sub =
298  QProjector<dim>::project_to_subface(q_base, face, sub);
299  const unsigned int child = GeometryInfo<dim>::child_cell_on_face(
301 
302  // On a certain face, we must
303  // compute the moments of ALL
304  // fine level functions with
305  // the coarse level weight
306  // functions belonging to
307  // that face. Due to the
308  // orthogonalization process
309  // when building the shape
310  // functions, these weights
311  // are equal to the
312  // corresponding shape
313  // functions.
314  for (unsigned int k = 0; k < n_face_points; ++k)
315  for (unsigned int i_child = 0; i_child < this->dofs_per_cell;
316  ++i_child)
317  for (unsigned int i_face = 0; i_face < this->dofs_per_face;
318  ++i_face)
319  {
320  // The quadrature
321  // weights on the
322  // subcell are NOT
323  // transformed, so we
324  // have to do it here.
325  this->restriction[iso][child](face * this->dofs_per_face +
326  i_face,
327  i_child) +=
328  Utilities::fixed_power<dim - 1>(.5) * q_sub.weight(k) *
329  cached_values_on_face(i_child, k) *
330  this->shape_value_component(
331  face * this->dofs_per_face + i_face,
332  q_sub.point(k),
334  }
335  }
336  }
337 
338  if (this->degree == 1)
339  return;
340 
341  // Create Legendre basis for the space D_xi Q_k. Here, we cannot
342  // use the shape functions
343  std::unique_ptr<AnisotropicPolynomials<dim>> polynomials[dim];
344  for (unsigned int dd = 0; dd < dim; ++dd)
345  {
346  std::vector<std::vector<Polynomials::Polynomial<double>>> poly(dim);
347  for (unsigned int d = 0; d < dim; ++d)
348  poly[d] =
350  poly[dd] =
352 
353  polynomials[dd] = std::make_unique<AnisotropicPolynomials<dim>>(poly);
354  }
355 
356  QGauss<dim> q_cell(this->degree);
357  const unsigned int start_cell_dofs =
359 
360  // Store shape values, since the
361  // evaluation suffers if not
362  // ordered by point
363  Table<3, double> cached_values_on_cell(this->dofs_per_cell,
364  q_cell.size(),
365  dim);
366  for (unsigned int k = 0; k < q_cell.size(); ++k)
367  for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
368  for (unsigned int d = 0; d < dim; ++d)
369  cached_values_on_cell(i, k, d) =
370  this->shape_value_component(i, q_cell.point(k), d);
371 
372  for (unsigned int child = 0; child < GeometryInfo<dim>::max_children_per_cell;
373  ++child)
374  {
375  Quadrature<dim> q_sub = QProjector<dim>::project_to_child(q_cell, child);
376 
377  for (unsigned int k = 0; k < q_sub.size(); ++k)
378  for (unsigned int i_child = 0; i_child < this->dofs_per_cell; ++i_child)
379  for (unsigned int d = 0; d < dim; ++d)
380  for (unsigned int i_weight = 0; i_weight < polynomials[d]->n();
381  ++i_weight)
382  {
383  this->restriction[iso][child](start_cell_dofs + i_weight * dim +
384  d,
385  i_child) +=
386  q_sub.weight(k) * cached_values_on_cell(i_child, k, d) *
387  polynomials[d]->compute_value(i_weight, q_sub.point(k));
388  }
389  }
390 }
391 
392 
393 
394 template <int dim>
395 std::vector<unsigned int>
397 {
398  // the element is face-based and we have
399  // (deg+1)^(dim-1) DoFs per face
400  unsigned int dofs_per_face = 1;
401  for (unsigned int d = 1; d < dim; ++d)
402  dofs_per_face *= deg + 1;
403 
404  // and then there are interior dofs
405  const unsigned int interior_dofs = dim * deg * dofs_per_face;
406 
407  std::vector<unsigned int> dpo(dim + 1);
408  dpo[dim - 1] = dofs_per_face;
409  dpo[dim] = interior_dofs;
410 
411  return dpo;
412 }
413 
414 
415 
416 template <int dim>
417 std::pair<Table<2, bool>, std::vector<unsigned int>>
419 {
420  Table<2, bool> constant_modes(dim, this->dofs_per_cell);
421  for (unsigned int d = 0; d < dim; ++d)
422  for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
423  constant_modes(d, i) = true;
424  std::vector<unsigned int> components;
425  for (unsigned int d = 0; d < dim; ++d)
426  components.push_back(d);
427  return std::pair<Table<2, bool>, std::vector<unsigned int>>(constant_modes,
428  components);
429 }
430 
431 
432 
433 //---------------------------------------------------------------------------
434 // Data field initialization
435 //---------------------------------------------------------------------------
436 
437 
438 template <int dim>
439 bool
440 FE_RaviartThomas<dim>::has_support_on_face(const unsigned int shape_index,
441  const unsigned int face_index) const
442 {
443  AssertIndexRange(shape_index, this->dofs_per_cell);
445 
446  // Return computed values if we
447  // know them easily. Otherwise, it
448  // is always safe to return true.
449  switch (this->degree)
450  {
451  case 1:
452  {
453  switch (dim)
454  {
455  case 2:
456  {
457  // only on the one
458  // non-adjacent face
459  // are the values
460  // actually zero. list
461  // these in a table
462  return (face_index !=
463  GeometryInfo<dim>::opposite_face[shape_index]);
464  }
465 
466  default:
467  return true;
468  }
469  }
470 
471  default: // other rt_order
472  return true;
473  }
474 
475  return true;
476 }
477 
478 
479 
480 template <int dim>
481 void
483  const std::vector<Vector<double>> &support_point_values,
484  std::vector<double> & nodal_values) const
485 {
486  Assert(support_point_values.size() == this->generalized_support_points.size(),
487  ExcDimensionMismatch(support_point_values.size(),
488  this->generalized_support_points.size()));
489  Assert(nodal_values.size() == this->dofs_per_cell,
490  ExcDimensionMismatch(nodal_values.size(), this->dofs_per_cell));
491  Assert(support_point_values[0].size() == this->n_components(),
492  ExcDimensionMismatch(support_point_values[0].size(),
493  this->n_components()));
494 
495  std::fill(nodal_values.begin(), nodal_values.end(), 0.);
496 
497  const unsigned int n_face_points = boundary_weights.size(0);
498  for (unsigned int face : GeometryInfo<dim>::face_indices())
499  for (unsigned int k = 0; k < n_face_points; ++k)
500  for (unsigned int i = 0; i < boundary_weights.size(1); ++i)
501  {
502  nodal_values[i + face * this->dofs_per_face] +=
503  boundary_weights(k, i) *
504  support_point_values[face * n_face_points + k](
506  }
507 
508  const unsigned int start_cell_dofs =
510  const unsigned int start_cell_points =
511  GeometryInfo<dim>::faces_per_cell * n_face_points;
512 
513  for (unsigned int k = 0; k < interior_weights.size(0); ++k)
514  for (unsigned int i = 0; i < interior_weights.size(1); ++i)
515  for (unsigned int d = 0; d < dim; ++d)
516  nodal_values[start_cell_dofs + i * dim + d] +=
517  interior_weights(k, i, d) *
518  support_point_values[k + start_cell_points](d);
519 }
520 
521 
522 
523 template <int dim>
524 std::size_t
526 {
527  Assert(false, ExcNotImplemented());
528  return 0;
529 }
530 
531 
532 
533 // explicit instantiations
534 #include "fe_raviart_thomas.inst"
535 
536 
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:621
static Quadrature< dim > project_to_child(const Quadrature< dim > &quadrature, const unsigned int child_no)
Definition: quadrature.cc:1064
std::vector< std::vector< FullMatrix< double > > > restriction
Definition: fe.h:2414
const unsigned int components
Definition: fe_base.h:292
std::vector< Point< dim > > generalized_support_points
Definition: fe.h:2465
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:2440
static void project_to_subface(const SubQuadrature &quadrature, const unsigned int face_no, const unsigned int subface_no, std::vector< Point< dim >> &q_points, const RefinementCase< dim - 1 > &ref_case=RefinementCase< dim - 1 >::isotropic_refinement)
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
Definition: polynomial.cc:877
#define AssertIndexRange(index, range)
Definition: exceptions.h:1649
const unsigned int degree
Definition: fe_base.h:298
const Point< dim > & point(const unsigned int i) const
void invert(const FullMatrix< number2 > &M)
STL namespace.
FullMatrix< double > inverse_node_matrix
FullMatrix< double > compute_node_matrix(const FiniteElement< dim, spacedim > &fe)
std::vector< Point< dim - 1 > > generalized_face_support_points
Definition: fe.h:2471
virtual std::size_t memory_consumption() const override
virtual std::string get_name() const override
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
friend class FE_RaviartThomas
std::vector< std::vector< FullMatrix< double > > > prolongation
Definition: fe.h:2428
void reinit(const TableIndices< N > &new_size, const bool omit_default_initialization=false)
#define Assert(cond, exc)
Definition: exceptions.h:1419
void initialize_support_points(const unsigned int rt_degree)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:362
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
unsigned int size() const
const unsigned int dofs_per_cell
Definition: fe_base.h:282
static Quadrature< dim > project_to_all_faces(const SubQuadrature &quadrature)
virtual double shape_value_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
virtual std::unique_ptr< FiniteElement< dim, dim > > clone() const override
static void project_to_face(const SubQuadrature &quadrature, const unsigned int face_no, std::vector< Point< dim >> &q_points)
size_type size(const unsigned int i) const
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
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const override
unsigned int n_components() const
static unsigned int child_cell_on_face(const RefinementCase< dim > &ref_case, const unsigned int face, const unsigned int subface, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false, const RefinementCase< dim - 1 > &face_refinement_case=RefinementCase< dim - 1 >::isotropic_refinement)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:361
const unsigned int dofs_per_face
Definition: fe_base.h:275
void reinit_restriction_and_prolongation_matrices(const bool isotropic_restriction_only=false, const bool isotropic_prolongation_only=false)
Definition: fe.cc:275
static ::ExceptionBase & ExcNotImplemented()
static std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
virtual std::pair< Table< 2, bool >, std::vector< unsigned int > > get_constant_modes() const override
Table< 3, double > interior_weights
double weight(const unsigned int i) const
Table< 2, double > boundary_weights
static ::ExceptionBase & ExcInternalError()
std::vector< MappingKind > mapping_kind