Reference documentation for deal.II version 9.5.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
fe_q_base.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2000 - 2023 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16
26
27#include <deal.II/fe/fe_dgp.h>
28#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>
89struct FE_Q_Base<xdim, xspacedim>::Implementation
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 if (q_deg > 1)
165 {
166 const unsigned int n = q_deg - 1;
167 const double step = 1. / q_deg;
168 // subface 0
169 for (unsigned int i = 1; i <= n; ++i)
170 constraint_points.push_back(
172 Point<dim - 1>(i * step), 0));
173 // subface 1
174 for (unsigned int i = 1; i <= n; ++i)
175 constraint_points.push_back(
177 Point<dim - 1>(i * step), 1));
178 }
179
180 // Now construct relation between destination (child) and source (mother)
181 // dofs.
182
183 fe.interface_constraints.TableBase<2, double>::reinit(
185
186 // use that the element evaluates to 1 at index 0 and along the line at
187 // zero
188 const std::vector<unsigned int> &index_map_inverse =
190 const std::vector<unsigned int> face_index_map =
193 fe.poly_space->compute_value(index_map_inverse[0], Point<dim>()) -
194 1.) < 1e-14,
196
197 for (unsigned int i = 0; i < constraint_points.size(); ++i)
198 for (unsigned int j = 0; j < q_deg + 1; ++j)
199 {
200 Point<dim> p;
201 p[0] = constraint_points[i](0);
202 fe.interface_constraints(i, face_index_map[j]) =
203 fe.poly_space->compute_value(index_map_inverse[j], p);
204
205 // if the value is small up to round-off, then simply set it to zero
206 // to avoid unwanted fill-in of the constraint matrices (which would
207 // then increase the number of other DoFs a constrained DoF would
208 // couple to)
209 if (std::fabs(fe.interface_constraints(i, face_index_map[j])) < 1e-13)
210 fe.interface_constraints(i, face_index_map[j]) = 0;
211 }
212 }
213
214
215 template <int spacedim>
216 static void
217 initialize_constraints(const std::vector<Point<1>> & /*points*/,
219 {
220 const unsigned int dim = 3;
221
222 unsigned int q_deg = fe.degree;
223 if (dynamic_cast<const TensorProductPolynomialsBubbles<dim> *>(
224 &fe.get_poly_space()) != nullptr)
225 q_deg = fe.degree - 1;
226
227 // For a detailed documentation of the interpolation see the
228 // FE_Q_Base<2>::initialize_constraints function.
229
230 // In the following the points x_i are constructed in the order as
231 // described in the documentation of the FiniteElement class (fe_base.h),
232 // i.e.
233 // *--15--4--16--*
234 // | | |
235 // 10 19 6 20 12
236 // | | |
237 // 1--7---0--8---2
238 // | | |
239 // 9 17 5 18 11
240 // | | |
241 // *--13--3--14--*
242 std::vector<Point<dim - 1>> constraint_points;
243
244 // Add midpoint
245 constraint_points.emplace_back(0.5, 0.5);
246
247 // Add midpoints of lines of "mother-face"
248 constraint_points.emplace_back(0, 0.5);
249 constraint_points.emplace_back(1, 0.5);
250 constraint_points.emplace_back(0.5, 0);
251 constraint_points.emplace_back(0.5, 1);
252
253 if (q_deg > 1)
254 {
255 const unsigned int n = q_deg - 1;
256 const double step = 1. / q_deg;
257 std::vector<Point<dim - 2>> line_support_points(n);
258 for (unsigned int i = 0; i < n; ++i)
259 line_support_points[i](0) = (i + 1) * step;
260 Quadrature<dim - 2> qline(line_support_points);
261
262 // auxiliary points in 2d
263 std::vector<Point<dim - 1>> p_line(n);
264
265 // Add nodes of lines interior in the "mother-face"
266
267 // line 5: use line 9
269 ReferenceCells::get_hypercube<dim - 1>(), qline, 0, 0, p_line);
270 for (unsigned int i = 0; i < n; ++i)
271 constraint_points.push_back(p_line[i] + Point<dim - 1>(0.5, 0));
272 // line 6: use line 10
274 ReferenceCells::get_hypercube<dim - 1>(), qline, 0, 1, p_line);
275 for (unsigned int i = 0; i < n; ++i)
276 constraint_points.push_back(p_line[i] + Point<dim - 1>(0.5, 0));
277 // line 7: use line 13
279 ReferenceCells::get_hypercube<dim - 1>(), qline, 2, 0, p_line);
280 for (unsigned int i = 0; i < n; ++i)
281 constraint_points.push_back(p_line[i] + Point<dim - 1>(0, 0.5));
282 // line 8: use line 14
284 ReferenceCells::get_hypercube<dim - 1>(), qline, 2, 1, p_line);
285 for (unsigned int i = 0; i < n; ++i)
286 constraint_points.push_back(p_line[i] + Point<dim - 1>(0, 0.5));
287
288 // DoFs on bordering lines lines 9-16
289 for (unsigned int face = 0;
290 face < GeometryInfo<dim - 1>::faces_per_cell;
291 ++face)
292 for (unsigned int subface = 0;
293 subface < GeometryInfo<dim - 1>::max_children_per_face;
294 ++subface)
295 {
297 ReferenceCells::get_hypercube<dim - 1>(),
298 qline,
299 face,
300 subface,
301 p_line);
302 constraint_points.insert(constraint_points.end(),
303 p_line.begin(),
304 p_line.end());
305 }
306
307 // Create constraints for interior nodes
308 std::vector<Point<dim - 1>> inner_points(n * n);
309 for (unsigned int i = 0, iy = 1; iy <= n; ++iy)
310 for (unsigned int ix = 1; ix <= n; ++ix)
311 inner_points[i++] = Point<dim - 1>(ix * step, iy * step);
312
313 // at the moment do this for isotropic face refinement only
314 for (unsigned int child = 0;
315 child < GeometryInfo<dim - 1>::max_children_per_cell;
316 ++child)
317 for (const auto &inner_point : inner_points)
318 constraint_points.push_back(
320 child));
321 }
322
323 // Now construct relation between destination (child) and source (mother)
324 // dofs.
325 const unsigned int pnts = (q_deg + 1) * (q_deg + 1);
326
327 // use that the element evaluates to 1 at index 0 and along the line at
328 // zero
329 const std::vector<unsigned int> &index_map_inverse =
331 const std::vector<unsigned int> face_index_map =
334 fe.poly_space->compute_value(index_map_inverse[0], Point<dim>()) -
335 1.) < 1e-14,
337
338 fe.interface_constraints.TableBase<2, double>::reinit(
340
341 for (unsigned int i = 0; i < constraint_points.size(); ++i)
342 {
343 const double interval = static_cast<double>(q_deg * 2);
344 bool mirror[dim - 1];
345 Point<dim> constraint_point;
346
347 // Eliminate FP errors in constraint points. Due to their origin, they
348 // must all be fractions of the unit interval. If we have polynomial
349 // degree 4, the refined element has 8 intervals. Hence the
350 // coordinates must be 0, 0.125, 0.25, 0.375 etc. Now the coordinates
351 // of the constraint points will be multiplied by the inverse of the
352 // interval size (in the example by 8). After that the coordinates
353 // must be integral numbers. Hence a normal truncation is performed
354 // and the coordinates will be scaled back. The equal treatment of all
355 // coordinates should eliminate any FP errors.
356 for (unsigned int k = 0; k < dim - 1; ++k)
357 {
358 const int coord_int =
359 static_cast<int>(constraint_points[i](k) * interval + 0.25);
360 constraint_point(k) = 1. * coord_int / interval;
361
362 // The following lines of code should eliminate the problems with
363 // the constraints object which appeared for P>=4. The
364 // AffineConstraints class complained about different constraints
365 // for the same entry: Actually, this
366 // difference could be attributed to FP errors, as it was in the
367 // range of 1.0e-16. These errors originate in the loss of
368 // symmetry in the FP approximation of the shape-functions.
369 // Considering a 3rd order shape function in 1d, we have
370 // N0(x)=N3(1-x) and N1(x)=N2(1-x). For higher order polynomials
371 // the FP approximations of the shape functions do not satisfy
372 // these equations any more! Thus in the following code
373 // everything is computed in the interval x \in [0..0.5], which is
374 // sufficient to express all values that could come out from a
375 // computation of any shape function in the full interval
376 // [0..1]. If x > 0.5 the computation is done for 1-x with the
377 // shape function N_{p-n} instead of N_n. Hence symmetry is
378 // preserved and everything works fine...
379 //
380 // For a different explanation of the problem, see the discussion
381 // in the FiniteElement class for constraint matrices in 3d.
382 mirror[k] = (constraint_point(k) > 0.5);
383 if (mirror[k])
384 constraint_point(k) = 1.0 - constraint_point(k);
385 }
386
387 for (unsigned int j = 0; j < pnts; ++j)
388 {
389 unsigned int indices[2] = {j % (q_deg + 1), j / (q_deg + 1)};
390
391 for (unsigned int k = 0; k < 2; ++k)
392 if (mirror[k])
393 indices[k] = q_deg - indices[k];
394
395 const unsigned int new_index =
396 indices[1] * (q_deg + 1) + indices[0];
397
398 fe.interface_constraints(i, face_index_map[j]) =
399 fe.poly_space->compute_value(index_map_inverse[new_index],
400 constraint_point);
401
402 // if the value is small up to round-off, then simply set it to
403 // zero to avoid unwanted fill-in of the constraint matrices
404 // (which would then increase the number of other DoFs a
405 // constrained DoF would couple to)
406 if (std::fabs(fe.interface_constraints(i, face_index_map[j])) <
407 1e-13)
408 fe.interface_constraints(i, face_index_map[j]) = 0;
409 }
410 }
411 }
412};
413
414
415
416template <int dim, int spacedim>
418 const ScalarPolynomialsBase<dim> &poly_space,
419 const FiniteElementData<dim> & fe_data,
420 const std::vector<bool> & restriction_is_additive_flags)
421 : FE_Poly<dim, spacedim>(
422 poly_space,
423 fe_data,
424 restriction_is_additive_flags,
425 std::vector<ComponentMask>(1, std::vector<bool>(1, true)))
426 , q_degree(dynamic_cast<const TensorProductPolynomialsBubbles<dim> *>(
427 &poly_space) != nullptr ?
428 this->degree - 1 :
429 this->degree)
430{}
431
432
433
434template <int dim, int spacedim>
435void
436FE_Q_Base<dim, spacedim>::initialize(const std::vector<Point<1>> &points)
437{
438 Assert(points[0][0] == 0,
439 ExcMessage("The first support point has to be zero."));
440 Assert(points.back()[0] == 1,
441 ExcMessage("The last support point has to be one."));
442
443 // distinguish q/q_dg0 case: need to be flexible enough to allow more
444 // degrees of freedom than there are FE_Q degrees of freedom for derived
445 // class FE_Q_DG0 that otherwise shares 95% of the code.
446 const unsigned int q_dofs_per_cell =
447 Utilities::fixed_power<dim>(q_degree + 1);
448 Assert(q_dofs_per_cell == this->n_dofs_per_cell() ||
449 q_dofs_per_cell + 1 == this->n_dofs_per_cell() ||
450 q_dofs_per_cell + dim == this->n_dofs_per_cell(),
452
453 [this, q_dofs_per_cell]() {
454 std::vector<unsigned int> renumber =
455 FETools::hierarchic_to_lexicographic_numbering<dim>(q_degree);
456 for (unsigned int i = q_dofs_per_cell; i < this->n_dofs_per_cell(); ++i)
457 renumber.push_back(i);
458 auto *tensor_poly_space_ptr =
459 dynamic_cast<TensorProductPolynomials<dim> *>(this->poly_space.get());
460 if (tensor_poly_space_ptr != nullptr)
461 {
462 tensor_poly_space_ptr->set_numbering(renumber);
463 return;
464 }
465 auto *tensor_piecewise_poly_space_ptr = dynamic_cast<
467 *>(this->poly_space.get());
468 if (tensor_piecewise_poly_space_ptr != nullptr)
469 {
470 tensor_piecewise_poly_space_ptr->set_numbering(renumber);
471 return;
472 }
473 auto *tensor_bubbles_poly_space_ptr =
475 this->poly_space.get());
476 if (tensor_bubbles_poly_space_ptr != nullptr)
477 {
478 tensor_bubbles_poly_space_ptr->set_numbering(renumber);
479 return;
480 }
481 auto *tensor_const_poly_space_ptr =
483 this->poly_space.get());
484 if (tensor_const_poly_space_ptr != nullptr)
485 {
486 tensor_const_poly_space_ptr->set_numbering(renumber);
487 return;
488 }
489 Assert(false, ExcNotImplemented());
490 }();
491
492 // Finally fill in support points on cell and face and initialize
493 // constraints. All of this can happen in parallel
495 tasks += Threads::new_task([&]() { initialize_unit_support_points(points); });
496 tasks +=
498 tasks += Threads::new_task([&]() { initialize_constraints(points); });
499 tasks +=
501 tasks.join_all();
502
503 // do not initialize embedding and restriction here. these matrices are
504 // initialized on demand in get_restriction_matrix and
505 // get_prolongation_matrix
506}
507
508
509
510template <int dim, int spacedim>
511void
513 const FiniteElement<dim, spacedim> &x_source_fe,
514 FullMatrix<double> & interpolation_matrix) const
515{
516 // go through the list of elements we can interpolate from
517 if (const FE_Q_Base<dim, spacedim> *source_fe =
518 dynamic_cast<const FE_Q_Base<dim, spacedim> *>(&x_source_fe))
519 {
520 // ok, source is a Q element, so we will be able to do the work
521 Assert(interpolation_matrix.m() == this->n_dofs_per_cell(),
522 ExcDimensionMismatch(interpolation_matrix.m(),
523 this->n_dofs_per_cell()));
524 Assert(interpolation_matrix.n() == x_source_fe.n_dofs_per_cell(),
525 ExcDimensionMismatch(interpolation_matrix.m(),
526 x_source_fe.n_dofs_per_cell()));
527
528 // only evaluate Q dofs
529 const unsigned int q_dofs_per_cell =
530 Utilities::fixed_power<dim>(q_degree + 1);
531 const unsigned int source_q_dofs_per_cell =
532 Utilities::fixed_power<dim>(source_fe->degree + 1);
533
534 // evaluation is simply done by evaluating the other FE's basis functions
535 // on the unit support points (FE_Q has the property that the cell
536 // interpolation matrix is a unit matrix, so no need to evaluate it and
537 // invert it)
538 for (unsigned int j = 0; j < q_dofs_per_cell; ++j)
539 {
540 // read in a point on this cell and evaluate the shape functions there
541 const Point<dim> p = this->unit_support_points[j];
542
543 // FE_Q element evaluates to 1 in unit support point and to zero in
544 // all other points by construction
545 Assert(std::abs(this->poly_space->compute_value(j, p) - 1.) < 1e-13,
547
548 for (unsigned int i = 0; i < source_q_dofs_per_cell; ++i)
549 interpolation_matrix(j, i) =
550 source_fe->poly_space->compute_value(i, p);
551 }
552
553 // for FE_Q_DG0, add one last row of identity
554 if (q_dofs_per_cell < this->n_dofs_per_cell())
555 {
556 AssertDimension(source_q_dofs_per_cell + 1,
557 source_fe->n_dofs_per_cell());
558 for (unsigned int i = 0; i < source_q_dofs_per_cell; ++i)
559 interpolation_matrix(q_dofs_per_cell, i) = 0.;
560 for (unsigned int j = 0; j < q_dofs_per_cell; ++j)
561 interpolation_matrix(j, source_q_dofs_per_cell) = 0.;
562 interpolation_matrix(q_dofs_per_cell, source_q_dofs_per_cell) = 1.;
563 }
564
565 // cut off very small values
566 const double eps = 2e-13 * q_degree * dim;
567 for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
568 for (unsigned int j = 0; j < source_fe->n_dofs_per_cell(); ++j)
569 if (std::fabs(interpolation_matrix(i, j)) < eps)
570 interpolation_matrix(i, j) = 0.;
571
572#ifdef DEBUG
573 // make sure that the row sum of each of the matrices is 1 at this
574 // point. this must be so since the shape functions sum up to 1
575 for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
576 {
577 double sum = 0.;
578 for (unsigned int j = 0; j < source_fe->n_dofs_per_cell(); ++j)
579 sum += interpolation_matrix(i, j);
580
581 Assert(std::fabs(sum - 1) < eps, ExcInternalError());
582 }
583#endif
584 }
585 else if (dynamic_cast<const FE_Nothing<dim> *>(&x_source_fe))
586 {
587 // the element we want to interpolate from is an FE_Nothing. this
588 // element represents a function that is constant zero and has no
589 // degrees of freedom, so the interpolation is simply a multiplication
590 // with a n_dofs x 0 matrix. there is nothing to do here
591
592 // we would like to verify that the number of rows and columns of
593 // the matrix equals this->n_dofs_per_cell() and zero. unfortunately,
594 // whenever we do FullMatrix::reinit(m,0), it sets both rows and
595 // columns to zero, instead of m and zero. thus, only test the
596 // number of columns
597 Assert(interpolation_matrix.n() == x_source_fe.n_dofs_per_cell(),
598 ExcDimensionMismatch(interpolation_matrix.m(),
599 x_source_fe.n_dofs_per_cell()));
600 }
601 else
603 false,
604 (typename FiniteElement<dim,
606}
607
608
609
610template <int dim, int spacedim>
611void
613 const FiniteElement<dim, spacedim> &source_fe,
614 FullMatrix<double> & interpolation_matrix,
615 const unsigned int face_no) const
616{
619 interpolation_matrix,
620 face_no);
621}
622
623
624
625template <int dim, int spacedim>
626void
628 const FiniteElement<dim, spacedim> &source_fe,
629 const unsigned int subface,
630 FullMatrix<double> & interpolation_matrix,
631 const unsigned int face_no) const
632{
633 Assert(interpolation_matrix.m() == source_fe.n_dofs_per_face(face_no),
634 ExcDimensionMismatch(interpolation_matrix.m(),
635 source_fe.n_dofs_per_face(face_no)));
636
637 Assert(source_fe.n_components() == this->n_components(),
638 ExcDimensionMismatch(source_fe.n_components(), this->n_components()));
639
640 if (source_fe.has_face_support_points(face_no))
641 {
642 // have this test in here since a table of size 2x0 reports its size as
643 // 0x0
644 Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
645 ExcDimensionMismatch(interpolation_matrix.n(),
646 this->n_dofs_per_face(face_no)));
647
648 // Make sure that the element for which the DoFs should be constrained
649 // is the one with the higher polynomial degree. Actually the procedure
650 // will work also if this assertion is not satisfied. But the matrices
651 // produced in that case might lead to problems in the hp-procedures,
652 // which use this method.
653 Assert(
654 this->n_dofs_per_face(face_no) <= source_fe.n_dofs_per_face(face_no),
655 (typename FiniteElement<dim,
657
658 // generate a point on this cell and evaluate the shape functions there
659 const Quadrature<dim - 1> quad_face_support(
660 source_fe.get_unit_face_support_points(face_no));
661
662 // Rule of thumb for FP accuracy, that can be expected for a given
663 // polynomial degree. This value is used to cut off values close to
664 // zero.
665 const double eps = 2e-13 * this->q_degree * std::max(dim - 1, 1);
666
667 // compute the interpolation matrix by simply taking the value at the
668 // support points.
669 // TODO: Verify that all faces are the same with respect to
670 // these support points. Furthermore, check if something has to
671 // be done for the face orientation flag in 3d.
672 const Quadrature<dim> subface_quadrature =
675 quad_face_support,
676 0) :
677 QProjector<dim>::project_to_subface(this->reference_cell(),
678 quad_face_support,
679 0,
680 subface);
681 for (unsigned int i = 0; i < source_fe.n_dofs_per_face(face_no); ++i)
682 {
683 const Point<dim> &p = subface_quadrature.point(i);
684
685 for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
686 {
687 double matrix_entry =
688 this->shape_value(this->face_to_cell_index(j, 0), p);
689
690 // Correct the interpolated value. I.e. if it is close to 1 or
691 // 0, make it exactly 1 or 0. Unfortunately, this is required to
692 // avoid problems with higher order elements.
693 if (std::fabs(matrix_entry - 1.0) < eps)
694 matrix_entry = 1.0;
695 if (std::fabs(matrix_entry) < eps)
696 matrix_entry = 0.0;
697
698 interpolation_matrix(i, j) = matrix_entry;
699 }
700 }
701
702#ifdef DEBUG
703 // make sure that the row sum of each of the matrices is 1 at this
704 // point. this must be so since the shape functions sum up to 1
705 for (unsigned int j = 0; j < source_fe.n_dofs_per_face(face_no); ++j)
706 {
707 double sum = 0.;
708
709 for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
710 sum += interpolation_matrix(j, i);
711
712 Assert(std::fabs(sum - 1) < eps, ExcInternalError());
713 }
714#endif
715 }
716 else if (dynamic_cast<const FE_Nothing<dim> *>(&source_fe) != nullptr)
717 {
718 // nothing to do here, the FE_Nothing has no degrees of freedom anyway
719 }
720 else
722 false,
723 (typename FiniteElement<dim,
725}
726
727
728
729template <int dim, int spacedim>
730bool
732{
733 return true;
734}
735
736
737
738template <int dim, int spacedim>
739std::vector<std::pair<unsigned int, unsigned int>>
741 const FiniteElement<dim, spacedim> &fe_other) const
742{
743 if (dynamic_cast<const FE_Q_Base<dim, spacedim> *>(&fe_other) != nullptr)
744 {
745 // there should be exactly one single DoF of each FE at a vertex, and they
746 // should have identical value
747 return {{0U, 0U}};
748 }
749 else if (dynamic_cast<const FE_SimplexP<dim, spacedim> *>(&fe_other) !=
750 nullptr)
751 {
752 // there should be exactly one single DoF of each FE at a vertex, and they
753 // should have identical value
754 return {{0U, 0U}};
755 }
756 else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
757 {
758 // the FE_Nothing has no degrees of freedom, so there are no
759 // equivalencies to be recorded
760 return {};
761 }
762 else if (fe_other.n_unique_faces() == 1 && fe_other.n_dofs_per_face(0) == 0)
763 {
764 // if the other element has no elements on faces at all,
765 // then it would be impossible to enforce any kind of
766 // continuity even if we knew exactly what kind of element
767 // we have -- simply because the other element declares
768 // that it is discontinuous because it has no DoFs on
769 // its faces. in that case, just state that we have no
770 // constraints to declare
771 return {};
772 }
773 else
774 {
775 Assert(false, ExcNotImplemented());
776 return {};
777 }
778}
779
780
781
782template <int dim, int spacedim>
783std::vector<std::pair<unsigned int, unsigned int>>
785 const FiniteElement<dim, spacedim> &fe_other) const
786{
787 // we can presently only compute these identities if both FEs are FE_Qs or
788 // if the other one is an FE_Nothing
789 if (const FE_Q_Base<dim, spacedim> *fe_q_other =
790 dynamic_cast<const FE_Q_Base<dim, spacedim> *>(&fe_other))
791 {
792 // dofs are located along lines, so two dofs are identical if they are
793 // located at identical positions. if we had only equidistant points, we
794 // could simply check for similarity like (i+1)*q == (j+1)*p, but we
795 // might have other support points (e.g. Gauss-Lobatto
796 // points). Therefore, read the points in unit_support_points for the
797 // first coordinate direction. We take the lexicographic ordering of the
798 // points in the first direction (i.e., x-direction), which we access
799 // between index 1 and p-1 (index 0 and p are vertex dofs).
800 const unsigned int p = this->degree;
801 const unsigned int q = fe_q_other->degree;
802
803 std::vector<std::pair<unsigned int, unsigned int>> identities;
804
805 const std::vector<unsigned int> &index_map_inverse =
807 const std::vector<unsigned int> &index_map_inverse_other =
808 fe_q_other->get_poly_space_numbering_inverse();
809
810 for (unsigned int i = 0; i < p - 1; ++i)
811 for (unsigned int j = 0; j < q - 1; ++j)
812 if (std::fabs(
813 this->unit_support_points[index_map_inverse[i + 1]][0] -
814 fe_q_other->unit_support_points[index_map_inverse_other[j + 1]]
815 [0]) < 1e-14)
816 identities.emplace_back(i, j);
817
818 return identities;
819 }
820 else if (const FE_SimplexP<dim, spacedim> *fe_p_other =
821 dynamic_cast<const FE_SimplexP<dim, spacedim> *>(&fe_other))
822 {
823 // DoFs are located along lines, so two dofs are identical if they are
824 // located at identical positions. If we had only equidistant points, we
825 // could simply check for similarity like (i+1)*q == (j+1)*p, but we
826 // might have other support points (e.g. Gauss-Lobatto
827 // points). Therefore, read the points in unit_support_points for the
828 // first coordinate direction. For FE_Q, we take the lexicographic
829 // ordering of the line support points in the first direction (i.e.,
830 // x-direction), which we access between index 1 and p-1 (index 0 and p
831 // are vertex dofs). For FE_SimplexP, they are currently hard-coded and we
832 // iterate over points on the first line which begin after the 3 vertex
833 // points in the complete list of unit support points
834
835 Assert(fe_p_other->degree <= 2, ExcNotImplemented());
836
837 const std::vector<unsigned int> &index_map_inverse_q =
839
840 std::vector<std::pair<unsigned int, unsigned int>> identities;
841
842 for (unsigned int i = 0; i < this->degree - 1; ++i)
843 for (unsigned int j = 0; j < fe_p_other->degree - 1; ++j)
844 if (std::fabs(
845 this->unit_support_points[index_map_inverse_q[i + 1]][0] -
846 fe_p_other->get_unit_support_points()[j + 3][0]) < 1e-14)
847 identities.emplace_back(i, j);
848
849 return identities;
850 }
851 else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
852 {
853 // the FE_Nothing has no degrees of freedom, so there are no
854 // equivalencies to be recorded
855 return {};
856 }
857 else if (fe_other.n_unique_faces() == 1 && fe_other.n_dofs_per_face(0) == 0)
858 {
859 // if the other element has no elements on faces at all,
860 // then it would be impossible to enforce any kind of
861 // continuity even if we knew exactly what kind of element
862 // we have -- simply because the other element declares
863 // that it is discontinuous because it has no DoFs on
864 // its faces. in that case, just state that we have no
865 // constraints to declare
866 return {};
867 }
868 else
869 {
870 Assert(false, ExcNotImplemented());
871 return {};
872 }
873}
874
875
876
877template <int dim, int spacedim>
878std::vector<std::pair<unsigned int, unsigned int>>
880 const FiniteElement<dim, spacedim> &fe_other,
881 const unsigned int) const
882{
883 // we can presently only compute these identities if both FEs are FE_Qs or
884 // if the other one is an FE_Nothing
885 if (const FE_Q_Base<dim, spacedim> *fe_q_other =
886 dynamic_cast<const FE_Q_Base<dim, spacedim> *>(&fe_other))
887 {
888 // this works exactly like the line case above, except that now we have
889 // to have two indices i1, i2 and j1, j2 to characterize the dofs on the
890 // face of each of the finite elements. since they are ordered
891 // lexicographically along the first line and we have a tensor product,
892 // the rest is rather straightforward
893 const unsigned int p = this->degree;
894 const unsigned int q = fe_q_other->degree;
895
896 std::vector<std::pair<unsigned int, unsigned int>> identities;
897
898 const std::vector<unsigned int> &index_map_inverse =
900 const std::vector<unsigned int> &index_map_inverse_other =
901 fe_q_other->get_poly_space_numbering_inverse();
902
903 for (unsigned int i1 = 0; i1 < p - 1; ++i1)
904 for (unsigned int i2 = 0; i2 < p - 1; ++i2)
905 for (unsigned int j1 = 0; j1 < q - 1; ++j1)
906 for (unsigned int j2 = 0; j2 < q - 1; ++j2)
907 if ((std::fabs(
908 this->unit_support_points[index_map_inverse[i1 + 1]][0] -
909 fe_q_other
910 ->unit_support_points[index_map_inverse_other[j1 + 1]]
911 [0]) < 1e-14) &&
912 (std::fabs(
913 this->unit_support_points[index_map_inverse[i2 + 1]][0] -
914 fe_q_other
915 ->unit_support_points[index_map_inverse_other[j2 + 1]]
916 [0]) < 1e-14))
917 identities.emplace_back(i1 * (p - 1) + i2, j1 * (q - 1) + j2);
918
919 return identities;
920 }
921 else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
922 {
923 // the FE_Nothing has no degrees of freedom, so there are no
924 // equivalencies to be recorded
925 return std::vector<std::pair<unsigned int, unsigned int>>();
926 }
927 else if (fe_other.n_unique_faces() == 1 && fe_other.n_dofs_per_face(0) == 0)
928 {
929 // if the other element has no elements on faces at all,
930 // then it would be impossible to enforce any kind of
931 // continuity even if we knew exactly what kind of element
932 // we have -- simply because the other element declares
933 // that it is discontinuous because it has no DoFs on
934 // its faces. in that case, just state that we have no
935 // constraints to declare
936 return std::vector<std::pair<unsigned int, unsigned int>>();
937 }
938 else
939 {
940 Assert(false, ExcNotImplemented());
941 return std::vector<std::pair<unsigned int, unsigned int>>();
942 }
943}
944
945
946
947//---------------------------------------------------------------------------
948// Auxiliary functions
949//---------------------------------------------------------------------------
950
951
952
953template <int dim, int spacedim>
954void
956 const std::vector<Point<1>> &points)
957{
958 const std::vector<unsigned int> &index_map_inverse =
960
961 // We can compute the support points by computing the tensor
962 // product of the 1d set of points. We could do this by hand, but it's
963 // easier to just re-use functionality that's already been implemented
964 // for quadrature formulas.
965 const Quadrature<1> support_1d(points);
966 const Quadrature<dim> support_quadrature(support_1d); // NOLINT
967
968 // The only thing we have to do is reorder the points from tensor
969 // product order to the order in which we enumerate DoFs on cells
970 this->unit_support_points.resize(support_quadrature.size());
971 for (unsigned int k = 0; k < support_quadrature.size(); ++k)
972 this->unit_support_points[index_map_inverse[k]] =
973 support_quadrature.point(k);
974}
975
976
977
978template <int dim, int spacedim>
979void
981 const std::vector<Point<1>> &points)
982{
983 // TODO: the implementation makes the assumption that all faces have the
984 // same number of dofs
986 const unsigned int face_no = 0;
987
988 this->unit_face_support_points[face_no].resize(
989 Utilities::fixed_power<dim - 1>(q_degree + 1));
990
991 // In 1d, there is only one 0-dimensional support point, so there is nothing
992 // more to be done.
993 if (dim == 1)
994 return;
995
996 // find renumbering of faces and assign from values of quadrature
997 const std::vector<unsigned int> face_index_map =
999
1000 // We can compute the support points by computing the tensor
1001 // product of the 1d set of points. We could do this by hand, but it's
1002 // easier to just re-use functionality that's already been implemented
1003 // for quadrature formulas.
1004 const Quadrature<1> support_1d(points);
1005 const Quadrature<dim - 1> support_quadrature(support_1d); // NOLINT
1006
1007 // The only thing we have to do is reorder the points from tensor
1008 // product order to the order in which we enumerate DoFs on cells
1009 this->unit_face_support_points[face_no].resize(support_quadrature.size());
1010 for (unsigned int k = 0; k < support_quadrature.size(); ++k)
1011 this->unit_face_support_points[face_no][face_index_map[k]] =
1012 support_quadrature.point(k);
1013}
1014
1015
1016
1017template <int dim, int spacedim>
1018void
1020{
1021 // for 1d and 2d, do nothing
1022 if (dim < 3)
1023 return;
1024
1025 // TODO: the implementation makes the assumption that all faces have the
1026 // same number of dofs
1027 AssertDimension(this->n_unique_faces(), 1);
1028 const unsigned int face_no = 0;
1029
1030 Assert(
1032 this->reference_cell().n_face_orientations(face_no) *
1033 this->n_dofs_per_quad(face_no),
1035
1036 const unsigned int n = q_degree - 1;
1037 Assert(n * n == this->n_dofs_per_quad(face_no), ExcInternalError());
1038
1039 // the dofs on a face are connected to a n x n matrix. for example, for
1040 // degree==4 we have the following dofs on a quad
1041
1042 // ___________
1043 // | |
1044 // | 6 7 8 |
1045 // | |
1046 // | 3 4 5 |
1047 // | |
1048 // | 0 1 2 |
1049 // |___________|
1050 //
1051 // we have dof_no=i+n*j with index i in x-direction and index j in
1052 // y-direction running from 0 to n-1. to extract i and j we can use
1053 // i=dof_no%n and j=dof_no/n. i and j can then be used to construct the
1054 // rotated and mirrored numbers.
1055
1056
1057 for (unsigned int local = 0; local < this->n_dofs_per_quad(face_no); ++local)
1058 // face support points are in lexicographic ordering with x running
1059 // fastest. invert that (y running fastest)
1060 {
1061 unsigned int i = local % n, j = local / n;
1062
1063 // face_orientation=false, face_flip=false, face_rotation=false
1065 0) =
1066 j + i * n - local;
1067 // face_orientation=false, face_flip=false, face_rotation=true
1069 1) =
1070 i + (n - 1 - j) * n - local;
1071 // face_orientation=false, face_flip=true, face_rotation=false
1073 2) =
1074 (n - 1 - j) + (n - 1 - i) * n - local;
1075 // face_orientation=false, face_flip=true, face_rotation=true
1077 3) =
1078 (n - 1 - i) + j * n - local;
1079 // face_orientation=true, face_flip=false, face_rotation=false
1081 4) = 0;
1082 // face_orientation=true, face_flip=false, face_rotation=true
1084 5) =
1085 j + (n - 1 - i) * n - local;
1086 // face_orientation=true, face_flip=true, face_rotation=false
1088 6) =
1089 (n - 1 - i) + (n - 1 - j) * n - local;
1090 // face_orientation=true, face_flip=true, face_rotation=true
1092 7) =
1093 (n - 1 - j) + i * n - local;
1094 }
1095
1096 // additionally initialize reordering of line dofs
1097 for (unsigned int i = 0; i < this->n_dofs_per_line(); ++i)
1099 this->n_dofs_per_line() - 1 - i - i;
1100}
1101
1102
1103
1104template <int dim, int spacedim>
1105unsigned int
1106FE_Q_Base<dim, spacedim>::face_to_cell_index(const unsigned int face_index,
1107 const unsigned int face,
1108 const bool face_orientation,
1109 const bool face_flip,
1110 const bool face_rotation) const
1111{
1112 AssertIndexRange(face_index, this->n_dofs_per_face(face));
1114
1115 // TODO: we could presumably solve the 3d case below using the
1116 // adjust_quad_dof_index_for_face_orientation_table field. for the
1117 // 2d case, we can't use adjust_line_dof_index_for_line_orientation_table
1118 // since that array is empty (presumably because we thought that
1119 // there are no flipped edges in 2d, but these can happen in
1120 // DoFTools::make_periodicity_constraints, for example). so we
1121 // would need to either fill this field, or rely on derived classes
1122 // implementing this function, as we currently do
1123
1124 // we need to distinguish between DoFs on vertices, lines and in 3d quads.
1125 // do so in a sequence of if-else statements
1126 if (face_index < this->get_first_face_line_index(face))
1127 // DoF is on a vertex
1128 {
1129 // get the number of the vertex on the face that corresponds to this DoF,
1130 // along with the number of the DoF on this vertex
1131 const unsigned int face_vertex = face_index / this->n_dofs_per_vertex();
1132 const unsigned int dof_index_on_vertex =
1133 face_index % this->n_dofs_per_vertex();
1134
1135 // then get the number of this vertex on the cell and translate
1136 // this to a DoF number on the cell
1138 face, face_vertex, face_orientation, face_flip, face_rotation) *
1139 this->n_dofs_per_vertex() +
1140 dof_index_on_vertex);
1141 }
1142 else if (face_index < this->get_first_face_quad_index(face))
1143 // DoF is on a face
1144 {
1145 // do the same kind of translation as before. we need to only consider
1146 // DoFs on the lines, i.e., ignoring those on the vertices
1147 const unsigned int index =
1148 face_index - this->get_first_face_line_index(face);
1149
1150 const unsigned int face_line = index / this->n_dofs_per_line();
1151 const unsigned int dof_index_on_line = index % this->n_dofs_per_line();
1152
1153 // we now also need to adjust the line index for the case of
1154 // face orientation, face flips, etc
1155 unsigned int adjusted_dof_index_on_line = 0;
1156 switch (dim)
1157 {
1158 case 1:
1159 Assert(false, ExcInternalError());
1160 break;
1161
1162 case 2:
1163 // in 2d, only face_flip has a meaning. if it is set, consider
1164 // dofs in reverse order
1165 if (face_flip == false)
1166 adjusted_dof_index_on_line = dof_index_on_line;
1167 else
1168 adjusted_dof_index_on_line =
1169 this->n_dofs_per_line() - dof_index_on_line - 1;
1170 break;
1171
1172 case 3:
1173 // in 3d, things are difficult. someone will have to think
1174 // about how this code here should look like, by drawing a bunch
1175 // of pictures of how all the faces can look like with the various
1176 // flips and rotations.
1177 //
1178 // that said, the Q2 case is easy enough to implement, as is the
1179 // case where everything is in standard orientation
1180 Assert((this->n_dofs_per_line() <= 1) ||
1181 ((face_orientation == true) && (face_flip == false) &&
1182 (face_rotation == false)),
1184 adjusted_dof_index_on_line = dof_index_on_line;
1185 break;
1186
1187 default:
1188 Assert(false, ExcInternalError());
1189 }
1190
1191 return (this->get_first_line_index() +
1193 face, face_line, face_orientation, face_flip, face_rotation) *
1194 this->n_dofs_per_line() +
1195 adjusted_dof_index_on_line);
1196 }
1197 else
1198 // DoF is on a quad
1199 {
1200 Assert(dim >= 3, ExcInternalError());
1201
1202 // ignore vertex and line dofs
1203 const unsigned int index =
1204 face_index - this->get_first_face_quad_index(face);
1205
1206 // the same is true here as above for the 3d case -- someone will
1207 // just have to draw a bunch of pictures. in the meantime,
1208 // we can implement the Q2 case in which it is simple
1209 Assert((this->n_dofs_per_quad(face) <= 1) ||
1210 ((face_orientation == true) && (face_flip == false) &&
1211 (face_rotation == false)),
1213 return (this->get_first_quad_index(face) + index);
1214 }
1215}
1216
1217
1218
1219template <int dim, int spacedim>
1220std::vector<unsigned int>
1222{
1223 using FEQ = FE_Q_Base<dim, spacedim>;
1224 AssertThrow(degree > 0, typename FEQ::ExcFEQCannotHaveDegree0());
1225 std::vector<unsigned int> dpo(dim + 1, 1U);
1226 for (unsigned int i = 1; i < dpo.size(); ++i)
1227 dpo[i] = dpo[i - 1] * (degree - 1);
1228 return dpo;
1229}
1230
1231
1232
1233template <int dim, int spacedim>
1234void
1236 const std::vector<Point<1>> &points)
1237{
1239}
1240
1241
1242
1243template <int dim, int spacedim>
1244const FullMatrix<double> &
1246 const unsigned int child,
1247 const RefinementCase<dim> &refinement_case) const
1248{
1249 AssertIndexRange(refinement_case,
1251 Assert(refinement_case != RefinementCase<dim>::no_refinement,
1252 ExcMessage(
1253 "Prolongation matrices are only available for refined cells!"));
1254 AssertIndexRange(child, GeometryInfo<dim>::n_children(refinement_case));
1255
1256 // initialization upon first request
1257 if (this->prolongation[refinement_case - 1][child].n() == 0)
1258 {
1259 std::lock_guard<std::mutex> lock(this->mutex);
1260
1261 // if matrix got updated while waiting for the lock
1262 if (this->prolongation[refinement_case - 1][child].n() ==
1263 this->n_dofs_per_cell())
1264 return this->prolongation[refinement_case - 1][child];
1265
1266 // distinguish q/q_dg0 case: only treat Q dofs first
1267 const unsigned int q_dofs_per_cell =
1268 Utilities::fixed_power<dim>(q_degree + 1);
1269
1270 // compute the interpolation matrices in much the same way as we do for
1271 // the constraints. it's actually simpler here, since we don't have this
1272 // weird renumbering stuff going on. The trick is again that we the
1273 // interpolation matrix is formed by a permutation of the indices of the
1274 // cell matrix. The value eps is used a threshold to decide when certain
1275 // evaluations of the Lagrange polynomials are zero or one.
1276 const double eps = 1e-15 * q_degree * dim;
1277
1278#ifdef DEBUG
1279 // in DEBUG mode, check that the evaluation of support points in the
1280 // current numbering gives the identity operation
1281 for (unsigned int i = 0; i < q_dofs_per_cell; ++i)
1282 {
1283 Assert(std::fabs(1. - this->poly_space->compute_value(
1284 i, this->unit_support_points[i])) < eps,
1285 ExcInternalError("The Lagrange polynomial does not evaluate "
1286 "to one or zero in a nodal point. "
1287 "This typically indicates that the "
1288 "polynomial interpolation is "
1289 "ill-conditioned such that round-off "
1290 "prevents the sum to be one."));
1291 for (unsigned int j = 0; j < q_dofs_per_cell; ++j)
1292 if (j != i)
1293 Assert(std::fabs(this->poly_space->compute_value(
1294 i, this->unit_support_points[j])) < eps,
1296 "The Lagrange polynomial does not evaluate "
1297 "to one or zero in a nodal point. "
1298 "This typically indicates that the "
1299 "polynomial interpolation is "
1300 "ill-conditioned such that round-off "
1301 "prevents the sum to be one."));
1302 }
1303#endif
1304
1305 // to efficiently evaluate the polynomial at the subcell, make use of
1306 // the tensor product structure of this element and only evaluate 1d
1307 // information from the polynomial. This makes the cost of this function
1308 // almost negligible also for high order elements
1309 const unsigned int dofs1d = q_degree + 1;
1310 std::vector<Table<2, double>> subcell_evaluations(
1311 dim, Table<2, double>(dofs1d, dofs1d));
1312
1313 const std::vector<unsigned int> &index_map_inverse =
1315
1316 // helper value: step size how to walk through diagonal and how many
1317 // points we have left apart from the first dimension
1318 unsigned int step_size_diag = 0;
1319 {
1320 unsigned int factor = 1;
1321 for (unsigned int d = 0; d < dim; ++d)
1322 {
1323 step_size_diag += factor;
1324 factor *= dofs1d;
1325 }
1326 }
1327
1328 FullMatrix<double> prolongate(this->n_dofs_per_cell(),
1329 this->n_dofs_per_cell());
1330
1331 // go through the points in diagonal to capture variation in all
1332 // directions simultaneously
1333 for (unsigned int j = 0; j < dofs1d; ++j)
1334 {
1335 const unsigned int diag_comp = index_map_inverse[j * step_size_diag];
1336 const Point<dim> p_subcell = this->unit_support_points[diag_comp];
1337 const Point<dim> p_cell =
1339 child,
1340 refinement_case);
1341 for (unsigned int i = 0; i < dofs1d; ++i)
1342 for (unsigned int d = 0; d < dim; ++d)
1343 {
1344 // evaluate along line where only x is different from zero
1346 point[0] = p_cell[d];
1347 const double cell_value =
1348 this->poly_space->compute_value(index_map_inverse[i], point);
1349
1350 // cut off values that are too small. note that we have here
1351 // Lagrange interpolation functions, so they should be zero at
1352 // almost all points, and one at the others, at least on the
1353 // subcells. so set them to their exact values
1354 //
1355 // the actual cut-off value is somewhat fuzzy, but it works
1356 // for 2e-13*degree*dim (see above), which is kind of
1357 // reasonable given that we compute the values of the
1358 // polynomials via an degree-step recursion and then multiply
1359 // the 1d-values. this gives us a linear growth in degree*dim,
1360 // times a small constant.
1361 //
1362 // the embedding matrix is given by applying the inverse of
1363 // the subcell matrix on the cell_interpolation matrix. since
1364 // the subcell matrix is actually only a permutation vector,
1365 // all we need to do is to switch the rows we write the data
1366 // into. moreover, cut off very small values here
1367 if (std::fabs(cell_value) < eps)
1368 subcell_evaluations[d](j, i) = 0;
1369 else
1370 subcell_evaluations[d](j, i) = cell_value;
1371 }
1372 }
1373
1374 // now expand from 1d info. block innermost dimension (x_0) in order to
1375 // avoid difficult checks at innermost loop
1376 unsigned int j_indices[dim];
1377 internal::FE_Q_Base::zero_indices<dim>(j_indices);
1378 for (unsigned int j = 0; j < q_dofs_per_cell; j += dofs1d)
1379 {
1380 unsigned int i_indices[dim];
1381 internal::FE_Q_Base::zero_indices<dim>(i_indices);
1382 for (unsigned int i = 0; i < q_dofs_per_cell; i += dofs1d)
1383 {
1384 double val_extra_dim = 1.;
1385 for (unsigned int d = 1; d < dim; ++d)
1386 val_extra_dim *=
1387 subcell_evaluations[d](j_indices[d - 1], i_indices[d - 1]);
1388
1389 // innermost sum where we actually compute. the same as
1390 // prolongate(j,i) = this->poly_space->compute_value (i, p_cell)
1391 for (unsigned int jj = 0; jj < dofs1d; ++jj)
1392 {
1393 const unsigned int j_ind = index_map_inverse[j + jj];
1394 for (unsigned int ii = 0; ii < dofs1d; ++ii)
1395 prolongate(j_ind, index_map_inverse[i + ii]) =
1396 val_extra_dim * subcell_evaluations[0](jj, ii);
1397 }
1398
1399 // update indices that denote the tensor product position. a bit
1400 // fuzzy and therefore not done for innermost x_0 direction
1401 internal::FE_Q_Base::increment_indices<dim>(i_indices, dofs1d);
1402 }
1403 Assert(i_indices[dim - 1] == 1, ExcInternalError());
1404 internal::FE_Q_Base::increment_indices<dim>(j_indices, dofs1d);
1405 }
1406
1407 // the discontinuous node is simply mapped on the discontinuous node on
1408 // the child element
1409 if (q_dofs_per_cell < this->n_dofs_per_cell())
1410 prolongate(q_dofs_per_cell, q_dofs_per_cell) = 1.;
1411
1412 // and make sure that the row sum is 1. this must be so since for this
1413 // element, the shape functions add up to one
1414#ifdef DEBUG
1415 for (unsigned int row = 0; row < this->n_dofs_per_cell(); ++row)
1416 {
1417 double sum = 0;
1418 for (unsigned int col = 0; col < this->n_dofs_per_cell(); ++col)
1419 sum += prolongate(row, col);
1420 Assert(std::fabs(sum - 1.) <
1421 std::max(eps, 5e-16 * std::sqrt(this->n_dofs_per_cell())),
1422 ExcInternalError("The entries in a row of the local "
1423 "prolongation matrix do not add to one. "
1424 "This typically indicates that the "
1425 "polynomial interpolation is "
1426 "ill-conditioned such that round-off "
1427 "prevents the sum to be one."));
1428 }
1429#endif
1430
1431 // swap matrices
1432 prolongate.swap(const_cast<FullMatrix<double> &>(
1433 this->prolongation[refinement_case - 1][child]));
1434 }
1435
1436 // finally return the matrix
1437 return this->prolongation[refinement_case - 1][child];
1438}
1439
1440
1441
1442template <int dim, int spacedim>
1443const FullMatrix<double> &
1445 const unsigned int child,
1446 const RefinementCase<dim> &refinement_case) const
1447{
1448 AssertIndexRange(refinement_case,
1450 Assert(refinement_case != RefinementCase<dim>::no_refinement,
1451 ExcMessage(
1452 "Restriction matrices are only available for refined cells!"));
1453 AssertIndexRange(child, GeometryInfo<dim>::n_children(refinement_case));
1454
1455 // initialization upon first request
1456 if (this->restriction[refinement_case - 1][child].n() == 0)
1457 {
1458 std::lock_guard<std::mutex> lock(this->mutex);
1459
1460 // if matrix got updated while waiting for the lock...
1461 if (this->restriction[refinement_case - 1][child].n() ==
1462 this->n_dofs_per_cell())
1463 return this->restriction[refinement_case - 1][child];
1464
1465 FullMatrix<double> my_restriction(this->n_dofs_per_cell(),
1466 this->n_dofs_per_cell());
1467 // distinguish q/q_dg0 case
1468 const unsigned int q_dofs_per_cell =
1469 Utilities::fixed_power<dim>(q_degree + 1);
1470
1471 // for Lagrange interpolation polynomials based on equidistant points,
1472 // construction of the restriction matrices is relatively simple. the
1473 // reason is that in this case the interpolation points on the mother
1474 // cell are always also interpolation points for some shape function on
1475 // one or the other child.
1476 //
1477 // in the general case with non-equidistant points, we need to actually
1478 // do an interpolation. thus, we take the interpolation points on the
1479 // mother cell and evaluate the shape functions of the child cell on
1480 // those points. it does not hurt in the equidistant case because then
1481 // simple one shape function evaluates to one and the others to zero.
1482 //
1483 // this element is non-additive in all its degrees of freedom by
1484 // default, which requires care in downstream use. fortunately, even the
1485 // interpolation on non-equidistant points is invariant under the
1486 // assumption that whenever a row makes a non-zero contribution to the
1487 // mother's residual, the correct value is interpolated.
1488
1489 const double eps = 1e-15 * q_degree * dim;
1490 const std::vector<unsigned int> &index_map_inverse =
1492
1493 const unsigned int dofs1d = q_degree + 1;
1494 std::vector<Tensor<1, dim>> evaluations1d(dofs1d);
1495
1496 my_restriction.reinit(this->n_dofs_per_cell(), this->n_dofs_per_cell());
1497
1498 for (unsigned int i = 0; i < q_dofs_per_cell; ++i)
1499 {
1500 unsigned int mother_dof = index_map_inverse[i];
1501 const Point<dim> p_cell = this->unit_support_points[mother_dof];
1502
1503 // check whether this interpolation point is inside this child cell
1504 const Point<dim> p_subcell =
1506 child,
1507 refinement_case);
1509 {
1510 // same logic as in initialize_embedding to evaluate the
1511 // polynomial faster than from the tensor product: since we
1512 // evaluate all polynomials, it is much faster to just compute
1513 // the 1d values for all polynomials before and then get the
1514 // dim-data.
1515 for (unsigned int j = 0; j < dofs1d; ++j)
1516 for (unsigned int d = 0; d < dim; ++d)
1517 {
1519 point[0] = p_subcell[d];
1520 evaluations1d[j][d] =
1521 this->poly_space->compute_value(index_map_inverse[j],
1522 point);
1523 }
1524 unsigned int j_indices[dim];
1525 internal::FE_Q_Base::zero_indices<dim>(j_indices);
1526#ifdef DEBUG
1527 double sum_check = 0;
1528#endif
1529 for (unsigned int j = 0; j < q_dofs_per_cell; j += dofs1d)
1530 {
1531 double val_extra_dim = 1.;
1532 for (unsigned int d = 1; d < dim; ++d)
1533 val_extra_dim *= evaluations1d[j_indices[d - 1]][d];
1534 for (unsigned int jj = 0; jj < dofs1d; ++jj)
1535 {
1536 // find the child shape function(s) corresponding to
1537 // this point. Usually this is just one function;
1538 // however, when we use FE_Q on arbitrary nodes a parent
1539 // support point might not be a child support point, and
1540 // then we will get more than one nonzero value per
1541 // row. Still, the values should sum up to 1
1542 const double val = val_extra_dim * evaluations1d[jj][0];
1543 const unsigned int child_dof = index_map_inverse[j + jj];
1544 if (std::fabs(val - 1.) < eps)
1545 my_restriction(mother_dof, child_dof) = 1.;
1546 else if (std::fabs(val) > eps)
1547 my_restriction(mother_dof, child_dof) = val;
1548#ifdef DEBUG
1549 sum_check += val;
1550#endif
1551 }
1552 internal::FE_Q_Base::increment_indices<dim>(j_indices,
1553 dofs1d);
1554 }
1555 Assert(std::fabs(sum_check - 1) <
1556 std::max(eps,
1557 5e-16 * std::sqrt(this->n_dofs_per_cell())),
1558 ExcInternalError("The entries in a row of the local "
1559 "restriction matrix do not add to one. "
1560 "This typically indicates that the "
1561 "polynomial interpolation is "
1562 "ill-conditioned such that round-off "
1563 "prevents the sum to be one."));
1564 }
1565
1566 // part for FE_Q_DG0
1567 if (q_dofs_per_cell < this->n_dofs_per_cell())
1568 my_restriction(this->n_dofs_per_cell() - 1,
1569 this->n_dofs_per_cell() - 1) =
1571 RefinementCase<dim>(refinement_case));
1572 }
1573
1574 // swap the just computed restriction matrix into the
1575 // element of the vector stored in the base class
1576 my_restriction.swap(const_cast<FullMatrix<double> &>(
1577 this->restriction[refinement_case - 1][child]));
1578 }
1579
1580 return this->restriction[refinement_case - 1][child];
1581}
1582
1583
1584
1585//---------------------------------------------------------------------------
1586// Data field initialization
1587//---------------------------------------------------------------------------
1588
1589
1590template <int dim, int spacedim>
1591bool
1593 const unsigned int shape_index,
1594 const unsigned int face_index) const
1595{
1596 AssertIndexRange(shape_index, this->n_dofs_per_cell());
1598
1599 // in 1d, things are simple. since there is only one degree of freedom per
1600 // vertex in this class, the first is on vertex 0 (==face 0 in some sense),
1601 // the second on face 1:
1602 if (dim == 1)
1603 return (((shape_index == 0) && (face_index == 0)) ||
1604 ((shape_index == 1) && (face_index == 1)));
1605
1606 // first, special-case interior shape functions, since they have no support
1607 // no-where on the boundary
1608 if (((dim == 2) &&
1609 (shape_index >= this->get_first_quad_index(0 /*first quad*/))) ||
1610 ((dim == 3) && (shape_index >= this->get_first_hex_index())))
1611 return false;
1612
1613 // let's see whether this is a vertex
1614 if (shape_index < this->get_first_line_index())
1615 {
1616 // for Q elements, there is one dof per vertex, so
1617 // shape_index==vertex_number. check whether this vertex is on the given
1618 // face. thus, for each face, give a list of vertices
1619 const unsigned int vertex_no = shape_index;
1622
1623 for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face; ++v)
1624 if (GeometryInfo<dim>::face_to_cell_vertices(face_index, v) ==
1625 vertex_no)
1626 return true;
1627
1628 return false;
1629 }
1630 else if (shape_index < this->get_first_quad_index(0 /*first quad*/))
1631 // ok, dof is on a line
1632 {
1633 const unsigned int line_index =
1634 (shape_index - this->get_first_line_index()) / this->n_dofs_per_line();
1637
1638 // in 2d, the line is the face, so get the line index
1639 if (dim == 2)
1640 return (line_index == face_index);
1641 else if (dim == 3)
1642 {
1643 // silence compiler warning
1644 const unsigned int lines_per_face =
1645 dim == 3 ? GeometryInfo<dim>::lines_per_face : 1;
1646 // see whether the given line is on the given face.
1647 for (unsigned int l = 0; l < lines_per_face; ++l)
1648 if (GeometryInfo<3>::face_to_cell_lines(face_index, l) ==
1649 line_index)
1650 return true;
1651
1652 return false;
1653 }
1654 else
1655 Assert(false, ExcNotImplemented());
1656 }
1657 else if (shape_index < this->get_first_hex_index())
1658 // dof is on a quad
1659 {
1660 const unsigned int quad_index =
1661 (shape_index - this->get_first_quad_index(0)) /
1662 this->n_dofs_per_quad(face_index); // this won't work
1663 Assert(static_cast<signed int>(quad_index) <
1664 static_cast<signed int>(GeometryInfo<dim>::quads_per_cell),
1666
1667 // in 2d, cell bubble are zero on all faces. but we have treated this
1668 // case above already
1669 Assert(dim != 2, ExcInternalError());
1670
1671 // in 3d, quad_index=face_index
1672 if (dim == 3)
1673 return (quad_index == face_index);
1674 else
1675 Assert(false, ExcNotImplemented());
1676 }
1677 else
1678 // dof on hex
1679 {
1680 // can only happen in 3d, but this case has already been covered above
1681 Assert(false, ExcNotImplemented());
1682 return false;
1683 }
1684
1685 // we should not have gotten here
1686 Assert(false, ExcInternalError());
1687 return false;
1688}
1689
1690
1691
1692template <int dim, int spacedim>
1693std::pair<Table<2, bool>, std::vector<unsigned int>>
1695{
1696 Table<2, bool> constant_modes(1, this->n_dofs_per_cell());
1697 // We here just care for the constant mode due to the polynomial space
1698 // without any enrichments
1699 // As there may be more constant modes derived classes may to implement this
1700 // themselves. An example for this is FE_Q_DG0.
1701 for (unsigned int i = 0; i < Utilities::fixed_power<dim>(q_degree + 1); ++i)
1702 constant_modes(0, i) = true;
1703 return std::pair<Table<2, bool>, std::vector<unsigned int>>(
1704 constant_modes, std::vector<unsigned int>(1, 0));
1705}
1706
1707
1708
1709// explicit instantiations
1710#include "fe_q_base.inst"
1711
std::vector< unsigned int > get_poly_space_numbering_inverse() const
const ScalarPolynomialsBase< dim > & get_poly_space() 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:533
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
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:346
virtual std::pair< Table< 2, bool >, std::vector< unsigned int > > get_constant_modes() const override
Threads::Mutex mutex
Definition fe_q_base.h:339
void initialize_quad_dof_index_permutation()
virtual unsigned int face_to_cell_index(const unsigned int face_dof_index, const unsigned int face, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false) 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)
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 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:453
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
unsigned int n_components() 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:2447
bool has_face_support_points(const unsigned int face_no=0) const
std::vector< std::vector< FullMatrix< double > > > restriction
Definition fe.h:2402
std::vector< Table< 2, int > > adjust_quad_dof_index_for_face_orientation_table
Definition fe.h:2476
const std::vector< Point< dim - 1 > > & get_unit_face_support_points(const unsigned int face_no=0) const
std::vector< int > adjust_line_dof_index_for_line_orientation_table
Definition fe.h:2491
TableIndices< 2 > interface_constraints_size() const
std::vector< Point< dim > > unit_support_points
Definition fe.h:2440
FullMatrix< double > interface_constraints
Definition fe.h:2428
std::vector< std::vector< FullMatrix< double > > > prolongation
Definition fe.h:2416
size_type n() const
size_type m() const
Definition point.h:112
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)
const Point< dim > & point(const unsigned int i) const
void set_numbering(const std::vector< unsigned int > &renumber)
void set_numbering(const std::vector< unsigned int > &renumber)
void set_numbering(const std::vector< unsigned int > &renumber)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#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)
std::vector< unsigned int > lexicographic_to_hierarchic_numbering(unsigned int degree)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:189
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)
static const unsigned int invalid_unsigned_int
Definition types.h:213
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:217
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 unsigned int n_children(const RefinementCase< dim > &refinement_case)
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)