Reference documentation for deal.II version 9.5.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\}}\)
Loading...
Searching...
No Matches
fe_raviart_thomas_nodal.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2003 - 2023 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>
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
39namespace
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
68template <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
126template <int dim>
127std::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
141template <int dim>
142std::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
155template <int dim>
156void
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
214template <int dim>
215bool
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
238template <int dim>
239void
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
293template <int dim>
294bool
296{
297 return true;
298}
299
300
301template <int dim>
302std::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
322template <int dim>
323std::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
377template <int dim>
378std::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
423template <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
462template <>
463void
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
473template <>
474void
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
486template <int dim>
487void
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 =
534 QProjector<dim>::project_to_face(this->reference_cell(),
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),
571 }
572#endif
573}
574
575
576template <int dim>
577void
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 =
624 QProjector<dim>::project_to_subface(this->reference_cell(),
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),
662 }
663#endif
664}
665
666
667
668template <int dim>
669const FullMatrix<double> &
671 const unsigned int child,
672 const RefinementCase<dim> &refinement_case) const
673{
674 AssertIndexRange(refinement_case,
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
726template <int dim>
727const FullMatrix<double> &
729 const unsigned int child,
730 const RefinementCase<dim> &refinement_case) const
731{
732 AssertIndexRange(refinement_case,
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 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 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 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:453
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:2402
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:2459
FullMatrix< double > interface_constraints
Definition fe.h:2428
std::vector< Point< dim > > generalized_support_points
Definition fe.h:2453
std::vector< std::vector< FullMatrix< double > > > prolongation
Definition fe.h:2416
size_type n() const
size_type m() const
Definition point.h:112
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)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
@ mapping_raviart_thomas
Definition mapping.h:133
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)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
STL namespace.