deal.II version GIT relicensing-2659-g040196caa3 2025-02-18 14:20: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_q_base.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2013 - 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
25
26#include <deal.II/fe/fe_dgp.h>
27#include <deal.II/fe/fe_dgq.h>
34#include <deal.II/fe/fe_tools.h>
36
37#include <memory>
38#include <sstream>
39#include <vector>
40
42
43
44namespace internal
45{
46 namespace FE_Q_Base
47 {
48 namespace
49 {
50 // in get_restriction_matrix() and get_prolongation_matrix(), want to undo
51 // tensorization on inner loops for performance reasons. this clears a
52 // dim-array
53 template <int dim>
54 inline void
55 zero_indices(unsigned int (&indices)[dim])
56 {
57 for (unsigned int d = 0; d < dim; ++d)
58 indices[d] = 0;
59 }
60
61
62
63 // in get_restriction_matrix() and get_prolongation_matrix(), want to undo
64 // tensorization on inner loops for performance reasons. this increments
65 // tensor product indices
66 template <int dim>
67 inline void
68 increment_indices(unsigned int (&indices)[dim], const unsigned int dofs1d)
69 {
70 ++indices[0];
71 for (unsigned int d = 0; d < dim - 1; ++d)
72 if (indices[d] == dofs1d)
73 {
74 indices[d] = 0;
75 indices[d + 1]++;
76 }
77 }
78 } // namespace
79 } // namespace FE_Q_Base
80} // namespace internal
81
82
83
88template <int xdim, int xspacedim>
90{
95 template <int spacedim>
96 static void
97 initialize_constraints(const std::vector<Point<1>> &,
99 {
100 // no constraints in 1d
101 }
102
103
104 template <int spacedim>
105 static void
106 initialize_constraints(const std::vector<Point<1>> & /*points*/,
108 {
109 const unsigned int dim = 2;
110
111 unsigned int q_deg = fe.degree;
112 if (dynamic_cast<const TensorProductPolynomialsBubbles<dim> *>(
113 &fe.get_poly_space()) != nullptr)
114 q_deg = fe.degree - 1;
115
116 // restricted to each face, the traces of the shape functions is an
117 // element of P_{k} (in 2d), or Q_{k} (in 3d), where k is the degree of
118 // the element. from this, we interpolate between mother and cell face.
119
120 // the interpolation process works as follows: on each subface, we want
121 // that finite element solutions from both sides coincide. i.e. if a and b
122 // are expansion coefficients for the shape functions from both sides, we
123 // seek a relation between a and b such that
124 // sum_j a_j phi^c_j(x) == sum_j b_j phi_j(x)
125 // for all points x on the interface. here, phi^c_j are the shape
126 // functions on the small cell on one side of the face, and phi_j those on
127 // the big cell on the other side. To get this relation, it suffices to
128 // look at a sufficient number of points for which this has to hold. if
129 // there are n functions, then we need n evaluation points, and we choose
130 // them equidistantly.
131 //
132 // we obtain the matrix system
133 // A a == B b
134 // where
135 // A_ij = phi^c_j(x_i)
136 // B_ij = phi_j(x_i)
137 // and the relation we are looking for is
138 // a = A^-1 B b
139 //
140 // for the special case of Lagrange interpolation polynomials, A_ij
141 // reduces to delta_ij, and
142 // a_i = B_ij b_j
143 // Hence, interface_constraints(i,j)=B_ij.
144 //
145 // for the general case, where we don't have Lagrange interpolation
146 // polynomials, this is a little more complicated. Then we would evaluate
147 // at a number of points and invert the interpolation matrix A.
148 //
149 // Note, that we build up these matrices for all subfaces at once, rather
150 // than considering them separately. the reason is that we finally will
151 // want to have them in this order anyway, as this is the format we need
152 // inside deal.II
153
154 // In the following the points x_i are constructed in following order
155 // (n=degree-1)
156 // *----------*---------*
157 // 1..n 0 n+1..2n
158 // i.e. first the midpoint of the line, then the support points on subface
159 // 0 and on subface 1
160 std::vector<Point<dim - 1>> constraint_points;
161 // Add midpoint
162 constraint_points.emplace_back(0.5);
163
164 // subface 0
165 for (unsigned int i = 1; i < q_deg; ++i)
166 constraint_points.push_back(
168 Point<dim - 1>(i / double(q_deg)), 0));
169 // subface 1
170 for (unsigned int i = 1; i < q_deg; ++i)
171 constraint_points.push_back(
173 Point<dim - 1>(i / double(q_deg)), 1));
174
175 // Now construct relation between destination (child) and source (mother)
176 // dofs.
177
178 fe.interface_constraints.TableBase<2, double>::reinit(
179 fe.interface_constraints_size());
180
181 // use that the element evaluates to 1 at index 0 and along the line at
182 // zero
183 const std::vector<unsigned int> &index_map_inverse =
184 fe.get_poly_space_numbering_inverse();
185 const std::vector<unsigned int> face_index_map =
188 fe.poly_space->compute_value(index_map_inverse[0], Point<dim>()) -
189 1.) < 1e-14,
191
192 for (unsigned int i = 0; i < constraint_points.size(); ++i)
193 for (unsigned int j = 0; j < q_deg + 1; ++j)
194 {
195 Point<dim> p;
196 p[0] = constraint_points[i][0];
197 fe.interface_constraints(i, face_index_map[j]) =
198 fe.poly_space->compute_value(index_map_inverse[j], p);
199
200 // if the value is small up to round-off, then simply set it to zero
201 // to avoid unwanted fill-in of the constraint matrices (which would
202 // then increase the number of other DoFs a constrained DoF would
203 // couple to)
204 if (std::fabs(fe.interface_constraints(i, face_index_map[j])) < 1e-13)
205 fe.interface_constraints(i, face_index_map[j]) = 0;
206 }
207 }
208
209
210 template <int spacedim>
211 static void
212 initialize_constraints(const std::vector<Point<1>> & /*points*/,
214 {
215 const unsigned int dim = 3;
216
217 unsigned int q_deg = fe.degree;
218 if (dynamic_cast<const TensorProductPolynomialsBubbles<dim> *>(
219 &fe.get_poly_space()) != nullptr)
220 q_deg = fe.degree - 1;
221
222 // For a detailed documentation of the interpolation see the
223 // FE_Q_Base<2>::initialize_constraints function.
224
225 // In the following the points x_i are constructed in the order as
226 // described in the documentation of the FiniteElement class (fe_base.h),
227 // i.e.
228 // *--15--4--16--*
229 // | | |
230 // 10 19 6 20 12
231 // | | |
232 // 1--7---0--8---2
233 // | | |
234 // 9 17 5 18 11
235 // | | |
236 // *--13--3--14--*
237 std::vector<Point<dim - 1>> constraint_points;
238
239 // Add midpoint
240 constraint_points.emplace_back(0.5, 0.5);
241
242 // Add midpoints of lines of "mother-face"
243 constraint_points.emplace_back(0, 0.5);
244 constraint_points.emplace_back(1, 0.5);
245 constraint_points.emplace_back(0.5, 0);
246 constraint_points.emplace_back(0.5, 1);
247
248 if (q_deg > 1)
249 {
250 const unsigned int n = q_deg - 1;
251 const double step = 1. / q_deg;
252 std::vector<Point<1>> line_support_points(n);
253 for (unsigned int i = 0; i < n; ++i)
254 line_support_points[i][0] = (i + 1) * step;
255 const Quadrature<1> qline(line_support_points);
256
257 // Add nodes of lines interior in the "mother-face"
258 auto get_points = [&](const unsigned int face_no,
259 const unsigned int subface_no) {
261 ReferenceCells::get_hypercube<2>(),
262 qline,
263 face_no,
267 .get_points();
268 };
269
270 // line 5: use line 9
271 for (const Point<dim - 1> &p : get_points(0, 0))
272 constraint_points.push_back(p + Point<dim - 1>(0.5, 0));
273 // line 6: use line 10
274 for (const Point<dim - 1> &p : get_points(0, 1))
275 constraint_points.push_back(p + Point<dim - 1>(0.5, 0));
276 // line 7: use line 13
277 for (const Point<dim - 1> &p : get_points(2, 0))
278 constraint_points.push_back(p + Point<dim - 1>(0, 0.5));
279 // line 8: use line 14
280 for (const Point<dim - 1> &p : get_points(2, 1))
281 constraint_points.push_back(p + Point<dim - 1>(0, 0.5));
282
283 // DoFs on bordering lines lines 9-16
284 for (unsigned int face = 0;
285 face < GeometryInfo<dim - 1>::faces_per_cell;
286 ++face)
287 for (unsigned int subface = 0;
288 subface < GeometryInfo<dim - 1>::max_children_per_face;
289 ++subface)
290 {
291 const auto p_line = get_points(face, subface);
293 p_line.begin(),
294 p_line.end());
295 }
296
297 // Create constraints for interior nodes
298 std::vector<Point<dim - 1>> inner_points(n * n);
299 for (unsigned int i = 0, iy = 1; iy <= n; ++iy)
300 for (unsigned int ix = 1; ix <= n; ++ix)
301 inner_points[i++] = Point<dim - 1>(ix * step, iy * step);
302
303 // at the moment do this for isotropic face refinement only
304 for (unsigned int child = 0;
305 child < GeometryInfo<dim - 1>::max_children_per_cell;
306 ++child)
307 for (const auto &inner_point : inner_points)
308 constraint_points.push_back(
310 child));
311 }
312
313 // Now construct relation between destination (child) and source (mother)
314 // dofs.
315 const unsigned int pnts = (q_deg + 1) * (q_deg + 1);
316
317 // use that the element evaluates to 1 at index 0 and along the line at
318 // zero
319 const std::vector<unsigned int> &index_map_inverse =
320 fe.get_poly_space_numbering_inverse();
321 const std::vector<unsigned int> face_index_map =
324 fe.poly_space->compute_value(index_map_inverse[0], Point<dim>()) -
325 1.) < 1e-14,
327
328 fe.interface_constraints.TableBase<2, double>::reinit(
329 fe.interface_constraints_size());
330
331 for (unsigned int i = 0; i < constraint_points.size(); ++i)
332 {
333 const double interval = static_cast<double>(q_deg * 2);
334 bool mirror[dim - 1];
336
337 // Eliminate FP errors in constraint points. Due to their origin, they
338 // must all be fractions of the unit interval. If we have polynomial
339 // degree 4, the refined element has 8 intervals. Hence the
340 // coordinates must be 0, 0.125, 0.25, 0.375 etc. Now the coordinates
341 // of the constraint points will be multiplied by the inverse of the
342 // interval size (in the example by 8). After that the coordinates
343 // must be integral numbers. Hence a normal truncation is performed
344 // and the coordinates will be scaled back. The equal treatment of all
345 // coordinates should eliminate any FP errors.
346 for (unsigned int k = 0; k < dim - 1; ++k)
347 {
348 const int coord_int =
349 static_cast<int>(constraint_points[i][k] * interval + 0.25);
350 constraint_point[k] = 1. * coord_int / interval;
351
352 // The following lines of code should eliminate the problems with
353 // the constraints object which appeared for P>=4. The
354 // AffineConstraints class complained about different constraints
355 // for the same entry: Actually, this
356 // difference could be attributed to FP errors, as it was in the
357 // range of 1.0e-16. These errors originate in the loss of
358 // symmetry in the FP approximation of the shape-functions.
359 // Considering a 3rd order shape function in 1d, we have
360 // N0(x)=N3(1-x) and N1(x)=N2(1-x). For higher order polynomials
361 // the FP approximations of the shape functions do not satisfy
362 // these equations any more! Thus in the following code
363 // everything is computed in the interval x \in [0..0.5], which is
364 // sufficient to express all values that could come out from a
365 // computation of any shape function in the full interval
366 // [0..1]. If x > 0.5 the computation is done for 1-x with the
367 // shape function N_{p-n} instead of N_n. Hence symmetry is
368 // preserved and everything works fine...
369 //
370 // For a different explanation of the problem, see the discussion
371 // in the FiniteElement class for constraint matrices in 3d.
372 mirror[k] = (constraint_point[k] > 0.5);
373 if (mirror[k])
375 }
376
377 for (unsigned int j = 0; j < pnts; ++j)
378 {
379 unsigned int indices[2] = {j % (q_deg + 1), j / (q_deg + 1)};
380
381 for (unsigned int k = 0; k < 2; ++k)
382 if (mirror[k])
383 indices[k] = q_deg - indices[k];
384
385 const unsigned int new_index =
386 indices[1] * (q_deg + 1) + indices[0];
387
388 fe.interface_constraints(i, face_index_map[j]) =
389 fe.poly_space->compute_value(index_map_inverse[new_index],
391
392 // if the value is small up to round-off, then simply set it to
393 // zero to avoid unwanted fill-in of the constraint matrices
394 // (which would then increase the number of other DoFs a
395 // constrained DoF would couple to)
396 if (std::fabs(fe.interface_constraints(i, face_index_map[j])) <
397 1e-13)
398 fe.interface_constraints(i, face_index_map[j]) = 0;
399 }
400 }
401 }
402};
403
404#ifndef DOXYGEN
405
406template <int dim, int spacedim>
408 const ScalarPolynomialsBase<dim> &poly_space,
409 const FiniteElementData<dim> &fe_data,
410 const std::vector<bool> &restriction_is_additive_flags)
411 : FE_Poly<dim, spacedim>(
412 poly_space,
413 fe_data,
414 restriction_is_additive_flags,
415 std::vector<ComponentMask>(1, ComponentMask(std::vector<bool>(1, true))))
417 &poly_space) != nullptr ?
418 this->degree - 1 :
419 this->degree)
420{}
421
422
423
424template <int dim, int spacedim>
425void
426FE_Q_Base<dim, spacedim>::initialize(const std::vector<Point<1>> &points)
427{
428 Assert(points[0][0] == 0,
429 ExcMessage("The first support point has to be zero."));
430 Assert(points.back()[0] == 1,
431 ExcMessage("The last support point has to be one."));
432
433 // distinguish q/q_dg0 case: need to be flexible enough to allow more
434 // degrees of freedom than there are FE_Q degrees of freedom for derived
435 // class FE_Q_DG0 that otherwise shares 95% of the code.
436 const unsigned int q_dofs_per_cell =
437 Utilities::fixed_power<dim>(q_degree + 1);
439 q_dofs_per_cell + 1 == this->n_dofs_per_cell() ||
440 q_dofs_per_cell + dim == this->n_dofs_per_cell(),
442
443 [this, q_dofs_per_cell]() {
444 std::vector<unsigned int> renumber =
445 FETools::hierarchic_to_lexicographic_numbering<dim>(q_degree);
446 for (unsigned int i = q_dofs_per_cell; i < this->n_dofs_per_cell(); ++i)
447 renumber.push_back(i);
449 dynamic_cast<TensorProductPolynomials<dim> *>(this->poly_space.get());
450 if (tensor_poly_space_ptr != nullptr)
451 {
452 tensor_poly_space_ptr->set_numbering(renumber);
453 return;
454 }
455 auto *tensor_piecewise_poly_space_ptr = dynamic_cast<
457 *>(this->poly_space.get());
458 if (tensor_piecewise_poly_space_ptr != nullptr)
459 {
460 tensor_piecewise_poly_space_ptr->set_numbering(renumber);
461 return;
462 }
465 this->poly_space.get());
466 if (tensor_bubbles_poly_space_ptr != nullptr)
467 {
468 tensor_bubbles_poly_space_ptr->set_numbering(renumber);
469 return;
470 }
473 this->poly_space.get());
474 if (tensor_const_poly_space_ptr != nullptr)
475 {
476 tensor_const_poly_space_ptr->set_numbering(renumber);
477 return;
478 }
480 }();
481
482 // Finally fill in support points on cell and face and initialize
483 // constraints. All of this can happen in parallel
485 tasks += Threads::new_task([&]() { initialize_unit_support_points(points); });
486 tasks +=
488 tasks += Threads::new_task([&]() { initialize_constraints(points); });
489 tasks +=
491 tasks.join_all();
492
493 // do not initialize embedding and restriction here. these matrices are
494 // initialized on demand in get_restriction_matrix and
495 // get_prolongation_matrix
496}
497
498
499
500template <int dim, int spacedim>
501void
505{
506 // go through the list of elements we can interpolate from
508 dynamic_cast<const FE_Q_Base<dim, spacedim> *>(&x_source_fe))
509 {
510 // ok, source is a Q element, so we will be able to do the work
511 Assert(interpolation_matrix.m() == this->n_dofs_per_cell(),
513 this->n_dofs_per_cell()));
514 Assert(interpolation_matrix.n() == x_source_fe.n_dofs_per_cell(),
516 x_source_fe.n_dofs_per_cell()));
517
518 // only evaluate Q dofs
519 const unsigned int q_dofs_per_cell =
520 Utilities::fixed_power<dim>(q_degree + 1);
521 const unsigned int source_q_dofs_per_cell =
522 Utilities::fixed_power<dim>(source_fe->degree + 1);
523
524 // evaluation is simply done by evaluating the other FE's basis functions
525 // on the unit support points (FE_Q has the property that the cell
526 // interpolation matrix is a unit matrix, so no need to evaluate it and
527 // invert it)
528 for (unsigned int j = 0; j < q_dofs_per_cell; ++j)
529 {
530 // read in a point on this cell and evaluate the shape functions there
531 const Point<dim> p = this->unit_support_points[j];
532
533 // FE_Q element evaluates to 1 in unit support point and to zero in
534 // all other points by construction
535 Assert(std::abs(this->poly_space->compute_value(j, p) - 1.) < 1e-13,
537
538 for (unsigned int i = 0; i < source_q_dofs_per_cell; ++i)
540 source_fe->poly_space->compute_value(i, p);
541 }
542
543 // for FE_Q_DG0, add one last row of identity
545 {
547 source_fe->n_dofs_per_cell());
548 for (unsigned int i = 0; i < source_q_dofs_per_cell; ++i)
550 for (unsigned int j = 0; j < q_dofs_per_cell; ++j)
553 }
554
555 // cut off very small values
556 const double eps = 2e-13 * q_degree * dim;
557 for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
558 for (unsigned int j = 0; j < source_fe->n_dofs_per_cell(); ++j)
559 if (std::fabs(interpolation_matrix(i, j)) < eps)
560 interpolation_matrix(i, j) = 0.;
561
562# ifdef DEBUG
563 // make sure that the row sum of each of the matrices is 1 at this
564 // point. this must be so since the shape functions sum up to 1
565 for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
566 {
567 double sum = 0.;
568 for (unsigned int j = 0; j < source_fe->n_dofs_per_cell(); ++j)
569 sum += interpolation_matrix(i, j);
570
571 Assert(std::fabs(sum - 1) < eps, ExcInternalError());
572 }
573# endif
574 }
575 else if (dynamic_cast<const FE_Nothing<dim> *>(&x_source_fe))
576 {
577 // the element we want to interpolate from is an FE_Nothing. this
578 // element represents a function that is constant zero and has no
579 // degrees of freedom, so the interpolation is simply a multiplication
580 // with a n_dofs x 0 matrix. there is nothing to do here
581
582 // we would like to verify that the number of rows and columns of
583 // the matrix equals this->n_dofs_per_cell() and zero. unfortunately,
584 // whenever we do FullMatrix::reinit(m,0), it sets both rows and
585 // columns to zero, instead of m and zero. thus, only test the
586 // number of columns
587 Assert(interpolation_matrix.n() == x_source_fe.n_dofs_per_cell(),
589 x_source_fe.n_dofs_per_cell()));
590 }
591 else
593 false,
594 (typename FiniteElement<dim,
596}
597
598
599
600template <int dim, int spacedim>
601void
605 const unsigned int face_no) const
606{
610 face_no);
611}
612
613
614
615template <int dim, int spacedim>
616void
619 const unsigned int subface,
621 const unsigned int face_no) const
622{
623 Assert(interpolation_matrix.m() == source_fe.n_dofs_per_face(face_no),
625 source_fe.n_dofs_per_face(face_no)));
626
627 Assert(source_fe.n_components() == this->n_components(),
628 ExcDimensionMismatch(source_fe.n_components(), this->n_components()));
629
630 if ((source_fe.has_face_support_points(face_no)) &&
631 (source_fe.n_dofs_per_face(face_no) > 0))
632 {
633 // have this test in here since a table of size 2x0 reports its size as
634 // 0x0
635 Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
637 this->n_dofs_per_face(face_no)));
638
639 // Make sure that the element for which the DoFs should be constrained
640 // is the one with the higher polynomial degree. Actually the procedure
641 // will work also if this assertion is not satisfied. But the matrices
642 // produced in that case might lead to problems in the hp-procedures,
643 // which use this method.
644 Assert(
645 this->n_dofs_per_face(face_no) <= source_fe.n_dofs_per_face(face_no),
646 (typename FiniteElement<dim,
648
649 // generate a point on this cell and evaluate the shape functions there
650 const Quadrature<dim - 1> quad_face_support(
651 source_fe.get_unit_face_support_points(face_no));
652
653 // Rule of thumb for FP accuracy, that can be expected for a given
654 // polynomial degree. This value is used to cut off values close to
655 // zero.
656 const double eps = 2e-13 * this->q_degree * std::max(dim - 1, 1);
657
658 // compute the interpolation matrix by simply taking the value at the
659 // support points.
660 // TODO: Verify that all faces are the same with respect to
661 // these support points. Furthermore, check if something has to
662 // be done for the face orientation flag in 3d.
666 this->reference_cell(),
668 0,
670 QProjector<dim>::project_to_subface(
673 0,
674 subface,
677 for (unsigned int i = 0; i < source_fe.n_dofs_per_face(face_no); ++i)
678 {
679 const Point<dim> &p = subface_quadrature.point(i);
680
681 for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
682 {
683 double matrix_entry =
684 this->shape_value(this->face_to_cell_index(j, 0), p);
685
686 // Correct the interpolated value. I.e. if it is close to 1 or
687 // 0, make it exactly 1 or 0. Unfortunately, this is required to
688 // avoid problems with higher order elements.
689 if (std::fabs(matrix_entry - 1.0) < eps)
690 matrix_entry = 1.0;
691 if (std::fabs(matrix_entry) < eps)
692 matrix_entry = 0.0;
693
695 }
696 }
697
698# ifdef DEBUG
699 // make sure that the row sum of each of the matrices is 1 at this
700 // point. this must be so since the shape functions sum up to 1
701 for (unsigned int j = 0; j < source_fe.n_dofs_per_face(face_no); ++j)
702 {
703 double sum = 0.;
704
705 for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
706 sum += interpolation_matrix(j, i);
707
708 Assert(std::fabs(sum - 1) < eps, ExcInternalError());
709 }
710# endif
711 }
712 else if (dynamic_cast<const FE_Nothing<dim> *>(&source_fe) != nullptr)
713 {
714 // nothing to do here, the FE_Nothing has no degrees of freedom anyway
715 }
716 else
718 false,
719 (typename FiniteElement<dim,
721}
722
723
724
725template <int dim, int spacedim>
726bool
728{
729 return true;
730}
731
732
733
734template <int dim, int spacedim>
735std::vector<std::pair<unsigned int, unsigned int>>
738{
739 if (dynamic_cast<const FE_Q_Base<dim, spacedim> *>(&fe_other) != nullptr)
740 {
741 // there should be exactly one single DoF of each FE at a vertex, and they
742 // should have identical value
743 return {{0U, 0U}};
744 }
745 else if (dynamic_cast<const FE_SimplexP<dim, spacedim> *>(&fe_other) !=
746 nullptr)
747 {
748 // there should be exactly one single DoF of each FE at a vertex, and they
749 // should have identical value
750 return {{0U, 0U}};
751 }
752 else if (dynamic_cast<const FE_Hermite<dim, spacedim> *>(&fe_other) !=
753 nullptr)
754 {
755 // FE_Hermite will usually have several degrees of freedom on
756 // each vertex, however only the first one will actually
757 // correspond to the shape value at the vertex, meaning it's
758 // the only one of interest for FE_Q_Base
759 return {{0U, 0U}};
760 }
761 else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
762 {
763 // the FE_Nothing has no degrees of freedom, so there are no
764 // equivalencies to be recorded
765 return {};
766 }
767 else if (fe_other.n_unique_faces() == 1 && fe_other.n_dofs_per_face(0) == 0)
768 {
769 // if the other element has no elements on faces at all,
770 // then it would be impossible to enforce any kind of
771 // continuity even if we knew exactly what kind of element
772 // we have -- simply because the other element declares
773 // that it is discontinuous because it has no DoFs on
774 // its faces. in that case, just state that we have no
775 // constraints to declare
776 return {};
777 }
778 else
779 {
781 return {};
782 }
783}
784
785
786
787template <int dim, int spacedim>
788std::vector<std::pair<unsigned int, unsigned int>>
791{
792 // we can presently only compute these identities if both FEs are FE_Qs or
793 // if the other one is an FE_Nothing
795 dynamic_cast<const FE_Q_Base<dim, spacedim> *>(&fe_other))
796 {
797 // dofs are located along lines, so two dofs are identical if they are
798 // located at identical positions. if we had only equidistant points, we
799 // could simply check for similarity like (i+1)*q == (j+1)*p, but we
800 // might have other support points (e.g. Gauss-Lobatto
801 // points). Therefore, read the points in unit_support_points for the
802 // first coordinate direction. We take the lexicographic ordering of the
803 // points in the first direction (i.e., x-direction), which we access
804 // between index 1 and p-1 (index 0 and p are vertex dofs).
805 const unsigned int p = this->degree;
806 const unsigned int q = fe_q_other->degree;
807
808 std::vector<std::pair<unsigned int, unsigned int>> identities;
809
810 const std::vector<unsigned int> &index_map_inverse =
812 const std::vector<unsigned int> &index_map_inverse_other =
813 fe_q_other->get_poly_space_numbering_inverse();
814
815 for (unsigned int i = 0; i < p - 1; ++i)
816 for (unsigned int j = 0; j < q - 1; ++j)
817 if (std::fabs(
818 this->unit_support_points[index_map_inverse[i + 1]][0] -
819 fe_q_other->unit_support_points[index_map_inverse_other[j + 1]]
820 [0]) < 1e-14)
821 identities.emplace_back(i, j);
822
823 return identities;
824 }
825 else if (const FE_SimplexP<dim, spacedim> *fe_p_other =
826 dynamic_cast<const FE_SimplexP<dim, spacedim> *>(&fe_other))
827 {
828 // DoFs are located along lines, so two dofs are identical if they are
829 // located at identical positions. If we had only equidistant points, we
830 // could simply check for similarity like (i+1)*q == (j+1)*p, but we
831 // might have other support points (e.g. Gauss-Lobatto
832 // points). Therefore, read the points in unit_support_points for the
833 // first coordinate direction. For FE_Q, we take the lexicographic
834 // ordering of the line support points in the first direction (i.e.,
835 // x-direction), which we access between index 1 and p-1 (index 0 and p
836 // are vertex dofs). For FE_SimplexP, they are currently hard-coded and we
837 // iterate over points on the first line which begin after the 3 vertex
838 // points in the complete list of unit support points
839
840 Assert(fe_p_other->degree <= 2, ExcNotImplemented());
841
842 const std::vector<unsigned int> &index_map_inverse_q =
844
845 std::vector<std::pair<unsigned int, unsigned int>> identities;
846
847 for (unsigned int i = 0; i < this->degree - 1; ++i)
848 for (unsigned int j = 0; j < fe_p_other->degree - 1; ++j)
849 if (std::fabs(
850 this->unit_support_points[index_map_inverse_q[i + 1]][0] -
851 fe_p_other->get_unit_support_points()[j + 3][0]) < 1e-14)
852 identities.emplace_back(i, j);
853
854 return identities;
855 }
856 else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
857 {
858 // the FE_Nothing has no degrees of freedom, so there are no
859 // equivalencies to be recorded
860 return {};
861 }
862 else if (fe_other.n_unique_faces() == 1 && fe_other.n_dofs_per_face(0) == 0)
863 {
864 // if the other element has no elements on faces at all,
865 // then it would be impossible to enforce any kind of
866 // continuity even if we knew exactly what kind of element
867 // we have -- simply because the other element declares
868 // that it is discontinuous because it has no DoFs on
869 // its faces. in that case, just state that we have no
870 // constraints to declare
871 return {};
872 }
873 else
874 {
876 return {};
877 }
878}
879
880
881
882template <int dim, int spacedim>
883std::vector<std::pair<unsigned int, unsigned int>>
886 const unsigned int) const
887{
888 // we can presently only compute these identities if both FEs are FE_Qs or
889 // if the other one is an FE_Nothing
891 dynamic_cast<const FE_Q_Base<dim, spacedim> *>(&fe_other))
892 {
893 // this works exactly like the line case above, except that now we have
894 // to have two indices i1, i2 and j1, j2 to characterize the dofs on the
895 // face of each of the finite elements. since they are ordered
896 // lexicographically along the first line and we have a tensor product,
897 // the rest is rather straightforward
898 const unsigned int p = this->degree;
899 const unsigned int q = fe_q_other->degree;
900
901 std::vector<std::pair<unsigned int, unsigned int>> identities;
902
903 const std::vector<unsigned int> &index_map_inverse =
905 const std::vector<unsigned int> &index_map_inverse_other =
906 fe_q_other->get_poly_space_numbering_inverse();
907
908 for (unsigned int i1 = 0; i1 < p - 1; ++i1)
909 for (unsigned int i2 = 0; i2 < p - 1; ++i2)
910 for (unsigned int j1 = 0; j1 < q - 1; ++j1)
911 for (unsigned int j2 = 0; j2 < q - 1; ++j2)
912 if ((std::fabs(
913 this->unit_support_points[index_map_inverse[i1 + 1]][0] -
916 [0]) < 1e-14) &&
917 (std::fabs(
918 this->unit_support_points[index_map_inverse[i2 + 1]][0] -
921 [0]) < 1e-14))
922 identities.emplace_back(i1 * (p - 1) + i2, j1 * (q - 1) + j2);
923
924 return identities;
925 }
926 else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
927 {
928 // the FE_Nothing has no degrees of freedom, so there are no
929 // equivalencies to be recorded
930 return std::vector<std::pair<unsigned int, unsigned int>>();
931 }
932 else if (fe_other.n_unique_faces() == 1 && fe_other.n_dofs_per_face(0) == 0)
933 {
934 // if the other element has no elements on faces at all,
935 // then it would be impossible to enforce any kind of
936 // continuity even if we knew exactly what kind of element
937 // we have -- simply because the other element declares
938 // that it is discontinuous because it has no DoFs on
939 // its faces. in that case, just state that we have no
940 // constraints to declare
941 return std::vector<std::pair<unsigned int, unsigned int>>();
942 }
943 else
944 {
946 return std::vector<std::pair<unsigned int, unsigned int>>();
947 }
948}
949
950
951
952//---------------------------------------------------------------------------
953// Auxiliary functions
954//---------------------------------------------------------------------------
955
956
957
958template <int dim, int spacedim>
959void
961 const std::vector<Point<1>> &points)
962{
963 const std::vector<unsigned int> &index_map_inverse =
965
966 // We can compute the support points by computing the tensor
967 // product of the 1d set of points. We could do this by hand, but it's
968 // easier to just re-use functionality that's already been implemented
969 // for quadrature formulas.
970 const Quadrature<1> support_1d(points);
971 const Quadrature<dim> support_quadrature(support_1d); // NOLINT
972
973 // The only thing we have to do is reorder the points from tensor
974 // product order to the order in which we enumerate DoFs on cells
975 this->unit_support_points.resize(support_quadrature.size());
976 for (unsigned int k = 0; k < support_quadrature.size(); ++k)
977 this->unit_support_points[index_map_inverse[k]] =
978 support_quadrature.point(k);
979}
980
981
982
983template <int dim, int spacedim>
984void
986 const std::vector<Point<1>> &points)
987{
988 // TODO: the implementation makes the assumption that all faces have the
989 // same number of dofs
991 const unsigned int face_no = 0;
992
993 this->unit_face_support_points[face_no].resize(
994 Utilities::fixed_power<dim - 1>(q_degree + 1));
995
996 // In 1d, there is only one 0-dimensional support point, so there is nothing
997 // more to be done.
998 if (dim == 1)
999 return;
1000
1001 // find renumbering of faces and assign from values of quadrature
1002 const std::vector<unsigned int> face_index_map =
1004
1005 // We can compute the support points by computing the tensor
1006 // product of the 1d set of points. We could do this by hand, but it's
1007 // easier to just re-use functionality that's already been implemented
1008 // for quadrature formulas.
1009 const Quadrature<1> support_1d(points);
1010 const Quadrature<dim - 1> support_quadrature(support_1d); // NOLINT
1011
1012 // The only thing we have to do is reorder the points from tensor
1013 // product order to the order in which we enumerate DoFs on cells
1014 this->unit_face_support_points[face_no].resize(support_quadrature.size());
1015 for (unsigned int k = 0; k < support_quadrature.size(); ++k)
1016 this->unit_face_support_points[face_no][face_index_map[k]] =
1017 support_quadrature.point(k);
1018}
1019
1020
1021
1022template <int dim, int spacedim>
1023void
1025{
1026 // initialize reordering of line dofs
1027 for (unsigned int i = 0; i < this->n_dofs_per_line(); ++i)
1029 this->n_dofs_per_line() - 1 - i - i;
1030
1031 // for 1d and 2d we can skip adjust_quad_dof_index_for_face_orientation_table
1032 if (dim < 3)
1033 return;
1034
1035 // TODO: the implementation makes the assumption that all faces have the
1036 // same number of dofs
1037 AssertDimension(this->n_unique_faces(), 1);
1038 const unsigned int face_no = 0;
1039
1040 Assert(
1042 this->reference_cell().n_face_orientations(face_no) *
1043 this->n_dofs_per_quad(face_no),
1045
1046 const unsigned int n = q_degree - 1;
1047 Assert(n * n == this->n_dofs_per_quad(face_no), ExcInternalError());
1048
1049 // the dofs on a face are connected to a n x n matrix. for example, for
1050 // degree==4 we have the following dofs on a quad
1051
1052 // ___________
1053 // | |
1054 // | 6 7 8 |
1055 // | |
1056 // | 3 4 5 |
1057 // | |
1058 // | 0 1 2 |
1059 // |___________|
1060 //
1061 // we have dof_no=i+n*j with index i in x-direction and index j in
1062 // y-direction running from 0 to n-1. to extract i and j we can use
1063 // i=dof_no%n and j=dof_no/n. i and j can then be used to construct the
1064 // rotated and mirrored numbers.
1065
1066
1067 for (unsigned int local = 0; local < this->n_dofs_per_quad(face_no); ++local)
1068 // face support points are in lexicographic ordering with x running
1069 // fastest. invert that (y running fastest)
1070 {
1071 unsigned int i = local % n, j = local / n;
1072
1073 // face_orientation=false, face_flip=false, face_rotation=false
1075 local, internal::combined_face_orientation(false, false, false)) =
1076 j + i * n - local;
1077 // face_orientation=false, face_flip=false, face_rotation=true
1079 local, internal::combined_face_orientation(false, true, false)) =
1080 i + (n - 1 - j) * n - local;
1081 // face_orientation=false, face_flip=true, face_rotation=false
1083 local, internal::combined_face_orientation(false, false, true)) =
1084 (n - 1 - j) + (n - 1 - i) * n - local;
1085 // face_orientation=false, face_flip=true, face_rotation=true
1087 local, internal::combined_face_orientation(false, true, true)) =
1088 (n - 1 - i) + j * n - local;
1089 // face_orientation=true, face_flip=false, face_rotation=false
1091 local, internal::combined_face_orientation(true, false, false)) = 0;
1092 // face_orientation=true, face_flip=false, face_rotation=true
1094 local, internal::combined_face_orientation(true, true, false)) =
1095 j + (n - 1 - i) * n - local;
1096 // face_orientation=true, face_flip=true, face_rotation=false
1098 local, internal::combined_face_orientation(true, false, true)) =
1099 (n - 1 - i) + (n - 1 - j) * n - local;
1100 // face_orientation=true, face_flip=true, face_rotation=true
1102 local, internal::combined_face_orientation(true, true, true)) =
1103 (n - 1 - j) + i * n - local;
1104 }
1105}
1106
1107
1108
1109template <int dim, int spacedim>
1110unsigned int
1112 const unsigned int face_index,
1113 const unsigned int face,
1115{
1116 AssertIndexRange(face_index, this->n_dofs_per_face(face));
1118
1119 // we need to distinguish between DoFs on vertices, lines and in 3d quads.
1120 // do so in a sequence of if-else statements
1122 // DoF is on a vertex
1123 {
1124 // get the number of the vertex on the face that corresponds to this DoF,
1125 // along with the number of the DoF on this vertex
1126 const unsigned int face_vertex = face_index / this->n_dofs_per_vertex();
1127 const unsigned int dof_index_on_vertex =
1128 face_index % this->n_dofs_per_vertex();
1129
1130 // then get the number of this vertex on the cell and translate
1131 // this to a DoF number on the cell
1134 this->n_dofs_per_vertex() +
1135 dof_index_on_vertex;
1136 }
1138 // DoF is on a face
1139 {
1140 // do the same kind of translation as before. we need to only consider
1141 // DoFs on the lines, i.e., ignoring those on the vertices
1142 const unsigned int index =
1143 face_index - this->get_first_face_line_index(face);
1144
1145 const unsigned int face_line = index / this->n_dofs_per_line();
1146 const unsigned int dof_index_on_line = index % this->n_dofs_per_line();
1147
1148 // we now also need to adjust the line index for the case of
1149 // face orientation, face flips, etc
1150 unsigned int adjusted_dof_index_on_line = 0;
1151 switch (dim)
1152 {
1153 case 1:
1155 break;
1156
1157 case 2:
1160 else
1162 this->n_dofs_per_line() - dof_index_on_line - 1;
1163 break;
1164
1165 case 3:
1166 // in 3d, things are difficult. someone will have to think
1167 // about how this code here should look like, by drawing a bunch
1168 // of pictures of how all the faces can look like with the various
1169 // flips and rotations.
1170 //
1171 // that said, the Q2 case is easy enough to implement, as is the
1172 // case where everything is in standard orientation
1173 Assert((this->n_dofs_per_line() <= 1) ||
1174 combined_orientation ==
1178 break;
1179
1180 default:
1182 }
1183
1184 return (this->get_first_line_index() +
1185 this->reference_cell().face_to_cell_lines(face,
1186 face_line,
1188 this->n_dofs_per_line() +
1189 adjusted_dof_index_on_line);
1190 }
1191 else
1192 // DoF is on a quad
1193 {
1194 Assert(dim >= 3, ExcInternalError());
1195
1196 // ignore vertex and line dofs
1197 const unsigned int index =
1198 face_index - this->get_first_face_quad_index(face);
1199
1200 // the same is true here as above for the 3d case -- someone will
1201 // just have to draw a bunch of pictures. in the meantime,
1202 // we can implement the Q2 case in which it is simple
1203 Assert((this->n_dofs_per_quad(face) <= 1) ||
1206 return (this->get_first_quad_index(face) + index);
1207 }
1208}
1209
1210
1211
1212template <int dim, int spacedim>
1213std::vector<unsigned int>
1215{
1217 AssertThrow(degree > 0, typename FEQ::ExcFEQCannotHaveDegree0());
1218 std::vector<unsigned int> dpo(dim + 1, 1U);
1219 for (unsigned int i = 1; i < dpo.size(); ++i)
1220 dpo[i] = dpo[i - 1] * (degree - 1);
1221 return dpo;
1222}
1223
1224
1225
1226template <int dim, int spacedim>
1227void
1229 const std::vector<Point<1>> &points)
1230{
1232}
1233
1234
1235
1236template <int dim, int spacedim>
1237const FullMatrix<double> &
1239 const unsigned int child,
1240 const RefinementCase<dim> &refinement_case) const
1241{
1242 AssertIndexRange(refinement_case,
1244 Assert(refinement_case != RefinementCase<dim>::no_refinement,
1245 ExcMessage(
1246 "Prolongation matrices are only available for refined cells!"));
1248 child, this->reference_cell().template n_children<dim>(refinement_case));
1249
1250 // initialization upon first request
1251 if (this->prolongation[refinement_case - 1][child].n() == 0)
1252 {
1253 std::lock_guard<std::mutex> lock(prolongation_matrix_mutex);
1254
1255 // if matrix got updated while waiting for the lock
1256 if (this->prolongation[refinement_case - 1][child].n() ==
1257 this->n_dofs_per_cell())
1258 return this->prolongation[refinement_case - 1][child];
1259
1260 // distinguish q/q_dg0 case: only treat Q dofs first
1261 const unsigned int q_dofs_per_cell =
1262 Utilities::fixed_power<dim>(q_degree + 1);
1263
1264 // compute the interpolation matrices in much the same way as we do for
1265 // the constraints. it's actually simpler here, since we don't have this
1266 // weird renumbering stuff going on. The trick is again that we the
1267 // interpolation matrix is formed by a permutation of the indices of the
1268 // cell matrix. The value eps is used a threshold to decide when certain
1269 // evaluations of the Lagrange polynomials are zero or one.
1270 const double eps = 1e-15 * q_degree * dim;
1271
1272# ifdef DEBUG
1273 // in DEBUG mode, check that the evaluation of support points in the
1274 // current numbering gives the identity operation
1275 for (unsigned int i = 0; i < q_dofs_per_cell; ++i)
1276 {
1277 Assert(std::fabs(1. - this->poly_space->compute_value(
1278 i, this->unit_support_points[i])) < eps,
1279 ExcInternalError("The Lagrange polynomial does not evaluate "
1280 "to one or zero in a nodal point. "
1281 "This typically indicates that the "
1282 "polynomial interpolation is "
1283 "ill-conditioned such that round-off "
1284 "prevents the sum to be one."));
1285 for (unsigned int j = 0; j < q_dofs_per_cell; ++j)
1286 if (j != i)
1287 Assert(std::fabs(this->poly_space->compute_value(
1288 i, this->unit_support_points[j])) < eps,
1290 "The Lagrange polynomial does not evaluate "
1291 "to one or zero in a nodal point. "
1292 "This typically indicates that the "
1293 "polynomial interpolation is "
1294 "ill-conditioned such that round-off "
1295 "prevents the sum to be one."));
1296 }
1297# endif
1298
1299 // to efficiently evaluate the polynomial at the subcell, make use of
1300 // the tensor product structure of this element and only evaluate 1d
1301 // information from the polynomial. This makes the cost of this function
1302 // almost negligible also for high order elements
1303 const unsigned int dofs1d = q_degree + 1;
1304 std::vector<Table<2, double>> subcell_evaluations(
1306
1307 const std::vector<unsigned int> &index_map_inverse =
1309
1310 // helper value: step size how to walk through diagonal and how many
1311 // points we have left apart from the first dimension
1312 unsigned int step_size_diag = 0;
1313 {
1314 unsigned int factor = 1;
1315 for (unsigned int d = 0; d < dim; ++d)
1316 {
1317 step_size_diag += factor;
1318 factor *= dofs1d;
1319 }
1320 }
1321
1322 FullMatrix<double> prolongate(this->n_dofs_per_cell(),
1323 this->n_dofs_per_cell());
1324
1325 // go through the points in diagonal to capture variation in all
1326 // directions simultaneously
1327 for (unsigned int j = 0; j < dofs1d; ++j)
1328 {
1329 const unsigned int diag_comp = index_map_inverse[j * step_size_diag];
1330 const Point<dim> p_subcell = this->unit_support_points[diag_comp];
1331 const Point<dim> p_cell =
1333 child,
1334 refinement_case);
1335 for (unsigned int i = 0; i < dofs1d; ++i)
1336 for (unsigned int d = 0; d < dim; ++d)
1337 {
1338 // evaluate along line where only x is different from zero
1340 point[0] = p_cell[d];
1341 const double cell_value =
1342 this->poly_space->compute_value(index_map_inverse[i], point);
1343
1344 // cut off values that are too small. note that we have here
1345 // Lagrange interpolation functions, so they should be zero at
1346 // almost all points, and one at the others, at least on the
1347 // subcells. so set them to their exact values
1348 //
1349 // the actual cut-off value is somewhat fuzzy, but it works
1350 // for 2e-13*degree*dim (see above), which is kind of
1351 // reasonable given that we compute the values of the
1352 // polynomials via an degree-step recursion and then multiply
1353 // the 1d-values. this gives us a linear growth in degree*dim,
1354 // times a small constant.
1355 //
1356 // the embedding matrix is given by applying the inverse of
1357 // the subcell matrix on the cell_interpolation matrix. since
1358 // the subcell matrix is actually only a permutation vector,
1359 // all we need to do is to switch the rows we write the data
1360 // into. moreover, cut off very small values here
1361 if (std::fabs(cell_value) < eps)
1362 subcell_evaluations[d](j, i) = 0;
1363 else
1365 }
1366 }
1367
1368 // now expand from 1d info. block innermost dimension (x_0) in order to
1369 // avoid difficult checks at innermost loop
1370 unsigned int j_indices[dim];
1371 internal::FE_Q_Base::zero_indices<dim>(j_indices);
1372 for (unsigned int j = 0; j < q_dofs_per_cell; j += dofs1d)
1373 {
1374 unsigned int i_indices[dim];
1375 internal::FE_Q_Base::zero_indices<dim>(i_indices);
1376 for (unsigned int i = 0; i < q_dofs_per_cell; i += dofs1d)
1377 {
1378 double val_extra_dim = 1.;
1379 for (unsigned int d = 1; d < dim; ++d)
1380 val_extra_dim *=
1381 subcell_evaluations[d](j_indices[d - 1], i_indices[d - 1]);
1382
1383 // innermost sum where we actually compute. the same as
1384 // prolongate(j,i) = this->poly_space->compute_value (i, p_cell)
1385 for (unsigned int jj = 0; jj < dofs1d; ++jj)
1386 {
1387 const unsigned int j_ind = index_map_inverse[j + jj];
1388 for (unsigned int ii = 0; ii < dofs1d; ++ii)
1389 prolongate(j_ind, index_map_inverse[i + ii]) =
1391 }
1392
1393 // update indices that denote the tensor product position. a bit
1394 // fuzzy and therefore not done for innermost x_0 direction
1395 internal::FE_Q_Base::increment_indices<dim>(i_indices, dofs1d);
1396 }
1397 Assert(i_indices[dim - 1] == 1, ExcInternalError());
1398 internal::FE_Q_Base::increment_indices<dim>(j_indices, dofs1d);
1399 }
1400
1401 // the discontinuous node is simply mapped on the discontinuous node on
1402 // the child element
1404 prolongate(q_dofs_per_cell, q_dofs_per_cell) = 1.;
1405
1406 // and make sure that the row sum is 1. this must be so since for this
1407 // element, the shape functions add up to one
1408# ifdef DEBUG
1409 for (unsigned int row = 0; row < this->n_dofs_per_cell(); ++row)
1410 {
1411 double sum = 0;
1412 for (unsigned int col = 0; col < this->n_dofs_per_cell(); ++col)
1413 sum += prolongate(row, col);
1414 Assert(std::fabs(sum - 1.) <
1415 std::max(eps, 5e-16 * std::sqrt(this->n_dofs_per_cell())),
1416 ExcInternalError("The entries in a row of the local "
1417 "prolongation matrix do not add to one. "
1418 "This typically indicates that the "
1419 "polynomial interpolation is "
1420 "ill-conditioned such that round-off "
1421 "prevents the sum to be one."));
1422 }
1423# endif
1424
1425 // move result into place
1426 const_cast<FullMatrix<double> &>(
1427 this->prolongation[refinement_case - 1][child]) = std::move(prolongate);
1428 }
1429
1430 // finally return the matrix
1431 return this->prolongation[refinement_case - 1][child];
1432}
1433
1434
1435
1436template <int dim, int spacedim>
1437const FullMatrix<double> &
1439 const unsigned int child,
1440 const RefinementCase<dim> &refinement_case) const
1441{
1442 AssertIndexRange(refinement_case,
1444 Assert(refinement_case != RefinementCase<dim>::no_refinement,
1445 ExcMessage(
1446 "Restriction matrices are only available for refined cells!"));
1448 child, this->reference_cell().template n_children<dim>(refinement_case));
1449
1450 // initialization upon first request
1451 if (this->restriction[refinement_case - 1][child].n() == 0)
1452 {
1453 std::lock_guard<std::mutex> lock(restriction_matrix_mutex);
1454
1455 // if matrix got updated while waiting for the lock...
1456 if (this->restriction[refinement_case - 1][child].n() ==
1457 this->n_dofs_per_cell())
1458 return this->restriction[refinement_case - 1][child];
1459
1461 this->n_dofs_per_cell());
1462 // distinguish q/q_dg0 case
1463 const unsigned int q_dofs_per_cell =
1464 Utilities::fixed_power<dim>(q_degree + 1);
1465
1466 // for Lagrange interpolation polynomials based on equidistant points,
1467 // construction of the restriction matrices is relatively simple. the
1468 // reason is that in this case the interpolation points on the mother
1469 // cell are always also interpolation points for some shape function on
1470 // one or the other child.
1471 //
1472 // in the general case with non-equidistant points, we need to actually
1473 // do an interpolation. thus, we take the interpolation points on the
1474 // mother cell and evaluate the shape functions of the child cell on
1475 // those points. it does not hurt in the equidistant case because then
1476 // simple one shape function evaluates to one and the others to zero.
1477 //
1478 // this element is non-additive in all its degrees of freedom by
1479 // default, which requires care in downstream use. fortunately, even the
1480 // interpolation on non-equidistant points is invariant under the
1481 // assumption that whenever a row makes a non-zero contribution to the
1482 // mother's residual, the correct value is interpolated.
1483
1484 const double eps = 1e-15 * q_degree * dim;
1485 const std::vector<unsigned int> &index_map_inverse =
1487
1488 const unsigned int dofs1d = q_degree + 1;
1489 std::vector<Tensor<1, dim>> evaluations1d(dofs1d);
1490
1491 my_restriction.reinit(this->n_dofs_per_cell(), this->n_dofs_per_cell());
1492
1493 for (unsigned int i = 0; i < q_dofs_per_cell; ++i)
1494 {
1495 unsigned int mother_dof = index_map_inverse[i];
1496 const Point<dim> p_cell = this->unit_support_points[mother_dof];
1497
1498 // check whether this interpolation point is inside this child cell
1499 const Point<dim> p_subcell =
1501 child,
1502 refinement_case);
1504 {
1505 // same logic as in initialize_embedding to evaluate the
1506 // polynomial faster than from the tensor product: since we
1507 // evaluate all polynomials, it is much faster to just compute
1508 // the 1d values for all polynomials before and then get the
1509 // dim-data.
1510 for (unsigned int j = 0; j < dofs1d; ++j)
1511 for (unsigned int d = 0; d < dim; ++d)
1512 {
1514 point[0] = p_subcell[d];
1515 evaluations1d[j][d] =
1516 this->poly_space->compute_value(index_map_inverse[j],
1517 point);
1518 }
1519 unsigned int j_indices[dim];
1520 internal::FE_Q_Base::zero_indices<dim>(j_indices);
1521 double sum_check = 0;
1522 for (unsigned int j = 0; j < q_dofs_per_cell; j += dofs1d)
1523 {
1524 double val_extra_dim = 1.;
1525 for (unsigned int d = 1; d < dim; ++d)
1527 for (unsigned int jj = 0; jj < dofs1d; ++jj)
1528 {
1529 // find the child shape function(s) corresponding to
1530 // this point. Usually this is just one function;
1531 // however, when we use FE_Q on arbitrary nodes a parent
1532 // support point might not be a child support point, and
1533 // then we will get more than one nonzero value per
1534 // row. Still, the values should sum up to 1
1535 const double val = val_extra_dim * evaluations1d[jj][0];
1536 const unsigned int child_dof = index_map_inverse[j + jj];
1537 if (std::fabs(val - 1.) < eps)
1539 else if (std::fabs(val) > eps)
1541 sum_check += val;
1542 }
1543 internal::FE_Q_Base::increment_indices<dim>(j_indices,
1544 dofs1d);
1545 }
1546 (void)sum_check;
1547 Assert(std::fabs(sum_check - 1.0) <
1548 std::max(eps,
1549 5e-16 * std::sqrt(this->n_dofs_per_cell())),
1550 ExcInternalError("The entries in a row of the local "
1551 "restriction matrix do not add to one. "
1552 "This typically indicates that the "
1553 "polynomial interpolation is "
1554 "ill-conditioned such that round-off "
1555 "prevents the sum to be one."));
1556 }
1557
1558 // part for FE_Q_DG0
1560 my_restriction(this->n_dofs_per_cell() - 1,
1561 this->n_dofs_per_cell() - 1) =
1562 1. / this->reference_cell().template n_children<dim>(
1563 RefinementCase<dim>(refinement_case));
1564 }
1565
1566 // move result into place
1567 const_cast<FullMatrix<double> &>(
1568 this->restriction[refinement_case - 1][child]) =
1569 std::move(my_restriction);
1570 }
1571
1572 return this->restriction[refinement_case - 1][child];
1573}
1574
1575
1576
1577//---------------------------------------------------------------------------
1578// Data field initialization
1579//---------------------------------------------------------------------------
1580
1581
1582template <int dim, int spacedim>
1583bool
1585 const unsigned int shape_index,
1586 const unsigned int face_index) const
1587{
1590
1591 // in 1d, things are simple. since there is only one degree of freedom per
1592 // vertex in this class, the first is on vertex 0 (==face 0 in some sense),
1593 // the second on face 1:
1594 if (dim == 1)
1595 return (((shape_index == 0) && (face_index == 0)) ||
1596 ((shape_index == 1) && (face_index == 1)));
1597
1598 // first, special-case interior shape functions, since they have no support
1599 // no-where on the boundary
1600 if (((dim == 2) &&
1601 (shape_index >= this->get_first_quad_index(0 /*first quad*/))) ||
1602 ((dim == 3) && (shape_index >= this->get_first_hex_index())))
1603 return false;
1604
1605 // let's see whether this is a vertex
1607 {
1608 // for Q elements, there is one dof per vertex, so
1609 // shape_index==vertex_number. check whether this vertex is on the given
1610 // face. thus, for each face, give a list of vertices
1611 const unsigned int vertex_no = shape_index;
1614
1615 for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face; ++v)
1616 if (GeometryInfo<dim>::face_to_cell_vertices(face_index, v) ==
1617 vertex_no)
1618 return true;
1619
1620 return false;
1621 }
1622 else if (shape_index < this->get_first_quad_index(0 /*first quad*/))
1623 // ok, dof is on a line
1624 {
1625 const unsigned int line_index =
1626 (shape_index - this->get_first_line_index()) / this->n_dofs_per_line();
1629
1630 // in 2d, the line is the face, so get the line index
1631 if constexpr (dim == 2)
1632 return (line_index == face_index);
1633 else if constexpr (dim == 3)
1634 {
1635 // see whether the given line is on the given face.
1636 for (unsigned int l = 0; l < GeometryInfo<3>::lines_per_face; ++l)
1637 if (GeometryInfo<3>::face_to_cell_lines(face_index, l) ==
1638 line_index)
1639 return true;
1640
1641 return false;
1642 }
1643 else
1645 }
1647 // dof is on a quad
1648 {
1649 const unsigned int quad_index =
1650 (shape_index - this->get_first_quad_index(0)) /
1651 this->n_dofs_per_quad(face_index); // this won't work
1653
1654 // in 2d, cell bubble are zero on all faces. but we have treated this
1655 // case above already
1656 Assert(dim != 2, ExcInternalError());
1657
1658 // in 3d, quad_index=face_index
1659 if (dim == 3)
1660 return (quad_index == face_index);
1661 else
1663 }
1664 else
1665 // dof on hex
1666 {
1667 // can only happen in 3d, but this case has already been covered above
1669 return false;
1670 }
1671
1672 // we should not have gotten here
1674 return false;
1675}
1676
1677
1678
1679template <int dim, int spacedim>
1680std::pair<Table<2, bool>, std::vector<unsigned int>>
1682{
1683 Table<2, bool> constant_modes(1, this->n_dofs_per_cell());
1684 // We here just care for the constant mode due to the polynomial space
1685 // without any enrichments
1686 // As there may be more constant modes derived classes may to implement this
1687 // themselves. An example for this is FE_Q_DG0.
1688 for (unsigned int i = 0; i < Utilities::fixed_power<dim>(q_degree + 1); ++i)
1689 constant_modes(0, i) = true;
1690 return std::pair<Table<2, bool>, std::vector<unsigned int>>(
1691 constant_modes, std::vector<unsigned int>(1, 0));
1692}
1693
1694#endif
1695
1696// explicit instantiations
1697#include "fe/fe_q_base.inst"
1698
std::vector< unsigned int > get_poly_space_numbering_inverse() const
virtual double shape_value(const unsigned int i, const Point< dim > &p) const override
const std::unique_ptr< ScalarPolynomialsBase< dim > > poly_space
Definition fe_poly.h:532
Threads::Mutex prolongation_matrix_mutex
Definition fe_q_base.h:333
void initialize_unit_face_support_points(const std::vector< Point< 1 > > &points)
virtual bool hp_constraints_are_implemented() const override
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_line_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
virtual void get_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix) const override
Threads::Mutex restriction_matrix_mutex
Definition fe_q_base.h:332
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_quad_dof_identities(const FiniteElement< dim, spacedim > &fe_other, const unsigned int face_no=0) const override
const unsigned int q_degree
Definition fe_q_base.h:340
virtual std::pair< Table< 2, bool >, std::vector< unsigned int > > get_constant_modes() const override
void initialize(const std::vector< Point< 1 > > &support_points_1d)
void initialize_unit_support_points(const std::vector< Point< 1 > > &points)
virtual void get_subface_interpolation_matrix(const FiniteElement< dim, spacedim > &source, const unsigned int subface, FullMatrix< double > &matrix, const unsigned int face_no=0) const override
virtual void get_face_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix, const unsigned int face_no=0) const override
static std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
void initialize_constraints(const std::vector< Point< 1 > > &points)
void initialize_dof_index_permutations()
virtual const FullMatrix< double > & get_restriction_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const override
virtual const FullMatrix< double > & get_prolongation_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const override
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const override
virtual unsigned int face_to_cell_index(const unsigned int face_dof_index, const unsigned int face, const types::geometric_orientation combined_orientation=numbers::default_geometric_orientation) const override
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_vertex_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
FE_Q_Base(const ScalarPolynomialsBase< dim > &poly_space, const FiniteElementData< dim > &fe_data, const std::vector< bool > &restriction_is_additive_flags)
unsigned int get_first_line_index() const
unsigned int n_dofs_per_vertex() const
const unsigned int degree
Definition fe_data.h:452
unsigned int n_dofs_per_cell() const
unsigned int n_dofs_per_line() const
unsigned int get_first_quad_index(const unsigned int quad_no=0) const
unsigned int n_dofs_per_face(unsigned int face_no=0, unsigned int child=0) const
ReferenceCell reference_cell() const
unsigned int get_first_face_line_index(const unsigned int face_no=0) const
unsigned int get_first_face_quad_index(const unsigned int face_no=0) const
unsigned int n_unique_faces() const
unsigned int n_dofs_per_quad(unsigned int face_no=0) const
unsigned int get_first_hex_index() const
std::vector< std::vector< Point< dim - 1 > > > unit_face_support_points
Definition fe.h:2568
std::vector< std::vector< FullMatrix< double > > > restriction
Definition fe.h:2523
std::vector< Table< 2, int > > adjust_quad_dof_index_for_face_orientation_table
Definition fe.h:2597
std::vector< int > adjust_line_dof_index_for_line_orientation_table
Definition fe.h:2610
std::vector< Point< dim > > unit_support_points
Definition fe.h:2561
std::vector< std::vector< FullMatrix< double > > > prolongation
Definition fe.h:2537
Definition point.h:113
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)
unsigned int face_to_cell_vertices(const unsigned int face, const unsigned int vertex, const types::geometric_orientation face_orientation) const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:518
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:519
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcInterpolationNotImplemented()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Task< RT > new_task(const std::function< RT()> &function)
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
std::vector< unsigned int > lexicographic_to_hierarchic_numbering(unsigned int degree)
constexpr char U
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:193
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
T sum(const T &t, const MPI_Comm mpi_communicator)
types::geometric_orientation combined_face_orientation(const bool face_orientation, const bool face_rotation, const bool face_flip)
constexpr unsigned int invalid_unsigned_int
Definition types.h:232
constexpr types::geometric_orientation default_geometric_orientation
Definition types.h:346
STL namespace.
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
static void initialize_constraints(const std::vector< Point< 1 > > &, FE_Q_Base< 3, spacedim > &fe)
Definition fe_q_base.cc:212
static void initialize_constraints(const std::vector< Point< 1 > > &, FE_Q_Base< 1, spacedim > &)
Definition fe_q_base.cc:97
static void initialize_constraints(const std::vector< Point< 1 > > &, FE_Q_Base< 2, spacedim > &fe)
Definition fe_q_base.cc:106
static Point< dim > cell_to_child_coordinates(const Point< dim > &p, const unsigned int child_index, const RefinementCase< dim > refine_case=RefinementCase< dim >::isotropic_refinement)
static Point< dim > child_to_cell_coordinates(const Point< dim > &p, const unsigned int child_index, const RefinementCase< dim > refine_case=RefinementCase< dim >::isotropic_refinement)