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