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