implemented codim 1 index set in 3D using side vectors

Note: There was a bug in UG/ug/gm/algebra.c, i.e. you have to update
UG as well.

[[Imported from SVN: r2951]]

grid/uggrid.hh

@@ -289,12 +289,10 @@ namespace Dune {
       {
         return size(codim,cube);
       }
-#if (dim==3)
       if (codim==1)
       {
         return size(1,cube)+size(1,simplex);
       }
-#endif
 
       std::cout << "dim=" << dim << " codim=" << codim << std::endl;
       DUNE_THROW(NotImplemented, "not implemented");

grid/uggrid/ugfunctions.hh

@@ -241,6 +241,16 @@ namespace Dune {
       return theElement->ge.levelIndex;
     }
 
+    //! Gets the level index of a UG sidevector
+    static int& levelIndex(UGVectorType<2>::T* theVector) {
+      return reinterpret_cast<int&>(theVector->index);
+    }
+
+    //! Gets the level index of a UG sidevector
+    static const int& levelIndex(const UGVectorType<2>::T* theVector) {
+      return reinterpret_cast<const int&>(theVector->index);
+    }
+
     //! Gets the level index of a UG edge
     static int& levelIndex(TargetType<1,2>::T* theEdge) {
       return theEdge->levelIndex;
@@ -275,6 +285,16 @@ namespace Dune {
       return theElement->ge.leafIndex;
     }
 
+    //! Gets the level index of a UG sidevector
+    static int& leafIndex(UGVectorType<2>::T* theVector) {
+      return reinterpret_cast<int &>(theVector->skip);
+    }
+
+    //! Gets the level index of a UG sidevector
+    static const int& leafIndex(const UGVectorType<2>::T* theVector) {
+      return reinterpret_cast<const int &>(theVector->skip);
+    }
+
     //! Gets the leaf index of a UG edge
     static int& leafIndex(TargetType<1,2>::T* theEdge) {
       return theEdge->leafIndex;
@@ -351,6 +371,12 @@ namespace Dune {
       return UG2d::GetEdge(nodei,nodej);
     }
 
+    //! access side vector from element (this is just a dummy to compile code also in 2d)
+    static UGVectorType<2>::T* SideVector (TargetType<0,2>::T* theElement, int i)
+    {
+      return NULL;
+    }
+
     //! \todo Please doc me!
     static TargetType<0,2>::T* EFather(TargetType<0,2>::T* theElement) {
       using UG2d::ELEMENT;

grid/uggrid/ugfunctions3d.hh

@@ -243,6 +243,16 @@ namespace Dune {
       return theElement->ge.levelIndex;
     }
 
+    //! Gets the level index of a UG sidevector
+    static int& levelIndex(UGVectorType<3>::T* theVector) {
+      return reinterpret_cast<int&>(theVector->index);
+    }
+
+    //! Gets the level index of a UG sidevector
+    static const int& levelIndex(const UGVectorType<3>::T* theVector) {
+      return reinterpret_cast<const int&>(theVector->index);
+    }
+
     //! Gets the level index of a UG edge
     static int& levelIndex(TargetType<2,3>::T* theEdge) {
       return theEdge->levelIndex;
@@ -277,6 +287,16 @@ namespace Dune {
       return theElement->ge.leafIndex;
     }
 
+    //! Gets the level index of a UG sidevector
+    static int& leafIndex(UGVectorType<3>::T* theVector) {
+      return reinterpret_cast<int &>(theVector->skip);
+    }
+
+    //! Gets the level index of a UG sidevector
+    static const int& leafIndex(const UGVectorType<3>::T* theVector) {
+      return reinterpret_cast<const int &>(theVector->skip);
+    }
+
     //! Gets the leaf index of a UG edge
     static int& leafIndex(TargetType<2,3>::T* theEdge) {
       return theEdge->leafIndex;
@@ -353,6 +373,14 @@ namespace Dune {
       return UG3d::GetEdge(nodei,nodej);
     }
 
+    //! access side vector from element
+    static UGVectorType<3>::T* SideVector (TargetType<0,3>::T* theElement, int i)
+    {
+      using UG3d::VECTOR;
+      using UG3d::svector_offset;
+      return SVECTOR(theElement,i);
+    }
+
     //! \todo Please doc me!
     static TargetType<0,3>::T* EFather(TargetType<0,3>::T* theElement) {
       using UG3d::ELEMENT;

grid/uggrid/uggrid.cc

@@ -416,7 +416,7 @@ inline int Dune::UGGrid < dim, dimworld >::size (int level, int codim) const
   {
     return this->levelIndexSet(level).size(0,simplex)+this->levelIndexSet(level).size(0,cube)
            +this->levelIndexSet(level).size(0,pyramid)+this->levelIndexSet(level).size(0,prism);
-  } else
+  }
   if(codim == dim)
   {
     return this->levelIndexSet(level).size(dim,cube);
@@ -425,12 +425,10 @@ inline int Dune::UGGrid < dim, dimworld >::size (int level, int codim) const
   {
     return this->levelIndexSet(level).size(dim-1,cube);
   }
-#if (dim==3)
   if (codim == 1)
   {
     return this->levelIndexSet(level).size(1,cube)+this->levelIndexSet(level).size(1,simplex);
   }
-#endif
   DUNE_THROW(GridError, "UGGrid<" << dim << ", " << dimworld
                                   << ">::size(int level, int codim) is only implemented"
                                   << " for codim==0 and codim==dim!");

grid/uggrid/uggridentity.cc

@@ -140,6 +140,30 @@ inline int UGGridEntity<0, dim, GridImp>::renumberVertex(int i) const {
 
 }
 template <int dim, class GridImp>
+inline int UGGridEntity<0, dim, GridImp>::renumberFace(int i) const {
+
+  if (geometry().type()==cube) {
+
+    // Dune numbers the vertices of a hexahedron and quadrilaterals differently than UG.
+    // The following two lines do the transformation
+    // The renumbering scheme is {0,1,3,2} for quadrilaterals, therefore, the
+    // following code works for 2d and 3d.
+    const int renumbering[6] = {4, 2, 1, 3, 0, 5};
+    return renumbering[i];
+
+  }
+  if (geometry().type()==simplex) {
+
+    // Dune numbers the vertices of a hexahedron and quadrilaterals differently than UG.
+    // The following two lines do the transformation
+    // The renumbering scheme is {0,1,3,2} for quadrilaterals, therefore, the
+    // following code works for 2d and 3d.
+    const int renumbering[4] = {1, 2, 0, 3};
+    return renumbering[i];
+  }
+  return i;
+}
+template <int dim, class GridImp>
 template <int cc>
 inline int UGGridEntity<0, dim, GridImp>::subIndex(int i) const
 {
@@ -155,9 +179,10 @@ inline int UGGridEntity<0, dim, GridImp>::subIndex(int i) const
     int b=ReferenceElements<double,dim>::general(geometry().type()).subentity(i,dim-1,1,dim);
     return UG_NS<dim>::levelIndex(UG_NS<dim>::GetEdge(UG_NS<dim>::Corner(target_,renumberVertex(a)),UG_NS<dim>::Corner(target_,renumberVertex(b))));
   }
-#if (dim==3)
-  // something for faces here
-#endif
+  if (cc==1)
+  {
+    return UG_NS<dim>::levelIndex(UG_NS<dim>::SideVector(target_,renumberFace(i)));
+  }
 
   DUNE_THROW(GridError, "UGGrid<" << dim << "," << dim << ">::subIndex isn't implemented for cc==" << cc );
 }
@@ -178,9 +203,10 @@ inline int UGGridEntity<0, dim, GridImp>::subLeafIndex(int i) const
     int b=ReferenceElements<double,dim>::general(geometry().type()).subentity(i,dim-1,1,dim);
     return UG_NS<dim>::leafIndex(UG_NS<dim>::GetEdge(UG_NS<dim>::Corner(target_,renumberVertex(a)),UG_NS<dim>::Corner(target_,renumberVertex(b))));
   }
-#if (dim==3)
-  // something for faces here
-#endif
+  if (cc==1)
+  {
+    return UG_NS<dim>::leafIndex(UG_NS<dim>::SideVector(target_,renumberFace(i)));
+  }
 
   DUNE_THROW(GridError, "UGGrid<" << dim << "," << dim << ">::subLeafIndex isn't implemented for cc==" << cc );
 }

grid/uggrid/uggridentity.hh

@@ -345,6 +345,7 @@ namespace Dune {
 
     //private:
     int renumberVertex(int i) const;
+    int renumberFace(int i) const;
 
     void setToTarget(typename TargetType<0,dim>::T* target, int level);
 

grid/uggrid/uggridgeometry.cc

@@ -96,7 +96,6 @@ struct UGGridGeometryPositionAccess<2,2>
 template< int mydim, int coorddim, class GridImp>
 inline GeometryType UGGridGeometry<mydim,coorddim,GridImp>::type() const
 {
-
   switch (mydim)
   {
   case 0 : return vertex;

grid/uggrid/uggridgeometry.hh

@@ -199,7 +199,7 @@ namespace Dune {
 
     /** \brief Default constructor */
     UGGridGeometry()
-    {}
+    {elementType_=simplex;}
 
     //! return the element type identifier (triangle or quadrilateral)
     GeometryType type () const {return elementType_;}

grid/uggrid/uggridindexsets.hh

@@ -70,13 +70,15 @@ namespace Dune {
         return numVertices_;
       if (codim==dim-1)
         return numEdges_;
+      if (codim==1)
+        return numTriFaces_+numQuadFaces_;
+      DUNE_THROW(NotImplemented, "wrong codim!");
     }
 
     //! get number of entities of given codim, type and on this level
     int size (int codim, GeometryType type) const
     {
       if (codim==0) {
-
         switch (type) {
         case simplex :
           return numSimplices_;
@@ -89,12 +91,24 @@ namespace Dune {
         default :
           return 0;
         }
-      } else if (codim==dim) {
+      }
+      if (codim==dim) {
         return numVertices_;
-      } else if (codim==dim-1) {
+      }
+      if (codim==dim-1) {
         return numEdges_;
-      } else
-        DUNE_THROW(NotImplemented, "Not yet implemented for this codim!");
+      }
+      if (codim==1) {
+        switch (type) {
+        case simplex :
+          return numTriFaces_;
+        case cube :
+          return numQuadFaces_;
+        default :
+          return 0;
+        }
+      }
+      DUNE_THROW(NotImplemented, "Wrong codim!");
     }
 
     /** \brief Deliver all geometry types used in this grid */
@@ -132,6 +146,30 @@ namespace Dune {
         return i;
     }
 
+    int renumberFace(GeometryType gt, int i) const
+    {
+
+      if (gt==cube) {
+
+        // Dune numbers the vertices of a hexahedron and quadrilaterals differently than UG.
+        // The following two lines do the transformation
+        // The renumbering scheme is {0,1,3,2} for quadrilaterals, therefore, the
+        // following code works for 2d and 3d.
+        const int renumbering[6] = {4, 2, 1, 3, 0, 5};
+        return renumbering[i];
+      }
+      if (gt==simplex) {
+
+        // Dune numbers the vertices of a hexahedron and quadrilaterals differently than UG.
+        // The following two lines do the transformation
+        // The renumbering scheme is {0,1,3,2} for quadrilaterals, therefore, the
+        // following code works for 2d and 3d.
+        const int renumbering[4] = {1, 2, 0, 3};
+        return renumbering[i];
+      }
+      return i;
+    }
+
     void update(const GridImp& grid, int level) {
 
       // Commit the index set to a specific level of a specific grid
@@ -146,19 +184,25 @@ namespace Dune {
       typename GridImp::Traits::template Codim<0>::LevelIterator eEndIt = grid_->template lend<0>(level_);
 
       for (; eIt!=eEndIt; ++eIt) {
+        typename TargetType<0,dim>::T* target_ = grid_->template getRealEntity<0>(*eIt).target_;
         // codim dim-1
         for (int i=0; i<eIt->template count<dim-1>(); i++)
         {
           GeometryType gt = eIt->geometry().type();
-          typename TargetType<0,dim>::T* target_ = grid_->template getRealEntity<0>(*eIt).target_;
           int a=ReferenceElements<double,dim>::general(gt).subentity(i,dim-1,0,dim);
           int b=ReferenceElements<double,dim>::general(gt).subentity(i,dim-1,1,dim);
           int& index = UG_NS<dim>::levelIndex(UG_NS<dim>::GetEdge(UG_NS<dim>::Corner(target_,renumberVertex(gt,a)),UG_NS<dim>::Corner(target_,renumberVertex(gt,b))));
           index = -1;
         }
         // codim 1 (faces): todo
-#if (dim==3)
-#endif
+        if (dim==3)
+          for (int i=0; i<eIt->template count<1>(); i++)
+          {
+            GeometryType gt = eIt->geometry().type();
+            int& index = UG_NS<dim>::levelIndex(UG_NS<dim>::SideVector(target_,renumberFace(gt,i)));
+            index = -1;
+          }
+
       }
 
       // ///////////////////////////////
@@ -169,6 +213,8 @@ namespace Dune {
       numPrisms_    = 0;
       numCubes_     = 0;
       numEdges_     = 0;
+      numTriFaces_  = 0;
+      numQuadFaces_ = 0;
 
       eIt    = grid_->template lbegin<0>(level_);
       eEndIt = grid_->template lend<0>(level_);
@@ -194,11 +240,12 @@ namespace Dune {
                                                           << ", which should never occur in a UGGrid!");
         }
 
+        typename TargetType<0,dim>::T* target_ = grid_->template getRealEntity<0>(*eIt).target_;
+
         // codim dim-1 (edges)
         for (int i=0; i<eIt->template count<dim-1>(); i++)
         {
           GeometryType gt = eIt->geometry().type();
-          typename TargetType<0,dim>::T* target_ = grid_->template getRealEntity<0>(*eIt).target_;
           int a=ReferenceElements<double,dim>::general(gt).subentity(i,dim-1,0,dim);
           int b=ReferenceElements<double,dim>::general(gt).subentity(i,dim-1,1,dim);
           int& index = UG_NS<dim>::levelIndex(UG_NS<dim>::GetEdge(UG_NS<dim>::Corner(target_,renumberVertex(gt,a)),UG_NS<dim>::Corner(target_,renumberVertex(gt,b))));
@@ -206,9 +253,24 @@ namespace Dune {
         }
 
         // codim 1 (faces): todo
-#if (dim==3)
-#endif
-
+        if (dim==3)
+          for (int i=0; i<eIt->template count<1>(); i++)
+          {
+            GeometryType gt = eIt->geometry().type();
+            int& index = UG_NS<dim>::levelIndex(UG_NS<dim>::SideVector(target_,renumberFace(gt,i)));
+            if (index<0)                       // not visited yet
+              switch (ReferenceElements<double,dim>::general(gt).type(i,1))
+              {
+              case simplex :
+                index = numTriFaces_++;
+                break;
+              case cube :
+                index = numQuadFaces_++;
+                break;
+              default :
+                DUNE_THROW(GridError, "wrong geometry type in face");
+              }
+          }
       }
 
       // Update the list of geometry types present
@@ -225,6 +287,15 @@ namespace Dune {
       myTypes_[dim-1].resize(0);
       myTypes_[dim-1].push_back(cube);
 
+      if (dim==3)
+      {
+        myTypes_[1].resize(0);
+        if (numTriFaces_ > 0)
+          myTypes_[1].push_back(simplex);
+        if (numQuadFaces_ > 0)
+          myTypes_[1].push_back(cube);
+      }
+
       // //////////////////////////////
       //   Init the vertex indices
       // //////////////////////////////
@@ -249,10 +320,8 @@ namespace Dune {
     int numCubes_;
     int numVertices_;
     int numEdges_;
-#if (dim==3)
     int numTriFaces_;
     int numQuadFaces_;
-#endif
 
     std::vector<GeometryType> myTypes_[dim+1];
   };
@@ -305,7 +374,6 @@ namespace Dune {
     int size (int codim, GeometryType type) const
     {
       if (codim==0) {
-
         switch (type) {
         case simplex :
           return numSimplices_;
@@ -318,12 +386,24 @@ namespace Dune {
         default :
           return 0;
         }
-      } else if (codim==dim) {
+      }
+      if (codim==dim) {
         return numVertices_;
-      } else if (codim==dim-1) {
+      }
+      if (codim==dim-1) {
         return numEdges_;
-      } else
-        DUNE_THROW(NotImplemented, "UGGridLeafIndexSet::size(codim,type) for codim neither 0 nor dim");
+      }
+      if (codim==1) {
+        switch (type) {
+        case simplex :
+          return numTriFaces_;
+        case cube :
+          return numQuadFaces_;
+        default :
+          return 0;
+        }
+      }
+      DUNE_THROW(NotImplemented, "Wrong codim!");
     }
 
     /** deliver all geometry types used in this grid */
@@ -361,6 +441,29 @@ namespace Dune {
         return i;
     }
 
+    int renumberFace(GeometryType gt, int i) const
+    {
+
+      if (gt==cube) {
+
+        // Dune numbers the vertices of a hexahedron and quadrilaterals differently than UG.
+        // The following two lines do the transformation
+        // The renumbering scheme is {0,1,3,2} for quadrilaterals, therefore, the
+        // following code works for 2d and 3d.
+        const int renumbering[6] = {4, 2, 1, 3, 0, 5};
+        return renumbering[i];
+      }
+      if (gt==simplex) {
+
+        // Dune numbers the vertices of a hexahedron and quadrilaterals differently than UG.
+        // The following two lines do the transformation
+        // The renumbering scheme is {0,1,3,2} for quadrilaterals, therefore, the
+        // following code works for 2d and 3d.
+        const int renumbering[4] = {1, 2, 0, 3};
+        return renumbering[i];
+      }
+      return i;
+    }
     void update() {
 
       // ///////////////////////////////////////////////////////////
@@ -376,19 +479,24 @@ namespace Dune {
 
       for (; eIt!=eEndIt; ++eIt)
       {
+        typename TargetType<0,dim>::T* target_ = grid_.template getRealEntity<0>(*eIt).target_;
         // codim dim-1
         for (int i=0; i<eIt->template count<dim-1>(); i++)
         {
           GeometryType gt = eIt->geometry().type();
-          typename TargetType<0,dim>::T* target_ = grid_.template getRealEntity<0>(*eIt).target_;
           int a=ReferenceElements<double,dim>::general(gt).subentity(i,dim-1,0,dim);
           int b=ReferenceElements<double,dim>::general(gt).subentity(i,dim-1,1,dim);
           int& index = UG_NS<dim>::leafIndex(UG_NS<dim>::GetEdge(UG_NS<dim>::Corner(target_,renumberVertex(gt,a)),UG_NS<dim>::Corner(target_,renumberVertex(gt,b))));
           index = -1;
         }
         // codim 1 (faces): todo
-#if (dim==3)
-#endif
+        if (dim==3)
+          for (int i=0; i<eIt->template count<1>(); i++)
+          {
+            GeometryType gt = eIt->geometry().type();
+            int& index = UG_NS<dim>::leafIndex(UG_NS<dim>::SideVector(target_,renumberFace(gt,i)));
+            index = -1;
+          }
       }
       // ///////////////////////////////
       //   Init the element indices
@@ -398,6 +506,8 @@ namespace Dune {
       numPrisms_    = 0;
       numCubes_     = 0;
       numEdges_     = 0;
+      numTriFaces_  = 0;
+      numQuadFaces_ = 0;
 
       eIt    = grid_.template leafbegin<0>();
       eEndIt = grid_.template leafend<0>();
@@ -422,11 +532,12 @@ namespace Dune {
                                                           << ", which should never occur in a UGGrid!");
         }
 
+        typename TargetType<0,dim>::T* target_ = grid_.template getRealEntity<0>(*eIt).target_;
+
         // codim dim-1 (edges)
         for (int i=0; i<eIt->template count<dim-1>(); i++)
         {
           GeometryType gt = eIt->geometry().type();
-          typename TargetType<0,dim>::T* target_ = grid_.template getRealEntity<0>(*eIt).target_;
           int a=ReferenceElements<double,dim>::general(gt).subentity(i,dim-1,0,dim);
           int b=ReferenceElements<double,dim>::general(gt).subentity(i,dim-1,1,dim);
           int& index = UG_NS<dim>::leafIndex(UG_NS<dim>::GetEdge(UG_NS<dim>::Corner(target_,renumberVertex(gt,a)),UG_NS<dim>::Corner(target_,renumberVertex(gt,b))));
@@ -434,8 +545,24 @@ namespace Dune {
         }
 
         // codim 1 (faces): todo
-#if (dim==3)
-#endif
+        if (dim==3)
+          for (int i=0; i<eIt->template count<1>(); i++)
+          {
+            GeometryType gt = eIt->geometry().type();
+            int& index = UG_NS<dim>::leafIndex(UG_NS<dim>::SideVector(target_,renumberFace(gt,i)));
+            if (index<0)                       // not visited yet
+              switch (ReferenceElements<double,dim>::general(gt).type(i,1))
+              {
+              case simplex :
+                index = numTriFaces_++;
+                break;
+              case cube :
+                index = numQuadFaces_++;
+                break;
+              default :
+                DUNE_THROW(GridError, "wrong geometry type in face");
+              }
+          }
 
       }
 
@@ -453,6 +580,15 @@ namespace Dune {
       myTypes_[dim-1].resize(0);
       myTypes_[dim-1].push_back(cube);
 
+      if (dim==3)
+      {
+        myTypes_[1].resize(0);
+        if (numTriFaces_ > 0)
+          myTypes_[1].push_back(simplex);
+        if (numQuadFaces_ > 0)
+          myTypes_[1].push_back(cube);
+      }
+
       // //////////////////////////////
       //   Init the vertex indices
       // //////////////////////////////
@@ -478,10 +614,8 @@ namespace Dune {
     int numCubes_;
     int numVertices_;
     int numEdges_;
-#if (dim==3)
     int numTriFaces_;
     int numQuadFaces_;
-#endif
 
     std::vector<GeometryType> myTypes_[dim+1];
   };

grid/uggrid/ugintersectionit.cc

@@ -65,7 +65,7 @@ inline const typename UGGridIntersectionIterator<GridImp>::LocalGeometry&
 UGGridIntersectionIterator<GridImp>::
 intersectionSelfLocal() const
 {
-  DUNE_THROW(NotImplemented, "intersection_self_local()");
+  //    DUNE_THROW(NotImplemented, "intersection_self_local()");
   return fakeNeigh_;
 }
 

grid/uggrid/ugtypes.hh

@@ -15,6 +15,7 @@ namespace UG2d {
   union element;
   struct node;
   struct edge;
+  struct vector;
 };
 
 namespace UG3d {
@@ -25,6 +26,7 @@ namespace UG3d {
   union element;
   struct node;
   struct edge;
+  struct vector;
 };
 
 
@@ -68,6 +70,27 @@ namespace Dune {
   /*****************************************************************/
   /*****************************************************************/
 
+  template <int dim>
+  class UGVectorType
+  {
+  public:
+    typedef void T;
+
+  };
+
+  template <>
+  class UGVectorType<3>
+  {
+  public:
+    typedef UG3d::vector T;
+  };
+
+  template <>
+  class UGVectorType<2>
+  {
+  public:
+    typedef UG2d::vector T;
+  };
 
   template <int codim, int dim>
   class TargetType