Reference documentation for deal.II version GIT d77e5ebb0a 2023-01-27 22:35: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_raviart_thomas.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2003 - 2021 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 
22 #include <deal.II/base/table.h>
23 #include <deal.II/base/utilities.h>
24 
26 
27 #include <deal.II/fe/fe.h>
29 #include <deal.II/fe/fe_tools.h>
30 #include <deal.II/fe/fe_values.h>
31 #include <deal.II/fe/mapping.h>
32 
33 #include <deal.II/grid/tria.h>
35 
36 #include <iostream>
37 #include <memory>
38 #include <sstream>
39 
40 
42 
43 
44 template <int dim>
46  : FE_PolyTensor<dim>(
47  PolynomialsRaviartThomas<dim>(deg + 1, deg),
49  dim,
50  deg + 1,
51  FiniteElementData<dim>::Hdiv),
52  std::vector<bool>(PolynomialsRaviartThomas<dim>::n_polynomials(deg + 1,
53  deg),
54  true),
55  std::vector<ComponentMask>(
56  PolynomialsRaviartThomas<dim>::n_polynomials(deg + 1, deg),
57  std::vector<bool>(dim, true)))
58 {
59  Assert(dim >= 2, ExcImpossibleInDim(dim));
60  const unsigned int n_dofs = this->n_dofs_per_cell();
61 
63  // First, initialize the
64  // generalized support points and
65  // quadrature weights, since they
66  // are required for interpolation.
68 
69  // Now compute the inverse node matrix, generating the correct
70  // basis functions from the raw ones. For a discussion of what
71  // exactly happens here, see FETools::compute_node_matrix.
73  this->inverse_node_matrix.reinit(n_dofs, n_dofs);
74  this->inverse_node_matrix.invert(M);
75  // From now on, the shape functions provided by FiniteElement::shape_value
76  // and similar functions will be the correct ones, not
77  // the raw shape functions from the polynomial space anymore.
78 
79  // Reinit the vectors of
80  // restriction and prolongation
81  // matrices to the right sizes.
82  // Restriction only for isotropic
83  // refinement
85  // Fill prolongation matrices with embedding operators
88 
89  // TODO: the implementation makes the assumption that all faces have the
90  // same number of dofs
91  AssertDimension(this->n_unique_faces(), 1);
92  const unsigned int face_no = 0;
93 
94  // TODO[TL]: for anisotropic refinement we will probably need a table of
95  // submatrices with an array for each refine case
97  for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_face; ++i)
98  face_embeddings[i].reinit(this->n_dofs_per_face(face_no),
99  this->n_dofs_per_face(face_no));
100  FETools::compute_face_embedding_matrices<dim, double>(*this,
101  face_embeddings,
102  0,
103  0);
104  this->interface_constraints.reinit((1 << (dim - 1)) *
105  this->n_dofs_per_face(face_no),
106  this->n_dofs_per_face(face_no));
107  unsigned int target_row = 0;
108  for (unsigned int d = 0; d < GeometryInfo<dim>::max_children_per_face; ++d)
109  for (unsigned int i = 0; i < face_embeddings[d].m(); ++i)
110  {
111  for (unsigned int j = 0; j < face_embeddings[d].n(); ++j)
112  this->interface_constraints(target_row, j) = face_embeddings[d](i, j);
113  ++target_row;
114  }
115 
116  // We need to initialize the dof permutation table and the one for the sign
117  // change.
119 }
120 
121 
122 
123 template <int dim>
124 std::string
126 {
127  // note that the
128  // FETools::get_fe_by_name
129  // function depends on the
130  // particular format of the string
131  // this function returns, so they
132  // have to be kept in synch
133 
134  // note that this->degree is the maximal
135  // polynomial degree and is thus one higher
136  // than the argument given to the
137  // constructor
138  std::ostringstream namebuf;
139  namebuf << "FE_RaviartThomas<" << dim << ">(" << this->degree - 1 << ")";
140 
141  return namebuf.str();
142 }
143 
144 
145 template <int dim>
146 std::unique_ptr<FiniteElement<dim, dim>>
148 {
149  return std::make_unique<FE_RaviartThomas<dim>>(*this);
150 }
151 
152 
153 //---------------------------------------------------------------------------
154 // Auxiliary and internal functions
155 //---------------------------------------------------------------------------
156 
157 
158 template <int dim>
159 void
161 {
162  QGauss<dim> cell_quadrature(deg + 1);
163  const unsigned int n_interior_points = (deg > 0) ? cell_quadrature.size() : 0;
164 
165  // TODO: the implementation makes the assumption that all faces have the
166  // same number of dofs
167  AssertDimension(this->n_unique_faces(), 1);
168  const unsigned int face_no = 0;
169 
170  unsigned int n_face_points = (dim > 1) ? 1 : 0;
171  // compute (deg+1)^(dim-1)
172  for (unsigned int d = 1; d < dim; ++d)
173  n_face_points *= deg + 1;
174 
175 
176  this->generalized_support_points.resize(
177  GeometryInfo<dim>::faces_per_cell * n_face_points + n_interior_points);
178  this->generalized_face_support_points[face_no].resize(n_face_points);
179 
180  // Number of the point being entered
181  unsigned int current = 0;
182 
183  if (dim > 1)
184  {
185  QGauss<dim - 1> face_points(deg + 1);
188 
189  boundary_weights.reinit(n_face_points, legendre.n());
190 
191  for (unsigned int k = 0; k < n_face_points; ++k)
192  {
193  this->generalized_face_support_points[face_no][k] =
194  face_points.point(k);
195  // Compute its quadrature
196  // contribution for each
197  // moment.
198  for (unsigned int i = 0; i < legendre.n(); ++i)
199  {
200  boundary_weights(k, i) =
201  face_points.weight(k) *
202  legendre.compute_value(i, face_points.point(k));
203  }
204  }
205 
206  Quadrature<dim> faces =
208  face_points);
209  for (; current < GeometryInfo<dim>::faces_per_cell * n_face_points;
210  ++current)
211  {
212  // Enter the support point
213  // into the vector
214  this->generalized_support_points[current] = faces.point(
215  current +
217  this->reference_cell(), 0, true, false, false, n_face_points));
218  }
219  }
220 
221  if (deg == 0)
222  return;
223 
224  // Create Legendre basis for the space D_xi Q_k
225  std::unique_ptr<AnisotropicPolynomials<dim>> polynomials[dim];
226  for (unsigned int dd = 0; dd < dim; ++dd)
227  {
228  std::vector<std::vector<Polynomials::Polynomial<double>>> poly(dim);
229  for (unsigned int d = 0; d < dim; ++d)
232 
233  polynomials[dd] = std::make_unique<AnisotropicPolynomials<dim>>(poly);
234  }
235 
236  interior_weights.reinit(
237  TableIndices<3>(n_interior_points, polynomials[0]->n(), dim));
238 
239  for (unsigned int k = 0; k < cell_quadrature.size(); ++k)
240  {
241  this->generalized_support_points[current++] = cell_quadrature.point(k);
242  for (unsigned int i = 0; i < polynomials[0]->n(); ++i)
243  for (unsigned int d = 0; d < dim; ++d)
244  interior_weights(k, i, d) =
245  cell_quadrature.weight(k) *
246  polynomials[d]->compute_value(i, cell_quadrature.point(k));
247  }
248 
249  Assert(current == this->generalized_support_points.size(),
250  ExcInternalError());
251 }
252 
253 
254 template <int dim>
255 void
257 {
258  // For 1D do nothing.
259  //
260  // TODO: For 2D we simply keep the legacy behavior for now. This should be
261  // changed in the future and can be taken care of by similar means as the 3D
262  // case below. The legacy behavior can be found in fe_poly_tensor.cc in the
263  // function internal::FE_PolyTensor::get_dof_sign_change_h_div(...)
264  if (dim < 3)
265  return;
266 
267  // TODO: the implementation makes the assumption that all faces have the
268  // same number of dofs
269  AssertDimension(this->n_unique_faces(), 1);
270  const unsigned int face_no = 0;
271 
272  Assert(
273  this->adjust_quad_dof_index_for_face_orientation_table[0].n_elements() ==
274  this->reference_cell().n_face_orientations(face_no) *
275  this->n_dofs_per_quad(face_no),
276  ExcInternalError());
277 
278  // The 3D RaviartThomas space has tensor_degree*tensor_degree face dofs
279  const unsigned int n = this->tensor_degree();
280  Assert(n * n == this->n_dofs_per_quad(face_no), ExcInternalError());
281 
282  // The vector of tables adjust_quad_dof_index_for_face_orientation_table
283  // contains offsets for local face_dofs and is being filled with zeros in
284  // fe.cc. We need to fill it with the correct values in case of non-standard,
285  // flipped (rotated by +180 degrees) or rotated (rotated by +90 degrees) faces
286 
287  // The dofs on a face are connected to a n x n matrix. for example, for
288  // tensor_degree==3 we have the following dofs on a quad:
289 
290  // ___________
291  // | |
292  // | 6 7 8 |
293  // | |
294  // | 3 4 5 |
295  // | |
296  // | 0 1 2 |
297  // |___________|
298  //
299  // We have dof_index=i+n*j with index i in x-direction and index j in
300  // y-direction running from 0 to n-1. to extract i and j we can use
301  // i=dof_index%n and j=dof_index/n. The indices i and j can then be used to
302  // compute the offset.
303 
304  // Example: if the switches are (true | true | true) that means we rotate the
305  // face first by + 90 degree(counterclockwise) then by another +180
306  // degrees but we do not flip it since the face has standard
307  // orientation. The flip axis is the diagonal from the lower left to the upper
308  // right corner of the face. With these flags the configuration above becomes:
309  // ___________
310  // | |
311  // | 2 5 8 |
312  // | |
313  // | 1 4 7 |
314  // | |
315  // | 0 3 6 |
316  // |___________|
317  //
318  // This is exactly what is realized by the formulas implemented below. Note
319  // that the necessity of a permuattion depends on the three flags.
320 
321  // There is also a pattern for the sign change of the mapped face_dofs
322  // depending on the switches. In the above example it would be
323  // ___________
324  // | |
325  // | + - + |
326  // | |
327  // | + - + |
328  // | |
329  // | + - + |
330  // |___________|
331  //
332 
333  for (unsigned int local_face_dof = 0;
334  local_face_dof < this->n_dofs_per_quad(face_no);
335  ++local_face_dof)
336  {
337  // Row and column
338  unsigned int i = local_face_dof % n;
339  unsigned int j = local_face_dof / n;
340 
341  // We have 8 cases that are all treated the same way. Note that the
342  // corresponding case to case_no is just its binary representation.
343  // The above example of (false | true | true) would be case_no=3
344  for (unsigned int case_no = 0; case_no < 8; ++case_no)
345  {
346  // Get the binary representation of the case
347  const bool face_orientation = Utilities::get_bit(case_no, 2);
348  const bool face_flip = Utilities::get_bit(case_no, 1);
349  const bool face_rotation = Utilities::get_bit(case_no, 0);
350 
351  if (((!face_orientation) && (!face_rotation)) ||
352  ((face_orientation) && (face_rotation)))
353  {
354  // We flip across the diagonal
355  // This is the local face dof offset
356  this->adjust_quad_dof_index_for_face_orientation_table[face_no](
357  local_face_dof, case_no) = j + i * n - local_face_dof;
358  }
359  else
360  {
361  // Offset is zero
362  this->adjust_quad_dof_index_for_face_orientation_table[face_no](
363  local_face_dof, case_no) = 0;
364  } // if face needs dof permutation
365 
366  // Get new local face_dof by adding offset
367  const unsigned int new_local_face_dof =
368  local_face_dof +
369  this->adjust_quad_dof_index_for_face_orientation_table[face_no](
370  local_face_dof, case_no);
371  // compute new row and column index
372  i = new_local_face_dof % n;
373  j = new_local_face_dof / n;
374 
375  /*
376  * Now compute if a sign change is necessary. This is done for the
377  * case of face_orientation==true
378  */
379  // flip = false, rotation=true
380  if (!face_flip && face_rotation)
381  {
382  this->adjust_quad_dof_sign_for_face_orientation_table[face_no](
383  local_face_dof, case_no) = ((j % 2) == 1);
384  }
385  // flip = true, rotation=false
386  else if (face_flip && !face_rotation)
387  {
388  // This case is symmetric (although row and column may be
389  // switched)
390  this->adjust_quad_dof_sign_for_face_orientation_table[face_no](
391  local_face_dof, case_no) = ((j % 2) == 1) != ((i % 2) == 1);
392  }
393  // flip = true, rotation=true
394  else if (face_flip && face_rotation)
395  {
396  this->adjust_quad_dof_sign_for_face_orientation_table[face_no](
397  local_face_dof, case_no) = ((i % 2) == 1);
398  }
399  /*
400  * flip = false, rotation=false => nothing to do
401  */
402 
403  /*
404  * If face_orientation==false the sign flip is exactly the opposite.
405  */
406  if (!face_orientation)
407  this->adjust_quad_dof_sign_for_face_orientation_table[face_no](
408  local_face_dof, case_no) =
409  !this->adjust_quad_dof_sign_for_face_orientation_table[face_no](
410  local_face_dof, case_no);
411  } // case_no
412  } // local_face_dof
413 }
414 
415 
416 template <>
417 void
419 {
420  // there is only one refinement case in 1d,
421  // which is the isotropic one (first index of
422  // the matrix array has to be 0)
423  for (unsigned int i = 0; i < GeometryInfo<1>::max_children_per_cell; ++i)
424  this->restriction[0][i].reinit(0, 0);
425 }
426 
427 
428 
429 // This function is the same Raviart-Thomas interpolation performed by
430 // interpolate. Still, we cannot use interpolate, since it was written
431 // for smooth functions. The functions interpolated here are not
432 // smooth, maybe even not continuous. Therefore, we must double the
433 // number of quadrature points in each direction in order to integrate
434 // only smooth functions.
435 
436 // Then again, the interpolated function is chosen such that the
437 // moments coincide with the function to be interpolated.
438 
439 template <int dim>
440 void
442 {
443  const unsigned int iso = RefinementCase<dim>::isotropic_refinement - 1;
444 
445  QGauss<dim - 1> q_base(this->degree);
446  const unsigned int n_face_points = q_base.size();
447  // First, compute interpolation on
448  // subfaces
449  for (unsigned int face : GeometryInfo<dim>::face_indices())
450  {
451  // The shape functions of the
452  // child cell are evaluated
453  // in the quadrature points
454  // of a full face.
455  Quadrature<dim> q_face =
456  QProjector<dim>::project_to_face(this->reference_cell(), q_base, face);
457  // Store shape values, since the
458  // evaluation suffers if not
459  // ordered by point
460  Table<2, double> cached_values_on_face(this->n_dofs_per_cell(),
461  q_face.size());
462  for (unsigned int k = 0; k < q_face.size(); ++k)
463  for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
464  cached_values_on_face(i, k) = this->shape_value_component(
465  i, q_face.point(k), GeometryInfo<dim>::unit_normal_direction[face]);
466 
467  for (unsigned int sub = 0; sub < GeometryInfo<dim>::max_children_per_face;
468  ++sub)
469  {
470  // The weight functions for
471  // the coarse face are
472  // evaluated on the subface
473  // only.
475  this->reference_cell(), q_base, face, sub);
476  const unsigned int child = GeometryInfo<dim>::child_cell_on_face(
478 
479  // On a certain face, we must
480  // compute the moments of ALL
481  // fine level functions with
482  // the coarse level weight
483  // functions belonging to
484  // that face. Due to the
485  // orthogonalization process
486  // when building the shape
487  // functions, these weights
488  // are equal to the
489  // corresponding shape
490  // functions.
491  for (unsigned int k = 0; k < n_face_points; ++k)
492  for (unsigned int i_child = 0; i_child < this->n_dofs_per_cell();
493  ++i_child)
494  for (unsigned int i_face = 0;
495  i_face < this->n_dofs_per_face(face);
496  ++i_face)
497  {
498  // The quadrature
499  // weights on the
500  // subcell are NOT
501  // transformed, so we
502  // have to do it here.
503  this->restriction[iso][child](
504  face * this->n_dofs_per_face(face) + i_face, i_child) +=
505  Utilities::fixed_power<dim - 1>(.5) * q_sub.weight(k) *
506  cached_values_on_face(i_child, k) *
507  this->shape_value_component(
508  face * this->n_dofs_per_face(face) + i_face,
509  q_sub.point(k),
511  }
512  }
513  }
514 
515  if (this->degree == 1)
516  return;
517 
518  // Create Legendre basis for the space D_xi Q_k. Here, we cannot
519  // use the shape functions
520  std::unique_ptr<AnisotropicPolynomials<dim>> polynomials[dim];
521  for (unsigned int dd = 0; dd < dim; ++dd)
522  {
523  std::vector<std::vector<Polynomials::Polynomial<double>>> poly(dim);
524  for (unsigned int d = 0; d < dim; ++d)
525  poly[d] =
527  poly[dd] =
529 
530  polynomials[dd] = std::make_unique<AnisotropicPolynomials<dim>>(poly);
531  }
532 
533  // TODO: the implementation makes the assumption that all faces have the
534  // same number of dofs
535  AssertDimension(this->n_unique_faces(), 1);
536  const unsigned int face_no = 0;
537 
538  QGauss<dim> q_cell(this->degree);
539  const unsigned int start_cell_dofs =
540  GeometryInfo<dim>::faces_per_cell * this->n_dofs_per_face(face_no);
541 
542  // Store shape values, since the
543  // evaluation suffers if not
544  // ordered by point
545  Table<3, double> cached_values_on_cell(this->n_dofs_per_cell(),
546  q_cell.size(),
547  dim);
548  for (unsigned int k = 0; k < q_cell.size(); ++k)
549  for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
550  for (unsigned int d = 0; d < dim; ++d)
551  cached_values_on_cell(i, k, d) =
552  this->shape_value_component(i, q_cell.point(k), d);
553 
554  for (unsigned int child = 0; child < GeometryInfo<dim>::max_children_per_cell;
555  ++child)
556  {
557  Quadrature<dim> q_sub =
559  q_cell,
560  child);
561 
562  for (unsigned int k = 0; k < q_sub.size(); ++k)
563  for (unsigned int i_child = 0; i_child < this->n_dofs_per_cell();
564  ++i_child)
565  for (unsigned int d = 0; d < dim; ++d)
566  for (unsigned int i_weight = 0; i_weight < polynomials[d]->n();
567  ++i_weight)
568  {
569  this->restriction[iso][child](start_cell_dofs + i_weight * dim +
570  d,
571  i_child) +=
572  q_sub.weight(k) * cached_values_on_cell(i_child, k, d) *
573  polynomials[d]->compute_value(i_weight, q_sub.point(k));
574  }
575  }
576 }
577 
578 
579 
580 template <int dim>
581 std::vector<unsigned int>
583 {
584  // the element is face-based and we have
585  // (deg+1)^(dim-1) DoFs per face
586  unsigned int dofs_per_face = 1;
587  for (unsigned int d = 1; d < dim; ++d)
588  dofs_per_face *= deg + 1;
589 
590  // and then there are interior dofs
591  const unsigned int interior_dofs = dim * deg * dofs_per_face;
592 
593  std::vector<unsigned int> dpo(dim + 1);
594  dpo[dim - 1] = dofs_per_face;
595  dpo[dim] = interior_dofs;
596 
597  return dpo;
598 }
599 
600 
601 
602 template <int dim>
603 std::pair<Table<2, bool>, std::vector<unsigned int>>
605 {
606  Table<2, bool> constant_modes(dim, this->n_dofs_per_cell());
607  for (unsigned int d = 0; d < dim; ++d)
608  for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
609  constant_modes(d, i) = true;
610  std::vector<unsigned int> components;
611  for (unsigned int d = 0; d < dim; ++d)
612  components.push_back(d);
613  return std::pair<Table<2, bool>, std::vector<unsigned int>>(constant_modes,
614  components);
615 }
616 
617 
618 
619 //---------------------------------------------------------------------------
620 // Data field initialization
621 //---------------------------------------------------------------------------
622 
623 
624 template <int dim>
625 bool
626 FE_RaviartThomas<dim>::has_support_on_face(const unsigned int shape_index,
627  const unsigned int face_index) const
628 {
629  AssertIndexRange(shape_index, this->n_dofs_per_cell());
631 
632  // Return computed values if we
633  // know them easily. Otherwise, it
634  // is always safe to return true.
635  switch (this->degree)
636  {
637  case 1:
638  {
639  switch (dim)
640  {
641  case 2:
642  {
643  // only on the one
644  // non-adjacent face
645  // are the values
646  // actually zero. list
647  // these in a table
648  return (face_index !=
649  GeometryInfo<dim>::opposite_face[shape_index]);
650  }
651 
652  default:
653  return true;
654  }
655  }
656 
657  default: // other rt_order
658  return true;
659  }
660 
661  return true;
662 }
663 
664 
665 
666 template <int dim>
667 void
669  const std::vector<Vector<double>> &support_point_values,
670  std::vector<double> & nodal_values) const
671 {
672  Assert(support_point_values.size() == this->generalized_support_points.size(),
673  ExcDimensionMismatch(support_point_values.size(),
674  this->generalized_support_points.size()));
675  Assert(nodal_values.size() == this->n_dofs_per_cell(),
676  ExcDimensionMismatch(nodal_values.size(), this->n_dofs_per_cell()));
677  Assert(support_point_values[0].size() == this->n_components(),
678  ExcDimensionMismatch(support_point_values[0].size(),
679  this->n_components()));
680 
681  std::fill(nodal_values.begin(), nodal_values.end(), 0.);
682 
683  const unsigned int n_face_points = boundary_weights.size(0);
684  for (unsigned int face : GeometryInfo<dim>::face_indices())
685  for (unsigned int k = 0; k < n_face_points; ++k)
686  for (unsigned int i = 0; i < boundary_weights.size(1); ++i)
687  {
688  nodal_values[i + face * this->n_dofs_per_face(face)] +=
689  boundary_weights(k, i) *
690  support_point_values[face * n_face_points + k](
692  }
693 
694  // TODO: the implementation makes the assumption that all faces have the
695  // same number of dofs
696  AssertDimension(this->n_unique_faces(), 1);
697  const unsigned int face_no = 0;
698 
699  const unsigned int start_cell_dofs =
700  GeometryInfo<dim>::faces_per_cell * this->n_dofs_per_face(face_no);
701  const unsigned int start_cell_points =
702  GeometryInfo<dim>::faces_per_cell * n_face_points;
703 
704  for (unsigned int k = 0; k < interior_weights.size(0); ++k)
705  for (unsigned int i = 0; i < interior_weights.size(1); ++i)
706  for (unsigned int d = 0; d < dim; ++d)
707  nodal_values[start_cell_dofs + i * dim + d] +=
708  interior_weights(k, i, d) *
709  support_point_values[k + start_cell_points](d);
710 }
711 
712 
713 
714 template <int dim>
715 std::size_t
717 {
718  Assert(false, ExcNotImplemented());
719  return 0;
720 }
721 
722 
723 
724 // explicit instantiations
725 #include "fe_raviart_thomas.inst"
726 
727 
FullMatrix< double > inverse_node_matrix
std::vector< MappingKind > mapping_kind
virtual std::size_t memory_consumption() const override
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
virtual std::pair< Table< 2, bool >, std::vector< unsigned int > > get_constant_modes() const override
void initialize_quad_dof_index_permutation_and_sign_change()
virtual std::unique_ptr< FiniteElement< dim, dim > > clone() const override
friend class FE_RaviartThomas
void initialize_support_points(const unsigned int rt_degree)
static std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
virtual std::string get_name() const override
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
void reinit_restriction_and_prolongation_matrices(const bool isotropic_restriction_only=false, const bool isotropic_prolongation_only=false)
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 std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
Definition: polynomial.cc:745
static void project_to_subface(const ReferenceCell &reference_cell, 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 Quadrature< dim > project_to_child(const ReferenceCell &reference_cell, const Quadrature< dim > &quadrature, const unsigned int child_no)
Definition: qprojector.cc:1132
static Quadrature< dim > project_to_all_faces(const ReferenceCell &reference_cell, const hp::QCollection< dim - 1 > &quadrature)
static void project_to_face(const ReferenceCell &reference_cell, const SubQuadrature &quadrature, const unsigned int face_no, std::vector< Point< dim >> &q_points)
double weight(const unsigned int i) const
const Point< dim > & point(const unsigned int i) const
unsigned int size() const
Definition: vector.h:109
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:461
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:462
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1509
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1703
#define AssertIndexRange(index, range)
Definition: exceptions.h:1768
@ mapping_raviart_thomas
Definition: mapping.h:128
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)
bool get_bit(const unsigned char number, const unsigned int n)
Definition: utilities.h:1549
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618
double legendre(unsigned int l, double x)
Definition: cmath.h:75
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)