Reference documentation for deal.II version 9.6.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// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2005 - 2024 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>(
70 PolynomialsRaviartThomas<dim>(degree + 1, degree),
71 FiniteElementData<dim>(get_rt_dpo_vector(dim, degree),
72 dim,
73 degree + 1,
74 FiniteElementData<dim>::Hdiv),
75 std::vector<bool>(1, false),
76 std::vector<ComponentMask>(
77 PolynomialsRaviartThomas<dim>::n_polynomials(degree + 1, degree),
78 ComponentMask(std::vector<bool>(dim, true))))
79{
80 Assert(dim >= 2, ExcImpossibleInDim(dim));
81
83
84 // First, initialize the generalized support points and quadrature weights,
85 // since they are required for interpolation.
90 this->n_dofs_per_cell());
91
92 const unsigned int face_no = 0;
93 if (dim > 1)
94 this->generalized_face_support_points[face_no] =
95 degree == 0 ? QGauss<dim - 1>(1).get_points() :
96 QGaussLobatto<dim - 1>(degree + 1).get_points();
97
99 for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_face; ++i)
100 face_embeddings[i].reinit(this->n_dofs_per_face(face_no),
101 this->n_dofs_per_face(face_no));
103 face_embeddings,
104 0,
105 0);
107 this->n_dofs_per_face(face_no),
108 this->n_dofs_per_face(face_no));
109 unsigned int target_row = 0;
110 for (unsigned int d = 0; d < GeometryInfo<dim>::max_children_per_face; ++d)
111 for (unsigned int i = 0; i < face_embeddings[d].m(); ++i)
112 {
113 for (unsigned int j = 0; j < face_embeddings[d].n(); ++j)
114 this->interface_constraints(target_row, j) = face_embeddings[d](i, j);
115 ++target_row;
116 }
117
118 // We need to initialize the dof permutation table and the one for the sign
119 // change.
121}
122
123
124
125template <int dim>
126std::string
128{
129 // note that the FETools::get_fe_by_name function depends on the particular
130 // format of the string this function returns, so they have to be kept in
131 // synch
132
133 // note that this->degree is the maximal polynomial degree and is thus one
134 // higher than the argument given to the constructor
135 return "FE_RaviartThomasNodal<" + std::to_string(dim) + ">(" +
136 std::to_string(this->degree - 1) + ")";
137}
138
139
140template <int dim>
141std::unique_ptr<FiniteElement<dim, dim>>
143{
144 return std::make_unique<FE_RaviartThomasNodal<dim>>(*this);
145}
146
147
148//---------------------------------------------------------------------------
149// Auxiliary and internal functions
150//---------------------------------------------------------------------------
151
152
153
154template <int dim>
155void
157 dim>::initialize_quad_dof_index_permutation_and_sign_change()
158{
159 // for 1d and 2d, do nothing
160 if (dim < 3)
161 return;
162
163 const unsigned int n = this->degree;
164 const unsigned int face_no = 0;
165 Assert(n * n == this->n_dofs_per_quad(face_no), ExcInternalError());
166 for (unsigned int local = 0; local < this->n_dofs_per_quad(face_no); ++local)
167 // face support points are in lexicographic ordering with x running
168 // fastest. invert that (y running fastest)
169 {
170 unsigned int i = local % n, j = local / n;
171
172 // face_orientation=false, face_rotation=false, face_flip=false
173 this->adjust_quad_dof_index_for_face_orientation_table[face_no](
174 local, internal::combined_face_orientation(false, false, false)) =
175 j + i * n - local;
176 // face_orientation=false, face_rotation=true, face_flip=false
177 this->adjust_quad_dof_index_for_face_orientation_table[face_no](
178 local, internal::combined_face_orientation(false, true, false)) =
179 i + (n - 1 - j) * n - local;
180 // face_orientation=false, face_rotation=false, face_flip=true
181 this->adjust_quad_dof_index_for_face_orientation_table[face_no](
182 local, internal::combined_face_orientation(false, false, true)) =
183 (n - 1 - j) + (n - 1 - i) * n - local;
184 // face_orientation=false, face_rotation=true, face_flip=true
185 this->adjust_quad_dof_index_for_face_orientation_table[face_no](
186 local, internal::combined_face_orientation(false, true, true)) =
187 (n - 1 - i) + j * n - local;
188 // face_orientation=true, face_rotation=false, face_flip=false
189 this->adjust_quad_dof_index_for_face_orientation_table[face_no](
190 local, internal::combined_face_orientation(true, false, false)) = 0;
191 // face_orientation=true, face_rotation=true, face_flip=false
192 this->adjust_quad_dof_index_for_face_orientation_table[face_no](
193 local, internal::combined_face_orientation(true, true, false)) =
194 j + (n - 1 - i) * n - local;
195 // face_orientation=true, face_rotation=false, face_flip=true
196 this->adjust_quad_dof_index_for_face_orientation_table[face_no](
197 local, internal::combined_face_orientation(true, false, true)) =
198 (n - 1 - i) + (n - 1 - j) * n - local;
199 // face_orientation=true, face_rotation=true, face_flip=true
200 this->adjust_quad_dof_index_for_face_orientation_table[face_no](
201 local, internal::combined_face_orientation(true, true, true)) =
202 (n - 1 - j) + i * n - local;
203
204 // for face_orientation == false, we need to switch the sign
205 for (const bool rotation : {false, true})
206 for (const bool flip : {false, true})
207 this->adjust_quad_dof_sign_for_face_orientation_table[face_no](
208 local, internal::combined_face_orientation(false, rotation, flip)) =
209 1;
210 }
211}
212
213
214
215template <int dim>
216bool
218 const unsigned int shape_index,
219 const unsigned int face_index) const
220{
221 AssertIndexRange(shape_index, this->n_dofs_per_cell());
223
224 // The first degrees of freedom are on the faces and each face has degree
225 // degrees.
226 const unsigned int support_face = shape_index / this->n_dofs_per_face();
227
228 // The only thing we know for sure is that shape functions with support on
229 // one face are zero on the opposite face.
230 if (support_face < GeometryInfo<dim>::faces_per_cell)
231 return (face_index != GeometryInfo<dim>::opposite_face[support_face]);
232
233 // In all other cases, return true, which is safe
234 return true;
235}
236
237
238
239template <int dim>
240void
243 const std::vector<Vector<double>> &support_point_values,
244 std::vector<double> &nodal_values) const
245{
246 Assert(support_point_values.size() == this->generalized_support_points.size(),
247 ExcDimensionMismatch(support_point_values.size(),
248 this->generalized_support_points.size()));
249 Assert(nodal_values.size() == this->n_dofs_per_cell(),
250 ExcDimensionMismatch(nodal_values.size(), this->n_dofs_per_cell()));
251 Assert(support_point_values[0].size() == this->n_components(),
252 ExcDimensionMismatch(support_point_values[0].size(),
253 this->n_components()));
254
255 // First do interpolation on faces. There, the component evaluated depends
256 // on the face direction and orientation.
257 unsigned int fbase = 0;
258 unsigned int f = 0;
259 for (; f < GeometryInfo<dim>::faces_per_cell;
260 ++f, fbase += this->n_dofs_per_face(f))
261 {
262 for (unsigned int i = 0; i < this->n_dofs_per_face(f); ++i)
263 {
264 nodal_values[fbase + i] = support_point_values[fbase + i](
266 }
267 }
268
269 // The remaining points form dim chunks, one for each component
270 const unsigned int istep = (this->n_dofs_per_cell() - fbase) / dim;
271 Assert((this->n_dofs_per_cell() - fbase) % dim == 0, ExcInternalError());
272
273 f = 0;
274 while (fbase < this->n_dofs_per_cell())
275 {
276 for (unsigned int i = 0; i < istep; ++i)
277 {
278 nodal_values[fbase + i] = support_point_values[fbase + i](f);
279 }
280 fbase += istep;
281 ++f;
282 }
283 Assert(fbase == this->n_dofs_per_cell(), ExcInternalError());
284}
285
286
287
288// TODO: There are tests that check that the following few functions don't
289// produce assertion failures, but none that actually check whether they do the
290// right thing. one example for such a test would be to project a function onto
291// an hp-space and make sure that the convergence order is correct with regard
292// to the lowest used polynomial degree
293
294template <int dim>
295bool
300
301
302template <int dim>
303std::vector<std::pair<unsigned int, unsigned int>>
305 const FiniteElement<dim> &fe_other) const
306{
307 // we can presently only compute these identities if both FEs are
308 // FE_RaviartThomasNodals or the other is FE_Nothing. In either case, no
309 // dofs are assigned on the vertex, so we shouldn't be getting here at all.
310 if (dynamic_cast<const FE_RaviartThomasNodal<dim> *>(&fe_other) != nullptr)
311 return std::vector<std::pair<unsigned int, unsigned int>>();
312 else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
313 return std::vector<std::pair<unsigned int, unsigned int>>();
314 else
315 {
317 return std::vector<std::pair<unsigned int, unsigned int>>();
318 }
319}
320
321
322
323template <int dim>
324std::vector<std::pair<unsigned int, unsigned int>>
326 const FiniteElement<dim> &fe_other) const
327{
328 // we can presently only compute these identities if both FEs are
329 // FE_RaviartThomasNodals or if the other one is FE_Nothing
330 if (const FE_RaviartThomasNodal<dim> *fe_q_other =
331 dynamic_cast<const FE_RaviartThomasNodal<dim> *>(&fe_other))
332 {
333 // dofs are located on faces; these are only lines in 2d
334 if (dim != 2)
335 return std::vector<std::pair<unsigned int, unsigned int>>();
336
337 // dofs are located along lines, so two dofs are identical only if in
338 // the following two cases (remember that the face support points are
339 // Gauss points):
340 // 1. this->degree = fe_q_other->degree,
341 // in the case, all the dofs on the line are identical
342 // 2. this->degree-1 and fe_q_other->degree-1
343 // are both even, i.e. this->dof_per_line and fe_q_other->dof_per_line
344 // are both odd, there exists only one point (the middle one) such
345 // that dofs are identical on this point
346 //
347 // to understand this, note that this->degree is the *maximal*
348 // polynomial degree, and is thus one higher than the argument given to
349 // the constructor
350 const unsigned int p = this->degree - 1;
351 const unsigned int q = fe_q_other->degree - 1;
352
353 std::vector<std::pair<unsigned int, unsigned int>> identities;
354
355 if (p == q)
356 for (unsigned int i = 0; i < p + 1; ++i)
357 identities.emplace_back(i, i);
358
359 else if (p % 2 == 0 && q % 2 == 0)
360 identities.emplace_back(p / 2, q / 2);
361
362 return identities;
363 }
364 else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
365 {
366 // the FE_Nothing has no degrees of freedom, so there are no
367 // equivalencies to be recorded
368 return std::vector<std::pair<unsigned int, unsigned int>>();
369 }
370 else
371 {
373 return std::vector<std::pair<unsigned int, unsigned int>>();
374 }
375}
376
377
378template <int dim>
379std::vector<std::pair<unsigned int, unsigned int>>
381 const FiniteElement<dim> &fe_other,
382 const unsigned int face_no) const
383{
384 // we can presently only compute these identities if both FEs are
385 // FE_RaviartThomasNodals or if the other one is FE_Nothing
386 if (const FE_RaviartThomasNodal<dim> *fe_q_other =
387 dynamic_cast<const FE_RaviartThomasNodal<dim> *>(&fe_other))
388 {
389 // dofs are located on faces; these are only quads in 3d
390 if (dim != 3)
391 return std::vector<std::pair<unsigned int, unsigned int>>();
392
393 // this works exactly like the line case above
394 const unsigned int p = this->n_dofs_per_quad(face_no);
395
396 AssertDimension(fe_q_other->n_unique_faces(), 1);
397 const unsigned int q = fe_q_other->n_dofs_per_quad(0);
398
399 std::vector<std::pair<unsigned int, unsigned int>> identities;
400
401 if (p == q)
402 for (unsigned int i = 0; i < p; ++i)
403 identities.emplace_back(i, i);
404
405 else if (p % 2 != 0 && q % 2 != 0)
406 identities.emplace_back(p / 2, q / 2);
407
408 return identities;
409 }
410 else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
411 {
412 // the FE_Nothing has no degrees of freedom, so there are no
413 // equivalencies to be recorded
414 return std::vector<std::pair<unsigned int, unsigned int>>();
415 }
416 else
417 {
419 return std::vector<std::pair<unsigned int, unsigned int>>();
420 }
421}
422
423
424template <int dim>
427 const FiniteElement<dim> &fe_other,
428 const unsigned int codim) const
429{
430 Assert(codim <= dim, ExcImpossibleInDim(dim));
431 (void)codim;
432
433 // vertex/line/face/cell domination
434 // --------------------------------
435 if (const FE_RaviartThomasNodal<dim> *fe_rt_nodal_other =
436 dynamic_cast<const FE_RaviartThomasNodal<dim> *>(&fe_other))
437 {
438 if (this->degree < fe_rt_nodal_other->degree)
440 else if (this->degree == fe_rt_nodal_other->degree)
442 else
444 }
445 else if (const FE_Nothing<dim> *fe_nothing =
446 dynamic_cast<const FE_Nothing<dim> *>(&fe_other))
447 {
448 if (fe_nothing->is_dominating())
450 else
451 // the FE_Nothing has no degrees of freedom and it is typically used
452 // in a context where we don't require any continuity along the
453 // interface
455 }
456
459}
460
461
462
463template <>
464void
466 const FiniteElement<1, 1> & /*x_source_fe*/,
467 FullMatrix<double> & /*interpolation_matrix*/,
468 const unsigned int) const
469{
470 Assert(false, ExcImpossibleInDim(1));
471}
472
473
474template <>
475void
477 const FiniteElement<1, 1> & /*x_source_fe*/,
478 const unsigned int /*subface*/,
479 FullMatrix<double> & /*interpolation_matrix*/,
480 const unsigned int) const
481{
482 Assert(false, ExcImpossibleInDim(1));
483}
484
485
486
487template <int dim>
488void
490 const FiniteElement<dim> &x_source_fe,
491 FullMatrix<double> &interpolation_matrix,
492 const unsigned int face_no) const
493{
494 // this is only implemented, if the
495 // source FE is also a
496 // RaviartThomasNodal element
497 AssertThrow((x_source_fe.get_name().find("FE_RaviartThomasNodal<") == 0) ||
498 (dynamic_cast<const FE_RaviartThomasNodal<dim> *>(
499 &x_source_fe) != nullptr),
501
502 Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
503 ExcDimensionMismatch(interpolation_matrix.n(),
504 this->n_dofs_per_face(face_no)));
505 Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
506 ExcDimensionMismatch(interpolation_matrix.m(),
507 x_source_fe.n_dofs_per_face(face_no)));
508
509 // ok, source is a RaviartThomasNodal element, so we will be able to do the
510 // work
511 const FE_RaviartThomasNodal<dim> &source_fe =
512 dynamic_cast<const FE_RaviartThomasNodal<dim> &>(x_source_fe);
513
514 // Make sure that the element for which the DoFs should be constrained is
515 // the one with the higher polynomial degree. Actually the procedure will
516 // work also if this assertion is not satisfied. But the matrices produced
517 // in that case might lead to problems in the hp-procedures, which use this
518 // method.
519 Assert(this->n_dofs_per_face(face_no) <= source_fe.n_dofs_per_face(face_no),
521
522 // generate a quadrature with the generalized support points. This is later
523 // based as a basis for the QProjector, which returns the support points on
524 // the face.
525 Quadrature<dim - 1> quad_face_support(
526 source_fe.generalized_face_support_points[face_no]);
527
528 // Rule of thumb for FP accuracy, that can be expected for a given
529 // polynomial degree. This value is used to cut off values close to zero.
530 double eps = 2e-13 * this->degree * (dim - 1);
531
532 // compute the interpolation matrix by simply taking the value at the
533 // support points.
534 const Quadrature<dim> face_projection =
535 QProjector<dim>::project_to_face(this->reference_cell(),
536 quad_face_support,
537 0);
538
539 for (unsigned int i = 0; i < source_fe.n_dofs_per_face(face_no); ++i)
540 {
541 const Point<dim> &p = face_projection.point(i);
542
543 for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
544 {
545 double matrix_entry =
546 this->shape_value_component(this->face_to_cell_index(j, 0), p, 0);
547
548 // Correct the interpolated value. I.e. if it is close to 1 or 0,
549 // make it exactly 1 or 0. Unfortunately, this is required to avoid
550 // problems with higher order elements.
551 if (std::fabs(matrix_entry - 1.0) < eps)
552 matrix_entry = 1.0;
553 if (std::fabs(matrix_entry) < eps)
554 matrix_entry = 0.0;
555
556 interpolation_matrix(i, j) = matrix_entry;
557 }
558 }
559
560#ifdef DEBUG
561 // make sure that the row sum of each of the matrices is 1 at this
562 // point. this must be so since the shape functions sum up to 1
563 for (unsigned int j = 0; j < source_fe.n_dofs_per_face(face_no); ++j)
564 {
565 double sum = 0.;
566
567 for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
568 sum += interpolation_matrix(j, i);
569
570 Assert(std::fabs(sum - 1) < 2e-13 * this->degree * (dim - 1),
572 }
573#endif
574}
575
576
577template <int dim>
578void
580 const FiniteElement<dim> &x_source_fe,
581 const unsigned int subface,
582 FullMatrix<double> &interpolation_matrix,
583 const unsigned int face_no) const
584{
585 // this is only implemented, if the source FE is also a RaviartThomasNodal
586 // element
587 AssertThrow((x_source_fe.get_name().find("FE_RaviartThomasNodal<") == 0) ||
588 (dynamic_cast<const FE_RaviartThomasNodal<dim> *>(
589 &x_source_fe) != nullptr),
591
592 Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
593 ExcDimensionMismatch(interpolation_matrix.n(),
594 this->n_dofs_per_face(face_no)));
595 Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
596 ExcDimensionMismatch(interpolation_matrix.m(),
597 x_source_fe.n_dofs_per_face(face_no)));
598
599 // ok, source is a RaviartThomasNodal element, so we will be able to do the
600 // work
601 const FE_RaviartThomasNodal<dim> &source_fe =
602 dynamic_cast<const FE_RaviartThomasNodal<dim> &>(x_source_fe);
603
604 // Make sure that the element for which the DoFs should be constrained is
605 // the one with the higher polynomial degree. Actually the procedure will
606 // work also if this assertion is not satisfied. But the matrices produced
607 // in that case might lead to problems in the hp-procedures, which use this
608 // method.
609 Assert(this->n_dofs_per_face(face_no) <= source_fe.n_dofs_per_face(face_no),
611
612 // generate a quadrature with the generalized support points. This is later
613 // based as a basis for the QProjector, which returns the support points on
614 // the face.
615 Quadrature<dim - 1> quad_face_support(
616 source_fe.generalized_face_support_points[face_no]);
617
618 // Rule of thumb for FP accuracy, that can be expected for a given
619 // polynomial degree. This value is used to cut off values close to zero.
620 double eps = 2e-13 * this->degree * (dim - 1);
621
622 // compute the interpolation matrix by simply taking the value at the
623 // support points.
624 const Quadrature<dim> subface_projection =
625 QProjector<dim>::project_to_subface(this->reference_cell(),
626 quad_face_support,
627 0,
628 subface);
629
630 for (unsigned int i = 0; i < source_fe.n_dofs_per_face(face_no); ++i)
631 {
632 const Point<dim> &p = subface_projection.point(i);
633
634 for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
635 {
636 double matrix_entry =
637 this->shape_value_component(this->face_to_cell_index(j, 0), p, 0);
638
639 // Correct the interpolated value. I.e. if it is close to 1 or 0,
640 // make it exactly 1 or 0. Unfortunately, this is required to avoid
641 // problems with higher order elements.
642 if (std::fabs(matrix_entry - 1.0) < eps)
643 matrix_entry = 1.0;
644 if (std::fabs(matrix_entry) < eps)
645 matrix_entry = 0.0;
646
647 interpolation_matrix(i, j) = matrix_entry;
648 }
649 }
650
651#ifdef DEBUG
652 // make sure that the row sum of each of the matrices is 1 at this
653 // point. this must be so since the shape functions sum up to 1
654 for (unsigned int j = 0; j < source_fe.n_dofs_per_face(face_no); ++j)
655 {
656 double sum = 0.;
657
658 for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
659 sum += interpolation_matrix(j, i);
660
661 Assert(std::fabs(sum - 1) < 2e-13 * this->degree * (dim - 1),
663 }
664#endif
665}
666
667
668
669template <int dim>
670const FullMatrix<double> &
672 const unsigned int child,
673 const RefinementCase<dim> &refinement_case) const
674{
675 AssertIndexRange(refinement_case,
679 "Prolongation matrices are only available for refined cells!"));
680 AssertIndexRange(child, GeometryInfo<dim>::n_children(refinement_case));
681
682 // initialization upon first request
683 if (this->prolongation[refinement_case - 1][child].n() == 0)
684 {
685 std::lock_guard<std::mutex> lock(prolongation_matrix_mutex);
686
687 // if matrix got updated while waiting for the lock
688 if (this->prolongation[refinement_case - 1][child].n() ==
689 this->n_dofs_per_cell())
690 return this->prolongation[refinement_case - 1][child];
691
692 // now do the work. need to get a non-const version of data in order to
693 // be able to modify them inside a const function
694 FE_RaviartThomasNodal<dim> &this_nonconst =
695 const_cast<FE_RaviartThomasNodal<dim> &>(*this);
696 if (refinement_case == RefinementCase<dim>::isotropic_refinement)
697 {
698 std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
700 isotropic_matrices.back().resize(
702 FullMatrix<double>(this->n_dofs_per_cell(),
703 this->n_dofs_per_cell()));
704 FETools::compute_embedding_matrices(*this, isotropic_matrices, true);
705 this_nonconst.prolongation[refinement_case - 1] =
706 std::move(isotropic_matrices.back());
707 }
708 else
709 {
710 // must compute both restriction and prolongation matrices because
711 // we only check for their size and the reinit call initializes them
712 // all
715 this_nonconst.prolongation);
717 this_nonconst.restriction);
718 }
719 }
720
721 // finally return the matrix
722 return this->prolongation[refinement_case - 1][child];
723}
724
725
726
727template <int dim>
728const FullMatrix<double> &
730 const unsigned int child,
731 const RefinementCase<dim> &refinement_case) const
732{
733 AssertIndexRange(refinement_case,
737 "Restriction matrices are only available for refined cells!"));
738 AssertIndexRange(child, GeometryInfo<dim>::n_children(refinement_case));
739
740 // initialization upon first request
741 if (this->restriction[refinement_case - 1][child].n() == 0)
742 {
743 std::lock_guard<std::mutex> lock(restriction_matrix_mutex);
744
745 // if matrix got updated while waiting for the lock...
746 if (this->restriction[refinement_case - 1][child].n() ==
747 this->n_dofs_per_cell())
748 return this->restriction[refinement_case - 1][child];
749
750 // now do the work. need to get a non-const version of data in order to
751 // be able to modify them inside a const function
752 FE_RaviartThomasNodal<dim> &this_nonconst =
753 const_cast<FE_RaviartThomasNodal<dim> &>(*this);
754 if (refinement_case == RefinementCase<dim>::isotropic_refinement)
755 {
756 std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
758 isotropic_matrices.back().resize(
760 FullMatrix<double>(this->n_dofs_per_cell(),
761 this->n_dofs_per_cell()));
762 FETools::compute_projection_matrices(*this, isotropic_matrices, true);
763 this_nonconst.restriction[refinement_case - 1] =
764 std::move(isotropic_matrices.back());
765 }
766 else
767 {
768 // must compute both restriction and prolongation matrices because
769 // we only check for their size and the reinit call initializes them
770 // all
773 this_nonconst.prolongation);
775 this_nonconst.restriction);
776 }
777 }
778
779 // finally return the matrix
780 return this->restriction[refinement_case - 1][child];
781}
782
783
784
785// explicit instantiations
786#include "fe_raviart_thomas_nodal.inst"
787
788
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: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:2413
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:2470
FullMatrix< double > interface_constraints
Definition fe.h:2439
std::vector< Point< dim > > generalized_support_points
Definition fe.h:2464
std::vector< std::vector< FullMatrix< double > > > prolongation
Definition fe.h:2427
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)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:503
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:504
#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:134
#define DEAL_II_NOT_IMPLEMENTED()
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_face_embedding_matrices(const FiniteElement< dim, spacedim > &fe, FullMatrix< number >(&matrices)[GeometryInfo< dim >::max_children_per_face], const unsigned int face_coarse, const unsigned int face_fine, 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)
unsigned char combined_face_orientation(const bool face_orientation, const bool face_rotation, const bool face_flip)
STL namespace.