]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Rather pervasive clean-ups. Move some functions to private sections, fix two bugs...
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 7 Mar 2001 13:57:11 +0000 (13:57 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 7 Mar 2001 13:57:11 +0000 (13:57 +0000)
git-svn-id: https://svn.dealii.org/trunk@4144 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/fe/fe_q.h
deal.II/deal.II/source/fe/fe_q.cc
deal.II/deal.II/source/fe/mat_q.1
deal.II/deal.II/source/fe/mat_q.2
deal.II/deal.II/source/fe/mat_q.3

index 4314c96e47d060e9cef3597a04e83c48fb35914e..90fd7ca6ce8e5e0af6cfdccba25e168643e69809 100644 (file)
 template <int dim> class TensorProductPolynomials;
 template <int dim> class MappingQ;
 
+
+
 /**
  * Tensor product elements based on equidistant support points.
+//TODO: Document node numbering etc. copy from old documentation 
  */
 template <int dim>
 class FE_Q : public FiniteElement<dim>
@@ -37,18 +40,6 @@ class FE_Q : public FiniteElement<dim>
                                      */
     ~FE_Q ();
 
-                                    /**
-                                     * Compute flags for initial
-                                     * update only.
-                                     */
-    virtual UpdateFlags update_once (UpdateFlags flags) const;
-  
-                                    /**
-                                     * Compute flags for update on
-                                     * each cell.
-                                     */
-    virtual UpdateFlags update_each (UpdateFlags flags) const;
-
                                     /**
                                      * Return the support points of
                                      * the trial functions on the
@@ -114,8 +105,7 @@ class FE_Q : public FiniteElement<dim>
                                      */
     virtual unsigned int memory_consumption () const;
 
-  protected:
-
+  protected:    
                                     /**
                                      * @p{clone} function instead of
                                      * a copy constructor.
@@ -178,6 +168,74 @@ class FE_Q : public FiniteElement<dim>
                            FEValuesData<dim>& data) const ;
 
   private:
+
+                                    /**
+                                     * Declare a nested class which
+                                     * will has static definitions of
+                                     * various matrices such as
+                                     * constraint and embedding
+                                     * matrices. The definition of
+                                     * the various static fields are
+                                     * in the files @p{mat_q.[123]}
+                                     * in the source directory.
+                                     */
+    struct Matrices {
+                                        /**
+                                         * Embedding matrices. For
+                                         * each element type (the
+                                         * first index) there are as
+                                         * many embedding matrices as
+                                         * there are children per
+                                         * cell. The first index
+                                         * starts with linear
+                                         * elements and goes up in
+                                         * polynomial degree. The
+                                         * array may grow in the
+                                         * future with the number of
+                                         * elements for which these
+                                         * matrices have been
+                                         * computed. If for some
+                                         * element, the matrices have
+                                         * not been computed then you
+                                         * may use the element
+                                         * nevertheless but can not
+                                         * access the respective
+                                         * fields.
+                                         */
+       static const double * const
+       embedding[][GeometryInfo<dim>::children_per_cell];
+
+                                        /**
+                                         * Number of elements (first
+                                         * index) the above field
+                                         * has. Equals the highest
+                                         * polynomial degree for
+                                         * which the embedding
+                                         * matrices have been
+                                         * computed.
+                                         */
+       static const unsigned int n_embedding_matrices;
+
+                                        /**
+                                         * As the
+                                         * @p{embedding_matrices}
+                                         * field, but for the
+                                         * interface constraints. One
+                                         * for each element for which
+                                         * it has been computed.
+                                         */
+       static const double * const constraint_matrices[];
+
+                                        /**
+                                         * Like
+                                         * @p{n_embedding_matrices},
+                                         * but for the number of
+                                         * interface constraint
+                                         * matrices.
+                                         */
+       static const unsigned int n_constraint_matrices;
+    };
+    
                                     /**
                                      * Only for internal use. Its
                                      * full name is
@@ -221,8 +279,8 @@ class FE_Q : public FiniteElement<dim>
                                      * @p{FiniteElementData} argument.
                                      */
     static void build_renumbering (const FiniteElementData<dim> &fe_data,
-                                  unsigned int degree,
-                                  std::vector<unsigned int>numbering);
+                                  const unsigned int            degree,
+                                  std::vector<unsigned int>    &numbering);
 
                                     /**
                                      * Map tensor product data to
@@ -235,18 +293,34 @@ class FE_Q : public FiniteElement<dim>
                                      *
                                      * It is used to compute the
                                      * order of face support points.
+                                     *
+                                     * In 1d, this function does
+                                     * nothing.
                                      */
     static void build_face_renumbering (const FiniteElementData<dim-1> &fe_data,
-                                       unsigned int degree,
-                                       std::vector<unsigned int>numbering);
+                                       const unsigned int              degree,
+                                       std::vector<unsigned int>      &numbering);
 
                                     /**
                                      * Compute support points.
                                      */
-    static void compute_support_points (std::vector<Point<dim> >& support_points,
-                                       unsigned int degree,
-                                       const std::vector<unsigned int>& renumber);
+    static void compute_support_points (const unsigned int               degree,
+                                       const std::vector<unsigned int> &renumber,
+                                       std::vector<Point<dim> >        &support_points);
 
+
+                                    /**
+                                     * Compute flags for initial
+                                     * update only.
+                                     */
+    virtual UpdateFlags update_once (UpdateFlags flags) const;
+  
+                                    /**
+                                     * Compute flags for update on
+                                     * each cell.
+                                     */
+    virtual UpdateFlags update_each (UpdateFlags flags) const;
+    
                                     /**
                                      * Degree of the polynomials.
                                      */  
index a1f585675a542e899b46eb196fa3a4666df41bda..29cb80fd3ab7de8a21d470b22e5ba9a86afc8786 100644 (file)
@@ -24,7 +24,8 @@
 #include <fe/fe_values.h>
 
 
-// Embedding matrices (produced by tests/fe/embedding with postprocessing
+// Include embedding matrices (produced by tests/fe/embedding with
+// postprocessing
 
 #if (deal_II_dimension == 1)
 #include "mat_q.1"
 #endif
 
 
+//TODO: move build_renumbering to mapping class
+
 template <int dim>
-FE_Q<dim>::FE_Q (unsigned int degree)
+FE_Q<dim>::FE_Q (const unsigned int degree)
                :
                FiniteElement<dim> (FiniteElementData<dim>(get_dpo_vector(degree),1),
                                    std::vector<bool> (1,false)),
