diff --git a/grid/uggrid.hh b/grid/uggrid.hh index 352de0e5f04f1fa41d4c86c027aeb6ab41964a3d..9e95ce97a08f86b6ec972aa03e5b9f88e244d6cd 100644 --- a/grid/uggrid.hh +++ b/grid/uggrid.hh @@ -178,8 +178,8 @@ namespace Dune { //! Iterator to first entity of given codim on level template<int codim> - UGGridLevelIterator<codim,All_Partition, const UGGrid<dim,dimworld> > lbegin (int level) const; - //typename Traits::template codim<cd>::template partition<pitype>::LevelIterator lbegin (int level) const; + //UGGridLevelIterator<codim,All_Partition, const UGGrid<dim,dimworld> > lbegin (int level) const; + typename Traits::template codim<codim>::LevelIterator lbegin (int level) const; //! one past the end on this level template<int codim> diff --git a/grid/uggrid/uggrid.cc b/grid/uggrid/uggrid.cc index 2197216ba53603e54ed032dd4d5d3a8348234d72..4bd6afcf8af6e0fb02321d4a56735d824075f94a 100644 --- a/grid/uggrid/uggrid.cc +++ b/grid/uggrid/uggrid.cc @@ -89,28 +89,28 @@ public: #endif #ifdef _3 -template <> -class UGGridLevelIteratorFactory<3,3,3,All_Partition> +template <class GridImp> +class UGGridLevelIteratorFactory<3,All_Partition,GridImp> { public: static inline - UGGridLevelIterator<3,3,3,All_Partition> getIterator(UGTypes<3>::GridType* theGrid, int level) { + UGGridLevelIterator<3,All_Partition,GridImp> getIterator(typename UGTypes<GridImp::dimension>::GridType* theGrid, int level) { - UGGridLevelIterator<3,3,3,All_Partition> it(level); + UGGridLevelIterator<3,All_Partition,GridImp> it(level); it.setToTarget(UG_NS<3>::PFirstNode(theGrid), level); return it; } }; -template <> -class UGGridLevelIteratorFactory<0,3,3,All_Partition> +template <class GridImp> +class UGGridLevelIteratorFactory<0,All_Partition,GridImp> { public: static inline - UGGridLevelIterator<0,3,3,All_Partition> getIterator(UGTypes<3>::GridType* theGrid, int level) { + UGGridLevelIterator<0,All_Partition,GridImp> getIterator(UGTypes<3>::GridType* theGrid, int level) { - UGGridLevelIterator<0,3,3,All_Partition> it(level); + UGGridLevelIterator<0,All_Partition,GridImp> it(level); it.setToTarget(UG_NS<3>::PFirstElement(theGrid), level); return it; @@ -293,7 +293,8 @@ inline int UGGrid < dim, dimworld >::maxlevel() const template<int dim, int dimworld> template<int codim> -inline UGGridLevelIterator<codim, All_Partition, const UGGrid<dim,dimworld> > +//inline UGGridLevelIterator<codim, All_Partition, const UGGrid<dim,dimworld> > +typename UGGrid<dim,dimworld>::Traits::template codim<codim>::LevelIterator UGGrid<dim, dimworld>::lbegin (int level) const { assert(multigrid_); @@ -489,11 +490,17 @@ template < int dim, int dimworld > void UGGrid < dim, dimworld >::globalRefine(int refCount) { // mark all entities for grid refinement - UGGridLevelIterator<0, All_Partition, UGGrid<dim,dimworld> > iIt = lbegin<0>(maxlevel()); - UGGridLevelIterator<0, All_Partition, UGGrid<dim,dimworld> > iEndIt = lend<0>(maxlevel()); + //UGGridLevelIterator<0, All_Partition, const UGGrid<dim,dimworld> > iIt = lbegin<0>(maxlevel()); + //UGGridLevelIterator<0, All_Partition, const UGGrid<dim,dimworld> > iEndIt = lend<0>(maxlevel()); + typename UGGrid<dim,dimworld>::Traits::template codim<0>::LevelIterator iIt = lbegin<0>(maxlevel()); + typename UGGrid<dim,dimworld>::Traits::template codim<0>::LevelIterator iEndIt = lend<0>(maxlevel()); + + +#if 0 for (; iIt!=iEndIt; ++iIt) iIt->mark(1); +#endif this->preAdapt(); adapt(); diff --git a/grid/uggrid/uggridgeometry.cc b/grid/uggrid/uggridgeometry.cc index 0698194596542170036fde8435983b4d32946069..f09a7593c943bd3af59dfc19c2ac1570d8309088 100644 --- a/grid/uggrid/uggridgeometry.cc +++ b/grid/uggrid/uggridgeometry.cc @@ -8,6 +8,49 @@ #include <algorithm> +template <int mydim, int coorddim> +struct UGGridGeometryPositionAccess +{}; + +template <> +struct UGGridGeometryPositionAccess<0,3> +{ + static inline + void get(TargetType<3,3>::T* target, + int i, + FieldVector<UGCtype, 3>& coord) { + const UG3d::VERTEX* vertex = target->myvertex; + + coord[0] = vertex->iv.x[0]; + coord[1] = vertex->iv.x[1]; + coord[2] = vertex->iv.x[2]; + } + +}; + + +template <> +struct UGGridGeometryPositionAccess<3,3> +{ + static inline + void get(TargetType<0,3>::T* target, + int i, + FieldVector<UGCtype, 3>& coord) { + if (UG_NS<3>::Tag(target) == UG3d::HEXAHEDRON) { + // Dune numbers the vertices of a hexahedron differently than UG. + // The following two lines do the transformation + const int renumbering[8] = {0, 1, 3, 2, 4, 5, 7, 6}; + i = renumbering[i]; + } + + UG3d::VERTEX* vertex = UG_NS<3>::Corner(target,i)->myvertex; + + for (int j=0; j<3; j++) + coord[j] = vertex->iv.x[j]; + } + +}; + template< int mydim, int coorddim, class GridImp> inline GeometryType UGGridGeometry<mydim,coorddim,GridImp>::type() const @@ -65,11 +108,11 @@ template<int mydim, int coorddim, class GridImp> inline const FieldVector<UGCtype, coorddim>& UGGridGeometry<mydim,coorddim,GridImp>:: operator [](int i) const { - std::cerr << "UGGridGeometry<" << mydim << "," << coorddim << ">::operator[]:\n" - "Default implementation, should not be called!\n"; + UGGridGeometryPositionAccess<mydim,coorddim>::get(target_, i, coord_[i]); return coord_[i]; } +#if 0 #ifdef _3 template<class GridImp> inline const FieldVector<UGCtype, 3>& UGGridGeometry<0,3, GridImp>:: @@ -87,23 +130,7 @@ operator [](int i) const template<class GridImp> inline const FieldVector<UGCtype, 3>& UGGridGeometry<3,3,GridImp>:: operator [](int i) const -{ - assert(0<=i && i<corners()); - - if (type()==hexahedron) { - // Dune numbers the vertices of a hexahedron differently than UG. - // The following two lines do the transformation - const int renumbering[8] = {0, 1, 3, 2, 4, 5, 7, 6}; - i = renumbering[i]; - } - - UG3d::VERTEX* vertex = UG_NS<3>::Corner(target_,i)->myvertex; - - for (int j=0; j<3; j++) - coord_[i][j] = vertex->iv.x[j]; - - return coord_[i]; -} +{} #endif #ifdef _2 @@ -137,6 +164,8 @@ operator [](int i) const #endif +#endif + template< int mydim, int coorddim, class GridImp> inline FieldVector<UGCtype, coorddim> UGGridGeometry<mydim,coorddim,GridImp>:: global(const FieldVector<UGCtype, mydim>& local) const diff --git a/grid/uggrid/uggridgeometry.hh b/grid/uggrid/uggridgeometry.hh index 96646e89278c93b5dc0205cb6d33dbfdcd3c9ab7..a8ddc8b11c9ffadf5836ae487ee642da8c13b003 100644 --- a/grid/uggrid/uggridgeometry.hh +++ b/grid/uggrid/uggridgeometry.hh @@ -19,7 +19,6 @@ namespace Dune { UGMakeableGeometry() : Geometry<mydim, coorddim, GridImp, UGGridGeometry>(UGGridGeometry<mydim, coorddim, GridImp>()) {}; - //void make (Mat<cdim,mydim+1,sgrid_ctype>& __As) { this->realGeometry.make(__As); } void setToTarget(typename TargetType<coorddim-mydim,coorddim>::T* target) { this->realGeometry.setToTarget(target); @@ -27,6 +26,26 @@ namespace Dune { }; + template<class GridImp> + class UGMakeableGeometry<2,3,GridImp> : public Geometry<2, 3, GridImp, UGGridGeometry> + { + public: + UGMakeableGeometry() : + Geometry<2, 3, GridImp, UGGridGeometry>(UGGridGeometry<2,3,GridImp>()) + {}; + //void make (Mat<cdim,mydim+1,sgrid_ctype>& __As) { this->realGeometry.make(__As); } + + void setCoords(int n, const FieldVector<UGCtype, 3>& pos) { + this->realGeometry.coord_[n] = pos; + } + + void setNumberOfCorners(int n) { + this->realGeometry.setNumberOfCorners(n); + } + }; + + + //********************************************************************** // // --UGGridElement @@ -161,6 +180,8 @@ namespace Dune { template <class GridImp_> friend class UGGridIntersectionIterator; + friend class UGMakeableGeometry<2,3,GridImp>; + public: /** \brief Default constructor */ @@ -174,7 +195,7 @@ namespace Dune { int corners () const {return (elementType_==triangle) ? 3 : 4;} //! access to coordinates of corners. Index is the number of the corner - const FieldVector<UGCtype, 3>& operator[] (int i) { + const FieldVector<UGCtype, 3>& operator[] (int i) const { return coord_[i]; } diff --git a/grid/uggrid/ugintersectionit.cc b/grid/uggrid/ugintersectionit.cc index efb43815c167b622a4c5a214b510e27075b88775..47110ddabb5ae29ac0a4d9c08b4fa42f4f992ee9 100644 --- a/grid/uggrid/ugintersectionit.cc +++ b/grid/uggrid/ugintersectionit.cc @@ -93,7 +93,7 @@ UGGridIntersectionIterator<GridImp>::boundary() const template<class GridImp> inline FieldVector<UGCtype, GridImp::dimensionworld>& -UGGridIntersectionIterator <GridImp>::unit_outer_normal () +UGGridIntersectionIterator <GridImp>::unitOuterNormal () const { // ////////////////////////////////////////////////////// // Implementation for 3D @@ -121,10 +121,6 @@ UGGridIntersectionIterator <GridImp>::unit_outer_normal () FieldVector<UGCtype, 3> ba = bPos - aPos; FieldVector<UGCtype, 3> ca = cPos - aPos; - // std::cout << "aPos: " << aPos << std::endl; - // std::cout << "bPos: " << bPos << std::endl; - // std::cout << "cPos: " << cPos << std::endl; - #define V3_VECTOR_PRODUCT(A,B,C) {(C)[0] = (A)[1]*(B)[2] - (A)[2]*(B)[1];\ (C)[1] = (A)[2]*(B)[0] - (A)[0]*(B)[2];\ (C)[2] = (A)[0]*(B)[1] - (A)[1]*(B)[0];} @@ -160,24 +156,24 @@ UGGridIntersectionIterator <GridImp>::unit_outer_normal () template<class GridImp> inline FieldVector<UGCtype, GridImp::dimensionworld>& UGGridIntersectionIterator < GridImp >:: -unit_outer_normal (const FieldVector<UGCtype, GridImp::dimension-1>& local) +unitOuterNormal (const FieldVector<UGCtype, GridImp::dimension-1>& local) const { - return unit_outer_normal(); + return unitOuterNormal(); } template< class GridImp> -inline UGGridGeometry< GridImp::dimension-1, GridImp::dimension-1,GridImp >& +inline typename UGGridIntersectionIterator<GridImp>::LocalGeometry& UGGridIntersectionIterator<GridImp>:: -intersectionSelfLocal() +intersectionSelfLocal() const { DUNE_THROW(NotImplemented, "intersection_self_local()"); return fakeNeigh_; } template< class GridImp> -inline UGGridGeometry< GridImp::dimension-1, GridImp::dimensionworld,GridImp>& +inline typename UGGridIntersectionIterator<GridImp>::Geometry& UGGridIntersectionIterator<GridImp>:: -intersectionGlobal() +intersectionGlobal() const { int numCornersOfSide = UG_NS<GridImp::dimensionworld>::Corners_Of_Side(center_, neighborCount_); @@ -188,8 +184,13 @@ intersectionGlobal() int cornerIdx = UG_NS<GridImp::dimensionworld>::Corner_Of_Side(center_, neighborCount_, i); typename TargetType<dim,dim>::T* node = UG_NS<GridImp::dimensionworld>::Corner(center_, cornerIdx); + // for (int j=0; j<GridImp::dimensionworld; j++) + // neighGlob_.coord_[i][j] = node->myvertex->iv.x[j]; + /** \todo Avoid the temporary */ + FieldVector<UGCtype, dimworld> tmp; for (int j=0; j<GridImp::dimensionworld; j++) - neighGlob_.coord_[i][j] = node->myvertex->iv.x[j]; + tmp[j] = node->myvertex->iv.x[j]; + neighGlob_.setCoords(i, tmp); } @@ -197,9 +198,9 @@ intersectionGlobal() } template< class GridImp> -inline UGGridGeometry<GridImp::dimension-1, GridImp::dimension,GridImp>& +inline typename UGGridIntersectionIterator<GridImp>::LocalGeometry& UGGridIntersectionIterator<GridImp>:: -intersectionNeighborLocal() +intersectionNeighborLocal() const { DUNE_THROW(NotImplemented, "intersection_neighbor_local()"); return fakeNeigh_; @@ -217,7 +218,7 @@ numberInSelf () const template< class GridImp> inline int UGGridIntersectionIterator<GridImp>:: -numberInNeighbor () +numberInNeighbor () const { const typename TargetType<0,GridImp::dimensionworld>::T* other = target(); diff --git a/grid/uggrid/ugintersectionit.hh b/grid/uggrid/ugintersectionit.hh index 4ffbecef324ec0a85c09ca5cfc92d1bd53f45db0..7f6761db0ba75c05af62e12b2ccc42937a2a5d9d 100644 --- a/grid/uggrid/ugintersectionit.hh +++ b/grid/uggrid/ugintersectionit.hh @@ -33,6 +33,10 @@ namespace Dune { friend class UGGridEntity<0,dim,GridImp>; public: + + typedef typename GridImp::template codim<1>::Geometry Geometry; + typedef typename GridImp::template codim<1>::LocalGeometry LocalGeometry; + //! prefix increment UGGridIntersectionIterator& operator ++(); @@ -66,37 +70,37 @@ namespace Dune { //! return unit outer normal, this should be dependent on local //! coordinates for higher order boundary - FieldVector<UGCtype, GridImp::dimensionworld>& unit_outer_normal (const FieldVector<UGCtype, GridImp::dimension-1>& local); + FieldVector<UGCtype, GridImp::dimensionworld>& unitOuterNormal (const FieldVector<UGCtype, GridImp::dimension-1>& local) const; //! return unit outer normal, if you know it is constant use this function instead - FieldVector<UGCtype, GridImp::dimensionworld>& unit_outer_normal (); + FieldVector<UGCtype, GridImp::dimensionworld>& unitOuterNormal () const; //! intersection of codimension 1 of this neighbor with element where //! iteration started. //! Here returned element is in LOCAL coordinates of the element //! where iteration started. - UGGridGeometry<GridImp::dimension-1,GridImp::dimension-1,GridImp>& intersectionSelfLocal (); + LocalGeometry& intersectionSelfLocal () const; //! intersection of codimension 1 of this neighbor with element where iteration started. //! Here returned element is in GLOBAL coordinates of the element where iteration started. - UGGridGeometry<GridImp::dimension-1,GridImp::dimensionworld,GridImp>& intersectionGlobal (); + Geometry& intersectionGlobal () const; //! local number of codim 1 entity in self where intersection is contained in int numberInSelf () const; //! intersection of codimension 1 of this neighbor with element where iteration started. //! Here returned element is in LOCAL coordinates of neighbor - UGGridGeometry<GridImp::dimension-1,GridImp::dimension,GridImp>& intersectionNeighborLocal (); + LocalGeometry& intersectionNeighborLocal () const; //! local number of codim 1 entity in neighbor where intersection is contained - int numberInNeighbor (); + int numberInNeighbor () const; //! return outer normal, this should be dependent on local //! coordinates for higher order boundary - FieldVector<UGCtype, dimworld>& outer_normal (const FieldVector<UGCtype, dim-1>& local); + FieldVector<UGCtype, dimworld>& outerNormal (const FieldVector<UGCtype, dim-1>& local) const; //! return unit outer normal, if you know it is constant use this function instead - FieldVector<UGCtype, dimworld>& outer_normal (); + FieldVector<UGCtype, dimworld>& outerNormal () const; private: //********************************************************** @@ -116,7 +120,7 @@ namespace Dune { UGGridEntity<0,dim,GridImp> virtualEntity_; //! vector storing the outer normal - FieldVector<UGCtype, dimworld> outerNormal_; + mutable FieldVector<UGCtype, dimworld> outerNormal_; //! pointer to element holding the self_local and self_global information. //! This element is created on demand. @@ -124,7 +128,7 @@ namespace Dune { //! pointer to element holding the neighbor_global and neighbor_local //! information. This element is created on demand. - UGGridGeometry<dim-1,dimworld,GridImp> neighGlob_; + mutable UGMakeableGeometry<dim-1,dimworld,GridImp> neighGlob_; //! BoundaryEntity UGGridBoundaryEntity<GridImp> boundaryEntity_;