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