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