]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
No longer preset the correct sizes of matrices in the FE base class, but force those...
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 30 Sep 2002 14:07:58 +0000 (14:07 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 30 Sep 2002 14:07:58 +0000 (14:07 +0000)
git-svn-id: https://svn.dealii.org/trunk@6557 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/fe/fe_base.h
deal.II/deal.II/source/fe/fe.cc
deal.II/deal.II/source/fe/fe_dgp.cc
deal.II/deal.II/source/fe/fe_dgq.cc
deal.II/deal.II/source/fe/fe_nedelec.cc
deal.II/deal.II/source/fe/fe_q.cc
deal.II/deal.II/source/fe/fe_system.cc

index 271eb473212ca27c662246f44da22f02a012e7c9..297e0846aced3e1eb56afc0becf106646de0f51f 100644 (file)
@@ -18,6 +18,7 @@
 #include <base/subscriptor.h>
 #include <base/point.h>
 #include <base/tensor.h>
+#include <base/table.h>
 #include <grid/geometry_info.h>
 #include <lac/full_matrix.h>
 #include <fe/fe_update_flags.h>
@@ -1106,6 +1107,21 @@ class FiniteElementBase : public Subscriptor,
                                      */
     FullMatrix<double> interface_constraints;
 
+                                     /**
+                                      * Return the size of interface
+                                      * constraint matrices. Since
+                                      * this is needed in every
+                                      * derived finite element class
+                                      * when initializing their size,
+                                      * it is placed into this
+                                      * function, to avoid having to
+                                      * recomputed the
+                                      * dimension-dependent size of
+                                      * these matrices each time.
+                                      */
+    TableIndices<2>
+    interface_constraints_size () const;
+    
                                     /**
                                      * Store what
                                      * @p{system_to_component_index}
index cc08011d62d40a1b8b638300ac8a9971f58520ea..9c59f139e5c0be89a8cffe35fb857f25d68c8c91 100644 (file)
@@ -27,6 +27,9 @@
 #include <numeric>
 
 
+//TODO[GK]: in all FEs, somehow make the handling of first_cell more intuitive. presently, this flag is cleared at the end of the fill_fe_*_values function, but this leads to confusion, if for example in the FESystems class, we ask for the present values of the update flags, using update_current, because the element thinks it is already on the second cell and the FESystem does not want to do some actions it is supposed to do. thus, in the respective FESystem function, we store these update_values once before calling the fill_* functions of the subelements, and reuse the stored values. if we would request their values only after initializing the base element (i.e. requesting them at the point _where we actually need them_!), this would fail, and we would not copy the shape_values, for example.
+//TODO[GK]: the solution would probably to not set and clear the first_cell flag in the FEs themselves, but to have a function which sets or clears this flag, and call it from the FEValues object, which should always know what the present state is. this way, we could be more sure that this flag is actually consistent with the present state, and among the various elements that are involved in a FESystem.
+
 // if necessary try to work around a bug in the IBM xlC compiler
 #ifdef XLC_WORK_AROUND_STD_BUG
 using namespace std;
@@ -151,43 +154,9 @@ FiniteElementBase (const FiniteElementData<dim> &fe_data,
              ExcInternalError());      
     };
   
-  for (unsigned int i=0; i<GeometryInfo<dim>::children_per_cell; ++i) 
-    {
-      restriction[i].reinit (this->dofs_per_cell, this->dofs_per_cell);
-      prolongation[i].reinit (this->dofs_per_cell, this->dofs_per_cell);
-    };
-
-                                  // first set sizes of some
-                                  // matrices. they will be filled by
-                                  // derived classes, or re-set to
-                                  // zero size of not defined
-  switch (dim)
-    {
-      case 1:
-       Assert ((interface_constraints.m() == 0) &&
-               (interface_constraints.n() == 0),
-               ExcInternalError());
-       break;
-           
-      case 2:
-       interface_constraints.reinit (this->dofs_per_vertex
-                                     +2*this->dofs_per_line,
-                                     this->dofs_per_face);
-       break;
-
-      case 3:
-       interface_constraints.reinit (5*this->dofs_per_vertex +
-                                     12*this->dofs_per_line  +
-                                     4*this->dofs_per_quad,
-                                     this->dofs_per_face);
-       break;
-
-      default:
-       Assert (false, ExcNotImplemented());
-    };
-           
-                                  // this is the default way, if
-                                  // there is only one component; if
+                                   // initialize some tables in the
+                                  // default way, i.e. if there is
+                                  // only one (vector-)component; if
                                   // there are several, then the
                                   // constructor of the derived class
                                   // needs to overwrite these arrays
@@ -313,6 +282,32 @@ FiniteElementBase<dim>::constraints () const
 
 
 
+template <int dim>
+TableIndices<2>
+FiniteElementBase<dim>::interface_constraints_size () const 
+{
+  switch (dim)
+    {
+      case 1:
+            return TableIndices<2> (0U, 0U);
+      case 2:
+            return TableIndices<2> (this->dofs_per_vertex +
+                                    2*this->dofs_per_line,
+                                    this->dofs_per_face);
+      case 3:
+            return TableIndices<2> (5*this->dofs_per_vertex +
+                                    12*this->dofs_per_line  +
+                                    4*this->dofs_per_quad,
+                                    this->dofs_per_face);
+      default:
+            Assert (false, ExcNotImplemented());
+    };
+  return TableIndices<2> (-1, -1);
+};
+
+                                   
+
+
 template <int dim>
 bool FiniteElementBase<dim>::operator == (const FiniteElementBase<dim> &f) const
 {
index b8331e09931cc75ee2d26cdf839c77d68b80f4e5..7018ccf84cdc1fb860044e0190d13a1bff62dc1d 100644 (file)
@@ -37,7 +37,11 @@ FE_DGP<dim>::FE_DGP (const unsigned int degree)
   if ((degree < Matrices::n_embedding_matrices) &&
       (Matrices::embedding[degree][0] != 0))
     for (unsigned int c=0; c<GeometryInfo<dim>::children_per_cell; ++c)
-      this->prolongation[c].fill (Matrices::embedding[degree][c]);
+      {
+        this->prolongation[c].reinit (this->dofs_per_cell,
+                                      this->dofs_per_cell);
+        this->prolongation[c].fill (Matrices::embedding[degree][c]);
+      }
   else
     for (unsigned int i=0; i<GeometryInfo<dim>::children_per_cell;++i)
       this->prolongation[i].reinit(0,0);
@@ -50,6 +54,7 @@ FE_DGP<dim>::FE_DGP (const unsigned int degree)
                                    //
                                    // if it were, then the following
                                    // snippet would be the right code
+//  this->restriction[i].reinit (this->dofs_per_cell, this->dofs_per_cell);
 //    if ((degree < Matrices::n_projection_matrices) &&
 //        (Matrices::projection_matrices[degree] != 0))
 //      {
@@ -59,10 +64,6 @@ FE_DGP<dim>::FE_DGP (const unsigned int degree)
 //                                  // matrix undefined, set size to zero
 //      for (unsigned int i=0;i<GeometryInfo<dim>::children_per_cell;++i)
 //        restriction[i].reinit(0, 0);
-                                   // since not implemented, set to
-                                   // "empty"
-  for (unsigned int i=0;i<GeometryInfo<dim>::children_per_cell;++i)
-    restriction[i].reinit(0, 0);
 
                                    // note further, that these
                                    // elements have neither support
index 6ddd18b9301929cc09d85749d228a8df20f84d57..9f7f38eb8a6778fe5ff874635bafdf2cd3b4a8bd 100644 (file)
@@ -57,14 +57,21 @@ FE_DGQ<dim>::FE_DGQ (const unsigned int degree)
   if ((degree < Matrices::n_embedding_matrices) &&
       (Matrices::embedding[degree] != 0))
     {
+      for (unsigned int i=0; i<GeometryInfo<dim>::children_per_cell; ++i)
+        this->prolongation[i].reinit (this->dofs_per_cell,
+                                      this->dofs_per_cell);
       this->prolongation[0].fill (Matrices::embedding[degree]);
       switch (dim)
        {
          case 1:
+          {
            this->prolongation[1].fill_permutation (this->prolongation[0],
                                                    right, right);
            break;
+          };
+          
          case 2:
+          {
            this->prolongation[1].fill_permutation (this->prolongation[0],
                                                    right, right);
            this->prolongation[2].fill_permutation (this->prolongation[1],
@@ -72,7 +79,10 @@ FE_DGQ<dim>::FE_DGQ (const unsigned int degree)
            this->prolongation[3].fill_permutation (this->prolongation[2],
                                                    right, right);
            break;
+          };
+          
          case 3:
+          {
            this->prolongation[1].fill_permutation (this->prolongation[0],
                                                    right, right);
            this->prolongation[5].fill_permutation (this->prolongation[1],
@@ -88,14 +98,16 @@ FE_DGQ<dim>::FE_DGQ (const unsigned int degree)
            this->prolongation[2].fill_permutation (this->prolongation[6],
                                                    top, top);
            break;
+          };
+          
          default:
            Assert (false, ExcNotImplemented());
        }
     }
   else
-                                    // matrix undefined, set size to zero
-    for (unsigned int i=0;i<GeometryInfo<dim>::children_per_cell;++i)
-      this->prolongation[i].reinit(0, 0);
+                                    // matrix undefined, leave matrix
+                                    // at size zero
+    ;
 
                                   // same as above: copy over matrix
                                   // from predefined values and
@@ -103,14 +115,21 @@ FE_DGQ<dim>::FE_DGQ (const unsigned int degree)
   if ((degree < Matrices::n_projection_matrices) &&
       (Matrices::projection_matrices[degree] != 0))
     {
+      for (unsigned int i=0; i<GeometryInfo<dim>::children_per_cell; ++i)
+        this->restriction[i].reinit (this->dofs_per_cell,
+                                     this->dofs_per_cell);
       this->restriction[0].fill (Matrices::projection_matrices[degree]);
       switch (dim)
        {
          case 1:
+          {
            this->restriction[1].fill_permutation (this->restriction[0],
                                                   right, right);
            break;
+          };
+          
          case 2:
+          {
            this->restriction[1].fill_permutation (this->restriction[0],
                                                   right, right);
            this->restriction[2].fill_permutation (this->restriction[1],
@@ -118,7 +137,10 @@ FE_DGQ<dim>::FE_DGQ (const unsigned int degree)
            this->restriction[3].fill_permutation (this->restriction[2],
                                                   right, right);
            break;
+          };
+          
          case 3:
+          {
            this->restriction[1].fill_permutation (this->restriction[0],
                                                   right, right);
            this->restriction[5].fill_permutation (this->restriction[1],
@@ -134,15 +156,16 @@ FE_DGQ<dim>::FE_DGQ (const unsigned int degree)
            this->restriction[2].fill_permutation (this->restriction[6],
                                                   top, top);
            break;
+          };
+          
          default:
            Assert (false, ExcNotImplemented());
        }
     }
   else
-                                    // matrix undefined, set size to zero
-    for (unsigned int i=0;i<GeometryInfo<dim>::children_per_cell;++i)
-      this->restriction[i].reinit(0, 0);
-
+                                    // matrix undefined, leave matrix
+                                    // at size zero
+    ;
   
                                   // finally fill in support points
   if (degree == 0)
index 22ada26dc7d6d788a09c82c8aab9b37691174520..079b01e0f18bf4bf3bead99b71ccd84c97810e1f 100644 (file)
@@ -12,8 +12,7 @@
 //----------------------------------------------------------------
 
 #include <base/quadrature.h>
-#include <base/polynomial.h>
-#include <base/tensor_product_polynomials.h>
+#include <base/table.h>
 #include <grid/tria.h>
 #include <grid/tria_iterator.h>
 #include <dofs/dof_accessor.h>
@@ -36,12 +35,14 @@ FE_Nedelec<dim>::FE_Nedelec (const unsigned int degree)
   Assert (dim >= 2, ExcNotUsefulInThisDimension());
   
                                   // copy constraint matrices if they
-                                  // are defined. otherwise set them
-                                  // to invalid size
+                                  // are defined. otherwise leave them
+                                  // at zero size
   if (degree<Matrices::n_constraint_matrices+1)
-    this->interface_constraints.fill (Matrices::constraint_matrices[degree-1]);
-  else
-    this->interface_constraints.reinit(0,0);
+    {
+      this->interface_constraints.
+        TableBase<2,double>::reinit (this->interface_constraints_size());
+      this->interface_constraints.fill (Matrices::constraint_matrices[degree-1]);
+    };
 
                                   // next copy over embedding
                                   // matrices if they are defined
@@ -50,6 +51,8 @@ FE_Nedelec<dim>::FE_Nedelec (const unsigned int degree)
     for (unsigned int c=0; c<GeometryInfo<dim>::children_per_cell; ++c)
       {
                                          // copy
+        this->prolongation[c].reinit (this->dofs_per_cell,
+                                      this->dofs_per_cell);
         this->prolongation[c].fill (Matrices::embedding[degree-1][c]);
                                          // and make sure that the row
                                          // sum is 0.5 (for usual
@@ -68,10 +71,7 @@ FE_Nedelec<dim>::FE_Nedelec (const unsigned int degree)
             Assert (std::fabs(sum-.5) < 1e-14,
                     ExcInternalError());
           };
-      }
-  else
-    for (unsigned int i=0; i<GeometryInfo<dim>::children_per_cell;++i)
-      this->prolongation[i].reinit(0,0);
+      };
 
                                   // then fill restriction
                                   // matrices. they are hardcoded for
@@ -165,6 +165,10 @@ FE_Nedelec<dim>::FE_Nedelec (const unsigned int degree)
                                                // always, in the
                                                // canonical direction
                                                // of lines
+              for (unsigned int c=0; c<GeometryInfo<dim>::children_per_cell; ++c)
+                this->restriction[c].reinit (this->dofs_per_cell,
+                                             this->dofs_per_cell);
+              
              this->restriction[0](0,0) = 2.;
              this->restriction[1](1,1) = 2.;
              this->restriction[3](2,2) = 2.;
@@ -177,17 +181,13 @@ FE_Nedelec<dim>::FE_Nedelec (const unsigned int degree)
            {
                                               // in case we don't
                                               // have the matrices
-                                              // (yet), set them to
-                                              // impossible
-                                              // values. this does
-                                              // not prevent the use
-                                              // of this FE, but will
+                                              // (yet), leave them
+                                              // empty. this does not
+                                              // prevent the use of
+                                              // this FE, but will
                                               // prevent the use of
                                               // these matrices
-             for (unsigned int i=0;
-                  i<GeometryInfo<dim>::children_per_cell;
-                  ++i)
-               this->restriction[i].reinit(0,0);
+              break;
            };
          };
        
@@ -206,6 +206,9 @@ FE_Nedelec<dim>::FE_Nedelec (const unsigned int degree)
                                               // cell to get at the
                                               // values of each of
                                               // the 12 lines
+              for (unsigned int c=0; c<GeometryInfo<dim>::children_per_cell; ++c)
+                this->restriction[c].reinit (this->dofs_per_cell,
+                                             this->dofs_per_cell);
              this->restriction[0](0,0) = 2.;
              this->restriction[0](3,3) = 2.;
              this->restriction[1](1,1) = 2.;
@@ -228,17 +231,13 @@ FE_Nedelec<dim>::FE_Nedelec (const unsigned int degree)
            {
                                               // in case we don't
                                               // have the matrices
-                                              // (yet), set them to
-                                              // impossible
-                                              // values. this does
-                                              // not prevent the use
-                                              // of this FE, but will
+                                              // (yet), leave them
+                                              // empty. this does not
+                                              // prevent the use of
+                                              // this FE, but will
                                               // prevent the use of
                                               // these matrices
-             for (unsigned int i=0;
-                  i<GeometryInfo<dim>::children_per_cell;
-                  ++i)
-               this->restriction[i].reinit(0,0);
+              break;
            };
          };
        
@@ -253,6 +252,14 @@ FE_Nedelec<dim>::FE_Nedelec (const unsigned int degree)
                                   // on cell and face
   initialize_unit_support_points ();
   initialize_unit_face_support_points ();
+
+                                   // then make
+                                   // system_to_component_table
+                                   // invalid, since this has no
+                                   // meaning for the present element
+  std::vector<std::pair<unsigned,unsigned> > tmp1, tmp2;
+  this->system_to_component_table.swap (tmp1);
+  this->face_system_to_component_table.swap (tmp2);
 };
 
 
index 845629644ec0a6e7fe6cb124d330e7d0cad3dd91..450063a9e568e65aca2d28709d19a63752c38c79 100644 (file)
@@ -49,22 +49,25 @@ FE_Q<dim>::FE_Q (const unsigned int degree)
     renumber_inverse[renumber[i]]=i;
 
                                   // copy constraint matrices if they
-                                  // are defined. otherwise set them
-                                  // to invalid size
+                                  // are defined. otherwise leave them
+                                  // at invalid size
   if ((dim>1) && (degree<Matrices::n_constraint_matrices+1))
-    this->interface_constraints.fill (Matrices::constraint_matrices[degree-1]);
-  else
-    this->interface_constraints.reinit(0,0);
+    {
+      this->interface_constraints.
+        TableBase<2,double>::reinit (this->interface_constraints_size());
+      this->interface_constraints.fill (Matrices::constraint_matrices[degree-1]);
+    };
 
                                   // next copy over embedding
                                   // matrices if they are defined
   if ((degree < Matrices::n_embedding_matrices+1) &&
       (Matrices::embedding[degree-1][0] != 0))
     for (unsigned int c=0; c<GeometryInfo<dim>::children_per_cell; ++c)
-      this->prolongation[c].fill (Matrices::embedding[degree-1][c]);
-  else
-    for (unsigned int i=0; i<GeometryInfo<dim>::children_per_cell;++i)
-      this->prolongation[i].reinit(0,0);
+      {
+        this->prolongation[c].reinit (this->dofs_per_cell,
+                                      this->dofs_per_cell);
+        this->prolongation[c].fill (Matrices::embedding[degree-1][c]);
+      };
 
                                   // then fill restriction
                                   // matrices. they are hardcoded for
@@ -77,548 +80,547 @@ FE_Q<dim>::FE_Q (const unsigned int degree)
                                   // them element-wise is more
                                   // expensive than just setting the
                                   // nonzero elements as done here
+  for (unsigned int c=0; c<GeometryInfo<dim>::children_per_cell; ++c)
+    this->restriction[c].reinit (this->dofs_per_cell, this->dofs_per_cell);
   switch (dim)
     {
       case 1:  // 1d
-    {
-      switch (degree)
-       {
-         case 1:
-           this->restriction[0](0,0) = 1;
-           this->restriction[1](1,1) = 1;
-           break;
-         case 2:
-           this->restriction[0](0,0) = 1;
-           this->restriction[0](2,1) = 1;
-           this->restriction[1](1,1) = 1;
-           this->restriction[1](1,1) = 1;
-           break;
-         case 3:
-           this->restriction[0](0,0) = 1;
-           this->restriction[0](2,3) = 1;
-           this->restriction[1](1,1) = 1;
-           this->restriction[1](3,2) = 1;
-           break;
-         case 4:
-           this->restriction[0](0,0) = 1;
-           this->restriction[0](2,3) = 1;
-           this->restriction[0](3,1) = 1;
-           this->restriction[1](1,1) = 1;
-           this->restriction[1](3,0) = 1;
-           this->restriction[1](4,3) = 1;
-           break;
+      {
+        switch (degree)
+          {
+            case 1:
+                  this->restriction[0](0,0) = 1;
+                  this->restriction[1](1,1) = 1;
+                  break;
+            case 2:
+                  this->restriction[0](0,0) = 1;
+                  this->restriction[0](2,1) = 1;
+                  this->restriction[1](1,1) = 1;
+                  this->restriction[1](1,1) = 1;
+                  break;
+            case 3:
+                  this->restriction[0](0,0) = 1;
+                  this->restriction[0](2,3) = 1;
+                  this->restriction[1](1,1) = 1;
+                  this->restriction[1](3,2) = 1;
+                  break;
+            case 4:
+                  this->restriction[0](0,0) = 1;
+                  this->restriction[0](2,3) = 1;
+                  this->restriction[0](3,1) = 1;
+                  this->restriction[1](1,1) = 1;
+                  this->restriction[1](3,0) = 1;
+                  this->restriction[1](4,3) = 1;
+                  break;
                  
-         default:
-       {
-                                          // in case we don't
-                                          // have the matrices
-                                          // (yet), set them to
-                                          // impossible
-                                          // values. this does
-                                          // not prevent the use
-                                          // of this FE, but will
-                                          // prevent the use of
-                                          // these matrices
-         for (unsigned int i=0;
-              i<GeometryInfo<dim>::children_per_cell;
-              ++i)
-           this->restriction[i].reinit(0,0);
-       };
-       }
-      break;
-    };
+            default:
+            {
+                                               // in case we don't
+                                               // have the matrices
+                                               // (yet), reset them to
+                                               // zero size. this does
+                                               // not prevent the use
+                                               // of this FE, but will
+                                               // prevent the use of
+                                               // these matrices
+              for (unsigned int i=0;
+                   i<GeometryInfo<dim>::children_per_cell;
+                   ++i)
+                this->restriction[i].reinit(0,0);
+            };
+          }
+        break;
+      };
       
       case 2:  // 2d
-    {
-      switch (degree)
-       {
-         case 1:
-           this->restriction[0](0,0) = 1;
-           this->restriction[1](1,1) = 1;
-           this->restriction[2](2,2) = 1;
-           this->restriction[3](3,3) = 1;
-           break;
-         case 2:
-           this->restriction[0](0,0) = 1;
-           this->restriction[0](4,1) = 1;
-           this->restriction[0](7,3) = 1;
-           this->restriction[0](8,2) = 1;
-           this->restriction[1](1,1) = 1;
-           this->restriction[1](4,0) = 1;
-           this->restriction[1](5,2) = 1;
-           this->restriction[1](8,3) = 1;
-           this->restriction[2](2,2) = 1;
-           this->restriction[2](5,1) = 1;
-           this->restriction[2](6,3) = 1;
-           this->restriction[2](8,0) = 1;
-           this->restriction[3](3,3) = 1;
-           this->restriction[3](6,2) = 1;
-           this->restriction[3](7,0) = 1;
-           this->restriction[3](8,1) = 1;
-           break;
-         case 3:
-           this->restriction[0](0,0) = 1;
-           this->restriction[0](4,5) = 1;
-           this->restriction[0](10,11) = 1;
-           this->restriction[0](12,15) = 1;
-           this->restriction[1](1,1) = 1;
-           this->restriction[1](5,4) = 1;
-           this->restriction[1](6,7) = 1;
-           this->restriction[1](13,14) = 1;
-           this->restriction[2](2,2) = 1;
-           this->restriction[2](7,6) = 1;
-           this->restriction[2](9,8) = 1;
-           this->restriction[2](15,12) = 1;
-           this->restriction[3](3,3) = 1;
-           this->restriction[3](8,9) = 1;
-           this->restriction[3](11,10) = 1;
-           this->restriction[3](14,13) = 1;
-           break;
-         case 4:
-           this->restriction[0](0,0) = 1;
-           this->restriction[0](4,5) = 1;
-           this->restriction[0](5,1) = 1;
-           this->restriction[0](13,14) = 1;
-           this->restriction[0](14,3) = 1;
-           this->restriction[0](16,20) = 1;
-           this->restriction[0](17,8) = 1;
-           this->restriction[0](19,11) = 1;
-           this->restriction[0](20,2) = 1;
-           this->restriction[1](1,1) = 1;
-           this->restriction[1](5,0) = 1;
-           this->restriction[1](6,5) = 1;
-           this->restriction[1](7,8) = 1;
-           this->restriction[1](8,2) = 1;
-           this->restriction[1](17,14) = 1;
-           this->restriction[1](18,20) = 1;
-           this->restriction[1](20,3) = 1;
-           this->restriction[1](21,11) = 1;
-           this->restriction[2](2,2) = 1;
-           this->restriction[2](8,1) = 1;
-           this->restriction[2](9,8) = 1;
-           this->restriction[2](11,3) = 1;
-           this->restriction[2](12,11) = 1;
-           this->restriction[2](20,0) = 1;
-           this->restriction[2](21,5) = 1;
-           this->restriction[2](23,14) = 1;
-           this->restriction[2](24,20) = 1;
-           this->restriction[3](3,3) = 1;
-           this->restriction[3](10,11) = 1;
-           this->restriction[3](11,2) = 1;
-           this->restriction[3](14,0) = 1;
-           this->restriction[3](15,14) = 1;
-           this->restriction[3](19,5) = 1;
-           this->restriction[3](20,1) = 1;
-           this->restriction[3](22,20) = 1;
-           this->restriction[3](23,8) = 1;
-           break;
-
-         default:
-       {
-                                          // in case we don't
-                                          // have the matrices
-                                          // (yet), set them to
-                                          // impossible
-                                          // values. this does
-                                          // not prevent the use
-                                          // of this FE, but will
-                                          // prevent the use of
-                                          // these matrices
-         for (unsigned int i=0;
-              i<GeometryInfo<dim>::children_per_cell;
-              ++i)
-           this->restriction[i].reinit(0,0);
-       };
-       }
-      break;
-    };
+      {
+        switch (degree)
+          {
+            case 1:
+                  this->restriction[0](0,0) = 1;
+                  this->restriction[1](1,1) = 1;
+                  this->restriction[2](2,2) = 1;
+                  this->restriction[3](3,3) = 1;
+                  break;
+            case 2:
+                  this->restriction[0](0,0) = 1;
+                  this->restriction[0](4,1) = 1;
+                  this->restriction[0](7,3) = 1;
+                  this->restriction[0](8,2) = 1;
+                  this->restriction[1](1,1) = 1;
+                  this->restriction[1](4,0) = 1;
+                  this->restriction[1](5,2) = 1;
+                  this->restriction[1](8,3) = 1;
+                  this->restriction[2](2,2) = 1;
+                  this->restriction[2](5,1) = 1;
+                  this->restriction[2](6,3) = 1;
+                  this->restriction[2](8,0) = 1;
+                  this->restriction[3](3,3) = 1;
+                  this->restriction[3](6,2) = 1;
+                  this->restriction[3](7,0) = 1;
+                  this->restriction[3](8,1) = 1;
+                  break;
+            case 3:
+                  this->restriction[0](0,0) = 1;
+                  this->restriction[0](4,5) = 1;
+                  this->restriction[0](10,11) = 1;
+                  this->restriction[0](12,15) = 1;
+                  this->restriction[1](1,1) = 1;
+                  this->restriction[1](5,4) = 1;
+                  this->restriction[1](6,7) = 1;
+                  this->restriction[1](13,14) = 1;
+                  this->restriction[2](2,2) = 1;
+                  this->restriction[2](7,6) = 1;
+                  this->restriction[2](9,8) = 1;
+                  this->restriction[2](15,12) = 1;
+                  this->restriction[3](3,3) = 1;
+                  this->restriction[3](8,9) = 1;
+                  this->restriction[3](11,10) = 1;
+                  this->restriction[3](14,13) = 1;
+                  break;
+            case 4:
+                  this->restriction[0](0,0) = 1;
+                  this->restriction[0](4,5) = 1;
+                  this->restriction[0](5,1) = 1;
+                  this->restriction[0](13,14) = 1;
+                  this->restriction[0](14,3) = 1;
+                  this->restriction[0](16,20) = 1;
+                  this->restriction[0](17,8) = 1;
+                  this->restriction[0](19,11) = 1;
+                  this->restriction[0](20,2) = 1;
+                  this->restriction[1](1,1) = 1;
+                  this->restriction[1](5,0) = 1;
+                  this->restriction[1](6,5) = 1;
+                  this->restriction[1](7,8) = 1;
+                  this->restriction[1](8,2) = 1;
+                  this->restriction[1](17,14) = 1;
+                  this->restriction[1](18,20) = 1;
+                  this->restriction[1](20,3) = 1;
+                  this->restriction[1](21,11) = 1;
+                  this->restriction[2](2,2) = 1;
+                  this->restriction[2](8,1) = 1;
+                  this->restriction[2](9,8) = 1;
+                  this->restriction[2](11,3) = 1;
+                  this->restriction[2](12,11) = 1;
+                  this->restriction[2](20,0) = 1;
+                  this->restriction[2](21,5) = 1;
+                  this->restriction[2](23,14) = 1;
+                  this->restriction[2](24,20) = 1;
+                  this->restriction[3](3,3) = 1;
+                  this->restriction[3](10,11) = 1;
+                  this->restriction[3](11,2) = 1;
+                  this->restriction[3](14,0) = 1;
+                  this->restriction[3](15,14) = 1;
+                  this->restriction[3](19,5) = 1;
+                  this->restriction[3](20,1) = 1;
+                  this->restriction[3](22,20) = 1;
+                  this->restriction[3](23,8) = 1;
+                  break;
+
+            default:
+            {
+                                               // in case we don't
+                                               // have the matrices
+                                               // (yet), reset them to
+                                               // zero size. this does
+                                               // not prevent the use
+                                               // of this FE, but will
+                                               // prevent the use of
+                                               // these matrices
+              for (unsigned int i=0;
+                   i<GeometryInfo<dim>::children_per_cell;
+                   ++i)
+                this->restriction[i].reinit(0,0);
+            };
+          }
+        break;
+      };
       
       case 3:  // 3d
-    {
-      switch (degree)
-       {
-         case 1:
-           this->restriction[0](0,0) = 1;
-           this->restriction[1](1,1) = 1;
-           this->restriction[2](2,2) = 1;
-           this->restriction[3](3,3) = 1;
-           this->restriction[4](4,4) = 1;
-           this->restriction[5](5,5) = 1;
-           this->restriction[6](6,6) = 1;
-           this->restriction[7](7,7) = 1;
-           break;
-         case 2:
-           this->restriction[0](0,0) = 1;
-           this->restriction[0](8,1) = 1;
-           this->restriction[0](11,3) = 1;
-           this->restriction[0](16,4) = 1;
-           this->restriction[0](20,2) = 1;
-           this->restriction[0](22,5) = 1;
-           this->restriction[0](25,7) = 1;
-           this->restriction[0](26,6) = 1;
-           this->restriction[1](1,1) = 1;
-           this->restriction[1](8,0) = 1;
-           this->restriction[1](9,2) = 1;
-           this->restriction[1](17,5) = 1;
-           this->restriction[1](20,3) = 1;
-           this->restriction[1](22,4) = 1;
-           this->restriction[1](23,6) = 1;
-           this->restriction[1](26,7) = 1;
-           this->restriction[2](2,2) = 1;
-           this->restriction[2](9,1) = 1;
-           this->restriction[2](10,3) = 1;
-           this->restriction[2](18,6) = 1;
-           this->restriction[2](20,0) = 1;
-           this->restriction[2](23,5) = 1;
-           this->restriction[2](24,7) = 1;
-           this->restriction[2](26,4) = 1;
-           this->restriction[3](3,3) = 1;
-           this->restriction[3](10,2) = 1;
-           this->restriction[3](11,0) = 1;
-           this->restriction[3](19,7) = 1;
-           this->restriction[3](20,1) = 1;
-           this->restriction[3](24,6) = 1;
-           this->restriction[3](25,4) = 1;
-           this->restriction[3](26,5) = 1;
-           this->restriction[4](4,4) = 1;
-           this->restriction[4](12,5) = 1;
-           this->restriction[4](15,7) = 1;
-           this->restriction[4](16,0) = 1;
-           this->restriction[4](21,6) = 1;
-           this->restriction[4](22,1) = 1;
-           this->restriction[4](25,3) = 1;
-           this->restriction[4](26,2) = 1;
-           this->restriction[5](5,5) = 1;
-           this->restriction[5](12,4) = 1;
-           this->restriction[5](13,6) = 1;
-           this->restriction[5](17,1) = 1;
-           this->restriction[5](21,7) = 1;
-           this->restriction[5](22,0) = 1;
-           this->restriction[5](23,2) = 1;
-           this->restriction[5](26,3) = 1;
-           this->restriction[6](6,6) = 1;
-           this->restriction[6](13,5) = 1;
-           this->restriction[6](14,7) = 1;
-           this->restriction[6](18,2) = 1;
-           this->restriction[6](21,4) = 1;
-           this->restriction[6](23,1) = 1;
-           this->restriction[6](24,3) = 1;
-           this->restriction[6](26,0) = 1;
-           this->restriction[7](7,7) = 1;
-           this->restriction[7](14,6) = 1;
-           this->restriction[7](15,4) = 1;
-           this->restriction[7](19,3) = 1;
-           this->restriction[7](21,5) = 1;
-           this->restriction[7](24,2) = 1;
-           this->restriction[7](25,0) = 1;
-           this->restriction[7](26,1) = 1;
-           break;
-         case 3:
-           this->restriction[0](0,0) = 1;
-           this->restriction[0](8,9) = 1;
-           this->restriction[0](14,15) = 1;
-           this->restriction[0](24,25) = 1;
-           this->restriction[0](32,35) = 1;
-           this->restriction[0](40,43) = 1;
-           this->restriction[0](52,55) = 1;
-           this->restriction[0](56,63) = 1;
-           this->restriction[1](1,1) = 1;
-           this->restriction[1](9,8) = 1;
-           this->restriction[1](10,11) = 1;
-           this->restriction[1](26,27) = 1;
-           this->restriction[1](33,34) = 1;
-           this->restriction[1](41,42) = 1;
-           this->restriction[1](44,47) = 1;
-           this->restriction[1](57,62) = 1;
-           this->restriction[2](2,2) = 1;
-           this->restriction[2](11,10) = 1;
-           this->restriction[2](13,12) = 1;
-           this->restriction[2](28,29) = 1;
-           this->restriction[2](35,32) = 1;
-           this->restriction[2](46,45) = 1;
-           this->restriction[2](49,50) = 1;
-           this->restriction[2](61,58) = 1;
-           this->restriction[3](3,3) = 1;
-           this->restriction[3](12,13) = 1;
-           this->restriction[3](15,14) = 1;
-           this->restriction[3](30,31) = 1;
-           this->restriction[3](34,33) = 1;
-           this->restriction[3](48,51) = 1;
-           this->restriction[3](54,53) = 1;
-           this->restriction[3](60,59) = 1;
-           this->restriction[4](4,4) = 1;
-           this->restriction[4](16,17) = 1;
-           this->restriction[4](22,23) = 1;
-           this->restriction[4](25,24) = 1;
-           this->restriction[4](36,39) = 1;
-           this->restriction[4](42,41) = 1;
-           this->restriction[4](53,54) = 1;
-           this->restriction[4](58,61) = 1;
-           this->restriction[5](5,5) = 1;
-           this->restriction[5](17,16) = 1;
-           this->restriction[5](18,19) = 1;
-           this->restriction[5](27,26) = 1;
-           this->restriction[5](37,38) = 1;
-           this->restriction[5](43,40) = 1;
-           this->restriction[5](45,46) = 1;
-           this->restriction[5](59,60) = 1;
-           this->restriction[6](6,6) = 1;
-           this->restriction[6](19,18) = 1;
-           this->restriction[6](21,20) = 1;
-           this->restriction[6](29,28) = 1;
-           this->restriction[6](39,36) = 1;
-           this->restriction[6](47,44) = 1;
-           this->restriction[6](51,48) = 1;
-           this->restriction[6](63,56) = 1;
-           this->restriction[7](7,7) = 1;
-           this->restriction[7](20,21) = 1;
-           this->restriction[7](23,22) = 1;
-           this->restriction[7](31,30) = 1;
-           this->restriction[7](38,37) = 1;
-           this->restriction[7](50,49) = 1;
-           this->restriction[7](55,52) = 1;
-           this->restriction[7](62,57) = 1;
-           break;
-         case 4:
-           this->restriction[0](0,0) = 1;
-           this->restriction[0](8,9) = 1;
-           this->restriction[0](9,1) = 1;
-           this->restriction[0](17,18) = 1;
-           this->restriction[0](18,3) = 1;
-           this->restriction[0](32,33) = 1;
-           this->restriction[0](33,4) = 1;
-           this->restriction[0](44,48) = 1;
-           this->restriction[0](45,12) = 1;
-           this->restriction[0](47,15) = 1;
-           this->restriction[0](48,2) = 1;
-           this->restriction[0](62,66) = 1;
-           this->restriction[0](63,36) = 1;
-           this->restriction[0](65,21) = 1;
-           this->restriction[0](66,5) = 1;
-           this->restriction[0](89,93) = 1;
-           this->restriction[0](90,30) = 1;
-           this->restriction[0](92,42) = 1;
-           this->restriction[0](93,7) = 1;
-           this->restriction[0](98,111) = 1;
-           this->restriction[0](99,75) = 1;
-           this->restriction[0](101,57) = 1;
-           this->restriction[0](102,24) = 1;
-           this->restriction[0](107,84) = 1;
-           this->restriction[0](108,39) = 1;
-           this->restriction[0](110,27) = 1;
-           this->restriction[0](111,6) = 1;
-           this->restriction[1](1,1) = 1;
-           this->restriction[1](9,0) = 1;
-           this->restriction[1](10,9) = 1;
-           this->restriction[1](11,12) = 1;
-           this->restriction[1](12,2) = 1;
-           this->restriction[1](35,36) = 1;
-           this->restriction[1](36,5) = 1;
-           this->restriction[1](45,18) = 1;
-           this->restriction[1](46,48) = 1;
-           this->restriction[1](48,3) = 1;
-           this->restriction[1](49,15) = 1;
-           this->restriction[1](63,33) = 1;
-           this->restriction[1](64,66) = 1;
-           this->restriction[1](66,4) = 1;
-           this->restriction[1](67,21) = 1;
-           this->restriction[1](71,75) = 1;
-           this->restriction[1](72,24) = 1;
-           this->restriction[1](74,39) = 1;
-           this->restriction[1](75,6) = 1;
-           this->restriction[1](99,93) = 1;
-           this->restriction[1](100,111) = 1;
-           this->restriction[1](102,30) = 1;
-           this->restriction[1](103,57) = 1;
-           this->restriction[1](108,42) = 1;
-           this->restriction[1](109,84) = 1;
-           this->restriction[1](111,7) = 1;
-           this->restriction[1](112,27) = 1;
-           this->restriction[2](2,2) = 1;
-           this->restriction[2](12,1) = 1;
-           this->restriction[2](13,12) = 1;
-           this->restriction[2](15,3) = 1;
-           this->restriction[2](16,15) = 1;
-           this->restriction[2](38,39) = 1;
-           this->restriction[2](39,6) = 1;
-           this->restriction[2](48,0) = 1;
-           this->restriction[2](49,9) = 1;
-           this->restriction[2](51,18) = 1;
-           this->restriction[2](52,48) = 1;
-           this->restriction[2](74,36) = 1;
-           this->restriction[2](75,5) = 1;
-           this->restriction[2](77,75) = 1;
-           this->restriction[2](78,24) = 1;
-           this->restriction[2](81,42) = 1;
-           this->restriction[2](82,84) = 1;
-           this->restriction[2](84,7) = 1;
-           this->restriction[2](85,27) = 1;
-           this->restriction[2](108,33) = 1;
-           this->restriction[2](109,66) = 1;
-           this->restriction[2](111,4) = 1;
-           this->restriction[2](112,21) = 1;
-           this->restriction[2](117,93) = 1;
-           this->restriction[2](118,111) = 1;
-           this->restriction[2](120,30) = 1;
-           this->restriction[2](121,57) = 1;
-           this->restriction[3](3,3) = 1;
-           this->restriction[3](14,15) = 1;
-           this->restriction[3](15,2) = 1;
-           this->restriction[3](18,0) = 1;
-           this->restriction[3](19,18) = 1;
-           this->restriction[3](41,42) = 1;
-           this->restriction[3](42,7) = 1;
-           this->restriction[3](47,9) = 1;
-           this->restriction[3](48,1) = 1;
-           this->restriction[3](50,48) = 1;
-           this->restriction[3](51,12) = 1;
-           this->restriction[3](80,84) = 1;
-           this->restriction[3](81,39) = 1;
-           this->restriction[3](83,27) = 1;
-           this->restriction[3](84,6) = 1;
-           this->restriction[3](92,33) = 1;
-           this->restriction[3](93,4) = 1;
-           this->restriction[3](95,93) = 1;
-           this->restriction[3](96,30) = 1;
-           this->restriction[3](107,66) = 1;
-           this->restriction[3](108,36) = 1;
-           this->restriction[3](110,21) = 1;
-           this->restriction[3](111,5) = 1;
-           this->restriction[3](116,111) = 1;
-           this->restriction[3](117,75) = 1;
-           this->restriction[3](119,57) = 1;
-           this->restriction[3](120,24) = 1;
-           this->restriction[4](4,4) = 1;
-           this->restriction[4](20,21) = 1;
-           this->restriction[4](21,5) = 1;
-           this->restriction[4](29,30) = 1;
-           this->restriction[4](30,7) = 1;
-           this->restriction[4](33,0) = 1;
-           this->restriction[4](34,33) = 1;
-           this->restriction[4](53,57) = 1;
-           this->restriction[4](54,24) = 1;
-           this->restriction[4](56,27) = 1;
-           this->restriction[4](57,6) = 1;
-           this->restriction[4](65,9) = 1;
-           this->restriction[4](66,1) = 1;
-           this->restriction[4](68,66) = 1;
-           this->restriction[4](69,36) = 1;
-           this->restriction[4](90,18) = 1;
-           this->restriction[4](91,93) = 1;
-           this->restriction[4](93,3) = 1;
-           this->restriction[4](94,42) = 1;
-           this->restriction[4](101,48) = 1;
-           this->restriction[4](102,12) = 1;
-           this->restriction[4](104,111) = 1;
-           this->restriction[4](105,75) = 1;
-           this->restriction[4](110,15) = 1;
-           this->restriction[4](111,2) = 1;
-           this->restriction[4](113,84) = 1;
-           this->restriction[4](114,39) = 1;
-           this->restriction[5](5,5) = 1;
-           this->restriction[5](21,4) = 1;
-           this->restriction[5](22,21) = 1;
-           this->restriction[5](23,24) = 1;
-           this->restriction[5](24,6) = 1;
-           this->restriction[5](36,1) = 1;
-           this->restriction[5](37,36) = 1;
-           this->restriction[5](54,30) = 1;
-           this->restriction[5](55,57) = 1;
-           this->restriction[5](57,7) = 1;
-           this->restriction[5](58,27) = 1;
-           this->restriction[5](66,0) = 1;
-           this->restriction[5](67,9) = 1;
-           this->restriction[5](69,33) = 1;
-           this->restriction[5](70,66) = 1;
-           this->restriction[5](72,12) = 1;
-           this->restriction[5](73,75) = 1;
-           this->restriction[5](75,2) = 1;
-           this->restriction[5](76,39) = 1;
-           this->restriction[5](102,18) = 1;
-           this->restriction[5](103,48) = 1;
-           this->restriction[5](105,93) = 1;
-           this->restriction[5](106,111) = 1;
-           this->restriction[5](111,3) = 1;
-           this->restriction[5](112,15) = 1;
-           this->restriction[5](114,42) = 1;
-           this->restriction[5](115,84) = 1;
-           this->restriction[6](6,6) = 1;
-           this->restriction[6](24,5) = 1;
-           this->restriction[6](25,24) = 1;
-           this->restriction[6](27,7) = 1;
-           this->restriction[6](28,27) = 1;
-           this->restriction[6](39,2) = 1;
-           this->restriction[6](40,39) = 1;
-           this->restriction[6](57,4) = 1;
-           this->restriction[6](58,21) = 1;
-           this->restriction[6](60,30) = 1;
-           this->restriction[6](61,57) = 1;
-           this->restriction[6](75,1) = 1;
-           this->restriction[6](76,36) = 1;
-           this->restriction[6](78,12) = 1;
-           this->restriction[6](79,75) = 1;
-           this->restriction[6](84,3) = 1;
-           this->restriction[6](85,15) = 1;
-           this->restriction[6](87,42) = 1;
-           this->restriction[6](88,84) = 1;
-           this->restriction[6](111,0) = 1;
-           this->restriction[6](112,9) = 1;
-           this->restriction[6](114,33) = 1;
-           this->restriction[6](115,66) = 1;
-           this->restriction[6](120,18) = 1;
-           this->restriction[6](121,48) = 1;
-           this->restriction[6](123,93) = 1;
-           this->restriction[6](124,111) = 1;
-           this->restriction[7](7,7) = 1;
-           this->restriction[7](26,27) = 1;
-           this->restriction[7](27,6) = 1;
-           this->restriction[7](30,4) = 1;
-           this->restriction[7](31,30) = 1;
-           this->restriction[7](42,3) = 1;
-           this->restriction[7](43,42) = 1;
-           this->restriction[7](56,21) = 1;
-           this->restriction[7](57,5) = 1;
-           this->restriction[7](59,57) = 1;
-           this->restriction[7](60,24) = 1;
-           this->restriction[7](83,15) = 1;
-           this->restriction[7](84,2) = 1;
-           this->restriction[7](86,84) = 1;
-           this->restriction[7](87,39) = 1;
-           this->restriction[7](93,0) = 1;
-           this->restriction[7](94,33) = 1;
-           this->restriction[7](96,18) = 1;
-           this->restriction[7](97,93) = 1;
-           this->restriction[7](110,9) = 1;
-           this->restriction[7](111,1) = 1;
-           this->restriction[7](113,66) = 1;
-           this->restriction[7](114,36) = 1;
-           this->restriction[7](119,48) = 1;
-           this->restriction[7](120,12) = 1;
-           this->restriction[7](122,111) = 1;
-           this->restriction[7](123,75) = 1;
-           break;
-         default:
-       {
-                                          // in case we don't
-                                          // have the matrices
-                                          // (yet), set them to
-                                          // impossible
-                                          // values. this does
-                                          // not prevent the use
-                                          // of this FE, but will
-                                          // prevent the use of
-                                          // these matrices
-         for (unsigned int i=0;
-              i<GeometryInfo<dim>::children_per_cell;
-              ++i)
-           this->restriction[i].reinit(0,0);
-       };
-       }
-      break;
-    };
+      {
+        switch (degree)
+          {
+            case 1:
+                  this->restriction[0](0,0) = 1;
+                  this->restriction[1](1,1) = 1;
+                  this->restriction[2](2,2) = 1;
+                  this->restriction[3](3,3) = 1;
+                  this->restriction[4](4,4) = 1;
+                  this->restriction[5](5,5) = 1;
+                  this->restriction[6](6,6) = 1;
+                  this->restriction[7](7,7) = 1;
+                  break;
+            case 2:
+                  this->restriction[0](0,0) = 1;
+                  this->restriction[0](8,1) = 1;
+                  this->restriction[0](11,3) = 1;
+                  this->restriction[0](16,4) = 1;
+                  this->restriction[0](20,2) = 1;
+                  this->restriction[0](22,5) = 1;
+                  this->restriction[0](25,7) = 1;
+                  this->restriction[0](26,6) = 1;
+                  this->restriction[1](1,1) = 1;
+                  this->restriction[1](8,0) = 1;
+                  this->restriction[1](9,2) = 1;
+                  this->restriction[1](17,5) = 1;
+                  this->restriction[1](20,3) = 1;
+                  this->restriction[1](22,4) = 1;
+                  this->restriction[1](23,6) = 1;
+                  this->restriction[1](26,7) = 1;
+                  this->restriction[2](2,2) = 1;
+                  this->restriction[2](9,1) = 1;
+                  this->restriction[2](10,3) = 1;
+                  this->restriction[2](18,6) = 1;
+                  this->restriction[2](20,0) = 1;
+                  this->restriction[2](23,5) = 1;
+                  this->restriction[2](24,7) = 1;
+                  this->restriction[2](26,4) = 1;
+                  this->restriction[3](3,3) = 1;
+                  this->restriction[3](10,2) = 1;
+                  this->restriction[3](11,0) = 1;
+                  this->restriction[3](19,7) = 1;
+                  this->restriction[3](20,1) = 1;
+                  this->restriction[3](24,6) = 1;
+                  this->restriction[3](25,4) = 1;
+                  this->restriction[3](26,5) = 1;
+                  this->restriction[4](4,4) = 1;
+                  this->restriction[4](12,5) = 1;
+                  this->restriction[4](15,7) = 1;
+                  this->restriction[4](16,0) = 1;
+                  this->restriction[4](21,6) = 1;
+                  this->restriction[4](22,1) = 1;
+                  this->restriction[4](25,3) = 1;
+                  this->restriction[4](26,2) = 1;
+                  this->restriction[5](5,5) = 1;
+                  this->restriction[5](12,4) = 1;
+                  this->restriction[5](13,6) = 1;
+                  this->restriction[5](17,1) = 1;
+                  this->restriction[5](21,7) = 1;
+                  this->restriction[5](22,0) = 1;
+                  this->restriction[5](23,2) = 1;
+                  this->restriction[5](26,3) = 1;
+                  this->restriction[6](6,6) = 1;
+                  this->restriction[6](13,5) = 1;
+                  this->restriction[6](14,7) = 1;
+                  this->restriction[6](18,2) = 1;
+                  this->restriction[6](21,4) = 1;
+                  this->restriction[6](23,1) = 1;
+                  this->restriction[6](24,3) = 1;
+                  this->restriction[6](26,0) = 1;
+                  this->restriction[7](7,7) = 1;
+                  this->restriction[7](14,6) = 1;
+                  this->restriction[7](15,4) = 1;
+                  this->restriction[7](19,3) = 1;
+                  this->restriction[7](21,5) = 1;
+                  this->restriction[7](24,2) = 1;
+                  this->restriction[7](25,0) = 1;
+                  this->restriction[7](26,1) = 1;
+                  break;
+            case 3:
+                  this->restriction[0](0,0) = 1;
+                  this->restriction[0](8,9) = 1;
+                  this->restriction[0](14,15) = 1;
+                  this->restriction[0](24,25) = 1;
+                  this->restriction[0](32,35) = 1;
+                  this->restriction[0](40,43) = 1;
+                  this->restriction[0](52,55) = 1;
+                  this->restriction[0](56,63) = 1;
+                  this->restriction[1](1,1) = 1;
+                  this->restriction[1](9,8) = 1;
+                  this->restriction[1](10,11) = 1;
+                  this->restriction[1](26,27) = 1;
+                  this->restriction[1](33,34) = 1;
+                  this->restriction[1](41,42) = 1;
+                  this->restriction[1](44,47) = 1;
+                  this->restriction[1](57,62) = 1;
+                  this->restriction[2](2,2) = 1;
+                  this->restriction[2](11,10) = 1;
+                  this->restriction[2](13,12) = 1;
+                  this->restriction[2](28,29) = 1;
+                  this->restriction[2](35,32) = 1;
+                  this->restriction[2](46,45) = 1;
+                  this->restriction[2](49,50) = 1;
+                  this->restriction[2](61,58) = 1;
+                  this->restriction[3](3,3) = 1;
+                  this->restriction[3](12,13) = 1;
+                  this->restriction[3](15,14) = 1;
+                  this->restriction[3](30,31) = 1;
+                  this->restriction[3](34,33) = 1;
+                  this->restriction[3](48,51) = 1;
+                  this->restriction[3](54,53) = 1;
+                  this->restriction[3](60,59) = 1;
+                  this->restriction[4](4,4) = 1;
+                  this->restriction[4](16,17) = 1;
+                  this->restriction[4](22,23) = 1;
+                  this->restriction[4](25,24) = 1;
+                  this->restriction[4](36,39) = 1;
+                  this->restriction[4](42,41) = 1;
+                  this->restriction[4](53,54) = 1;
+                  this->restriction[4](58,61) = 1;
+                  this->restriction[5](5,5) = 1;
+                  this->restriction[5](17,16) = 1;
+                  this->restriction[5](18,19) = 1;
+                  this->restriction[5](27,26) = 1;
+                  this->restriction[5](37,38) = 1;
+                  this->restriction[5](43,40) = 1;
+                  this->restriction[5](45,46) = 1;
+                  this->restriction[5](59,60) = 1;
+                  this->restriction[6](6,6) = 1;
+                  this->restriction[6](19,18) = 1;
+                  this->restriction[6](21,20) = 1;
+                  this->restriction[6](29,28) = 1;
+                  this->restriction[6](39,36) = 1;
+                  this->restriction[6](47,44) = 1;
+                  this->restriction[6](51,48) = 1;
+                  this->restriction[6](63,56) = 1;
+                  this->restriction[7](7,7) = 1;
+                  this->restriction[7](20,21) = 1;
+                  this->restriction[7](23,22) = 1;
+                  this->restriction[7](31,30) = 1;
+                  this->restriction[7](38,37) = 1;
+                  this->restriction[7](50,49) = 1;
+                  this->restriction[7](55,52) = 1;
+                  this->restriction[7](62,57) = 1;
+                  break;
+            case 4:
+                  this->restriction[0](0,0) = 1;
+                  this->restriction[0](8,9) = 1;
+                  this->restriction[0](9,1) = 1;
+                  this->restriction[0](17,18) = 1;
+                  this->restriction[0](18,3) = 1;
+                  this->restriction[0](32,33) = 1;
+                  this->restriction[0](33,4) = 1;
+                  this->restriction[0](44,48) = 1;
+                  this->restriction[0](45,12) = 1;
+                  this->restriction[0](47,15) = 1;
+                  this->restriction[0](48,2) = 1;
+                  this->restriction[0](62,66) = 1;
+                  this->restriction[0](63,36) = 1;
+                  this->restriction[0](65,21) = 1;
+                  this->restriction[0](66,5) = 1;
+                  this->restriction[0](89,93) = 1;
+                  this->restriction[0](90,30) = 1;
+                  this->restriction[0](92,42) = 1;
+                  this->restriction[0](93,7) = 1;
+                  this->restriction[0](98,111) = 1;
+                  this->restriction[0](99,75) = 1;
+                  this->restriction[0](101,57) = 1;
+                  this->restriction[0](102,24) = 1;
+                  this->restriction[0](107,84) = 1;
+                  this->restriction[0](108,39) = 1;
+                  this->restriction[0](110,27) = 1;
+                  this->restriction[0](111,6) = 1;
+                  this->restriction[1](1,1) = 1;
+                  this->restriction[1](9,0) = 1;
+                  this->restriction[1](10,9) = 1;
+                  this->restriction[1](11,12) = 1;
+                  this->restriction[1](12,2) = 1;
+                  this->restriction[1](35,36) = 1;
+                  this->restriction[1](36,5) = 1;
+                  this->restriction[1](45,18) = 1;
+                  this->restriction[1](46,48) = 1;
+                  this->restriction[1](48,3) = 1;
+                  this->restriction[1](49,15) = 1;
+                  this->restriction[1](63,33) = 1;
+                  this->restriction[1](64,66) = 1;
+                  this->restriction[1](66,4) = 1;
+                  this->restriction[1](67,21) = 1;
+                  this->restriction[1](71,75) = 1;
+                  this->restriction[1](72,24) = 1;
+                  this->restriction[1](74,39) = 1;
+                  this->restriction[1](75,6) = 1;
+                  this->restriction[1](99,93) = 1;
+                  this->restriction[1](100,111) = 1;
+                  this->restriction[1](102,30) = 1;
+                  this->restriction[1](103,57) = 1;
+                  this->restriction[1](108,42) = 1;
+                  this->restriction[1](109,84) = 1;
+                  this->restriction[1](111,7) = 1;
+                  this->restriction[1](112,27) = 1;
+                  this->restriction[2](2,2) = 1;
+                  this->restriction[2](12,1) = 1;
+                  this->restriction[2](13,12) = 1;
+                  this->restriction[2](15,3) = 1;
+                  this->restriction[2](16,15) = 1;
+                  this->restriction[2](38,39) = 1;
+                  this->restriction[2](39,6) = 1;
+                  this->restriction[2](48,0) = 1;
+                  this->restriction[2](49,9) = 1;
+                  this->restriction[2](51,18) = 1;
+                  this->restriction[2](52,48) = 1;
+                  this->restriction[2](74,36) = 1;
+                  this->restriction[2](75,5) = 1;
+                  this->restriction[2](77,75) = 1;
+                  this->restriction[2](78,24) = 1;
+                  this->restriction[2](81,42) = 1;
+                  this->restriction[2](82,84) = 1;
+                  this->restriction[2](84,7) = 1;
+                  this->restriction[2](85,27) = 1;
+                  this->restriction[2](108,33) = 1;
+                  this->restriction[2](109,66) = 1;
+                  this->restriction[2](111,4) = 1;
+                  this->restriction[2](112,21) = 1;
+                  this->restriction[2](117,93) = 1;
+                  this->restriction[2](118,111) = 1;
+                  this->restriction[2](120,30) = 1;
+                  this->restriction[2](121,57) = 1;
+                  this->restriction[3](3,3) = 1;
+                  this->restriction[3](14,15) = 1;
+                  this->restriction[3](15,2) = 1;
+                  this->restriction[3](18,0) = 1;
+                  this->restriction[3](19,18) = 1;
+                  this->restriction[3](41,42) = 1;
+                  this->restriction[3](42,7) = 1;
+                  this->restriction[3](47,9) = 1;
+                  this->restriction[3](48,1) = 1;
+                  this->restriction[3](50,48) = 1;
+                  this->restriction[3](51,12) = 1;
+                  this->restriction[3](80,84) = 1;
+                  this->restriction[3](81,39) = 1;
+                  this->restriction[3](83,27) = 1;
+                  this->restriction[3](84,6) = 1;
+                  this->restriction[3](92,33) = 1;
+                  this->restriction[3](93,4) = 1;
+                  this->restriction[3](95,93) = 1;
+                  this->restriction[3](96,30) = 1;
+                  this->restriction[3](107,66) = 1;
+                  this->restriction[3](108,36) = 1;
+                  this->restriction[3](110,21) = 1;
+                  this->restriction[3](111,5) = 1;
+                  this->restriction[3](116,111) = 1;
+                  this->restriction[3](117,75) = 1;
+                  this->restriction[3](119,57) = 1;
+                  this->restriction[3](120,24) = 1;
+                  this->restriction[4](4,4) = 1;
+                  this->restriction[4](20,21) = 1;
+                  this->restriction[4](21,5) = 1;
+                  this->restriction[4](29,30) = 1;
+                  this->restriction[4](30,7) = 1;
+                  this->restriction[4](33,0) = 1;
+                  this->restriction[4](34,33) = 1;
+                  this->restriction[4](53,57) = 1;
+                  this->restriction[4](54,24) = 1;
+                  this->restriction[4](56,27) = 1;
+                  this->restriction[4](57,6) = 1;
+                  this->restriction[4](65,9) = 1;
+                  this->restriction[4](66,1) = 1;
+                  this->restriction[4](68,66) = 1;
+                  this->restriction[4](69,36) = 1;
+                  this->restriction[4](90,18) = 1;
+                  this->restriction[4](91,93) = 1;
+                  this->restriction[4](93,3) = 1;
+                  this->restriction[4](94,42) = 1;
+                  this->restriction[4](101,48) = 1;
+                  this->restriction[4](102,12) = 1;
+                  this->restriction[4](104,111) = 1;
+                  this->restriction[4](105,75) = 1;
+                  this->restriction[4](110,15) = 1;
+                  this->restriction[4](111,2) = 1;
+                  this->restriction[4](113,84) = 1;
+                  this->restriction[4](114,39) = 1;
+                  this->restriction[5](5,5) = 1;
+                  this->restriction[5](21,4) = 1;
+                  this->restriction[5](22,21) = 1;
+                  this->restriction[5](23,24) = 1;
+                  this->restriction[5](24,6) = 1;
+                  this->restriction[5](36,1) = 1;
+                  this->restriction[5](37,36) = 1;
+                  this->restriction[5](54,30) = 1;
+                  this->restriction[5](55,57) = 1;
+                  this->restriction[5](57,7) = 1;
+                  this->restriction[5](58,27) = 1;
+                  this->restriction[5](66,0) = 1;
+                  this->restriction[5](67,9) = 1;
+                  this->restriction[5](69,33) = 1;
+                  this->restriction[5](70,66) = 1;
+                  this->restriction[5](72,12) = 1;
+                  this->restriction[5](73,75) = 1;
+                  this->restriction[5](75,2) = 1;
+                  this->restriction[5](76,39) = 1;
+                  this->restriction[5](102,18) = 1;
+                  this->restriction[5](103,48) = 1;
+                  this->restriction[5](105,93) = 1;
+                  this->restriction[5](106,111) = 1;
+                  this->restriction[5](111,3) = 1;
+                  this->restriction[5](112,15) = 1;
+                  this->restriction[5](114,42) = 1;
+                  this->restriction[5](115,84) = 1;
+                  this->restriction[6](6,6) = 1;
+                  this->restriction[6](24,5) = 1;
+                  this->restriction[6](25,24) = 1;
+                  this->restriction[6](27,7) = 1;
+                  this->restriction[6](28,27) = 1;
+                  this->restriction[6](39,2) = 1;
+                  this->restriction[6](40,39) = 1;
+                  this->restriction[6](57,4) = 1;
+                  this->restriction[6](58,21) = 1;
+                  this->restriction[6](60,30) = 1;
+                  this->restriction[6](61,57) = 1;
+                  this->restriction[6](75,1) = 1;
+                  this->restriction[6](76,36) = 1;
+                  this->restriction[6](78,12) = 1;
+                  this->restriction[6](79,75) = 1;
+                  this->restriction[6](84,3) = 1;
+                  this->restriction[6](85,15) = 1;
+                  this->restriction[6](87,42) = 1;
+                  this->restriction[6](88,84) = 1;
+                  this->restriction[6](111,0) = 1;
+                  this->restriction[6](112,9) = 1;
+                  this->restriction[6](114,33) = 1;
+                  this->restriction[6](115,66) = 1;
+                  this->restriction[6](120,18) = 1;
+                  this->restriction[6](121,48) = 1;
+                  this->restriction[6](123,93) = 1;
+                  this->restriction[6](124,111) = 1;
+                  this->restriction[7](7,7) = 1;
+                  this->restriction[7](26,27) = 1;
+                  this->restriction[7](27,6) = 1;
+                  this->restriction[7](30,4) = 1;
+                  this->restriction[7](31,30) = 1;
+                  this->restriction[7](42,3) = 1;
+                  this->restriction[7](43,42) = 1;
+                  this->restriction[7](56,21) = 1;
+                  this->restriction[7](57,5) = 1;
+                  this->restriction[7](59,57) = 1;
+                  this->restriction[7](60,24) = 1;
+                  this->restriction[7](83,15) = 1;
+                  this->restriction[7](84,2) = 1;
+                  this->restriction[7](86,84) = 1;
+                  this->restriction[7](87,39) = 1;
+                  this->restriction[7](93,0) = 1;
+                  this->restriction[7](94,33) = 1;
+                  this->restriction[7](96,18) = 1;
+                  this->restriction[7](97,93) = 1;
+                  this->restriction[7](110,9) = 1;
+                  this->restriction[7](111,1) = 1;
+                  this->restriction[7](113,66) = 1;
+                  this->restriction[7](114,36) = 1;
+                  this->restriction[7](119,48) = 1;
+                  this->restriction[7](120,12) = 1;
+                  this->restriction[7](122,111) = 1;
+                  this->restriction[7](123,75) = 1;
+                  break;
+            default:
+            {
+                                               // in case we don't
+                                               // have the matrices
+                                               // (yet), reset them to
+                                               // zero size. this does
+                                               // not prevent the use
+                                               // of this FE, but will
+                                               // prevent the use of
+                                               // these matrices
+              for (unsigned int i=0;
+                   i<GeometryInfo<dim>::children_per_cell;
+                   ++i)
+                this->restriction[i].reinit(0,0);
+            };
+          }
+        break;
+      };
       
       default:
-       Assert (false,ExcNotImplemented());
+            Assert (false, ExcNotImplemented());
     }
 
                                   // finally fill in support points
index a890debab6e3001ab22af51fd34ddfb5b644db3f..0579ba33a5627761ef0d0174eec3eb951d03b397 100644 (file)
@@ -1215,6 +1215,9 @@ FESystem<dim>::build_face_tables()
 template <int dim>
 void FESystem<dim>::build_interface_constraints () 
 {
+  this->interface_constraints.
+    TableBase<2,double>::reinit (this->interface_constraints_size());
+  
                                   // the layout of the constraints
                                   // matrix is described in the
                                   // FiniteElement class. you may
@@ -1482,15 +1485,17 @@ void FESystem<dim>::initialize ()
        do_prolongation = false;
     }
   
-                                  // if we encountered void matrices,
-                                  // disable them for the composite
-                                  // element as well
-  if (!do_restriction)
+                                  // if we did not encounter void
+                                  // matrices, initialize the
+                                  // respective matrix sizes
+  if (do_restriction)
     for (unsigned int i=0;i<GeometryInfo<dim>::children_per_cell;++i)
-      this->restriction[i].reinit(0,0);
-  if (!do_prolongation)
+      this->restriction[i].reinit(this->dofs_per_cell,
+                                  this->dofs_per_cell);
+  if (do_prolongation)
     for (unsigned int i=0;i<GeometryInfo<dim>::children_per_cell;++i)
-      this->prolongation[i].reinit(0,0);
+      this->prolongation[i].reinit(this->dofs_per_cell,
+                                   this->dofs_per_cell);
        
                                   // distribute the matrices of the
                                   // base finite elements to the

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.