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