deal.II version GIT relicensing-3718-gc13d52846a 2025-07-09 19:30:01+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\}}\)
Loading...
Searching...
No Matches
fe_raviart_thomas_nodal.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2005 - 2025 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15
20
22
23#include <deal.II/fe/fe.h>
26#include <deal.II/fe/fe_tools.h>
27#include <deal.II/fe/mapping.h>
28
29#include <memory>
30#include <sstream>
31
32
34
35
36// ---------------- polynomial class for FE_RaviartThomasNodal ---------------
37
38namespace
39{
40 // Return a vector of "dofs per object" where the components of the returned
41 // vector refer to:
42 // 0 = vertex
43 // 1 = edge
44 // 2 = face (which is a cell in 2d)
45 // 3 = cell
46 std::vector<unsigned int>
47 get_rt_dpo_vector(const unsigned int dim, const unsigned int degree)
48 {
49 std::vector<unsigned int> dpo(dim + 1);
50 dpo[0] = 0;
51 dpo[1] = 0;
52 unsigned int dofs_per_face = 1;
53 for (unsigned int d = 1; d < dim; ++d)
54 dofs_per_face *= (degree + 1);
55
56 dpo[dim - 1] = dofs_per_face;
57 dpo[dim] = dim * degree * dofs_per_face;
58
59 return dpo;
60 }
61} // namespace
62
63
64
65// --------------------- actual implementation of element --------------------
66
67template <int dim>
69 : FE_PolyTensor<dim>(
71 degree + 1,
72 degree,
73 FE_RaviartThomas<dim>::get_lexicographic_numbering(degree)),
74 FiniteElementData<dim>(get_rt_dpo_vector(dim, degree),
75 dim,
76 degree + 1,
77 FiniteElementData<dim>::Hdiv),
78 std::vector<bool>(1, false),
79 std::vector<ComponentMask>(
80 PolynomialsVectorAnisotropic<dim>::n_polynomials(degree + 1, degree),
81 ComponentMask(std::vector<bool>(dim, true))))
82{
83 Assert(dim >= 2, ExcImpossibleInDim(dim));
84
86
87 const std::vector<unsigned int> numbering =
89
90 // First, initialize the generalized support points and quadrature weights,
91 // since they are required for interpolation.
96 this->n_dofs_per_cell());
97
98 const unsigned int face_no = 0;
99 if (dim > 1)
100 this->generalized_face_support_points[face_no] =
101 degree == 0 ? QGauss<dim - 1>(1).get_points() :
102 QGaussLobatto<dim - 1>(degree + 1).get_points();
103
105 for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_face; ++i)
106 face_embeddings[i].reinit(this->n_dofs_per_face(face_no),
107 this->n_dofs_per_face(face_no));
108 FETools::compute_face_embedding_matrices<dim, double>(*this,
109 face_embeddings,
110 0,
111 0);
113 this->n_dofs_per_face(face_no),
114 this->n_dofs_per_face(face_no));
115 unsigned int target_row = 0;
116 for (unsigned int d = 0; d < GeometryInfo<dim>::max_children_per_face; ++d)
117 for (unsigned int i = 0; i < face_embeddings[d].m(); ++i)
118 {
119 for (unsigned int j = 0; j < face_embeddings[d].n(); ++j)
120 this->interface_constraints(target_row, j) = face_embeddings[d](i, j);
121 ++target_row;
122 }
123
124 // We need to initialize the dof permutation table and the one for the sign
125 // change.
127}
128
129
130
131template <int dim>
132std::string
134{
135 // note that the FETools::get_fe_by_name function depends on the particular
136 // format of the string this function returns, so they have to be kept in
137 // synch
138
139 // note that this->degree is the maximal polynomial degree and is thus one
140 // higher than the argument given to the constructor
141 return "FE_RaviartThomasNodal<" + std::to_string(dim) + ">(" +
142 std::to_string(this->degree - 1) + ")";
143}
144
145
146template <int dim>
147std::unique_ptr<FiniteElement<dim, dim>>
149{
150 return std::make_unique<FE_RaviartThomasNodal<dim>>(*this);
151}
152
153
154//---------------------------------------------------------------------------
155// Auxiliary and internal functions
156//---------------------------------------------------------------------------
157
158
159
160template <int dim>
161void
163 dim>::initialize_quad_dof_index_permutation_and_sign_change()
164{
165 // for 1d and 2d, do nothing
166 if (dim < 3)
167 return;
168
169 const unsigned int n = this->degree;
170 const unsigned int face_no = 0;
171 Assert(n * n == this->n_dofs_per_quad(face_no), ExcInternalError());
172 for (unsigned int local = 0; local < this->n_dofs_per_quad(face_no); ++local)
173 // face support points are in lexicographic ordering with x running
174 // fastest. invert that (y running fastest)
175 {
176 unsigned int i = local % n, j = local / n;
177
178 // face_orientation=false, face_rotation=false, face_flip=false
179 this->adjust_quad_dof_index_for_face_orientation_table[face_no](
180 local, internal::combined_face_orientation(false, false, false)) =
181 j + i * n - local;
182 // face_orientation=false, face_rotation=true, face_flip=false
183 this->adjust_quad_dof_index_for_face_orientation_table[face_no](
184 local, internal::combined_face_orientation(false, true, false)) =
185 i + (n - 1 - j) * n - local;
186 // face_orientation=false, face_rotation=false, face_flip=true
187 this->adjust_quad_dof_index_for_face_orientation_table[face_no](
188 local, internal::combined_face_orientation(false, false, true)) =
189 (n - 1 - j) + (n - 1 - i) * n - local;
190 // face_orientation=false, face_rotation=true, face_flip=true
191 this->adjust_quad_dof_index_for_face_orientation_table[face_no](
192 local, internal::combined_face_orientation(false, true, true)) =
193 (n - 1 - i) + j * n - local;
194 // face_orientation=true, face_rotation=false, face_flip=false
195 this->adjust_quad_dof_index_for_face_orientation_table[face_no](
196 local, internal::combined_face_orientation(true, false, false)) = 0;
197 // face_orientation=true, face_rotation=true, face_flip=false
198 this->adjust_quad_dof_index_for_face_orientation_table[face_no](
199 local, internal::combined_face_orientation(true, true, false)) =
200 j + (n - 1 - i) * n - local;
201 // face_orientation=true, face_rotation=false, face_flip=true
202 this->adjust_quad_dof_index_for_face_orientation_table[face_no](
203 local, internal::combined_face_orientation(true, false, true)) =
204 (n - 1 - i) + (n - 1 - j) * n - local;
205 // face_orientation=true, face_rotation=true, face_flip=true
206 this->adjust_quad_dof_index_for_face_orientation_table[face_no](
207 local, internal::combined_face_orientation(true, true, true)) =
208 (n - 1 - j) + i * n - local;
209
210 // for face_orientation == false, we need to switch the sign
211 for (const bool rotation : {false, true})
212 for (const bool flip : {false, true})
213 this->adjust_quad_dof_sign_for_face_orientation_table[face_no](
214 local, internal::combined_face_orientation(false, rotation, flip)) =
215 1;
216 }
217}
218
219
220
221template <int dim>
222bool
224 const unsigned int shape_index,
225 const unsigned int face_index) const
226{
227 AssertIndexRange(shape_index, this->n_dofs_per_cell());
229
230 // The first degrees of freedom are on the faces and each face has degree
231 // degrees.
232 const unsigned int support_face = shape_index / this->n_dofs_per_face();
233
234 // The only thing we know for sure is that shape functions with support on
235 // one face are zero on the opposite face.
236 if (support_face < GeometryInfo<dim>::faces_per_cell)
237 return (face_index != GeometryInfo<dim>::opposite_face[support_face]);
238
239 // In all other cases, return true, which is safe
240 return true;
241}
242
243
244
245template <int dim>
246void
249 const std::vector<Vector<double>> &support_point_values,
250 std::vector<double> &nodal_values) const
251{
252 Assert(support_point_values.size() == this->generalized_support_points.size(),
253 ExcDimensionMismatch(support_point_values.size(),
254 this->generalized_support_points.size()));
255 Assert(nodal_values.size() == this->n_dofs_per_cell(),
256 ExcDimensionMismatch(nodal_values.size(), this->n_dofs_per_cell()));
257 Assert(support_point_values[0].size() == this->n_components(),
258 ExcDimensionMismatch(support_point_values[0].size(),
259 this->n_components()));
260
261 // First do interpolation on faces. There, the component evaluated depends
262 // on the face direction and orientation.
263 unsigned int fbase = 0;
264 unsigned int f = 0;
265 for (; f < GeometryInfo<dim>::faces_per_cell;
266 ++f, fbase += this->n_dofs_per_face(f))
267 {
268 for (unsigned int i = 0; i < this->n_dofs_per_face(f); ++i)
269 {
270 nodal_values[fbase + i] = support_point_values[fbase + i](
272 }
273 }
274
275 // The remaining points form dim chunks, one for each component
276 const unsigned int istep = (this->n_dofs_per_cell() - fbase) / dim;
277 Assert((this->n_dofs_per_cell() - fbase) % dim == 0, ExcInternalError());
278
279 f = 0;
280 while (fbase < this->n_dofs_per_cell())
281 {
282 for (unsigned int i = 0; i < istep; ++i)
283 {
284 nodal_values[fbase + i] = support_point_values[fbase + i](f);
285 }
286 fbase += istep;
287 ++f;
288 }
289 Assert(fbase == this->n_dofs_per_cell(), ExcInternalError());
290}
291
292
293
294// TODO: There are tests that check that the following few functions don't
295// produce assertion failures, but none that actually check whether they do the
296// right thing. one example for such a test would be to project a function onto
297// an hp-space and make sure that the convergence order is correct with regard
298// to the lowest used polynomial degree
299
300template <int dim>
301bool
306
307
308template <int dim>
309std::vector<std::pair<unsigned int, unsigned int>>
311 const FiniteElement<dim> &fe_other) const
312{
313 // we can presently only compute these identities if both FEs are
314 // FE_RaviartThomasNodals or the other is FE_Nothing. In either case, no
315 // dofs are assigned on the vertex, so we shouldn't be getting here at all.
316 if (dynamic_cast<const FE_RaviartThomasNodal<dim> *>(&fe_other) != nullptr)
317 return std::vector<std::pair<unsigned int, unsigned int>>();
318 else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
319 return std::vector<std::pair<unsigned int, unsigned int>>();
320 else
321 {
323 return std::vector<std::pair<unsigned int, unsigned int>>();
324 }
325}
326
327
328
329template <int dim>
330std::vector<std::pair<unsigned int, unsigned int>>
332 const FiniteElement<dim> &fe_other) const
333{
334 // we can presently only compute these identities if both FEs are
335 // FE_RaviartThomasNodals or if the other one is FE_Nothing
336 if (const FE_RaviartThomasNodal<dim> *fe_q_other =
337 dynamic_cast<const FE_RaviartThomasNodal<dim> *>(&fe_other))
338 {
339 // dofs are located on faces; these are only lines in 2d
340 if (dim != 2)
341 return std::vector<std::pair<unsigned int, unsigned int>>();
342
343 // dofs are located along lines, so two dofs are identical only if in
344 // the following two cases (remember that the face support points are
345 // Gauss points):
346 // 1. this->degree = fe_q_other->degree,
347 // in the case, all the dofs on the line are identical
348 // 2. this->degree-1 and fe_q_other->degree-1
349 // are both even, i.e. this->dof_per_line and fe_q_other->dof_per_line
350 // are both odd, there exists only one point (the middle one) such
351 // that dofs are identical on this point
352 //
353 // to understand this, note that this->degree is the *maximal*
354 // polynomial degree, and is thus one higher than the argument given to
355 // the constructor
356 const unsigned int p = this->degree - 1;
357 const unsigned int q = fe_q_other->degree - 1;
358
359 std::vector<std::pair<unsigned int, unsigned int>> identities;
360
361 if (p == q)
362 for (unsigned int i = 0; i < p + 1; ++i)
363 identities.emplace_back(i, i);
364
365 else if (p % 2 == 0 && q % 2 == 0)
366 identities.emplace_back(p / 2, q / 2);
367
368 return identities;
369 }
370 else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
371 {
372 // the FE_Nothing has no degrees of freedom, so there are no
373 // equivalencies to be recorded
374 return std::vector<std::pair<unsigned int, unsigned int>>();
375 }
376 else
377 {
379 return std::vector<std::pair<unsigned int, unsigned int>>();
380 }
381}
382
383
384template <int dim>
385std::vector<std::pair<unsigned int, unsigned int>>
387 const FiniteElement<dim> &fe_other,
388 const unsigned int face_no) const
389{
390 // we can presently only compute these identities if both FEs are
391 // FE_RaviartThomasNodals or if the other one is FE_Nothing
392 if (const FE_RaviartThomasNodal<dim> *fe_q_other =
393 dynamic_cast<const FE_RaviartThomasNodal<dim> *>(&fe_other))
394 {
395 // dofs are located on faces; these are only quads in 3d
396 if (dim != 3)
397 return std::vector<std::pair<unsigned int, unsigned int>>();
398
399 // this works exactly like the line case above
400 const unsigned int p = this->n_dofs_per_quad(face_no);
401
402 AssertDimension(fe_q_other->n_unique_faces(), 1);
403 const unsigned int q = fe_q_other->n_dofs_per_quad(0);
404
405 std::vector<std::pair<unsigned int, unsigned int>> identities;
406
407 if (p == q)
408 for (unsigned int i = 0; i < p; ++i)
409 identities.emplace_back(i, i);
410
411 else if (p % 2 != 0 && q % 2 != 0)
412 identities.emplace_back(p / 2, q / 2);
413
414 return identities;
415 }
416 else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
417 {
418 // the FE_Nothing has no degrees of freedom, so there are no
419 // equivalencies to be recorded
420 return std::vector<std::pair<unsigned int, unsigned int>>();
421 }
422 else
423 {
425 return std::vector<std::pair<unsigned int, unsigned int>>();
426 }
427}
428
429
430template <int dim>
433 const FiniteElement<dim> &fe_other,
434 const unsigned int codim) const
435{
436 Assert(codim <= dim, ExcImpossibleInDim(dim));
437 (void)codim;
438
439 // vertex/line/face/cell domination
440 // --------------------------------
441 if (const FE_RaviartThomasNodal<dim> *fe_rt_nodal_other =
442 dynamic_cast<const FE_RaviartThomasNodal<dim> *>(&fe_other))
443 {
444 if (this->degree < fe_rt_nodal_other->degree)
446 else if (this->degree == fe_rt_nodal_other->degree)
448 else
450 }
451 else if (const FE_Nothing<dim> *fe_nothing =
452 dynamic_cast<const FE_Nothing<dim> *>(&fe_other))
453 {
454 if (fe_nothing->is_dominating())
456 else
457 // the FE_Nothing has no degrees of freedom and it is typically used
458 // in a context where we don't require any continuity along the
459 // interface
461 }
462
465}
466
467
468
469template <>
470void
472 const FiniteElement<1, 1> & /*x_source_fe*/,
473 FullMatrix<double> & /*interpolation_matrix*/,
474 const unsigned int) const
475{
476 Assert(false, ExcImpossibleInDim(1));
477}
478
479
480template <>
481void
483 const FiniteElement<1, 1> & /*x_source_fe*/,
484 const unsigned int /*subface*/,
485 FullMatrix<double> & /*interpolation_matrix*/,
486 const unsigned int) const
487{
488 Assert(false, ExcImpossibleInDim(1));
489}
490
491
492
493template <int dim>
494void
496 const FiniteElement<dim> &x_source_fe,
497 FullMatrix<double> &interpolation_matrix,
498 const unsigned int face_no) const
499{
500 // this is only implemented, if the
501 // source FE is also a
502 // RaviartThomasNodal element
503 AssertThrow((x_source_fe.get_name().find("FE_RaviartThomasNodal<") == 0) ||
504 (dynamic_cast<const FE_RaviartThomasNodal<dim> *>(
505 &x_source_fe) != nullptr),
507
508 Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
509 ExcDimensionMismatch(interpolation_matrix.n(),
510 this->n_dofs_per_face(face_no)));
511 Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
512 ExcDimensionMismatch(interpolation_matrix.m(),
513 x_source_fe.n_dofs_per_face(face_no)));
514
515 // ok, source is a RaviartThomasNodal element, so we will be able to do the
516 // work
517 const FE_RaviartThomasNodal<dim> &source_fe =
518 dynamic_cast<const FE_RaviartThomasNodal<dim> &>(x_source_fe);
519
520 // Make sure that the element for which the DoFs should be constrained is
521 // the one with the higher polynomial degree. Actually the procedure will
522 // work also if this assertion is not satisfied. But the matrices produced
523 // in that case might lead to problems in the hp-procedures, which use this
524 // method.
525 Assert(this->n_dofs_per_face(face_no) <= source_fe.n_dofs_per_face(face_no),
527
528 // generate a quadrature with the generalized support points. This is later
529 // based as a basis for the QProjector, which returns the support points on
530 // the face.
531 Quadrature<dim - 1> quad_face_support(
532 source_fe.generalized_face_support_points[face_no]);
533
534 // Rule of thumb for FP accuracy, that can be expected for a given
535 // polynomial degree. This value is used to cut off values close to zero.
536 double eps = 2e-13 * this->degree * (dim - 1);
537
538 // compute the interpolation matrix by simply taking the value at the
539 // support points.
540 const Quadrature<dim> face_projection =
541 QProjector<dim>::project_to_face(this->reference_cell(),
542 quad_face_support,
543 0,
545
546 for (unsigned int i = 0; i < source_fe.n_dofs_per_face(face_no); ++i)
547 {
548 const Point<dim> &p = face_projection.point(i);
549
550 for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
551 {
552 double matrix_entry =
553 this->shape_value_component(this->face_to_cell_index(j, 0), p, 0);
554
555 // Correct the interpolated value. I.e. if it is close to 1 or 0,
556 // make it exactly 1 or 0. Unfortunately, this is required to avoid
557 // problems with higher order elements.
558 if (std::fabs(matrix_entry - 1.0) < eps)
559 matrix_entry = 1.0;
560 if (std::fabs(matrix_entry) < eps)
561 matrix_entry = 0.0;
562
563 interpolation_matrix(i, j) = matrix_entry;
564 }
565 }
566
567 if constexpr (running_in_debug_mode())
568 {
569 // make sure that the row sum of each of the matrices is 1 at this
570 // point. this must be so since the shape functions sum up to 1
571 for (unsigned int j = 0; j < source_fe.n_dofs_per_face(face_no); ++j)
572 {
573 double sum = 0.;
574
575 for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
576 sum += interpolation_matrix(j, i);
577
578 Assert(std::fabs(sum - 1) < 2e-13 * this->degree * (dim - 1),
580 }
581 }
582}
583
584
585template <int dim>
586void
588 const FiniteElement<dim> &x_source_fe,
589 const unsigned int subface,
590 FullMatrix<double> &interpolation_matrix,
591 const unsigned int face_no) const
592{
593 // this is only implemented, if the source FE is also a RaviartThomasNodal
594 // element
595 AssertThrow((x_source_fe.get_name().find("FE_RaviartThomasNodal<") == 0) ||
596 (dynamic_cast<const FE_RaviartThomasNodal<dim> *>(
597 &x_source_fe) != nullptr),
599
600 Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
601 ExcDimensionMismatch(interpolation_matrix.n(),
602 this->n_dofs_per_face(face_no)));
603 Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
604 ExcDimensionMismatch(interpolation_matrix.m(),
605 x_source_fe.n_dofs_per_face(face_no)));
606
607 // ok, source is a RaviartThomasNodal element, so we will be able to do the
608 // work
609 const FE_RaviartThomasNodal<dim> &source_fe =
610 dynamic_cast<const FE_RaviartThomasNodal<dim> &>(x_source_fe);
611
612 // Make sure that the element for which the DoFs should be constrained is
613 // the one with the higher polynomial degree. Actually the procedure will
614 // work also if this assertion is not satisfied. But the matrices produced
615 // in that case might lead to problems in the hp-procedures, which use this
616 // method.
617 Assert(this->n_dofs_per_face(face_no) <= source_fe.n_dofs_per_face(face_no),
619
620 // generate a quadrature with the generalized support points. This is later
621 // based as a basis for the QProjector, which returns the support points on
622 // the face.
623 Quadrature<dim - 1> quad_face_support(
624 source_fe.generalized_face_support_points[face_no]);
625
626 // Rule of thumb for FP accuracy, that can be expected for a given
627 // polynomial degree. This value is used to cut off values close to zero.
628 double eps = 2e-13 * this->degree * (dim - 1);
629
630 // compute the interpolation matrix by simply taking the value at the
631 // support points.
632 const Quadrature<dim> subface_projection =
634 this->reference_cell(),
635 quad_face_support,
636 0,
637 subface,
640
641 for (unsigned int i = 0; i < source_fe.n_dofs_per_face(face_no); ++i)
642 {
643 const Point<dim> &p = subface_projection.point(i);
644
645 for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
646 {
647 double matrix_entry =
648 this->shape_value_component(this->face_to_cell_index(j, 0), p, 0);
649
650 // Correct the interpolated value. I.e. if it is close to 1 or 0,
651 // make it exactly 1 or 0. Unfortunately, this is required to avoid
652 // problems with higher order elements.
653 if (std::fabs(matrix_entry - 1.0) < eps)
654 matrix_entry = 1.0;
655 if (std::fabs(matrix_entry) < eps)
656 matrix_entry = 0.0;
657
658 interpolation_matrix(i, j) = matrix_entry;
659 }
660 }
661
662 if constexpr (running_in_debug_mode())
663 {
664 // make sure that the row sum of each of the matrices is 1 at this
665 // point. this must be so since the shape functions sum up to 1
666 for (unsigned int j = 0; j < source_fe.n_dofs_per_face(face_no); ++j)
667 {
668 double sum = 0.;
669
670 for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
671 sum += interpolation_matrix(j, i);
672
673 Assert(std::fabs(sum - 1) < 2e-13 * this->degree * (dim - 1),
675 }
676 }
677}
678
679
680
681template <int dim>
682const FullMatrix<double> &
684 const unsigned int child,
685 const RefinementCase<dim> &refinement_case) const
686{
687 AssertIndexRange(refinement_case,
691 "Prolongation matrices are only available for refined cells!"));
693 child, this->reference_cell().template n_children<dim>(refinement_case));
694
695 // initialization upon first request
696 if (this->prolongation[refinement_case - 1][child].n() == 0)
697 {
698 std::lock_guard<std::mutex> lock(prolongation_matrix_mutex);
699
700 // if matrix got updated while waiting for the lock
701 if (this->prolongation[refinement_case - 1][child].n() ==
702 this->n_dofs_per_cell())
703 return this->prolongation[refinement_case - 1][child];
704
705 // now do the work. need to get a non-const version of data in order to
706 // be able to modify them inside a const function
707 FE_RaviartThomasNodal<dim> &this_nonconst =
708 const_cast<FE_RaviartThomasNodal<dim> &>(*this);
709 if (refinement_case == RefinementCase<dim>::isotropic_refinement)
710 {
711 std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
713 isotropic_matrices.back().resize(
714 this->reference_cell().template n_children<dim>(
715 RefinementCase<dim>(refinement_case)),
716 FullMatrix<double>(this->n_dofs_per_cell(),
717 this->n_dofs_per_cell()));
718 FETools::compute_embedding_matrices(*this, isotropic_matrices, true);
719 this_nonconst.prolongation[refinement_case - 1] =
720 std::move(isotropic_matrices.back());
721 }
722 else
723 {
724 // must compute both restriction and prolongation matrices because
725 // we only check for their size and the reinit call initializes them
726 // all
729 this_nonconst.prolongation);
731 this_nonconst.restriction);
732 }
733 }
734
735 // finally return the matrix
736 return this->prolongation[refinement_case - 1][child];
737}
738
739
740
741template <int dim>
742const FullMatrix<double> &
744 const unsigned int child,
745 const RefinementCase<dim> &refinement_case) const
746{
747 AssertIndexRange(refinement_case,
751 "Restriction matrices are only available for refined cells!"));
753 child, this->reference_cell().template n_children<dim>(refinement_case));
754
755 // initialization upon first request
756 if (this->restriction[refinement_case - 1][child].n() == 0)
757 {
758 std::lock_guard<std::mutex> lock(restriction_matrix_mutex);
759
760 // if matrix got updated while waiting for the lock...
761 if (this->restriction[refinement_case - 1][child].n() ==
762 this->n_dofs_per_cell())
763 return this->restriction[refinement_case - 1][child];
764
765 // now do the work. need to get a non-const version of data in order to
766 // be able to modify them inside a const function
767 FE_RaviartThomasNodal<dim> &this_nonconst =
768 const_cast<FE_RaviartThomasNodal<dim> &>(*this);
769 if (refinement_case == RefinementCase<dim>::isotropic_refinement)
770 {
771 std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
773 isotropic_matrices.back().resize(
774 this->reference_cell().template n_children<dim>(
775 RefinementCase<dim>(refinement_case)),
776 FullMatrix<double>(this->n_dofs_per_cell(),
777 this->n_dofs_per_cell()));
778 FETools::compute_projection_matrices(*this, isotropic_matrices, true);
779 this_nonconst.restriction[refinement_case - 1] =
780 std::move(isotropic_matrices.back());
781 }
782 else
783 {
784 // must compute both restriction and prolongation matrices because
785 // we only check for their size and the reinit call initializes them
786 // all
789 this_nonconst.prolongation);
791 this_nonconst.restriction);
792 }
793 }
794
795 // finally return the matrix
796 return this->restriction[refinement_case - 1][child];
797}
798
799
800
801// explicit instantiations
802#include "fe/fe_raviart_thomas_nodal.inst"
803
804
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)
static std::vector< unsigned int > get_lexicographic_numbering(const unsigned int degree)
const unsigned int degree
Definition fe_data.h:452
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:2524
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:2581
FullMatrix< double > interface_constraints
Definition fe.h:2550
std::vector< Point< dim > > generalized_support_points
Definition fe.h:2575
std::vector< std::vector< FullMatrix< double > > > prolongation
Definition fe.h:2538
size_type n() const
size_type m() const
Definition point.h:113
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:35
constexpr bool running_in_debug_mode()
Definition config.h:73
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
#define DEAL_II_NOT_IMPLEMENTED()
#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:136
std::size_t size
Definition mpi.cc:749
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)
types::geometric_orientation combined_face_orientation(const bool face_orientation, const bool face_rotation, const bool face_flip)
constexpr types::geometric_orientation default_geometric_orientation
Definition types.h:352
STL namespace.