@@ -49,539 +52,537 @@ FE_Q<dim>::FE_Q (unsigned int degree)
                face_renumber(dofs_per_face, 0),
                poly(0)
 {
+                                  // Q0 elements cannot be
+                                  // continuous, use FE_DGQ<dim>(0)
+                                  // instead
+  Assert (degree>0, ExcNotImplemented());
   std::vector<LagrangeEquidistant> v;
   for (unsigned int i=0;i<=degree;++i)
     v.push_back(LagrangeEquidistant(degree,i));
   
   poly = new TensorProductPolynomials<dim> (v);
-  build_renumbering (*this, degree, renumber);
 
-#if (deal_II_dimension > 1)
+                                  // do some internal book-keeping
+  build_renumbering (*this, degree, renumber);
+#if deal_II_dimension > 1
   build_face_renumbering (FiniteElementData<dim-1>(FE_Q<dim-1>::get_dpo_vector(degree),1),
                          degree,
                          face_renumber);
 #endif
 
-  Assert (degree <= 10, ExcNotImplemented());
-  
-  if ((dim>1) && (constraint_matrices[degree-1] != 0))
-    interface_constraints.fill (constraint_matrices[degree-1]);
+                                  // copy constraint matrices if they
+                                  // are defined. otherwise set them
+                                  // to invalid size
+  if ((dim>1) && (degree<Matrices::n_constraint_matrices+1))
+    interface_constraints.fill (Matrices::constraint_matrices[degree-1]);
   else
     interface_constraints.reinit(0,0);
-  
-  if (embedding0[degree-1] != 0)
-    {
-      prolongation[0].fill (embedding0[degree-1]);
-      prolongation[1].fill (embedding1[degree-1]);
-      if (dim>1)
-       {
-         prolongation[2].fill (embedding2[degree-1]);
-         prolongation[3].fill (embedding3[degree-1]);
-       }
-      if (dim>2)
-       {
-         prolongation[4].fill (embedding4[degree-1]);
-         prolongation[5].fill (embedding5[degree-1]);      
-         prolongation[6].fill (embedding6[degree-1]);
-         prolongation[7].fill (embedding7[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)
+      prolongation[c].fill (Matrices::embedding[degree-1][c]);
   else
     for (unsigned int i=0; i<GeometryInfo<dim>::children_per_cell;++i)
       prolongation[i].reinit(0,0);
-  
+
+                                  // then fill restriction
+                                  // matrices. they are hardcoded for
+                                  // the first few elements
   switch (dim)
     {
-    case 1:
-      switch (degree)
-       {
-       case 1:
-         restriction[0](0,0) = 1;
-         restriction[1](1,1) = 1;
-         break;
-       case 2:
-         restriction[0](0,0) = 1;
-         restriction[0](2,1) = 1;
-         restriction[1](1,1) = 1;
-         restriction[1](1,1) = 1;
-         break;
-       case 3:
-         restriction[0](0,0) = 1;
-         restriction[0](2,3) = 1;
-         restriction[1](1,1) = 1;
-         restriction[1](3,2) = 1;
-         break;
-       case 4:
-         restriction[0](0,0) = 1;
-         restriction[0](2,3) = 1;
-         restriction[0](3,1) = 1;
-         restriction[1](1,1) = 1;
-         restriction[1](3,0) = 1;
-         restriction[1](4,3) = 1;
-         break;
-       default:
-         for (unsigned int i=0; i<GeometryInfo<dim>::children_per_cell;++i)
-           restriction[i].reinit(0,0);
-       }
-      break;
-    case 2:
-      switch (degree)
-       {
-       case 1:
-         restriction[0](0,0) = 1;
-         restriction[1](1,1) = 1;
-         restriction[2](2,2) = 1;
-         restriction[3](3,3) = 1;
-         break;
-       case 2:
-         restriction[0](0,0) = 1;
-         restriction[0](4,1) = 1;
-         restriction[0](7,3) = 1;
-         restriction[0](8,2) = 1;
-         restriction[1](1,1) = 1;
-         restriction[1](4,0) = 1;
-         restriction[1](5,2) = 1;
-         restriction[1](8,3) = 1;
-         restriction[2](2,2) = 1;
-         restriction[2](5,1) = 1;
-         restriction[2](6,3) = 1;
-         restriction[2](8,0) = 1;
-         restriction[3](3,3) = 1;
-         restriction[3](6,2) = 1;
-         restriction[3](7,0) = 1;
-         restriction[3](8,1) = 1;
-         break;
-       case 3:
-         restriction[0](0,0) = 1;
-         restriction[0](4,5) = 1;
-         restriction[0](10,11) = 1;
-         restriction[0](12,15) = 1;
-         restriction[1](1,1) = 1;
-         restriction[1](5,4) = 1;
-         restriction[1](6,7) = 1;
-         restriction[1](13,14) = 1;
-         restriction[2](2,2) = 1;
-         restriction[2](7,6) = 1;
-         restriction[2](9,8) = 1;
-         restriction[2](15,12) = 1;
-         restriction[3](3,3) = 1;
-         restriction[3](8,9) = 1;
-         restriction[3](11,10) = 1;
-         restriction[3](14,13) = 1;
-         break;
-       case 4:
-         restriction[0](0,0) = 1;
-         restriction[0](4,5) = 1;
-         restriction[0](5,1) = 1;
-         restriction[0](13,14) = 1;
-         restriction[0](14,3) = 1;
-         restriction[0](16,20) = 1;
-         restriction[0](17,8) = 1;
-         restriction[0](19,11) = 1;
-         restriction[0](20,2) = 1;
-         restriction[1](1,1) = 1;
-         restriction[1](5,0) = 1;
-         restriction[1](6,5) = 1;
-         restriction[1](7,8) = 1;
-         restriction[1](8,2) = 1;
-         restriction[1](17,14) = 1;
-         restriction[1](18,20) = 1;
-         restriction[1](20,3) = 1;
-         restriction[1](21,11) = 1;
-         restriction[2](2,2) = 1;
-         restriction[2](8,1) = 1;
-         restriction[2](9,8) = 1;
-         restriction[2](11,3) = 1;
-         restriction[2](12,11) = 1;
-         restriction[2](20,0) = 1;
-         restriction[2](21,5) = 1;
-         restriction[2](23,14) = 1;
-         restriction[2](24,20) = 1;
-         restriction[3](3,3) = 1;
-         restriction[3](10,11) = 1;
-         restriction[3](11,2) = 1;
-         restriction[3](14,0) = 1;
-         restriction[3](15,14) = 1;
-         restriction[3](19,5) = 1;
-         restriction[3](20,1) = 1;
-         restriction[3](22,20) = 1;
-         restriction[3](23,8) = 1;
-         break;
-       default:
-         for (unsigned int i=0; i<GeometryInfo<dim>::children_per_cell;++i)
-           restriction[i].reinit(0,0);
-       }
-      break;
-    case 3:
-      switch (degree)
-       {
-       case 1:
-         restriction[0](0,0) = 1;
-         restriction[1](1,1) = 1;
-         restriction[2](2,2) = 1;
-         restriction[3](3,3) = 1;
-         restriction[4](4,4) = 1;
-         restriction[5](5,5) = 1;
-         restriction[6](6,6) = 1;
-         restriction[7](7,7) = 1;
-         break;
-       case 2:
-         restriction[0](0,0) = 1;
-         restriction[0](8,1) = 1;
-         restriction[0](11,3) = 1;
-         restriction[0](16,4) = 1;
-         restriction[0](20,2) = 1;
-         restriction[0](22,5) = 1;
-         restriction[0](25,7) = 1;
-         restriction[0](26,6) = 1;
-         restriction[1](1,1) = 1;
-         restriction[1](8,0) = 1;
-         restriction[1](9,2) = 1;
-         restriction[1](17,5) = 1;
-         restriction[1](20,3) = 1;
-         restriction[1](22,4) = 1;
-         restriction[1](23,6) = 1;
-         restriction[1](26,7) = 1;
-         restriction[2](2,2) = 1;
-         restriction[2](9,1) = 1;
-         restriction[2](10,3) = 1;
-         restriction[2](18,6) = 1;
-         restriction[2](20,0) = 1;
-         restriction[2](23,5) = 1;
-         restriction[2](24,7) = 1;
-         restriction[2](26,4) = 1;
-         restriction[3](3,3) = 1;
-         restriction[3](10,2) = 1;
-         restriction[3](11,0) = 1;
-         restriction[3](19,7) = 1;
-         restriction[3](20,1) = 1;
-         restriction[3](24,6) = 1;
-         restriction[3](25,4) = 1;
-         restriction[3](26,5) = 1;
-         restriction[4](4,4) = 1;
-         restriction[4](12,5) = 1;
-         restriction[4](15,7) = 1;
-         restriction[4](16,0) = 1;
-         restriction[4](21,6) = 1;
-         restriction[4](22,1) = 1;
-         restriction[4](25,3) = 1;
-         restriction[4](26,2) = 1;
-         restriction[5](5,5) = 1;
-         restriction[5](12,4) = 1;
-         restriction[5](13,6) = 1;
-         restriction[5](17,1) = 1;
-         restriction[5](21,7) = 1;
-         restriction[5](22,0) = 1;
-         restriction[5](23,2) = 1;
-         restriction[5](26,3) = 1;
-         restriction[6](6,6) = 1;
-         restriction[6](13,5) = 1;
-         restriction[6](14,7) = 1;
-         restriction[6](18,2) = 1;
-         restriction[6](21,4) = 1;
-         restriction[6](23,1) = 1;
-         restriction[6](24,3) = 1;
-         restriction[6](26,0) = 1;
-         restriction[7](7,7) = 1;
-         restriction[7](14,6) = 1;
-         restriction[7](15,4) = 1;
-         restriction[7](19,3) = 1;
-         restriction[7](21,5) = 1;
-         restriction[7](24,2) = 1;
-         restriction[7](25,0) = 1;
-         restriction[7](26,1) = 1;
-         break;
-       case 3:
-         restriction[0](0,0) = 1;
-         restriction[0](8,9) = 1;
-         restriction[0](14,15) = 1;
-         restriction[0](24,25) = 1;
-         restriction[0](32,35) = 1;
-         restriction[0](40,43) = 1;
-         restriction[0](52,55) = 1;
-         restriction[0](56,63) = 1;
-         restriction[1](1,1) = 1;
-         restriction[1](9,8) = 1;
-         restriction[1](10,11) = 1;
-         restriction[1](26,27) = 1;
-         restriction[1](33,34) = 1;
-         restriction[1](41,42) = 1;
-         restriction[1](44,47) = 1;
-         restriction[1](57,62) = 1;
-         restriction[2](2,2) = 1;
-         restriction[2](11,10) = 1;
-         restriction[2](13,12) = 1;
-         restriction[2](28,29) = 1;
-         restriction[2](35,32) = 1;
-         restriction[2](46,45) = 1;
-         restriction[2](49,50) = 1;
-         restriction[2](61,58) = 1;
-         restriction[3](3,3) = 1;
-         restriction[3](12,13) = 1;
-         restriction[3](15,14) = 1;
-         restriction[3](30,31) = 1;
-         restriction[3](34,33) = 1;
-         restriction[3](48,51) = 1;
-         restriction[3](54,53) = 1;
-         restriction[3](60,59) = 1;
-         restriction[4](4,4) = 1;
-         restriction[4](16,17) = 1;
-         restriction[4](22,23) = 1;
-         restriction[4](25,24) = 1;
-         restriction[4](36,39) = 1;
-         restriction[4](42,41) = 1;
-         restriction[4](53,54) = 1;
-         restriction[4](58,61) = 1;
-         restriction[5](5,5) = 1;
-         restriction[5](17,16) = 1;
-         restriction[5](18,19) = 1;
-         restriction[5](27,26) = 1;
-         restriction[5](37,38) = 1;
-         restriction[5](43,40) = 1;
-         restriction[5](45,46) = 1;
-         restriction[5](59,60) = 1;
-         restriction[6](6,6) = 1;
-         restriction[6](19,18) = 1;
-         restriction[6](21,20) = 1;
-         restriction[6](29,28) = 1;
-         restriction[6](39,36) = 1;
-         restriction[6](47,44) = 1;
-         restriction[6](51,48) = 1;
-         restriction[6](63,56) = 1;
-         restriction[7](7,7) = 1;
-         restriction[7](20,21) = 1;
-         restriction[7](23,22) = 1;
-         restriction[7](31,30) = 1;
-         restriction[7](38,37) = 1;
-         restriction[7](50,49) = 1;
-         restriction[7](55,52) = 1;
-         restriction[7](62,57) = 1;
-         break;
-       case 4:
-         restriction[0](0,0) = 1;
-         restriction[0](8,9) = 1;
-         restriction[0](9,1) = 1;
-         restriction[0](17,18) = 1;
-         restriction[0](18,3) = 1;
-         restriction[0](32,33) = 1;
-         restriction[0](33,4) = 1;
-         restriction[0](44,48) = 1;
-         restriction[0](45,12) = 1;
-         restriction[0](47,15) = 1;
-         restriction[0](48,2) = 1;
-         restriction[0](62,66) = 1;
-         restriction[0](63,36) = 1;
-         restriction[0](65,21) = 1;
-         restriction[0](66,5) = 1;
-         restriction[0](89,93) = 1;
-         restriction[0](90,30) = 1;
-         restriction[0](92,42) = 1;
-         restriction[0](93,7) = 1;
-         restriction[0](98,111) = 1;
-         restriction[0](99,75) = 1;
-         restriction[0](101,57) = 1;
-         restriction[0](102,24) = 1;
-         restriction[0](107,84) = 1;
-         restriction[0](108,39) = 1;
-         restriction[0](110,27) = 1;
-         restriction[0](111,6) = 1;
-         restriction[1](1,1) = 1;
-         restriction[1](9,0) = 1;
-         restriction[1](10,9) = 1;
-         restriction[1](11,12) = 1;
-         restriction[1](12,2) = 1;
-         restriction[1](35,36) = 1;
-         restriction[1](36,5) = 1;
-         restriction[1](45,18) = 1;
-         restriction[1](46,48) = 1;
-         restriction[1](48,3) = 1;
-         restriction[1](49,15) = 1;
-         restriction[1](63,33) = 1;
-         restriction[1](64,66) = 1;
-         restriction[1](66,4) = 1;
-         restriction[1](67,21) = 1;
-         restriction[1](71,75) = 1;
-         restriction[1](72,24) = 1;
-         restriction[1](74,39) = 1;
-         restriction[1](75,6) = 1;
-         restriction[1](99,93) = 1;
-         restriction[1](100,111) = 1;
-         restriction[1](102,30) = 1;
-         restriction[1](103,57) = 1;
-         restriction[1](108,42) = 1;
-         restriction[1](109,84) = 1;
-         restriction[1](111,7) = 1;
-         restriction[1](112,27) = 1;
-         restriction[2](2,2) = 1;
-         restriction[2](12,1) = 1;
-         restriction[2](13,12) = 1;
-         restriction[2](15,3) = 1;
-         restriction[2](16,15) = 1;
-         restriction[2](38,39) = 1;
-         restriction[2](39,6) = 1;
-         restriction[2](48,0) = 1;
-         restriction[2](49,9) = 1;
-         restriction[2](51,18) = 1;
-         restriction[2](52,48) = 1;
-         restriction[2](74,36) = 1;
-         restriction[2](75,5) = 1;
-         restriction[2](77,75) = 1;
-         restriction[2](78,24) = 1;
-         restriction[2](81,42) = 1;
-         restriction[2](82,84) = 1;
-         restriction[2](84,7) = 1;
-         restriction[2](85,27) = 1;
-         restriction[2](108,33) = 1;
-         restriction[2](109,66) = 1;
-         restriction[2](111,4) = 1;
-         restriction[2](112,21) = 1;
-         restriction[2](117,93) = 1;
-         restriction[2](118,111) = 1;
-         restriction[2](120,30) = 1;
-         restriction[2](121,57) = 1;
-         restriction[3](3,3) = 1;
-         restriction[3](14,15) = 1;
-         restriction[3](15,2) = 1;
-         restriction[3](18,0) = 1;
-         restriction[3](19,18) = 1;
-         restriction[3](41,42) = 1;
-         restriction[3](42,7) = 1;
-         restriction[3](47,9) = 1;
-         restriction[3](48,1) = 1;
-         restriction[3](50,48) = 1;
-         restriction[3](51,12) = 1;
-         restriction[3](80,84) = 1;
-         restriction[3](81,39) = 1;
-         restriction[3](83,27) = 1;
-         restriction[3](84,6) = 1;
-         restriction[3](92,33) = 1;
-         restriction[3](93,4) = 1;
-         restriction[3](95,93) = 1;
-         restriction[3](96,30) = 1;
-         restriction[3](107,66) = 1;
-         restriction[3](108,36) = 1;
-         restriction[3](110,21) = 1;
-         restriction[3](111,5) = 1;
-         restriction[3](116,111) = 1;
-         restriction[3](117,75) = 1;
-         restriction[3](119,57) = 1;
-         restriction[3](120,24) = 1;
-         restriction[4](4,4) = 1;
-         restriction[4](20,21) = 1;
-         restriction[4](21,5) = 1;
-         restriction[4](29,30) = 1;
-         restriction[4](30,7) = 1;
-         restriction[4](33,0) = 1;
-         restriction[4](34,33) = 1;
-         restriction[4](53,57) = 1;
-         restriction[4](54,24) = 1;
-         restriction[4](56,27) = 1;
-         restriction[4](57,6) = 1;
-         restriction[4](65,9) = 1;
-         restriction[4](66,1) = 1;
-         restriction[4](68,66) = 1;
-         restriction[4](69,36) = 1;
-         restriction[4](90,18) = 1;
-         restriction[4](91,93) = 1;
-         restriction[4](93,3) = 1;
-         restriction[4](94,42) = 1;
-         restriction[4](101,48) = 1;
-         restriction[4](102,12) = 1;
-         restriction[4](104,111) = 1;
-         restriction[4](105,75) = 1;
-         restriction[4](110,15) = 1;
-         restriction[4](111,2) = 1;
-         restriction[4](113,84) = 1;
-         restriction[4](114,39) = 1;
-         restriction[5](5,5) = 1;
-         restriction[5](21,4) = 1;
-         restriction[5](22,21) = 1;
-         restriction[5](23,24) = 1;
-         restriction[5](24,6) = 1;
-         restriction[5](36,1) = 1;
-         restriction[5](37,36) = 1;
-         restriction[5](54,30) = 1;
-         restriction[5](55,57) = 1;
-         restriction[5](57,7) = 1;
-         restriction[5](58,27) = 1;
-         restriction[5](66,0) = 1;
-         restriction[5](67,9) = 1;
-         restriction[5](69,33) = 1;
-         restriction[5](70,66) = 1;
-         restriction[5](72,12) = 1;
-         restriction[5](73,75) = 1;
-         restriction[5](75,2) = 1;
-         restriction[5](76,39) = 1;
-         restriction[5](102,18) = 1;
-         restriction[5](103,48) = 1;
-         restriction[5](105,93) = 1;
-         restriction[5](106,111) = 1;
-         restriction[5](111,3) = 1;
-         restriction[5](112,15) = 1;
-         restriction[5](114,42) = 1;
-         restriction[5](115,84) = 1;
-         restriction[6](6,6) = 1;
-         restriction[6](24,5) = 1;
-         restriction[6](25,24) = 1;
-         restriction[6](27,7) = 1;
-         restriction[6](28,27) = 1;
-         restriction[6](39,2) = 1;
-         restriction[6](40,39) = 1;
-         restriction[6](57,4) = 1;
-         restriction[6](58,21) = 1;
-         restriction[6](60,30) = 1;
-         restriction[6](61,57) = 1;
-         restriction[6](75,1) = 1;
-         restriction[6](76,36) = 1;
-         restriction[6](78,12) = 1;
-         restriction[6](79,75) = 1;
-         restriction[6](84,3) = 1;
-         restriction[6](85,15) = 1;
-         restriction[6](87,42) = 1;
-         restriction[6](88,84) = 1;
-         restriction[6](111,0) = 1;
-         restriction[6](112,9) = 1;
-         restriction[6](114,33) = 1;
-         restriction[6](115,66) = 1;
-         restriction[6](120,18) = 1;
-         restriction[6](121,48) = 1;
-         restriction[6](123,93) = 1;
-         restriction[6](124,111) = 1;
-         restriction[7](7,7) = 1;
-         restriction[7](26,27) = 1;
-         restriction[7](27,6) = 1;
-         restriction[7](30,4) = 1;
-         restriction[7](31,30) = 1;
-         restriction[7](42,3) = 1;
-         restriction[7](43,42) = 1;
-         restriction[7](56,21) = 1;
-         restriction[7](57,5) = 1;
-         restriction[7](59,57) = 1;
-         restriction[7](60,24) = 1;
-         restriction[7](83,15) = 1;
-         restriction[7](84,2) = 1;
-         restriction[7](86,84) = 1;
-         restriction[7](87,39) = 1;
-         restriction[7](93,0) = 1;
-         restriction[7](94,33) = 1;
-         restriction[7](96,18) = 1;
-         restriction[7](97,93) = 1;
-         restriction[7](110,9) = 1;
-         restriction[7](111,1) = 1;
-         restriction[7](113,66) = 1;
-         restriction[7](114,36) = 1;
-         restriction[7](119,48) = 1;
-         restriction[7](120,12) = 1;
-         restriction[7](122,111) = 1;
-         restriction[7](123,75) = 1;
-         break;
-       default:
-         for (unsigned int i=0; i<GeometryInfo<dim>::children_per_cell;++i)
-           restriction[i].reinit(0,0);
-       }
-      break;
-    default:
-      Assert (false,ExcNotImplemented());
+      case 1:
+           switch (degree)
+             {
+               case 1:
+                     restriction[0](0,0) = 1;
+                     restriction[1](1,1) = 1;
+                     break;
+               case 2:
+                     restriction[0](0,0) = 1;
+                     restriction[0](2,1) = 1;
+                     restriction[1](1,1) = 1;
+                     restriction[1](1,1) = 1;
+                     break;
+               case 3:
+                     restriction[0](0,0) = 1;
+                     restriction[0](2,3) = 1;
+                     restriction[1](1,1) = 1;
+                     restriction[1](3,2) = 1;
+                     break;
+               case 4:
+                     restriction[0](0,0) = 1;
+                     restriction[0](2,3) = 1;
+                     restriction[0](3,1) = 1;
+                     restriction[1](1,1) = 1;
+                     restriction[1](3,0) = 1;
+                     restriction[1](4,3) = 1;
+                     break;
+               default:
+                     for (unsigned int i=0; i<GeometryInfo<dim>::children_per_cell;++i)
+                       restriction[i].reinit(0,0);
+             }
+           break;
+      case 2:
+           switch (degree)
+             {
+               case 1:
+                     restriction[0](0,0) = 1;
+                     restriction[1](1,1) = 1;
+                     restriction[2](2,2) = 1;
+                     restriction[3](3,3) = 1;
+                     break;
+               case 2:
+                     restriction[0](0,0) = 1;
+                     restriction[0](4,1) = 1;
+                     restriction[0](7,3) = 1;
+                     restriction[0](8,2) = 1;
+                     restriction[1](1,1) = 1;
+                     restriction[1](4,0) = 1;
+                     restriction[1](5,2) = 1;
+                     restriction[1](8,3) = 1;
+                     restriction[2](2,2) = 1;
+                     restriction[2](5,1) = 1;
+                     restriction[2](6,3) = 1;
+                     restriction[2](8,0) = 1;
+                     restriction[3](3,3) = 1;
+                     restriction[3](6,2) = 1;
+                     restriction[3](7,0) = 1;
+                     restriction[3](8,1) = 1;
+                     break;
+               case 3:
+                     restriction[0](0,0) = 1;
+                     restriction[0](4,5) = 1;
+                     restriction[0](10,11) = 1;
+                     restriction[0](12,15) = 1;
+                     restriction[1](1,1) = 1;
+                     restriction[1](5,4) = 1;
+                     restriction[1](6,7) = 1;
+                     restriction[1](13,14) = 1;
+                     restriction[2](2,2) = 1;
+                     restriction[2](7,6) = 1;
+                     restriction[2](9,8) = 1;
+                     restriction[2](15,12) = 1;
+                     restriction[3](3,3) = 1;
+                     restriction[3](8,9) = 1;
+                     restriction[3](11,10) = 1;
+                     restriction[3](14,13) = 1;
+                     break;
+               case 4:
+                     restriction[0](0,0) = 1;
+                     restriction[0](4,5) = 1;
+                     restriction[0](5,1) = 1;
+                     restriction[0](13,14) = 1;
+                     restriction[0](14,3) = 1;
+                     restriction[0](16,20) = 1;
+                     restriction[0](17,8) = 1;
+                     restriction[0](19,11) = 1;
+                     restriction[0](20,2) = 1;
+                     restriction[1](1,1) = 1;
+                     restriction[1](5,0) = 1;
+                     restriction[1](6,5) = 1;
+                     restriction[1](7,8) = 1;
+                     restriction[1](8,2) = 1;
+                     restriction[1](17,14) = 1;
+                     restriction[1](18,20) = 1;
+                     restriction[1](20,3) = 1;
+                     restriction[1](21,11) = 1;
+                     restriction[2](2,2) = 1;
+                     restriction[2](8,1) = 1;
+                     restriction[2](9,8) = 1;
+                     restriction[2](11,3) = 1;
+                     restriction[2](12,11) = 1;
+                     restriction[2](20,0) = 1;
+                     restriction[2](21,5) = 1;
+                     restriction[2](23,14) = 1;
+                     restriction[2](24,20) = 1;
+                     restriction[3](3,3) = 1;
+                     restriction[3](10,11) = 1;
+                     restriction[3](11,2) = 1;
+                     restriction[3](14,0) = 1;
+                     restriction[3](15,14) = 1;
+                     restriction[3](19,5) = 1;
+                     restriction[3](20,1) = 1;
+                     restriction[3](22,20) = 1;
+                     restriction[3](23,8) = 1;
+                     break;
+               default:
+                     for (unsigned int i=0; i<GeometryInfo<dim>::children_per_cell;++i)
+                       restriction[i].reinit(0,0);
+             }
+           break;
+      case 3:
+           switch (degree)
+             {
+               case 1:
+                     restriction[0](0,0) = 1;
+                     restriction[1](1,1) = 1;
+                     restriction[2](2,2) = 1;
+                     restriction[3](3,3) = 1;
+                     restriction[4](4,4) = 1;
+                     restriction[5](5,5) = 1;
+                     restriction[6](6,6) = 1;
+                     restriction[7](7,7) = 1;
+                     break;
+               case 2:
+                     restriction[0](0,0) = 1;
+                     restriction[0](8,1) = 1;
+                     restriction[0](11,3) = 1;
+                     restriction[0](16,4) = 1;
+                     restriction[0](20,2) = 1;
+                     restriction[0](22,5) = 1;
+                     restriction[0](25,7) = 1;
+                     restriction[0](26,6) = 1;
+                     restriction[1](1,1) = 1;
+                     restriction[1](8,0) = 1;
+                     restriction[1](9,2) = 1;
+                     restriction[1](17,5) = 1;
+                     restriction[1](20,3) = 1;
+                     restriction[1](22,4) = 1;
+                     restriction[1](23,6) = 1;
+                     restriction[1](26,7) = 1;
+                     restriction[2](2,2) = 1;
+                     restriction[2](9,1) = 1;
+                     restriction[2](10,3) = 1;
+                     restriction[2](18,6) = 1;
+                     restriction[2](20,0) = 1;
+                     restriction[2](23,5) = 1;
+                     restriction[2](24,7) = 1;
+                     restriction[2](26,4) = 1;
+                     restriction[3](3,3) = 1;
+                     restriction[3](10,2) = 1;
+                     restriction[3](11,0) = 1;
+                     restriction[3](19,7) = 1;
+                     restriction[3](20,1) = 1;
+                     restriction[3](24,6) = 1;
+                     restriction[3](25,4) = 1;
+                     restriction[3](26,5) = 1;
+                     restriction[4](4,4) = 1;
+                     restriction[4](12,5) = 1;
+                     restriction[4](15,7) = 1;
+                     restriction[4](16,0) = 1;
+                     restriction[4](21,6) = 1;
+                     restriction[4](22,1) = 1;
+                     restriction[4](25,3) = 1;
+                     restriction[4](26,2) = 1;
+                     restriction[5](5,5) = 1;
+                     restriction[5](12,4) = 1;
+                     restriction[5](13,6) = 1;
+                     restriction[5](17,1) = 1;
+                     restriction[5](21,7) = 1;
+                     restriction[5](22,0) = 1;
+                     restriction[5](23,2) = 1;
+                     restriction[5](26,3) = 1;
+                     restriction[6](6,6) = 1;
+                     restriction[6](13,5) = 1;
+                     restriction[6](14,7) = 1;
+                     restriction[6](18,2) = 1;
+                     restriction[6](21,4) = 1;
+                     restriction[6](23,1) = 1;
+                     restriction[6](24,3) = 1;
+                     restriction[6](26,0) = 1;
+                     restriction[7](7,7) = 1;
+                     restriction[7](14,6) = 1;
+                     restriction[7](15,4) = 1;
+                     restriction[7](19,3) = 1;
+                     restriction[7](21,5) = 1;
+                     restriction[7](24,2) = 1;
+                     restriction[7](25,0) = 1;
+                     restriction[7](26,1) = 1;
+                     break;
+               case 3:
+                     restriction[0](0,0) = 1;
+                     restriction[0](8,9) = 1;
+                     restriction[0](14,15) = 1;
+                     restriction[0](24,25) = 1;
+                     restriction[0](32,35) = 1;
+                     restriction[0](40,43) = 1;
+                     restriction[0](52,55) = 1;
+                     restriction[0](56,63) = 1;
+                     restriction[1](1,1) = 1;
+                     restriction[1](9,8) = 1;
+                     restriction[1](10,11) = 1;
+                     restriction[1](26,27) = 1;
+                     restriction[1](33,34) = 1;
+                     restriction[1](41,42) = 1;
+                     restriction[1](44,47) = 1;
+                     restriction[1](57,62) = 1;
+                     restriction[2](2,2) = 1;
+                     restriction[2](11,10) = 1;
+                     restriction[2](13,12) = 1;
+                     restriction[2](28,29) = 1;
+                     restriction[2](35,32) = 1;
+                     restriction[2](46,45) = 1;
+                     restriction[2](49,50) = 1;
+                     restriction[2](61,58) = 1;
+                     restriction[3](3,3) = 1;
+                     restriction[3](12,13) = 1;
+                     restriction[3](15,14) = 1;
+                     restriction[3](30,31) = 1;
+                     restriction[3](34,33) = 1;
+                     restriction[3](48,51) = 1;
+                     restriction[3](54,53) = 1;
+                     restriction[3](60,59) = 1;
+                     restriction[4](4,4) = 1;
+                     restriction[4](16,17) = 1;
+                     restriction[4](22,23) = 1;
+                     restriction[4](25,24) = 1;
+                     restriction[4](36,39) = 1;
+                     restriction[4](42,41) = 1;
+                     restriction[4](53,54) = 1;
+                     restriction[4](58,61) = 1;
+                     restriction[5](5,5) = 1;
+                     restriction[5](17,16) = 1;
+                     restriction[5](18,19) = 1;
+                     restriction[5](27,26) = 1;
+                     restriction[5](37,38) = 1;
+                     restriction[5](43,40) = 1;
+                     restriction[5](45,46) = 1;
+                     restriction[5](59,60) = 1;
+                     restriction[6](6,6) = 1;
+                     restriction[6](19,18) = 1;
+                     restriction[6](21,20) = 1;
+                     restriction[6](29,28) = 1;
+                     restriction[6](39,36) = 1;
+                     restriction[6](47,44) = 1;
+                     restriction[6](51,48) = 1;
+                     restriction[6](63,56) = 1;
+                     restriction[7](7,7) = 1;
+                     restriction[7](20,21) = 1;
+                     restriction[7](23,22) = 1;
+                     restriction[7](31,30) = 1;
+                     restriction[7](38,37) = 1;
+                     restriction[7](50,49) = 1;
+                     restriction[7](55,52) = 1;
+                     restriction[7](62,57) = 1;
+                     break;
+               case 4:
+                     restriction[0](0,0) = 1;
+                     restriction[0](8,9) = 1;
+                     restriction[0](9,1) = 1;
+                     restriction[0](17,18) = 1;
+                     restriction[0](18,3) = 1;
+                     restriction[0](32,33) = 1;
+                     restriction[0](33,4) = 1;
+                     restriction[0](44,48) = 1;
+                     restriction[0](45,12) = 1;
+                     restriction[0](47,15) = 1;
+                     restriction[0](48,2) = 1;
+                     restriction[0](62,66) = 1;
+                     restriction[0](63,36) = 1;
+                     restriction[0](65,21) = 1;
+                     restriction[0](66,5) = 1;
+                     restriction[0](89,93) = 1;
+                     restriction[0](90,30) = 1;
+                     restriction[0](92,42) = 1;
+                     restriction[0](93,7) = 1;
+                     restriction[0](98,111) = 1;
+                     restriction[0](99,75) = 1;
+                     restriction[0](101,57) = 1;
+                     restriction[0](102,24) = 1;
+                     restriction[0](107,84) = 1;
+                     restriction[0](108,39) = 1;
+                     restriction[0](110,27) = 1;
+                     restriction[0](111,6) = 1;
+                     restriction[1](1,1) = 1;
+                     restriction[1](9,0) = 1;
+                     restriction[1](10,9) = 1;
+                     restriction[1](11,12) = 1;
+                     restriction[1](12,2) = 1;
+                     restriction[1](35,36) = 1;
+                     restriction[1](36,5) = 1;
+                     restriction[1](45,18) = 1;
+                     restriction[1](46,48) = 1;
+                     restriction[1](48,3) = 1;
+                     restriction[1](49,15) = 1;
+                     restriction[1](63,33) = 1;
+                     restriction[1](64,66) = 1;
+                     restriction[1](66,4) = 1;
+                     restriction[1](67,21) = 1;
+                     restriction[1](71,75) = 1;
+                     restriction[1](72,24) = 1;
+                     restriction[1](74,39) = 1;
+                     restriction[1](75,6) = 1;
+                     restriction[1](99,93) = 1;
+                     restriction[1](100,111) = 1;
+                     restriction[1](102,30) = 1;
+                     restriction[1](103,57) = 1;
+                     restriction[1](108,42) = 1;
+                     restriction[1](109,84) = 1;
+                     restriction[1](111,7) = 1;
+                     restriction[1](112,27) = 1;
+                     restriction[2](2,2) = 1;
+                     restriction[2](12,1) = 1;
+                     restriction[2](13,12) = 1;
+                     restriction[2](15,3) = 1;
+                     restriction[2](16,15) = 1;
+                     restriction[2](38,39) = 1;
+                     restriction[2](39,6) = 1;
+                     restriction[2](48,0) = 1;
+                     restriction[2](49,9) = 1;
+                     restriction[2](51,18) = 1;
+                     restriction[2](52,48) = 1;
+                     restriction[2](74,36) = 1;
+                     restriction[2](75,5) = 1;
+                     restriction[2](77,75) = 1;
+                     restriction[2](78,24) = 1;
+                     restriction[2](81,42) = 1;
+                     restriction[2](82,84) = 1;
+                     restriction[2](84,7) = 1;
+                     restriction[2](85,27) = 1;
+                     restriction[2](108,33) = 1;
+                     restriction[2](109,66) = 1;
+                     restriction[2](111,4) = 1;
+                     restriction[2](112,21) = 1;
+                     restriction[2](117,93) = 1;
+                     restriction[2](118,111) = 1;
+                     restriction[2](120,30) = 1;
+                     restriction[2](121,57) = 1;
+                     restriction[3](3,3) = 1;
+                     restriction[3](14,15) = 1;
+                     restriction[3](15,2) = 1;
+                     restriction[3](18,0) = 1;
+                     restriction[3](19,18) = 1;
+                     restriction[3](41,42) = 1;
+                     restriction[3](42,7) = 1;
+                     restriction[3](47,9) = 1;
+                     restriction[3](48,1) = 1;
+                     restriction[3](50,48) = 1;
+                     restriction[3](51,12) = 1;
+                     restriction[3](80,84) = 1;
+                     restriction[3](81,39) = 1;
+                     restriction[3](83,27) = 1;
+                     restriction[3](84,6) = 1;
+                     restriction[3](92,33) = 1;
+                     restriction[3](93,4) = 1;
+                     restriction[3](95,93) = 1;
+                     restriction[3](96,30) = 1;
+                     restriction[3](107,66) = 1;
+                     restriction[3](108,36) = 1;
+                     restriction[3](110,21) = 1;
+                     restriction[3](111,5) = 1;
+                     restriction[3](116,111) = 1;
+                     restriction[3](117,75) = 1;
+                     restriction[3](119,57) = 1;
+                     restriction[3](120,24) = 1;
+                     restriction[4](4,4) = 1;
+                     restriction[4](20,21) = 1;
+                     restriction[4](21,5) = 1;
+                     restriction[4](29,30) = 1;
+                     restriction[4](30,7) = 1;
+                     restriction[4](33,0) = 1;
+                     restriction[4](34,33) = 1;
+                     restriction[4](53,57) = 1;
+                     restriction[4](54,24) = 1;
+                     restriction[4](56,27) = 1;
+                     restriction[4](57,6) = 1;
+                     restriction[4](65,9) = 1;
+                     restriction[4](66,1) = 1;
+                     restriction[4](68,66) = 1;
+                     restriction[4](69,36) = 1;
+                     restriction[4](90,18) = 1;
+                     restriction[4](91,93) = 1;
+                     restriction[4](93,3) = 1;
+                     restriction[4](94,42) = 1;
+                     restriction[4](101,48) = 1;
+                     restriction[4](102,12) = 1;
+                     restriction[4](104,111) = 1;
+                     restriction[4](105,75) = 1;
+                     restriction[4](110,15) = 1;
+                     restriction[4](111,2) = 1;
+                     restriction[4](113,84) = 1;
+                     restriction[4](114,39) = 1;
+                     restriction[5](5,5) = 1;
+                     restriction[5](21,4) = 1;
+                     restriction[5](22,21) = 1;
+                     restriction[5](23,24) = 1;
+                     restriction[5](24,6) = 1;
+                     restriction[5](36,1) = 1;
+                     restriction[5](37,36) = 1;
+                     restriction[5](54,30) = 1;
+                     restriction[5](55,57) = 1;
+                     restriction[5](57,7) = 1;
+                     restriction[5](58,27) = 1;
+                     restriction[5](66,0) = 1;
+                     restriction[5](67,9) = 1;
+                     restriction[5](69,33) = 1;
+                     restriction[5](70,66) = 1;
+                     restriction[5](72,12) = 1;
+                     restriction[5](73,75) = 1;
+                     restriction[5](75,2) = 1;
+                     restriction[5](76,39) = 1;
+                     restriction[5](102,18) = 1;
+                     restriction[5](103,48) = 1;
+                     restriction[5](105,93) = 1;
+                     restriction[5](106,111) = 1;
+                     restriction[5](111,3) = 1;
+                     restriction[5](112,15) = 1;
+                     restriction[5](114,42) = 1;
+                     restriction[5](115,84) = 1;
+                     restriction[6](6,6) = 1;
+                     restriction[6](24,5) = 1;
+                     restriction[6](25,24) = 1;
+                     restriction[6](27,7) = 1;
+                     restriction[6](28,27) = 1;
+                     restriction[6](39,2) = 1;
+                     restriction[6](40,39) = 1;
+                     restriction[6](57,4) = 1;
+                     restriction[6](58,21) = 1;
+                     restriction[6](60,30) = 1;
+                     restriction[6](61,57) = 1;
+                     restriction[6](75,1) = 1;
+                     restriction[6](76,36) = 1;
+                     restriction[6](78,12) = 1;
+                     restriction[6](79,75) = 1;
+                     restriction[6](84,3) = 1;
+                     restriction[6](85,15) = 1;
+                     restriction[6](87,42) = 1;
+                     restriction[6](88,84) = 1;
+                     restriction[6](111,0) = 1;
+                     restriction[6](112,9) = 1;
+                     restriction[6](114,33) = 1;
+                     restriction[6](115,66) = 1;
+                     restriction[6](120,18) = 1;
+                     restriction[6](121,48) = 1;
+                     restriction[6](123,93) = 1;
+                     restriction[6](124,111) = 1;
+                     restriction[7](7,7) = 1;
+                     restriction[7](26,27) = 1;
+                     restriction[7](27,6) = 1;
+                     restriction[7](30,4) = 1;
+                     restriction[7](31,30) = 1;
+                     restriction[7](42,3) = 1;
+                     restriction[7](43,42) = 1;
+                     restriction[7](56,21) = 1;
+                     restriction[7](57,5) = 1;
+                     restriction[7](59,57) = 1;
+                     restriction[7](60,24) = 1;
+                     restriction[7](83,15) = 1;
+                     restriction[7](84,2) = 1;
+                     restriction[7](86,84) = 1;
+                     restriction[7](87,39) = 1;
+                     restriction[7](93,0) = 1;
+                     restriction[7](94,33) = 1;
+                     restriction[7](96,18) = 1;
+                     restriction[7](97,93) = 1;
+                     restriction[7](110,9) = 1;
+                     restriction[7](111,1) = 1;
+                     restriction[7](113,66) = 1;
+                     restriction[7](114,36) = 1;
+                     restriction[7](119,48) = 1;
+                     restriction[7](120,12) = 1;
+                     restriction[7](122,111) = 1;
+                     restriction[7](123,75) = 1;
+                     break;
+               default:
+                     for (unsigned int i=0; i<GeometryInfo<dim>::children_per_cell;++i)
+                       restriction[i].reinit(0,0);
+             }
+           break;
+      default:
+           Assert (false,ExcNotImplemented());
     }
 }
 
@@ -596,6 +597,7 @@ FE_Q<dim>::~FE_Q ()
 }
 
 
+
 template <int dim>
 FiniteElement<dim> *
 FE_Q<dim>::clone() const
@@ -604,19 +606,21 @@ FE_Q<dim>::clone() const
 }
 
 
+
 template <int dim>
 void
 FE_Q<dim>::get_unit_support_points (std::vector<Point<dim> > &points) const
 {
-  compute_support_points (points, degree, renumber);
+  compute_support_points (degree, renumber, points);
 }
 
+
     
 template <int dim>
 void
 FE_Q<dim>::get_unit_face_support_points (std::vector<Point<dim-1> > &points) const
 {
-  FE_Q<dim-1>::compute_support_points (points, degree, face_renumber);
+  FE_Q<dim-1>::compute_support_points (degree, face_renumber, points);
 }
 
 
@@ -630,31 +634,34 @@ FE_Q<1>::get_unit_face_support_points (std::vector<Point<0> > &points) const
 }
 
 #endif
