]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Unify the way the various finite elements are initialized in their constructors,...
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 23 Apr 2003 23:56:24 +0000 (23:56 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 23 Apr 2003 23:56:24 +0000 (23:56 +0000)
git-svn-id: https://svn.dealii.org/trunk@7464 0785d39b-7218-0410-832d-ea1e28bc413d

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

index a7553e53be1d1a2b2945d5e9d0817f58ee62c6c9..e1d1472ca87096d1578e4e0b6fd91df069b4b688 100644 (file)
@@ -447,6 +447,28 @@ class FE_Nedelec : public FiniteElement<dim>
                                      */
     static std::vector<unsigned int> get_dpo_vector(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 54d59b934bc5ab997975c4d341d8e58b2ab014e4..816bb73f6878de5291cd538bdeccd2a7caa4bc78 100644 (file)
@@ -560,10 +560,9 @@ class FE_Q_Hierarchical : public FiniteElement<dim>
                                      * proper design.
                                      */
     static
-    void
+    std::vector<unsigned int>
     lexicographic_to_hierarchic_numbering (const FiniteElementData<dim> &fe_data,
-                                          const unsigned int            degree,
-                                          std::vector<unsigned int>    &numbering);
+                                          const unsigned int            degree);
 
                                     /**
                                      * This is an analogon to the
@@ -571,10 +570,32 @@ class FE_Q_Hierarchical : public FiniteElement<dim>
                                      * on faces.
                                      */
     static
-    void
-    face_lexicographic_to_hierarchic_numbering (const unsigned int         degree,
-                                               std::vector<unsigned int> &numbering);
+    std::vector<unsigned int>
+    face_lexicographic_to_hierarchic_numbering (const unsigned int degree);
 
+                                    /**
+                                     * Initialize two auxiliary
+                                     * fields that will be used in
+                                     * setting up the various
+                                     * matrices in the constructor.
+                                     */
+    void build_dofs_cell (std::vector<FullMatrix<double> > &dofs_cell,
+                         std::vector<FullMatrix<double> > &dofs_subcell) const;
+    
+                                    /**
+                                     * Initialize the hanging node
+                                     * constraints matrices. Called
+                                     * from the constructor.
+                                     */
+    void initialize_constraints (const std::vector<FullMatrix<double> > &dofs_subcell);
+
+                                    /**
+                                     * Initialize the embedding
+                                     * matrices. Called from the
+                                     * constructor.
+                                     */
+    void initialize_embedding_and_restriction (const std::vector<FullMatrix<double> > &dofs_cell,
+                                              const std::vector<FullMatrix<double> > &dofs_subcell);
 
                                     /**
                                      * Initialize the
@@ -672,7 +693,7 @@ class FE_Q_Hierarchical : public FiniteElement<dim>
                                      * Mapping from lexicographic to
                                      * shape function numbering.
                                      */
-    std::vector<unsigned int> renumber;
+    const std::vector<unsigned int> renumber;
 
                                     /**
                                      * Inverse renumber
@@ -680,32 +701,13 @@ class FE_Q_Hierarchical : public FiniteElement<dim>
                                      * shape function numbering to
                                      * lexicographic numbering.
                                      */
-    std::vector<unsigned int> renumber_inverse;
+    const std::vector<unsigned int> renumber_inverse;
              
                                     /**
                                      * Mapping from lexicographic to
                                      * shape function numbering on first face.
                                      */
-    std::vector<unsigned int> face_renumber;
-
-                                     /**
-                                      * The matrix @p{dofs_cell} contains the 
-                                     * values of the linear functionals of 
-                                      * the master 1d cell applied to the 
-                                     * shape functions of the two 1d subcells.
-                                     * The matrix @p{dofs_subcell} constains
-                                     * the values of the linear functionals 
-                                     * on each 1d subcell applied to the 
-                                     * shape functions on the master 1d 
-                                     * subcell. 
-                                     * We use @p{dofs_cell} and 
-                                     * @p{dofs_subcell} to compute the 
-                                     * @p{prolongation}, @p{restriction} and 
-                                     * @p{interface_constraints} matrices 
-                                     * for all dimensions.
-                                     */
-    std::vector<FullMatrix<double> > dofs_cell;
-    std::vector<FullMatrix<double> > dofs_subcell;
+    const std::vector<unsigned int> face_renumber;
 
                                     /**
                                      * Pointer to the tensor
@@ -782,8 +784,11 @@ class FE_Q_Hierarchical : public FiniteElement<dim>
 
 /* -------------- declaration of explicit specializations ------------- */
 
-template <> void FE_Q_Hierarchical<1>::initialize_unit_face_support_points ();
-template <> void FE_Q_Hierarchical<1>::face_lexicographic_to_hierarchic_numbering (const unsigned int,
-                                                                                  std::vector<unsigned int>&);
+template <>
+void FE_Q_Hierarchical<1>::initialize_unit_face_support_points ();
+
+template <>
+std::vector<unsigned int>
+FE_Q_Hierarchical<1>::face_lexicographic_to_hierarchic_numbering (const unsigned int);
 
 #endif
index f93bc4e22fc3d3c73696f6283482f8fc4c726bd2..37bf020088880217d178493692d920d5a2fae7db 100644 (file)
@@ -34,219 +34,13 @@ FE_Nedelec<dim>::FE_Nedelec (const unsigned int degree)
 {
   Assert (dim >= 2, ExcNotUsefulInThisDimension());
   
-                                  // copy constraint matrices if they
-                                  // are defined. otherwise leave them
-                                  // at zero size
-  if (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]);
-    };
-
-                                  // 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)
-      {
-                                         // copy
-        this->prolongation[c].reinit (this->dofs_per_cell,
-                                      this->dofs_per_cell);
-        this->prolongation[c].fill (Matrices::embedding[degree-1][c]);
-                                         // and make sure that the row
-                                         // sum is 0.5 (for usual
-                                         // elements, the row sum must
-                                         // be 1, but here the shape
-                                         // function is multiplied by
-                                         // the inverse of the
-                                         // Jacobian, which introduces
-                                         // a factor of 1/2 when going
-                                         // from mother to child)
-        for (unsigned int row=0; row<this->dofs_per_cell; ++row)
-          {
-            double sum = 0;
-            for (unsigned int col=0; col<this->dofs_per_cell; ++col)
-              sum += this->prolongation[c](row,col);
-            Assert (std::fabs(sum-.5) < 1e-14,
-                    ExcInternalError());
-          };
-      };
-
-                                  // then fill restriction
-                                  // matrices. they are hardcoded for
-                                  // the first few elements
-  switch (dim)
-    {
-      case 2:   // 2d
-      {
-       switch (degree)
-         {
-           case 1:
-           {
-                                               // this is a strange
-                                               // element, since it is
-                                               // both additive and
-                                               // then it is also
-                                               // not. ideally, we
-                                               // would like to have
-                                               // the value of the
-                                               // shape function on
-                                               // the coarse line to
-                                               // be the mean value of
-                                               // that on the two
-                                               // child ones. thus,
-                                               // one should make it
-                                               // additive. however,
-                                               // additivity only
-                                               // works if an element
-                                               // does not have any
-                                               // continuity
-                                               // requirements, since
-                                               // otherwise degrees of
-                                               // freedom are shared
-                                               // between adjacent
-                                               // elements, and when
-                                               // we make the element
-                                               // additive, that would
-                                               // mean that we end up
-                                               // adding up
-                                               // contributions not
-                                               // only from the child
-                                               // cells of this cell,
-                                               // but also from the
-                                               // child cells of the
-                                               // neighbor, and since
-                                               // we cannot know
-                                               // whether there even
-                                               // exists a neighbor we
-                                               // cannot simply make
-                                               // the element
-                                               // additive.
-                                              //
-                                               // so, until someone
-                                               // comes along with a
-                                               // better alternative,
-                                               // we do the following:
-                                               // make the element
-                                               // non-additive, and
-                                               // simply pick the
-                                               // value of one of the
-                                               // child lines for the
-                                               // value of the mother
-                                               // line (note that we
-                                               // have to multiply by
-                                               // two, since the shape
-                                               // functions scale with
-                                               // the inverse
-                                               // Jacobian). we thus
-                                               // throw away the
-                                               // information of one
-                                               // of the child lines,
-                                               // but there seems to
-                                               // be no other way than
-                                               // that...
-                                               //
-                                               // note: to make things
-                                               // consistent, and
-                                               // restriction
-                                               // independent of the
-                                               // order in which we
-                                               // travel across the
-                                               // cells of the coarse
-                                               // grid, we have to
-                                               // make sure that we
-                                               // take the same small
-                                               // line when visiting
-                                               // its two neighbors,
-                                               // to get the value for
-                                               // the mother line. we
-                                               // take the first line
-                                               // always, in the
-                                               // canonical direction
-                                               // of lines
-              for (unsigned int c=0; c<GeometryInfo<dim>::children_per_cell; ++c)
-                this->restriction[c].reinit (this->dofs_per_cell,
-                                             this->dofs_per_cell);
-              
-             this->restriction[0](0,0) = 2.;
-             this->restriction[1](1,1) = 2.;
-             this->restriction[3](2,2) = 2.;
-             this->restriction[0](3,3) = 2.;
-
-             break;
-           };
-           
-           default:
-           {
-                                              // in case we don't
-                                              // have the matrices
-                                              // (yet), leave them
-                                              // empty. this does not
-                                              // prevent the use of
-                                              // this FE, but will
-                                              // prevent the use of
-                                              // these matrices
-              break;
-           };
-         };
-       
-       break;
-      };
-
-
-      case 3:   // 3d
-      {
-       switch (degree)
-         {
-           case 1:
-           {
-                                              // same principle as in
-                                              // 2d, take one child
-                                              // cell to get at the
-                                              // values of each of
-                                              // the 12 lines
-              for (unsigned int c=0; c<GeometryInfo<dim>::children_per_cell; ++c)
-                this->restriction[c].reinit (this->dofs_per_cell,
-                                             this->dofs_per_cell);
-             this->restriction[0](0,0) = 2.;
-             this->restriction[0](3,3) = 2.;
-             this->restriction[1](1,1) = 2.;
-             this->restriction[3](2,2) = 2.;
-              
-             this->restriction[4](4,4) = 2.;
-             this->restriction[4](7,7) = 2.;
-             this->restriction[5](5,5) = 2.;
-             this->restriction[7](6,6) = 2.;
-              
-             this->restriction[0](8,8) = 2.;
-             this->restriction[1](9,9) = 2.;
-             this->restriction[2](10,10) = 2.;
-             this->restriction[3](11,11) = 2.;
-              
-             break;
-           };
-           
-           default:
-           {
-                                              // in case we don't
-                                              // have the matrices
-                                              // (yet), leave them
-                                              // empty. this does not
-                                              // prevent the use of
-                                              // this FE, but will
-                                              // prevent the use of
-                                              // these matrices
-              break;
-           };
-         };
-       
-       break;
-      };
-      
-      default:
-           Assert (false,ExcNotImplemented());
-    }
+                                  // 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
@@ -708,6 +502,236 @@ FE_Nedelec<3>::shape_grad_grad_component (const unsigned int i,
 
 
 
+template <int dim>
+void
+FE_Nedelec<dim>::initialize_constraints ()
+{
+                                  // copy constraint matrices if they
+                                  // are defined. otherwise leave
+                                  // them at zero size
+  if (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_Nedelec<dim>::initialize_embedding ()
+{
+  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)
+      {
+                                         // copy
+        this->prolongation[c].reinit (this->dofs_per_cell,
+                                      this->dofs_per_cell);
+        this->prolongation[c].fill (Matrices::embedding[degree-1][c]);
+                                         // and make sure that the row
+                                         // sum is 0.5 (for usual
+                                         // elements, the row sum must
+                                         // be 1, but here the shape
+                                         // function is multiplied by
+                                         // the inverse of the
+                                         // Jacobian, which introduces
+                                         // a factor of 1/2 when going
+                                         // from mother to child)
+        for (unsigned int row=0; row<this->dofs_per_cell; ++row)
+          {
+            double sum = 0;
+            for (unsigned int col=0; col<this->dofs_per_cell; ++col)
+              sum += this->prolongation[c](row,col);
+            Assert (std::fabs(sum-.5) < 1e-14,
+                    ExcInternalError());
+          };
+      };
+}
+
+
+
+template <int dim>
+void
+FE_Nedelec<dim>::initialize_restriction ()
+{
+  switch (dim)
+    {
+      case 2:   // 2d
+      {
+       switch (degree)
+         {
+           case 1:
+           {
+                                               // this is a strange
+                                               // element, since it is
+                                               // both additive and
+                                               // then it is also
+                                               // not. ideally, we
+                                               // would like to have
+                                               // the value of the
+                                               // shape function on
+                                               // the coarse line to
+                                               // be the mean value of
+                                               // that on the two
+                                               // child ones. thus,
+                                               // one should make it
+                                               // additive. however,
+                                               // additivity only
+                                               // works if an element
+                                               // does not have any
+                                               // continuity
+                                               // requirements, since
+                                               // otherwise degrees of
+                                               // freedom are shared
+                                               // between adjacent
+                                               // elements, and when
+                                               // we make the element
+                                               // additive, that would
+                                               // mean that we end up
+                                               // adding up
+                                               // contributions not
+                                               // only from the child
+                                               // cells of this cell,
+                                               // but also from the
+                                               // child cells of the
+                                               // neighbor, and since
+                                               // we cannot know
+                                               // whether there even
+                                               // exists a neighbor we
+                                               // cannot simply make
+                                               // the element
+                                               // additive.
+                                              //
+                                               // so, until someone
+                                               // comes along with a
+                                               // better alternative,
+                                               // we do the following:
+                                               // make the element
+                                               // non-additive, and
+                                               // simply pick the
+                                               // value of one of the
+                                               // child lines for the
+                                               // value of the mother
+                                               // line (note that we
+                                               // have to multiply by
+                                               // two, since the shape
+                                               // functions scale with
+                                               // the inverse
+                                               // Jacobian). we thus
+                                               // throw away the
+                                               // information of one
+                                               // of the child lines,
+                                               // but there seems to
+                                               // be no other way than
+                                               // that...
+                                               //
+                                               // note: to make things
+                                               // consistent, and
+                                               // restriction
+                                               // independent of the
+                                               // order in which we
+                                               // travel across the
+                                               // cells of the coarse
+                                               // grid, we have to
+                                               // make sure that we
+                                               // take the same small
+                                               // line when visiting
+                                               // its two neighbors,
+                                               // to get the value for
+                                               // the mother line. we
+                                               // take the first line
+                                               // always, in the
+                                               // canonical direction
+                                               // of lines
+              for (unsigned int c=0; c<GeometryInfo<dim>::children_per_cell; ++c)
+                this->restriction[c].reinit (this->dofs_per_cell,
+                                             this->dofs_per_cell);
+              
+             this->restriction[0](0,0) = 2.;
+             this->restriction[1](1,1) = 2.;
+             this->restriction[3](2,2) = 2.;
+             this->restriction[0](3,3) = 2.;
+
+             break;
+           };
+           
+           default:
+           {
+                                              // in case we don't
+                                              // have the matrices
+                                              // (yet), leave them
+                                              // empty. this does not
+                                              // prevent the use of
+                                              // this FE, but will
+                                              // prevent the use of
+                                              // these matrices
+              break;
+           };
+         };
+       
+       break;
+      };
+
+
+      case 3:   // 3d
+      {
+       switch (degree)
+         {
+           case 1:
+           {
+                                              // same principle as in
+                                              // 2d, take one child
+                                              // cell to get at the
+                                              // values of each of
+                                              // the 12 lines
+              for (unsigned int c=0; c<GeometryInfo<dim>::children_per_cell; ++c)
+                this->restriction[c].reinit (this->dofs_per_cell,
+                                             this->dofs_per_cell);
+             this->restriction[0](0,0) = 2.;
+             this->restriction[0](3,3) = 2.;
+             this->restriction[1](1,1) = 2.;
+             this->restriction[3](2,2) = 2.;
+              
+             this->restriction[4](4,4) = 2.;
+             this->restriction[4](7,7) = 2.;
+             this->restriction[5](5,5) = 2.;
+             this->restriction[7](6,6) = 2.;
+              
+             this->restriction[0](8,8) = 2.;
+             this->restriction[1](9,9) = 2.;
+             this->restriction[2](10,10) = 2.;
+             this->restriction[3](11,11) = 2.;
+              
+             break;
+           };
+           
+           default:
+           {
+                                              // in case we don't
+                                              // have the matrices
+                                              // (yet), leave them
+                                              // empty. this does not
+                                              // prevent the use of
+                                              // this FE, but will
+                                              // prevent the use of
+                                              // these matrices
+              break;
+           };
+         };
+       
+       break;
+      };
+      
+      default:
+           Assert (false,ExcNotImplemented());
+    }
+}
+
+
+
 template <int dim>
 void FE_Nedelec<dim>::initialize_unit_support_points ()
 {
index 8a85752502da8824d183e8a3df4e438e0de94c50..aa166ffcc4b924b40118705f618888784995942e 100644 (file)
@@ -522,7 +522,18 @@ FE_Q<dim>::initialize_embedding ()
         this->prolongation[c].reinit (this->dofs_per_cell,
                                       this->dofs_per_cell);
         this->prolongation[c].fill (Matrices::embedding[degree-1][c]);
-      };
+
+                                         // and make sure that the row
+                                         // sum is 1
+        for (unsigned int row=0; row<this->dofs_per_cell; ++row)
+          {
+            double sum = 0;
+            for (unsigned int col=0; col<this->dofs_per_cell; ++col)
+              sum += this->prolongation[c](row,col);
+            Assert (std::fabs(sum-1.) < 1e-14,
+                    ExcInternalError());
+          };
+      };  
 }
 
 
