Reference documentation for deal.II version 9.4.0
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>
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 template <int dim>
42 class PolynomialsRaviartThomasNodal : public TensorPolynomialsBase<dim>
43 {
44 public:
45 PolynomialsRaviartThomasNodal(const unsigned int degree);
46
54 void
55 evaluate(const Point<dim> & unit_point,
56 std::vector<Tensor<1, dim>> &values,
57 std::vector<Tensor<2, dim>> &grads,
58 std::vector<Tensor<3, dim>> &grad_grads,
59 std::vector<Tensor<4, dim>> &third_derivatives,
60 std::vector<Tensor<5, dim>> &fourth_derivatives) const override;
61
65 std::string
66 name() const override;
67
73 static unsigned int
74 n_polynomials(const unsigned int degree);
75
76 const std::vector<unsigned int> &
77 get_renumbering() const;
78
82 virtual std::unique_ptr<TensorPolynomialsBase<dim>>
83 clone() const override;
84
92 std::vector<Point<dim>>
93 get_polynomial_support_points() const;
94
95 private:
99 const unsigned int degree;
100
105 const AnisotropicPolynomials<dim> polynomial_space;
106
110 std::vector<unsigned int> lexicographic_to_hierarchic;
111
116 std::vector<unsigned int> hierarchic_to_lexicographic;
117
121 std::array<std::vector<unsigned int>, dim> renumber_aniso;
122 };
123
124
125
126 // Create nodal Raviart-Thomas polynomials as the tensor product of Lagrange
127 // polynomials on Gauss-Lobatto points of degree + 2 points in the
128 // continuous direction and degree + 1 points in the discontinuous
129 // directions (we could also choose Lagrange polynomials on Gauss points but
130 // those are slightly more expensive to handle in classes).
131 std::vector<std::vector<Polynomials::Polynomial<double>>>
132 create_rt_polynomials(const unsigned int dim, const unsigned int degree)
133 {
134 std::vector<std::vector<Polynomials::Polynomial<double>>> pols(dim);
136 QGaussLobatto<1>(degree + 2).get_points());
137 if (degree > 0)
138 for (unsigned int d = 1; d < dim; ++d)
140 QGaussLobatto<1>(degree + 1).get_points());
141 else
142 for (unsigned int d = 1; d < dim; ++d)
144 QMidpoint<1>().get_points());
145
146 return pols;
147 }
148
149
150 template <int dim>
151 PolynomialsRaviartThomasNodal<dim>::PolynomialsRaviartThomasNodal(
152 const unsigned int degree)
153 : TensorPolynomialsBase<dim>(degree, n_polynomials(degree))
154 , degree(degree)
155 , polynomial_space(create_rt_polynomials(dim, degree))
156 {
157 // create renumbering of the unknowns from the lexicographic order to the
158 // actual order required by the finite element class with unknowns on
159 // faces placed first
160 const unsigned int n_pols = polynomial_space.n();
161 lexicographic_to_hierarchic =
162 internal::get_lexicographic_numbering_rt_nodal<dim>(degree + 1);
163
164 hierarchic_to_lexicographic =
165 Utilities::invert_permutation(lexicographic_to_hierarchic);
166
167 // since we only store an anisotropic polynomial for the first component,
168 // we set up a second numbering to switch out the actual coordinate
169 // direction
170 renumber_aniso[0].resize(n_pols);
171 for (unsigned int i = 0; i < n_pols; ++i)
172 renumber_aniso[0][i] = i;
173 if (dim == 2)
174 {
175 // switch x and y component (i and j loops)
176 renumber_aniso[1].resize(n_pols);
177 for (unsigned int j = 0; j < degree + 2; ++j)
178 for (unsigned int i = 0; i < degree + 1; ++i)
179 renumber_aniso[1][j * (degree + 1) + i] = j + i * (degree + 2);
180 }
181 if (dim == 3)
182 {
183 // switch x, y, and z component (i, j, k) -> (j, k, i)
184 renumber_aniso[1].resize(n_pols);
185 for (unsigned int k = 0; k < degree + 1; ++k)
186 for (unsigned int j = 0; j < degree + 2; ++j)
187 for (unsigned int i = 0; i < degree + 1; ++i)
188 renumber_aniso[1][(k * (degree + 2) + j) * (degree + 1) + i] =
189 j + k * (degree + 2) + i * (degree + 2) * (degree + 1);
190
191 // switch x, y, and z component (i, j, k) -> (k, i, j)
192 renumber_aniso[2].resize(n_pols);
193 for (unsigned int k = 0; k < degree + 2; ++k)
194 for (unsigned int j = 0; j < degree + 1; ++j)
195 for (unsigned int i = 0; i < degree + 1; ++i)
196 renumber_aniso[2][(k * (degree + 1) + j) * (degree + 1) + i] =
197 k + i * (degree + 2) + j * (degree + 2) * (degree + 1);
198 }
199 }
200
201
202
203 template <int dim>
204 void
205 PolynomialsRaviartThomasNodal<dim>::evaluate(
206 const Point<dim> & unit_point,
207 std::vector<Tensor<1, dim>> &values,
208 std::vector<Tensor<2, dim>> &grads,
209 std::vector<Tensor<3, dim>> &grad_grads,
210 std::vector<Tensor<4, dim>> &third_derivatives,
211 std::vector<Tensor<5, dim>> &fourth_derivatives) const
212 {
213 Assert(values.size() == this->n() || values.size() == 0,
214 ExcDimensionMismatch(values.size(), this->n()));
215 Assert(grads.size() == this->n() || grads.size() == 0,
216 ExcDimensionMismatch(grads.size(), this->n()));
217 Assert(grad_grads.size() == this->n() || grad_grads.size() == 0,
218 ExcDimensionMismatch(grad_grads.size(), this->n()));
219 Assert(third_derivatives.size() == this->n() ||
220 third_derivatives.size() == 0,
221 ExcDimensionMismatch(third_derivatives.size(), this->n()));
222 Assert(fourth_derivatives.size() == this->n() ||
223 fourth_derivatives.size() == 0,
224 ExcDimensionMismatch(fourth_derivatives.size(), this->n()));
225
226 std::vector<double> p_values;
227 std::vector<Tensor<1, dim>> p_grads;
228 std::vector<Tensor<2, dim>> p_grad_grads;
229 std::vector<Tensor<3, dim>> p_third_derivatives;
230 std::vector<Tensor<4, dim>> p_fourth_derivatives;
231
232 const unsigned int n_sub = polynomial_space.n();
233 p_values.resize((values.size() == 0) ? 0 : n_sub);
234 p_grads.resize((grads.size() == 0) ? 0 : n_sub);
235 p_grad_grads.resize((grad_grads.size() == 0) ? 0 : n_sub);
236 p_third_derivatives.resize((third_derivatives.size() == 0) ? 0 : n_sub);
237 p_fourth_derivatives.resize((fourth_derivatives.size() == 0) ? 0 : n_sub);
238
239 for (unsigned int d = 0; d < dim; ++d)
240 {
241 // First we copy the point. The polynomial space for component d
242 // consists of polynomials of degree k in x_d and degree k+1 in the
243 // other variables. in order to simplify this, we use the same
244 // AnisotropicPolynomial space and simply rotate the coordinates
245 // through all directions.
246 Point<dim> p;
247 for (unsigned int c = 0; c < dim; ++c)
248 p(c) = unit_point((c + d) % dim);
249
250 polynomial_space.evaluate(p,
251 p_values,
252 p_grads,
253 p_grad_grads,
254 p_third_derivatives,
255 p_fourth_derivatives);
256
257 for (unsigned int i = 0; i < p_values.size(); ++i)
258 values[lexicographic_to_hierarchic[i + d * n_sub]][d] =
259 p_values[renumber_aniso[d][i]];
260
261 for (unsigned int i = 0; i < p_grads.size(); ++i)
262 for (unsigned int d1 = 0; d1 < dim; ++d1)
263 grads[lexicographic_to_hierarchic[i + d * n_sub]][d]
264 [(d1 + d) % dim] = p_grads[renumber_aniso[d][i]][d1];
265
266 for (unsigned int i = 0; i < p_grad_grads.size(); ++i)
267 for (unsigned int d1 = 0; d1 < dim; ++d1)
268 for (unsigned int d2 = 0; d2 < dim; ++d2)
269 grad_grads[lexicographic_to_hierarchic[i + d * n_sub]][d]
270 [(d1 + d) % dim][(d2 + d) % dim] =
271 p_grad_grads[renumber_aniso[d][i]][d1][d2];
272
273 for (unsigned int i = 0; i < p_third_derivatives.size(); ++i)
274 for (unsigned int d1 = 0; d1 < dim; ++d1)
275 for (unsigned int d2 = 0; d2 < dim; ++d2)
276 for (unsigned int d3 = 0; d3 < dim; ++d3)
277 third_derivatives[lexicographic_to_hierarchic[i + d * n_sub]][d]
278 [(d1 + d) % dim][(d2 + d) % dim]
279 [(d3 + d) % dim] =
280 p_third_derivatives[renumber_aniso[d][i]][d1]
281 [d2][d3];
282
283 for (unsigned int i = 0; i < p_fourth_derivatives.size(); ++i)
284 for (unsigned int d1 = 0; d1 < dim; ++d1)
285 for (unsigned int d2 = 0; d2 < dim; ++d2)
286 for (unsigned int d3 = 0; d3 < dim; ++d3)
287 for (unsigned int d4 = 0; d4 < dim; ++d4)
288 fourth_derivatives[lexicographic_to_hierarchic[i + d * n_sub]]
289 [d][(d1 + d) % dim][(d2 + d) % dim]
290 [(d3 + d) % dim][(d4 + d) % dim] =
291 p_fourth_derivatives[renumber_aniso[d][i]]
292 [d1][d2][d3][d4];
293 }
294 }
295
296
297
298 template <int dim>
299 std::string
300 PolynomialsRaviartThomasNodal<dim>::name() const
301 {
302 return "PolynomialsRaviartThomasNodal";
303 }
304
305
306
307 template <int dim>
308 unsigned int
309 PolynomialsRaviartThomasNodal<dim>::n_polynomials(unsigned int degree)
310 {
311 return dim * (degree + 2) * Utilities::pow(degree + 1, dim - 1);
312 }
313
314
315
316 template <int dim>
317 std::unique_ptr<TensorPolynomialsBase<dim>>
318 PolynomialsRaviartThomasNodal<dim>::clone() const
319 {
320 return std::make_unique<PolynomialsRaviartThomasNodal<dim>>(*this);
321 }
322
323
324
325 template <int dim>
326 std::vector<Point<dim>>
327 PolynomialsRaviartThomasNodal<dim>::get_polynomial_support_points() const
328 {
329 Assert(dim > 0 && dim <= 3, ExcImpossibleInDim(dim));
330 const Quadrature<1> low(
331 degree == 0 ? static_cast<Quadrature<1>>(QMidpoint<1>()) :
332 static_cast<Quadrature<1>>(QGaussLobatto<1>(degree + 1)));
333 const QGaussLobatto<1> high(degree + 2);
334 const QAnisotropic<dim> quad =
335 (dim == 1 ? QAnisotropic<dim>(high) :
336 (dim == 2 ? QAnisotropic<dim>(high, low) :
337 QAnisotropic<dim>(high, low, low)));
338
339 const unsigned int n_sub = polynomial_space.n();
340 std::vector<Point<dim>> points(dim * n_sub);
341 points.resize(n_polynomials(degree));
342 for (unsigned int d = 0; d < dim; ++d)
343 for (unsigned int i = 0; i < n_sub; ++i)
344 for (unsigned int c = 0; c < dim; ++c)
345 points[lexicographic_to_hierarchic[i + d * n_sub]][(c + d) % dim] =
346 quad.point(renumber_aniso[d][i])[c];
347 return points;
348 }
349
350
351
352 // Return a vector of "dofs per object" where the components of the returned
353 // vector refer to:
354 // 0 = vertex
355 // 1 = edge
356 // 2 = face (which is a cell in 2D)
357 // 3 = cell
358 std::vector<unsigned int>
359 get_rt_dpo_vector(const unsigned int dim, const unsigned int degree)
360 {
361 std::vector<unsigned int> dpo(dim + 1);
362 dpo[0] = 0;
363 dpo[1] = 0;
364 unsigned int dofs_per_face = 1;
365 for (unsigned int d = 1; d < dim; ++d)
366 dofs_per_face *= (degree + 1);
367
368 dpo[dim - 1] = dofs_per_face;
369 dpo[dim] = dim * degree * dofs_per_face;
370
371 return dpo;
372 }
373} // namespace
374
375namespace internal
376{
377 template <int dim>
378 std::vector<unsigned int>
379 get_lexicographic_numbering_rt_nodal(const unsigned int degree)
380 {
381 const unsigned int n_dofs_face = Utilities::pow(degree, dim - 1);
382 std::vector<unsigned int> lexicographic_numbering;
383 // component 1
384 for (unsigned int j = 0; j < n_dofs_face; j++)
385 {
386 lexicographic_numbering.push_back(j);
387 for (unsigned int i = n_dofs_face * 2 * dim;
388 i < n_dofs_face * 2 * dim + degree - 1;
389 i++)
390 lexicographic_numbering.push_back(i + j * (degree - 1));
391 lexicographic_numbering.push_back(n_dofs_face + j);
392 }
393
394 // component 2
395 unsigned int layers = (dim == 3) ? degree : 1;
396 for (unsigned int k = 0; k < layers; k++)
397 {
398 unsigned int k_add = k * degree;
399 for (unsigned int j = n_dofs_face * 2; j < n_dofs_face * 2 + degree;
400 j++)
401 lexicographic_numbering.push_back(j + k_add);
402
403 for (unsigned int i = n_dofs_face * (2 * dim + (degree - 1));
404 i < n_dofs_face * (2 * dim + (degree - 1)) + (degree - 1) * degree;
405 i++)
406 {
407 lexicographic_numbering.push_back(i + k_add * (degree - 1));
408 }
409 for (unsigned int j = n_dofs_face * 3; j < n_dofs_face * 3 + degree;
410 j++)
411 lexicographic_numbering.push_back(j + k_add);
412 }
413
414 // component 3
415 if (dim == 3)
416 {
417 for (unsigned int i = 4 * n_dofs_face; i < 5 * n_dofs_face; i++)
418 lexicographic_numbering.push_back(i);
419 for (unsigned int i = 6 * n_dofs_face + n_dofs_face * 2 * (degree - 1);
420 i < 6 * n_dofs_face + n_dofs_face * 3 * (degree - 1);
421 i++)
422 lexicographic_numbering.push_back(i);
423 for (unsigned int i = 5 * n_dofs_face; i < 6 * n_dofs_face; i++)
424 lexicographic_numbering.push_back(i);
425 }
426
427 return lexicographic_numbering;
428 }
429} // namespace internal
430
431
432
433// --------------------- actual implementation of element --------------------
434
435template <int dim>
437 : FE_PolyTensor<dim>(PolynomialsRaviartThomasNodal<dim>(degree),
438 FiniteElementData<dim>(get_rt_dpo_vector(dim, degree),
439 dim,
440 degree + 1,
441 FiniteElementData<dim>::Hdiv),
442 std::vector<bool>(1, false),
443 std::vector<ComponentMask>(
444 PolynomialsRaviartThomasNodal<dim>::n_polynomials(
445 degree),
446 std::vector<bool>(dim, true)))
447{
448 Assert(dim >= 2, ExcImpossibleInDim(dim));
449
451
452 // First, initialize the generalized support points and quadrature weights,
453 // since they are required for interpolation.
455 PolynomialsRaviartThomasNodal<dim>(degree).get_polynomial_support_points();
457 this->n_dofs_per_cell());
458
459 const unsigned int face_no = 0;
460 if (dim > 1)
461 this->generalized_face_support_points[face_no] =
462 degree == 0 ? QGauss<dim - 1>(1).get_points() :
463 QGaussLobatto<dim - 1>(degree + 1).get_points();
464
466 for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_face; ++i)
467 face_embeddings[i].reinit(this->n_dofs_per_face(face_no),
468 this->n_dofs_per_face(face_no));
469 FETools::compute_face_embedding_matrices<dim, double>(*this,
470 face_embeddings,
471 0,
472 0);
474 this->n_dofs_per_face(face_no),
475 this->n_dofs_per_face(face_no));
476 unsigned int target_row = 0;
477 for (unsigned int d = 0; d < GeometryInfo<dim>::max_children_per_face; ++d)
478 for (unsigned int i = 0; i < face_embeddings[d].m(); ++i)
479 {
480 for (unsigned int j = 0; j < face_embeddings[d].n(); ++j)
481 this->interface_constraints(target_row, j) = face_embeddings[d](i, j);
482 ++target_row;
483 }
484
485 // We need to initialize the dof permutation table and the one for the sign
486 // change.
488}
489
490
491
492template <int dim>
493std::string
495{
496 // note that the FETools::get_fe_by_name function depends on the particular
497 // format of the string this function returns, so they have to be kept in
498 // synch
499
500 // note that this->degree is the maximal polynomial degree and is thus one
501 // higher than the argument given to the constructor
502 return "FE_RaviartThomasNodal<" + std::to_string(dim) + ">(" +
503 std::to_string(this->degree - 1) + ")";
504}
505
506
507template <int dim>
508std::unique_ptr<FiniteElement<dim, dim>>
510{
511 return std::make_unique<FE_RaviartThomasNodal<dim>>(*this);
512}
513
514
515//---------------------------------------------------------------------------
516// Auxiliary and internal functions
517//---------------------------------------------------------------------------
518
519
520
521template <int dim>
522void
524 dim>::initialize_quad_dof_index_permutation_and_sign_change()
525{
526 // for 1D and 2D, do nothing
527 if (dim < 3)
528 return;
529
530 const unsigned int n = this->degree;
531 const unsigned int face_no = 0;
532 Assert(n * n == this->n_dofs_per_quad(face_no), ExcInternalError());
533 for (unsigned int local = 0; local < this->n_dofs_per_quad(face_no); ++local)
534 // face support points are in lexicographic ordering with x running
535 // fastest. invert that (y running fastest)
536 {
537 unsigned int i = local % n, j = local / n;
538
539 // face_orientation=false, face_flip=false, face_rotation=false
540 this->adjust_quad_dof_index_for_face_orientation_table[face_no](local,
541 0) =
542 j + i * n - local;
543 // face_orientation=false, face_flip=false, face_rotation=true
544 this->adjust_quad_dof_index_for_face_orientation_table[face_no](local,
545 1) =
546 i + (n - 1 - j) * n - local;
547 // face_orientation=false, face_flip=true, face_rotation=false
548 this->adjust_quad_dof_index_for_face_orientation_table[face_no](local,
549 2) =
550 (n - 1 - j) + (n - 1 - i) * n - local;
551 // face_orientation=false, face_flip=true, face_rotation=true
552 this->adjust_quad_dof_index_for_face_orientation_table[face_no](local,
553 3) =
554 (n - 1 - i) + j * n - local;
555 // face_orientation=true, face_flip=false, face_rotation=false
556 this->adjust_quad_dof_index_for_face_orientation_table[face_no](local,
557 4) = 0;
558 // face_orientation=true, face_flip=false, face_rotation=true
559 this->adjust_quad_dof_index_for_face_orientation_table[face_no](local,
560 5) =
561 j + (n - 1 - i) * n - local;
562 // face_orientation=true, face_flip=true, face_rotation=false
563 this->adjust_quad_dof_index_for_face_orientation_table[face_no](local,
564 6) =
565 (n - 1 - i) + (n - 1 - j) * n - local;
566 // face_orientation=true, face_flip=true, face_rotation=true
567 this->adjust_quad_dof_index_for_face_orientation_table[face_no](local,
568 7) =
569 (n - 1 - j) + i * n - local;
570
571 // for face_orientation == false, we need to switch the sign
572 for (unsigned int i = 0; i < 4; ++i)
573 this->adjust_quad_dof_sign_for_face_orientation_table[face_no](local,
574 i) = 1;
575 }
576}
577
578
579
580template <int dim>
581bool
583 const unsigned int shape_index,
584 const unsigned int face_index) const
585{
586 AssertIndexRange(shape_index, this->n_dofs_per_cell());
588
589 // The first degrees of freedom are on the faces and each face has degree
590 // degrees.
591 const unsigned int support_face = shape_index / this->n_dofs_per_face();
592
593 // The only thing we know for sure is that shape functions with support on
594 // one face are zero on the opposite face.
595 if (support_face < GeometryInfo<dim>::faces_per_cell)
596 return (face_index != GeometryInfo<dim>::opposite_face[support_face]);
597
598 // In all other cases, return true, which is safe
599 return true;
600}
601
602
603
604template <int dim>
605void
608 const std::vector<Vector<double>> &support_point_values,
609 std::vector<double> & nodal_values) const
610{
611 Assert(support_point_values.size() == this->generalized_support_points.size(),
612 ExcDimensionMismatch(support_point_values.size(),
613 this->generalized_support_points.size()));
614 Assert(nodal_values.size() == this->n_dofs_per_cell(),
615 ExcDimensionMismatch(nodal_values.size(), this->n_dofs_per_cell()));
616 Assert(support_point_values[0].size() == this->n_components(),
617 ExcDimensionMismatch(support_point_values[0].size(),
618 this->n_components()));
619
620 // First do interpolation on faces. There, the component evaluated depends
621 // on the face direction and orientation.
622 unsigned int fbase = 0;
623 unsigned int f = 0;
624 for (; f < GeometryInfo<dim>::faces_per_cell;
625 ++f, fbase += this->n_dofs_per_face(f))
626 {
627 for (unsigned int i = 0; i < this->n_dofs_per_face(f); ++i)
628 {
629 nodal_values[fbase + i] = support_point_values[fbase + i](
631 }
632 }
633
634 // The remaining points form dim chunks, one for each component
635 const unsigned int istep = (this->n_dofs_per_cell() - fbase) / dim;
636 Assert((this->n_dofs_per_cell() - fbase) % dim == 0, ExcInternalError());
637
638 f = 0;
639 while (fbase < this->n_dofs_per_cell())
640 {
641 for (unsigned int i = 0; i < istep; ++i)
642 {
643 nodal_values[fbase + i] = support_point_values[fbase + i](f);
644 }
645 fbase += istep;
646 ++f;
647 }
648 Assert(fbase == this->n_dofs_per_cell(), ExcInternalError());
649}
650
651
652
653// TODO: There are tests that check that the following few functions don't
654// produce assertion failures, but none that actually check whether they do the
655// right thing. one example for such a test would be to project a function onto
656// an hp-space and make sure that the convergence order is correct with regard
657// to the lowest used polynomial degree
658
659template <int dim>
660bool
662{
663 return true;
664}
665
666
667template <int dim>
668std::vector<std::pair<unsigned int, unsigned int>>
670 const FiniteElement<dim> &fe_other) const
671{
672 // we can presently only compute these identities if both FEs are
673 // FE_RaviartThomasNodals or the other is FE_Nothing. In either case, no
674 // dofs are assigned on the vertex, so we shouldn't be getting here at all.
675 if (dynamic_cast<const FE_RaviartThomasNodal<dim> *>(&fe_other) != nullptr)
676 return std::vector<std::pair<unsigned int, unsigned int>>();
677 else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
678 return std::vector<std::pair<unsigned int, unsigned int>>();
679 else
680 {
681 Assert(false, ExcNotImplemented());
682 return std::vector<std::pair<unsigned int, unsigned int>>();
683 }
684}
685
686
687
688template <int dim>
689std::vector<std::pair<unsigned int, unsigned int>>
691 const FiniteElement<dim> &fe_other) const
692{
693 // we can presently only compute these identities if both FEs are
694 // FE_RaviartThomasNodals or if the other one is FE_Nothing
695 if (const FE_RaviartThomasNodal<dim> *fe_q_other =
696 dynamic_cast<const FE_RaviartThomasNodal<dim> *>(&fe_other))
697 {
698 // dofs are located on faces; these are only lines in 2d
699 if (dim != 2)
700 return std::vector<std::pair<unsigned int, unsigned int>>();
701
702 // dofs are located along lines, so two dofs are identical only if in
703 // the following two cases (remember that the face support points are
704 // Gauss points):
705 // 1. this->degree = fe_q_other->degree,
706 // in the case, all the dofs on the line are identical
707 // 2. this->degree-1 and fe_q_other->degree-1
708 // are both even, i.e. this->dof_per_line and fe_q_other->dof_per_line
709 // are both odd, there exists only one point (the middle one) such
710 // that dofs are identical on this point
711 //
712 // to understand this, note that this->degree is the *maximal*
713 // polynomial degree, and is thus one higher than the argument given to
714 // the constructor
715 const unsigned int p = this->degree - 1;
716 const unsigned int q = fe_q_other->degree - 1;
717
718 std::vector<std::pair<unsigned int, unsigned int>> identities;
719
720 if (p == q)
721 for (unsigned int i = 0; i < p + 1; ++i)
722 identities.emplace_back(i, i);
723
724 else if (p % 2 == 0 && q % 2 == 0)
725 identities.emplace_back(p / 2, q / 2);
726
727 return identities;
728 }
729 else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
730 {
731 // the FE_Nothing has no degrees of freedom, so there are no
732 // equivalencies to be recorded
733 return std::vector<std::pair<unsigned int, unsigned int>>();
734 }
735 else
736 {
737 Assert(false, ExcNotImplemented());
738 return std::vector<std::pair<unsigned int, unsigned int>>();
739 }
740}
741
742
743template <int dim>
744std::vector<std::pair<unsigned int, unsigned int>>
746 const FiniteElement<dim> &fe_other,
747 const unsigned int face_no) const
748{
749 // we can presently only compute these identities if both FEs are
750 // FE_RaviartThomasNodals or if the other one is FE_Nothing
751 if (const FE_RaviartThomasNodal<dim> *fe_q_other =
752 dynamic_cast<const FE_RaviartThomasNodal<dim> *>(&fe_other))
753 {
754 // dofs are located on faces; these are only quads in 3d
755 if (dim != 3)
756 return std::vector<std::pair<unsigned int, unsigned int>>();
757
758 // this works exactly like the line case above
759 const unsigned int p = this->n_dofs_per_quad(face_no);
760
761 AssertDimension(fe_q_other->n_unique_faces(), 1);
762 const unsigned int q = fe_q_other->n_dofs_per_quad(0);
763
764 std::vector<std::pair<unsigned int, unsigned int>> identities;
765
766 if (p == q)
767 for (unsigned int i = 0; i < p; ++i)
768 identities.emplace_back(i, i);
769
770 else if (p % 2 != 0 && q % 2 != 0)
771 identities.emplace_back(p / 2, q / 2);
772
773 return identities;
774 }
775 else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
776 {
777 // the FE_Nothing has no degrees of freedom, so there are no
778 // equivalencies to be recorded
779 return std::vector<std::pair<unsigned int, unsigned int>>();
780 }
781 else
782 {
783 Assert(false, ExcNotImplemented());
784 return std::vector<std::pair<unsigned int, unsigned int>>();
785 }
786}
787
788
789template <int dim>
792 const FiniteElement<dim> &fe_other,
793 const unsigned int codim) const
794{
795 Assert(codim <= dim, ExcImpossibleInDim(dim));
796 (void)codim;
797
798 // vertex/line/face/cell domination
799 // --------------------------------
800 if (const FE_RaviartThomasNodal<dim> *fe_rt_nodal_other =
801 dynamic_cast<const FE_RaviartThomasNodal<dim> *>(&fe_other))
802 {
803 if (this->degree < fe_rt_nodal_other->degree)
805 else if (this->degree == fe_rt_nodal_other->degree)
807 else
809 }
810 else if (const FE_Nothing<dim> *fe_nothing =
811 dynamic_cast<const FE_Nothing<dim> *>(&fe_other))
812 {
813 if (fe_nothing->is_dominating())
815 else
816 // the FE_Nothing has no degrees of freedom and it is typically used
817 // in a context where we don't require any continuity along the
818 // interface
820 }
821
822 Assert(false, ExcNotImplemented());
824}
825
826
827
828template <>
829void
831 const FiniteElement<1, 1> & /*x_source_fe*/,
832 FullMatrix<double> & /*interpolation_matrix*/,
833 const unsigned int) const
834{
835 Assert(false, ExcImpossibleInDim(1));
836}
837
838
839template <>
840void
842 const FiniteElement<1, 1> & /*x_source_fe*/,
843 const unsigned int /*subface*/,
844 FullMatrix<double> & /*interpolation_matrix*/,
845 const unsigned int) const
846{
847 Assert(false, ExcImpossibleInDim(1));
848}
849
850
851
852template <int dim>
853void
855 const FiniteElement<dim> &x_source_fe,
856 FullMatrix<double> & interpolation_matrix,
857 const unsigned int face_no) const
858{
859 // this is only implemented, if the
860 // source FE is also a
861 // RaviartThomasNodal element
862 AssertThrow((x_source_fe.get_name().find("FE_RaviartThomasNodal<") == 0) ||
863 (dynamic_cast<const FE_RaviartThomasNodal<dim> *>(
864 &x_source_fe) != nullptr),
866
867 Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
868 ExcDimensionMismatch(interpolation_matrix.n(),
869 this->n_dofs_per_face(face_no)));
870 Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
871 ExcDimensionMismatch(interpolation_matrix.m(),
872 x_source_fe.n_dofs_per_face(face_no)));
873
874 // ok, source is a RaviartThomasNodal element, so we will be able to do the
875 // work
876 const FE_RaviartThomasNodal<dim> &source_fe =
877 dynamic_cast<const FE_RaviartThomasNodal<dim> &>(x_source_fe);
878
879 // Make sure that the element for which the DoFs should be constrained is
880 // the one with the higher polynomial degree. Actually the procedure will
881 // work also if this assertion is not satisfied. But the matrices produced
882 // in that case might lead to problems in the hp-procedures, which use this
883 // method.
884 Assert(this->n_dofs_per_face(face_no) <= source_fe.n_dofs_per_face(face_no),
886
887 // generate a quadrature with the generalized support points. This is later
888 // based as a basis for the QProjector, which returns the support points on
889 // the face.
890 Quadrature<dim - 1> quad_face_support(
891 source_fe.generalized_face_support_points[face_no]);
892
893 // Rule of thumb for FP accuracy, that can be expected for a given
894 // polynomial degree. This value is used to cut off values close to zero.
895 double eps = 2e-13 * this->degree * (dim - 1);
896
897 // compute the interpolation matrix by simply taking the value at the
898 // support points.
899 const Quadrature<dim> face_projection =
901 quad_face_support,
902 0);
903
904 for (unsigned int i = 0; i < source_fe.n_dofs_per_face(face_no); ++i)
905 {
906 const Point<dim> &p = face_projection.point(i);
907
908 for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
909 {
910 double matrix_entry =
911 this->shape_value_component(this->face_to_cell_index(j, 0), p, 0);
912
913 // Correct the interpolated value. I.e. if it is close to 1 or 0,
914 // make it exactly 1 or 0. Unfortunately, this is required to avoid
915 // problems with higher order elements.
916 if (std::fabs(matrix_entry - 1.0) < eps)
917 matrix_entry = 1.0;
918 if (std::fabs(matrix_entry) < eps)
919 matrix_entry = 0.0;
920
921 interpolation_matrix(i, j) = matrix_entry;
922 }
923 }
924
925#ifdef DEBUG
926 // make sure that the row sum of each of the matrices is 1 at this
927 // point. this must be so since the shape functions sum up to 1
928 for (unsigned int j = 0; j < source_fe.n_dofs_per_face(face_no); ++j)
929 {
930 double sum = 0.;
931
932 for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
933 sum += interpolation_matrix(j, i);
934
935 Assert(std::fabs(sum - 1) < 2e-13 * this->degree * (dim - 1),
937 }
938#endif
939}
940
941
942template <int dim>
943void
945 const FiniteElement<dim> &x_source_fe,
946 const unsigned int subface,
947 FullMatrix<double> & interpolation_matrix,
948 const unsigned int face_no) const
949{
950 // this is only implemented, if the source FE is also a RaviartThomasNodal
951 // element
952 AssertThrow((x_source_fe.get_name().find("FE_RaviartThomasNodal<") == 0) ||
953 (dynamic_cast<const FE_RaviartThomasNodal<dim> *>(
954 &x_source_fe) != nullptr),
956
957 Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
958 ExcDimensionMismatch(interpolation_matrix.n(),
959 this->n_dofs_per_face(face_no)));
960 Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
961 ExcDimensionMismatch(interpolation_matrix.m(),
962 x_source_fe.n_dofs_per_face(face_no)));
963
964 // ok, source is a RaviartThomasNodal element, so we will be able to do the
965 // work
966 const FE_RaviartThomasNodal<dim> &source_fe =
967 dynamic_cast<const FE_RaviartThomasNodal<dim> &>(x_source_fe);
968
969 // Make sure that the element for which the DoFs should be constrained is
970 // the one with the higher polynomial degree. Actually the procedure will
971 // work also if this assertion is not satisfied. But the matrices produced
972 // in that case might lead to problems in the hp-procedures, which use this
973 // method.
974 Assert(this->n_dofs_per_face(face_no) <= source_fe.n_dofs_per_face(face_no),
976
977 // generate a quadrature with the generalized support points. This is later
978 // based as a basis for the QProjector, which returns the support points on
979 // the face.
980 Quadrature<dim - 1> quad_face_support(
981 source_fe.generalized_face_support_points[face_no]);
982
983 // Rule of thumb for FP accuracy, that can be expected for a given
984 // polynomial degree. This value is used to cut off values close to zero.
985 double eps = 2e-13 * this->degree * (dim - 1);
986
987 // compute the interpolation matrix by simply taking the value at the
988 // support points.
989 const Quadrature<dim> subface_projection =
991 quad_face_support,
992 0,
993 subface);
994
995 for (unsigned int i = 0; i < source_fe.n_dofs_per_face(face_no); ++i)
996 {
997 const Point<dim> &p = subface_projection.point(i);
998
999 for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
1000 {
1001 double matrix_entry =
1002 this->shape_value_component(this->face_to_cell_index(j, 0), p, 0);
1003
1004 // Correct the interpolated value. I.e. if it is close to 1 or 0,
1005 // make it exactly 1 or 0. Unfortunately, this is required to avoid
1006 // problems with higher order elements.
1007 if (std::fabs(matrix_entry - 1.0) < eps)
1008 matrix_entry = 1.0;
1009 if (std::fabs(matrix_entry) < eps)
1010 matrix_entry = 0.0;
1011
1012 interpolation_matrix(i, j) = matrix_entry;
1013 }
1014 }
1015
1016#ifdef DEBUG
1017 // make sure that the row sum of each of the matrices is 1 at this
1018 // point. this must be so since the shape functions sum up to 1
1019 for (unsigned int j = 0; j < source_fe.n_dofs_per_face(face_no); ++j)
1020 {
1021 double sum = 0.;
1022
1023 for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
1024 sum += interpolation_matrix(j, i);
1025
1026 Assert(std::fabs(sum - 1) < 2e-13 * this->degree * (dim - 1),
1028 }
1029#endif
1030}
1031
1032
1033
1034template <int dim>
1035const FullMatrix<double> &
1037 const unsigned int child,
1038 const RefinementCase<dim> &refinement_case) const
1039{
1040 AssertIndexRange(refinement_case,
1042 Assert(refinement_case != RefinementCase<dim>::no_refinement,
1043 ExcMessage(
1044 "Prolongation matrices are only available for refined cells!"));
1045 AssertIndexRange(child, GeometryInfo<dim>::n_children(refinement_case));
1046
1047 // initialization upon first request
1048 if (this->prolongation[refinement_case - 1][child].n() == 0)
1049 {
1050 std::lock_guard<std::mutex> lock(this->mutex);
1051
1052 // if matrix got updated while waiting for the lock
1053 if (this->prolongation[refinement_case - 1][child].n() ==
1054 this->n_dofs_per_cell())
1055 return this->prolongation[refinement_case - 1][child];
1056
1057 // now do the work. need to get a non-const version of data in order to
1058 // be able to modify them inside a const function
1059 FE_RaviartThomasNodal<dim> &this_nonconst =
1060 const_cast<FE_RaviartThomasNodal<dim> &>(*this);
1061 if (refinement_case == RefinementCase<dim>::isotropic_refinement)
1062 {
1063 std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
1065 isotropic_matrices.back().resize(
1067 FullMatrix<double>(this->n_dofs_per_cell(),
1068 this->n_dofs_per_cell()));
1069 FETools::compute_embedding_matrices(*this, isotropic_matrices, true);
1070 this_nonconst.prolongation[refinement_case - 1].swap(
1071 isotropic_matrices.back());
1072 }
1073 else
1074 {
1075 // must compute both restriction and prolongation matrices because
1076 // we only check for their size and the reinit call initializes them
1077 // all
1080 this_nonconst.prolongation);
1082 this_nonconst.restriction);
1083 }
1084 }
1085
1086 // finally return the matrix
1087 return this->prolongation[refinement_case - 1][child];
1088}
1089
1090
1091
1092template <int dim>
1093const FullMatrix<double> &
1095 const unsigned int child,
1096 const RefinementCase<dim> &refinement_case) const
1097{
1098 AssertIndexRange(refinement_case,
1100 Assert(refinement_case != RefinementCase<dim>::no_refinement,
1101 ExcMessage(
1102 "Restriction matrices are only available for refined cells!"));
1103 AssertIndexRange(child, GeometryInfo<dim>::n_children(refinement_case));
1104
1105 // initialization upon first request
1106 if (this->restriction[refinement_case - 1][child].n() == 0)
1107 {
1108 std::lock_guard<std::mutex> lock(this->mutex);
1109
1110 // if matrix got updated while waiting for the lock...
1111 if (this->restriction[refinement_case - 1][child].n() ==
1112 this->n_dofs_per_cell())
1113 return this->restriction[refinement_case - 1][child];
1114
1115 // now do the work. need to get a non-const version of data in order to
1116 // be able to modify them inside a const function
1117 FE_RaviartThomasNodal<dim> &this_nonconst =
1118 const_cast<FE_RaviartThomasNodal<dim> &>(*this);
1119 if (refinement_case == RefinementCase<dim>::isotropic_refinement)
1120 {
1121 std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
1123 isotropic_matrices.back().resize(
1125 FullMatrix<double>(this->n_dofs_per_cell(),
1126 this->n_dofs_per_cell()));
1127 FETools::compute_projection_matrices(*this, isotropic_matrices, true);
1128 this_nonconst.restriction[refinement_case - 1].swap(
1129 isotropic_matrices.back());
1130 }
1131 else
1132 {
1133 // must compute both restriction and prolongation matrices because
1134 // we only check for their size and the reinit call initializes them
1135 // all
1138 this_nonconst.prolongation);
1140 this_nonconst.restriction);
1141 }
1142 }
1143
1144 // finally return the matrix
1145 return this->restriction[refinement_case - 1][child];
1146}
1147
1148
1149
1150// explicit instantiations
1151#include "fe_raviart_thomas_nodal.inst"
1152
1153
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: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: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
static void project_to_subface(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 SubQuadrature &quadrature, const unsigned int face_no, std::vector< Point< dim > > &q_points)
const Point< dim > & point(const unsigned int i) const
SymmetricTensor< rank, dim, Number > sum(const SymmetricTensor< rank, dim, Number > &local, const MPI_Comm &mpi_communicator)
virtual std::string name() const =0
virtual void evaluate(const Point< dim > &unit_point, std::vector< Tensor< 1, dim > > &values, std::vector< Tensor< 2, dim > > &grads, std::vector< Tensor< 3, dim > > &grad_grads, std::vector< Tensor< 4, dim > > &third_derivatives, std::vector< Tensor< 5, dim > > &fourth_derivatives) const =0
virtual std::unique_ptr< TensorPolynomialsBase< dim > > clone() const =0
unsigned int degree() const
Definition: vector.h:109
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
Definition: exceptions.h:1473
std::string to_string(const T &t)
Definition: patterns.h:2403
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
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)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
std::vector< Polynomial< double > > generate_complete_Lagrange_basis(const std::vector< Point< 1 > > &points)
Definition: polynomial.cc:702
constexpr T pow(const T base, const int iexp)
Definition: utilities.h:462
std::vector< Integer > invert_permutation(const std::vector< Integer > &permutation)
Definition: utilities.h:1813
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618
std::vector< unsigned int > get_lexicographic_numbering_rt_nodal(const unsigned int degree)
STL namespace.