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