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