index 30980f23998eac82ffcade12270320dd47655e28..a5106aa7c36293c046c7a42d5cb38987ad2ac972 100644 (file)
 
 #include <cmath>
 
+namespace 
+{
+  std::vector<unsigned int>
+  invert_numbering (const std::vector<unsigned int> &in)
+  {
+    std::vector<unsigned int> out (in.size());
+    for (unsigned int i=0; i<in.size(); ++i)
+      out[in[i]]=i;
+    return out;
+  }
+}
+
 
 
 template <int dim>
@@ -34,284 +46,44 @@ FE_Q_Hierarchical<dim>::FE_Q_Hierarchical (const unsigned int degree)
                                                        false),
                                    std::vector<std::vector<bool> >(FiniteElementData<dim>(get_dpo_vector(degree),1).dofs_per_cell,
                                                                    std::vector<bool>(1,true))),
-                                                                     degree(degree),
-                                                                     renumber(this->dofs_per_cell, 0),
-                                                                     renumber_inverse(this->dofs_per_cell, 0),
-                                                                     face_renumber(this->dofs_per_face, 0),
-                                                                     polynomial_space(Polynomials::Hierarchical<double>::generate_complete_basis(degree))
+               degree(degree),
+               renumber(lexicographic_to_hierarchic_numbering (*this, degree)),
+               renumber_inverse(invert_numbering(renumber)),
+               face_renumber(face_lexicographic_to_hierarchic_numbering (degree)),
+               polynomial_space(Polynomials::Hierarchical<double>::generate_complete_basis(degree))
 {
-                                  // do some internal book-keeping on
-                                  // cells and faces. if in 1d, the
-                                  // face function is empty
-  lexicographic_to_hierarchic_numbering (*this, degree, renumber);
-  face_lexicographic_to_hierarchic_numbering (degree, face_renumber);
-  
-                                  // build inverse of renumbering
-                                  // vector
-  for (unsigned int i=0; i<this->dofs_per_cell; ++i)
-    renumber_inverse[renumber[i]]=i;
-
-                                  // build the dofs_subcell, dofs_cell 
-                                  // matrices for use in building prolongation,
-                                   // restriction, and constraint matrices.
-  const unsigned int dofs_1d = 2*this->dofs_per_vertex + this->dofs_per_line;
-
-  for (unsigned int c=0; c<GeometryInfo<1>::children_per_cell; ++c)
-    {
-      dofs_cell.push_back (FullMatrix<double> (dofs_1d,dofs_1d) );
-      dofs_subcell.push_back (FullMatrix<double> (dofs_1d,dofs_1d) );
-
-      for (unsigned int j=0; j<dofs_1d; ++j)
-       {
-         for (unsigned int k=0; k<dofs_1d; ++k)
-           {   
-                                              // upper diagonal block
-             if ((j<=1) && (k<=1))
-               {
-                 if (((c==0) && (j==0) && (k==0)) || 
-                     ((c==1) && (j==1) && (k==1)))
-                   dofs_cell[c](j,k) = 1.;
-                 else 
-                   dofs_cell[c](j,k) = 0.;
-         
-                 if      (((c==0) && (j==1)) || ((c==1) && (j==0)))
-                   dofs_subcell[c](j,k) = .5;
-                 else if (((c==0) && (k==0)) || ((c==1) && (k==1)))
-                   dofs_subcell[c](j,k) = 1.;
-                 else
-                   dofs_subcell[c](j,k) = 0.;
-               }
-                                              // upper right block
-             else if ((j<=1) && (k>=2))
-               {
-                 if (((c==0) && (j==1) && ((k % 2)==0)) ||
-                     ((c==1) && (j==0) && ((k % 2)==0)))
-                   dofs_subcell[c](j,k) = -1.;
-               }
-                               // lower diagonal block
-             else if ((j>=2) && (k>=2) && (j<=k))
-               {
-                 double factor = 1.;
-                 for (unsigned int i=1; i<=j;++i)
-                   factor *= ((double) (k-i+1))/((double) i);
-                 if (c==0)
-                   {
-                     dofs_subcell[c](j,k) = ((k+j) % 2 == 0) ? 
-                                            std::pow(.5,static_cast<double>(k))*factor :
-                                            -std::pow(.5,static_cast<double>(k))*factor;
-                     dofs_cell[c](j,k) = std::pow(2.,static_cast<double>(j))*factor;
-                   }
-                 else
-                   {
-                     dofs_subcell[c](j,k) = std::pow(.5,static_cast<double>(k))*factor;
-                     dofs_cell[c](j,k) = ((k+j) % 2 == 0) ? 
-                                         std::pow(2.,static_cast<double>(j))*factor :
-                                            -std::pow(2.,static_cast<double>(j))*factor;
-                   }
-               }
-           }
-       }
-    }
-                                  // fill constraint matrices
-  if (dim==2 || dim==3)
-    {
-      this->interface_constraints.reinit ( (dim==2) ? 1 + 2*(degree-1) : 
-                                          5 + 12*(degree-1) + 4*(degree-1)*(degree-1),
-                                          (dim==2) ? (degree+1) : 
-                                          (degree+1)*(degree+1) );
-      switch (dim)
-       {
-         case 2:
-                                                // vertex node
-               for (unsigned int i=0; i<dofs_1d; ++i)
-                 this->interface_constraints(0,i) = dofs_subcell[0](1,i); 
-                                                // edge nodes
-               for (unsigned int c=0; c<GeometryInfo<1>::children_per_cell; ++c)
-                 for (unsigned int i=0; i<dofs_1d; ++i)
-                   for (unsigned int j=2; j<dofs_1d; ++j)
-                     this->interface_constraints(1 + c*(degree-1) + j - 2,i) = 
-                       dofs_subcell[c](j,i);
-               break;
-         case 3:
-               for (unsigned int i=0; i<dofs_1d * dofs_1d; i++)
-                 {
-                                                    // center vertex node        
-                   this->interface_constraints(0,face_renumber[i]) = 
-                     dofs_subcell[0](1,i % dofs_1d) * 
-                     dofs_subcell[0](1,(i - (i % dofs_1d)) / dofs_1d);
-
-                                                    // boundary vertex nodes
-                   this->interface_constraints(1,face_renumber[i]) = 
-                     dofs_subcell[0](1, i % dofs_1d) * 
-                     dofs_subcell[0](0, (i - (i % dofs_1d)) / dofs_1d);
-                   this->interface_constraints(2,face_renumber[i]) = 
-                     dofs_subcell[1](1, i % dofs_1d) * 
-                     dofs_subcell[0](1, (i - (i % dofs_1d)) / dofs_1d);
-                   this->interface_constraints(3,face_renumber[i]) = 
-                     dofs_subcell[1](0, i % dofs_1d) * 
-                     dofs_subcell[1](1, (i - (i % dofs_1d)) / dofs_1d);
-                   this->interface_constraints(4,face_renumber[i]) = 
-                     dofs_subcell[0](0, i % dofs_1d) * 
-                     dofs_subcell[0](1, (i - (i % dofs_1d)) / dofs_1d);
-         
-                                                    // interior edges
-                   for (unsigned int j=0; j<(degree-1); j++)
-                     {
-                       this->interface_constraints(5 + j,face_renumber[i]) = 
-                         dofs_subcell[0](1, i % dofs_1d) * 
-                         dofs_subcell[0](2 + j, (i - (i % dofs_1d)) / dofs_1d);
-                       this->interface_constraints(5 + (degree-1) + j,face_renumber[i]) = 
-                         dofs_subcell[1](2 + j, i % dofs_1d) * 
-                         dofs_subcell[0](1, (i - (i % dofs_1d)) / dofs_1d);
-                       this->interface_constraints(5 + 2*(degree-1) + j,face_renumber[i]) = 
-                         dofs_subcell[0](1,i % dofs_1d) * 
-                         dofs_subcell[1](2 + j, (i - (i % dofs_1d)) / dofs_1d);
-                       this->interface_constraints(5 + 3*(degree-1) + j,face_renumber[i]) = 
-                         dofs_subcell[0](2 + j,i % dofs_1d) * 
-                         dofs_subcell[1](0, (i - (i % dofs_1d)) / dofs_1d);
-                     }
-
-                                                    // boundary edges
-                   for (unsigned int j=0; j<(degree-1); j++)
-                     {
-                                                        // bottom edge 
-                       this->interface_constraints(5 + 4*(degree-1) + j,face_renumber[i]) =
-                         dofs_subcell[0](2 + j, i % dofs_1d) * 
-                         dofs_subcell[0](0,     (i - (i % dofs_1d)) / dofs_1d);
-                       this->interface_constraints(5 + 4*(degree-1) + (degree-1) + j,face_renumber[i]) =
-                         dofs_subcell[1](2 + j, i % dofs_1d) * 
-                         dofs_subcell[0](0,     (i - (i % dofs_1d)) / dofs_1d);
-                                                        // right edge
-                       this->interface_constraints(5 + 4*(degree-1) + 2*(degree-1) + j,face_renumber[i]) =
-                         dofs_subcell[1](1,     i % dofs_1d) * 
-                         dofs_subcell[0](2 + j, (i - (i % dofs_1d)) / dofs_1d);
-                       this->interface_constraints(5 + 4*(degree-1) + 3*(degree-1) + j,face_renumber[i]) =
-                         dofs_subcell[1](1,     i % dofs_1d) * 
-                         dofs_subcell[1](2 + j, (i - (i % dofs_1d)) / dofs_1d);
-                                                        // top edge
-                       this->interface_constraints(5 + 4*(degree-1) + 4*(degree-1) + j,face_renumber[i]) =
-                         dofs_subcell[0](2 + j, i % dofs_1d) * 
-                         dofs_subcell[1](1,     (i - (i % dofs_1d)) / dofs_1d);
-                       this->interface_constraints(5 + 4*(degree-1) + 5*(degree-1) + j,face_renumber[i]) =
-                         dofs_subcell[1](2 + j, i % dofs_1d) * 
-                         dofs_subcell[1](1,     (i - (i % dofs_1d)) / dofs_1d);
-                                                        // left edge
-                       this->interface_constraints(5 + 4*(degree-1) + 6*(degree-1) + j,face_renumber[i]) =
-                         dofs_subcell[0](0,     i % dofs_1d) * 
-                         dofs_subcell[0](2 + j, (i - (i % dofs_1d)) / dofs_1d);
-                       this->interface_constraints(5 + 4*(degree-1) + 7*(degree-1) + j,face_renumber[i]) =
-                         dofs_subcell[0](0,     i % dofs_1d) * 
-                         dofs_subcell[1](2 + j, (i - (i % dofs_1d)) / dofs_1d);
-                     }
-
-                                                    // interior faces
-                   for (unsigned int j=0; j<(degree-1); j++)
-                     for (unsigned int k=0; k<(degree-1); k++)
-                       {
-                                                          // subcell 0
-                         this->interface_constraints(5 + 12*(degree-1) + j + k*(degree-1),face_renumber[i]) =
-                           dofs_subcell[0](2 + j, i % dofs_1d) * 
-                           dofs_subcell[0](2 + k, (i - (i % dofs_1d)) / dofs_1d);
-                                                          // subcell 1
-                         this->interface_constraints(5 + 12*(degree-1) + j + k*(degree-1) + (degree-1)*(degree-1),face_renumber[i]) =
-                           dofs_subcell[1](2 + j, i % dofs_1d) * 
-                           dofs_subcell[0](2 + k, (i - (i % dofs_1d)) / dofs_1d);
-                                                          // subcell 2
-                         this->interface_constraints(5 + 12*(degree-1) + j + k*(degree-1) + 2*(degree-1)*(degree-1),face_renumber[i]) =
-                           dofs_subcell[1](2 + j, i % dofs_1d) * 
-                           dofs_subcell[1](2 + k, (i - (i % dofs_1d)) / dofs_1d);
-                                                          // subcell 3
-                         this->interface_constraints(5 + 12*(degree-1) + j + k*(degree-1) + 3*(degree-1)*(degree-1),face_renumber[i]) =
-                           dofs_subcell[0](2 + j, i % dofs_1d) * 
-                           dofs_subcell[1](2 + k, (i - (i % dofs_1d)) / dofs_1d);
-                       }
-                 }
-               break;
-       }
-    };
-
-                                  // fill prolongation and restriction 
-                                   // matrices
-  if (dim==1)
-    {
-      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 (dofs_subcell[c]);
+                                  // The matrix @p{dofs_cell} contains the 
+                                  // values of the linear functionals of 
+                                  // the master 1d cell applied to the 
+                                  // shape functions of the two 1d subcells.
+                                  // The matrix @p{dofs_subcell} constains
+                                  // the values of the linear functionals 
+                                  // on each 1d subcell applied to the 
+                                  // shape functions on the master 1d 
+                                  // subcell. 
+                                  // We use @p{dofs_cell} and 
+                                  // @p{dofs_subcell} to compute the 
+                                  // @p{prolongation}, @p{restriction} and 
+                                  // @p{interface_constraints} matrices 
+                                  // for all dimensions.
+  std::vector<FullMatrix<double> >
+    dofs_cell (GeometryInfo<1>::children_per_cell,
+              FullMatrix<double> (2*this->dofs_per_vertex + this->dofs_per_line,
+                                  2*this->dofs_per_vertex + this->dofs_per_line));
+  std::vector<FullMatrix<double> >
+    dofs_subcell (GeometryInfo<1>::children_per_cell,
+                 FullMatrix<double> (2*this->dofs_per_vertex + this->dofs_per_line,
+                                     2*this->dofs_per_vertex + this->dofs_per_line));
+                                  // build these fields, as they are
+                                  // needed as auxiliary fields later
+                                  // on
+  build_dofs_cell (dofs_cell, dofs_subcell);
+
+                                  // then use them to initialize
+                                  // other fields
+  initialize_constraints (dofs_subcell);
+  initialize_embedding_and_restriction (dofs_cell, dofs_subcell);
 
-         this->restriction[c].reinit (this->dofs_per_cell,this->dofs_per_cell);
-         this->restriction[c].fill (dofs_cell[c]);
-       }
-    }
-  else if (dim==2 || dim==3)
-    {
-      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->restriction[c].reinit (this->dofs_per_cell,this->dofs_per_cell);
-       }
-                                       // j loops over dofs in the subcell. 
-                                       // These are the rows in the 
-                                       // embedding matrix.
-      for (unsigned int j=0; j<this->dofs_per_cell; ++j)
-       {
-                                          // i loops over the dofs in the master 
-                                          // cell. These are the columns in 
-                                          // the embedding matrix.
-         for (unsigned int i=0; i<this->dofs_per_cell; ++i)
-           {
-             switch (dim)
-               {
-                 case 2:
-                       for (unsigned int c=0; c<GeometryInfo<2>::children_per_cell; ++c)
-                         {
-                           unsigned int c0 = ((c==1) || (c==2)) ? 1 : 0;
-                           unsigned int c1 = ((c==2) || (c==3)) ? 1 : 0;
-
-                           this->prolongation[c](j,i) = 
-                             dofs_subcell[c0](renumber_inverse[j] % dofs_1d,
-                                              renumber_inverse[i] % dofs_1d) *
-                             dofs_subcell[c1]((renumber_inverse[j] - (renumber_inverse[j] % dofs_1d)) / dofs_1d,
-                                              (renumber_inverse[i] - (renumber_inverse[i] % dofs_1d)) / dofs_1d);
-
-                           this->restriction[c](j,i) = 
-                             dofs_cell[c0](renumber_inverse[j] % dofs_1d,
-                                           renumber_inverse[i] % dofs_1d) *
-                             dofs_cell[c1]((renumber_inverse[j] - (renumber_inverse[j] % dofs_1d)) / dofs_1d,
-                                           (renumber_inverse[i] - (renumber_inverse[i] % dofs_1d)) / dofs_1d);
-                         }
-                       break;
-
-                 case 3:
-                       for (unsigned int c=0; c<GeometryInfo<3>::children_per_cell; ++c)
-                         {
-                           unsigned int c0 = ((c==1) || (c==2) || (c==5) || (c==6)) ? 1 : 0;
-                           unsigned int c1 = ((c==4) || (c==5) || (c==6) || (c==7)) ? 1 : 0;
-                           unsigned int c2 = ((c==2) || (c==3) || (c==6) || (c==7)) ? 1 : 0;
-
-                           this->prolongation[c](j,i) = 
-                             dofs_subcell[c0](renumber_inverse[j] % dofs_1d,
-                                              renumber_inverse[i] % dofs_1d) *
-                             dofs_subcell[c1](((renumber_inverse[j] - (renumber_inverse[j] % dofs_1d)) / dofs_1d) % dofs_1d,
-                                              ((renumber_inverse[i] - (renumber_inverse[i] % dofs_1d)) / dofs_1d) % dofs_1d) *
-                             dofs_subcell[c2](((renumber_inverse[j] - (renumber_inverse[j] % dofs_1d)) / dofs_1d - (((renumber_inverse[j] - (renumber_inverse[j] % dofs_1d)) / dofs_1d ) % dofs_1d)) / dofs_1d,
-                                              ((renumber_inverse[i] - (renumber_inverse[i] % dofs_1d)) / dofs_1d - (((renumber_inverse[i] - (renumber_inverse[i] % dofs_1d)) / dofs_1d ) % dofs_1d)) / dofs_1d);
-
-                           this->restriction[c](j,i) = 
-                             dofs_cell[c0](renumber_inverse[j] % dofs_1d,
-                                           renumber_inverse[i] % dofs_1d) *
-                             dofs_cell[c1](((renumber_inverse[j] - (renumber_inverse[j] % dofs_1d)) / dofs_1d) % dofs_1d,
-                                           ((renumber_inverse[i] - (renumber_inverse[i] % dofs_1d)) / dofs_1d) % dofs_1d) *
-                             dofs_cell[c2](((renumber_inverse[j] - (renumber_inverse[j] % dofs_1d)) / dofs_1d - (((renumber_inverse[j] - (renumber_inverse[j] % dofs_1d)) / dofs_1d ) % dofs_1d)) / dofs_1d,
-                                           ((renumber_inverse[i] - (renumber_inverse[i] % dofs_1d)) / dofs_1d - (((renumber_inverse[i] - (renumber_inverse[i] % dofs_1d)) / dofs_1d ) % dofs_1d)) / dofs_1d);
-                         }
-                       break;
-               }
-           }
-       }
-    }
-  else
-    Assert (false, ExcNotImplemented());
                                   // finally fill in support points
                                   // on cell and face
   initialize_unit_support_points ();
