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