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