@@ -405,6 +177,298 @@ FE_Q_Hierarchical<dim>::shape_grad_grad_component (const unsigned int i,
 //----------------------------------------------------------------------
 
 
+template <int dim>
+void
+FE_Q_Hierarchical<dim>::build_dofs_cell (std::vector<FullMatrix<double> > &dofs_cell,
+                                        std::vector<FullMatrix<double> > &dofs_subcell) const
+{
+  const unsigned int dofs_1d = 2*this->dofs_per_vertex + this->dofs_per_line;
+
+  for (unsigned int c=0; c<GeometryInfo<1>::children_per_cell; ++c)
+    for (unsigned int j=0; j<dofs_1d; ++j)
+      for (unsigned int k=0; k<dofs_1d; ++k)
+       {
+                                          // upper diagonal block
+         if ((j<=1) && (k<=1))
+           {
+             if (((c==0) && (j==0) && (k==0)) || 
+                 ((c==1) && (j==1) && (k==1)))
+               dofs_cell[c](j,k) = 1.;
+             else 
+               dofs_cell[c](j,k) = 0.;
+             
+             if      (((c==0) && (j==1)) || ((c==1) && (j==0)))
+               dofs_subcell[c](j,k) = .5;
+             else if (((c==0) && (k==0)) || ((c==1) && (k==1)))
+               dofs_subcell[c](j,k) = 1.;
+             else
+               dofs_subcell[c](j,k) = 0.;
+           }
+                                          // upper right block
+         else if ((j<=1) && (k>=2))
+           {
+             if (((c==0) && (j==1) && ((k % 2)==0)) ||
+                 ((c==1) && (j==0) && ((k % 2)==0)))
+               dofs_subcell[c](j,k) = -1.;
+           }
+                               // lower diagonal block
+         else if ((j>=2) && (k>=2) && (j<=k))
+           {
+             double factor = 1.;
+             for (unsigned int i=1; i<=j;++i)
+               factor *= ((double) (k-i+1))/((double) i);
+             if (c==0)
+               {
+                 dofs_subcell[c](j,k) = ((k+j) % 2 == 0) ? 
+                                        std::pow(.5,static_cast<double>(k))*factor :
+                                        -std::pow(.5,static_cast<double>(k))*factor;
+                 dofs_cell[c](j,k) = std::pow(2.,static_cast<double>(j))*factor;
+               }
+             else
+               {
+                 dofs_subcell[c](j,k) = std::pow(.5,static_cast<double>(k))*factor;
+                 dofs_cell[c](j,k) = ((k+j) % 2 == 0) ? 
+                                     std::pow(2.,static_cast<double>(j))*factor :
+                                        -std::pow(2.,static_cast<double>(j))*factor;
+               }
+           }
+       }
+}
+
+
+
+template <int dim>
+void
+FE_Q_Hierarchical<dim>::
+initialize_constraints (const std::vector<FullMatrix<double> > &dofs_subcell)
+{
+  const unsigned int dofs_1d = 2*this->dofs_per_vertex + this->dofs_per_line;
+
+  this->interface_constraints
+    .TableBase<2,double>::reinit (this->interface_constraints_size());
+
+  switch (dim)
+    {
+      case 1:
+      {
+                                        // no constraints in 1d
+       break;
+      }
+       
+      case 2:
+      {
+                                        // vertex node
+       for (unsigned int i=0; i<dofs_1d; ++i)
+         this->interface_constraints(0,i) = dofs_subcell[0](1,i); 
+                                        // edge nodes
+       for (unsigned int c=0; c<GeometryInfo<1>::children_per_cell; ++c)
+         for (unsigned int i=0; i<dofs_1d; ++i)
+           for (unsigned int j=2; j<dofs_1d; ++j)
+             this->interface_constraints(1 + c*(degree-1) + j - 2,i) = 
+               dofs_subcell[c](j,i);
+       break;
+      }
+       
+      case 3:
+      {
+       for (unsigned int i=0; i<dofs_1d * dofs_1d; i++)
+         {
+                                            // center vertex node        
+           this->interface_constraints(0,face_renumber[i]) = 
+             dofs_subcell[0](1,i % dofs_1d) * 
+             dofs_subcell[0](1,(i - (i % dofs_1d)) / dofs_1d);
+
+                                                // boundary vertex nodes
+           this->interface_constraints(1,face_renumber[i]) = 
+             dofs_subcell[0](1, i % dofs_1d) * 
+             dofs_subcell[0](0, (i - (i % dofs_1d)) / dofs_1d);
+           this->interface_constraints(2,face_renumber[i]) = 
+             dofs_subcell[1](1, i % dofs_1d) * 
+             dofs_subcell[0](1, (i - (i % dofs_1d)) / dofs_1d);
+           this->interface_constraints(3,face_renumber[i]) = 
+             dofs_subcell[1](0, i % dofs_1d) * 
+             dofs_subcell[1](1, (i - (i % dofs_1d)) / dofs_1d);
+           this->interface_constraints(4,face_renumber[i]) = 
+             dofs_subcell[0](0, i % dofs_1d) * 
+             dofs_subcell[0](1, (i - (i % dofs_1d)) / dofs_1d);
+         
+                                                // interior edges
+           for (unsigned int j=0; j<(degree-1); j++)
+             {
+               this->interface_constraints(5 + j,face_renumber[i]) = 
+                 dofs_subcell[0](1, i % dofs_1d) * 
+                 dofs_subcell[0](2 + j, (i - (i % dofs_1d)) / dofs_1d);
+               this->interface_constraints(5 + (degree-1) + j,face_renumber[i]) = 
+                 dofs_subcell[1](2 + j, i % dofs_1d) * 
+                 dofs_subcell[0](1, (i - (i % dofs_1d)) / dofs_1d);
+               this->interface_constraints(5 + 2*(degree-1) + j,face_renumber[i]) = 
+                 dofs_subcell[0](1,i % dofs_1d) * 
+                 dofs_subcell[1](2 + j, (i - (i % dofs_1d)) / dofs_1d);
+               this->interface_constraints(5 + 3*(degree-1) + j,face_renumber[i]) = 
+                 dofs_subcell[0](2 + j,i % dofs_1d) * 
+                 dofs_subcell[1](0, (i - (i % dofs_1d)) / dofs_1d);
+             }
+
+                                            // boundary edges
+           for (unsigned int j=0; j<(degree-1); j++)
+             {
+                                                // bottom edge 
+               this->interface_constraints(5 + 4*(degree-1) + j,face_renumber[i]) =
+                 dofs_subcell[0](2 + j, i % dofs_1d) * 
+                 dofs_subcell[0](0,     (i - (i % dofs_1d)) / dofs_1d);
+               this->interface_constraints(5 + 4*(degree-1) + (degree-1) + j,face_renumber[i]) =
+                 dofs_subcell[1](2 + j, i % dofs_1d) * 
+                 dofs_subcell[0](0,     (i - (i % dofs_1d)) / dofs_1d);
+                                                // right edge
+               this->interface_constraints(5 + 4*(degree-1) + 2*(degree-1) + j,face_renumber[i]) =
+                 dofs_subcell[1](1,     i % dofs_1d) * 
+                 dofs_subcell[0](2 + j, (i - (i % dofs_1d)) / dofs_1d);
+               this->interface_constraints(5 + 4*(degree-1) + 3*(degree-1) + j,face_renumber[i]) =
+                 dofs_subcell[1](1,     i % dofs_1d) * 
+                 dofs_subcell[1](2 + j, (i - (i % dofs_1d)) / dofs_1d);
+                                                // top edge
+               this->interface_constraints(5 + 4*(degree-1) + 4*(degree-1) + j,face_renumber[i]) =
+                 dofs_subcell[0](2 + j, i % dofs_1d) * 
+                 dofs_subcell[1](1,     (i - (i % dofs_1d)) / dofs_1d);
+               this->interface_constraints(5 + 4*(degree-1) + 5*(degree-1) + j,face_renumber[i]) =
+                 dofs_subcell[1](2 + j, i % dofs_1d) * 
+                 dofs_subcell[1](1,     (i - (i % dofs_1d)) / dofs_1d);
+                                                // left edge
+               this->interface_constraints(5 + 4*(degree-1) + 6*(degree-1) + j,face_renumber[i]) =
+                 dofs_subcell[0](0,     i % dofs_1d) * 
+                 dofs_subcell[0](2 + j, (i - (i % dofs_1d)) / dofs_1d);
+               this->interface_constraints(5 + 4*(degree-1) + 7*(degree-1) + j,face_renumber[i]) =
+                 dofs_subcell[0](0,     i % dofs_1d) * 
+                 dofs_subcell[1](2 + j, (i - (i % dofs_1d)) / dofs_1d);
+             }
+
+                                            // interior faces
+           for (unsigned int j=0; j<(degree-1); j++)
+             for (unsigned int k=0; k<(degree-1); k++)
+               {
+                                                  // subcell 0
+                 this->interface_constraints(5 + 12*(degree-1) + j + k*(degree-1),face_renumber[i]) =
+                   dofs_subcell[0](2 + j, i % dofs_1d) * 
+                   dofs_subcell[0](2 + k, (i - (i % dofs_1d)) / dofs_1d);
+                                                  // subcell 1
+                 this->interface_constraints(5 + 12*(degree-1) + j + k*(degree-1) + (degree-1)*(degree-1),face_renumber[i]) =
+                   dofs_subcell[1](2 + j, i % dofs_1d) * 
+                   dofs_subcell[0](2 + k, (i - (i % dofs_1d)) / dofs_1d);
+                                                  // subcell 2
+                 this->interface_constraints(5 + 12*(degree-1) + j + k*(degree-1) + 2*(degree-1)*(degree-1),face_renumber[i]) =
+                   dofs_subcell[1](2 + j, i % dofs_1d) * 
+                   dofs_subcell[1](2 + k, (i - (i % dofs_1d)) / dofs_1d);
+                                                  // subcell 3
+                 this->interface_constraints(5 + 12*(degree-1) + j + k*(degree-1) + 3*(degree-1)*(degree-1),face_renumber[i]) =
+                   dofs_subcell[0](2 + j, i % dofs_1d) * 
+                   dofs_subcell[1](2 + k, (i - (i % dofs_1d)) / dofs_1d);
+               }
+         }
+       break;
+      }
+       
+      default:
+           Assert (false, ExcNotImplemented());
+    }
+}
+
+
+
+template <int dim>
+void
+FE_Q_Hierarchical<dim>::
+initialize_embedding_and_restriction (const std::vector<FullMatrix<double> > &dofs_cell,
+                                     const std::vector<FullMatrix<double> > &dofs_subcell)
+{
+  const unsigned int dofs_1d = 2*this->dofs_per_vertex + this->dofs_per_line;
+
+  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->restriction[c].reinit (this->dofs_per_cell, this->dofs_per_cell);
+    }
+
+                                  // the 1d case is particularly
+                                  // simple, so special case it:
+  if (dim==1)
+    {
+      for (unsigned int c=0; c<GeometryInfo<dim>::children_per_cell; ++c)
+       {
+         this->prolongation[c].fill (dofs_subcell[c]);
+         this->restriction[c].fill (dofs_cell[c]);
+       }
+      return;
+    }
+
+                                  // for higher dimensions, things
+                                  // are a little more tricky:
+  
+                                  // j loops over dofs in the
+                                  // subcell.  These are the rows in
+                                  // the embedding matrix.
+                                  //
+                                  // i loops over the dofs in the
+                                  // master cell. These are the
+                                  // columns in the embedding matrix.
+  for (unsigned int j=0; j<this->dofs_per_cell; ++j)
+    for (unsigned int i=0; i<this->dofs_per_cell; ++i)
+      switch (dim)
+       {
+         case 2:
+         {
+           for (unsigned int c=0; c<GeometryInfo<2>::children_per_cell; ++c)
+             {
+               unsigned int c0 = ((c==1) || (c==2)) ? 1 : 0;
+               unsigned int c1 = ((c==2) || (c==3)) ? 1 : 0;
+
+               this->prolongation[c](j,i) = 
+                 dofs_subcell[c0](renumber_inverse[j] % dofs_1d,
+                                  renumber_inverse[i] % dofs_1d) *
+                 dofs_subcell[c1]((renumber_inverse[j] - (renumber_inverse[j] % dofs_1d)) / dofs_1d,
+                                  (renumber_inverse[i] - (renumber_inverse[i] % dofs_1d)) / dofs_1d);
+
+               this->restriction[c](j,i) = 
+                 dofs_cell[c0](renumber_inverse[j] % dofs_1d,
+                               renumber_inverse[i] % dofs_1d) *
+                 dofs_cell[c1]((renumber_inverse[j] - (renumber_inverse[j] % dofs_1d)) / dofs_1d,
+                               (renumber_inverse[i] - (renumber_inverse[i] % dofs_1d)) / dofs_1d);
+             }
+           break;
+         }
+            
+         case 3:
+         {
+           for (unsigned int c=0; c<GeometryInfo<3>::children_per_cell; ++c)
+             {
+               unsigned int c0 = ((c==1) || (c==2) || (c==5) || (c==6)) ? 1 : 0;
+               unsigned int c1 = ((c==4) || (c==5) || (c==6) || (c==7)) ? 1 : 0;
+               unsigned int c2 = ((c==2) || (c==3) || (c==6) || (c==7)) ? 1 : 0;
+
+               this->prolongation[c](j,i) = 
+                 dofs_subcell[c0](renumber_inverse[j] % dofs_1d,
+                                  renumber_inverse[i] % dofs_1d) *
+                 dofs_subcell[c1](((renumber_inverse[j] - (renumber_inverse[j] % dofs_1d)) / dofs_1d) % dofs_1d,
+                                  ((renumber_inverse[i] - (renumber_inverse[i] % dofs_1d)) / dofs_1d) % dofs_1d) *
+                 dofs_subcell[c2](((renumber_inverse[j] - (renumber_inverse[j] % dofs_1d)) / dofs_1d - (((renumber_inverse[j] - (renumber_inverse[j] % dofs_1d)) / dofs_1d ) % dofs_1d)) / dofs_1d,
+                                  ((renumber_inverse[i] - (renumber_inverse[i] % dofs_1d)) / dofs_1d - (((renumber_inverse[i] - (renumber_inverse[i] % dofs_1d)) / dofs_1d ) % dofs_1d)) / dofs_1d);
+
+               this->restriction[c](j,i) = 
+                 dofs_cell[c0](renumber_inverse[j] % dofs_1d,
+                               renumber_inverse[i] % dofs_1d) *
+                 dofs_cell[c1](((renumber_inverse[j] - (renumber_inverse[j] % dofs_1d)) / dofs_1d) % dofs_1d,
+                               ((renumber_inverse[i] - (renumber_inverse[i] % dofs_1d)) / dofs_1d) % dofs_1d) *
+                 dofs_cell[c2](((renumber_inverse[j] - (renumber_inverse[j] % dofs_1d)) / dofs_1d - (((renumber_inverse[j] - (renumber_inverse[j] % dofs_1d)) / dofs_1d ) % dofs_1d)) / dofs_1d,
+                               ((renumber_inverse[i] - (renumber_inverse[i] % dofs_1d)) / dofs_1d - (((renumber_inverse[i] - (renumber_inverse[i] % dofs_1d)) / dofs_1d ) % dofs_1d)) / dofs_1d);
+             }
+           break;
+         }
+
+         default:
+               Assert (false, ExcNotImplemented());
+       }
+}
+
+
 
 template <int dim>
 void FE_Q_Hierarchical<dim>::initialize_unit_support_points ()
