]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Single out the functions that initialize matrices from the
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 23 Apr 2003 22:59:17 +0000 (22:59 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 23 Apr 2003 22:59:17 +0000 (22:59 +0000)
constructor, to make the whole process of initializing finite element
classes a little more transparent. Also makes it simple to replace
precomputed matrices with something that computes these values on the fly.

git-svn-id: https://svn.dealii.org/trunk@7463 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/fe/fe_q.h
deal.II/deal.II/source/fe/fe_q.cc

index cb4ec781fd06577ecda010ba9a9d58d5ea76211b..218949d220d0dadcad8b64c004d35d3096785e0d 100644 (file)
@@ -624,6 +624,9 @@ class FE_Q : public FiniteElement<dim>
                                      * chosen to weigh the simplicity
                                      * aspect a little more than
                                      * proper design.
+                                     *
+                                     * This function is called from
+                                     * the constructor.
                                      */
     static
     std::vector<unsigned int>
@@ -633,12 +636,34 @@ class FE_Q : public FiniteElement<dim>
                                     /**
                                      * This is an analogon to the
                                      * previous function, but working
-                                     * on faces.
+                                     * on faces. Called from the
+                                     * constructor.
                                      */
     static
     std::vector<unsigned int>
     face_lexicographic_to_hierarchic_numbering (const unsigned int degree);
 
+                                    /**
+                                     * Initialize the hanging node
+                                     * constraints matrices. Called
+                                     * from the constructor.
+                                     */
+    void initialize_constraints ();
+
+                                    /**
+                                     * Initialize the embedding
+                                     * matrices. Called from the
+                                     * constructor.
+                                     */
+    void initialize_embedding ();
+
+                                    /**
+                                     * Initialize the restriction
+                                     * matrices. Called from the
+                                     * constructor.
+                                     */
+    void initialize_restriction ();
+
                                     /**
                                      * Initialize the
                                      * @p{unit_support_points} field
index b39b40161601626005f13878a0899a97d2d169ec..8a85752502da8824d183e8a3df4e438e0de94c50 100644 (file)
@@ -51,779 +51,213 @@ FE_Q<dim>::FE_Q (const unsigned int degree)
                face_renumber(face_lexicographic_to_hierarchic_numbering (degree)),
                polynomial_space(Polynomials::LagrangeEquidistant::generate_complete_basis(degree))
 {
-                                  // copy constraint matrices if they
-                                  // are defined. otherwise leave them
-                                  // at invalid size
-  if ((dim>1) && (degree<Matrices::n_constraint_matrices+1))
+  
+                                  // copy constraint and embedding
+                                  // matrices if they are
+                                  // defined. otherwise leave them at
+                                  // invalid size
+  initialize_constraints ();
+  initialize_embedding ();
+  initialize_restriction ();
+
+                                  // finally fill in support points
+                                  // on cell and face
+  initialize_unit_support_points ();
+  initialize_unit_face_support_points ();
+}
+
+
+
+template <int dim>
+FiniteElement<dim> *
+FE_Q<dim>::clone() const
+{
+  return new FE_Q<dim>(degree);
+}
+
+
+
+template <int dim>
+double
+FE_Q<dim>::shape_value (const unsigned int i,
+                       const Point<dim> &p) const
+{
+  Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
+  return polynomial_space.compute_value(renumber_inverse[i], p);
+}
+
+
+template <int dim>
+double
+FE_Q<dim>::shape_value_component (const unsigned int i,
+                                 const Point<dim> &p,
+                                 const unsigned int component) const
+{
+  Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
+  Assert (component == 0, ExcIndexRange (component, 0, 1));
+  return polynomial_space.compute_value(renumber_inverse[i], p);
+}
+
+
+
+template <int dim>
+Tensor<1,dim>
+FE_Q<dim>::shape_grad (const unsigned int i,
+                      const Point<dim> &p) const
+{
+  Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
+  return polynomial_space.compute_grad(renumber_inverse[i], p);
+}
+
+
+
+template <int dim>
+Tensor<1,dim>
+FE_Q<dim>::shape_grad_component (const unsigned int i,
+                                const Point<dim> &p,
+                                const unsigned int component) const
+{
+  Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
+  Assert (component == 0, ExcIndexRange (component, 0, 1));
+  return polynomial_space.compute_grad(renumber_inverse[i], p);
+}
+
+
+
+template <int dim>
+Tensor<2,dim>
+FE_Q<dim>::shape_grad_grad (const unsigned int i,
+                           const Point<dim> &p) const
+{
+  Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
+  return polynomial_space.compute_grad_grad(renumber_inverse[i], p);
+}
+
+
+
+template <int dim>
+Tensor<2,dim>
+FE_Q<dim>::shape_grad_grad_component (const unsigned int i,
+                                     const Point<dim> &p,
+                                     const unsigned int component) const
+{
+  Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
+  Assert (component == 0, ExcIndexRange (component, 0, 1));
+  return polynomial_space.compute_grad_grad(renumber_inverse[i], p);
+}
+
+
+//----------------------------------------------------------------------
+// Auxiliary functions
+//----------------------------------------------------------------------
+
+
+
+template <int dim>
+void FE_Q<dim>::initialize_unit_support_points ()
+{
+                                  // number of points: (degree+1)^dim
+  unsigned int n = degree+1;
+  for (unsigned int i=1; i<dim; ++i)
+    n *= degree+1;
+  
+  this->unit_support_points.resize(n);
+  
+  const double step = 1./degree;
+  Point<dim> p;
+  
+  unsigned int k=0;
+  for (unsigned int iz=0; iz <= ((dim>2) ? degree : 0) ; ++iz)
+    for (unsigned int iy=0; iy <= ((dim>1) ? degree : 0) ; ++iy)
+      for (unsigned int ix=0; ix<=degree; ++ix)
+       {
+         p(0) = ix * step;
+         if (dim>1)
+           p(1) = iy * step;
+         if (dim>2)
+           p(2) = iz * step;
+         
+         this->unit_support_points[renumber[k++]] = p;
+       };
+}
+
+
+#if deal_II_dimension == 1
+
+template <>
+void FE_Q<1>::initialize_unit_face_support_points ()
+{
+                                  // no faces in 1d, so nothing to do
+}
+
+#endif
+
+
+template <int dim>
+void FE_Q<dim>::initialize_unit_face_support_points ()
+{
+  const unsigned int codim = dim-1;
+  
+                                  // number of points: (degree+1)^codim
+  unsigned int n = degree+1;
+  for (unsigned int i=1; i<codim; ++i)
+    n *= degree+1;
+  
+  this->unit_face_support_points.resize(n);
+  
+  const double step = 1./degree;
+  Point<codim> p;
+  
+  unsigned int k=0;
+  for (unsigned int iz=0; iz <= ((codim>2) ? degree : 0) ; ++iz)
+    for (unsigned int iy=0; iy <= ((codim>1) ? degree : 0) ; ++iy)
+      for (unsigned int ix=0; ix<=degree; ++ix)
+       {
+         p(0) = ix * step;
+         if (codim>1)
+           p(1) = iy * step;
+         if (codim>2)
+           p(2) = iz * step;
+         
+         this->unit_face_support_points[face_renumber[k++]] = p;
+       };
+}
+
+
+
+template <int dim>
+std::vector<unsigned int>
+FE_Q<dim>::get_dpo_vector(const unsigned int deg)
+{
+  std::vector<unsigned int> dpo(dim+1, static_cast<unsigned int>(1));
+  for (unsigned int i=1; i<dpo.size(); ++i)
+    dpo[i]=dpo[i-1]*(deg-1);
+  return dpo;
+}
+
+
+
+template <int dim>
+std::vector<unsigned int>
+FE_Q<dim>::lexicographic_to_hierarchic_numbering (const FiniteElementData<dim> &fe_data,
+                                                 const unsigned int            degree)
+{
+  std::vector<unsigned int> renumber (fe_data.dofs_per_cell);
+  
+  const unsigned int n = degree+1;
+
+
+  if (degree == 0)
     {
-      this->interface_constraints.
-        TableBase<2,double>::reinit (this->interface_constraints_size());
-      this->interface_constraints.fill (Matrices::constraint_matrices[degree-1]);
+      Assert ((fe_data.dofs_per_vertex == 0) &&
+             ((fe_data.dofs_per_line == 0) || (dim == 1)) &&
+             ((fe_data.dofs_per_quad == 0) || (dim == 2)) &&
+             ((fe_data.dofs_per_hex == 0)  || (dim == 3)),
+             ExcInternalError());
+      renumber[0] = 0;
     };
 
-                                  // next copy over embedding
-                                  // matrices if they are defined
-  if ((degree < Matrices::n_embedding_matrices+1) &&
-      (Matrices::embedding[degree-1][0] != 0))
-    for (unsigned int c=0; c<GeometryInfo<dim>::children_per_cell; ++c)
-      {
-        this->prolongation[c].reinit (this->dofs_per_cell,
-                                      this->dofs_per_cell);
-        this->prolongation[c].fill (Matrices::embedding[degree-1][c]);
-      };
-
-                                  // then fill restriction
-                                  // matrices. they are hardcoded for
-                                  // the first few elements. in
-                                  // contrast to the other matrices,
-                                  // these are not stored in the
-                                  // files fe_q_[123]d.cc, since they
-                                  // contain only a rather small
-                                  // number of zeros, and storing
-                                  // them element-wise is more
-                                  // expensive than just setting the
-                                  // nonzero elements as done here
-  for (unsigned int c=0; c<GeometryInfo<dim>::children_per_cell; ++c)
-    this->restriction[c].reinit (this->dofs_per_cell, this->dofs_per_cell);
-  switch (dim)
-    {
-      case 1:  // 1d
-      {
-        switch (degree)
-          {
-            case 1:
-                  this->restriction[0](0,0) = 1;
-                  this->restriction[1](1,1) = 1;
-                  break;
-            case 2:
-                  this->restriction[0](0,0) = 1;
-                  this->restriction[0](2,1) = 1;
-                  this->restriction[1](1,1) = 1;
-                  this->restriction[1](1,1) = 1;
-                  break;
-            case 3:
-                  this->restriction[0](0,0) = 1;
-                  this->restriction[0](2,3) = 1;
-                  this->restriction[1](1,1) = 1;
-                  this->restriction[1](3,2) = 1;
-                  break;
-            case 4:
-                  this->restriction[0](0,0) = 1;
-                  this->restriction[0](2,3) = 1;
-                  this->restriction[0](3,1) = 1;
-                  this->restriction[1](1,1) = 1;
-                  this->restriction[1](3,0) = 1;
-                  this->restriction[1](4,3) = 1;
-                  break;
-                 
-            default:
-            {
-                                               // in case we don't
-                                               // have the matrices
-                                               // (yet), reset them to
-                                               // zero size. this does
-                                               // not prevent the use
-                                               // of this FE, but will
-                                               // prevent the use of
-                                               // these matrices
-              for (unsigned int i=0;
-                   i<GeometryInfo<dim>::children_per_cell;
-                   ++i)
-                this->restriction[i].reinit(0,0);
-            };
-          }
-        break;
-      };
-      
-      case 2:  // 2d
-      {
-        switch (degree)
-          {
-            case 1:
-                  this->restriction[0](0,0) = 1;
-                  this->restriction[1](1,1) = 1;
-                  this->restriction[2](2,2) = 1;
-                  this->restriction[3](3,3) = 1;
-                  break;
-            case 2:
-                  this->restriction[0](0,0) = 1;
-                  this->restriction[0](4,1) = 1;
-                  this->restriction[0](7,3) = 1;
-                  this->restriction[0](8,2) = 1;
-                  this->restriction[1](1,1) = 1;
-                  this->restriction[1](4,0) = 1;
-                  this->restriction[1](5,2) = 1;
-                  this->restriction[1](8,3) = 1;
-                  this->restriction[2](2,2) = 1;
-                  this->restriction[2](5,1) = 1;
-                  this->restriction[2](6,3) = 1;
-                  this->restriction[2](8,0) = 1;
-                  this->restriction[3](3,3) = 1;
-                  this->restriction[3](6,2) = 1;
-                  this->restriction[3](7,0) = 1;
-                  this->restriction[3](8,1) = 1;
-                  break;
-            case 3:
-                  this->restriction[0](0,0) = 1;
-                  this->restriction[0](4,5) = 1;
-                  this->restriction[0](10,11) = 1;
-                  this->restriction[0](12,15) = 1;
-                  this->restriction[1](1,1) = 1;
-                  this->restriction[1](5,4) = 1;
-                  this->restriction[1](6,7) = 1;
-                  this->restriction[1](13,14) = 1;
-                  this->restriction[2](2,2) = 1;
-                  this->restriction[2](7,6) = 1;
-                  this->restriction[2](9,8) = 1;
-                  this->restriction[2](15,12) = 1;
-                  this->restriction[3](3,3) = 1;
-                  this->restriction[3](8,9) = 1;
-                  this->restriction[3](11,10) = 1;
-                  this->restriction[3](14,13) = 1;
-                  break;
-            case 4:
-                  this->restriction[0](0,0) = 1;
-                  this->restriction[0](4,5) = 1;
-                  this->restriction[0](5,1) = 1;
-                  this->restriction[0](13,14) = 1;
-                  this->restriction[0](14,3) = 1;
-                  this->restriction[0](16,20) = 1;
-                  this->restriction[0](17,8) = 1;
-                  this->restriction[0](19,11) = 1;
-                  this->restriction[0](20,2) = 1;
-                  this->restriction[1](1,1) = 1;
-                  this->restriction[1](5,0) = 1;
-                  this->restriction[1](6,5) = 1;
-                  this->restriction[1](7,8) = 1;
-                  this->restriction[1](8,2) = 1;
-                  this->restriction[1](17,14) = 1;
-                  this->restriction[1](18,20) = 1;
-                  this->restriction[1](20,3) = 1;
-                  this->restriction[1](21,11) = 1;
-                  this->restriction[2](2,2) = 1;
-                  this->restriction[2](8,1) = 1;
-                  this->restriction[2](9,8) = 1;
-                  this->restriction[2](11,3) = 1;
-                  this->restriction[2](12,11) = 1;
-                  this->restriction[2](20,0) = 1;
-                  this->restriction[2](21,5) = 1;
-                  this->restriction[2](23,14) = 1;
-                  this->restriction[2](24,20) = 1;
-                  this->restriction[3](3,3) = 1;
-                  this->restriction[3](10,11) = 1;
-                  this->restriction[3](11,2) = 1;
-                  this->restriction[3](14,0) = 1;
-                  this->restriction[3](15,14) = 1;
-                  this->restriction[3](19,5) = 1;
-                  this->restriction[3](20,1) = 1;
-                  this->restriction[3](22,20) = 1;
-                  this->restriction[3](23,8) = 1;
-                  break;
-
-            default:
-            {
-                                               // in case we don't
-                                               // have the matrices
-                                               // (yet), reset them to
-                                               // zero size. this does
-                                               // not prevent the use
-                                               // of this FE, but will
-                                               // prevent the use of
-                                               // these matrices
-              for (unsigned int i=0;
-                   i<GeometryInfo<dim>::children_per_cell;
-                   ++i)
-                this->restriction[i].reinit(0,0);
-            };
-          }
-        break;
-      };
-      
-      case 3:  // 3d
-      {
-        switch (degree)
-          {
-            case 1:
-                  this->restriction[0](0,0) = 1;
-                  this->restriction[1](1,1) = 1;
-                  this->restriction[2](2,2) = 1;
-                  this->restriction[3](3,3) = 1;
-                  this->restriction[4](4,4) = 1;
-                  this->restriction[5](5,5) = 1;
-                  this->restriction[6](6,6) = 1;
-                  this->restriction[7](7,7) = 1;
-                  break;
-            case 2:
-                  this->restriction[0](0,0) = 1;
-                  this->restriction[0](8,1) = 1;
-                  this->restriction[0](11,3) = 1;
-                  this->restriction[0](16,4) = 1;
-                  this->restriction[0](20,2) = 1;
-                  this->restriction[0](22,5) = 1;
-                  this->restriction[0](25,7) = 1;
-                  this->restriction[0](26,6) = 1;
-                  this->restriction[1](1,1) = 1;
-                  this->restriction[1](8,0) = 1;
-                  this->restriction[1](9,2) = 1;
-                  this->restriction[1](17,5) = 1;
-                  this->restriction[1](20,3) = 1;
-                  this->restriction[1](22,4) = 1;
-                  this->restriction[1](23,6) = 1;
-                  this->restriction[1](26,7) = 1;
-                  this->restriction[2](2,2) = 1;
-                  this->restriction[2](9,1) = 1;
-                  this->restriction[2](10,3) = 1;
-                  this->restriction[2](18,6) = 1;
-                  this->restriction[2](20,0) = 1;
-                  this->restriction[2](23,5) = 1;
-                  this->restriction[2](24,7) = 1;
-                  this->restriction[2](26,4) = 1;
-                  this->restriction[3](3,3) = 1;
-                  this->restriction[3](10,2) = 1;
-                  this->restriction[3](11,0) = 1;
-                  this->restriction[3](19,7) = 1;
-                  this->restriction[3](20,1) = 1;
-                  this->restriction[3](24,6) = 1;
-                  this->restriction[3](25,4) = 1;
-                  this->restriction[3](26,5) = 1;
-                  this->restriction[4](4,4) = 1;
-                  this->restriction[4](12,5) = 1;
-                  this->restriction[4](15,7) = 1;
-                  this->restriction[4](16,0) = 1;
-                  this->restriction[4](21,6) = 1;
-                  this->restriction[4](22,1) = 1;
-                  this->restriction[4](25,3) = 1;
-                  this->restriction[4](26,2) = 1;
-                  this->restriction[5](5,5) = 1;
-                  this->restriction[5](12,4) = 1;
-                  this->restriction[5](13,6) = 1;
-                  this->restriction[5](17,1) = 1;
-                  this->restriction[5](21,7) = 1;
-                  this->restriction[5](22,0) = 1;
-                  this->restriction[5](23,2) = 1;
-                  this->restriction[5](26,3) = 1;
-                  this->restriction[6](6,6) = 1;
-                  this->restriction[6](13,5) = 1;
-                  this->restriction[6](14,7) = 1;
-                  this->restriction[6](18,2) = 1;
-                  this->restriction[6](21,4) = 1;
-                  this->restriction[6](23,1) = 1;
-                  this->restriction[6](24,3) = 1;
-                  this->restriction[6](26,0) = 1;
-                  this->restriction[7](7,7) = 1;
-                  this->restriction[7](14,6) = 1;
-                  this->restriction[7](15,4) = 1;
-                  this->restriction[7](19,3) = 1;
-                  this->restriction[7](21,5) = 1;
-                  this->restriction[7](24,2) = 1;
-                  this->restriction[7](25,0) = 1;
-                  this->restriction[7](26,1) = 1;
-                  break;
-            case 3:
-                  this->restriction[0](0,0) = 1;
-                  this->restriction[0](8,9) = 1;
-                  this->restriction[0](14,15) = 1;
-                  this->restriction[0](24,25) = 1;
-                  this->restriction[0](32,35) = 1;
-                  this->restriction[0](40,43) = 1;
-                  this->restriction[0](52,55) = 1;
-                  this->restriction[0](56,63) = 1;
-                  this->restriction[1](1,1) = 1;
-                  this->restriction[1](9,8) = 1;
-                  this->restriction[1](10,11) = 1;
-                  this->restriction[1](26,27) = 1;
-                  this->restriction[1](33,34) = 1;
-                  this->restriction[1](41,42) = 1;
-                  this->restriction[1](44,47) = 1;
-                  this->restriction[1](57,62) = 1;
-                  this->restriction[2](2,2) = 1;
-                  this->restriction[2](11,10) = 1;
-                  this->restriction[2](13,12) = 1;
-                  this->restriction[2](28,29) = 1;
-                  this->restriction[2](35,32) = 1;
-                  this->restriction[2](46,45) = 1;
-                  this->restriction[2](49,50) = 1;
-                  this->restriction[2](61,58) = 1;
-                  this->restriction[3](3,3) = 1;
-                  this->restriction[3](12,13) = 1;
-                  this->restriction[3](15,14) = 1;
-                  this->restriction[3](30,31) = 1;
-                  this->restriction[3](34,33) = 1;
-                  this->restriction[3](48,51) = 1;
-                  this->restriction[3](54,53) = 1;
-                  this->restriction[3](60,59) = 1;
-                  this->restriction[4](4,4) = 1;
-                  this->restriction[4](16,17) = 1;
-                  this->restriction[4](22,23) = 1;
-                  this->restriction[4](25,24) = 1;
-                  this->restriction[4](36,39) = 1;
-                  this->restriction[4](42,41) = 1;
-                  this->restriction[4](53,54) = 1;
-                  this->restriction[4](58,61) = 1;
-                  this->restriction[5](5,5) = 1;
-                  this->restriction[5](17,16) = 1;
-                  this->restriction[5](18,19) = 1;
-                  this->restriction[5](27,26) = 1;
-                  this->restriction[5](37,38) = 1;
-                  this->restriction[5](43,40) = 1;
-                  this->restriction[5](45,46) = 1;
-                  this->restriction[5](59,60) = 1;
-                  this->restriction[6](6,6) = 1;
-                  this->restriction[6](19,18) = 1;
-                  this->restriction[6](21,20) = 1;
-                  this->restriction[6](29,28) = 1;
-                  this->restriction[6](39,36) = 1;
-                  this->restriction[6](47,44) = 1;
-                  this->restriction[6](51,48) = 1;
-                  this->restriction[6](63,56) = 1;
-                  this->restriction[7](7,7) = 1;
-                  this->restriction[7](20,21) = 1;
-                  this->restriction[7](23,22) = 1;
-                  this->restriction[7](31,30) = 1;
-                  this->restriction[7](38,37) = 1;
-                  this->restriction[7](50,49) = 1;
-                  this->restriction[7](55,52) = 1;
-                  this->restriction[7](62,57) = 1;
-                  break;
-            case 4:
-                  this->restriction[0](0,0) = 1;
-                  this->restriction[0](8,9) = 1;
-                  this->restriction[0](9,1) = 1;
-                  this->restriction[0](17,18) = 1;
-                  this->restriction[0](18,3) = 1;
-                  this->restriction[0](32,33) = 1;
-                  this->restriction[0](33,4) = 1;
-                  this->restriction[0](44,48) = 1;
-                  this->restriction[0](45,12) = 1;
-                  this->restriction[0](47,15) = 1;
-                  this->restriction[0](48,2) = 1;
-                  this->restriction[0](62,66) = 1;
-                  this->restriction[0](63,36) = 1;
-                  this->restriction[0](65,21) = 1;
-                  this->restriction[0](66,5) = 1;
-                  this->restriction[0](89,93) = 1;
-                  this->restriction[0](90,30) = 1;
-                  this->restriction[0](92,42) = 1;
-                  this->restriction[0](93,7) = 1;
-                  this->restriction[0](98,111) = 1;
-                  this->restriction[0](99,75) = 1;
-                  this->restriction[0](101,57) = 1;
-                  this->restriction[0](102,24) = 1;
-                  this->restriction[0](107,84) = 1;
-                  this->restriction[0](108,39) = 1;
-                  this->restriction[0](110,27) = 1;
-                  this->restriction[0](111,6) = 1;
-                  this->restriction[1](1,1) = 1;
-                  this->restriction[1](9,0) = 1;
-                  this->restriction[1](10,9) = 1;
-                  this->restriction[1](11,12) = 1;
-                  this->restriction[1](12,2) = 1;
-                  this->restriction[1](35,36) = 1;
-                  this->restriction[1](36,5) = 1;
-                  this->restriction[1](45,18) = 1;
-                  this->restriction[1](46,48) = 1;
-                  this->restriction[1](48,3) = 1;
-                  this->restriction[1](49,15) = 1;
-                  this->restriction[1](63,33) = 1;
-                  this->restriction[1](64,66) = 1;
-                  this->restriction[1](66,4) = 1;
-                  this->restriction[1](67,21) = 1;
-                  this->restriction[1](71,75) = 1;
-                  this->restriction[1](72,24) = 1;
-                  this->restriction[1](74,39) = 1;
-                  this->restriction[1](75,6) = 1;
-                  this->restriction[1](99,93) = 1;
-                  this->restriction[1](100,111) = 1;
-                  this->restriction[1](102,30) = 1;
-                  this->restriction[1](103,57) = 1;
-                  this->restriction[1](108,42) = 1;
-                  this->restriction[1](109,84) = 1;
-                  this->restriction[1](111,7) = 1;
-                  this->restriction[1](112,27) = 1;
-                  this->restriction[2](2,2) = 1;
-                  this->restriction[2](12,1) = 1;
-                  this->restriction[2](13,12) = 1;
-                  this->restriction[2](15,3) = 1;
-                  this->restriction[2](16,15) = 1;
-                  this->restriction[2](38,39) = 1;
-                  this->restriction[2](39,6) = 1;
-                  this->restriction[2](48,0) = 1;
-                  this->restriction[2](49,9) = 1;
-                  this->restriction[2](51,18) = 1;
-                  this->restriction[2](52,48) = 1;
-                  this->restriction[2](74,36) = 1;
-                  this->restriction[2](75,5) = 1;
-                  this->restriction[2](77,75) = 1;
-                  this->restriction[2](78,24) = 1;
-                  this->restriction[2](81,42) = 1;
-                  this->restriction[2](82,84) = 1;
-                  this->restriction[2](84,7) = 1;
-                  this->restriction[2](85,27) = 1;
-                  this->restriction[2](108,33) = 1;
-                  this->restriction[2](109,66) = 1;
-                  this->restriction[2](111,4) = 1;
-                  this->restriction[2](112,21) = 1;
-                  this->restriction[2](117,93) = 1;
-                  this->restriction[2](118,111) = 1;
-                  this->restriction[2](120,30) = 1;
-                  this->restriction[2](121,57) = 1;
-                  this->restriction[3](3,3) = 1;
-                  this->restriction[3](14,15) = 1;
-                  this->restriction[3](15,2) = 1;
-                  this->restriction[3](18,0) = 1;
-                  this->restriction[3](19,18) = 1;
-                  this->restriction[3](41,42) = 1;
-                  this->restriction[3](42,7) = 1;
-                  this->restriction[3](47,9) = 1;
-                  this->restriction[3](48,1) = 1;
-                  this->restriction[3](50,48) = 1;
-                  this->restriction[3](51,12) = 1;
-                  this->restriction[3](80,84) = 1;
-                  this->restriction[3](81,39) = 1;
-                  this->restriction[3](83,27) = 1;
-                  this->restriction[3](84,6) = 1;
-                  this->restriction[3](92,33) = 1;
-                  this->restriction[3](93,4) = 1;
-                  this->restriction[3](95,93) = 1;
-                  this->restriction[3](96,30) = 1;
-                  this->restriction[3](107,66) = 1;
-                  this->restriction[3](108,36) = 1;
-                  this->restriction[3](110,21) = 1;
-                  this->restriction[3](111,5) = 1;
-                  this->restriction[3](116,111) = 1;
-                  this->restriction[3](117,75) = 1;
-                  this->restriction[3](119,57) = 1;
-                  this->restriction[3](120,24) = 1;
-                  this->restriction[4](4,4) = 1;
-                  this->restriction[4](20,21) = 1;
-                  this->restriction[4](21,5) = 1;
-                  this->restriction[4](29,30) = 1;
-                  this->restriction[4](30,7) = 1;
-                  this->restriction[4](33,0) = 1;
-                  this->restriction[4](34,33) = 1;
-                  this->restriction[4](53,57) = 1;
-                  this->restriction[4](54,24) = 1;
-                  this->restriction[4](56,27) = 1;
-                  this->restriction[4](57,6) = 1;
-                  this->restriction[4](65,9) = 1;
-                  this->restriction[4](66,1) = 1;
-                  this->restriction[4](68,66) = 1;
-                  this->restriction[4](69,36) = 1;
-                  this->restriction[4](90,18) = 1;
-                  this->restriction[4](91,93) = 1;
-                  this->restriction[4](93,3) = 1;
-                  this->restriction[4](94,42) = 1;
-                  this->restriction[4](101,48) = 1;
-                  this->restriction[4](102,12) = 1;
-                  this->restriction[4](104,111) = 1;
-                  this->restriction[4](105,75) = 1;
-                  this->restriction[4](110,15) = 1;
-                  this->restriction[4](111,2) = 1;
-                  this->restriction[4](113,84) = 1;
-                  this->restriction[4](114,39) = 1;
-                  this->restriction[5](5,5) = 1;
-                  this->restriction[5](21,4) = 1;
-                  this->restriction[5](22,21) = 1;
-                  this->restriction[5](23,24) = 1;
-                  this->restriction[5](24,6) = 1;
-                  this->restriction[5](36,1) = 1;
-                  this->restriction[5](37,36) = 1;
-                  this->restriction[5](54,30) = 1;
-                  this->restriction[5](55,57) = 1;
-                  this->restriction[5](57,7) = 1;
-                  this->restriction[5](58,27) = 1;
-                  this->restriction[5](66,0) = 1;
-                  this->restriction[5](67,9) = 1;
-                  this->restriction[5](69,33) = 1;
-                  this->restriction[5](70,66) = 1;
-                  this->restriction[5](72,12) = 1;
-                  this->restriction[5](73,75) = 1;
-                  this->restriction[5](75,2) = 1;
-                  this->restriction[5](76,39) = 1;
-                  this->restriction[5](102,18) = 1;
-                  this->restriction[5](103,48) = 1;
-                  this->restriction[5](105,93) = 1;
-                  this->restriction[5](106,111) = 1;
-                  this->restriction[5](111,3) = 1;
-                  this->restriction[5](112,15) = 1;
-                  this->restriction[5](114,42) = 1;
-                  this->restriction[5](115,84) = 1;
-                  this->restriction[6](6,6) = 1;
-                  this->restriction[6](24,5) = 1;
-                  this->restriction[6](25,24) = 1;
-                  this->restriction[6](27,7) = 1;
-                  this->restriction[6](28,27) = 1;
-                  this->restriction[6](39,2) = 1;
-                  this->restriction[6](40,39) = 1;
-                  this->restriction[6](57,4) = 1;
-                  this->restriction[6](58,21) = 1;
-                  this->restriction[6](60,30) = 1;
-                  this->restriction[6](61,57) = 1;
-                  this->restriction[6](75,1) = 1;
-                  this->restriction[6](76,36) = 1;
-                  this->restriction[6](78,12) = 1;
-                  this->restriction[6](79,75) = 1;
-                  this->restriction[6](84,3) = 1;
-                  this->restriction[6](85,15) = 1;
-                  this->restriction[6](87,42) = 1;
-                  this->restriction[6](88,84) = 1;
-                  this->restriction[6](111,0) = 1;
-                  this->restriction[6](112,9) = 1;
-                  this->restriction[6](114,33) = 1;
-                  this->restriction[6](115,66) = 1;
-                  this->restriction[6](120,18) = 1;
-                  this->restriction[6](121,48) = 1;
-                  this->restriction[6](123,93) = 1;
-                  this->restriction[6](124,111) = 1;
-                  this->restriction[7](7,7) = 1;
-                  this->restriction[7](26,27) = 1;
-                  this->restriction[7](27,6) = 1;
-                  this->restriction[7](30,4) = 1;
-                  this->restriction[7](31,30) = 1;
-                  this->restriction[7](42,3) = 1;
-                  this->restriction[7](43,42) = 1;
-                  this->restriction[7](56,21) = 1;
-                  this->restriction[7](57,5) = 1;
-                  this->restriction[7](59,57) = 1;
-                  this->restriction[7](60,24) = 1;
-                  this->restriction[7](83,15) = 1;
-                  this->restriction[7](84,2) = 1;
-                  this->restriction[7](86,84) = 1;
-                  this->restriction[7](87,39) = 1;
-                  this->restriction[7](93,0) = 1;
-                  this->restriction[7](94,33) = 1;
-                  this->restriction[7](96,18) = 1;
-                  this->restriction[7](97,93) = 1;
-                  this->restriction[7](110,9) = 1;
-                  this->restriction[7](111,1) = 1;
-                  this->restriction[7](113,66) = 1;
-                  this->restriction[7](114,36) = 1;
-                  this->restriction[7](119,48) = 1;
-                  this->restriction[7](120,12) = 1;
-                  this->restriction[7](122,111) = 1;
-                  this->restriction[7](123,75) = 1;
-                  break;
-            default:
-            {
-                                               // in case we don't
-                                               // have the matrices
-                                               // (yet), reset them to
-                                               // zero size. this does
-                                               // not prevent the use
-                                               // of this FE, but will
-                                               // prevent the use of
-                                               // these matrices
-              for (unsigned int i=0;
-                   i<GeometryInfo<dim>::children_per_cell;
-                   ++i)
-                this->restriction[i].reinit(0,0);
-            };
-          }
-        break;
-      };
-      
-      default:
-            Assert (false, ExcNotImplemented());
-    }
-
-                                  // finally fill in support points
-                                  // on cell and face
-  initialize_unit_support_points ();
-  initialize_unit_face_support_points ();
-}
-
-
-
-template <int dim>
-FiniteElement<dim> *
-FE_Q<dim>::clone() const
-{
-  return new FE_Q<dim>(degree);
-}
-
-
-
-template <int dim>
-double
-FE_Q<dim>::shape_value (const unsigned int i,
-                       const Point<dim> &p) const
-{
-  Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
-  return polynomial_space.compute_value(renumber_inverse[i], p);
-}
-
-
-template <int dim>
-double
-FE_Q<dim>::shape_value_component (const unsigned int i,
-                                 const Point<dim> &p,
-                                 const unsigned int component) const
-{
-  Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
-  Assert (component == 0, ExcIndexRange (component, 0, 1));
-  return polynomial_space.compute_value(renumber_inverse[i], p);
-}
-
-
-
-template <int dim>
-Tensor<1,dim>
-FE_Q<dim>::shape_grad (const unsigned int i,
-                      const Point<dim> &p) const
-{
-  Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
-  return polynomial_space.compute_grad(renumber_inverse[i], p);
-}
-
-
-
-template <int dim>
-Tensor<1,dim>
-FE_Q<dim>::shape_grad_component (const unsigned int i,
-                                const Point<dim> &p,
-                                const unsigned int component) const
-{
-  Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
-  Assert (component == 0, ExcIndexRange (component, 0, 1));
-  return polynomial_space.compute_grad(renumber_inverse[i], p);
-}
-
-
-
-template <int dim>
-Tensor<2,dim>
-FE_Q<dim>::shape_grad_grad (const unsigned int i,
-                           const Point<dim> &p) const
-{
-  Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
-  return polynomial_space.compute_grad_grad(renumber_inverse[i], p);
-}
-
-
-
-template <int dim>
-Tensor<2,dim>
-FE_Q<dim>::shape_grad_grad_component (const unsigned int i,
-                                     const Point<dim> &p,
-                                     const unsigned int component) const
-{
-  Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
-  Assert (component == 0, ExcIndexRange (component, 0, 1));
-  return polynomial_space.compute_grad_grad(renumber_inverse[i], p);
-}
-
-
-//----------------------------------------------------------------------
-// Auxiliary functions
-//----------------------------------------------------------------------
-
-
-
-template <int dim>
-void FE_Q<dim>::initialize_unit_support_points ()
-{
-                                  // number of points: (degree+1)^dim
-  unsigned int n = degree+1;
-  for (unsigned int i=1; i<dim; ++i)
-    n *= degree+1;
-  
-  this->unit_support_points.resize(n);
-  
-  const double step = 1./degree;
-  Point<dim> p;
-  
-  unsigned int k=0;
-  for (unsigned int iz=0; iz <= ((dim>2) ? degree : 0) ; ++iz)
-    for (unsigned int iy=0; iy <= ((dim>1) ? degree : 0) ; ++iy)
-      for (unsigned int ix=0; ix<=degree; ++ix)
-       {
-         p(0) = ix * step;
-         if (dim>1)
-           p(1) = iy * step;
-         if (dim>2)
-           p(2) = iz * step;
-         
-         this->unit_support_points[renumber[k++]] = p;
-       };
-}
-
-
-#if deal_II_dimension == 1
-
-template <>
-void FE_Q<1>::initialize_unit_face_support_points ()
-{
-                                  // no faces in 1d, so nothing to do
-}
-
-#endif
-
-
-template <int dim>
-void FE_Q<dim>::initialize_unit_face_support_points ()
-{
-  const unsigned int codim = dim-1;
-  
-                                  // number of points: (degree+1)^codim
-  unsigned int n = degree+1;
-  for (unsigned int i=1; i<codim; ++i)
-    n *= degree+1;
-  
-  this->unit_face_support_points.resize(n);
-  
-  const double step = 1./degree;
-  Point<codim> p;
-  
-  unsigned int k=0;
-  for (unsigned int iz=0; iz <= ((codim>2) ? degree : 0) ; ++iz)
-    for (unsigned int iy=0; iy <= ((codim>1) ? degree : 0) ; ++iy)
-      for (unsigned int ix=0; ix<=degree; ++ix)
-       {
-         p(0) = ix * step;
-         if (codim>1)
-           p(1) = iy * step;
-         if (codim>2)
-           p(2) = iz * step;
-         
-         this->unit_face_support_points[face_renumber[k++]] = p;
-       };
-}
-
-
-
-template <int dim>
-std::vector<unsigned int>
-FE_Q<dim>::get_dpo_vector(const unsigned int deg)
-{
-  std::vector<unsigned int> dpo(dim+1, static_cast<unsigned int>(1));
-  for (unsigned int i=1; i<dpo.size(); ++i)
-    dpo[i]=dpo[i-1]*(deg-1);
-  return dpo;
-}
-
-
-
-template <int dim>
-std::vector<unsigned int>
-FE_Q<dim>::lexicographic_to_hierarchic_numbering (const FiniteElementData<dim> &fe_data,
-                                                 const unsigned int            degree)
-{
-  std::vector<unsigned int> renumber (fe_data.dofs_per_cell);
-  
-  const unsigned int n = degree+1;
-
-
-  if (degree == 0)
-    {
-      Assert ((fe_data.dofs_per_vertex == 0) &&
-             ((fe_data.dofs_per_line == 0) || (dim == 1)) &&
-             ((fe_data.dofs_per_quad == 0) || (dim == 2)) &&
-             ((fe_data.dofs_per_hex == 0)  || (dim == 3)),
-             ExcInternalError());
-      renumber[0] = 0;
-    };
-
-  if (degree > 0)
+  if (degree > 0)
     {
       Assert (fe_data.dofs_per_vertex == 1, ExcInternalError());
       for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
@@ -1015,46 +449,644 @@ FE_Q<dim>::lexicographic_to_hierarchic_numbering (const FiniteElementData<dim> &
              }
        }
 
-      if (GeometryInfo<dim>::hexes_per_cell > 0)
-       for (int i=0; i<static_cast<signed int>(GeometryInfo<dim>::hexes_per_cell); ++i)
-         {
-           unsigned int index = fe_data.first_hex_index;
-           
-           for (unsigned int jz = 1; jz<degree; jz++)
-             for (unsigned int jy = 1; jy<degree; jy++)
-               for (unsigned int jx = 1; jx<degree; jx++)
-                 {
-                   const unsigned int tensorindex = jx + jy*n + jz*n*n;
-                   Assert (tensorindex<renumber.size(), ExcInternalError());
-                   renumber[tensorindex]=index++;
-                 }  
-         } 
+      if (GeometryInfo<dim>::hexes_per_cell > 0)
+       for (int i=0; i<static_cast<signed int>(GeometryInfo<dim>::hexes_per_cell); ++i)
+         {
+           unsigned int index = fe_data.first_hex_index;
+           
+           for (unsigned int jz = 1; jz<degree; jz++)
+             for (unsigned int jy = 1; jy<degree; jy++)
+               for (unsigned int jx = 1; jx<degree; jx++)
+                 {
+                   const unsigned int tensorindex = jx + jy*n + jz*n*n;
+                   Assert (tensorindex<renumber.size(), ExcInternalError());
+                   renumber[tensorindex]=index++;
+                 }  
+         } 
+    }
+
+  return renumber;
+}
+
+
+
+template <int dim>
+std::vector<unsigned int>
+FE_Q<dim>::face_lexicographic_to_hierarchic_numbering (const unsigned int degree)
+{
+  const FiniteElementData<dim-1> fe_data(FE_Q<dim-1>::get_dpo_vector(degree),1);
+  return FE_Q<dim-1>::lexicographic_to_hierarchic_numbering (fe_data, degree); 
+}
+
+
+#if (deal_II_dimension == 1)
+
+template <>
+std::vector<unsigned int>
+FE_Q<1>::face_lexicographic_to_hierarchic_numbering (const unsigned int)
+{
+  return std::vector<unsigned int>();
+}
+
+#endif
+
+
+
+template <int dim>
+void
+FE_Q<dim>::initialize_constraints ()
+{  
+                                  // copy constraint matrices if they
+                                  // are defined. otherwise leave them
+                                  // at invalid size
+  if ((dim > 1) && (degree < Matrices::n_constraint_matrices+1))
+    {
+      this->interface_constraints.
+        TableBase<2,double>::reinit (this->interface_constraints_size());
+      this->interface_constraints.fill (Matrices::constraint_matrices[degree-1]);
+    };
+}
+
+
+
+template <int dim>
+void
+FE_Q<dim>::initialize_embedding ()
+{
+                                  // copy over embedding matrices if
+                                  // they are defined
+  if ((degree < Matrices::n_embedding_matrices+1) &&
+      (Matrices::embedding[degree-1][0] != 0))
+    for (unsigned int c=0; c<GeometryInfo<dim>::children_per_cell; ++c)
+      {
+        this->prolongation[c].reinit (this->dofs_per_cell,
+                                      this->dofs_per_cell);
+        this->prolongation[c].fill (Matrices::embedding[degree-1][c]);
+      };
+}
+
+
+
+template <int dim>
+void
+FE_Q<dim>::initialize_restriction ()
+{
+  
+                                  // then fill restriction
+                                  // matrices. they are hardcoded for
+                                  // the first few elements. in
+                                  // contrast to the other matrices,
+                                  // these are not stored in the
+                                  // files fe_q_[123]d.cc, since they
+                                  // contain only a rather small
+                                  // number of zeros, and storing
+                                  // them element-wise is more
+                                  // expensive than just setting the
+                                  // nonzero elements as done here
+  for (unsigned int c=0; c<GeometryInfo<dim>::children_per_cell; ++c)
+    this->restriction[c].reinit (this->dofs_per_cell, this->dofs_per_cell);
+  switch (dim)
+    {
+      case 1:  // 1d
+      {
+        switch (degree)
+          {
+            case 1:
+                  this->restriction[0](0,0) = 1;
+                  this->restriction[1](1,1) = 1;
+                  break;
+            case 2:
+                  this->restriction[0](0,0) = 1;
+                  this->restriction[0](2,1) = 1;
+                  this->restriction[1](1,1) = 1;
+                  this->restriction[1](1,1) = 1;
+                  break;
+            case 3:
+                  this->restriction[0](0,0) = 1;
+                  this->restriction[0](2,3) = 1;
+                  this->restriction[1](1,1) = 1;
+                  this->restriction[1](3,2) = 1;
+                  break;
+            case 4:
+                  this->restriction[0](0,0) = 1;
+                  this->restriction[0](2,3) = 1;
+                  this->restriction[0](3,1) = 1;
+                  this->restriction[1](1,1) = 1;
+                  this->restriction[1](3,0) = 1;
+                  this->restriction[1](4,3) = 1;
+                  break;
+                 
+            default:
+            {
+                                               // in case we don't
+                                               // have the matrices
+                                               // (yet), reset them to
+                                               // zero size. this does
+                                               // not prevent the use
+                                               // of this FE, but will
+                                               // prevent the use of
+                                               // these matrices
+              for (unsigned int i=0;
+                   i<GeometryInfo<dim>::children_per_cell;
+                   ++i)
+                this->restriction[i].reinit(0,0);
+            };
+          }
+        break;
+      };
+      
+      case 2:  // 2d
+      {
+        switch (degree)
+          {
+            case 1:
+                  this->restriction[0](0,0) = 1;
+                  this->restriction[1](1,1) = 1;
+                  this->restriction[2](2,2) = 1;
+                  this->restriction[3](3,3) = 1;
+                  break;
+            case 2:
+                  this->restriction[0](0,0) = 1;
+                  this->restriction[0](4,1) = 1;
+                  this->restriction[0](7,3) = 1;
+                  this->restriction[0](8,2) = 1;
+                  this->restriction[1](1,1) = 1;
+                  this->restriction[1](4,0) = 1;
+                  this->restriction[1](5,2) = 1;
+                  this->restriction[1](8,3) = 1;
+                  this->restriction[2](2,2) = 1;
+                  this->restriction[2](5,1) = 1;
+                  this->restriction[2](6,3) = 1;
+                  this->restriction[2](8,0) = 1;
+                  this->restriction[3](3,3) = 1;
+                  this->restriction[3](6,2) = 1;
+                  this->restriction[3](7,0) = 1;
+                  this->restriction[3](8,1) = 1;
+                  break;
+            case 3:
+                  this->restriction[0](0,0) = 1;
+                  this->restriction[0](4,5) = 1;
+                  this->restriction[0](10,11) = 1;
+                  this->restriction[0](12,15) = 1;
+                  this->restriction[1](1,1) = 1;
+                  this->restriction[1](5,4) = 1;
+                  this->restriction[1](6,7) = 1;
+                  this->restriction[1](13,14) = 1;
+                  this->restriction[2](2,2) = 1;
+                  this->restriction[2](7,6) = 1;
+                  this->restriction[2](9,8) = 1;
+                  this->restriction[2](15,12) = 1;
+                  this->restriction[3](3,3) = 1;
+                  this->restriction[3](8,9) = 1;
+                  this->restriction[3](11,10) = 1;
+                  this->restriction[3](14,13) = 1;
+                  break;
+            case 4:
+                  this->restriction[0](0,0) = 1;
+                  this->restriction[0](4,5) = 1;
+                  this->restriction[0](5,1) = 1;
+                  this->restriction[0](13,14) = 1;
+                  this->restriction[0](14,3) = 1;
+                  this->restriction[0](16,20) = 1;
+                  this->restriction[0](17,8) = 1;
+                  this->restriction[0](19,11) = 1;
+                  this->restriction[0](20,2) = 1;
+                  this->restriction[1](1,1) = 1;
+                  this->restriction[1](5,0) = 1;
+                  this->restriction[1](6,5) = 1;
+                  this->restriction[1](7,8) = 1;
+                  this->restriction[1](8,2) = 1;
+                  this->restriction[1](17,14) = 1;
+                  this->restriction[1](18,20) = 1;
+                  this->restriction[1](20,3) = 1;
+                  this->restriction[1](21,11) = 1;
+                  this->restriction[2](2,2) = 1;
+                  this->restriction[2](8,1) = 1;
+                  this->restriction[2](9,8) = 1;
+                  this->restriction[2](11,3) = 1;
+                  this->restriction[2](12,11) = 1;
+                  this->restriction[2](20,0) = 1;
+                  this->restriction[2](21,5) = 1;
+                  this->restriction[2](23,14) = 1;
+                  this->restriction[2](24,20) = 1;
+                  this->restriction[3](3,3) = 1;
+                  this->restriction[3](10,11) = 1;
+                  this->restriction[3](11,2) = 1;
+                  this->restriction[3](14,0) = 1;
+                  this->restriction[3](15,14) = 1;
+                  this->restriction[3](19,5) = 1;
+                  this->restriction[3](20,1) = 1;
+                  this->restriction[3](22,20) = 1;
+                  this->restriction[3](23,8) = 1;
+                  break;
+
+            default:
+            {
+                                               // in case we don't
+                                               // have the matrices
+                                               // (yet), reset them to
+                                               // zero size. this does
+                                               // not prevent the use
+                                               // of this FE, but will
+                                               // prevent the use of
+                                               // these matrices
+              for (unsigned int i=0;
+                   i<GeometryInfo<dim>::children_per_cell;
+                   ++i)
+                this->restriction[i].reinit(0,0);
+            };
+          }
+        break;
+      };
+      
+      case 3:  // 3d
+      {
+        switch (degree)
+          {
+            case 1:
+                  this->restriction[0](0,0) = 1;
+                  this->restriction[1](1,1) = 1;
+                  this->restriction[2](2,2) = 1;
+                  this->restriction[3](3,3) = 1;
+                  this->restriction[4](4,4) = 1;
+                  this->restriction[5](5,5) = 1;
+                  this->restriction[6](6,6) = 1;
+                  this->restriction[7](7,7) = 1;
+                  break;
+            case 2:
+                  this->restriction[0](0,0) = 1;
+                  this->restriction[0](8,1) = 1;
+                  this->restriction[0](11,3) = 1;
+                  this->restriction[0](16,4) = 1;
+                  this->restriction[0](20,2) = 1;
+                  this->restriction[0](22,5) = 1;
+                  this->restriction[0](25,7) = 1;
+                  this->restriction[0](26,6) = 1;
+                  this->restriction[1](1,1) = 1;
+                  this->restriction[1](8,0) = 1;
+                  this->restriction[1](9,2) = 1;
+                  this->restriction[1](17,5) = 1;
+                  this->restriction[1](20,3) = 1;
+                  this->restriction[1](22,4) = 1;
+                  this->restriction[1](23,6) = 1;
+                  this->restriction[1](26,7) = 1;
+                  this->restriction[2](2,2) = 1;
+                  this->restriction[2](9,1) = 1;
+                  this->restriction[2](10,3) = 1;
+                  this->restriction[2](18,6) = 1;
+                  this->restriction[2](20,0) = 1;
+                  this->restriction[2](23,5) = 1;
+                  this->restriction[2](24,7) = 1;
+                  this->restriction[2](26,4) = 1;
+                  this->restriction[3](3,3) = 1;
+                  this->restriction[3](10,2) = 1;
+                  this->restriction[3](11,0) = 1;
+                  this->restriction[3](19,7) = 1;
+                  this->restriction[3](20,1) = 1;
+                  this->restriction[3](24,6) = 1;
+                  this->restriction[3](25,4) = 1;
+                  this->restriction[3](26,5) = 1;
+                  this->restriction[4](4,4) = 1;
+                  this->restriction[4](12,5) = 1;
+                  this->restriction[4](15,7) = 1;
+                  this->restriction[4](16,0) = 1;
+                  this->restriction[4](21,6) = 1;
+                  this->restriction[4](22,1) = 1;
+                  this->restriction[4](25,3) = 1;
+                  this->restriction[4](26,2) = 1;
+                  this->restriction[5](5,5) = 1;
+                  this->restriction[5](12,4) = 1;
+                  this->restriction[5](13,6) = 1;
+                  this->restriction[5](17,1) = 1;
+                  this->restriction[5](21,7) = 1;
+                  this->restriction[5](22,0) = 1;
+                  this->restriction[5](23,2) = 1;
+                  this->restriction[5](26,3) = 1;
+                  this->restriction[6](6,6) = 1;
+                  this->restriction[6](13,5) = 1;
+                  this->restriction[6](14,7) = 1;
+                  this->restriction[6](18,2) = 1;
+                  this->restriction[6](21,4) = 1;
+                  this->restriction[6](23,1) = 1;
+                  this->restriction[6](24,3) = 1;
+                  this->restriction[6](26,0) = 1;
+                  this->restriction[7](7,7) = 1;
+                  this->restriction[7](14,6) = 1;
+                  this->restriction[7](15,4) = 1;
+                  this->restriction[7](19,3) = 1;
+                  this->restriction[7](21,5) = 1;
+                  this->restriction[7](24,2) = 1;
+                  this->restriction[7](25,0) = 1;
+                  this->restriction[7](26,1) = 1;
+                  break;
+            case 3:
+                  this->restriction[0](0,0) = 1;
+                  this->restriction[0](8,9) = 1;
+                  this->restriction[0](14,15) = 1;
+                  this->restriction[0](24,25) = 1;
+                  this->restriction[0](32,35) = 1;
+                  this->restriction[0](40,43) = 1;
+                  this->restriction[0](52,55) = 1;
+                  this->restriction[0](56,63) = 1;
+                  this->restriction[1](1,1) = 1;
+                  this->restriction[1](9,8) = 1;
+                  this->restriction[1](10,11) = 1;
+                  this->restriction[1](26,27) = 1;
+                  this->restriction[1](33,34) = 1;
+                  this->restriction[1](41,42) = 1;
+                  this->restriction[1](44,47) = 1;
+                  this->restriction[1](57,62) = 1;
+                  this->restriction[2](2,2) = 1;
+                  this->restriction[2](11,10) = 1;
+                  this->restriction[2](13,12) = 1;
+                  this->restriction[2](28,29) = 1;
+                  this->restriction[2](35,32) = 1;
+                  this->restriction[2](46,45) = 1;
+                  this->restriction[2](49,50) = 1;
+                  this->restriction[2](61,58) = 1;
+                  this->restriction[3](3,3) = 1;
+                  this->restriction[3](12,13) = 1;
+                  this->restriction[3](15,14) = 1;
+                  this->restriction[3](30,31) = 1;
+                  this->restriction[3](34,33) = 1;
+                  this->restriction[3](48,51) = 1;
+                  this->restriction[3](54,53) = 1;
+                  this->restriction[3](60,59) = 1;
+                  this->restriction[4](4,4) = 1;
+                  this->restriction[4](16,17) = 1;
+                  this->restriction[4](22,23) = 1;
+                  this->restriction[4](25,24) = 1;
+                  this->restriction[4](36,39) = 1;
+                  this->restriction[4](42,41) = 1;
+                  this->restriction[4](53,54) = 1;
+                  this->restriction[4](58,61) = 1;
+                  this->restriction[5](5,5) = 1;
+                  this->restriction[5](17,16) = 1;
+                  this->restriction[5](18,19) = 1;
+                  this->restriction[5](27,26) = 1;
+                  this->restriction[5](37,38) = 1;
+                  this->restriction[5](43,40) = 1;
+                  this->restriction[5](45,46) = 1;
+                  this->restriction[5](59,60) = 1;
+                  this->restriction[6](6,6) = 1;
+                  this->restriction[6](19,18) = 1;
+                  this->restriction[6](21,20) = 1;
+                  this->restriction[6](29,28) = 1;
+                  this->restriction[6](39,36) = 1;
+                  this->restriction[6](47,44) = 1;
+                  this->restriction[6](51,48) = 1;
+                  this->restriction[6](63,56) = 1;
+                  this->restriction[7](7,7) = 1;
+                  this->restriction[7](20,21) = 1;
+                  this->restriction[7](23,22) = 1;
+                  this->restriction[7](31,30) = 1;
+                  this->restriction[7](38,37) = 1;
+                  this->restriction[7](50,49) = 1;
+                  this->restriction[7](55,52) = 1;
+                  this->restriction[7](62,57) = 1;
+                  break;
+            case 4:
+                  this->restriction[0](0,0) = 1;
+                  this->restriction[0](8,9) = 1;
+                  this->restriction[0](9,1) = 1;
+                  this->restriction[0](17,18) = 1;
+                  this->restriction[0](18,3) = 1;
+                  this->restriction[0](32,33) = 1;
+                  this->restriction[0](33,4) = 1;
+                  this->restriction[0](44,48) = 1;
+                  this->restriction[0](45,12) = 1;
+                  this->restriction[0](47,15) = 1;
+                  this->restriction[0](48,2) = 1;
+                  this->restriction[0](62,66) = 1;
+                  this->restriction[0](63,36) = 1;
+                  this->restriction[0](65,21) = 1;
+                  this->restriction[0](66,5) = 1;
+                  this->restriction[0](89,93) = 1;
+                  this->restriction[0](90,30) = 1;
+                  this->restriction[0](92,42) = 1;
+                  this->restriction[0](93,7) = 1;
+                  this->restriction[0](98,111) = 1;
+                  this->restriction[0](99,75) = 1;
+                  this->restriction[0](101,57) = 1;
+                  this->restriction[0](102,24) = 1;
+                  this->restriction[0](107,84) = 1;
+                  this->restriction[0](108,39) = 1;
+                  this->restriction[0](110,27) = 1;
+                  this->restriction[0](111,6) = 1;
+                  this->restriction[1](1,1) = 1;
+                  this->restriction[1](9,0) = 1;
+                  this->restriction[1](10,9) = 1;
+                  this->restriction[1](11,12) = 1;
+                  this->restriction[1](12,2) = 1;
+                  this->restriction[1](35,36) = 1;
+                  this->restriction[1](36,5) = 1;
+                  this->restriction[1](45,18) = 1;
+                  this->restriction[1](46,48) = 1;
+                  this->restriction[1](48,3) = 1;
+                  this->restriction[1](49,15) = 1;
+                  this->restriction[1](63,33) = 1;
+                  this->restriction[1](64,66) = 1;
+                  this->restriction[1](66,4) = 1;
+                  this->restriction[1](67,21) = 1;
+                  this->restriction[1](71,75) = 1;
+                  this->restriction[1](72,24) = 1;
+                  this->restriction[1](74,39) = 1;
+                  this->restriction[1](75,6) = 1;
+                  this->restriction[1](99,93) = 1;
+                  this->restriction[1](100,111) = 1;
+                  this->restriction[1](102,30) = 1;
+                  this->restriction[1](103,57) = 1;
+                  this->restriction[1](108,42) = 1;
+                  this->restriction[1](109,84) = 1;
+                  this->restriction[1](111,7) = 1;
+                  this->restriction[1](112,27) = 1;
+                  this->restriction[2](2,2) = 1;
+                  this->restriction[2](12,1) = 1;
+                  this->restriction[2](13,12) = 1;
+                  this->restriction[2](15,3) = 1;
+                  this->restriction[2](16,15) = 1;
+                  this->restriction[2](38,39) = 1;
+                  this->restriction[2](39,6) = 1;
+                  this->restriction[2](48,0) = 1;
+                  this->restriction[2](49,9) = 1;
+                  this->restriction[2](51,18) = 1;
+                  this->restriction[2](52,48) = 1;
+                  this->restriction[2](74,36) = 1;
+                  this->restriction[2](75,5) = 1;
+                  this->restriction[2](77,75) = 1;
+                  this->restriction[2](78,24) = 1;
+                  this->restriction[2](81,42) = 1;
+                  this->restriction[2](82,84) = 1;
+                  this->restriction[2](84,7) = 1;
+                  this->restriction[2](85,27) = 1;
+                  this->restriction[2](108,33) = 1;
+                  this->restriction[2](109,66) = 1;
+                  this->restriction[2](111,4) = 1;
+                  this->restriction[2](112,21) = 1;
+                  this->restriction[2](117,93) = 1;
+                  this->restriction[2](118,111) = 1;
+                  this->restriction[2](120,30) = 1;
+                  this->restriction[2](121,57) = 1;
+                  this->restriction[3](3,3) = 1;
+                  this->restriction[3](14,15) = 1;
+                  this->restriction[3](15,2) = 1;
+                  this->restriction[3](18,0) = 1;
+                  this->restriction[3](19,18) = 1;
+                  this->restriction[3](41,42) = 1;
+                  this->restriction[3](42,7) = 1;
+                  this->restriction[3](47,9) = 1;
+                  this->restriction[3](48,1) = 1;
+                  this->restriction[3](50,48) = 1;
+                  this->restriction[3](51,12) = 1;
+                  this->restriction[3](80,84) = 1;
+                  this->restriction[3](81,39) = 1;
+                  this->restriction[3](83,27) = 1;
+                  this->restriction[3](84,6) = 1;
+                  this->restriction[3](92,33) = 1;
+                  this->restriction[3](93,4) = 1;
+                  this->restriction[3](95,93) = 1;
+                  this->restriction[3](96,30) = 1;
+                  this->restriction[3](107,66) = 1;
+                  this->restriction[3](108,36) = 1;
+                  this->restriction[3](110,21) = 1;
+                  this->restriction[3](111,5) = 1;
+                  this->restriction[3](116,111) = 1;
+                  this->restriction[3](117,75) = 1;
+                  this->restriction[3](119,57) = 1;
+                  this->restriction[3](120,24) = 1;
+                  this->restriction[4](4,4) = 1;
+                  this->restriction[4](20,21) = 1;
+                  this->restriction[4](21,5) = 1;
+                  this->restriction[4](29,30) = 1;
+                  this->restriction[4](30,7) = 1;
+                  this->restriction[4](33,0) = 1;
+                  this->restriction[4](34,33) = 1;
+                  this->restriction[4](53,57) = 1;
+                  this->restriction[4](54,24) = 1;
+                  this->restriction[4](56,27) = 1;
+                  this->restriction[4](57,6) = 1;
+                  this->restriction[4](65,9) = 1;
+                  this->restriction[4](66,1) = 1;
+                  this->restriction[4](68,66) = 1;
+                  this->restriction[4](69,36) = 1;
+                  this->restriction[4](90,18) = 1;
+                  this->restriction[4](91,93) = 1;
+                  this->restriction[4](93,3) = 1;
+                  this->restriction[4](94,42) = 1;
+                  this->restriction[4](101,48) = 1;
+                  this->restriction[4](102,12) = 1;
+                  this->restriction[4](104,111) = 1;
+                  this->restriction[4](105,75) = 1;
+                  this->restriction[4](110,15) = 1;
+                  this->restriction[4](111,2) = 1;
+                  this->restriction[4](113,84) = 1;
+                  this->restriction[4](114,39) = 1;
+                  this->restriction[5](5,5) = 1;
+                  this->restriction[5](21,4) = 1;
+                  this->restriction[5](22,21) = 1;
+                  this->restriction[5](23,24) = 1;
+                  this->restriction[5](24,6) = 1;
+                  this->restriction[5](36,1) = 1;
+                  this->restriction[5](37,36) = 1;
+                  this->restriction[5](54,30) = 1;
+                  this->restriction[5](55,57) = 1;
+                  this->restriction[5](57,7) = 1;
+                  this->restriction[5](58,27) = 1;
+                  this->restriction[5](66,0) = 1;
+                  this->restriction[5](67,9) = 1;
+                  this->restriction[5](69,33) = 1;
+                  this->restriction[5](70,66) = 1;
+                  this->restriction[5](72,12) = 1;
+                  this->restriction[5](73,75) = 1;
+                  this->restriction[5](75,2) = 1;
+                  this->restriction[5](76,39) = 1;
+                  this->restriction[5](102,18) = 1;
+                  this->restriction[5](103,48) = 1;
+                  this->restriction[5](105,93) = 1;
+                  this->restriction[5](106,111) = 1;
+                  this->restriction[5](111,3) = 1;
+                  this->restriction[5](112,15) = 1;
+                  this->restriction[5](114,42) = 1;
+                  this->restriction[5](115,84) = 1;
+                  this->restriction[6](6,6) = 1;
+                  this->restriction[6](24,5) = 1;
+                  this->restriction[6](25,24) = 1;
+                  this->restriction[6](27,7) = 1;
+                  this->restriction[6](28,27) = 1;
+                  this->restriction[6](39,2) = 1;
+                  this->restriction[6](40,39) = 1;
+                  this->restriction[6](57,4) = 1;
+                  this->restriction[6](58,21) = 1;
+                  this->restriction[6](60,30) = 1;
+                  this->restriction[6](61,57) = 1;
+                  this->restriction[6](75,1) = 1;
+                  this->restriction[6](76,36) = 1;
+                  this->restriction[6](78,12) = 1;
+                  this->restriction[6](79,75) = 1;
+                  this->restriction[6](84,3) = 1;
+                  this->restriction[6](85,15) = 1;
+                  this->restriction[6](87,42) = 1;
+                  this->restriction[6](88,84) = 1;
+                  this->restriction[6](111,0) = 1;
+                  this->restriction[6](112,9) = 1;
+                  this->restriction[6](114,33) = 1;
+                  this->restriction[6](115,66) = 1;
+                  this->restriction[6](120,18) = 1;
+                  this->restriction[6](121,48) = 1;
+                  this->restriction[6](123,93) = 1;
+                  this->restriction[6](124,111) = 1;
+                  this->restriction[7](7,7) = 1;
+                  this->restriction[7](26,27) = 1;
+                  this->restriction[7](27,6) = 1;
+                  this->restriction[7](30,4) = 1;
+                  this->restriction[7](31,30) = 1;
+                  this->restriction[7](42,3) = 1;
+                  this->restriction[7](43,42) = 1;
+                  this->restriction[7](56,21) = 1;
+                  this->restriction[7](57,5) = 1;
+                  this->restriction[7](59,57) = 1;
+                  this->restriction[7](60,24) = 1;
+                  this->restriction[7](83,15) = 1;
+                  this->restriction[7](84,2) = 1;
+                  this->restriction[7](86,84) = 1;
+                  this->restriction[7](87,39) = 1;
+                  this->restriction[7](93,0) = 1;
+                  this->restriction[7](94,33) = 1;
+                  this->restriction[7](96,18) = 1;
+                  this->restriction[7](97,93) = 1;
+                  this->restriction[7](110,9) = 1;
+                  this->restriction[7](111,1) = 1;
+                  this->restriction[7](113,66) = 1;
+                  this->restriction[7](114,36) = 1;
+                  this->restriction[7](119,48) = 1;
+                  this->restriction[7](120,12) = 1;
+                  this->restriction[7](122,111) = 1;
+                  this->restriction[7](123,75) = 1;
+                  break;
+            default:
+            {
+                                               // in case we don't
+                                               // have the matrices
+                                               // (yet), reset them to
+                                               // zero size. this does
+                                               // not prevent the use
+                                               // of this FE, but will
+                                               // prevent the use of
+                                               // these matrices
+              for (unsigned int i=0;
+                   i<GeometryInfo<dim>::children_per_cell;
+                   ++i)
+                this->restriction[i].reinit(0,0);
+            };
+          }
+        break;
+      };
+      
+      default:
+            Assert (false, ExcNotImplemented());
     }
-
-  return renumber;
-}
-
-
-
-template <int dim>
-std::vector<unsigned int>
-FE_Q<dim>::face_lexicographic_to_hierarchic_numbering (const unsigned int degree)
-{
-  const FiniteElementData<dim-1> fe_data(FE_Q<dim-1>::get_dpo_vector(degree),1);
-  return FE_Q<dim-1>::lexicographic_to_hierarchic_numbering (fe_data, degree); 
-}
-
-
-#if (deal_II_dimension == 1)
-
-template <>
-std::vector<unsigned int>
-FE_Q<1>::face_lexicographic_to_hierarchic_numbering (const unsigned int)
-{
-  return std::vector<unsigned int>();
 }
 
-#endif
 
 
 template <int dim>

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.