Reference documentation for deal.II version GIT b6bf1e606d 2022-08-11 15: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_raviart_thomas_nodal.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2003 - 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 
21 
23 
24 #include <deal.II/fe/fe.h>
25 #include <deal.II/fe/fe_nothing.h>
27 #include <deal.II/fe/fe_tools.h>
28 #include <deal.II/fe/mapping.h>
29 
30 #include <memory>
31 #include <sstream>
32 
33 
35 
36 
37 // ---------------- polynomial class for FE_RaviartThomasNodal ---------------
38 
39 namespace
40 {
41  // Return a vector of "dofs per object" where the components of the returned
42  // vector refer to:
43  // 0 = vertex
44  // 1 = edge
45  // 2 = face (which is a cell in 2D)
46  // 3 = cell
47  std::vector<unsigned int>
48  get_rt_dpo_vector(const unsigned int dim, const unsigned int degree)
49  {
50  std::vector<unsigned int> dpo(dim + 1);
51  dpo[0] = 0;
52  dpo[1] = 0;
53  unsigned int dofs_per_face = 1;
54  for (unsigned int d = 1; d < dim; ++d)
55  dofs_per_face *= (degree + 1);
56 
57  dpo[dim - 1] = dofs_per_face;
58  dpo[dim] = dim * degree * dofs_per_face;
59 
60  return dpo;
61  }
62 } // namespace
63 
64 
65 
66 // --------------------- actual implementation of element --------------------
67 
68 template <int dim>
70  : FE_PolyTensor<dim>(
71  PolynomialsRaviartThomas<dim>(degree + 1, degree),
72  FiniteElementData<dim>(get_rt_dpo_vector(dim, degree),
73  dim,
74  degree + 1,
75  FiniteElementData<dim>::Hdiv),
76  std::vector<bool>(1, false),
77  std::vector<ComponentMask>(
78  PolynomialsRaviartThomas<dim>::n_polynomials(degree + 1, degree),
79  std::vector<bool>(dim, true)))
80 {
81  Assert(dim >= 2, ExcImpossibleInDim(dim));
82 
84 
85  // First, initialize the generalized support points and quadrature weights,
86  // since they are required for interpolation.
91  this->n_dofs_per_cell());
92 
93  const unsigned int face_no = 0;
94  if (dim > 1)
95  this->generalized_face_support_points[face_no] =
96  degree == 0 ? QGauss<dim - 1>(1).get_points() :
97  QGaussLobatto<dim - 1>(degree + 1).get_points();
98 
100  for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_face; ++i)
101  face_embeddings[i].reinit(this->n_dofs_per_face(face_no),
102  this->n_dofs_per_face(face_no));
103  FETools::compute_face_embedding_matrices<dim, double>(*this,
104  face_embeddings,
105  0,
106  0);
108  this->n_dofs_per_face(face_no),
109  this->n_dofs_per_face(face_no));
110  unsigned int target_row = 0;
111  for (unsigned int d = 0; d < GeometryInfo<dim>::max_children_per_face; ++d)
112  for (unsigned int i = 0; i < face_embeddings[d].m(); ++i)
113  {
114  for (unsigned int j = 0; j < face_embeddings[d].n(); ++j)
115  this->interface_constraints(target_row, j) = face_embeddings[d](i, j);
116  ++target_row;
117  }
118 
119  // We need to initialize the dof permutation table and the one for the sign
120  // change.
122 }
123 
124 
125 
126 template <int dim>
127 std::string
129 {
130  // note that the FETools::get_fe_by_name function depends on the particular
131  // format of the string this function returns, so they have to be kept in
132  // synch
133 
134  // note that this->degree is the maximal polynomial degree and is thus one
135  // higher than the argument given to the constructor
136  return "FE_RaviartThomasNodal<" + std::to_string(dim) + ">(" +
137  std::to_string(this->degree - 1) + ")";
138 }
139 
140 
141 template <int dim>
142 std::unique_ptr<FiniteElement<dim, dim>>
144 {
145  return std::make_unique<FE_RaviartThomasNodal<dim>>(*this);
146 }
147 
148 
149 //---------------------------------------------------------------------------
150 // Auxiliary and internal functions
151 //---------------------------------------------------------------------------
152 
153 
154 
155 template <int dim>
156 void
158  dim>::initialize_quad_dof_index_permutation_and_sign_change()
159 {
160  // for 1D and 2D, do nothing
161  if (dim < 3)
162  return;
163 
164  const unsigned int n = this->degree;
165  const unsigned int face_no = 0;
166  Assert(n * n == this->n_dofs_per_quad(face_no), ExcInternalError());
167  for (unsigned int local = 0; local < this->n_dofs_per_quad(face_no); ++local)
168  // face support points are in lexicographic ordering with x running
169  // fastest. invert that (y running fastest)
170  {
171  unsigned int i = local % n, j = local / n;
172 
173  // face_orientation=false, face_flip=false, face_rotation=false
174  this->adjust_quad_dof_index_for_face_orientation_table[face_no](local,
175  0) =
176  j + i * n - local;
177  // face_orientation=false, face_flip=false, face_rotation=true
178  this->adjust_quad_dof_index_for_face_orientation_table[face_no](local,
179  1) =
180  i + (n - 1 - j) * n - local;
181  // face_orientation=false, face_flip=true, face_rotation=false
182  this->adjust_quad_dof_index_for_face_orientation_table[face_no](local,
183  2) =
184  (n - 1 - j) + (n - 1 - i) * n - local;
185  // face_orientation=false, face_flip=true, face_rotation=true
186  this->adjust_quad_dof_index_for_face_orientation_table[face_no](local,
187  3) =
188  (n - 1 - i) + j * n - local;
189  // face_orientation=true, face_flip=false, face_rotation=false
190  this->adjust_quad_dof_index_for_face_orientation_table[face_no](local,
191  4) = 0;
192  // face_orientation=true, face_flip=false, face_rotation=true
193  this->adjust_quad_dof_index_for_face_orientation_table[face_no](local,
194  5) =
195  j + (n - 1 - i) * n - local;
196  // face_orientation=true, face_flip=true, face_rotation=false
197  this->adjust_quad_dof_index_for_face_orientation_table[face_no](local,
198  6) =
199  (n - 1 - i) + (n - 1 - j) * n - local;
200  // face_orientation=true, face_flip=true, face_rotation=true
201  this->adjust_quad_dof_index_for_face_orientation_table[face_no](local,
202  7) =
203  (n - 1 - j) + i * n - local;
204 
205  // for face_orientation == false, we need to switch the sign
206  for (unsigned int i = 0; i < 4; ++i)
207  this->adjust_quad_dof_sign_for_face_orientation_table[face_no](local,
208  i) = 1;
209  }
210 }
211 
212 
213 
214 template <int dim>
215 bool
217  const unsigned int shape_index,
218  const unsigned int face_index) const
219 {
220  AssertIndexRange(shape_index, this->n_dofs_per_cell());
222 
223  // The first degrees of freedom are on the faces and each face has degree
224  // degrees.
225  const unsigned int support_face = shape_index / this->n_dofs_per_face();
226 
227  // The only thing we know for sure is that shape functions with support on
228  // one face are zero on the opposite face.
229  if (support_face < GeometryInfo<dim>::faces_per_cell)
230  return (face_index != GeometryInfo<dim>::opposite_face[support_face]);
231 
232  // In all other cases, return true, which is safe
233  return true;
234 }
235 
236 
237 
238 template <int dim>
239 void
242  const std::vector<Vector<double>> &support_point_values,
243  std::vector<double> & nodal_values) const
244 {
245  Assert(support_point_values.size() == this->generalized_support_points.size(),
246  ExcDimensionMismatch(support_point_values.size(),
247  this->generalized_support_points.size()));
248  Assert(nodal_values.size() == this->n_dofs_per_cell(),
249  ExcDimensionMismatch(nodal_values.size(), this->n_dofs_per_cell()));
250  Assert(support_point_values[0].size() == this->n_components(),
251  ExcDimensionMismatch(support_point_values[0].size(),
252  this->n_components()));
253 
254  // First do interpolation on faces. There, the component evaluated depends
255  // on the face direction and orientation.
256  unsigned int fbase = 0;
257  unsigned int f = 0;
258  for (; f < GeometryInfo<dim>::faces_per_cell;
259  ++f, fbase += this->n_dofs_per_face(f))
260  {
261  for (unsigned int i = 0; i < this->n_dofs_per_face(f); ++i)
262  {
263  nodal_values[fbase + i] = support_point_values[fbase + i](
265  }
266  }
267 
268  // The remaining points form dim chunks, one for each component
269  const unsigned int istep = (this->n_dofs_per_cell() - fbase) / dim;
270  Assert((this->n_dofs_per_cell() - fbase) % dim == 0, ExcInternalError());
271 
272  f = 0;
273  while (fbase < this->n_dofs_per_cell())
274  {
275  for (unsigned int i = 0; i < istep; ++i)
276  {
277  nodal_values[fbase + i] = support_point_values[fbase + i](f);
278  }
279  fbase += istep;
280  ++f;
281  }
282  Assert(fbase == this->n_dofs_per_cell(), ExcInternalError());
283 }
284 
285 
286 
287 // TODO: There are tests that check that the following few functions don't
288 // produce assertion failures, but none that actually check whether they do the
289 // right thing. one example for such a test would be to project a function onto
290 // an hp-space and make sure that the convergence order is correct with regard
291 // to the lowest used polynomial degree
292 
293 template <int dim>
294 bool
296 {
297  return true;
298 }
299 
300 
301 template <int dim>
302 std::vector<std::pair<unsigned int, unsigned int>>
304  const FiniteElement<dim> &fe_other) const
305 {
306  // we can presently only compute these identities if both FEs are
307  // FE_RaviartThomasNodals or the other is FE_Nothing. In either case, no
308  // dofs are assigned on the vertex, so we shouldn't be getting here at all.
309  if (dynamic_cast<const FE_RaviartThomasNodal<dim> *>(&fe_other) != nullptr)
310  return std::vector<std::pair<unsigned int, unsigned int>>();
311  else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
312  return std::vector<std::pair<unsigned int, unsigned int>>();
313  else
314  {
315  Assert(false, ExcNotImplemented());
316  return std::vector<std::pair<unsigned int, unsigned int>>();
317  }
318 }
319 
320 
321 
322 template <int dim>
323 std::vector<std::pair<unsigned int, unsigned int>>
325  const FiniteElement<dim> &fe_other) const
326 {
327  // we can presently only compute these identities if both FEs are
328  // FE_RaviartThomasNodals or if the other one is FE_Nothing
329  if (const FE_RaviartThomasNodal<dim> *fe_q_other =
330  dynamic_cast<const FE_RaviartThomasNodal<dim> *>(&fe_other))
331  {
332  // dofs are located on faces; these are only lines in 2d
333  if (dim != 2)
334  return std::vector<std::pair<unsigned int, unsigned int>>();
335 
336  // dofs are located along lines, so two dofs are identical only if in
337  // the following two cases (remember that the face support points are
338  // Gauss points):
339  // 1. this->degree = fe_q_other->degree,
340  // in the case, all the dofs on the line are identical
341  // 2. this->degree-1 and fe_q_other->degree-1
342  // are both even, i.e. this->dof_per_line and fe_q_other->dof_per_line
343  // are both odd, there exists only one point (the middle one) such
344  // that dofs are identical on this point
345  //
346  // to understand this, note that this->degree is the *maximal*
347  // polynomial degree, and is thus one higher than the argument given to
348  // the constructor
349  const unsigned int p = this->degree - 1;
350  const unsigned int q = fe_q_other->degree - 1;
351 
352  std::vector<std::pair<unsigned int, unsigned int>> identities;
353 
354  if (p == q)
355  for (unsigned int i = 0; i < p + 1; ++i)
356  identities.emplace_back(i, i);
357 
358  else if (p % 2 == 0 && q % 2 == 0)
359  identities.emplace_back(p / 2, q / 2);
360 
361  return identities;
362  }
363  else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
364  {
365  // the FE_Nothing has no degrees of freedom, so there are no
366  // equivalencies to be recorded
367  return std::vector<std::pair<unsigned int, unsigned int>>();
368  }
369  else
370  {
371  Assert(false, ExcNotImplemented());
372  return std::vector<std::pair<unsigned int, unsigned int>>();
373  }
374 }
375 
376 
377 template <int dim>
378 std::vector<std::pair<unsigned int, unsigned int>>
380  const FiniteElement<dim> &fe_other,
381  const unsigned int face_no) const
382 {
383  // we can presently only compute these identities if both FEs are
384  // FE_RaviartThomasNodals or if the other one is FE_Nothing
385  if (const FE_RaviartThomasNodal<dim> *fe_q_other =
386  dynamic_cast<const FE_RaviartThomasNodal<dim> *>(&fe_other))
387  {
388  // dofs are located on faces; these are only quads in 3d
389  if (dim != 3)
390  return std::vector<std::pair<unsigned int, unsigned int>>();
391 
392  // this works exactly like the line case above
393  const unsigned int p = this->n_dofs_per_quad(face_no);
394 
395  AssertDimension(fe_q_other->n_unique_faces(), 1);
396  const unsigned int q = fe_q_other->n_dofs_per_quad(0);
397 
398  std::vector<std::pair<unsigned int, unsigned int>> identities;
399 
400  if (p == q)
401  for (unsigned int i = 0; i < p; ++i)
402  identities.emplace_back(i, i);
403 
404  else if (p % 2 != 0 && q % 2 != 0)
405  identities.emplace_back(p / 2, q / 2);
406 
407  return identities;
408  }
409  else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
410  {
411  // the FE_Nothing has no degrees of freedom, so there are no
412  // equivalencies to be recorded
413  return std::vector<std::pair<unsigned int, unsigned int>>();
414  }
415  else
416  {
417  Assert(false, ExcNotImplemented());
418  return std::vector<std::pair<unsigned int, unsigned int>>();
419  }
420 }
421 
422 
423 template <int dim>
426  const FiniteElement<dim> &fe_other,
427  const unsigned int codim) const
428 {
429  Assert(codim <= dim, ExcImpossibleInDim(dim));
430  (void)codim;
431 
432  // vertex/line/face/cell domination
433  // --------------------------------
434  if (const FE_RaviartThomasNodal<dim> *fe_rt_nodal_other =
435  dynamic_cast<const FE_RaviartThomasNodal<dim> *>(&fe_other))
436  {
437  if (this->degree < fe_rt_nodal_other->degree)
439  else if (this->degree == fe_rt_nodal_other->degree)
441  else
443  }
444  else if (const FE_Nothing<dim> *fe_nothing =
445  dynamic_cast<const FE_Nothing<dim> *>(&fe_other))
446  {
447  if (fe_nothing->is_dominating())
449  else
450  // the FE_Nothing has no degrees of freedom and it is typically used
451  // in a context where we don't require any continuity along the
452  // interface
454  }
455 
456  Assert(false, ExcNotImplemented());
458 }
459 
460 
461 
462 template <>
463 void
465  const FiniteElement<1, 1> & /*x_source_fe*/,
466  FullMatrix<double> & /*interpolation_matrix*/,
467  const unsigned int) const
468 {
469  Assert(false, ExcImpossibleInDim(1));
470 }
471 
472 
473 template <>
474 void
476  const FiniteElement<1, 1> & /*x_source_fe*/,
477  const unsigned int /*subface*/,
478  FullMatrix<double> & /*interpolation_matrix*/,
479  const unsigned int) const
480 {
481  Assert(false, ExcImpossibleInDim(1));
482 }
483 
484 
485 
486 template <int dim>
487 void
489  const FiniteElement<dim> &x_source_fe,
490  FullMatrix<double> & interpolation_matrix,
491  const unsigned int face_no) const
492 {
493  // this is only implemented, if the
494  // source FE is also a
495  // RaviartThomasNodal element
496  AssertThrow((x_source_fe.get_name().find("FE_RaviartThomasNodal<") == 0) ||
497  (dynamic_cast<const FE_RaviartThomasNodal<dim> *>(
498  &x_source_fe) != nullptr),
500 
501  Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
502  ExcDimensionMismatch(interpolation_matrix.n(),
503  this->n_dofs_per_face(face_no)));
504  Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
505  ExcDimensionMismatch(interpolation_matrix.m(),
506  x_source_fe.n_dofs_per_face(face_no)));
507 
508  // ok, source is a RaviartThomasNodal element, so we will be able to do the
509  // work
510  const FE_RaviartThomasNodal<dim> &source_fe =
511  dynamic_cast<const FE_RaviartThomasNodal<dim> &>(x_source_fe);
512 
513  // Make sure that the element for which the DoFs should be constrained is
514  // the one with the higher polynomial degree. Actually the procedure will
515  // work also if this assertion is not satisfied. But the matrices produced
516  // in that case might lead to problems in the hp-procedures, which use this
517  // method.
518  Assert(this->n_dofs_per_face(face_no) <= source_fe.n_dofs_per_face(face_no),
520 
521  // generate a quadrature with the generalized support points. This is later
522  // based as a basis for the QProjector, which returns the support points on
523  // the face.
524  Quadrature<dim - 1> quad_face_support(
525  source_fe.generalized_face_support_points[face_no]);
526 
527  // Rule of thumb for FP accuracy, that can be expected for a given
528  // polynomial degree. This value is used to cut off values close to zero.
529  double eps = 2e-13 * this->degree * (dim - 1);
530 
531  // compute the interpolation matrix by simply taking the value at the
532  // support points.
533  const Quadrature<dim> face_projection =
535  quad_face_support,
536  0);
537 
538  for (unsigned int i = 0; i < source_fe.n_dofs_per_face(face_no); ++i)
539  {
540  const Point<dim> &p = face_projection.point(i);
541 
542  for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
543  {
544  double matrix_entry =
545  this->shape_value_component(this->face_to_cell_index(j, 0), p, 0);
546 
547  // Correct the interpolated value. I.e. if it is close to 1 or 0,
548  // make it exactly 1 or 0. Unfortunately, this is required to avoid
549  // problems with higher order elements.
550  if (std::fabs(matrix_entry - 1.0) < eps)
551  matrix_entry = 1.0;
552  if (std::fabs(matrix_entry) < eps)
553  matrix_entry = 0.0;
554 
555  interpolation_matrix(i, j) = matrix_entry;
556  }
557  }
558 
559 #ifdef DEBUG
560  // make sure that the row sum of each of the matrices is 1 at this
561  // point. this must be so since the shape functions sum up to 1
562  for (unsigned int j = 0; j < source_fe.n_dofs_per_face(face_no); ++j)
563  {
564  double sum = 0.;
565 
566  for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
567  sum += interpolation_matrix(j, i);
568 
569  Assert(std::fabs(sum - 1) < 2e-13 * this->degree * (dim - 1),
570  ExcInternalError());
571  }
572 #endif
573 }
574 
575 
576 template <int dim>
577 void
579  const FiniteElement<dim> &x_source_fe,
580  const unsigned int subface,
581  FullMatrix<double> & interpolation_matrix,
582  const unsigned int face_no) const
583 {
584  // this is only implemented, if the source FE is also a RaviartThomasNodal
585  // element
586  AssertThrow((x_source_fe.get_name().find("FE_RaviartThomasNodal<") == 0) ||
587  (dynamic_cast<const FE_RaviartThomasNodal<dim> *>(
588  &x_source_fe) != nullptr),
590 
591  Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
592  ExcDimensionMismatch(interpolation_matrix.n(),
593  this->n_dofs_per_face(face_no)));
594  Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
595  ExcDimensionMismatch(interpolation_matrix.m(),
596  x_source_fe.n_dofs_per_face(face_no)));
597 
598  // ok, source is a RaviartThomasNodal element, so we will be able to do the
599  // work
600  const FE_RaviartThomasNodal<dim> &source_fe =
601  dynamic_cast<const FE_RaviartThomasNodal<dim> &>(x_source_fe);
602 
603  // Make sure that the element for which the DoFs should be constrained is
604  // the one with the higher polynomial degree. Actually the procedure will
605  // work also if this assertion is not satisfied. But the matrices produced
606  // in that case might lead to problems in the hp-procedures, which use this
607  // method.
608  Assert(this->n_dofs_per_face(face_no) <= source_fe.n_dofs_per_face(face_no),
610 
611  // generate a quadrature with the generalized support points. This is later
612  // based as a basis for the QProjector, which returns the support points on
613  // the face.
614  Quadrature<dim - 1> quad_face_support(
615  source_fe.generalized_face_support_points[face_no]);
616 
617  // Rule of thumb for FP accuracy, that can be expected for a given
618  // polynomial degree. This value is used to cut off values close to zero.
619  double eps = 2e-13 * this->degree * (dim - 1);
620 
621  // compute the interpolation matrix by simply taking the value at the
622  // support points.
623  const Quadrature<dim> subface_projection =
625  quad_face_support,
626  0,
627  subface);
628 
629  for (unsigned int i = 0; i < source_fe.n_dofs_per_face(face_no); ++i)
630  {
631  const Point<dim> &p = subface_projection.point(i);
632 
633  for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
634  {
635  double matrix_entry =
636  this->shape_value_component(this->face_to_cell_index(j, 0), p, 0);
637 
638  // Correct the interpolated value. I.e. if it is close to 1 or 0,
639  // make it exactly 1 or 0. Unfortunately, this is required to avoid
640  // problems with higher order elements.
641  if (std::fabs(matrix_entry - 1.0) < eps)
642  matrix_entry = 1.0;
643  if (std::fabs(matrix_entry) < eps)
644  matrix_entry = 0.0;
645 
646  interpolation_matrix(i, j) = matrix_entry;
647  }
648  }
649 
650 #ifdef DEBUG
651  // make sure that the row sum of each of the matrices is 1 at this
652  // point. this must be so since the shape functions sum up to 1
653  for (unsigned int j = 0; j < source_fe.n_dofs_per_face(face_no); ++j)
654  {
655  double sum = 0.;
656 
657  for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
658  sum += interpolation_matrix(j, i);
659 
660  Assert(std::fabs(sum - 1) < 2e-13 * this->degree * (dim - 1),
661  ExcInternalError());
662  }
663 #endif
664 }
665 
666 
667 
668 template <int dim>
669 const FullMatrix<double> &
671  const unsigned int child,
672  const RefinementCase<dim> &refinement_case) const
673 {
674  AssertIndexRange(refinement_case,
676  Assert(refinement_case != RefinementCase<dim>::no_refinement,
677  ExcMessage(
678  "Prolongation matrices are only available for refined cells!"));
679  AssertIndexRange(child, GeometryInfo<dim>::n_children(refinement_case));
680 
681  // initialization upon first request
682  if (this->prolongation[refinement_case - 1][child].n() == 0)
683  {
684  std::lock_guard<std::mutex> lock(this->mutex);
685 
686  // if matrix got updated while waiting for the lock
687  if (this->prolongation[refinement_case - 1][child].n() ==
688  this->n_dofs_per_cell())
689  return this->prolongation[refinement_case - 1][child];
690 
691  // now do the work. need to get a non-const version of data in order to
692  // be able to modify them inside a const function
693  FE_RaviartThomasNodal<dim> &this_nonconst =
694  const_cast<FE_RaviartThomasNodal<dim> &>(*this);
695  if (refinement_case == RefinementCase<dim>::isotropic_refinement)
696  {
697  std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
699  isotropic_matrices.back().resize(
701  FullMatrix<double>(this->n_dofs_per_cell(),
702  this->n_dofs_per_cell()));
703  FETools::compute_embedding_matrices(*this, isotropic_matrices, true);
704  this_nonconst.prolongation[refinement_case - 1].swap(
705  isotropic_matrices.back());
706  }
707  else
708  {
709  // must compute both restriction and prolongation matrices because
710  // we only check for their size and the reinit call initializes them
711  // all
714  this_nonconst.prolongation);
716  this_nonconst.restriction);
717  }
718  }
719 
720  // finally return the matrix
721  return this->prolongation[refinement_case - 1][child];
722 }
723 
724 
725 
726 template <int dim>
727 const FullMatrix<double> &
729  const unsigned int child,
730  const RefinementCase<dim> &refinement_case) const
731 {
732  AssertIndexRange(refinement_case,
734  Assert(refinement_case != RefinementCase<dim>::no_refinement,
735  ExcMessage(
736  "Restriction matrices are only available for refined cells!"));
737  AssertIndexRange(child, GeometryInfo<dim>::n_children(refinement_case));
738 
739  // initialization upon first request
740  if (this->restriction[refinement_case - 1][child].n() == 0)
741  {
742  std::lock_guard<std::mutex> lock(this->mutex);
743 
744  // if matrix got updated while waiting for the lock...
745  if (this->restriction[refinement_case - 1][child].n() ==
746  this->n_dofs_per_cell())
747  return this->restriction[refinement_case - 1][child];
748 
749  // now do the work. need to get a non-const version of data in order to
750  // be able to modify them inside a const function
751  FE_RaviartThomasNodal<dim> &this_nonconst =
752  const_cast<FE_RaviartThomasNodal<dim> &>(*this);
753  if (refinement_case == RefinementCase<dim>::isotropic_refinement)
754  {
755  std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
757  isotropic_matrices.back().resize(
759  FullMatrix<double>(this->n_dofs_per_cell(),
760  this->n_dofs_per_cell()));
761  FETools::compute_projection_matrices(*this, isotropic_matrices, true);
762  this_nonconst.restriction[refinement_case - 1].swap(
763  isotropic_matrices.back());
764  }
765  else
766  {
767  // must compute both restriction and prolongation matrices because
768  // we only check for their size and the reinit call initializes them
769  // all
772  this_nonconst.prolongation);
774  this_nonconst.restriction);
775  }
776  }
777 
778  // finally return the matrix
779  return this->restriction[refinement_case - 1][child];
780 }
781 
782 
783 
784 // explicit instantiations
785 #include "fe_raviart_thomas_nodal.inst"
786 
787 
std::vector< MappingKind > mapping_kind
virtual std::unique_ptr< FiniteElement< dim, dim > > clone() const override
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const override
virtual const FullMatrix< double > & get_restriction_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const override
void initialize_quad_dof_index_permutation_and_sign_change()
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_line_dof_identities(const FiniteElement< dim > &fe_other) const override
virtual void get_face_interpolation_matrix(const FiniteElement< dim > &source, FullMatrix< double > &matrix, const unsigned int face_no=0) const override
virtual std::string get_name() 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 hp_constraints_are_implemented() const override
virtual void get_subface_interpolation_matrix(const FiniteElement< dim > &source, const unsigned int subface, FullMatrix< double > &matrix, const unsigned int face_no=0) const override
virtual FiniteElementDomination::Domination compare_for_domination(const FiniteElement< dim > &fe_other, const unsigned int codim=0) const override final
virtual const FullMatrix< double > & get_prolongation_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const override
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_vertex_dof_identities(const FiniteElement< dim > &fe_other) const override
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_quad_dof_identities(const FiniteElement< dim > &fe_other, const unsigned int face_no=0) const override
FE_RaviartThomasNodal(const unsigned int p)
const unsigned int degree
Definition: fe_data.h:449
unsigned int n_dofs_per_face(unsigned int face_no=0, unsigned int child=0) const
virtual std::string get_name() const =0
std::vector< std::vector< FullMatrix< double > > > restriction
Definition: fe.h:2386
void reinit_restriction_and_prolongation_matrices(const bool isotropic_restriction_only=false, const bool isotropic_prolongation_only=false)
std::vector< std::vector< Point< dim - 1 > > > generalized_face_support_points
Definition: fe.h:2443
FullMatrix< double > interface_constraints
Definition: fe.h:2412
std::vector< Point< dim > > generalized_support_points
Definition: fe.h:2437
std::vector< std::vector< FullMatrix< double > > > prolongation
Definition: fe.h:2400
size_type n() const
size_type m() const
Definition: point.h:111
std::vector< Point< dim > > get_polynomial_support_points() const
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 void project_to_face(const ReferenceCell reference_cell, const SubQuadrature &quadrature, const unsigned int face_no, std::vector< Point< dim >> &q_points)
SymmetricTensor< rank, dim, Number > sum(const SymmetricTensor< rank, dim, Number > &local, const MPI_Comm &mpi_communicator)
Definition: vector.h:110
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
@ mapping_raviart_thomas
Definition: mapping.h:127
Expression fabs(const Expression &x)
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)
void compute_projection_matrices(const FiniteElement< dim, spacedim > &fe, std::vector< std::vector< FullMatrix< number >>> &matrices, const bool isotropic_only=false)
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
std::string to_string(const T &t)
Definition: patterns.h:2393
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618