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