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