@@ -548,12 +612,13 @@ FE_Q_Hierarchical<dim>::get_dpo_vector(const unsigned int deg)
 
 
 template <int dim>
-void
-FE_Q_Hierarchical<dim>::lexicographic_to_hierarchic_numbering (
-  const FiniteElementData<dim> &fe_data,
-  const unsigned int            degree,
-  std::vector<unsigned int>    &renumber)
+std::vector<unsigned int>
+FE_Q_Hierarchical<dim>::
+lexicographic_to_hierarchic_numbering (const FiniteElementData<dim> &fe_data,
+                                      const unsigned int            degree)
 {
+  std::vector<unsigned int> renumber (fe_data.dofs_per_cell, 0);
+  
   const unsigned int n = degree+1;
 
 
@@ -624,7 +689,8 @@ FE_Q_Hierarchical<dim>::lexicographic_to_hierarchic_numbering (
              default:
                    Assert(false, ExcNotImplemented());
            }
-      
+
+         Assert (index < renumber.size(), ExcInternalError());
          renumber[index] = i;
        }
     };
@@ -709,6 +775,7 @@ FE_Q_Hierarchical<dim>::lexicographic_to_hierarchic_numbering (
          for (unsigned int jx = 2; jx<=degree ;++jx)
            {
              unsigned int tensorindex = tensorstart + jx * incr;
+             Assert (tensorindex < renumber.size(), ExcInternalError());
              renumber[tensorindex] = index++;
            }
        }
