Reference documentation for deal.II version GIT d9d8a449a2 2022-08-17 08:45: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_hierarchical.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2002 - 2022 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 
17 #include <deal.II/fe/fe_dgq.h>
18 #include <deal.II/fe/fe_nothing.h>
20 
21 #include <cmath>
22 #include <memory>
23 #include <sstream>
24 
25 // TODO: implement the adjust_quad_dof_index_for_face_orientation_table and
26 // adjust_line_dof_index_for_line_orientation_table fields, and write tests
27 // similar to bits/face_orientation_and_fe_q_*
28 
29 
31 
32 namespace internal
33 {
35  {
36  namespace
37  {
41  inline std::vector<unsigned int>
42  invert_numbering(const std::vector<unsigned int> &in)
43  {
44  std::vector<unsigned int> out(in.size());
45  for (unsigned int i = 0; i < in.size(); ++i)
46  {
47  AssertIndexRange(in[i], out.size());
48  out[in[i]] = i;
49  }
50  return out;
51  }
52  } // namespace
53  } // namespace FE_Q_Hierarchical
54 } // namespace internal
55 
56 
57 
58 template <int dim>
61  Polynomials::Hierarchical::generate_complete_basis(degree)),
62  FiniteElementData<dim>(get_dpo_vector(degree),
63  1,
64  degree,
65  FiniteElementData<dim>::H1),
66  std::vector<bool>(
67  FiniteElementData<dim>(get_dpo_vector(degree), 1, degree)
68  .n_dofs_per_cell(),
69  false),
70  std::vector<ComponentMask>(
71  FiniteElementData<dim>(get_dpo_vector(degree), 1, degree)
72  .n_dofs_per_cell(),
73  std::vector<bool>(1, true)))
74  , face_renumber(face_fe_q_hierarchical_to_hierarchic_numbering(degree))
75 {
76  TensorProductPolynomials<dim> *poly_space_derived_ptr =
77  dynamic_cast<TensorProductPolynomials<dim> *>(this->poly_space.get());
78  poly_space_derived_ptr->set_numbering(
80 
81  // The matrix @p{dofs_cell} contains the
82  // values of the linear functionals of
83  // the parent 1d cell applied to the
84  // shape functions of the two 1d subcells.
85  // The matrix @p{dofs_subcell} contains
86  // the values of the linear functionals
87  // on each 1d subcell applied to the
88  // shape functions on the parent 1d
89  // subcell.
90  // We use @p{dofs_cell} and
91  // @p{dofs_subcell} to compute the
92  // @p{prolongation}, @p{restriction} and
93  // @p{interface_constraints} matrices
94  // for all dimensions.
95  std::vector<FullMatrix<double>> dofs_cell(
98  2 * this->n_dofs_per_vertex() +
99  this->n_dofs_per_line()));
100  std::vector<FullMatrix<double>> dofs_subcell(
103  2 * this->n_dofs_per_vertex() +
104  this->n_dofs_per_line()));
105  // build these fields, as they are
106  // needed as auxiliary fields later
107  // on
108  build_dofs_cell(dofs_cell, dofs_subcell);
109 
110  // then use them to initialize
111  // other fields
112  initialize_constraints(dofs_subcell);
113  initialize_embedding_and_restriction(dofs_cell, dofs_subcell);
114 
115  // finally fill in support points
116  // on cell and face
119 }
120 
121 
122 
123 template <int dim>
124 std::string
126 {
127  // note that the
128  // FETools::get_fe_by_name
129  // function depends on the
130  // particular format of the string
131  // this function returns, so they
132  // have to be kept in synch
133 
134  std::ostringstream namebuf;
135  namebuf << "FE_Q_Hierarchical<" << dim << ">(" << this->degree << ")";
136 
137  return namebuf.str();
138 }
139 
140 
141 
142 template <int dim>
143 std::unique_ptr<FiniteElement<dim, dim>>
145 {
146  return std::make_unique<FE_Q_Hierarchical<dim>>(*this);
147 }
148 
149 
150 
151 template <int dim>
152 void
154  const FiniteElement<dim> &source,
155  FullMatrix<double> & matrix) const
156 {
157  // support interpolation between FE_Q_Hierarchical only.
158  if (const FE_Q_Hierarchical<dim> *source_fe =
159  dynamic_cast<const FE_Q_Hierarchical<dim> *>(&source))
160  {
161  // ok, source is a Q_Hierarchical element, so we will be able to do the
162  // work
163  Assert(matrix.m() == this->n_dofs_per_cell(),
164  ExcDimensionMismatch(matrix.m(), this->n_dofs_per_cell()));
165  Assert(matrix.n() == source.n_dofs_per_cell(),
166  ExcDimensionMismatch(matrix.m(), source_fe->n_dofs_per_cell()));
167 
168  // Recall that DoFs are renumbered in the following order:
169  // vertices, lines, quads, hexes.
170  // As we deal with hierarchical FE, interpolation matrix is rather easy:
171  // it has 1 on pairs of dofs which are the same.
172  // To get those use get_embedding_dofs();
173 
174  matrix = 0.;
175 
176  // distinguish between the case when we interpolate to a richer element
177  if (this->n_dofs_per_cell() >= source_fe->n_dofs_per_cell())
178  {
179  const std::vector<unsigned int> dof_map =
180  this->get_embedding_dofs(source_fe->degree);
181  for (unsigned int j = 0; j < dof_map.size(); ++j)
182  matrix[dof_map[j]][j] = 1.;
183  }
184  // and when just truncate higher modes.
185  else
186  {
187  const std::vector<unsigned int> dof_map =
188  source_fe->get_embedding_dofs(this->degree);
189  for (unsigned int j = 0; j < dof_map.size(); ++j)
190  matrix[j][dof_map[j]] = 1.;
191  }
192  }
193  else
194  {
195  AssertThrow(
196  false,
197  ::ExcMessage(
198  "Interpolation is supported only between FE_Q_Hierarchical"));
199  }
200 }
201 
202 template <int dim>
203 const FullMatrix<double> &
205  const unsigned int child,
206  const RefinementCase<dim> &refinement_case) const
207 {
208  Assert(
210  ExcMessage(
211  "Prolongation matrices are only available for isotropic refinement!"));
212 
213  AssertIndexRange(child, GeometryInfo<dim>::n_children(refinement_case));
214 
215  return this->prolongation[refinement_case - 1][child];
216 }
217 
218 
219 template <int dim>
220 bool
222 {
223  return true;
224 }
225 
226 
227 template <int dim>
228 std::vector<std::pair<unsigned int, unsigned int>>
230  const FiniteElement<dim> &fe_other) const
231 {
232  // we can presently only compute
233  // these identities if both FEs are
234  // FE_Q_Hierarchicals or if the other
235  // one is an FE_Nothing. in the first
236  // case, there should be exactly one
237  // single DoF of each FE at a
238  // vertex, and they should have
239  // identical value
240  if (dynamic_cast<const FE_Q_Hierarchical<dim> *>(&fe_other) != nullptr)
241  {
242  return std::vector<std::pair<unsigned int, unsigned int>>(
243  1, std::make_pair(0U, 0U));
244  }
245  else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
246  {
247  // the FE_Nothing has no
248  // degrees of freedom, so there
249  // are no equivalencies to be
250  // recorded
251  return std::vector<std::pair<unsigned int, unsigned int>>();
252  }
253  else
254  {
255  Assert(false, ExcNotImplemented());
256  return std::vector<std::pair<unsigned int, unsigned int>>();
257  }
258 }
259 
260 template <int dim>
261 std::vector<std::pair<unsigned int, unsigned int>>
263  const FiniteElement<dim> &fe_other) const
264 {
265  // we can presently only compute
266  // these identities if both FEs are
267  // FE_Q_Hierarchicals or if the other
268  // one is an FE_Nothing.
269  if (dynamic_cast<const FE_Q_Hierarchical<dim> *>(&fe_other) != nullptr)
270  {
271  const unsigned int this_dpl = this->n_dofs_per_line();
272  const unsigned int other_dpl = fe_other.n_dofs_per_line();
273 
274  // we deal with hierarchical 1d polynomials where dofs are enumerated
275  // increasingly. Thus we return a vector of pairs for the first N-1, where
276  // N is minimum number of dofs_per_line for each FE_Q_Hierarchical.
277  std::vector<std::pair<unsigned int, unsigned int>> res;
278  for (unsigned int i = 0; i < std::min(this_dpl, other_dpl); ++i)
279  res.emplace_back(i, i);
280 
281  return res;
282  }
283  else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
284  {
285  // the FE_Nothing has no
286  // degrees of freedom, so there
287  // are no equivalencies to be
288  // recorded
289  return std::vector<std::pair<unsigned int, unsigned int>>();
290  }
291  else
292  {
293  Assert(false, ExcNotImplemented());
294  return std::vector<std::pair<unsigned int, unsigned int>>();
295  }
296 }
297 
298 template <int dim>
299 std::vector<std::pair<unsigned int, unsigned int>>
301  const FiniteElement<dim> &fe_other,
302  const unsigned int) const
303 {
304  // we can presently only compute
305  // these identities if both FEs are
306  // FE_Q_Hierarchicals or if the other
307  // one is an FE_Nothing.
308  if (dynamic_cast<const FE_Q_Hierarchical<dim> *>(&fe_other) != nullptr)
309  {
310  // TODO: the implementation makes the assumption that all faces have the
311  // same number of dofs
312  AssertDimension(this->n_unique_faces(), 1);
313  const unsigned int face_no = 0;
314 
315  const unsigned int this_dpq = this->n_dofs_per_quad(face_no);
316  const unsigned int other_dpq = fe_other.n_dofs_per_quad(face_no);
317 
318  // we deal with hierarchical 1d polynomials where dofs are enumerated
319  // increasingly. Thus we return a vector of pairs for the first N-1, where
320  // N is minimum number of dofs_per_line for each FE_Q_Hierarchical.
321  std::vector<std::pair<unsigned int, unsigned int>> res;
322  for (unsigned int i = 0; i < std::min(this_dpq, other_dpq); ++i)
323  res.emplace_back(i, i);
324 
325  return res;
326  }
327  else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
328  {
329  // the FE_Nothing has no
330  // degrees of freedom, so there
331  // are no equivalencies to be
332  // recorded
333  return std::vector<std::pair<unsigned int, unsigned int>>();
334  }
335  else
336  {
337  Assert(false, ExcNotImplemented());
338  return std::vector<std::pair<unsigned int, unsigned int>>();
339  }
340 }
341 
342 template <int dim>
345  const FiniteElement<dim> &fe_other,
346  const unsigned int codim) const
347 {
348  Assert(codim <= dim, ExcImpossibleInDim(dim));
349 
350  // vertex/line/face domination
351  // (if fe_other is derived from FE_DGQ)
352  // ------------------------------------
353  if (codim > 0)
354  if (dynamic_cast<const FE_DGQ<dim> *>(&fe_other) != nullptr)
355  // there are no requirements between continuous and discontinuous elements
357 
358  // vertex/line/face domination
359  // (if fe_other is not derived from FE_DGQ)
360  // & cell domination
361  // ----------------------------------------
362  if (const FE_Q_Hierarchical<dim> *fe_hierarchical_other =
363  dynamic_cast<const FE_Q_Hierarchical<dim> *>(&fe_other))
364  {
365  if (this->degree < fe_hierarchical_other->degree)
367  else if (this->degree == fe_hierarchical_other->degree)
369  else
371  }
372  else if (const FE_Nothing<dim> *fe_nothing =
373  dynamic_cast<const FE_Nothing<dim> *>(&fe_other))
374  {
375  if (fe_nothing->is_dominating())
377  else
378  // the FE_Nothing has no degrees of freedom and it is typically used
379  // in a context where we don't require any continuity along the
380  // interface
382  }
383 
384  Assert(false, ExcNotImplemented());
386 }
387 
388 
389 //---------------------------------------------------------------------------
390 // Auxiliary functions
391 //---------------------------------------------------------------------------
392 
393 
394 template <int dim>
395 void
397  std::vector<FullMatrix<double>> &dofs_cell,
398  std::vector<FullMatrix<double>> &dofs_subcell) const
399 {
400  const unsigned int dofs_1d =
401  2 * this->n_dofs_per_vertex() + this->n_dofs_per_line();
402 
403  // The dofs_subcell matrices are transposed
404  // (4.19), (4.21) and (4.27),(4.28),(4.30) in
405  // Demkowicz, Oden, Rachowicz, Hardy, CMAMAE 77, 79-112, 1989
406  // so that
407  // DoFs_c(j) = dofs_subcell[c](j,k) dofs_cell(k)
408 
409  // TODO: The dofs_subcell shall differ by a factor 2^p due to shape functions
410  // defined on [0,1] instead of [-1,1]. However, that does not seem to be
411  // the case. Perhaps this factor is added later on in auxiliary functions
412  // which use these matrices.
413 
414  // dofs_cells[0](j,k):
415  // 0 1 | 2 3 4.
416  // 0 1 0 | .
417  // 1 0 0 | .
418  // -------------------
419  // 2 \ .
420  // 3 \ 2^k * k! / (k-j)!j!
421  // 4 \ .
422 
423  // dofs_subcells[0](j,k):
424  // 0 1 | 2 3 4 5 6 .
425  // 0 1 0 | .
426  // 1 1/2 1/2 | -1 0 -1 0 -1.
427  // -------------------------------
428  // 2 \ .
429  // 3 \ (-1)^(k+j)/ 2^k * k!/(k-j)!j!
430  // 4 \ .
431 
432  // dofs_cells[1](j,k):
433  // 0 1 | 2 3 4.
434  // 0 0 0 | .
435  // 1 0 1 | .
436  // -------------------
437  // 2 \ .
438  // 3 \ (-1)^(k+j) * 2^k * k!/(k-j)!j!
439  // 4 \.
440 
441  // dofs_subcells[1](j,k):
442  // 0 1 | 2 3 4 5 6 .
443  // 0 1/2 1/2 | -1 0 -1 0 -1.
444  // 1 0 1 | .
445  // -------------------------------
446  // 2 \ .
447  // 3 \ 1/ 2^k * k!/(k-j)!j!
448  // 4 .
449 
450  for (unsigned int c = 0; c < GeometryInfo<1>::max_children_per_cell; ++c)
451  for (unsigned int j = 0; j < dofs_1d; ++j)
452  for (unsigned int k = 0; k < dofs_1d; ++k)
453  {
454  // upper diagonal block
455  if ((j <= 1) && (k <= 1))
456  {
457  if (((c == 0) && (j == 0) && (k == 0)) ||
458  ((c == 1) && (j == 1) && (k == 1)))
459  dofs_cell[c](j, k) = 1.;
460  else
461  dofs_cell[c](j, k) = 0.;
462 
463  if (((c == 0) && (j == 1)) || ((c == 1) && (j == 0)))
464  dofs_subcell[c](j, k) = .5;
465  else if (((c == 0) && (k == 0)) || ((c == 1) && (k == 1)))
466  dofs_subcell[c](j, k) = 1.;
467  else
468  dofs_subcell[c](j, k) = 0.;
469  }
470  // upper right block
471  else if ((j <= 1) && (k >= 2))
472  {
473  if (((c == 0) && (j == 1) && ((k % 2) == 0)) ||
474  ((c == 1) && (j == 0) && ((k % 2) == 0)))
475  dofs_subcell[c](j, k) = -1.;
476  }
477  // upper diagonal block
478  else if ((j >= 2) && (k >= 2) && (j <= k))
479  {
480  double factor = 1.;
481  for (unsigned int i = 1; i <= j; ++i)
482  factor *=
483  (static_cast<double>(k - i + 1)) / (static_cast<double>(i));
484  // factor == k * (k-1) * ... * (k-j+1) / j! = k! / (k-j)! / j!
485  if (c == 0)
486  {
487  dofs_subcell[c](j, k) =
488  ((k + j) % 2 == 0) ?
489  std::pow(.5, static_cast<double>(k)) * factor :
490  -std::pow(.5, static_cast<double>(k)) * factor;
491  dofs_cell[c](j, k) =
492  std::pow(2., static_cast<double>(j)) * factor;
493  }
494  else
495  {
496  dofs_subcell[c](j, k) =
497  std::pow(.5, static_cast<double>(k)) * factor;
498  dofs_cell[c](j, k) =
499  ((k + j) % 2 == 0) ?
500  std::pow(2., static_cast<double>(j)) * factor :
501  -std::pow(2., static_cast<double>(j)) * factor;
502  }
503  }
504  }
505 }
506 
507 
508 
509 template <int dim>
510 void
512  const std::vector<FullMatrix<double>> &dofs_subcell)
513 {
514  const unsigned int dofs_1d =
515  2 * this->n_dofs_per_vertex() + this->n_dofs_per_line();
516 
517  this->interface_constraints.TableBase<2, double>::reinit(
518  this->interface_constraints_size());
519 
520  switch (dim)
521  {
522  case 1:
523  {
524  // no constraints in 1d
525  break;
526  }
527 
528  case 2:
529  {
530  // vertex node
531  for (unsigned int i = 0; i < dofs_1d; ++i)
532  this->interface_constraints(0, i) = dofs_subcell[0](1, i);
533  // edge nodes
534  for (unsigned int c = 0; c < GeometryInfo<1>::max_children_per_cell;
535  ++c)
536  for (unsigned int i = 0; i < dofs_1d; ++i)
537  for (unsigned int j = 2; j < dofs_1d; ++j)
538  this->interface_constraints(1 + c * (this->degree - 1) + j - 2,
539  i) = dofs_subcell[c](j, i);
540  break;
541  }
542 
543  case 3:
544  {
545  for (unsigned int i = 0; i < dofs_1d * dofs_1d; ++i)
546  {
547  // center vertex node
548  this->interface_constraints(0, face_renumber[i]) =
549  dofs_subcell[0](1, i % dofs_1d) *
550  dofs_subcell[0](1, (i - (i % dofs_1d)) / dofs_1d);
551 
552  // boundary vertex nodes
553  this->interface_constraints(1, face_renumber[i]) =
554  dofs_subcell[0](0, i % dofs_1d) *
555  dofs_subcell[0](1, (i - (i % dofs_1d)) / dofs_1d);
556  this->interface_constraints(2, face_renumber[i]) =
557  dofs_subcell[1](1, i % dofs_1d) *
558  dofs_subcell[0](1, (i - (i % dofs_1d)) / dofs_1d);
559  this->interface_constraints(3, face_renumber[i]) =
560  dofs_subcell[0](1, i % dofs_1d) *
561  dofs_subcell[0](0, (i - (i % dofs_1d)) / dofs_1d);
562  this->interface_constraints(4, face_renumber[i]) =
563  dofs_subcell[1](0, i % dofs_1d) *
564  dofs_subcell[1](1, (i - (i % dofs_1d)) / dofs_1d);
565 
566  // interior edges
567  for (unsigned int j = 0; j < (this->degree - 1); ++j)
568  {
569  this->interface_constraints(5 + j, face_renumber[i]) =
570  dofs_subcell[0](1, i % dofs_1d) *
571  dofs_subcell[0](2 + j, (i - (i % dofs_1d)) / dofs_1d);
572  this->interface_constraints(5 + (this->degree - 1) + j,
573  face_renumber[i]) =
574  dofs_subcell[0](1, i % dofs_1d) *
575  dofs_subcell[1](2 + j, (i - (i % dofs_1d)) / dofs_1d);
576  this->interface_constraints(5 + 2 * (this->degree - 1) + j,
577  face_renumber[i]) =
578  dofs_subcell[0](2 + j, i % dofs_1d) *
579  dofs_subcell[1](0, (i - (i % dofs_1d)) / dofs_1d);
580  this->interface_constraints(5 + 3 * (this->degree - 1) + j,
581  face_renumber[i]) =
582  dofs_subcell[1](2 + j, i % dofs_1d) *
583  dofs_subcell[0](1, (i - (i % dofs_1d)) / dofs_1d);
584  }
585 
586  // boundary edges
587  for (unsigned int j = 0; j < (this->degree - 1); ++j)
588  {
589  // left edge
590  this->interface_constraints(5 + 4 * (this->degree - 1) + j,
591  face_renumber[i]) =
592  dofs_subcell[0](0, i % dofs_1d) *
593  dofs_subcell[0](2 + j, (i - (i % dofs_1d)) / dofs_1d);
594  this->interface_constraints(5 + 4 * (this->degree - 1) +
595  (this->degree - 1) + j,
596  face_renumber[i]) =
597  dofs_subcell[0](0, i % dofs_1d) *
598  dofs_subcell[1](2 + j, (i - (i % dofs_1d)) / dofs_1d);
599  // right edge
600  this->interface_constraints(5 + 4 * (this->degree - 1) +
601  2 * (this->degree - 1) + j,
602  face_renumber[i]) =
603  dofs_subcell[1](1, i % dofs_1d) *
604  dofs_subcell[0](2 + j, (i - (i % dofs_1d)) / dofs_1d);
605  this->interface_constraints(5 + 4 * (this->degree - 1) +
606  3 * (this->degree - 1) + j,
607  face_renumber[i]) =
608  dofs_subcell[1](1, i % dofs_1d) *
609  dofs_subcell[1](2 + j, (i - (i % dofs_1d)) / dofs_1d);
610  // bottom edge
611  this->interface_constraints(5 + 4 * (this->degree - 1) +
612  4 * (this->degree - 1) + j,
613  face_renumber[i]) =
614  dofs_subcell[0](2 + j, i % dofs_1d) *
615  dofs_subcell[0](0, (i - (i % dofs_1d)) / dofs_1d);
616  this->interface_constraints(5 + 4 * (this->degree - 1) +
617  5 * (this->degree - 1) + j,
618  face_renumber[i]) =
619  dofs_subcell[1](2 + j, i % dofs_1d) *
620  dofs_subcell[0](0, (i - (i % dofs_1d)) / dofs_1d);
621  // top edge
622  this->interface_constraints(5 + 4 * (this->degree - 1) +
623  6 * (this->degree - 1) + j,
624  face_renumber[i]) =
625  dofs_subcell[0](2 + j, i % dofs_1d) *
626  dofs_subcell[1](1, (i - (i % dofs_1d)) / dofs_1d);
627  this->interface_constraints(5 + 4 * (this->degree - 1) +
628  7 * (this->degree - 1) + j,
629  face_renumber[i]) =
630  dofs_subcell[1](2 + j, i % dofs_1d) *
631  dofs_subcell[1](1, (i - (i % dofs_1d)) / dofs_1d);
632  }
633 
634  // interior faces
635  for (unsigned int j = 0; j < (this->degree - 1); ++j)
636  for (unsigned int k = 0; k < (this->degree - 1); ++k)
637  {
638  // subcell 0
639  this->interface_constraints(5 + 12 * (this->degree - 1) +
640  j + k * (this->degree - 1),
641  face_renumber[i]) =
642  dofs_subcell[0](2 + j, i % dofs_1d) *
643  dofs_subcell[0](2 + k, (i - (i % dofs_1d)) / dofs_1d);
644  // subcell 1
645  this->interface_constraints(5 + 12 * (this->degree - 1) +
646  j + k * (this->degree - 1) +
647  (this->degree - 1) *
648  (this->degree - 1),
649  face_renumber[i]) =
650  dofs_subcell[1](2 + j, i % dofs_1d) *
651  dofs_subcell[0](2 + k, (i - (i % dofs_1d)) / dofs_1d);
652  // subcell 2
653  this->interface_constraints(5 + 12 * (this->degree - 1) +
654  j + k * (this->degree - 1) +
655  2 * (this->degree - 1) *
656  (this->degree - 1),
657  face_renumber[i]) =
658  dofs_subcell[0](2 + j, i % dofs_1d) *
659  dofs_subcell[1](2 + k, (i - (i % dofs_1d)) / dofs_1d);
660  // subcell 3
661  this->interface_constraints(5 + 12 * (this->degree - 1) +
662  j + k * (this->degree - 1) +
663  3 * (this->degree - 1) *
664  (this->degree - 1),
665  face_renumber[i]) =
666  dofs_subcell[1](2 + j, i % dofs_1d) *
667  dofs_subcell[1](2 + k, (i - (i % dofs_1d)) / dofs_1d);
668  }
669  }
670  break;
671  }
672 
673  default:
674  Assert(false, ExcNotImplemented());
675  }
676 }
677 
678 
679 
680 template <int dim>
681 void
683  const std::vector<FullMatrix<double>> &dofs_cell,
684  const std::vector<FullMatrix<double>> &dofs_subcell)
685 {
686  unsigned int iso = RefinementCase<dim>::isotropic_refinement - 1;
687 
688  const unsigned int dofs_1d =
689  2 * this->n_dofs_per_vertex() + this->n_dofs_per_line();
690  TensorProductPolynomials<dim> *poly_space_derived_ptr =
691  dynamic_cast<TensorProductPolynomials<dim> *>(this->poly_space.get());
692  const std::vector<unsigned int> &renumber =
693  poly_space_derived_ptr->get_numbering();
694 
695  for (unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell; ++c)
696  {
697  this->prolongation[iso][c].reinit(this->n_dofs_per_cell(),
698  this->n_dofs_per_cell());
699  this->restriction[iso][c].reinit(this->n_dofs_per_cell(),
700  this->n_dofs_per_cell());
701  }
702 
703  // the 1d case is particularly
704  // simple, so special case it:
705  if (dim == 1)
706  {
707  for (unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
708  ++c)
709  {
710  this->prolongation[iso][c].fill(dofs_subcell[c]);
711  this->restriction[iso][c].fill(dofs_cell[c]);
712  }
713  return;
714  }
715 
716  // for higher dimensions, things
717  // are a little more tricky:
718 
719  // j loops over dofs in the
720  // subcell. These are the rows in
721  // the embedding matrix.
722  //
723  // i loops over the dofs in the
724  // parent cell. These are the
725  // columns in the embedding matrix.
726  for (unsigned int j = 0; j < this->n_dofs_per_cell(); ++j)
727  for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
728  switch (dim)
729  {
730  case 2:
731  {
732  for (unsigned int c = 0;
733  c < GeometryInfo<2>::max_children_per_cell;
734  ++c)
735  {
736  // left/right line: 0/1
737  const unsigned int c0 = c % 2;
738  // bottom/top line: 0/1
739  const unsigned int c1 = c / 2;
740 
741  this->prolongation[iso][c](j, i) =
742  dofs_subcell[c0](renumber[j] % dofs_1d,
743  renumber[i] % dofs_1d) *
744  dofs_subcell[c1]((renumber[j] - (renumber[j] % dofs_1d)) /
745  dofs_1d,
746  (renumber[i] - (renumber[i] % dofs_1d)) /
747  dofs_1d);
748 
749  this->restriction[iso][c](j, i) =
750  dofs_cell[c0](renumber[j] % dofs_1d,
751  renumber[i] % dofs_1d) *
752  dofs_cell[c1]((renumber[j] - (renumber[j] % dofs_1d)) /
753  dofs_1d,
754  (renumber[i] - (renumber[i] % dofs_1d)) /
755  dofs_1d);
756  }
757  break;
758  }
759 
760  case 3:
761  {
762  for (unsigned int c = 0;
763  c < GeometryInfo<3>::max_children_per_cell;
764  ++c)
765  {
766  // left/right face: 0/1
767  const unsigned int c0 = c % 2;
768  // front/back face: 0/1
769  const unsigned int c1 = (c % 4) / 2;
770  // bottom/top face: 0/1
771  const unsigned int c2 = c / 4;
772 
773  this->prolongation[iso][c](j, i) =
774  dofs_subcell[c0](renumber[j] % dofs_1d,
775  renumber[i] % dofs_1d) *
776  dofs_subcell[c1](
777  ((renumber[j] - (renumber[j] % dofs_1d)) / dofs_1d) %
778  dofs_1d,
779  ((renumber[i] - (renumber[i] % dofs_1d)) / dofs_1d) %
780  dofs_1d) *
781  dofs_subcell[c2](
782  ((renumber[j] - (renumber[j] % dofs_1d)) / dofs_1d -
783  (((renumber[j] - (renumber[j] % dofs_1d)) / dofs_1d) %
784  dofs_1d)) /
785  dofs_1d,
786  ((renumber[i] - (renumber[i] % dofs_1d)) / dofs_1d -
787  (((renumber[i] - (renumber[i] % dofs_1d)) / dofs_1d) %
788  dofs_1d)) /
789  dofs_1d);
790 
791  this->restriction[iso][c](j, i) =
792  dofs_cell[c0](renumber[j] % dofs_1d,
793  renumber[i] % dofs_1d) *
794  dofs_cell[c1](
795  ((renumber[j] - (renumber[j] % dofs_1d)) / dofs_1d) %
796  dofs_1d,
797  ((renumber[i] - (renumber[i] % dofs_1d)) / dofs_1d) %
798  dofs_1d) *
799  dofs_cell[c2](
800  ((renumber[j] - (renumber[j] % dofs_1d)) / dofs_1d -
801  (((renumber[j] - (renumber[j] % dofs_1d)) / dofs_1d) %
802  dofs_1d)) /
803  dofs_1d,
804  ((renumber[i] - (renumber[i] % dofs_1d)) / dofs_1d -
805  (((renumber[i] - (renumber[i] % dofs_1d)) / dofs_1d) %
806  dofs_1d)) /
807  dofs_1d);
808  }
809  break;
810  }
811 
812  default:
813  Assert(false, ExcNotImplemented());
814  }
815 }
816 
817 
818 
819 template <int dim>
820 void
822 {
823  // number of points: (degree+1)^dim
824  unsigned int n = this->degree + 1;
825  for (unsigned int i = 1; i < dim; ++i)
826  n *= this->degree + 1;
827 
828  this->generalized_support_points.resize(n);
829 
830  TensorProductPolynomials<dim> *poly_space_derived_ptr =
831  dynamic_cast<TensorProductPolynomials<dim> *>(this->poly_space.get());
832  const std::vector<unsigned int> &index_map_inverse =
833  poly_space_derived_ptr->get_numbering_inverse();
834 
835  Point<dim> p;
836  // the method of numbering allows
837  // each dof to be associated with a
838  // support point. There is
839  // only one support point per
840  // vertex, line, quad, hex, etc.
841  //
842  // note, however, that the support
843  // points thus associated with
844  // shape functions are not unique:
845  // the linear shape functions are
846  // associated with the vertices,
847  // but all others are associated
848  // with either line, quad, or hex
849  // midpoints, and there may be
850  // multiple shape functions
851  // associated with them. there
852  // really is no other useful
853  // numbering, since the
854  // hierarchical shape functions do
855  // not vanish at all-but-one
856  // interpolation points (like the
857  // Lagrange functions used in
858  // FE_Q), so there's not much we
859  // can do here.
860 
861  // TODO shouldn't we just at least make support points unique,
862  // even though the delta property is not satisfied for this FE?
863  unsigned int k = 0;
864  for (unsigned int iz = 0; iz <= ((dim > 2) ? this->degree : 0); ++iz)
865  for (unsigned int iy = 0; iy <= ((dim > 1) ? this->degree : 0); ++iy)
866  for (unsigned int ix = 0; ix <= this->degree; ++ix)
867  {
868  if (ix == 0)
869  p(0) = 0.;
870  else if (ix == 1)
871  p(0) = 1.;
872  else
873  p(0) = .5;
874  if (dim > 1)
875  {
876  if (iy == 0)
877  p(1) = 0.;
878  else if (iy == 1)
879  p(1) = 1.;
880  else
881  p(1) = .5;
882  }
883  if (dim > 2)
884  {
885  if (iz == 0)
886  p(2) = 0.;
887  else if (iz == 1)
888  p(2) = 1.;
889  else
890  p(2) = .5;
891  }
892  this->generalized_support_points[index_map_inverse[k++]] = p;
893  }
894 }
895 
896 
897 
898 template <>
899 void
901 {
902  // no faces in 1d, so nothing to do
903 }
904 
905 
906 template <>
907 void
909  const FiniteElement<1, 1> & /*x_source_fe*/,
910  FullMatrix<double> & /*interpolation_matrix*/,
911  const unsigned int) const
912 {
913  Assert(false, ExcImpossibleInDim(1));
914 }
915 
916 
917 template <>
918 void
920  const FiniteElement<1, 1> & /*x_source_fe*/,
921  const unsigned int /*subface*/,
922  FullMatrix<double> & /*interpolation_matrix*/,
923  const unsigned int) const
924 {
925  Assert(false, ExcImpossibleInDim(1));
926 }
927 
928 
929 
930 template <int dim>
931 void
933  const FiniteElement<dim> &x_source_fe,
934  FullMatrix<double> & interpolation_matrix,
935  const unsigned int face_no) const
936 {
937  // this is only implemented, if the
938  // source FE is also a
939  // Q_Hierarchical element
940  using FEQHierarchical = FE_Q_Hierarchical<dim>;
941  AssertThrow((x_source_fe.get_name().find("FE_Q_Hierarchical<") == 0) ||
942  (dynamic_cast<const FEQHierarchical *>(&x_source_fe) !=
943  nullptr),
945 
946  Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
947  ExcDimensionMismatch(interpolation_matrix.n(),
948  this->n_dofs_per_face(face_no)));
949  Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
950  ExcDimensionMismatch(interpolation_matrix.m(),
951  x_source_fe.n_dofs_per_face(face_no)));
952 
953  // ok, source is a Q_Hierarchical element, so
954  // we will be able to do the work
955  const FE_Q_Hierarchical<dim> &source_fe =
956  dynamic_cast<const FE_Q_Hierarchical<dim> &>(x_source_fe);
957  (void)source_fe;
958 
959  // Make sure, that the element,
960  // for which the DoFs should be
961  // constrained is the one with
962  // the higher polynomial degree.
963  // Actually the procedure will work
964  // also if this assertion is not
965  // satisfied. But the matrices
966  // produced in that case might
967  // lead to problems in the
968  // hp-procedures, which use this
969  // method.
970  Assert(this->n_dofs_per_face(face_no) <= source_fe.n_dofs_per_face(face_no),
972  interpolation_matrix = 0;
973 
974  switch (dim)
975  {
976  case 2:
977  {
978  // In 2-dimension the constraints are trivial.
979  // First this->dofs_per_face DoFs of the constrained
980  // element are made equal to the current (dominating)
981  // element, which corresponds to 1 on diagonal of the matrix.
982  // DoFs which correspond to higher polynomials
983  // are zeroed (zero rows in the matrix).
984  for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
985  interpolation_matrix(i, i) = 1;
986 
987  break;
988  }
989 
990  case 3:
991  {
992  for (unsigned int i = 0; i < GeometryInfo<3>::vertices_per_face; ++i)
993  interpolation_matrix(i, i) = 1;
994 
995  for (unsigned int i = 0; i < this->degree - 1; ++i)
996  {
997  for (unsigned int j = 0; j < GeometryInfo<3>::lines_per_face; ++j)
998  interpolation_matrix(i + j * (x_source_fe.degree - 1) +
1000  i + j * (this->degree - 1) +
1002 
1003  for (unsigned int j = 0; j < this->degree - 1; ++j)
1004  interpolation_matrix((i + GeometryInfo<3>::lines_per_face) *
1005  (x_source_fe.degree - 1) +
1008  (this->degree - 1) +
1010  1;
1011  }
1012  }
1013  }
1014 }
1015 
1016 
1017 
1018 template <int dim>
1019 void
1021  const FiniteElement<dim> &x_source_fe,
1022  const unsigned int subface,
1023  FullMatrix<double> & interpolation_matrix,
1024  const unsigned int face_no) const
1025 {
1026  // this is only implemented, if the
1027  // source FE is also a
1028  // Q_Hierarchical element
1029  using FEQHierarchical = FE_Q_Hierarchical<dim>;
1030  AssertThrow((x_source_fe.get_name().find("FE_Q_Hierarchical<") == 0) ||
1031  (dynamic_cast<const FEQHierarchical *>(&x_source_fe) !=
1032  nullptr),
1034 
1035  Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
1036  ExcDimensionMismatch(interpolation_matrix.n(),
1037  this->n_dofs_per_face(face_no)));
1038  Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
1039  ExcDimensionMismatch(interpolation_matrix.m(),
1040  x_source_fe.n_dofs_per_face(face_no)));
1041 
1042  // ok, source is a Q_Hierarchical element, so
1043  // we will be able to do the work
1044  const FE_Q_Hierarchical<dim> &source_fe =
1045  dynamic_cast<const FE_Q_Hierarchical<dim> &>(x_source_fe);
1046 
1047  // Make sure, that the element,
1048  // for which the DoFs should be
1049  // constrained is the one with
1050  // the higher polynomial degree.
1051  // Actually the procedure will work
1052  // also if this assertion is not
1053  // satisfied. But the matrices
1054  // produced in that case might
1055  // lead to problems in the
1056  // hp-procedures, which use this
1057  // method.
1058  Assert(this->n_dofs_per_face(face_no) <= source_fe.n_dofs_per_face(face_no),
1060 
1061  switch (dim)
1062  {
1063  case 2:
1064  {
1065  switch (subface)
1066  {
1067  case 0:
1068  {
1069  interpolation_matrix(0, 0) = 1.0;
1070  interpolation_matrix(1, 0) = 0.5;
1071  interpolation_matrix(1, 1) = 0.5;
1072 
1073  for (unsigned int dof = 2;
1074  dof < this->n_dofs_per_face(face_no);)
1075  {
1076  interpolation_matrix(1, dof) = -1.0;
1077  dof = dof + 2;
1078  }
1079 
1080  int factorial_i = 1;
1081 
1082  for (unsigned int i = 2; i < this->n_dofs_per_face(face_no);
1083  ++i)
1084  {
1085  interpolation_matrix(i, i) = std::pow(0.5, i);
1086  factorial_i *= i;
1087  int factorial_j = factorial_i;
1088  int factorial_ij = 1;
1089 
1090  for (unsigned int j = i + 1;
1091  j < this->n_dofs_per_face(face_no);
1092  ++j)
1093  {
1094  factorial_ij *= j - i;
1095  factorial_j *= j;
1096 
1097  if (((i + j) & 1) != 0u)
1098  interpolation_matrix(i, j) =
1099  -1.0 * std::pow(0.5, j) * factorial_j /
1100  (factorial_i * factorial_ij);
1101 
1102  else
1103  interpolation_matrix(i, j) =
1104  std::pow(0.5, j) * factorial_j /
1105  (factorial_i * factorial_ij);
1106  }
1107  }
1108 
1109  break;
1110  }
1111 
1112  case 1:
1113  {
1114  interpolation_matrix(0, 0) = 0.5;
1115  interpolation_matrix(0, 1) = 0.5;
1116 
1117  for (unsigned int dof = 2;
1118  dof < this->n_dofs_per_face(face_no);)
1119  {
1120  interpolation_matrix(0, dof) = -1.0;
1121  dof = dof + 2;
1122  }
1123 
1124  interpolation_matrix(1, 1) = 1.0;
1125 
1126  int factorial_i = 1;
1127 
1128  for (unsigned int i = 2; i < this->n_dofs_per_face(face_no);
1129  ++i)
1130  {
1131  interpolation_matrix(i, i) = std::pow(0.5, i);
1132  factorial_i *= i;
1133  int factorial_j = factorial_i;
1134  int factorial_ij = 1;
1135 
1136  for (unsigned int j = i + 1;
1137  j < this->n_dofs_per_face(face_no);
1138  ++j)
1139  {
1140  factorial_ij *= j - i;
1141  factorial_j *= j;
1142  interpolation_matrix(i, j) =
1143  std::pow(0.5, j) * factorial_j /
1144  (factorial_i * factorial_ij);
1145  }
1146  }
1147  }
1148  }
1149 
1150  break;
1151  }
1152 
1153  case 3:
1154  {
1155  switch (subface)
1156  {
1157  case 0:
1158  {
1159  interpolation_matrix(0, 0) = 1.0;
1160  interpolation_matrix(1, 0) = 0.5;
1161  interpolation_matrix(1, 1) = 0.5;
1162  interpolation_matrix(2, 0) = 0.5;
1163  interpolation_matrix(2, 2) = 0.5;
1164 
1165  for (unsigned int i = 0; i < this->degree - 1;)
1166  {
1167  for (unsigned int line = 0;
1168  line < GeometryInfo<3>::lines_per_face;
1169  ++line)
1170  interpolation_matrix(3,
1171  i + line * (this->degree - 1) +
1172  4) = -0.5;
1173 
1174  for (unsigned int j = 0; j < this->degree - 1;)
1175  {
1176  interpolation_matrix(3,
1177  i + (j + 4) * this->degree - j) =
1178  1.0;
1179  j = j + 2;
1180  }
1181 
1182  interpolation_matrix(1, i + 2 * (this->degree + 1)) =
1183  -1.0;
1184  interpolation_matrix(2, i + 4) = -1.0;
1185  i = i + 2;
1186  }
1187 
1188  for (unsigned int vertex = 0;
1189  vertex < GeometryInfo<3>::vertices_per_face;
1190  ++vertex)
1191  interpolation_matrix(3, vertex) = 0.25;
1192 
1193  int factorial_i = 1;
1194 
1195  for (unsigned int i = 2; i <= this->degree; ++i)
1196  {
1197  double tmp = std::pow(0.5, i);
1198  interpolation_matrix(i + 2, i + 2) = tmp;
1199  interpolation_matrix(i + 2 * source_fe.degree,
1200  i + 2 * this->degree) = tmp;
1201  tmp *= 0.5;
1202  interpolation_matrix(i + source_fe.degree + 1, i + 2) =
1203  tmp;
1204  interpolation_matrix(i + source_fe.degree + 1,
1205  i + this->degree + 1) = tmp;
1206  interpolation_matrix(i + 3 * source_fe.degree - 1,
1207  i + 2 * this->degree) = tmp;
1208  interpolation_matrix(i + 3 * source_fe.degree - 1,
1209  i + 3 * this->degree - 1) = tmp;
1210  tmp *= -2.0;
1211 
1212  for (unsigned int j = 0; j < this->degree - 1;)
1213  {
1214  interpolation_matrix(i + source_fe.degree + 1,
1215  (i + 2) * this->degree + j + 2 -
1216  i) = tmp;
1217  interpolation_matrix(i + 3 * source_fe.degree - 1,
1218  i + (j + 4) * this->degree - j -
1219  2) = tmp;
1220  j = j + 2;
1221  }
1222 
1223  int factorial_k = 1;
1224 
1225  for (unsigned int j = 2; j <= this->degree; ++j)
1226  {
1227  interpolation_matrix(i + (j + 2) * source_fe.degree -
1228  j,
1229  i + (j + 2) * this->degree - j) =
1230  std::pow(0.5, i + j);
1231  factorial_k *= j;
1232  int factorial_kl = 1;
1233  int factorial_l = factorial_k;
1234 
1235  for (unsigned int k = j + 1; k < this->degree; ++k)
1236  {
1237  factorial_kl *= k - j;
1238  factorial_l *= k;
1239 
1240  if (((j + k) & 1) != 0u)
1241  interpolation_matrix(
1242  i + (j + 2) * source_fe.degree - j,
1243  i + (k + 2) * this->degree - k) =
1244  -1.0 * std::pow(0.5, i + k) * factorial_l /
1245  (factorial_k * factorial_kl);
1246 
1247  else
1248  interpolation_matrix(
1249  i + (j + 2) * source_fe.degree - j,
1250  i + (k + 2) * this->degree - k) =
1251  std::pow(0.5, i + k) * factorial_l /
1252  (factorial_k * factorial_kl);
1253  }
1254  }
1255 
1256  factorial_i *= i;
1257  int factorial_j = factorial_i;
1258  int factorial_ij = 1;
1259 
1260  for (unsigned int j = i + 1; j <= this->degree; ++j)
1261  {
1262  factorial_ij *= j - i;
1263  factorial_j *= j;
1264 
1265  if (((i + j) & 1) != 0u)
1266  {
1267  tmp = -1.0 * std::pow(0.5, j) * factorial_j /
1268  (factorial_i * factorial_ij);
1269  interpolation_matrix(i + 2, j + 2) = tmp;
1270  interpolation_matrix(i + 2 * source_fe.degree,
1271  j + 2 * this->degree) = tmp;
1272  factorial_k = 1;
1273 
1274  for (unsigned int k = 2; k <= this->degree; ++k)
1275  {
1276  interpolation_matrix(
1277  i + (k + 2) * source_fe.degree - k,
1278  j + (k + 2) * this->degree - k) =
1279  tmp * std::pow(0.5, k);
1280  factorial_k *= k;
1281  int factorial_l = factorial_k;
1282  int factorial_kl = 1;
1283 
1284  for (unsigned int l = k + 1;
1285  l <= this->degree;
1286  ++l)
1287  {
1288  factorial_kl *= l - k;
1289  factorial_l *= l;
1290 
1291  if (((k + l) & 1) != 0u)
1292  interpolation_matrix(
1293  i + (k + 2) * source_fe.degree - k,
1294  j + (l + 2) * this->degree - l) =
1295  -1.0 * tmp * std::pow(0.5, l) *
1296  factorial_l /
1297  (factorial_k * factorial_kl);
1298 
1299  else
1300  interpolation_matrix(
1301  i + (k + 2) * source_fe.degree - k,
1302  j + (l + 2) * this->degree - l) =
1303  tmp * std::pow(0.5, l) * factorial_l /
1304  (factorial_k * factorial_kl);
1305  }
1306  }
1307 
1308  tmp *= 0.5;
1309  interpolation_matrix(i + source_fe.degree + 1,
1310  j + 2) = tmp;
1311  interpolation_matrix(i + source_fe.degree + 1,
1312  j + this->degree + 1) = tmp;
1313  interpolation_matrix(i + 3 * source_fe.degree - 1,
1314  j + 2 * this->degree) = tmp;
1315  interpolation_matrix(i + 3 * source_fe.degree - 1,
1316  j + 3 * this->degree - 1) =
1317  tmp;
1318  tmp *= -2.0;
1319 
1320  for (unsigned int k = 0; k < this->degree - 1;)
1321  {
1322  interpolation_matrix(i + source_fe.degree + 1,
1323  (j + 2) * this->degree +
1324  k + 2 - j) = tmp;
1325  interpolation_matrix(
1326  i + 3 * source_fe.degree - 1,
1327  j + (k + 4) * this->degree - k - 2) = tmp;
1328  k = k + 2;
1329  }
1330  }
1331  else
1332  {
1333  tmp = std::pow(0.5, j) * factorial_j /
1334  (factorial_i * factorial_ij);
1335  interpolation_matrix(i + 2, j + 2) = tmp;
1336  interpolation_matrix(i + 2 * source_fe.degree,
1337  j + 2 * this->degree) = tmp;
1338  factorial_k = 1;
1339 
1340  for (unsigned int k = 2; k <= this->degree; ++k)
1341  {
1342  interpolation_matrix(
1343  i + (k + 2) * source_fe.degree - k,
1344  j + (k + 2) * this->degree - k) =
1345  tmp * std::pow(0.5, k);
1346  factorial_k *= k;
1347  int factorial_l = factorial_k;
1348  int factorial_kl = 1;
1349 
1350  for (unsigned int l = k + 1;
1351  l <= this->degree;
1352  ++l)
1353  {
1354  factorial_kl *= l - k;
1355  factorial_l *= l;
1356 
1357  if (((k + l) & 1) != 0u)
1358  interpolation_matrix(
1359  i + (k + 2) * source_fe.degree - k,
1360  j + (l + 2) * this->degree - l) =
1361  -1.0 * tmp * std::pow(0.5, l) *
1362  factorial_l /
1363  (factorial_k * factorial_kl);
1364 
1365  else
1366  interpolation_matrix(
1367  i + (k + 2) * source_fe.degree - k,
1368  j + (l + 2) * this->degree - l) =
1369  tmp * std::pow(0.5, l) * factorial_l /
1370  (factorial_k * factorial_kl);
1371  }
1372  }
1373 
1374  tmp *= 0.5;
1375  interpolation_matrix(i + source_fe.degree + 1,
1376  j + 2) = tmp;
1377  interpolation_matrix(i + source_fe.degree + 1,
1378  j + this->degree + 1) = tmp;
1379  interpolation_matrix(i + 3 * source_fe.degree - 1,
1380  j + 2 * this->degree) = tmp;
1381  interpolation_matrix(i + 3 * source_fe.degree - 1,
1382  j + 3 * this->degree - 1) =
1383  tmp;
1384  tmp *= -2.0;
1385 
1386  for (unsigned int k = 0; k < this->degree - 1;)
1387  {
1388  interpolation_matrix(i + source_fe.degree + 1,
1389  (j + 2) * this->degree +
1390  k + 2 - j) = tmp;
1391  interpolation_matrix(
1392  i + 3 * source_fe.degree - 1,
1393  j + (k + 4) * this->degree - k - 2) = tmp;
1394  k = k + 2;
1395  }
1396  }
1397  }
1398  }
1399 
1400  break;
1401  }
1402 
1403  case 1:
1404  {
1405  interpolation_matrix(0, 0) = 0.5;
1406  interpolation_matrix(0, 1) = 0.5;
1407  interpolation_matrix(1, 1) = 1.0;
1408  interpolation_matrix(3, 1) = 0.5;
1409  interpolation_matrix(3, 3) = 0.5;
1410 
1411  for (unsigned int i = 0; i < this->degree - 1;)
1412  {
1413  for (unsigned int line = 0;
1414  line < GeometryInfo<3>::lines_per_face;
1415  ++line)
1416  interpolation_matrix(2,
1417  i + line * (this->degree - 1) +
1418  4) = -0.5;
1419 
1420  for (unsigned int j = 0; j < this->degree - 1;)
1421  {
1422  interpolation_matrix(2,
1423  i + (j + 4) * this->degree - j) =
1424  1.0;
1425  j = j + 2;
1426  }
1427 
1428  interpolation_matrix(0, i + 2 * (this->degree + 1)) =
1429  -1.0;
1430  interpolation_matrix(3, i + this->degree + 3) = -1.0;
1431  i = i + 2;
1432  }
1433 
1434  for (unsigned int vertex = 0;
1435  vertex < GeometryInfo<3>::vertices_per_face;
1436  ++vertex)
1437  interpolation_matrix(2, vertex) = 0.25;
1438 
1439  int factorial_i = 1;
1440 
1441  for (unsigned int i = 2; i <= this->degree; ++i)
1442  {
1443  double tmp = std::pow(0.5, i + 1);
1444  interpolation_matrix(i + 2, i + 2) = tmp;
1445  interpolation_matrix(i + 2, i + this->degree + 1) = tmp;
1446  interpolation_matrix(i + 3 * source_fe.degree - 1,
1447  i + 2 * this->degree) = tmp;
1448  interpolation_matrix(i + 3 * source_fe.degree - 1,
1449  i + 3 * this->degree - 1) = tmp;
1450  tmp *= -2.0;
1451 
1452  for (unsigned int j = 0; j < this->degree - 1;)
1453  {
1454  interpolation_matrix(i + 2,
1455  j + (i + 2) * this->degree + 2 -
1456  i) = tmp;
1457  interpolation_matrix(i + 3 * source_fe.degree - 1,
1458  i + (j + 4) * this->degree - j -
1459  2) = tmp;
1460  j = j + 2;
1461  }
1462 
1463  tmp *= -1.0;
1464  interpolation_matrix(i + source_fe.degree + 1,
1465  i + this->degree + 1) = tmp;
1466  interpolation_matrix(i + 2 * source_fe.degree,
1467  i + 2 * this->degree) = tmp;
1468  factorial_i *= i;
1469  int factorial_j = factorial_i;
1470  int factorial_ij = 1;
1471 
1472  for (unsigned int j = i + 1; j <= this->degree; ++j)
1473  {
1474  factorial_ij *= j - i;
1475  factorial_j *= j;
1476  tmp = std::pow(0.5, j) * factorial_j /
1477  (factorial_i * factorial_ij);
1478  interpolation_matrix(i + 2 * source_fe.degree,
1479  j + 2 * this->degree) = tmp;
1480  int factorial_k = 1;
1481 
1482  for (unsigned int k = 2; k <= this->degree; ++k)
1483  {
1484  interpolation_matrix(
1485  i + (k + 2) * source_fe.degree - k,
1486  j + (k + 2) * this->degree - k) =
1487  tmp * std::pow(0.5, k);
1488  factorial_k *= k;
1489  int factorial_l = factorial_k;
1490  int factorial_kl = 1;
1491 
1492  for (unsigned int l = k + 1; l <= this->degree;
1493  ++l)
1494  {
1495  factorial_kl *= l - k;
1496  factorial_l *= l;
1497 
1498  if (((k + l) & 1) != 0u)
1499  interpolation_matrix(
1500  i + (k + 2) * source_fe.degree - k,
1501  j + (l + 2) * this->degree - l) =
1502  -1.0 * tmp * std::pow(0.5, l) *
1503  factorial_l /
1504  (factorial_k * factorial_kl);
1505 
1506  else
1507  interpolation_matrix(
1508  i + (k + 2) * source_fe.degree - k,
1509  j + (l + 2) * this->degree - l) =
1510  tmp * std::pow(0.5, l) * factorial_l /
1511  (factorial_k * factorial_kl);
1512  }
1513  }
1514 
1515  tmp *= -1.0;
1516 
1517  for (unsigned int k = 0; k < this->degree - 1;)
1518  {
1519  interpolation_matrix(i + 3 * source_fe.degree - 1,
1520  j + (k + 4) * this->degree -
1521  k - 2) = tmp;
1522  k = k + 2;
1523  }
1524 
1525  tmp *= -0.5;
1526  interpolation_matrix(i + 3 * source_fe.degree - 1,
1527  j + 2 * this->degree) = tmp;
1528  interpolation_matrix(i + 3 * source_fe.degree - 1,
1529  j + 3 * this->degree - 1) = tmp;
1530 
1531  if (((i + j) & 1) != 0u)
1532  tmp *= -1.0;
1533 
1534  interpolation_matrix(i + 2, j + 2) = tmp;
1535  interpolation_matrix(i + 2, j + this->degree + 1) =
1536  tmp;
1537  interpolation_matrix(i + source_fe.degree + 1,
1538  j + this->degree + 1) =
1539  2.0 * tmp;
1540  tmp *= -2.0;
1541 
1542  for (unsigned int k = 0; k < this->degree - 1;)
1543  {
1544  interpolation_matrix(i + 2,
1545  k + (j + 2) * this->degree +
1546  2 - j) = tmp;
1547  k = k + 2;
1548  }
1549  }
1550 
1551  int factorial_k = 1;
1552 
1553  for (unsigned int j = 2; j <= this->degree; ++j)
1554  {
1555  interpolation_matrix(i + (j + 2) * source_fe.degree -
1556  j,
1557  i + (j + 2) * this->degree - j) =
1558  std::pow(0.5, i + j);
1559  factorial_k *= j;
1560  int factorial_l = factorial_k;
1561  int factorial_kl = 1;
1562 
1563  for (unsigned int k = j + 1; k <= this->degree; ++k)
1564  {
1565  factorial_kl *= k - j;
1566  factorial_l *= k;
1567 
1568  if (((j + k) & 1) != 0u)
1569  interpolation_matrix(
1570  i + (j + 2) * source_fe.degree - j,
1571  i + (k + 2) * this->degree - k) =
1572  -1.0 * std::pow(0.5, i + k) * factorial_l /
1573  (factorial_k * factorial_kl);
1574 
1575  else
1576  interpolation_matrix(
1577  i + (j + 2) * source_fe.degree - j,
1578  i + (k + 2) * this->degree - k) =
1579  std::pow(0.5, i + k) * factorial_l /
1580  (factorial_k * factorial_kl);
1581  }
1582  }
1583  }
1584 
1585  break;
1586  }
1587 
1588  case 2:
1589  {
1590  interpolation_matrix(0, 0) = 0.5;
1591  interpolation_matrix(0, 2) = 0.5;
1592  interpolation_matrix(2, 2) = 1.0;
1593  interpolation_matrix(3, 2) = 0.5;
1594  interpolation_matrix(3, 3) = 0.5;
1595 
1596  for (unsigned int i = 0; i < this->degree - 1;)
1597  {
1598  for (unsigned int line = 0;
1599  line < GeometryInfo<3>::lines_per_face;
1600  ++line)
1601  interpolation_matrix(1,
1602  i + line * (this->degree - 1) +
1603  4) = -0.5;
1604 
1605  for (unsigned int j = 0; j < this->degree - 1;)
1606  {
1607  interpolation_matrix(1,
1608  i + (j + 4) * this->degree - j) =
1609  1.0;
1610  j = j + 2;
1611  }
1612 
1613  interpolation_matrix(0, i + 4) = -1.0;
1614  interpolation_matrix(3, i + 3 * this->degree + 1) = -1.0;
1615  i = i + 2;
1616  }
1617 
1618  for (unsigned int vertex = 0;
1619  vertex < GeometryInfo<3>::vertices_per_face;
1620  ++vertex)
1621  interpolation_matrix(1, vertex) = 0.25;
1622 
1623  int factorial_i = 1;
1624 
1625  for (unsigned int i = 2; i <= this->degree; ++i)
1626  {
1627  double tmp = std::pow(0.5, i);
1628  interpolation_matrix(i + 2, i + 2) = tmp;
1629  interpolation_matrix(i + 3 * source_fe.degree - 1,
1630  i + 3 * this->degree - 1) = tmp;
1631  tmp *= 0.5;
1632  interpolation_matrix(i + source_fe.degree + 1, i + 2) =
1633  tmp;
1634  interpolation_matrix(i + source_fe.degree + 1,
1635  i + this->degree + 1) = tmp;
1636  interpolation_matrix(i + 2 * source_fe.degree,
1637  i + 2 * this->degree) = tmp;
1638  interpolation_matrix(i + 2 * source_fe.degree,
1639  i + 3 * this->degree - 1) = tmp;
1640  tmp *= -2.0;
1641 
1642  for (unsigned int j = 0; j < this->degree - 1;)
1643  {
1644  interpolation_matrix(i + source_fe.degree + 1,
1645  j + (i + 2) * this->degree + 2 -
1646  i) = tmp;
1647  interpolation_matrix(i + 2 * source_fe.degree,
1648  i + (j + 4) * this->degree - j -
1649  2) = tmp;
1650  j = j + 2;
1651  }
1652 
1653  int factorial_k = 1;
1654 
1655  for (unsigned int j = 2; j <= this->degree; ++j)
1656  {
1657  interpolation_matrix(i + (j + 2) * source_fe.degree -
1658  j,
1659  i + (j + 2) * this->degree - j) =
1660  std::pow(0.5, i + j);
1661  factorial_k *= j;
1662  int factorial_kl = 1;
1663  int factorial_l = factorial_k;
1664 
1665  for (unsigned int k = j + 1; k <= this->degree; ++k)
1666  {
1667  factorial_kl *= k - j;
1668  factorial_l *= k;
1669  interpolation_matrix(
1670  i + (j + 2) * source_fe.degree - j,
1671  i + (k + 2) * this->degree - k) =
1672  std::pow(0.5, i + k) * factorial_l /
1673  (factorial_k * factorial_kl);
1674  }
1675  }
1676 
1677  factorial_i *= i;
1678  int factorial_j = factorial_i;
1679  int factorial_ij = 1;
1680 
1681  for (unsigned int j = i + 1; j <= this->degree; ++j)
1682  {
1683  factorial_ij *= j - i;
1684  factorial_j *= j;
1685  tmp = std::pow(0.5, j) * factorial_j /
1686  (factorial_i * factorial_ij);
1687  interpolation_matrix(i + 2, j + 2) = tmp;
1688  tmp *= -1.0;
1689 
1690  for (unsigned int k = 0; k < this->degree - 1;)
1691  {
1692  interpolation_matrix(i + source_fe.degree + 1,
1693  k + (j + 2) * this->degree +
1694  2 - j) = tmp;
1695  k = k + 2;
1696  }
1697 
1698  tmp *= -0.5;
1699  interpolation_matrix(i + source_fe.degree + 1,
1700  j + 2) = tmp;
1701  interpolation_matrix(i + source_fe.degree + 1,
1702  j + this->degree + 1) = tmp;
1703 
1704  if (((i + j) & 1) != 0u)
1705  tmp *= -1.0;
1706 
1707  interpolation_matrix(i + 2 * source_fe.degree,
1708  j + 2 * this->degree) = tmp;
1709  interpolation_matrix(i + 2 * source_fe.degree,
1710  j + 3 * this->degree - 1) = tmp;
1711  tmp *= 2.0;
1712  interpolation_matrix(i + 3 * source_fe.degree - 1,
1713  j + 3 * this->degree - 1) = tmp;
1714  int factorial_k = 1;
1715 
1716  for (unsigned int k = 2; k <= this->degree; ++k)
1717  {
1718  interpolation_matrix(
1719  i + (k + 2) * source_fe.degree - k,
1720  j + (k + 2) * this->degree - k) =
1721  tmp * std::pow(0.5, k);
1722  factorial_k *= k;
1723  int factorial_l = factorial_k;
1724  int factorial_kl = 1;
1725 
1726  for (unsigned int l = k + 1; l <= this->degree;
1727  ++l)
1728  {
1729  factorial_kl *= l - k;
1730  factorial_l *= l;
1731  interpolation_matrix(
1732  i + (k + 2) * source_fe.degree - k,
1733  j + (l + 2) * this->degree - l) =
1734  tmp * std::pow(0.5, l) * factorial_l /
1735  (factorial_k * factorial_kl);
1736  }
1737  }
1738 
1739  tmp *= -1.0;
1740 
1741  for (unsigned int k = 0; k < this->degree - 1;)
1742  {
1743  interpolation_matrix(i + 2 * source_fe.degree,
1744  j + (k + 4) * this->degree -
1745  k - 2) = tmp;
1746  k = k + 2;
1747  }
1748  }
1749  }
1750 
1751  break;
1752  }
1753 
1754  case 3:
1755  {
1756  for (unsigned int vertex = 0;
1757  vertex < GeometryInfo<3>::vertices_per_face;
1758  ++vertex)
1759  interpolation_matrix(0, vertex) = 0.25;
1760 
1761  for (unsigned int i = 0; i < this->degree - 1;)
1762  {
1763  for (unsigned int line = 0;
1764  line < GeometryInfo<3>::lines_per_face;
1765  ++line)
1766  interpolation_matrix(0,
1767  i + line * (this->degree - 1) +
1768  4) = -0.5;
1769 
1770  for (unsigned int j = 0; j < this->degree - 1;)
1771  {
1772  interpolation_matrix(0,
1773  i + (j + 4) * this->degree - j) =
1774  1.0;
1775  j = j + 2;
1776  }
1777 
1778  interpolation_matrix(1, i + 4) = -1.0;
1779  interpolation_matrix(2, i + 3 * this->degree + 1) = -1.0;
1780  i = i + 2;
1781  }
1782 
1783  interpolation_matrix(1, 0) = 0.5;
1784  interpolation_matrix(1, 1) = 0.5;
1785  interpolation_matrix(2, 2) = 0.5;
1786  interpolation_matrix(2, 3) = 0.5;
1787  interpolation_matrix(3, 3) = 1.0;
1788 
1789  int factorial_i = 1;
1790 
1791  for (unsigned int i = 2; i <= this->degree; ++i)
1792  {
1793  double tmp = std::pow(0.5, i + 1);
1794  interpolation_matrix(i + 2, i + 2) = tmp;
1795  interpolation_matrix(i + 2, i + this->degree + 1) = tmp;
1796  interpolation_matrix(i + 2 * source_fe.degree,
1797  i + 2 * this->degree) = tmp;
1798  interpolation_matrix(i + 2 * source_fe.degree,
1799  i + 3 * this->degree - 1) = tmp;
1800  tmp *= -2.0;
1801 
1802  for (unsigned int j = 0; j < this->degree - 1;)
1803  {
1804  interpolation_matrix(i + 2,
1805  j + (i + 2) * this->degree + 2 -
1806  i) = tmp;
1807  interpolation_matrix(i + 2 * source_fe.degree,
1808  i + (j + 4) * this->degree - 2) =
1809  tmp;
1810  j = j + 2;
1811  }
1812 
1813  tmp *= -1.0;
1814  interpolation_matrix(i + source_fe.degree + 1,
1815  i + this->degree + 1) = tmp;
1816  interpolation_matrix(i + 3 * source_fe.degree - 1,
1817  i + 3 * this->degree - 1) = tmp;
1818  int factorial_k = 1;
1819 
1820  for (unsigned int j = 2; j <= this->degree; ++j)
1821  {
1822  interpolation_matrix(i + (j + 2) * source_fe.degree -
1823  j,
1824  i + (j + 2) * this->degree - j) =
1825  std::pow(0.5, i + j);
1826  factorial_k *= j;
1827  int factorial_l = factorial_k;
1828  int factorial_kl = 1;
1829 
1830  for (unsigned int k = j + 1; k <= this->degree; ++k)
1831  {
1832  factorial_kl *= k - j;
1833  factorial_l *= k;
1834  interpolation_matrix(
1835  i + (j + 2) * source_fe.degree - j,
1836  i + (k + 2) * this->degree - k) =
1837  std::pow(0.5, i + k) * factorial_l /
1838  (factorial_k * factorial_kl);
1839  }
1840  }
1841 
1842  factorial_i *= i;
1843  int factorial_j = factorial_i;
1844  int factorial_ij = 1;
1845 
1846  for (unsigned int j = i + 1; j <= this->degree; ++j)
1847  {
1848  factorial_ij *= j - i;
1849  factorial_j *= j;
1850  tmp = std::pow(0.5, j + 1) * factorial_j /
1851  (factorial_i * factorial_ij);
1852  interpolation_matrix(i + 2, j + 2) = tmp;
1853  interpolation_matrix(i + 2, j + this->degree + 1) =
1854  tmp;
1855  interpolation_matrix(i + 2 * source_fe.degree,
1856  j + 2 * this->degree) = tmp;
1857  interpolation_matrix(i + 2 * source_fe.degree,
1858  j + 3 * this->degree - 1) = tmp;
1859  tmp *= 2.0;
1860  interpolation_matrix(i + source_fe.degree + 1,
1861  j + this->degree + 1) = tmp;
1862  interpolation_matrix(i + 3 * source_fe.degree - 1,
1863  j + 3 * this->degree - 1) = tmp;
1864  int factorial_k = 1;
1865 
1866  for (unsigned int k = 2; k <= this->degree; ++k)
1867  {
1868  interpolation_matrix(
1869  i + (k + 2) * source_fe.degree - k,
1870  j + (k + 2) * this->degree - k) =
1871  tmp * std::pow(0.5, k);
1872  factorial_k *= k;
1873  int factorial_l = factorial_k;
1874  int factorial_kl = 1;
1875 
1876  for (unsigned int l = k + 1; l <= this->degree;
1877  ++l)
1878  {
1879  factorial_kl *= l - k;
1880  factorial_l *= l;
1881  interpolation_matrix(
1882  i + (k + 2) * source_fe.degree - k,
1883  j + (l + 2) * this->degree - l) =
1884  tmp * std::pow(0.5, l) * factorial_l /
1885  (factorial_k * factorial_kl);
1886  }
1887  }
1888 
1889  tmp *= -1.0;
1890 
1891  for (unsigned int k = 0; k < this->degree - 1;)
1892  {
1893  interpolation_matrix(i + 2,
1894  k + (j + 2) * this->degree +
1895  2 - j) = tmp;
1896  interpolation_matrix(i + 2 * source_fe.degree,
1897  j + (k + 4) * this->degree -
1898  2) = tmp;
1899  k = k + 2;
1900  }
1901  }
1902  }
1903  }
1904  }
1905  }
1906  }
1907 }
1908 
1909 
1910 
1911 template <int dim>
1912 void
1914 {
1915  const unsigned int codim = dim - 1;
1916 
1917  // TODO: the implementation makes the assumption that all faces have the
1918  // same number of dofs
1919  AssertDimension(this->n_unique_faces(), 1);
1920  const unsigned int face_no = 0;
1921 
1922  // number of points: (degree+1)^codim
1923  unsigned int n = this->degree + 1;
1924  for (unsigned int i = 1; i < codim; ++i)
1925  n *= this->degree + 1;
1926 
1927  this->generalized_face_support_points[face_no].resize(n);
1928 
1929  Point<codim> p;
1930 
1931  unsigned int k = 0;
1932  for (unsigned int iz = 0; iz <= ((codim > 2) ? this->degree : 0); ++iz)
1933  for (unsigned int iy = 0; iy <= ((codim > 1) ? this->degree : 0); ++iy)
1934  for (unsigned int ix = 0; ix <= this->degree; ++ix)
1935  {
1936  if (ix == 0)
1937  p(0) = 0.;
1938  else if (ix == 1)
1939  p(0) = 1.;
1940  else
1941  p(0) = .5;
1942  if (codim > 1)
1943  {
1944  if (iy == 0)
1945  p(1) = 0.;
1946  else if (iy == 1)
1947  p(1) = 1.;
1948  else
1949  p(1) = .5;
1950  }
1951  if (codim > 2)
1952  {
1953  if (iz == 0)
1954  p(2) = 0.;
1955  else if (iz == 1)
1956  p(2) = 1.;
1957  else
1958  p(2) = .5;
1959  }
1960  this->generalized_face_support_points[face_no][face_renumber[k++]] =
1961  p;
1962  }
1963 }
1964 
1965 
1966 // we use same dpo_vector as FE_Q
1967 template <int dim>
1968 std::vector<unsigned int>
1970 {
1971  std::vector<unsigned int> dpo(dim + 1, 1U);
1972  for (unsigned int i = 1; i < dpo.size(); ++i)
1973  dpo[i] = dpo[i - 1] * (deg - 1);
1974  return dpo;
1975 }
1976 
1977 
1978 
1979 template <int dim>
1980 std::vector<unsigned int>
1982  const FiniteElementData<dim> &fe)
1983 {
1984  Assert(fe.n_components() == 1, ExcInternalError());
1985  std::vector<unsigned int> h2l(fe.n_dofs_per_cell());
1986 
1987  // polynomial degree
1988  const unsigned int degree = fe.n_dofs_per_line() + 1;
1989  // number of grid points in each
1990  // direction
1991  const unsigned int n = degree + 1;
1992 
1993  // the following lines of code are
1994  // somewhat odd, due to the way the
1995  // hierarchic numbering is
1996  // organized. if someone would
1997  // really want to understand these
1998  // lines, you better draw some
1999  // pictures where you indicate the
2000  // indices and orders of vertices,
2001  // lines, etc, along with the
2002  // numbers of the degrees of
2003  // freedom in hierarchical and
2004  // lexicographical order
2005  switch (dim)
2006  {
2007  case 1:
2008  {
2009  for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
2010  h2l[i] = i;
2011 
2012  break;
2013  }
2014 
2015  case 2:
2016  {
2017  // Example: degree=3
2018  //
2019  // hierarchical numbering:
2020  // 2 10 11 3
2021  // 5 14 15 7
2022  // 4 12 13 6
2023  // 0 8 9 1
2024  //
2025  // fe_q_hierarchical numbering:
2026  // 4 6 7 5
2027  // 12 14 15 13
2028  // 8 10 11 9
2029  // 0 2 3 1
2030  unsigned int next_index = 0;
2031  // first the four vertices
2032  h2l[next_index++] = 0;
2033  h2l[next_index++] = 1;
2034  h2l[next_index++] = n;
2035  h2l[next_index++] = n + 1;
2036  // left line
2037  for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2038  h2l[next_index++] = (2 + i) * n;
2039  // right line
2040  for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2041  h2l[next_index++] = (2 + i) * n + 1;
2042  // bottom line
2043  for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2044  h2l[next_index++] = 2 + i;
2045  // top line
2046  for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2047  h2l[next_index++] = n + 2 + i;
2048  // inside quad
2049  Assert(fe.n_dofs_per_quad(0 /*only one quad in 2D*/) ==
2050  fe.n_dofs_per_line() * fe.n_dofs_per_line(),
2051  ExcInternalError());
2052  for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2053  for (unsigned int j = 0; j < fe.n_dofs_per_line(); ++j)
2054  h2l[next_index++] = (2 + i) * n + 2 + j;
2055 
2056  Assert(next_index == fe.n_dofs_per_cell(), ExcInternalError());
2057 
2058  break;
2059  }
2060 
2061  case 3:
2062  {
2063  unsigned int next_index = 0;
2064  const unsigned int n2 = n * n;
2065  // first the eight vertices
2066  // bottom face, lexicographic
2067  h2l[next_index++] = 0;
2068  h2l[next_index++] = 1;
2069  h2l[next_index++] = n;
2070  h2l[next_index++] = n + 1;
2071  // top face, lexicographic
2072  h2l[next_index++] = n2;
2073  h2l[next_index++] = n2 + 1;
2074  h2l[next_index++] = n2 + n;
2075  h2l[next_index++] = n2 + n + 1;
2076 
2077  // now the lines
2078  // bottom face
2079  for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2080  h2l[next_index++] = (2 + i) * n;
2081  for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2082  h2l[next_index++] = (2 + i) * n + 1;
2083  for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2084  h2l[next_index++] = 2 + i;
2085  for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2086  h2l[next_index++] = n + 2 + i;
2087  // top face
2088  for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2089  h2l[next_index++] = n2 + (2 + i) * n;
2090  for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2091  h2l[next_index++] = n2 + (2 + i) * n + 1;
2092  for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2093  h2l[next_index++] = n2 + 2 + i;
2094  for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2095  h2l[next_index++] = n2 + n + 2 + i;
2096  // lines in z-direction
2097  for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2098  h2l[next_index++] = (2 + i) * n2;
2099  for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2100  h2l[next_index++] = (2 + i) * n2 + 1;
2101  for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2102  h2l[next_index++] = (2 + i) * n2 + n;
2103  for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2104  h2l[next_index++] = (2 + i) * n2 + n + 1;
2105 
2106  // TODO: the implementation makes the assumption that all faces have
2107  // the same number of dofs
2109  const unsigned int face_no = 0;
2110  (void)face_no;
2111 
2112  // inside quads
2113  Assert(fe.n_dofs_per_quad(face_no) ==
2114  fe.n_dofs_per_line() * fe.n_dofs_per_line(),
2115  ExcInternalError());
2116  // left face
2117  for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2118  for (unsigned int j = 0; j < fe.n_dofs_per_line(); ++j)
2119  h2l[next_index++] = (2 + i) * n2 + (2 + j) * n;
2120  // right face
2121  for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2122  for (unsigned int j = 0; j < fe.n_dofs_per_line(); ++j)
2123  h2l[next_index++] = (2 + i) * n2 + (2 + j) * n + 1;
2124  // front face
2125  for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2126  for (unsigned int j = 0; j < fe.n_dofs_per_line(); ++j)
2127  h2l[next_index++] = (2 + i) * n2 + 2 + j;
2128  // back face
2129  for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2130  for (unsigned int j = 0; j < fe.n_dofs_per_line(); ++j)
2131  h2l[next_index++] = (2 + i) * n2 + n + 2 + j;
2132  // bottom face
2133  for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2134  for (unsigned int j = 0; j < fe.n_dofs_per_line(); ++j)
2135  h2l[next_index++] = (2 + i) * n + 2 + j;
2136  // top face
2137  for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2138  for (unsigned int j = 0; j < fe.n_dofs_per_line(); ++j)
2139  h2l[next_index++] = n2 + (2 + i) * n + 2 + j;
2140 
2141  // inside hex
2142  Assert(fe.n_dofs_per_hex() ==
2143  fe.n_dofs_per_quad(face_no) * fe.n_dofs_per_line(),
2144  ExcInternalError());
2145  for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2146  for (unsigned int j = 0; j < fe.n_dofs_per_line(); ++j)
2147  for (unsigned int k = 0; k < fe.n_dofs_per_line(); ++k)
2148  h2l[next_index++] = (2 + i) * n2 + (2 + j) * n + 2 + k;
2149 
2150  Assert(next_index == fe.n_dofs_per_cell(), ExcInternalError());
2151 
2152  break;
2153  }
2154 
2155  default:
2156  Assert(false, ExcNotImplemented());
2157  }
2158  return h2l;
2159 }
2160 
2161 
2162 template <int dim>
2163 std::vector<unsigned int>
2165  const unsigned int degree)
2166 {
2167  FiniteElementData<dim - 1> fe_data(
2168  FE_Q_Hierarchical<dim - 1>::get_dpo_vector(degree), 1, degree);
2169  return internal::FE_Q_Hierarchical::invert_numbering(
2171  fe_data));
2172 }
2173 
2174 
2175 
2176 template <>
2177 std::vector<unsigned int>
2179  const unsigned int)
2180 {
2181  return std::vector<unsigned int>();
2182 }
2183 
2184 
2185 template <>
2186 bool
2187 FE_Q_Hierarchical<1>::has_support_on_face(const unsigned int shape_index,
2188  const unsigned int face_index) const
2189 {
2190  AssertIndexRange(shape_index, this->n_dofs_per_cell());
2192 
2193 
2194  // in 1d, things are simple. since
2195  // there is only one degree of
2196  // freedom per vertex in this
2197  // class, the first is on vertex 0
2198  // (==face 0 in some sense), the
2199  // second on face 1:
2200  return (((shape_index == 0) && (face_index == 0)) ||
2201  ((shape_index == 1) && (face_index == 1)));
2202 }
2203 
2204 
2205 
2206 template <int dim>
2207 bool
2208 FE_Q_Hierarchical<dim>::has_support_on_face(const unsigned int shape_index,
2209  const unsigned int face_index) const
2210 {
2211  AssertIndexRange(shape_index, this->n_dofs_per_cell());
2213 
2214  // first, special-case interior
2215  // shape functions, since they
2216  // have no support no-where on
2217  // the boundary
2218  if (((dim == 2) && (shape_index >=
2219  this->get_first_quad_index(0 /*only one quad in 2D*/))) ||
2220  ((dim == 3) && (shape_index >= this->get_first_hex_index())))
2221  return false;
2222 
2223  // let's see whether this is a
2224  // vertex
2225  if (shape_index < this->get_first_line_index())
2226  {
2227  // for Q elements, there is
2228  // one dof per vertex, so
2229  // shape_index==vertex_number. check
2230  // whether this vertex is
2231  // on the given face.
2232  const unsigned int vertex_no = shape_index;
2234  ExcInternalError());
2235  for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_face; ++i)
2236  if (GeometryInfo<dim>::face_to_cell_vertices(face_index, i) ==
2237  vertex_no)
2238  return true;
2239  return false;
2240  }
2241  else if (shape_index < this->get_first_quad_index(0))
2242  // ok, dof is on a line
2243  {
2244  const unsigned int line_index =
2245  (shape_index - this->get_first_line_index()) / this->n_dofs_per_line();
2247  ExcInternalError());
2248 
2249  for (unsigned int i = 0; i < GeometryInfo<dim>::lines_per_face; ++i)
2250  if (GeometryInfo<dim>::face_to_cell_lines(face_index, i) == line_index)
2251  return true;
2252  return false;
2253  }
2254  else if (shape_index < this->get_first_hex_index())
2255  // dof is on a quad
2256  {
2257  const unsigned int quad_index =
2258  (shape_index - this->get_first_quad_index(0 /*first quad*/)) /
2259  this->n_dofs_per_quad(face_index);
2260  Assert(static_cast<signed int>(quad_index) <
2261  static_cast<signed int>(GeometryInfo<dim>::quads_per_cell),
2262  ExcInternalError());
2263 
2264  // in 2d, cell bubble are
2265  // zero on all faces. but
2266  // we have treated this
2267  // case above already
2268  Assert(dim != 2, ExcInternalError());
2269 
2270  // in 3d,
2271  // quad_index=face_index
2272  if (dim == 3)
2273  return (quad_index == face_index);
2274  else
2275  Assert(false, ExcNotImplemented());
2276  }
2277  else
2278  // dof on hex
2279  {
2280  // can only happen in 3d, but
2281  // this case has already been
2282  // covered above
2283  Assert(false, ExcNotImplemented());
2284  return false;
2285  }
2286 
2287  // we should not have gotten here
2288  Assert(false, ExcInternalError());
2289  return false;
2290 }
2291 
2292 
2293 
2294 template <int dim>
2295 std::vector<unsigned int>
2296 FE_Q_Hierarchical<dim>::get_embedding_dofs(const unsigned int sub_degree) const
2297 {
2298  Assert((sub_degree > 0) && (sub_degree <= this->degree),
2299  ExcIndexRange(sub_degree, 1, this->degree));
2300 
2301  if (dim == 1)
2302  {
2303  std::vector<unsigned int> embedding_dofs(sub_degree + 1);
2304  for (unsigned int i = 0; i < (sub_degree + 1); ++i)
2305  embedding_dofs[i] = i;
2306 
2307  return embedding_dofs;
2308  }
2309 
2310  if (sub_degree == 1)
2311  {
2312  std::vector<unsigned int> embedding_dofs(
2314  for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
2315  embedding_dofs[i] = i;
2316 
2317  return embedding_dofs;
2318  }
2319  else if (sub_degree == this->degree)
2320  {
2321  std::vector<unsigned int> embedding_dofs(this->n_dofs_per_cell());
2322  for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
2323  embedding_dofs[i] = i;
2324 
2325  return embedding_dofs;
2326  }
2327 
2328  if ((dim == 2) || (dim == 3))
2329  {
2330  std::vector<unsigned int> embedding_dofs(
2331  (dim == 2) ? (sub_degree + 1) * (sub_degree + 1) :
2332  (sub_degree + 1) * (sub_degree + 1) * (sub_degree + 1));
2333 
2334  for (unsigned int i = 0;
2335  i < ((dim == 2) ?
2336  (sub_degree + 1) * (sub_degree + 1) :
2337  (sub_degree + 1) * (sub_degree + 1) * (sub_degree + 1));
2338  ++i)
2339  {
2340  // vertex
2342  embedding_dofs[i] = i;
2343  // line
2344  else if (i < (GeometryInfo<dim>::vertices_per_cell +
2345  GeometryInfo<dim>::lines_per_cell * (sub_degree - 1)))
2346  {
2347  const unsigned int j =
2348  (i - GeometryInfo<dim>::vertices_per_cell) % (sub_degree - 1);
2349  const unsigned int line =
2351  (sub_degree - 1);
2352 
2353  embedding_dofs[i] = GeometryInfo<dim>::vertices_per_cell +
2354  line * (this->degree - 1) + j;
2355  }
2356  // quad
2357  else if (i < (GeometryInfo<dim>::vertices_per_cell +
2358  GeometryInfo<dim>::lines_per_cell * (sub_degree - 1)) +
2359  GeometryInfo<dim>::quads_per_cell * (sub_degree - 1) *
2360  (sub_degree - 1))
2361  {
2362  const unsigned int j =
2364  GeometryInfo<dim>::lines_per_cell * (sub_degree - 1)) %
2365  (sub_degree - 1);
2366  const unsigned int k =
2368  GeometryInfo<dim>::lines_per_cell * (sub_degree - 1) - j) /
2369  (sub_degree - 1)) %
2370  (sub_degree - 1);
2371  const unsigned int face =
2373  GeometryInfo<dim>::lines_per_cell * (sub_degree - 1) -
2374  k * (sub_degree - 1) - j) /
2375  ((sub_degree - 1) * (sub_degree - 1));
2376 
2377  embedding_dofs[i] =
2379  GeometryInfo<dim>::lines_per_cell * (this->degree - 1) +
2380  face * (this->degree - 1) * (this->degree - 1) +
2381  k * (this->degree - 1) + j;
2382  }
2383  // hex
2384  else if (i < (GeometryInfo<dim>::vertices_per_cell +
2385  GeometryInfo<dim>::lines_per_cell * (sub_degree - 1)) +
2386  GeometryInfo<dim>::quads_per_cell * (sub_degree - 1) *
2387  (sub_degree - 1) +
2388  GeometryInfo<dim>::hexes_per_cell * (sub_degree - 1) *
2389  (sub_degree - 1) * (sub_degree - 1))
2390  {
2391  const unsigned int j =
2393  GeometryInfo<dim>::lines_per_cell * (sub_degree - 1) -
2394  GeometryInfo<dim>::quads_per_cell * (sub_degree - 1) *
2395  (sub_degree - 1)) %
2396  (sub_degree - 1);
2397  const unsigned int k =
2399  GeometryInfo<dim>::lines_per_cell * (sub_degree - 1) -
2400  GeometryInfo<dim>::quads_per_cell * (sub_degree - 1) *
2401  (sub_degree - 1) -
2402  j) /
2403  (sub_degree - 1)) %
2404  (sub_degree - 1);
2405  const unsigned int l =
2407  GeometryInfo<dim>::lines_per_cell * (sub_degree - 1) -
2408  GeometryInfo<dim>::quads_per_cell * (sub_degree - 1) *
2409  (sub_degree - 1) -
2410  j - k * (sub_degree - 1)) /
2411  ((sub_degree - 1) * (sub_degree - 1));
2412 
2413  embedding_dofs[i] =
2415  GeometryInfo<dim>::lines_per_cell * (this->degree - 1) +
2416  GeometryInfo<dim>::quads_per_cell * (this->degree - 1) *
2417  (this->degree - 1) +
2418  l * (this->degree - 1) * (this->degree - 1) +
2419  k * (this->degree - 1) + j;
2420  }
2421  }
2422 
2423  return embedding_dofs;
2424  }
2425  else
2426  {
2427  Assert(false, ExcNotImplemented());
2428  return std::vector<unsigned int>();
2429  }
2430 }
2431 
2432 
2433 
2434 template <int dim>
2435 std::pair<Table<2, bool>, std::vector<unsigned int>>
2437 {
2438  Table<2, bool> constant_modes(1, this->n_dofs_per_cell());
2439  for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
2440  constant_modes(0, i) = true;
2441  for (unsigned int i = GeometryInfo<dim>::vertices_per_cell;
2442  i < this->n_dofs_per_cell();
2443  ++i)
2444  constant_modes(0, i) = false;
2445  return std::pair<Table<2, bool>, std::vector<unsigned int>>(
2446  constant_modes, std::vector<unsigned int>(1, 0));
2447 }
2448 
2449 
2450 
2451 template <int dim>
2452 std::size_t
2454 {
2455  Assert(false, ExcNotImplemented());
2456  return 0;
2457 }
2458 
2459 
2460 
2461 // explicit instantiations
2462 #include "fe_q_hierarchical.inst"
2463 
2464 
Definition: fe_dgq.h:113
const std::unique_ptr< ScalarPolynomialsBase< dim > > poly_space
Definition: fe_poly.h:536
virtual std::unique_ptr< FiniteElement< dim, dim > > clone() const override
virtual void get_interpolation_matrix(const FiniteElement< dim > &source, FullMatrix< double > &matrix) const override
virtual std::pair< Table< 2, bool >, std::vector< unsigned int > > get_constant_modes() const override
static std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
friend class FE_Q_Hierarchical
std::vector< unsigned int > get_embedding_dofs(const unsigned int sub_degree) const
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_quad_dof_identities(const FiniteElement< dim > &fe_other, const unsigned int face_no=0) const override
void build_dofs_cell(std::vector< FullMatrix< double >> &dofs_cell, std::vector< FullMatrix< double >> &dofs_subcell) const
void initialize_embedding_and_restriction(const std::vector< FullMatrix< double >> &dofs_cell, const std::vector< FullMatrix< double >> &dofs_subcell)
void initialize_generalized_support_points()
virtual FiniteElementDomination::Domination compare_for_domination(const FiniteElement< dim > &fe_other, const unsigned int codim=0) const override final
static std::vector< unsigned int > hierarchic_to_fe_q_hierarchical_numbering(const FiniteElementData< dim > &fe)
virtual const FullMatrix< double > & get_prolongation_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const override
virtual void get_subface_interpolation_matrix(const FiniteElement< dim > &source, const unsigned int subface, 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 > &fe_other) const override
static std::vector< unsigned int > face_fe_q_hierarchical_to_hierarchic_numbering(const unsigned int degree)
void initialize_generalized_face_support_points()
virtual std::string get_name() const override
virtual void get_face_interpolation_matrix(const FiniteElement< dim > &source, FullMatrix< double > &matrix, const unsigned int face_no=0) const override
virtual std::size_t memory_consumption() const override
void initialize_constraints(const std::vector< FullMatrix< double >> &dofs_subcell)
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const override
virtual bool hp_constraints_are_implemented() const override
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_line_dof_identities(const FiniteElement< dim > &fe_other) const override
unsigned int n_dofs_per_vertex() const
const unsigned int degree
Definition: fe_data.h:449
unsigned int n_dofs_per_cell() const
unsigned int n_dofs_per_line() const
unsigned int n_dofs_per_face(unsigned int face_no=0, unsigned int child=0) const
unsigned int n_components() const
unsigned int n_unique_faces() const
unsigned int n_dofs_per_quad(unsigned int face_no=0) const
unsigned int n_dofs_per_hex() const
virtual std::string get_name() const =0
size_type n() const
size_type m() const
Definition: point.h:111
const std::vector< unsigned int > & get_numbering_inverse() const
const std::vector< unsigned int > & get_numbering() const
void set_numbering(const std::vector< unsigned int > &renumber)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
static ::ExceptionBase & ExcIndexRange(std::size_t arg1, std::size_t arg2, std::size_t arg3)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
@ matrix
Contents is actually a matrix.
static const char U
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)