implemented geometryInFather correctly.

Also removed not needed getCoordVec method.

[[Imported from SVN: r3291]]

grid/alu3dgrid/entity_imp.cc

@@ -134,6 +134,7 @@ namespace Dune {
   inline const typename ALU3dGridEntity<cd,dim,GridImp>::Geometry &
   ALU3dGridEntity<cd,dim,GridImp>:: geometry() const
   {
+    assert( (cd == 1) ? (face_ >= 0) : 1 );
     if(!builtgeometry_) builtgeometry_ = geo_.buildGeom(*item_, twist_, face_ );
     return geo_;
   }
@@ -330,12 +331,10 @@ namespace Dune {
   inline const typename ALU3dGridEntity<0,dim,GridImp>::Geometry &
   ALU3dGridEntity<0,dim,GridImp> :: geometryInFather () const
   {
-    const Geometry & vati   = (*this->father()).geometry();
-    const Geometry & myself = this->geometry();
-
-    for(int i=0; i<vati.corners(); i++)
-      geoInFather_.getCoordVec(i) = vati.local( myself[i] );
+    const typename GridImp::template Codim<0> ::
+    EntityPointer ep = father();
 
+    geoInFather_.buildGeomInFather( (*ep).geometry() , geometry() );
     return geoInFather_;
   }
 
@@ -684,6 +683,7 @@ namespace Dune {
       , done_ (false)
   {
     assert( entity_ );
+    if( codim == 1 ) assert( face >= 0 );
     (*entity_).setElement( const_cast<MyHElementType &> (item), twist , face );
   }
 

grid/alu3dgrid/geometry.hh

@@ -65,12 +65,12 @@ namespace Dune {
     //! build geometry out of different ALU3dGrid Geometrys
     //! ItemType are HElementType, HFaceType, HEdgeType and VertexType
     template <class ItemType>
-    bool buildGeom(const ItemType & item, int twist = 0, int face = -1 )
+    bool buildGeom(const ItemType & item, int twist = 0 , int face = -1 )
     {
       return this->realGeometry.buildGeom(item, twist,face);
     }
 
-    bool buildGeom(const ALU3DSPACE HFaceType& item, int twist, int face = -1 ) {
+    bool buildGeom(const ALU3DSPACE HFaceType& item, int twist, int face ) {
       return this->realGeometry.buildGeom(item, twist, face );
     }
 
@@ -84,18 +84,18 @@ namespace Dune {
       return this->realGeometry.buildGhost(ghost);
     }
 
-    // print real entity for debugging
-    void print (std::ostream& ss) const
+    //! build geometry of local coordinates relative to father
+    template <class GeometryType>
+    bool buildGeomInFather(const GeometryType &fatherGeom , const GeometryType & myGeom)
     {
-      this->realGeometry.print(ss);
+      return this->realGeometry.buildGeomInFather(fatherGeom,myGeom);
     }
 
-    // for changing the coordinates of one element
-    FieldVector<alu3d_ctype, coorddim> & getCoordVec (int i)
+    // print real entity for debugging
+    void print (std::ostream& ss) const
     {
-      return this->realGeometry.getCoordVec(i);
+      this->realGeometry.print(ss);
     }
-
   };
 
   //! ALU3dGridGeometry
@@ -167,13 +167,14 @@ namespace Dune {
     //! build ghost out of internal boundary segment
     bool buildGhost(const PLLBndFaceType & ghost);
 
+    //! build geometry of local coordinates relative to father
+    template <class GeometryType>
+    bool buildGeomInFather(const GeometryType &fatherGeom , const GeometryType & myGeom);
+
     //! print internal data
     //! no interface method
     void print (std::ostream& ss) const;
 
-    // for changing the coordinates of one element
-    FieldVector<alu3d_ctype, cdim> & getCoordVec (int i);
-
   private:
     //! calculates the vertex index in the reference element out of a face index
     //! and a local vertex index
@@ -273,14 +274,22 @@ namespace Dune {
     //! build ghost out of internal boundary segment
     bool buildGhost(const PLLBndFaceType & ghost);
 
+    //! build geometry of local coordinates relative to father
+    template <class GeometryType>
+    bool buildGeomInFather(const GeometryType &fatherGeom , const GeometryType & myGeom);
+
     //! print internal data
     //! no interface method
     void print (std::ostream& ss) const;
 
-    // for changing the coordinates of one element
-    FieldVector<alu3d_ctype, cdim> & getCoordVec (int i);
 
   private:
+    // create triMap from coordinates
+    void buildMapping();
+
+    // delete triMap
+    void removeMapping();
+
     //! the vertex coordinates
     mutable FieldMatrix<alu3d_ctype, Power_m_p<2,mydim>::power, cdim> coord_;
     //mutable FieldVector<alu3d_ctype, mydim> tmp1_;

grid/alu3dgrid/geometry_imp.cc

@@ -95,6 +95,23 @@ namespace Dune {
     }
   }
 
+  // built Geometry
+  template <int mydim, int cdim>
+  template <class GeometryType>
+  inline bool
+  ALU3dGridGeometry<mydim, cdim, const ALU3dGrid<3, 3, tetra> >::
+  buildGeomInFather(const GeometryType &fatherGeom , const GeometryType & myGeom)
+  {
+    // reset flags, because mappings need to be calculated again
+    builtinverse_ = builtA_ = builtDetDF_ = false;
+
+    // compute the local coordinates in father refelem
+    for(int i=0; i < myGeom.corners() ; i++)
+      coord_[i] = fatherGeom.local( myGeom[i] );
+
+    return true;
+  }
+
   template <>
   inline bool ALU3dGridGeometry<3,3, const ALU3dGrid<3,3,tetra> > ::
   buildGeom(const IMPLElementType & item, int, int)
@@ -262,14 +279,6 @@ namespace Dune {
     return coord_[i];
   }
 
-  template<int mydim, int cdim>
-  inline FieldVector<alu3d_ctype, cdim>&
-  ALU3dGridGeometry<mydim,cdim,const ALU3dGrid<3, 3, tetra> > :: getCoordVec (int i)
-  {
-    assert((i>=0) && (i < mydim+1));
-    return coord_[i];
-  }
-
   //   G L O B A L   - - -
 
   // dim = 1,2,3 dimworld = 3
@@ -414,14 +423,6 @@ namespace Dune {
     return Power_m_p<2,mydim>::power;
   }
 
-  template<int mydim, int cdim>
-  inline FieldVector<alu3d_ctype, cdim>&
-  ALU3dGridGeometry<mydim,cdim,const ALU3dGrid<3, 3, hexa> > :: getCoordVec (int i)
-  {
-    assert((i>=0) && (i < corners()));
-    return coord_[i];
-  }
-
   template <int mydim, int cdim>
   const FieldVector<alu3d_ctype, cdim>&
   ALU3dGridGeometry<mydim, cdim, const ALU3dGrid<3, 3, hexa> >::
@@ -516,6 +517,48 @@ namespace Dune {
     ss << "} \n";
   }
 
+
+  template <int mydim, int cdim>
+  inline void
+  ALU3dGridGeometry<mydim, cdim, const ALU3dGrid<3, 3, hexa> >::
+  buildMapping()
+  {
+    assert( mydim == 3 );
+    assert( cdim  == 3 );
+    assert( triMap_ == 0 );
+    triMap_ = new TrilinearMapping(coord_[0], coord_[1], coord_[2], coord_[3],
+                                   coord_[4], coord_[5], coord_[6], coord_[7]);
+  }
+
+  template <int mydim, int cdim>
+  inline void
+  ALU3dGridGeometry<mydim, cdim, const ALU3dGrid<3, 3, hexa> >::
+  removeMapping()
+  {
+    if(triMap_)
+    {
+      delete triMap_;
+      triMap_ = 0;
+    }
+  }
+
+  // built Geometry
+  template <int mydim, int cdim>
+  template <class GeometryType>
+  inline bool
+  ALU3dGridGeometry<mydim, cdim, const ALU3dGrid<3, 3, hexa> >::
+  buildGeomInFather(const GeometryType &fatherGeom , const GeometryType & myGeom)
+  {
+    removeMapping();
+
+    // compute the local coordinates in father refelem
+    for(int i=0; i < myGeom.corners() ; i++)
+      coord_[i] = fatherGeom.local( myGeom[i] );
+
+    buildMapping();
+    return true;
+  }
+
   template <>
   inline bool
   ALU3dGridGeometry<3, 3, const ALU3dGrid<3, 3, hexa> >::
@@ -530,11 +573,9 @@ namespace Dune {
         coord_[i][j] = p[j];
       }
     }
-
-    if(triMap_) delete triMap_;
-    triMap_ = new TrilinearMapping(coord_[0], coord_[1], coord_[2], coord_[3],
-                                   coord_[4], coord_[5], coord_[6], coord_[7]);
-
+    // delete old mapping
+    removeMapping();
+    buildMapping();
     return true;
   }
 
@@ -605,6 +646,9 @@ namespace Dune {
 
     const GEOFaceType& face = static_cast<const GEOFaceType&> (item);
 
+    if( faceNum >= 6 ) std::cout << "\n" << faceNum << " wrong face num (too big) \n";
+    if( faceNum < 0 ) std::cout << "\n" << faceNum << " wrong face num \n";
+    // for all vertices of this face
     for (int i = 0; i < 4; ++i)
     {
       // Transform Dune index to ALU index and apply twist