@@ -758,6 +825,7 @@ FE_Q_Hierarchical<dim>::lexicographic_to_hierarchic_numbering (
              {
                unsigned int tensorindex = tensorstart
                                           + jx * incx + jy * incy;
+               Assert (tensorindex < renumber.size(), ExcInternalError());
                renumber[tensorindex] = index++;
              }
        }
@@ -772,34 +840,36 @@ FE_Q_Hierarchical<dim>::lexicographic_to_hierarchic_numbering (
                for (unsigned int jx = 2; 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>
-void
-FE_Q_Hierarchical<dim>::face_lexicographic_to_hierarchic_numbering (
-  const unsigned int         degree,
-  std::vector<unsigned int> &numbering)
+std::vector<unsigned int>
+FE_Q_Hierarchical<dim>::
+face_lexicographic_to_hierarchic_numbering (const unsigned int degree)
 {
   FiniteElementData<dim-1> fe_data(FE_Q_Hierarchical<dim-1>::get_dpo_vector(degree),1);
-  FE_Q_Hierarchical<dim-1>::lexicographic_to_hierarchic_numbering (fe_data,
-                                                                  degree,
-                                                                  numbering); 
+  return FE_Q_Hierarchical<dim-1>::lexicographic_to_hierarchic_numbering (fe_data,
+                                                                         degree); 
 }
 
 
 #if (deal_II_dimension == 1)
 
 template <>
-void
-FE_Q_Hierarchical<1>::face_lexicographic_to_hierarchic_numbering (const unsigned int,
-                                                                 std::vector<unsigned int>&)
-{}
+std::vector<unsigned int>
+FE_Q_Hierarchical<1>::face_lexicographic_to_hierarchic_numbering (const unsigned int)
+{
+  return std::vector<unsigned int> ();
+}
 
 #endif
 

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.