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