+
+
+
 //----------------------------------------------------------------------
-// Auxilliary functions
+// Auxiliary functions
 //----------------------------------------------------------------------
 
 
 template <int dim>
 void
-FE_Q<dim>::compute_support_points (std::vector<Point<dim> >& support_points,
-                                  unsigned int degree,
-                                  const std::vector<unsigned int>& renumber)
+FE_Q<dim>::compute_support_points (const unsigned int               degree,
+                                  const std::vector<unsigned int> &renumber,
+                                  std::vector<Point<dim> >        &support_points)
 {
                                   // number of points: (degree+1)^dim
-  unsigned int n= degree+1;
-  for (unsigned int i=1;i<dim;++i)
+  unsigned int n = degree+1;
+  for (unsigned int i=1; i<dim; ++i)
     n *= degree+1;
 
   support_points.resize(n);
   
-  double step = 1./degree;
+  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)
+  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)
@@ -669,7 +676,7 @@ FE_Q<dim>::compute_support_points (std::vector<Point<dim> >& support_points,
 
 template <int dim>
 std::vector<unsigned int>
-FE_Q<dim>::get_dpo_vector(unsigned int deg)
+FE_Q<dim>::get_dpo_vector(const unsigned int deg)
 {
   std::vector<unsigned int> dpo(dim+1, 1);
   for (unsigned int i=1; i<dpo.size(); ++i)
@@ -682,8 +689,8 @@ FE_Q<dim>::get_dpo_vector(unsigned int deg)
 template <int dim>
 void
 FE_Q<dim>::build_renumbering (const FiniteElementData<dim> &fe_data,
-                             unsigned int degree,
-                             std::vector<unsigned int> &renumber)
+                             const unsigned int            degree,
+                             std::vector<unsigned int>    &renumber)
 {
   const unsigned int n = degree+1;
   
@@ -702,58 +709,59 @@ FE_Q<dim>::build_renumbering (const FiniteElementData<dim> &fe_data,
                                         // sophisticated.
        switch (dim)
          {
-         case 1:
-           if (i==1)
-             index += degree;
-           break;
-         case 2:
-           switch (i)
-             {
-             case 1:
-               index += degree;
-               break;
-             case 3:
-               index += n*degree;
-               break;
-             case 2:
-               index += n*degree+degree;
-               break;
-             }
-           break;
-         case 3:
-           switch (i)
-             {
-             case 1:
-               index += degree;
-               break;
-             case 4:
-               index += n*degree;
-               break;
-             case 5:
-               index += n*degree+degree;
-               break;
-             case 3:
-               index += n*n*degree;
-               break;
-             case 2:
-               index += n*n*degree + degree;
-               break;
-             case 7:
-               index += n*n*degree + n*degree;
-               break;
-             case 6:
-               index += n*n*degree + n*degree+degree;
-               break;
-             }
-           break;
-           
-         default:
-           Assert(false, ExcNotImplemented());
+           case 1:
+                 if (i==1)
+                   index += degree;
+                 break;
+           case 2:
+                 switch (i)
+                   {
+                     case 1:
+                           index += degree;
+                           break;
+                     case 3:
+                           index += n*degree;
+                           break;
+                     case 2:
+                           index += n*degree+degree;
+                           break;
+                   }
+                 break;
+           case 3:
+                 switch (i)
+                   {
+                     case 1:
+                           index += degree;
+                           break;
+                     case 4:
+                           index += n*degree;
+                           break;
+                     case 5:
+                           index += n*degree+degree;
+                           break;
+                     case 3:
+                           index += n*n*degree;
+                           break;
+                     case 2:
+                           index += n*n*degree + degree;
+                           break;
+                     case 7:
+                           index += n*n*degree + n*degree;
+                           break;
+                     case 6:
+                           index += n*n*degree + n*degree+degree;
+                           break;
+                   }
+                 break;
+                 
+           default:
+                 Assert(false, ExcNotImplemented());
          }
        
        renumber[index] = i;
       }
   else
+                                    // degree == 0
     renumber[0] = 0;
   
                                   // Lines and higher
@@ -769,62 +777,62 @@ FE_Q<dim>::build_renumbering (const FiniteElementData<dim> &fe_data,
                                           // numbering scheme.
          switch (i+100*dim)
            {
-             // lines in x-direction
-           case 100:
-           case 200: case 202:
-           case 300: case 302: case 304: case 306:
-             incr = 1;
-             break;
-             // lines in y-direction
-           case 201: case 203:
-           case 308: case 309: case 310: case 311:
-             incr = n;
-             break;
-             // lines in z-direction
-           case 301: case 303: case 305: case 307:
-             incr = n*n;
-             break;
-           default:
-             Assert(false, ExcNotImplemented());
+                                              // lines in x-direction
+             case 100:
+             case 200: case 202:
+             case 300: case 302: case 304: case 306:
+                   incr = 1;
+                   break;
+                                                    // lines in y-direction
+             case 201: case 203:
+             case 308: case 309: case 310: case 311:
+                   incr = n;
+                   break;
+                                                    // lines in z-direction
+             case 301: case 303: case 305: case 307:
+                   incr = n*n;
+                   break;
+             default:
+                   Assert(false, ExcNotImplemented());
            }
          switch (i+100*dim)
            {
-             // x=y=z=0
-           case 100:
-           case 200: case 203:
-           case 300: case 303: case 308:
-             tensorstart = 0;
-             break;
-             // x=1 y=z=0
-           case 201:
-           case 301: case 309:
-             tensorstart = degree;
-             break;
-             // y=1 x=z=0
-           case 202:
-           case 304: case 307:
-             tensorstart = n*degree;
-             break;
-             // x=z=1 y=0
-           case 310:
-             tensorstart = n*n*degree+degree;
-             break;
-             // z=1 x=y=0
-           case 302: case 311:
-             tensorstart = n*n*degree;
-             break;
-             // x=y=1 z=0
-           case 305:
-             tensorstart = n*degree+degree;
-             break;
-             // y=z=1 x=0
-           case 306:
-             tensorstart = n*n*n-n;
-             break;
-           default:
-             Assert(false, ExcNotImplemented());             
+                                              // x=y=z=0
+             case 100:
+             case 200: case 203:
+             case 300: case 303: case 308:
+                   tensorstart = 0;
+                   break;
+                                                    // x=1 y=z=0
+             case 201:
+             case 301: case 309:
+                   tensorstart = degree;
+                   break;
+                                                    // y=1 x=z=0
+             case 202:
+             case 304: case 307:
+                   tensorstart = n*degree;
+                   break;
+                                                    // x=z=1 y=0
+             case 310:
+                   tensorstart = n*n*degree+degree;
+                   break;
+                                                    // z=1 x=y=0
+             case 302: case 311:
+                   tensorstart = n*n*degree;
+                   break;
+                                                    // x=y=1 z=0
+             case 305:
+                   tensorstart = n*degree+degree;
+                   break;
+                                                    // y=z=1 x=0
+             case 306:
+                   tensorstart = n*n*n-n;
+                   break;
+             default:
+                   Assert(false, ExcNotImplemented());       
            }
-
+         
          for (unsigned int jx = 1; jx<degree ;++jx)
            {
              unsigned int tensorindex = tensorstart + jx * incr;
@@ -841,29 +849,29 @@ FE_Q<dim>::build_renumbering (const FiniteElementData<dim> &fe_data,
          switch (i)
            {
              case 0:
-               tensorstart = 0; incx = 1;
-               if (dim==2)
-                 incy = n;
-               else
-                 incy = n*n;
-               break;
+                   tensorstart = 0; incx = 1;
+                   if (dim==2)
+                     incy = n;
+                   else
+                     incy = n*n;
+                   break;
              case 1:
-               tensorstart = n*degree; incx = 1; incy = n*n;
-               break;
+                   tensorstart = n*degree; incx = 1; incy = n*n;
+                   break;
              case 2:
-               tensorstart = 0; incx = 1; incy = n;
-               break;
+                   tensorstart = 0; incx = 1; incy = n;
+                   break;
              case 3:
-               tensorstart = degree; incx = n; incy = n*n;
-               break;
+                   tensorstart = degree; incx = n; incy = n*n;
+                   break;
              case 4:
-               tensorstart = n*n*degree; incx = 1; incy = n;
-               break;
+                   tensorstart = n*n*degree; incx = 1; incy = n;
+                   break;
              case 5:
-               tensorstart = 0; incx = n; incy = n*n;
-               break;
-           default:
-             Assert(false, ExcNotImplemented());             
+                   tensorstart = 0; incx = n; incy = n*n;
+                   break;
+             default:
+                   Assert(false, ExcNotImplemented());       
            }
          
          for (unsigned int jy = 1; jy<degree; jy++)
@@ -883,8 +891,7 @@ FE_Q<dim>::build_renumbering (const FiniteElementData<dim> &fe_data,
            for (unsigned int jy = 1; jy<degree; jy++)
              for (unsigned int jx = 1; jx<degree; jx++)
                {
-                 unsigned int tensorindex =
-                   jx + jy * n + jz * n*n;
+                 const unsigned int tensorindex = jx + jy*n + jz*n*n;
                  renumber[tensorindex]=index++;
                }  
        }
@@ -893,29 +900,32 @@ FE_Q<dim>::build_renumbering (const FiniteElementData<dim> &fe_data,
 }
 
 
+
 template <int dim>
 void
 FE_Q<dim>::build_face_renumbering (const FiniteElementData<dim-1> &fe_data,
-                                  unsigned int degree,
-                                  std::vector<unsigned int>numbering)
+                                  const unsigned int              degree,
+                                  std::vector<unsigned int>      &numbering)
 {
   FE_Q<dim-1>::build_renumbering (fe_data, degree, numbering); 
 }
 
+
 #if (deal_II_dimension == 1)
 
 template <>
 void
 FE_Q<1>::build_face_renumbering (const FiniteElementData<0> &,
-                                unsigned int,
+                                const unsigned int,
                                 std::vector<unsigned int>&)
 {}
+
 #endif
 
 
 template <int dim>
 UpdateFlags
-FE_Q<dim>::update_once (UpdateFlags flags) const
+FE_Q<dim>::update_once (const UpdateFlags flags) const
 {
   UpdateFlags out = update_default;
 
@@ -926,6 +936,7 @@ FE_Q<dim>::update_once (UpdateFlags flags) const
 }
 
 
+
 template <int dim>
 UpdateFlags
 FE_Q<dim>::update_each (UpdateFlags flags) const
@@ -941,13 +952,15 @@ FE_Q<dim>::update_each (UpdateFlags flags) const
 }
 
 
+
 //----------------------------------------------------------------------
 // Data field initialization
 //----------------------------------------------------------------------
+
 template <int dim>
 Mapping<dim>::InternalDataBase*
-FE_Q<dim>::get_data (const UpdateFlags update_flags,
-                    const Mapping<dim>mapping,
+FE_Q<dim>::get_data (const UpdateFlags      update_flags,
+                    const Mapping<dim>    &mapping,
                     const Quadrature<dim> &quadrature) const
 {
   InternalData* data = new InternalData;
@@ -960,29 +973,30 @@ FE_Q<dim>::get_data (const UpdateFlags update_flags,
   data->update_flags = data->update_once | data->update_each;
 
   const UpdateFlags flags(data->update_flags);
+  const unsigned int n_q_points = quadrature.n_quadrature_points;
   
   if (flags & update_values)
     {
       values.resize (dofs_per_cell);
       data->shape_values.resize(dofs_per_cell,
-                               std::vector<double>(quadrature.n_quadrature_points));
+                               std::vector<double>(n_q_points));
     }
 
   if (flags & update_gradients)
     {
       grads.resize (dofs_per_cell);
       data->shape_gradients.resize(dofs_per_cell,
-                                  std::vector<Tensor<1,dim> >(quadrature.n_quadrature_points));
+                                  std::vector<Tensor<1,dim> >(n_q_points));
     }
 
   if (flags & update_second_derivatives)
     data->initialize (this, mapping, quadrature);
   
   if (flags & (update_values | update_gradients))
-    for (unsigned int i=0;i<quadrature.n_quadrature_points;++i)
+    for (unsigned int i=0; i<n_q_points; ++i)
       {
        poly->compute(quadrature.point(i), values, grads, grad_grads);
-       for (unsigned int k=0;k<dofs_per_cell; ++k)
+       for (unsigned int k=0; k<dofs_per_cell; ++k)
          {
            if (flags & update_values)
              data->shape_values[renumber[k]][i] = values[k];
@@ -995,36 +1009,40 @@ FE_Q<dim>::get_data (const UpdateFlags update_flags,
 
 
 
+
 //----------------------------------------------------------------------
 // Fill data of FEValues
 //----------------------------------------------------------------------
+
 template <int dim>
 void
 FE_Q<dim>::fill_fe_values (const Mapping<dim>                   &mapping,
                           const DoFHandler<dim>::cell_iterator &cell,
                           const Quadrature<dim>                &quadrature,
-                          Mapping<dim>::InternalDataBase      &mapping_data,
-                          Mapping<dim>::InternalDataBase      &fedata,
+                          Mapping<dim>::InternalDataBase       &mapping_data,
+                          Mapping<dim>::InternalDataBase       &fedata,
                           FEValuesData<dim>                    &data) const
 {
+                                  // convert data object to internal
+                                  // data for this class. fails with
+                                  // an exception if that is not
+                                  // possible
   InternalData &fe_data = dynamic_cast<InternalData &> (fedata);
   
   const UpdateFlags flags(fe_data.current_update_flags());
 
-  for (unsigned int k=0;k<dofs_per_cell;++k)
+  for (unsigned int k=0; k<dofs_per_cell; ++k)
     {
-      for (unsigned int i=0;i<quadrature.n_quadrature_points;++i)
-       {
-         if (flags & update_values)
-           data.shape_values(k,i)
-             = fe_data.shape_values[k][i];
-       }
+      for (unsigned int i=0; i<quadrature.n_quadrature_points; ++i)
+       if (flags & update_values)
+         data.shape_values(k,i) = fe_data.shape_values[k][i];
       
       if (flags & update_gradients)
        mapping.transform_covariant(data.shape_gradients[k],
                                    fe_data.shape_gradients[k],
                                    mapping_data, 0);
     }
+
   if (flags & update_second_derivatives)
     compute_2nd (mapping, cell, 0, mapping_data, fe_data, data);
   
@@ -1035,35 +1053,41 @@ FE_Q<dim>::fill_fe_values (const Mapping<dim>                   &mapping,
 
 template <int dim>
 void
-FE_Q<dim>::fill_fe_face_values (const Mapping<dim> &mapping,
+FE_Q<dim>::fill_fe_face_values (const Mapping<dim>                   &mapping,
                                const DoFHandler<dim>::cell_iterator &cell,
-                               const unsigned int face,
-                               const Quadrature<dim-1>quadrature,
-                               Mapping<dim>::InternalDataBasemapping_data,
-                               Mapping<dim>::InternalDataBasefedata,
-                               FEValuesData<dim>data) const
+                               const unsigned int                    face,
+                               const Quadrature<dim-1>              &quadrature,
+                               Mapping<dim>::InternalDataBase       &mapping_data,
+                               Mapping<dim>::InternalDataBase       &fedata,
+                               FEValuesData<dim>                    &data) const
 {
+                                  // convert data object to internal
+                                  // data for this class. fails with
+                                  // an exception if that is not
+                                  // possible
   InternalData &fe_data = dynamic_cast<InternalData &> (fedata);
-  unsigned int offset = face * quadrature.n_quadrature_points;
+
+                                  // offset determines which data set
+                                  // to take (all data sets for all
+                                  // faces are stored contiguously)
+  const unsigned int offset = face * quadrature.n_quadrature_points;
   
   const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
 
-  for (unsigned int k=0;k<dofs_per_cell;++k)
+  for (unsigned int k=0; k<dofs_per_cell; ++k)
     {
       for (unsigned int i=0;i<quadrature.n_quadrature_points;++i)
-       {
-         if (flags & update_values)
-           data.shape_values(k,i)
-             = fe_data.shape_values[k][i+offset];
-       }
+       if (flags & update_values)
+         data.shape_values(k,i) = fe_data.shape_values[k][i+offset];
       
       if (flags & update_gradients)
        mapping.transform_covariant(data.shape_gradients[k],
                                    fe_data.shape_gradients[k],
                                    mapping_data, offset);
-      if (flags & update_second_derivatives)
-       compute_2nd (mapping, cell, offset, mapping_data, fe_data, data);
     }
+
+  if (flags & update_second_derivatives)
+    compute_2nd (mapping, cell, offset, mapping_data, fe_data, data);
   
   fe_data.first_cell = false;
 }
@@ -1072,38 +1096,44 @@ FE_Q<dim>::fill_fe_face_values (const Mapping<dim> &mapping,
 
 template <int dim>
 void
-FE_Q<dim>::fill_fe_subface_values (const Mapping<dim> &mapping,
+FE_Q<dim>::fill_fe_subface_values (const Mapping<dim>                   &mapping,
                                   const DoFHandler<dim>::cell_iterator &cell,
-                                  const unsigned int face,
-                                  const unsigned int subface,
-                                  const Quadrature<dim-1>quadrature,
-                                  Mapping<dim>::InternalDataBasemapping_data,
-                                  Mapping<dim>::InternalDataBasefedata,
-                                  FEValuesData<dim>data) const
+                                  const unsigned int                    face,
+                                  const unsigned int                    subface,
+                                  const Quadrature<dim-1>              &quadrature,
+                                  Mapping<dim>::InternalDataBase       &mapping_data,
+                                  Mapping<dim>::InternalDataBase       &fedata,
+                                  FEValuesData<dim>                    &data) const
 {
+                                  // convert data object to internal
+                                  // data for this class. fails with
+                                  // an exception if that is not
+                                  // possible
   InternalData &fe_data = dynamic_cast<InternalData &> (fedata);
-  unsigned int offset = (face * GeometryInfo<dim>::subfaces_per_face
-                        + subface) * quadrature.n_quadrature_points;
+
+                                  // offset determines which data set
+                                  // to take (all data sets for all
+                                  // sub-faces are stored contiguously)
+  const unsigned int offset = (face * GeometryInfo<dim>::subfaces_per_face + subface)
+                             * quadrature.n_quadrature_points;
 
   const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
 
-  for (unsigned int k=0;k<dofs_per_cell;++k)
+  for (unsigned int k=0; k<dofs_per_cell; ++k)
     {
       for (unsigned int i=0;i<quadrature.n_quadrature_points;++i)
-       {
-         if (flags & update_values)
-           data.shape_values(k,i)
-             = fe_data.shape_values[k][i+offset];
-       }
+       if (flags & update_values)
+         data.shape_values(k,i) = fe_data.shape_values[k][i+offset];
       
       if (flags & update_gradients)
        mapping.transform_covariant(data.shape_gradients[k],
                                    fe_data.shape_gradients[k],
                                    mapping_data, offset);
-      if (flags & update_second_derivatives)
-       compute_2nd (mapping, cell, offset, mapping_data, fe_data, data);
     }
   
+  if (flags & update_second_derivatives)
+    compute_2nd (mapping, cell, offset, mapping_data, fe_data, data);
+  
   fe_data.first_cell = false;
 }
 
index c6c848b5b541defa0b0fba4e4bcd3ec1f0339473..5f8ab47dcb4deab5ffaed660b270e969ecb8a724 100644 (file)
 
 // Transfer matrices for finite elements
 
-static double q1_into_q1_refined_0[] =
+static const double q1_into_q1_refined_0[] =
 {
   1., 0.,
   13.5/27., 13.5/27.,
 };
 
-static double q1_into_q1_refined_1[] =
+static const double q1_into_q1_refined_1[] =
 {
  13.5/27., 13.5/27.,
  0., 1.,
 };
 
-static double q2_into_q2_refined_0[] =
+static const double q2_into_q2_refined_0[] =
 {
   1., 0., 0.,
   0., 0., 1.,
  10.125/27., -3.375/27., 20.25/27.,
 };
 
-static double q2_into_q2_refined_1[] =
+static const double q2_into_q2_refined_1[] =
 {
  0., 0., 1.,
  0., 1., 0.,
  -3.375/27., 10.125/27., 20.25/27.,
 };
 
-static double q3_into_q3_refined_0[] =
+static const double q3_into_q3_refined_0[] =
 {
  1., 0., 0., 0.,
  -1.6875/27., -1.6875/27., 15.1875/27., 15.1875/27.,
@@ -47,7 +47,7 @@ static double q3_into_q3_refined_0[] =
  0., 0., 1., 0.
 };
 
-static double q3_into_q3_refined_1[] =
+static const double q3_into_q3_refined_1[] =
 {
  -1.6875/27., -1.6875/27., 15.1875/27., 15.1875/27.,
  0., 1., 0., 0.,
@@ -55,7 +55,7 @@ static double q3_into_q3_refined_1[] =
  1.6875/27., 8.4375/27., -8.4375/27., 25.3125/27.,
 };
 
-static double q4_into_q4_refined_0[] =
+static const double q4_into_q4_refined_0[] =
 {
   1., 0., 0., 0., 0.,
   0., 0., 0., 1., 0.,
@@ -64,7 +64,7 @@ static double q4_into_q4_refined_0[] =
  -1.0546875/27., 0.6328125/27., 12.65625/27., 18.984375/27., -4.21875/27.,
 };
 
-static double q4_into_q4_refined_1[] =
+static const double q4_into_q4_refined_1[] =
 {
   0., 0., 0., 1., 0.,
   0., 1., 0., 0., 0.,
@@ -73,42 +73,27 @@ static double q4_into_q4_refined_1[] =
  -1.0546875/27., 7.3828125/27., 5.90625/27., -14.765625/27., 29.53125/27.,
 };
 
-static double * embedding0[] =
-{
-  q1_into_q1_refined_0,
-  q2_into_q2_refined_0,
-  q3_into_q3_refined_0,
-  q4_into_q4_refined_0,
-  0,
-  0,
-  0,
-  0,
-  0,
-  0
-};
 
-static double * embedding1[] =
+
+// embedding matrices
+
+
+const unsigned int  FE_Q<1>::Matrices::n_embedding_matrices = 4;
+
+const double * const
+FE_Q<1>::Matrices::embedding[][GeometryInfo<1>::children_per_cell] =
 {
-  q1_into_q1_refined_1,
-  q2_into_q2_refined_1,
-  q3_into_q3_refined_1,
-  q4_into_q4_refined_1,
-  0,
-  0,
-  0,
-  0,
-  0,
-  0
+  {q1_into_q1_refined_0, q1_into_q1_refined_1},
+  {q2_into_q2_refined_0, q2_into_q2_refined_1},
+  {q3_into_q3_refined_0, q3_into_q3_refined_1},
+  {q4_into_q4_refined_0, q4_into_q4_refined_1},
 };
 
-static double** embedding2;
-static double** embedding3;
-static double** embedding4;
-static double** embedding5;
-static double** embedding6;
-static double** embedding7;
 
-// No constrants in 1d
+// No constraints in 1d
+const unsigned int 
+FE_Q<1>::Matrices::n_constraint_matrices = 0;
 
-static double ** constraint_matrices = 0;
+const double * const
+FE_Q<1>::Matrices::constraint_matrices[] = { 0 };
 
index 089a416c63919cf77a7229806663169ac4a75df9..c52d9cc7bd0fc2e1737bf45d38d4961a5ff0137d 100644 (file)
@@ -13,7 +13,7 @@
 
 // Transfer matrices for finite elements
 
-static double q1_into_q1_refined_0[] =
+static const double q1_into_q1_refined_0[] =
 {
  1., 0., 0., 0.,
  13.5/27., 13.5/27., 0., 0.,
@@ -21,7 +21,7 @@ static double q1_into_q1_refined_0[] =
  13.5/27., 0., 0., 13.5/27.,
 };
 
-static double q1_into_q1_refined_1[] =
+static const double q1_into_q1_refined_1[] =
 {
  13.5/27., 13.5/27., 0., 0.,
  0., 1., 0., 0.,
@@ -29,7 +29,7 @@ static double q1_into_q1_refined_1[] =
  6.75/27., 6.75/27., 6.75/27., 6.75/27.,
 };
 
-static double q1_into_q1_refined_2[] =
+static const double q1_into_q1_refined_2[] =
 {
  6.75/27., 6.75/27., 6.75/27., 6.75/27.,
  0., 13.5/27., 13.5/27., 0.,
@@ -37,7 +37,7 @@ static double q1_into_q1_refined_2[] =
  0., 0., 13.5/27., 13.5/27.,
 };
 
-static double q1_into_q1_refined_3[] =
+static const double q1_into_q1_refined_3[] =
 {
  13.5/27., 0., 0., 13.5/27.,
  6.75/27., 6.75/27., 6.75/27., 6.75/27.,
@@ -45,7 +45,7 @@ static double q1_into_q1_refined_3[] =
  0., 0., 0., 1.,
 };
 
-static double q2_into_q2_refined_0[] =
+static const double q2_into_q2_refined_0[] =
 {
  1., 0., 0., 0., 0., 0., 0., 0., 0.,
  0., 0., 0., 0., 1., 0., 0., 0., 0.,
@@ -58,7 +58,7 @@ static double q2_into_q2_refined_0[] =
  3.796875/27., -1.265625/27., 0.421875/27., -1.265625/27., 7.59375/27., -2.53125/27., -2.53125/27., 7.59375/27., 15.1875/27.,
 };
 
-static double q2_into_q2_refined_1[] =
+static const double q2_into_q2_refined_1[] =
 {
  0., 0., 0., 0., 1., 0., 0., 0., 0.,
  0., 1., 0., 0., 0., 0., 0., 0., 0.,
@@ -71,7 +71,7 @@ static double q2_into_q2_refined_1[] =
  -1.265625/27., 3.796875/27., -1.265625/27., 0.421875/27., 7.59375/27., 7.59375/27., -2.53125/27., -2.53125/27., 15.1875/27.,
 };
 
-static double q2_into_q2_refined_2[] =
+static const double q2_into_q2_refined_2[] =
 {
  0., 0., 0., 0., 0., 0., 0., 0., 1.,
  0., 0., 0., 0., 0., 1., 0., 0., 0.,
@@ -84,7 +84,7 @@ static double q2_into_q2_refined_2[] =
  0.421875/27., -1.265625/27., 3.796875/27., -1.265625/27., -2.53125/27., 7.59375/27., 7.59375/27., -2.53125/27., 15.1875/27.,
 };
 
-static double q2_into_q2_refined_3[] =
+static const double q2_into_q2_refined_3[] =
 {
  0., 0., 0., 0., 0., 0., 0., 1., 0.,
  0., 0., 0., 0., 0., 0., 0., 0., 1.,
@@ -97,7 +97,7 @@ static double q2_into_q2_refined_3[] =
  -1.265625/27., 0.421875/27., -1.265625/27., 3.796875/27., -2.53125/27., -2.53125/27., 7.59375/27., 7.59375/27., 15.1875/27.,
 };
 
-static double q3_into_q3_refined_0[] =
+static const double q3_into_q3_refined_0[] =
 {
  1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
  -1.6875/27., -1.6875/27., 0., 0., 15.1875/27., 15.1875/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
@@ -117,7 +117,7 @@ static double q3_into_q3_refined_0[] =
  0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0.,
 };
 
-static double q3_into_q3_refined_1[] =
+static const double q3_into_q3_refined_1[] =
 {
  -1.6875/27., -1.6875/27., 0., 0., 15.1875/27., 15.1875/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
  0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
@@ -137,7 +137,7 @@ static double q3_into_q3_refined_1[] =
  0., 0., 0., 0., 0., 0., 8.4375/27., 0., 0., 0., 1.6875/27., 0., -8.4375/27., 25.3125/27., 0., 0.,
 };
 
-static double q3_into_q3_refined_2[] =
+static const double q3_into_q3_refined_2[] =
 {
  0.10546875/27., 0.10546875/27., 0.10546875/27., 0.10546875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., 8.54296875/27., 8.54296875/27., 8.54296875/27., 8.54296875/27.,
  0., -1.6875/27., -1.6875/27., 0., 0., 0., 15.1875/27., 15.1875/27., 0., 0., 0., 0., 0., 0., 0., 0.,
@@ -157,7 +157,7 @@ static double q3_into_q3_refined_2[] =
  0.10546875/27., 0.52734375/27., 2.63671875/27., 0.52734375/27., -0.52734375/27., 1.58203125/27., -2.63671875/27., 7.91015625/27., -2.63671875/27., 7.91015625/27., -0.52734375/27., 1.58203125/27., 2.63671875/27., -7.91015625/27., -7.91015625/27., 23.73046875/27.,
 };
 
-static double q3_into_q3_refined_3[] =
+static const double q3_into_q3_refined_3[] =
 {
  -1.6875/27., 0., 0., -1.6875/27., 0., 0., 0., 0., 0., 0., 15.1875/27., 15.1875/27., 0., 0., 0., 0.,
  0.10546875/27., 0.10546875/27., 0.10546875/27., 0.10546875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., 8.54296875/27., 8.54296875/27., 8.54296875/27., 8.54296875/27.,
@@ -178,82 +178,34 @@ static double q3_into_q3_refined_3[] =
 };
 
 
-static double * embedding0[] =
-{
-  q1_into_q1_refined_0,
-  q2_into_q2_refined_0,
-  q3_into_q3_refined_0,
-  0,
-  0,
-  0,
-  0,
-  0,
-  0,
-  0
-};
+// embedding matrices
 
-static double * embedding1[] =
-{
-  q1_into_q1_refined_1,
-  q2_into_q2_refined_1,
-  q3_into_q3_refined_1,
-  0,
-  0,
-  0,
-  0,
-  0,
-  0,
-  0
-};
+const unsigned int FE_Q<2>::Matrices::n_embedding_matrices = 3;
 
-static double * embedding2[] =
+const double * const 
+FE_Q<2>::Matrices::embedding[][GeometryInfo<2>::children_per_cell] =
 {
-  q1_into_q1_refined_2,
-  q2_into_q2_refined_2,
-  q3_into_q3_refined_2,
-  0,
-  0,
-  0,
-  0,
-  0,
-  0,
-  0
+  { q1_into_q1_refined_0, q1_into_q1_refined_1, q1_into_q1_refined_2, q1_into_q1_refined_3 },
+  { q2_into_q2_refined_0, q2_into_q2_refined_1, q2_into_q2_refined_2, q2_into_q2_refined_3 },
+  { q3_into_q3_refined_0, q3_into_q3_refined_1, q3_into_q3_refined_2, q3_into_q3_refined_3 }
 };
 
-static double * embedding3[] =
-{
-  q1_into_q1_refined_3,
-  q2_into_q2_refined_3,
-  q3_into_q3_refined_3,
-  0,
-  0,
-  0,
-  0,
-  0,
-  0,
-  0
-};
 
-static double** embedding4;
-static double** embedding5;
-static double** embedding6;
-static double** embedding7;
 
 // Constraint matrices taken from Wolfgangs old version
-
-static double constraint_q1[] =
+static const double constraint_q1[] =
 {
   .5, .5
 };
 
-static double constraint_q2[] =
+static const double constraint_q2[] =
 {
   0., 0., 1.,
   .375, -.125, .75,
   -.125, .375, .75
 };
 
-static double constraint_q3[] =
+static const double constraint_q3[] =
 {
   -.0625, -.0625, .5625, .5625,
   .3125, .0625, .9375, -.3125,
@@ -262,7 +214,7 @@ static double constraint_q3[] =
   .0625, .3125, -.3125, 0.9375
 };
 
-static double constraint_q4[] =
+static const double constraint_q4[] =
 {
   0., 0., 0., 1., 0.,
   0.2734375, -0.0390625, 1.09375, -0.546875,  0.21875,
@@ -273,20 +225,15 @@ static double constraint_q4[] =
   -0.0390625, 0.2734375, 0.21875, -0.546875, 1.09375
 };
 
-static double* constraint_matrices[] =
+const double * const 
+FE_Q<2>::Matrices::constraint_matrices[] =
 {
   constraint_q1,
   constraint_q2,
   constraint_q3,
   constraint_q4,
-  0,
-  0,
-  0,
-  0,
-  0,
-  0,
-  0,
-  0,
-  0
 };
 
+const unsigned int 
+FE_Q<2>::Matrices::n_constraint_matrices = 4;
+
index 38791acf43ce3ad58050aa2cd14ead1c346d7b94..6159077b43465d24b72c235085f06b2a816ff47f 100644 (file)
@@ -14,7 +14,7 @@
 // Transfer matrices for finite elements
 
 
-static double q1_into_q1_refined_0[] =
+static const double q1_into_q1_refined_0[] =
 {
  1., 0., 0., 0., 0., 0., 0., 0.,
  13.5/27., 13.5/27., 0., 0., 0., 0., 0., 0.,
@@ -26,7 +26,7 @@ static double q1_into_q1_refined_0[] =
  6.75/27., 0., 0., 6.75/27., 6.75/27., 0., 0., 6.75/27.,
 };
 
-static double q1_into_q1_refined_1[] =
+static const double q1_into_q1_refined_1[] =
 {
  13.5/27., 13.5/27., 0., 0., 0., 0., 0., 0.,
  0., 1., 0., 0., 0., 0., 0., 0.,
@@ -38,7 +38,7 @@ static double q1_into_q1_refined_1[] =
  3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27.,
 };
 
-static double q1_into_q1_refined_2[] =
+static const double q1_into_q1_refined_2[] =
 {
  6.75/27., 6.75/27., 6.75/27., 6.75/27., 0., 0., 0., 0.,
  0., 13.5/27., 13.5/27., 0., 0., 0., 0., 0.,
@@ -50,7 +50,7 @@ static double q1_into_q1_refined_2[] =
  0., 0., 6.75/27., 6.75/27., 0., 0., 6.75/27., 6.75/27.,
 };
 
-static double q1_into_q1_refined_3[] =
+static const double q1_into_q1_refined_3[] =
 {
  13.5/27., 0., 0., 13.5/27., 0., 0., 0., 0.,
  6.75/27., 6.75/27., 6.75/27., 6.75/27., 0., 0., 0., 0.,
@@ -62,7 +62,7 @@ static double q1_into_q1_refined_3[] =
  0., 0., 0., 13.5/27., 0., 0., 0., 13.5/27.,
 };
 
-static double q1_into_q1_refined_4[] =
+static const double q1_into_q1_refined_4[] =
 {
  13.5/27., 0., 0., 0., 13.5/27., 0., 0., 0.,
  6.75/27., 6.75/27., 0., 0., 6.75/27., 6.75/27., 0., 0.,
@@ -74,7 +74,7 @@ static double q1_into_q1_refined_4[] =
  0., 0., 0., 0., 13.5/27., 0., 0., 13.5/27.,
 };
 
-static double q1_into_q1_refined_5[] =
+static const double q1_into_q1_refined_5[] =
 {
  6.75/27., 6.75/27., 0., 0., 6.75/27., 6.75/27., 0., 0.,
  0., 13.5/27., 0., 0., 0., 13.5/27., 0., 0.,
@@ -86,7 +86,7 @@ static double q1_into_q1_refined_5[] =
  0., 0., 0., 0., 6.75/27., 6.75/27., 6.75/27., 6.75/27.,
 };
 
-static double q1_into_q1_refined_6[] =
+static const double q1_into_q1_refined_6[] =
 {
  3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27.,
  0., 6.75/27., 6.75/27., 0., 0., 6.75/27., 6.75/27., 0.,
@@ -98,7 +98,7 @@ static double q1_into_q1_refined_6[] =
  0., 0., 0., 0., 0., 0., 13.5/27., 13.5/27.,
 };
 
-static double q1_into_q1_refined_7[] =
+static const double q1_into_q1_refined_7[] =
 {
  6.75/27., 0., 0., 6.75/27., 6.75/27., 0., 0., 6.75/27.,
  3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27.,
@@ -110,7 +110,7 @@ static double q1_into_q1_refined_7[] =
  0., 0., 0., 0., 0., 0., 0., 1.,
 };
 
-static double q2_into_q2_refined_0[] =
+static const double q2_into_q2_refined_0[] =
 {
  1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
  0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
@@ -141,7 +141,7 @@ static double q2_into_q2_refined_0[] =
  1.423828125/27., -0.474609375/27., 0.158203125/27., -0.474609375/27., -0.474609375/27., 0.158203125/27., -0.052734375/27., 0.158203125/27., 2.84765625/27., -0.94921875/27., -0.94921875/27., 2.84765625/27., -0.94921875/27., 0.31640625/27., 0.31640625/27., -0.94921875/27., 2.84765625/27., -0.94921875/27., 0.31640625/27., -0.94921875/27., 5.6953125/27., -1.8984375/27., 5.6953125/27., -1.8984375/27., -1.8984375/27., 5.6953125/27., 11.390625/27.,
 };
 
-static double q2_into_q2_refined_1[] =
+static const double q2_into_q2_refined_1[] =
 {
  0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
  0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
@@ -172,7 +172,7 @@ static double q2_into_q2_refined_1[] =
  -0.474609375/27., 1.423828125/27., -0.474609375/27., 0.158203125/27., 0.158203125/27., -0.474609375/27., 0.158203125/27., -0.052734375/27., 2.84765625/27., 2.84765625/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., 0.31640625/27., 0.31640625/27., -0.94921875/27., 2.84765625/27., -0.94921875/27., 0.31640625/27., 5.6953125/27., -1.8984375/27., 5.6953125/27., 5.6953125/27., -1.8984375/27., -1.8984375/27., 11.390625/27.,
 };
 
-static double q2_into_q2_refined_2[] =
+static const double q2_into_q2_refined_2[] =
 {
  0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0.,
  0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
@@ -203,7 +203,7 @@ static double q2_into_q2_refined_2[] =
  0.158203125/27., -0.474609375/27., 1.423828125/27., -0.474609375/27., -0.052734375/27., 0.158203125/27., -0.474609375/27., 0.158203125/27., -0.94921875/27., 2.84765625/27., 2.84765625/27., -0.94921875/27., 0.31640625/27., -0.94921875/27., -0.94921875/27., 0.31640625/27., 0.31640625/27., -0.94921875/27., 2.84765625/27., -0.94921875/27., 5.6953125/27., -1.8984375/27., -1.8984375/27., 5.6953125/27., 5.6953125/27., -1.8984375/27., 11.390625/27.,
 };
 
-static double q2_into_q2_refined_3[] =
+static const double q2_into_q2_refined_3[] =
 {
  0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
  0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0.,
@@ -234,7 +234,7 @@ static double q2_into_q2_refined_3[] =
  -0.474609375/27., 0.158203125/27., -0.474609375/27., 1.423828125/27., 0.158203125/27., -0.052734375/27., 0.158203125/27., -0.474609375/27., -0.94921875/27., -0.94921875/27., 2.84765625/27., 2.84765625/27., 0.31640625/27., 0.31640625/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., 0.31640625/27., -0.94921875/27., 2.84765625/27., 5.6953125/27., -1.8984375/27., -1.8984375/27., -1.8984375/27., 5.6953125/27., 5.6953125/27., 11.390625/27.,
 };
 
-static double q2_into_q2_refined_4[] =
+static const double q2_into_q2_refined_4[] =
 {
  0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
  0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0.,
@@ -265,7 +265,7 @@ static double q2_into_q2_refined_4[] =
  -0.474609375/27., 0.158203125/27., -0.052734375/27., 0.158203125/27., 1.423828125/27., -0.474609375/27., 0.158203125/27., -0.474609375/27., -0.94921875/27., 0.31640625/27., 0.31640625/27., -0.94921875/27., 2.84765625/27., -0.94921875/27., -0.94921875/27., 2.84765625/27., 2.84765625/27., -0.94921875/27., 0.31640625/27., -0.94921875/27., -1.8984375/27., 5.6953125/27., 5.6953125/27., -1.8984375/27., -1.8984375/27., 5.6953125/27., 11.390625/27.,
 };
 
-static double q2_into_q2_refined_5[] =
+static const double q2_into_q2_refined_5[] =
 {
  0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0.,
  0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
@@ -296,7 +296,7 @@ static double q2_into_q2_refined_5[] =
  0.158203125/27., -0.474609375/27., 0.158203125/27., -0.052734375/27., -0.474609375/27., 1.423828125/27., -0.474609375/27., 0.158203125/27., -0.94921875/27., -0.94921875/27., 0.31640625/27., 0.31640625/27., 2.84765625/27., 2.84765625/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., 2.84765625/27., -0.94921875/27., 0.31640625/27., -1.8984375/27., 5.6953125/27., 5.6953125/27., 5.6953125/27., -1.8984375/27., -1.8984375/27., 11.390625/27.,
 };
 
-static double q2_into_q2_refined_6[] =
+static const double q2_into_q2_refined_6[] =
 {
  0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.,
  0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0.,
@@ -327,7 +327,7 @@ static double q2_into_q2_refined_6[] =
  -0.052734375/27., 0.158203125/27., -0.474609375/27., 0.158203125/27., 0.158203125/27., -0.474609375/27., 1.423828125/27., -0.474609375/27., 0.31640625/27., -0.94921875/27., -0.94921875/27., 0.31640625/27., -0.94921875/27., 2.84765625/27., 2.84765625/27., -0.94921875/27., 0.31640625/27., -0.94921875/27., 2.84765625/27., -0.94921875/27., -1.8984375/27., 5.6953125/27., -1.8984375/27., 5.6953125/27., 5.6953125/27., -1.8984375/27., 11.390625/27.,
 };
 
-static double q2_into_q2_refined_7[] =
+static const double q2_into_q2_refined_7[] =
 {
  0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.,
  0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.,
@@ -360,146 +360,27 @@ static double q2_into_q2_refined_7[] =
 
 
 
+// embedding matrices
 
-static double * embedding0[] =
-{
-  q1_into_q1_refined_0,
-  q2_into_q2_refined_0,
-//    q3_into_q3_refined_0,
-//    q4_into_q4_refined_0,
-  0,
-  0,
-  0,
-  0,
-  0,
-  0,
-  0,
-  0,
-  0
-};
-
-static double * embedding1[] =
-{
-  q1_into_q1_refined_1,
-  q2_into_q2_refined_1,
-//    q3_into_q3_refined_1,
-//    q4_into_q4_refined_1,
-  0,
-  0,
-  0,
-  0,
-  0,
-  0,
-  0,
-  0,
-  0
-};
+const unsigned int 
+FE_Q<3>::Matrices::n_embedding_matrices = 2;
 
-static double * embedding2[] =
+const double * const 
+FE_Q<3>::Matrices::embedding[][GeometryInfo<3>::children_per_cell] =
 {
-  q1_into_q1_refined_2,
-  q2_into_q2_refined_2,
-//    q3_into_q3_refined_2,
-//    q4_into_q4_refined_2,
-  0,
-  0,
-  0,
-  0,
-  0,
-  0,
-  0,
-  0,
-  0
+  { q1_into_q1_refined_0, q1_into_q1_refined_1, q1_into_q1_refined_2, q1_into_q1_refined_3,
+    q1_into_q1_refined_4, q1_into_q1_refined_5, q1_into_q1_refined_6, q1_into_q1_refined_7 },
+  { q2_into_q2_refined_0, q2_into_q2_refined_1, q2_into_q2_refined_2, q2_into_q2_refined_3,
+    q2_into_q2_refined_4, q2_into_q2_refined_5, q2_into_q2_refined_6, q2_into_q2_refined_7 }
 };
 
-static double * embedding3[] =
-{
-  q1_into_q1_refined_3,
-  q2_into_q2_refined_3,
-//    q3_into_q3_refined_3,
-//    q4_into_q4_refined_3,
-  0,
-  0,
-  0,
-  0,
-  0,
-  0,
-  0,
-  0,
-  0
-};
-
-static double * embedding4[] =
-{
-  q1_into_q1_refined_4,
-  q2_into_q2_refined_4,
-//    q3_into_q3_refined_4,
-//    q4_into_q4_refined_4,
-  0,
-  0,
-  0,
-  0,
-  0,
-  0,
-  0,
-  0,
-  0
-};
-
-static double * embedding5[] =
-{
-  q1_into_q1_refined_5,
-  q2_into_q2_refined_5,
-//    q3_into_q3_refined_5,
-//    q4_into_q4_refined_5,
-  0,
-  0,
-  0,
-  0,
-  0,
-  0,
-  0,
-  0,
-  0
-};
-
-static double * embedding6[] =
-{
-  q1_into_q1_refined_6,
-  q2_into_q2_refined_6,
-//    q3_into_q3_refined_6,
-//    q4_into_q4_refined_6,
-  0,
-  0,
-  0,
-  0,
-  0,
-  0,
-  0,
-  0,
-  0
-};
-
-static double * embedding7[] =
-{
-  q1_into_q1_refined_7,
-  q2_into_q2_refined_7,
-//    q3_into_q3_refined_7,
-//    q4_into_q4_refined_7,
-  0,
-  0,
-  0,
-  0,
-  0,
-  0,
-  0,
-  0,
-  0
-};
 
 
+// Constraint matrices taken from Wolfgangs old version
+const unsigned int 
+FE_Q<3>::Matrices::n_constraint_matrices = 2;
 
-static double constraint_q1[] =
+static const double constraint_q1[] =
 {
   .25,.25,.25,.25,
   .5,.5,0.,0.,
@@ -508,7 +389,7 @@ static double constraint_q1[] =
   .5,0.,0.,.5
 };
 
-static double constraint_q2[] =
+static const double constraint_q2[] =
 {
   0,0,0,0,0,0,0,0,1,
   0,0,0,0,1,0,0,0,0,
@@ -533,20 +414,10 @@ static double constraint_q2[] =
   -.046875,.015625,-.046875,.140625,-.09375,-.09375,.28125,.28125,.5625
 };
 
-static double* constraint_matrices[] =
+const double * const 
+FE_Q<3>::Matrices::constraint_matrices[] =
 {
   constraint_q1,
-  constraint_q2,
-  0,
-  0,
-  0,
-  0,
-  0,
-  0,
-  0,
-  0,
-  0
-//  constraint_q3,
-//  constraint_q4
+  constraint_q2
 };
 

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.