Reference documentation for deal.II version 8.4.2
fe_q_hierarchical.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2002 - 2015 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_unit_support_points ();
103  initialize_unit_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_from_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  const unsigned int degree=this->degree;
506 
507  this->interface_constraints
508  .TableBase<2,double>::reinit (this->interface_constraints_size());
509 
510  switch (dim)
511  {
512  case 1:
513  {
514  // no constraints in 1d
515  break;
516  }
517 
518  case 2:
519  {
520  // vertex node
521  for (unsigned int i=0; i<dofs_1d; ++i)
522  this->interface_constraints(0,i) = dofs_subcell[0](1,i);
523  // edge nodes
524  for (unsigned int c=0; c<GeometryInfo<1>::max_children_per_cell; ++c)
525  for (unsigned int i=0; i<dofs_1d; ++i)
526  for (unsigned int j=2; j<dofs_1d; ++j)
527  this->interface_constraints(1 + c*(degree-1) + j - 2,i) =
528  dofs_subcell[c](j,i);
529  break;
530  }
531 
532  case 3:
533  {
534  for (unsigned int i=0; i<dofs_1d * dofs_1d; i++)
535  {
536  // center vertex node
537  this->interface_constraints(0,face_renumber[i]) =
538  dofs_subcell[0](1,i % dofs_1d) *
539  dofs_subcell[0](1,(i - (i % dofs_1d)) / dofs_1d);
540 
541  // boundary vertex nodes
542  this->interface_constraints(1,face_renumber[i]) =
543  dofs_subcell[0](0, i % dofs_1d) *
544  dofs_subcell[0](1, (i - (i % dofs_1d)) / dofs_1d);
545  this->interface_constraints(2,face_renumber[i]) =
546  dofs_subcell[1](1, i % dofs_1d) *
547  dofs_subcell[0](1, (i - (i % dofs_1d)) / dofs_1d);
548  this->interface_constraints(3,face_renumber[i]) =
549  dofs_subcell[0](1, i % dofs_1d) *
550  dofs_subcell[0](0, (i - (i % dofs_1d)) / dofs_1d);
551  this->interface_constraints(4,face_renumber[i]) =
552  dofs_subcell[1](0, i % dofs_1d) *
553  dofs_subcell[1](1, (i - (i % dofs_1d)) / dofs_1d);
554 
555  // interior edges
556  for (unsigned int j=0; j<(degree-1); j++)
557  {
558  this->interface_constraints(5 + j,face_renumber[i]) =
559  dofs_subcell[0](1, i % dofs_1d) *
560  dofs_subcell[0](2 + j, (i - (i % dofs_1d)) / dofs_1d);
561  this->interface_constraints(5 + (degree-1) + j,face_renumber[i]) =
562  dofs_subcell[0](1,i % dofs_1d) *
563  dofs_subcell[1](2 + j, (i - (i % dofs_1d)) / dofs_1d);
564  this->interface_constraints(5 + 2*(degree-1) + j,face_renumber[i]) =
565  dofs_subcell[0](2 + j,i % dofs_1d) *
566  dofs_subcell[1](0, (i - (i % dofs_1d)) / dofs_1d);
567  this->interface_constraints(5 + 3*(degree-1) + j,face_renumber[i]) =
568  dofs_subcell[1](2 + j, i % dofs_1d) *
569  dofs_subcell[0](1, (i - (i % dofs_1d)) / dofs_1d);
570  }
571 
572  // boundary edges
573  for (unsigned int j=0; j<(degree-1); j++)
574  {
575  // left edge
576  this->interface_constraints(5 + 4*(degree-1) + j,face_renumber[i]) =
577  dofs_subcell[0](0, i % dofs_1d) *
578  dofs_subcell[0](2 + j, (i - (i % dofs_1d)) / dofs_1d);
579  this->interface_constraints(5 + 4*(degree-1) + (degree-1) + j,face_renumber[i]) =
580  dofs_subcell[0](0, i % dofs_1d) *
581  dofs_subcell[1](2 + j, (i - (i % dofs_1d)) / dofs_1d);
582  // right edge
583  this->interface_constraints(5 + 4*(degree-1) + 2*(degree-1) + j,face_renumber[i]) =
584  dofs_subcell[1](1, i % dofs_1d) *
585  dofs_subcell[0](2 + j, (i - (i % dofs_1d)) / dofs_1d);
586  this->interface_constraints(5 + 4*(degree-1) + 3*(degree-1) + j,face_renumber[i]) =
587  dofs_subcell[1](1, i % dofs_1d) *
588  dofs_subcell[1](2 + j, (i - (i % dofs_1d)) / dofs_1d);
589  // bottom edge
590  this->interface_constraints(5 + 4*(degree-1) + 4*(degree-1) + j,face_renumber[i]) =
591  dofs_subcell[0](2 + j, i % dofs_1d) *
592  dofs_subcell[0](0, (i - (i % dofs_1d)) / dofs_1d);
593  this->interface_constraints(5 + 4*(degree-1) + 5*(degree-1) + j,face_renumber[i]) =
594  dofs_subcell[1](2 + j, i % dofs_1d) *
595  dofs_subcell[0](0, (i - (i % dofs_1d)) / dofs_1d);
596  // top edge
597  this->interface_constraints(5 + 4*(degree-1) + 6*(degree-1) + j,face_renumber[i]) =
598  dofs_subcell[0](2 + j, i % dofs_1d) *
599  dofs_subcell[1](1, (i - (i % dofs_1d)) / dofs_1d);
600  this->interface_constraints(5 + 4*(degree-1) + 7*(degree-1) + j,face_renumber[i]) =
601  dofs_subcell[1](2 + j, i % dofs_1d) *
602  dofs_subcell[1](1, (i - (i % dofs_1d)) / dofs_1d);
603  }
604 
605  // interior faces
606  for (unsigned int j=0; j<(degree-1); j++)
607  for (unsigned int k=0; k<(degree-1); k++)
608  {
609  // subcell 0
610  this->interface_constraints(5 + 12*(degree-1) + j + k*(degree-1),face_renumber[i]) =
611  dofs_subcell[0](2 + j, i % dofs_1d) *
612  dofs_subcell[0](2 + k, (i - (i % dofs_1d)) / dofs_1d);
613  // subcell 1
614  this->interface_constraints(5 + 12*(degree-1) + j + k*(degree-1) + (degree-1)*(degree-1),face_renumber[i]) =
615  dofs_subcell[1](2 + j, i % dofs_1d) *
616  dofs_subcell[0](2 + k, (i - (i % dofs_1d)) / dofs_1d);
617  // subcell 2
618  this->interface_constraints(5 + 12*(degree-1) + j + k*(degree-1) + 2*(degree-1)*(degree-1),face_renumber[i]) =
619  dofs_subcell[0](2 + j, i % dofs_1d) *
620  dofs_subcell[1](2 + k, (i - (i % dofs_1d)) / dofs_1d);
621  // subcell 3
622  this->interface_constraints(5 + 12*(degree-1) + j + k*(degree-1) + 3*(degree-1)*(degree-1),face_renumber[i]) =
623  dofs_subcell[1](2 + j, i % dofs_1d) *
624  dofs_subcell[1](2 + k, (i - (i % dofs_1d)) / dofs_1d);
625  }
626  }
627  break;
628  }
629 
630  default:
631  Assert (false, ExcNotImplemented());
632  }
633 }
634 
635 
636 
637 template <int dim>
638 void
639 FE_Q_Hierarchical<dim>::
640 initialize_embedding_and_restriction (const std::vector<FullMatrix<double> > &dofs_cell,
641  const std::vector<FullMatrix<double> > &dofs_subcell)
642 {
643  unsigned int iso=RefinementCase<dim>::isotropic_refinement-1;
644 
645  const unsigned int dofs_1d = 2*this->dofs_per_vertex + this->dofs_per_line;
646  const std::vector<unsigned int> &renumber=
647  this->poly_space.get_numbering();
648 
649  for (unsigned int c=0; c<GeometryInfo<dim>::max_children_per_cell; ++c)
650  {
651  this->prolongation[iso][c].reinit (this->dofs_per_cell, this->dofs_per_cell);
652  this->restriction[iso][c].reinit (this->dofs_per_cell, this->dofs_per_cell);
653  }
654 
655  // the 1d case is particularly
656  // simple, so special case it:
657  if (dim==1)
658  {
659  for (unsigned int c=0; c<GeometryInfo<dim>::max_children_per_cell; ++c)
660  {
661  this->prolongation[iso][c].fill (dofs_subcell[c]);
662  this->restriction[iso][c].fill (dofs_cell[c]);
663  }
664  return;
665  }
666 
667  // for higher dimensions, things
668  // are a little more tricky:
669 
670  // j loops over dofs in the
671  // subcell. These are the rows in
672  // the embedding matrix.
673  //
674  // i loops over the dofs in the
675  // master cell. These are the
676  // columns in the embedding matrix.
677  for (unsigned int j=0; j<this->dofs_per_cell; ++j)
678  for (unsigned int i=0; i<this->dofs_per_cell; ++i)
679  switch (dim)
680  {
681  case 2:
682  {
683  for (unsigned int c=0; c<GeometryInfo<2>::max_children_per_cell; ++c)
684  {
685  // left/right line: 0/1
686  const unsigned int c0 = c%2;
687  // bottom/top line: 0/1
688  const unsigned int c1 = c/2;
689 
690  this->prolongation[iso][c](j,i) =
691  dofs_subcell[c0](renumber[j] % dofs_1d,
692  renumber[i] % dofs_1d) *
693  dofs_subcell[c1]((renumber[j] - (renumber[j] % dofs_1d)) / dofs_1d,
694  (renumber[i] - (renumber[i] % dofs_1d)) / dofs_1d);
695 
696  this->restriction[iso][c](j,i) =
697  dofs_cell[c0](renumber[j] % dofs_1d,
698  renumber[i] % dofs_1d) *
699  dofs_cell[c1]((renumber[j] - (renumber[j] % dofs_1d)) / dofs_1d,
700  (renumber[i] - (renumber[i] % dofs_1d)) / dofs_1d);
701  }
702  break;
703  }
704 
705  case 3:
706  {
707  for (unsigned int c=0; c<GeometryInfo<3>::max_children_per_cell; ++c)
708  {
709  // left/right face: 0/1
710  const unsigned int c0 = c%2;
711  // front/back face: 0/1
712  const unsigned int c1 = (c%4)/2;
713  // bottom/top face: 0/1
714  const unsigned int c2 = c/4;
715 
716  this->prolongation[iso][c](j,i) =
717  dofs_subcell[c0](renumber[j] % dofs_1d,
718  renumber[i] % dofs_1d) *
719  dofs_subcell[c1](((renumber[j] - (renumber[j] % dofs_1d)) / dofs_1d) % dofs_1d,
720  ((renumber[i] - (renumber[i] % dofs_1d)) / dofs_1d) % dofs_1d) *
721  dofs_subcell[c2](((renumber[j] - (renumber[j] % dofs_1d)) / dofs_1d - (((renumber[j] - (renumber[j] % dofs_1d)) / dofs_1d ) % dofs_1d)) / dofs_1d,
722  ((renumber[i] - (renumber[i] % dofs_1d)) / dofs_1d - (((renumber[i] - (renumber[i] % dofs_1d)) / dofs_1d ) % dofs_1d)) / dofs_1d);
723 
724  this->restriction[iso][c](j,i) =
725  dofs_cell[c0](renumber[j] % dofs_1d,
726  renumber[i] % dofs_1d) *
727  dofs_cell[c1](((renumber[j] - (renumber[j] % dofs_1d)) / dofs_1d) % dofs_1d,
728  ((renumber[i] - (renumber[i] % dofs_1d)) / dofs_1d) % dofs_1d) *
729  dofs_cell[c2](((renumber[j] - (renumber[j] % dofs_1d)) / dofs_1d - (((renumber[j] - (renumber[j] % dofs_1d)) / dofs_1d ) % dofs_1d)) / dofs_1d,
730  ((renumber[i] - (renumber[i] % dofs_1d)) / dofs_1d - (((renumber[i] - (renumber[i] % dofs_1d)) / dofs_1d ) % dofs_1d)) / dofs_1d);
731  }
732  break;
733  }
734 
735  default:
736  Assert (false, ExcNotImplemented());
737  }
738 }
739 
740 
741 
742 template <int dim>
743 void FE_Q_Hierarchical<dim>::initialize_unit_support_points ()
744 {
745  // number of points: (degree+1)^dim
746  unsigned int n = this->degree+1;
747  for (unsigned int i=1; i<dim; ++i)
748  n *= this->degree+1;
749 
750  this->unit_support_points.resize(n);
751 
752  const std::vector<unsigned int> &index_map_inverse=
753  this->poly_space.get_numbering_inverse();
754 
755  Point<dim> p;
756  // the method of numbering allows
757  // each dof to be associated with a
758  // support point. There is
759  // only one support point per
760  // vertex, line, quad, hex, etc.
761  //
762  // note, however, that the support
763  // points thus associated with
764  // shape functions are not unique:
765  // the linear shape functions are
766  // associated with the vertices,
767  // but all others are associated
768  // with either line, quad, or hex
769  // midpoints, and there may be
770  // multiple shape functions
771  // associated with them. there
772  // really is no other useful
773  // numbering, since the
774  // hierarchical shape functions do
775  // not vanish at all-but-one
776  // interpolation points (like the
777  // Lagrange functions used in
778  // FE_Q), so there's not much we
779  // can do here.
780 
781  // TODO shouldn't we just at least make support points unique,
782  // even though the delta property is not satisfied for this FE?
783  unsigned int k=0;
784  for (unsigned int iz=0; iz <= ((dim>2) ? this->degree : 0) ; ++iz)
785  for (unsigned int iy=0; iy <= ((dim>1) ? this->degree : 0) ; ++iy)
786  for (unsigned int ix=0; ix<=this->degree; ++ix)
787  {
788  if (ix==0)
789  p(0) = 0.;
790  else if (ix==1)
791  p(0) = 1.;
792  else
793  p(0) = .5;
794  if (dim>1)
795  {
796  if (iy==0)
797  p(1) = 0.;
798  else if (iy==1)
799  p(1) = 1.;
800  else
801  p(1) = .5;
802  }
803  if (dim>2)
804  {
805  if (iz==0)
806  p(2) = 0.;
807  else if (iz==1)
808  p(2) = 1.;
809  else
810  p(2) = .5;
811  }
812  this->unit_support_points[index_map_inverse[k++]] = p;
813  };
814 }
815 
816 
817 
818 template <>
819 void FE_Q_Hierarchical<1>::initialize_unit_face_support_points ()
820 {
821  // no faces in 1d, so nothing to do
822 }
823 
824 
825 template <>
826 void FE_Q_Hierarchical<1>::
827 get_face_interpolation_matrix (const FiniteElement<1,1> &/*x_source_fe*/,
828  FullMatrix<double> &/*interpolation_matrix*/) const
829 {
830  Assert (false, ExcImpossibleInDim(1));
831 }
832 
833 
834 template <>
835 void
836 FE_Q_Hierarchical<1>::
837 get_subface_interpolation_matrix (const FiniteElement<1,1> &/*x_source_fe*/,
838  const unsigned int /*subface*/,
839  FullMatrix<double> &/*interpolation_matrix*/) const
840 {
841  Assert (false, ExcImpossibleInDim(1));
842 }
843 
844 
845 
846 template <int dim>
847 void
848 FE_Q_Hierarchical<dim>::
849 get_face_interpolation_matrix (const FiniteElement<dim> &x_source_fe,
850  FullMatrix<double> &interpolation_matrix) const
851 {
852  // this is only implemented, if the
853  // source FE is also a
854  // Q_Hierarchical element
855  typedef FE_Q_Hierarchical<dim> FEQHierarchical;
856  typedef FiniteElement<dim> FEL;
857  AssertThrow ((x_source_fe.get_name().find ("FE_Q_Hierarchical<") == 0)
858  ||
859  (dynamic_cast<const FEQHierarchical *>(&x_source_fe) != 0),
860  typename FEL::
861  ExcInterpolationNotImplemented());
862 
863  Assert (interpolation_matrix.n() == this->dofs_per_face,
864  ExcDimensionMismatch (interpolation_matrix.n(),
865  this->dofs_per_face));
866  Assert (interpolation_matrix.m() == x_source_fe.dofs_per_face,
867  ExcDimensionMismatch (interpolation_matrix.m(),
868  x_source_fe.dofs_per_face));
869 
870  // ok, source is a Q_Hierarchical element, so
871  // we will be able to do the work
872  const FE_Q_Hierarchical<dim> &source_fe
873  = dynamic_cast<const FE_Q_Hierarchical<dim>&>(x_source_fe);
874  (void)source_fe;
875 
876  // Make sure, that the element,
877  // for which the DoFs should be
878  // constrained is the one with
879  // the higher polynomial degree.
880  // Actually the procedure will work
881  // also if this assertion is not
882  // satisfied. But the matrices
883  // produced in that case might
884  // lead to problems in the
885  // hp procedures, which use this
886  // method.
887  Assert (this->dofs_per_face <= source_fe.dofs_per_face,
888  typename FEL::
889  ExcInterpolationNotImplemented ());
890  interpolation_matrix = 0;
891 
892  switch (dim)
893  {
894  case 2:
895  {
896  // In 2-dimension the constraints are trivial.
897  // First this->dofs_per_face DoFs of the constrained
898  // element are made equal to the current (dominating)
899  // element, which corresponds to 1 on diagonal of the matrix.
900  // DoFs which correspond to higher polynomials
901  // are zeroed (zero rows in the matrix).
902  for (unsigned int i = 0; i < this->dofs_per_face; ++i)
903  interpolation_matrix (i, i) = 1;
904 
905  break;
906  }
907 
908  case 3:
909  {
910  for (unsigned int i = 0; i < GeometryInfo<3>::vertices_per_face; ++i)
911  interpolation_matrix (i, i) = 1;
912 
913  for (unsigned int i = 0; i < this->degree - 1; ++i)
914  {
915  for (unsigned int j = 0; j < GeometryInfo<3>::lines_per_face; ++j)
916  interpolation_matrix (
917  i + j * (x_source_fe.degree - 1) + GeometryInfo<3>::vertices_per_face,
918  i + j * (this->degree - 1) + GeometryInfo<3>::vertices_per_face) = 1;
919 
920  for (unsigned int j = 0; j < this->degree - 1; ++j)
921  interpolation_matrix (
922  (i + GeometryInfo<3>::lines_per_face) * (x_source_fe.degree - 1) + j
923  + GeometryInfo<3>::vertices_per_face,
924  (i + GeometryInfo<3>::lines_per_face) * (this->degree - 1) + j
925  + GeometryInfo<3>::vertices_per_face) = 1;
926  }
927  }
928  }
929 }
930 
931 
932 
933 template <int dim>
934 void
935 FE_Q_Hierarchical<dim>::
936 get_subface_interpolation_matrix (const FiniteElement<dim> &x_source_fe,
937  const unsigned int subface,
938  FullMatrix<double> &interpolation_matrix) const
939 {
940  // this is only implemented, if the
941  // source FE is also a
942  // Q_Hierarchical element
943  typedef FE_Q_Hierarchical<dim> FEQHierarchical;
944  typedef FiniteElement<dim> FEL;
945  AssertThrow ((x_source_fe.get_name().find ("FE_Q_Hierarchical<") == 0)
946  ||
947  (dynamic_cast<const FEQHierarchical *>(&x_source_fe) != 0),
948  typename FEL::
949  ExcInterpolationNotImplemented());
950 
951  Assert (interpolation_matrix.n() == this->dofs_per_face,
952  ExcDimensionMismatch (interpolation_matrix.n(),
953  this->dofs_per_face));
954  Assert (interpolation_matrix.m() == x_source_fe.dofs_per_face,
955  ExcDimensionMismatch (interpolation_matrix.m(),
956  x_source_fe.dofs_per_face));
957 
958  // ok, source is a Q_Hierarchical element, so
959  // we will be able to do the work
960  const FE_Q_Hierarchical<dim> &source_fe
961  = dynamic_cast<const FE_Q_Hierarchical<dim>&>(x_source_fe);
962 
963  // Make sure, that the element,
964  // for which the DoFs should be
965  // constrained is the one with
966  // the higher polynomial degree.
967  // Actually the procedure will work
968  // also if this assertion is not
969  // satisfied. But the matrices
970  // produced in that case might
971  // lead to problems in the
972  // hp procedures, which use this
973  // method.
974  Assert (this->dofs_per_face <= source_fe.dofs_per_face,
975  typename FEL::
976  ExcInterpolationNotImplemented ());
977 
978  switch (dim)
979  {
980  case 2:
981  {
982  switch (subface)
983  {
984  case 0:
985  {
986  interpolation_matrix (0, 0) = 1.0;
987  interpolation_matrix (1, 0) = 0.5;
988  interpolation_matrix (1, 1) = 0.5;
989 
990  for (unsigned int dof = 2; dof < this->dofs_per_face;)
991  {
992  interpolation_matrix (1, dof) = -1.0;
993  dof = dof + 2;
994  }
995 
996  int factorial_i = 1;
997  int factorial_ij;
998  int factorial_j;
999 
1000  for (int i = 2; i < (int) this->dofs_per_face; ++i)
1001  {
1002  interpolation_matrix (i, i) = std::pow (0.5, i);
1003  factorial_i *= i;
1004  factorial_j = factorial_i;
1005  factorial_ij = 1;
1006 
1007  for (int j = i + 1; j < (int) this->dofs_per_face; ++j)
1008  {
1009  factorial_ij *= j - i;
1010  factorial_j *= j;
1011 
1012  if ((i + j) & 1)
1013  interpolation_matrix (i, j)
1014  = -1.0 * std::pow (0.5, j) *
1015  factorial_j / (factorial_i * factorial_ij);
1016 
1017  else
1018  interpolation_matrix (i, j)
1019  = std::pow (0.5, j) *
1020  factorial_j / (factorial_i * factorial_ij);
1021  }
1022  }
1023 
1024  break;
1025  }
1026 
1027  case 1:
1028  {
1029  interpolation_matrix (0, 0) = 0.5;
1030  interpolation_matrix (0, 1) = 0.5;
1031 
1032  for (unsigned int dof = 2; dof < this->dofs_per_face;)
1033  {
1034  interpolation_matrix (0, dof) = -1.0;
1035  dof = dof + 2;
1036  }
1037 
1038  interpolation_matrix (1, 1) = 1.0;
1039 
1040  int factorial_i = 1;
1041  int factorial_ij;
1042  int factorial_j;
1043 
1044  for (int i = 2; i < (int) this->dofs_per_face; ++i)
1045  {
1046  interpolation_matrix (i, i) = std::pow (0.5, i);
1047  factorial_i *= i;
1048  factorial_j = factorial_i;
1049  factorial_ij = 1;
1050 
1051  for (int j = i + 1; j < (int) this->dofs_per_face; ++j)
1052  {
1053  factorial_ij *= j - i;
1054  factorial_j *= j;
1055  interpolation_matrix (i, j)
1056  = std::pow (0.5, j) * factorial_j / (factorial_i * factorial_ij);
1057  }
1058  }
1059  }
1060  }
1061 
1062  break;
1063  }
1064 
1065  case 3:
1066  {
1067  switch (subface)
1068  {
1069  case 0:
1070  {
1071  interpolation_matrix (0, 0) = 1.0;
1072  interpolation_matrix (1, 0) = 0.5;
1073  interpolation_matrix (1, 1) = 0.5;
1074  interpolation_matrix (2, 0) = 0.5;
1075  interpolation_matrix (2, 2) = 0.5;
1076 
1077  for (unsigned int i = 0; i < this->degree - 1;)
1078  {
1079  for (unsigned int line = 0; line < GeometryInfo<3>::lines_per_face; ++line)
1080  interpolation_matrix (3, i + line * (this->degree - 1) + 4) = -0.5;
1081 
1082  for (unsigned int j = 0; j < this->degree - 1;)
1083  {
1084  interpolation_matrix (3, i + (j + 4) * this->degree - j) = 1.0;
1085  j = j + 2;
1086  }
1087 
1088  interpolation_matrix (1, i + 2 * (this->degree + 1)) = -1.0;
1089  interpolation_matrix (2, i + 4) = -1.0;
1090  i = i + 2;
1091  }
1092 
1093  for (unsigned int vertex = 0; vertex < GeometryInfo<3>::vertices_per_face; ++vertex)
1094  interpolation_matrix (3, vertex) = 0.25;
1095 
1096  int factorial_i = 1;
1097  int factorial_ij;
1098  int factorial_j;
1099  int factorial_k;
1100  int factorial_kl;
1101  int factorial_l;
1102 
1103  for (int i = 2; i <= (int) this->degree; ++i)
1104  {
1105  double tmp = std::pow (0.5, i);
1106  interpolation_matrix (i + 2, i + 2) = tmp;
1107  interpolation_matrix (i + 2 * source_fe.degree, i + 2 * this->degree) = tmp;
1108  tmp *= 0.5;
1109  interpolation_matrix (i + source_fe.degree + 1, i + 2) = tmp;
1110  interpolation_matrix (i + source_fe.degree + 1, i + this->degree + 1) = tmp;
1111  interpolation_matrix (i + 3 * source_fe.degree - 1, i + 2 * this->degree) = tmp;
1112  interpolation_matrix (i + 3 * source_fe.degree - 1, i + 3 * this->degree - 1) = tmp;
1113  tmp *= -2.0;
1114 
1115  for (unsigned int j = 0; j < this->degree - 1;)
1116  {
1117  interpolation_matrix (i + source_fe.degree + 1, (i + 2) * this->degree + j + 2 - i) = tmp;
1118  interpolation_matrix (i + 3 * source_fe.degree - 1, i + (j + 4) * this->degree - j - 2) = tmp;
1119  j = j + 2;
1120  }
1121 
1122  factorial_k = 1;
1123 
1124  for (int j = 2; j <= (int) this->degree; ++j)
1125  {
1126  interpolation_matrix (i + (j + 2) * source_fe.degree - j, i + (j + 2) * this->degree - j) = std::pow (0.5, i + j);
1127  factorial_k *= j;
1128  factorial_kl = 1;
1129  factorial_l = factorial_k;
1130 
1131  for (int k = j + 1; k < (int) this->degree; ++k)
1132  {
1133  factorial_kl *= k - j;
1134  factorial_l *= k;
1135 
1136  if ((j + k) & 1)
1137  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);
1138 
1139  else
1140  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);
1141  }
1142  }
1143 
1144  factorial_i *= i;
1145  factorial_j = factorial_i;
1146  factorial_ij = 1;
1147 
1148  for (int j = i + 1; j <= (int) this->degree; ++j)
1149  {
1150  factorial_ij *= j - i;
1151  factorial_j *= j;
1152 
1153  if ((i + j) & 1)
1154  {
1155  tmp = -1.0 * std::pow (0.5, j) * factorial_j / (factorial_i * factorial_ij);
1156  interpolation_matrix (i + 2, j + 2) = tmp;
1157  interpolation_matrix (i + 2 * source_fe.degree, j + 2 * this->degree) = tmp;
1158  factorial_k = 1;
1159 
1160  for (int k = 2; k <= (int) this->degree; ++k)
1161  {
1162  interpolation_matrix (i + (k + 2) * source_fe.degree - k, j + (k + 2) * this->degree - k) = tmp * std::pow (0.5, k);
1163  factorial_k *= k;
1164  factorial_l = factorial_k;
1165  factorial_kl = 1;
1166 
1167  for (int l = k + 1; l <= (int) this->degree; ++l)
1168  {
1169  factorial_kl *= l - k;
1170  factorial_l *= l;
1171 
1172  if ((k + l) & 1)
1173  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);
1174 
1175  else
1176  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);
1177  }
1178  }
1179 
1180  tmp *= 0.5;
1181  interpolation_matrix (i + source_fe.degree + 1, j + 2) = tmp;
1182  interpolation_matrix (i + source_fe.degree + 1, j + this->degree + 1) = tmp;
1183  interpolation_matrix (i + 3 * source_fe.degree - 1, j + 2 * this->degree) = tmp;
1184  interpolation_matrix (i + 3 * source_fe.degree - 1, j + 3 * this->degree - 1) = tmp;
1185  tmp *= -2.0;
1186 
1187  for (unsigned int k = 0; k < this->degree - 1;)
1188  {
1189  interpolation_matrix (i + source_fe.degree + 1, (j + 2) * this->degree + k + 2 - j) = tmp;
1190  interpolation_matrix (i + 3 * source_fe.degree - 1, j + (k + 4) * this->degree - k - 2) = tmp;
1191  k = k + 2;
1192  }
1193  }
1194  else
1195  {
1196  tmp = std::pow (0.5, j) * factorial_j / (factorial_i * factorial_ij);
1197  interpolation_matrix (i + 2, j + 2) = tmp;
1198  interpolation_matrix (i + 2 * source_fe.degree, j + 2 * this->degree) = tmp;
1199  factorial_k = 1;
1200 
1201  for (int k = 2; k <= (int) this->degree; ++k)
1202  {
1203  interpolation_matrix (i + (k + 2) * source_fe.degree - k, j + (k + 2) * this->degree - k) = tmp * std::pow (0.5, k);
1204  factorial_k *= k;
1205  factorial_l = factorial_k;
1206  factorial_kl = 1;
1207 
1208  for (int l = k + 1; l <= (int) this->degree; ++l)
1209  {
1210  factorial_kl *= l - k;
1211  factorial_l *= l;
1212 
1213  if ((k + l) & 1)
1214  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);
1215 
1216  else
1217  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);
1218  }
1219  }
1220 
1221  tmp *= 0.5;
1222  interpolation_matrix (i + source_fe.degree + 1, j + 2) = tmp;
1223  interpolation_matrix (i + source_fe.degree + 1, j + this->degree + 1) = tmp;
1224  interpolation_matrix (i + 3 * source_fe.degree - 1, j + 2 * this->degree) = tmp;
1225  interpolation_matrix (i + 3 * source_fe.degree - 1, j + 3 * this->degree - 1) = tmp;
1226  tmp *= -2.0;
1227 
1228  for (unsigned int k = 0; k < this->degree - 1;)
1229  {
1230  interpolation_matrix (i + source_fe.degree + 1, (j + 2) * this->degree + k + 2 - j) = tmp;
1231  interpolation_matrix (i + 3 * source_fe.degree - 1, j + (k + 4) * this->degree - k - 2) = tmp;
1232  k = k + 2;
1233  }
1234  }
1235  }
1236  }
1237 
1238  break;
1239  }
1240 
1241  case 1:
1242  {
1243  interpolation_matrix (0, 0) = 0.5;
1244  interpolation_matrix (0, 1) = 0.5;
1245  interpolation_matrix (1, 1) = 1.0;
1246  interpolation_matrix (3, 1) = 0.5;
1247  interpolation_matrix (3, 3) = 0.5;
1248 
1249  for (unsigned int i = 0; i < this->degree - 1;)
1250  {
1251  for (unsigned int line = 0; line < GeometryInfo<3>::lines_per_face; ++line)
1252  interpolation_matrix (2, i + line * (this->degree - 1) + 4) = -0.5;
1253 
1254  for (unsigned int j = 0; j < this->degree - 1;)
1255  {
1256  interpolation_matrix (2, i + (j + 4) * this->degree - j) = 1.0;
1257  j = j + 2;
1258  }
1259 
1260  interpolation_matrix (0, i + 2 * (this->degree + 1)) = -1.0;
1261  interpolation_matrix (3, i + this->degree + 3) = -1.0;
1262  i = i + 2;
1263  }
1264 
1265  for (unsigned int vertex = 0; vertex < GeometryInfo<3>::vertices_per_face; ++vertex)
1266  interpolation_matrix (2, vertex) = 0.25;
1267 
1268  int factorial_i = 1;
1269  int factorial_ij;
1270  int factorial_j;
1271  int factorial_k;
1272  int factorial_kl;
1273  int factorial_l;
1274 
1275  for (int i = 2; i <= (int) this->degree; ++i)
1276  {
1277  double tmp = std::pow (0.5, i + 1);
1278  interpolation_matrix (i + 2, i + 2) = tmp;
1279  interpolation_matrix (i + 2, i + this->degree + 1) = tmp;
1280  interpolation_matrix (i + 3 * source_fe.degree - 1, i + 2 * this->degree) = tmp;
1281  interpolation_matrix (i + 3 * source_fe.degree - 1, i + 3 * this->degree - 1) = tmp;
1282  tmp *= -2.0;
1283 
1284  for (unsigned int j = 0; j < this->degree - 1;)
1285  {
1286  interpolation_matrix (i + 2, j + (i + 2) * this->degree + 2 - i) = tmp;
1287  interpolation_matrix (i + 3 * source_fe.degree - 1, i + (j + 4) * this->degree - j - 2) = tmp;
1288  j = j + 2;
1289  }
1290 
1291  tmp *= - 1.0;
1292  interpolation_matrix (i + source_fe.degree + 1, i + this->degree + 1) = tmp;
1293  interpolation_matrix (i + 2 * source_fe.degree, i + 2 * this->degree) = tmp;
1294  factorial_i *= i;
1295  factorial_j = factorial_i;
1296  factorial_ij = 1;
1297 
1298  for (int j = i + 1; j <= (int) this->degree; ++j)
1299  {
1300  factorial_ij *= j - i;
1301  factorial_j *= j;
1302  tmp = std::pow (0.5, j) * factorial_j / (factorial_i * factorial_ij);
1303  interpolation_matrix (i + 2 * source_fe.degree, j + 2 * this->degree) = tmp;
1304  factorial_k = 1;
1305 
1306  for (int k = 2; k <= (int) this->degree; ++k)
1307  {
1308  interpolation_matrix (i + (k + 2) * source_fe.degree - k, j + (k + 2) * this->degree - k) = tmp * std::pow (0.5, k);
1309  factorial_k *= k;
1310  factorial_l = factorial_k;
1311  factorial_kl = 1;
1312 
1313  for (int l = k + 1; l <= (int) this->degree; ++l)
1314  {
1315  factorial_kl *= l - k;
1316  factorial_l *= l;
1317 
1318  if ((k + l) & 1)
1319  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);
1320 
1321  else
1322  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);
1323  }
1324  }
1325 
1326  tmp *= -1.0;
1327 
1328  for (unsigned int k = 0; k < this->degree - 1;)
1329  {
1330  interpolation_matrix (i + 3 * source_fe.degree - 1, j + (k + 4) * this->degree - k - 2) = tmp;
1331  k = k + 2;
1332  }
1333 
1334  tmp *= -0.5;
1335  interpolation_matrix (i + 3 * source_fe.degree - 1, j + 2 * this->degree) = tmp;
1336  interpolation_matrix (i + 3 * source_fe.degree - 1, j + 3 * this->degree - 1) = tmp;
1337 
1338  if ((i + j) & 1)
1339  tmp *= -1.0;
1340 
1341  interpolation_matrix (i + 2, j + 2) = tmp;
1342  interpolation_matrix (i + 2, j + this->degree + 1) = tmp;
1343  interpolation_matrix (i + source_fe.degree + 1, j + this->degree + 1) = 2.0 * tmp;
1344  tmp *= -2.0;
1345 
1346  for (unsigned int k = 0; k < this->degree - 1;)
1347  {
1348  interpolation_matrix (i + 2, k + (j + 2) * this->degree + 2 - j) = tmp;
1349  k = k + 2;
1350  }
1351  }
1352 
1353  factorial_k = 1;
1354 
1355  for (int j = 2; j <= (int) this->degree; ++j)
1356  {
1357  interpolation_matrix (i + (j + 2) * source_fe.degree - j, i + (j + 2) * this->degree - j) = std::pow (0.5, i + j);
1358  factorial_k *= j;
1359  factorial_l = factorial_k;
1360  factorial_kl = 1;
1361 
1362  for (int k = j + 1; k <= (int) this->degree; ++k)
1363  {
1364  factorial_kl *= k - j;
1365  factorial_l *= k;
1366 
1367  if ((j + k) & 1)
1368  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);
1369 
1370  else
1371  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);
1372  }
1373  }
1374  }
1375 
1376  break;
1377  }
1378 
1379  case 2:
1380  {
1381  interpolation_matrix (0, 0) = 0.5;
1382  interpolation_matrix (0, 2) = 0.5;
1383  interpolation_matrix (2, 2) = 1.0;
1384  interpolation_matrix (3, 2) = 0.5;
1385  interpolation_matrix (3, 3) = 0.5;
1386 
1387  for (unsigned int i = 0; i < this->degree - 1;)
1388  {
1389  for (unsigned int line = 0; line < GeometryInfo<3>::lines_per_face; ++line)
1390  interpolation_matrix (1, i + line * (this->degree - 1) + 4) = -0.5;
1391 
1392  for (unsigned int j = 0; j < this->degree - 1;)
1393  {
1394  interpolation_matrix (1, i + (j + 4) * this->degree - j) = 1.0;
1395  j = j + 2;
1396  }
1397 
1398  interpolation_matrix (0, i + 4) = -1.0;
1399  interpolation_matrix (3, i + 3 * this->degree + 1) = -1.0;
1400  i = i + 2;
1401  }
1402 
1403  for (unsigned int vertex = 0; vertex < GeometryInfo<3>::vertices_per_face; ++vertex)
1404  interpolation_matrix (1, vertex) = 0.25;
1405 
1406  int factorial_i = 1;
1407  int factorial_ij;
1408  int factorial_j;
1409  int factorial_k;
1410  int factorial_kl;
1411  int factorial_l;
1412 
1413  for (int i = 2; i <= (int) this->degree; ++i)
1414  {
1415  double tmp = std::pow (0.5, i);
1416  interpolation_matrix (i + 2, i + 2) = tmp;
1417  interpolation_matrix (i + 3 * source_fe.degree - 1, i + 3 * this->degree - 1) = tmp;
1418  tmp *= 0.5;
1419  interpolation_matrix (i + source_fe.degree + 1, i + 2) = tmp;
1420  interpolation_matrix (i + source_fe.degree + 1, i + this->degree + 1) = tmp;
1421  interpolation_matrix (i + 2 * source_fe.degree, i + 2 * this->degree) = tmp;
1422  interpolation_matrix (i + 2 * source_fe.degree, i + 3 * this->degree - 1) = tmp;
1423  tmp *= -2.0;
1424 
1425  for (unsigned int j = 0; j < this->degree - 1;)
1426  {
1427  interpolation_matrix (i + source_fe.degree + 1, j + (i + 2) * this->degree + 2 - i) = tmp;
1428  interpolation_matrix (i + 2 * source_fe.degree, i + (j + 4) * this->degree - j - 2) = tmp;
1429  j = j + 2;
1430  }
1431 
1432  factorial_k = 1;
1433 
1434  for (int j = 2; j <= (int) this->degree; ++j)
1435  {
1436  interpolation_matrix (i + (j + 2) * source_fe.degree - j, i + (j + 2) * this->degree - j) = std::pow (0.5, i + j);
1437  factorial_k *= j;
1438  factorial_l = factorial_k;
1439  factorial_kl = 1;
1440 
1441  for (int k = j + 1; k <= (int) this->degree; ++k)
1442  {
1443  factorial_kl *= k - j;
1444  factorial_l *= k;
1445  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);
1446  }
1447  }
1448 
1449  factorial_i *= i;
1450  factorial_j = factorial_i;
1451  factorial_ij = 1;
1452 
1453  for (int j = i + 1; j <= (int) this->degree; ++j)
1454  {
1455  factorial_ij *= j - i;
1456  factorial_j *= j;
1457  tmp = std::pow (0.5, j) * factorial_j / (factorial_i * factorial_ij);
1458  interpolation_matrix (i + 2, j + 2) = tmp;
1459  tmp *= -1.0;
1460 
1461  for (unsigned int k = 0; k < this->degree - 1;)
1462  {
1463  interpolation_matrix (i + source_fe.degree + 1, k + (j + 2) * this->degree + 2 - j) = tmp;
1464  k = k + 2;
1465  }
1466 
1467  tmp *= -0.5;
1468  interpolation_matrix (i + source_fe.degree + 1, j + 2) = tmp;
1469  interpolation_matrix (i + source_fe.degree + 1, j + this->degree + 1) = tmp;
1470 
1471  if ((i + j) & 1)
1472  tmp *= -1.0;
1473 
1474  interpolation_matrix (i + 2 * source_fe.degree, j + 2 * this->degree) = tmp;
1475  interpolation_matrix (i + 2 * source_fe.degree, j + 3 * this->degree - 1) = tmp;
1476  tmp *= 2.0;
1477  interpolation_matrix (i + 3 * source_fe.degree - 1, j + 3 * this->degree - 1) = tmp;
1478  factorial_k = 1;
1479 
1480  for (int k = 2; k <= (int) this->degree; ++k)
1481  {
1482  interpolation_matrix (i + (k + 2) * source_fe.degree - k, j + (k + 2) * this->degree - k) = tmp * std::pow (0.5, k);
1483  factorial_k *= k;
1484  factorial_l = factorial_k;
1485  factorial_kl = 1;
1486 
1487  for (int l = k + 1; l <= (int) this->degree; ++l)
1488  {
1489  factorial_kl *= l - k;
1490  factorial_l *= l;
1491  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);
1492  }
1493  }
1494 
1495  tmp *= -1.0;
1496 
1497  for (unsigned int k = 0; k < this->degree - 1;)
1498  {
1499  interpolation_matrix (i + 2 * source_fe.degree, j + (k + 4) * this->degree - k - 2) = tmp;
1500  k = k + 2;
1501  }
1502  }
1503  }
1504 
1505  break;
1506  }
1507 
1508  case 3:
1509  {
1510  for (unsigned int vertex = 0; vertex < GeometryInfo<3>::vertices_per_face; ++vertex)
1511  interpolation_matrix (0, vertex) = 0.25;
1512 
1513  for (unsigned int i = 0; i < this->degree - 1;)
1514  {
1515  for (unsigned int line = 0; line < GeometryInfo<3>::lines_per_face; ++line)
1516  interpolation_matrix (0, i + line * (this->degree - 1) + 4) = -0.5;
1517 
1518  for (unsigned int j = 0; j < this->degree - 1;)
1519  {
1520  interpolation_matrix (0, i + (j + 4) * this->degree - j) = 1.0;
1521  j = j + 2;
1522  }
1523 
1524  interpolation_matrix (1, i + 4) = -1.0;
1525  interpolation_matrix (2, i + 3 * this->degree + 1) = -1.0;
1526  i = i + 2;
1527  }
1528 
1529  interpolation_matrix (1, 0) = 0.5;
1530  interpolation_matrix (1, 1) = 0.5;
1531  interpolation_matrix (2, 2) = 0.5;
1532  interpolation_matrix (2, 3) = 0.5;
1533  interpolation_matrix (3, 3) = 1.0;
1534 
1535  int factorial_i = 1;
1536  int factorial_ij;
1537  int factorial_j;
1538  int factorial_k;
1539  int factorial_kl;
1540  int factorial_l;
1541 
1542  for (int i = 2; i <= (int) this->degree; ++i)
1543  {
1544  double tmp = std::pow (0.5, i + 1);
1545  interpolation_matrix (i + 2, i + 2) = tmp;
1546  interpolation_matrix (i + 2, i + this->degree + 1) = tmp;
1547  interpolation_matrix (i + 2 * source_fe.degree, i + 2 * this->degree) = tmp;
1548  interpolation_matrix (i + 2 * source_fe.degree, i + 3 * this->degree - 1) = tmp;
1549  tmp *= -2.0;
1550 
1551  for (unsigned int j = 0; j < this->degree - 1;)
1552  {
1553  interpolation_matrix (i + 2, j + (i + 2) * this->degree + 2 - i) = tmp;
1554  interpolation_matrix (i + 2 * source_fe.degree, i + (j + 4) * this->degree - 2) = tmp;
1555  j = j + 2;
1556  }
1557 
1558  tmp *= -1.0;
1559  interpolation_matrix (i + source_fe.degree + 1, i + this->degree + 1) = tmp;
1560  interpolation_matrix (i + 3 * source_fe.degree - 1, i + 3 * this->degree - 1) = tmp;
1561  factorial_k = 1;
1562 
1563  for (int j = 2; j <= (int) this->degree; ++j)
1564  {
1565  interpolation_matrix (i + (j + 2) * source_fe.degree - j, i + (j + 2) * this->degree - j) = std::pow (0.5, i + j);
1566  factorial_k *= j;
1567  factorial_l = factorial_k;
1568  factorial_kl = 1;
1569 
1570  for (int k = j + 1; k <= (int) this->degree; ++k)
1571  {
1572  factorial_kl *= k - j;
1573  factorial_l *= k;
1574  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);
1575  }
1576  }
1577 
1578  factorial_i *= i;
1579  factorial_j = factorial_i;
1580  factorial_ij = 1;
1581 
1582  for (int j = i + 1; j <= (int) this->degree; ++j)
1583  {
1584  factorial_ij *= j - i;
1585  factorial_j *= j;
1586  tmp = std::pow (0.5, j + 1) * factorial_j / (factorial_i * factorial_ij);
1587  interpolation_matrix (i + 2, j + 2) = tmp;
1588  interpolation_matrix (i + 2, j + this->degree + 1) = tmp;
1589  interpolation_matrix (i + 2 * source_fe.degree, j + 2 * this->degree) = tmp;
1590  interpolation_matrix (i + 2 * source_fe.degree, j + 3 * this->degree - 1) = tmp;
1591  tmp *= 2.0;
1592  interpolation_matrix (i + source_fe.degree + 1, j + this->degree + 1) = tmp;
1593  interpolation_matrix (i + 3 * source_fe.degree - 1, j + 3 * this->degree - 1) = tmp;
1594  factorial_k = 1;
1595 
1596  for (int k = 2; k <= (int) this->degree; ++k)
1597  {
1598  interpolation_matrix (i + (k + 2) * source_fe.degree - k, j + (k + 2) * this->degree - k) = tmp * std::pow (0.5, k);
1599  factorial_k *= k;
1600  factorial_l = factorial_k;
1601  factorial_kl = 1;
1602 
1603  for (int l = k + 1; l <= (int) this->degree; ++l)
1604  {
1605  factorial_kl *= l - k;
1606  factorial_l *= l;
1607  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);
1608  }
1609  }
1610 
1611  tmp *= -1.0;
1612 
1613  for (unsigned int k = 0; k < this->degree - 1;)
1614  {
1615  interpolation_matrix (i + 2, k + (j + 2) * this->degree + 2 - j) = tmp;
1616  interpolation_matrix (i + 2 * source_fe.degree, j + (k + 4) * this->degree - 2) = tmp;
1617  k = k + 2;
1618  }
1619  }
1620  }
1621  }
1622  }
1623  }
1624  }
1625 }
1626 
1627 
1628 
1629 template <int dim>
1630 void FE_Q_Hierarchical<dim>::initialize_unit_face_support_points ()
1631 {
1632  const unsigned int codim = dim-1;
1633 
1634  // number of points: (degree+1)^codim
1635  unsigned int n = this->degree+1;
1636  for (unsigned int i=1; i<codim; ++i)
1637  n *= this->degree+1;
1638 
1639  this->unit_face_support_points.resize(n);
1640 
1641  Point<codim> p;
1642 
1643  unsigned int k=0;
1644  for (unsigned int iz=0; iz <= ((codim>2) ? this->degree : 0) ; ++iz)
1645  for (unsigned int iy=0; iy <= ((codim>1) ? this->degree : 0) ; ++iy)
1646  for (unsigned int ix=0; ix<=this->degree; ++ix)
1647  {
1648  if (ix==0)
1649  p(0) = 0.;
1650  else if (ix==1)
1651  p(0) = 1.;
1652  else
1653  p(0) = .5;
1654  if (codim>1)
1655  {
1656  if (iy==0)
1657  p(1) = 0.;
1658  else if (iy==1)
1659  p(1) = 1.;
1660  else
1661  p(1) = .5;
1662  }
1663  if (codim>2)
1664  {
1665  if (iz==0)
1666  p(2) = 0.;
1667  else if (iz==1)
1668  p(2) = 1.;
1669  else
1670  p(2) = .5;
1671  }
1672  this->unit_face_support_points[face_renumber[k++]] = p;
1673  };
1674 }
1675 
1676 
1677 // we use same dpo_vector as FE_Q
1678 template <int dim>
1679 std::vector<unsigned int>
1680 FE_Q_Hierarchical<dim>::get_dpo_vector(const unsigned int deg)
1681 {
1682  std::vector<unsigned int> dpo(dim+1, 1U);
1683  for (unsigned int i=1; i<dpo.size(); ++i)
1684  dpo[i]=dpo[i-1]*(deg-1);
1685  return dpo;
1686 }
1687 
1688 
1689 
1690 template <int dim>
1691 std::vector<unsigned int>
1692 FE_Q_Hierarchical<dim>::
1693 hierarchic_to_fe_q_hierarchical_numbering (const FiniteElementData<dim> &fe)
1694 {
1695  Assert (fe.n_components() == 1, ExcInternalError());
1696  std::vector<unsigned int> h2l(fe.dofs_per_cell);
1697 
1698  // polynomial degree
1699  const unsigned int degree = fe.dofs_per_line+1;
1700  // number of grid points in each
1701  // direction
1702  const unsigned int n = degree+1;
1703 
1704  // the following lines of code are
1705  // somewhat odd, due to the way the
1706  // hierarchic numbering is
1707  // organized. if someone would
1708  // really want to understand these
1709  // lines, you better draw some
1710  // pictures where you indicate the
1711  // indices and orders of vertices,
1712  // lines, etc, along with the
1713  // numbers of the degrees of
1714  // freedom in hierarchical and
1715  // lexicographical order
1716  switch (dim)
1717  {
1718  case 1:
1719  {
1720  for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
1721  h2l[i] = i;
1722 
1723  break;
1724  }
1725 
1726  case 2:
1727  {
1728  // Example: degree=3
1729  //
1730  // hierarchical numbering:
1731  // 2 10 11 3
1732  // 5 14 15 7
1733  // 4 12 13 6
1734  // 0 8 9 1
1735  //
1736  // fe_q_hierarchical numbering:
1737  // 4 6 7 5
1738  // 12 14 15 13
1739  // 8 10 11 9
1740  // 0 2 3 1
1741  unsigned int next_index = 0;
1742  // first the four vertices
1743  h2l[next_index++] = 0;
1744  h2l[next_index++] = 1;
1745  h2l[next_index++] = n;
1746  h2l[next_index++] = n+1;
1747  // left line
1748  for (unsigned int i=0; i<fe.dofs_per_line; ++i)
1749  h2l[next_index++] = (2+i)*n;
1750  // right line
1751  for (unsigned int i=0; i<fe.dofs_per_line; ++i)
1752  h2l[next_index++] = (2+i)*n+1;
1753  // bottom line
1754  for (unsigned int i=0; i<fe.dofs_per_line; ++i)
1755  h2l[next_index++] = 2+i;
1756  // top line
1757  for (unsigned int i=0; i<fe.dofs_per_line; ++i)
1758  h2l[next_index++] = n+2+i;
1759  // inside quad
1760  Assert (fe.dofs_per_quad == fe.dofs_per_line*fe.dofs_per_line,
1761  ExcInternalError());
1762  for (unsigned int i=0; i<fe.dofs_per_line; ++i)
1763  for (unsigned int j=0; j<fe.dofs_per_line; ++j)
1764  h2l[next_index++] = (2+i)*n+2+j;
1765 
1766  Assert (next_index == fe.dofs_per_cell, ExcInternalError());
1767 
1768  break;
1769  }
1770 
1771  case 3:
1772  {
1773  unsigned int next_index = 0;
1774  const unsigned int n2=n*n;
1775  // first the eight vertices
1776  // bottom face, lexicographic
1777  h2l[next_index++] = 0;
1778  h2l[next_index++] = 1;
1779  h2l[next_index++] = n;
1780  h2l[next_index++] = n+1;
1781  // top face, lexicographic
1782  h2l[next_index++] = n2;
1783  h2l[next_index++] = n2+1;
1784  h2l[next_index++] = n2+n;
1785  h2l[next_index++] = n2+n+1;
1786 
1787  // now the lines
1788  // bottom face
1789  for (unsigned int i=0; i<fe.dofs_per_line; ++i)
1790  h2l[next_index++] = (2+i)*n;
1791  for (unsigned int i=0; i<fe.dofs_per_line; ++i)
1792  h2l[next_index++] = (2+i)*n+1;
1793  for (unsigned int i=0; i<fe.dofs_per_line; ++i)
1794  h2l[next_index++] = 2+i;
1795  for (unsigned int i=0; i<fe.dofs_per_line; ++i)
1796  h2l[next_index++] = n+2+i;
1797  // top face
1798  for (unsigned int i=0; i<fe.dofs_per_line; ++i)
1799  h2l[next_index++] = n2+(2+i)*n;
1800  for (unsigned int i=0; i<fe.dofs_per_line; ++i)
1801  h2l[next_index++] = n2+(2+i)*n+1;
1802  for (unsigned int i=0; i<fe.dofs_per_line; ++i)
1803  h2l[next_index++] = n2+2+i;
1804  for (unsigned int i=0; i<fe.dofs_per_line; ++i)
1805  h2l[next_index++] = n2+n+2+i;
1806  // lines in z-direction
1807  for (unsigned int i=0; i<fe.dofs_per_line; ++i)
1808  h2l[next_index++] = (2+i)*n2;
1809  for (unsigned int i=0; i<fe.dofs_per_line; ++i)
1810  h2l[next_index++] = (2+i)*n2+1;
1811  for (unsigned int i=0; i<fe.dofs_per_line; ++i)
1812  h2l[next_index++] = (2+i)*n2+n;
1813  for (unsigned int i=0; i<fe.dofs_per_line; ++i)
1814  h2l[next_index++] = (2+i)*n2+n+1;
1815 
1816  // inside quads
1817  Assert (fe.dofs_per_quad == fe.dofs_per_line*fe.dofs_per_line,
1818  ExcInternalError());
1819  // left face
1820  for (unsigned int i=0; i<fe.dofs_per_line; ++i)
1821  for (unsigned int j=0; j<fe.dofs_per_line; ++j)
1822  h2l[next_index++] = (2+i)*n2+(2+j)*n;
1823  // right face
1824  for (unsigned int i=0; i<fe.dofs_per_line; ++i)
1825  for (unsigned int j=0; j<fe.dofs_per_line; ++j)
1826  h2l[next_index++] = (2+i)*n2+(2+j)*n+1;
1827  // front face
1828  for (unsigned int i=0; i<fe.dofs_per_line; ++i)
1829  for (unsigned int j=0; j<fe.dofs_per_line; ++j)
1830  h2l[next_index++] = (2+i)*n2+2+j;
1831  // back face
1832  for (unsigned int i=0; i<fe.dofs_per_line; ++i)
1833  for (unsigned int j=0; j<fe.dofs_per_line; ++j)
1834  h2l[next_index++] = (2+i)*n2+n+2+j;
1835  // bottom face
1836  for (unsigned int i=0; i<fe.dofs_per_line; ++i)
1837  for (unsigned int j=0; j<fe.dofs_per_line; ++j)
1838  h2l[next_index++] = (2+i)*n+2+j;
1839  // top face
1840  for (unsigned int i=0; i<fe.dofs_per_line; ++i)
1841  for (unsigned int j=0; j<fe.dofs_per_line; ++j)
1842  h2l[next_index++] = n2+(2+i)*n+2+j;
1843 
1844  // inside hex
1845  Assert (fe.dofs_per_hex == fe.dofs_per_quad*fe.dofs_per_line,
1846  ExcInternalError());
1847  for (unsigned int i=0; i<fe.dofs_per_line; ++i)
1848  for (unsigned int j=0; j<fe.dofs_per_line; ++j)
1849  for (unsigned int k=0; k<fe.dofs_per_line; ++k)
1850  h2l[next_index++] = (2+i)*n2+(2+j)*n+2+k;
1851 
1852  Assert (next_index == fe.dofs_per_cell, ExcInternalError());
1853 
1854  break;
1855  }
1856 
1857  default:
1858  Assert (false, ExcNotImplemented());
1859  }
1860  return h2l;
1861 }
1862 
1863 
1864 template <int dim>
1865 std::vector<unsigned int>
1866 FE_Q_Hierarchical<dim>::
1867 face_fe_q_hierarchical_to_hierarchic_numbering (const unsigned int degree)
1868 {
1869  FiniteElementData<dim-1> fe_data(FE_Q_Hierarchical<dim-1>::get_dpo_vector(degree),1,degree);
1870  return invert_numbering(FE_Q_Hierarchical<dim-1>::
1871  hierarchic_to_fe_q_hierarchical_numbering (fe_data));
1872 }
1873 
1874 
1875 
1876 template <>
1877 std::vector<unsigned int>
1878 FE_Q_Hierarchical<1>::face_fe_q_hierarchical_to_hierarchic_numbering (const unsigned int)
1879 {
1880  return std::vector<unsigned int> ();
1881 }
1882 
1883 
1884 template <>
1885 bool
1886 FE_Q_Hierarchical<1>::has_support_on_face (const unsigned int shape_index,
1887  const unsigned int face_index) const
1888 {
1889  Assert (shape_index < this->dofs_per_cell,
1890  ExcIndexRange (shape_index, 0, this->dofs_per_cell));
1891  Assert (face_index < GeometryInfo<1>::faces_per_cell,
1892  ExcIndexRange (face_index, 0, GeometryInfo<1>::faces_per_cell));
1893 
1894 
1895  // in 1d, things are simple. since
1896  // there is only one degree of
1897  // freedom per vertex in this
1898  // class, the first is on vertex 0
1899  // (==face 0 in some sense), the
1900  // second on face 1:
1901  return (((shape_index == 0) && (face_index == 0)) ||
1902  ((shape_index == 1) && (face_index == 1)));
1903 }
1904 
1905 
1906 
1907 
1908 template <int dim>
1909 bool
1910 FE_Q_Hierarchical<dim>::has_support_on_face (const unsigned int shape_index,
1911  const unsigned int face_index) const
1912 {
1913  Assert (shape_index < this->dofs_per_cell,
1914  ExcIndexRange (shape_index, 0, this->dofs_per_cell));
1915  Assert (face_index < GeometryInfo<dim>::faces_per_cell,
1916  ExcIndexRange (face_index, 0, GeometryInfo<dim>::faces_per_cell));
1917 
1918  // first, special-case interior
1919  // shape functions, since they
1920  // have no support no-where on
1921  // the boundary
1922  if (((dim==2) && (shape_index>=this->first_quad_index))
1923  ||
1924  ((dim==3) && (shape_index>=this->first_hex_index)))
1925  return false;
1926 
1927  // let's see whether this is a
1928  // vertex
1929  if (shape_index < this->first_line_index)
1930  {
1931  // for Q elements, there is
1932  // one dof per vertex, so
1933  // shape_index==vertex_number. check
1934  // whether this vertex is
1935  // on the given face.
1936  const unsigned int vertex_no = shape_index;
1937  Assert (vertex_no < GeometryInfo<dim>::vertices_per_cell,
1938  ExcInternalError());
1939  for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_face; ++i)
1940  if (GeometryInfo<dim>::face_to_cell_vertices(face_index,i) == vertex_no)
1941  return true;
1942  return false;
1943  }
1944  else if (shape_index < this->first_quad_index)
1945  // ok, dof is on a line
1946  {
1947  const unsigned int line_index
1948  = (shape_index - this->first_line_index) / this->dofs_per_line;
1949  Assert (line_index < GeometryInfo<dim>::lines_per_cell,
1950  ExcInternalError());
1951 
1952  for (unsigned int i=0; i<GeometryInfo<dim>::lines_per_face; ++i)
1953  if (GeometryInfo<dim>::face_to_cell_lines(face_index,i) == line_index)
1954  return true;
1955  return false;
1956  }
1957  else if (shape_index < this->first_hex_index)
1958  // dof is on a quad
1959  {
1960  const unsigned int quad_index
1961  = (shape_index - this->first_quad_index) / this->dofs_per_quad;
1962  Assert (static_cast<signed int>(quad_index) <
1963  static_cast<signed int>(GeometryInfo<dim>::quads_per_cell),
1964  ExcInternalError());
1965 
1966  // in 2d, cell bubble are
1967  // zero on all faces. but
1968  // we have treated this
1969  // case above already
1970  Assert (dim != 2, ExcInternalError());
1971 
1972  // in 3d,
1973  // quad_index=face_index
1974  if (dim == 3)
1975  return (quad_index == face_index);
1976  else
1977  Assert (false, ExcNotImplemented());
1978  }
1979  else
1980  // dof on hex
1981  {
1982  // can only happen in 3d, but
1983  // this case has already been
1984  // covered above
1985  Assert (false, ExcNotImplemented());
1986  return false;
1987  }
1988 
1989  // we should not have gotten here
1990  Assert (false, ExcInternalError());
1991  return false;
1992 }
1993 
1994 
1995 
1996 template <int dim>
1997 std::vector<unsigned int>
1998 FE_Q_Hierarchical<dim>::get_embedding_dofs (const unsigned int sub_degree) const
1999 {
2000  Assert ((sub_degree>0) && (sub_degree<=this->degree),
2001  ExcIndexRange(sub_degree, 1, this->degree));
2002 
2003  if (dim==1)
2004  {
2005  std::vector<unsigned int> embedding_dofs (sub_degree+1);
2006  for (unsigned int i=0; i<(sub_degree+1); ++i)
2007  embedding_dofs[i] = i;
2008 
2009  return embedding_dofs;
2010  }
2011 
2012  if (sub_degree==1)
2013  {
2014  std::vector<unsigned int> embedding_dofs (GeometryInfo<dim>::vertices_per_cell);
2015  for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
2016  embedding_dofs[i] = i;
2017 
2018  return embedding_dofs;
2019  }
2020  else if (sub_degree==this->degree)
2021  {
2022  std::vector<unsigned int> embedding_dofs (this->dofs_per_cell);
2023  for (unsigned int i=0; i<this->dofs_per_cell; ++i)
2024  embedding_dofs[i] = i;
2025 
2026  return embedding_dofs;
2027  }
2028 
2029  if ((dim==2) || (dim==3))
2030  {
2031  std::vector<unsigned int> embedding_dofs ( (dim==2) ?
2032  (sub_degree+1) * (sub_degree+1) :
2033  (sub_degree+1) * (sub_degree+1) * (sub_degree+1) );
2034 
2035  for (unsigned int i=0; i<( (dim==2) ?
2036  (sub_degree+1) * (sub_degree+1) :
2037  (sub_degree+1) * (sub_degree+1) * (sub_degree+1) ); ++i)
2038  {
2039  // vertex
2040  if (i<GeometryInfo<dim>::vertices_per_cell)
2041  embedding_dofs[i] = i;
2042  // line
2043  else if (i<(GeometryInfo<dim>::vertices_per_cell +
2044  GeometryInfo<dim>::lines_per_cell * (sub_degree-1)))
2045  {
2046  const unsigned int j = (i - GeometryInfo<dim>::vertices_per_cell) %
2047  (sub_degree-1);
2048  const unsigned int line = (i - GeometryInfo<dim>::vertices_per_cell - j) / (sub_degree-1);
2049 
2050  embedding_dofs[i] = GeometryInfo<dim>::vertices_per_cell +
2051  line * (this->degree-1) + j;
2052  }
2053  // quad
2054  else if (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))
2057  {
2058  const unsigned int j = (i - GeometryInfo<dim>::vertices_per_cell -
2059  GeometryInfo<dim>::lines_per_cell * (sub_degree-1)) % (sub_degree-1);
2060  const unsigned int k = ( (i - GeometryInfo<dim>::vertices_per_cell -
2061  GeometryInfo<dim>::lines_per_cell * (sub_degree-1) - j) / (sub_degree-1) ) % (sub_degree-1);
2062  const unsigned int face = (i - GeometryInfo<dim>::vertices_per_cell -
2063  GeometryInfo<dim>::lines_per_cell * (sub_degree-1) - k * (sub_degree-1) - j) / ( (sub_degree-1) * (sub_degree-1) );
2064 
2065  embedding_dofs[i] = GeometryInfo<dim>::vertices_per_cell +
2066  GeometryInfo<dim>::lines_per_cell * (this->degree-1) +
2067  face * (this->degree-1) * (this->degree-1) +
2068  k * (this->degree-1) + j;
2069  }
2070  // hex
2071  else if (i<(GeometryInfo<dim>::vertices_per_cell +
2072  GeometryInfo<dim>::lines_per_cell * (sub_degree-1)) +
2073  GeometryInfo<dim>::quads_per_cell * (sub_degree-1) * (sub_degree-1) +
2074  GeometryInfo<dim>::hexes_per_cell * (sub_degree-1) * (sub_degree-1) * (sub_degree-1))
2075  {
2076  const unsigned int j = (i - GeometryInfo<dim>::vertices_per_cell -
2077  GeometryInfo<dim>::lines_per_cell * (sub_degree-1) -
2078  GeometryInfo<dim>::quads_per_cell * (sub_degree-1) * (sub_degree-1) ) % (sub_degree-1);
2079  const unsigned int k = ( (i - GeometryInfo<dim>::vertices_per_cell -
2080  GeometryInfo<dim>::lines_per_cell * (sub_degree-1) -
2081  GeometryInfo<dim>::quads_per_cell * (sub_degree-1) * (sub_degree-1) - j) / (sub_degree-1) ) % (sub_degree-1);
2082  const unsigned int l = (i - GeometryInfo<dim>::vertices_per_cell -
2083  GeometryInfo<dim>::lines_per_cell * (sub_degree-1) -
2084  GeometryInfo<dim>::quads_per_cell * (sub_degree-1) * (sub_degree-1) - j - k * (sub_degree-1)) / ( (sub_degree-1) * (sub_degree-1) );
2085 
2086  embedding_dofs[i] = GeometryInfo<dim>::vertices_per_cell +
2087  GeometryInfo<dim>::lines_per_cell * (this->degree-1) +
2088  GeometryInfo<dim>::quads_per_cell * (this->degree-1) * (this->degree-1) +
2089  l * (this->degree-1) * (this->degree-1) + k * (this->degree-1) + j;
2090  }
2091  }
2092 
2093  return embedding_dofs;
2094  }
2095  else
2096  {
2097  Assert(false, ExcNotImplemented ());
2098  return std::vector<unsigned int> ();
2099  }
2100 }
2101 
2102 
2103 
2104 template <int dim>
2105 std::pair<Table<2,bool>, std::vector<unsigned int> >
2106 FE_Q_Hierarchical<dim>::get_constant_modes () const
2107 {
2108  Table<2,bool> constant_modes(1, this->dofs_per_cell);
2109  for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
2110  constant_modes(0,i) = true;
2111  for (unsigned int i=GeometryInfo<dim>::vertices_per_cell; i<this->dofs_per_cell; ++i)
2112  constant_modes(0,i) = false;
2113  return std::pair<Table<2,bool>, std::vector<unsigned int> >
2114  (constant_modes, std::vector<unsigned int>(1, 0));
2115 }
2116 
2117 
2118 
2119 template <int dim>
2120 std::size_t
2121 FE_Q_Hierarchical<dim>::memory_consumption () const
2122 {
2123  Assert (false, ExcNotImplemented ());
2124  return 0;
2125 }
2126 
2127 
2128 
2129 // explicit instantiations
2130 #include "fe_q_hierarchical.inst"
2131 
2132 
2133 DEAL_II_NAMESPACE_CLOSE
FiniteElementData::first_hex_index
const unsigned int first_hex_index
Definition: fe_base.h:259
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
FiniteElement< dim, dim >::restriction
std::vector< std::vector< FullMatrix< double > > > restriction
Definition: fe.h:2106
FiniteElement< dim, dim >::interface_constraints
FullMatrix< double > interface_constraints
Definition: fe.h:2132
FE_Q_Hierarchical::hp_constraints_are_implemented
virtual bool hp_constraints_are_implemented() const
Definition: fe_q_hierarchical.cc:232
StandardExceptions::ExcMessage
::ExceptionBase & ExcMessage(std::string arg1)
FiniteElementData::dofs_per_quad
const unsigned int dofs_per_quad
Definition: fe_base.h:238
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:640
FiniteElementData::degree
const unsigned int degree
Definition: fe_base.h:299
FE_Q_Hierarchical::initialize_unit_support_points
void initialize_unit_support_points()
Definition: fe_q_hierarchical.cc:743
TensorProductPolynomials::set_numbering
void set_numbering(const std::vector< unsigned int > &renumber)
Definition: tensor_product_polynomials.cc:110
FE_Q_Hierarchical::clone
virtual FiniteElement< dim > * clone() const
Definition: fe_q_hierarchical.cc:162
std
STL namespace.
FE_Q_Hierarchical::face_renumber
const std::vector< unsigned int > face_renumber
Definition: fe_q_hierarchical.h:809
AssertThrow
#define AssertThrow(cond, exc)
Definition: exceptions.h:358
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
FiniteElementData::dofs_per_line
const unsigned int dofs_per_line
Definition: fe_base.h:232
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:1998
FE_Q_Hierarchical::get_name
virtual std::string get_name() const
Definition: fe_q_hierarchical.cc:110
FiniteElement< dim, dim >::unit_support_points
std::vector< Point< dim > > unit_support_points
Definition: fe.h:2144
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:1910
FiniteElementData::first_quad_index
const unsigned int first_quad_index
Definition: fe_base.h:254
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:1867
FullMatrix::n
size_type n() const
FiniteElement< dim, dim >::prolongation
std::vector< std::vector< FullMatrix< double > > > prolongation
Definition: fe.h:2120
FiniteElementData::dofs_per_hex
const unsigned int dofs_per_hex
Definition: fe_base.h:244
Assert
#define Assert(cond, exc)
Definition: exceptions.h:294
FiniteElement< dim, dim >::unit_face_support_points
std::vector< Point< dim-1 > > unit_face_support_points
Definition: fe.h:2151
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:2106
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:544
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:936
TensorProductPolynomials::get_numbering
const std::vector< unsigned int > & get_numbering() const
FiniteElementData::dofs_per_cell
const unsigned int dofs_per_cell
Definition: fe_base.h:283
FE_Q_Hierarchical::FE_Q_Hierarchical
friend class FE_Q_Hierarchical
Definition: fe_q_hierarchical.h:817
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:1693
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
FE_Q_Hierarchical::interpolate
virtual void interpolate(std::vector< double > &local_dofs, const std::vector< double > &values) const
Definition: fe_q_hierarchical.cc:128
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:276
TensorProductPolynomials::get_numbering_inverse
const std::vector< unsigned int > & get_numbering_inverse() const
FiniteElementData::first_line_index
const unsigned int first_line_index
Definition: fe_base.h:249
FE_Q_Hierarchical::initialize_unit_face_support_points
void initialize_unit_face_support_points()
Definition: fe_q_hierarchical.cc:1630
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:849
FiniteElementData::dofs_per_vertex
const unsigned int dofs_per_vertex
Definition: fe_base.h:226
FE_Q_Hierarchical::memory_consumption
virtual std::size_t memory_consumption() const
Definition: fe_q_hierarchical.cc:2121
FE_Q_Hierarchical::get_dpo_vector
static std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
Definition: fe_q_hierarchical.cc:1680
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
FiniteElement< dim, dim >::interface_constraints_size
TableIndices< 2 > interface_constraints_size() const
Definition: fe.cc:824
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
Generated by   doxygen 1.8.12