1 #ifndef DUNE_ALU3DGRID_HSFC_HH
2 #define DUNE_ALU3DGRID_HSFC_HH
3 
4 #include <string>
5 #include <sstream>
6 #include <fstream>
7 #include <vector>
8 
9 #include <dune/common/exceptions.hh>
10 
11 #include <dune/common/parallel/mpihelper.hh>
12 #include <dune/common/parallel/communication.hh>
13 #include <dune/common/parallel/mpicommunication.hh>
14 
15 #include <dune/alugrid/impl/parallel/zcurve.hh>
16 #if HAVE_ZOLTAN
17 #include <dune/alugrid/impl/parallel/aluzoltan.hh>
18 
19 extern "C" {
20   extern double Zoltan_HSFC_InvHilbert3d (Zoltan_Struct *zz, double *coord);
21   extern double Zoltan_HSFC_InvHilbert2d (Zoltan_Struct *zz, double *coord);
22 }
23 #endif
24 
25 namespace ALUGridSFC {
26 
27   template <class Coordinate>
28   class ZoltanSpaceFillingCurveOrdering
29   {
30     // type of communicator
31     typedef Dune :: CollectiveCommunication< typename Dune :: MPIHelper :: MPICommunicator >
32         CollectiveCommunication ;
33 
34 #if ! HAVE_ZOLTAN
35     typedef void                      Zoltan_Struct;
36     typedef CollectiveCommunication   Zoltan;
37 #endif
38 
39     // type of Zoltan HSFC ordering function
40     typedef double zoltan_hsfc_inv_t(Zoltan_Struct *zz, double *coord);
41 
42     static const int dimension = Coordinate::dimension;
43 
44     Coordinate lower_;
45     Coordinate length_;
46 
47     zoltan_hsfc_inv_t* hsfcInv_;
48 
49     mutable Zoltan zz_;
50 
51   public:
ZoltanSpaceFillingCurveOrdering(const Coordinate & lower,const Coordinate & upper,const CollectiveCommunication & comm=CollectiveCommunication (Dune::MPIHelper::getCommunicator ()))52     ZoltanSpaceFillingCurveOrdering( const Coordinate& lower,
53                                      const Coordinate& upper,
54                                      const CollectiveCommunication& comm =
55                                      CollectiveCommunication( Dune::MPIHelper::getCommunicator() ) )
56       : lower_( lower ),
57         length_( upper ),
58 #if HAVE_ZOLTAN
59         hsfcInv_( dimension == 3 ? Zoltan_HSFC_InvHilbert3d : Zoltan_HSFC_InvHilbert2d ),
60 #else
61         hsfcInv_( 0 ),
62 #endif
63         zz_( comm )
64     {
65       // compute length
66       length_ -= lower_;
67     }
68 
69     // return unique hilbert index in interval [0,1] given an element's center
hilbertIndex(const Coordinate & point) const70     double hilbertIndex( const Coordinate& point ) const
71     {
72 #if HAVE_ZOLTAN
73       assert( point.size() == (unsigned int)dimension );
74 
75       Coordinate center( 0 ) ;
76       // scale center into [0,1]^d box which is needed by Zoltan_HSFC_InvHilbert{2,3}d
77       for( int d=0; d<dimension; ++d )
78         center[ d ] = (point[ d ] - lower_[ d ]) / length_[ d ];
79 
80       // return hsfc index in interval [0,1]
81       return hsfcInv_( zz_.Get_C_Handle(), &center[ 0 ] );
82 #else
83       DUNE_THROW(Dune::SystemError,"Zoltan not found, cannot use Zoltan's Hilbert curve");
84       return 0.0;
85 #endif
86     }
87 
88     // return unique hilbert index in interval [0,1] given an element's center
index(const Coordinate & point) const89     double index( const Coordinate& point ) const
90     {
91       return hilbertIndex( point );
92     }
93   };
94 
95   template< class GridView >
printSpaceFillingCurve(const GridView & view,std::string name="sfc",const bool vtk=false)96   void printSpaceFillingCurve ( const GridView& view, std::string name = "sfc", const bool vtk = false  )
97   {
98     typedef typename GridView :: template Codim< 0 > :: Iterator  Iterator ;
99     typedef typename Iterator :: Entity :: Geometry :: GlobalCoordinate GlobalCoordinate ;
100 
101     std::vector< GlobalCoordinate > vertices;
102     vertices.reserve( view.indexSet().size( 0 ) );
103 
104     const Iterator endit = view.template end< 0 > ();
105     for(Iterator it = view.template begin< 0 > (); it != endit; ++ it )
106     {
107       GlobalCoordinate center = (*it).geometry().center();
108       vertices.push_back( center );
109     }
110 
111     std::stringstream gnufilename;
112     gnufilename << name << ".gnu";
113     if( view.grid().comm().size() > 1 )
114       gnufilename << "." << view.grid().comm().rank();
115 
116     std::ofstream gnuFile ( gnufilename.str() );
117     if( gnuFile )
118     {
119       for( size_t i=0; i<vertices.size(); ++i )
120       {
121         gnuFile << vertices[ i ] << std::endl;
122       }
123       gnuFile.close();
124     }
125 
126     if( vtk )
127     {
128       std::stringstream vtkfilename;
129       vtkfilename << name << ".vtk";
130       if( view.grid().comm().size() > 1 )
131         vtkfilename << "." << view.grid().comm().rank();
132       std::ofstream vtkFile ( vtkfilename.str() );
133 
134       if( vtkFile )
135       {
136         vtkFile << "# vtk DataFile Version 1.0" << std::endl;
137         vtkFile << "Line representation of vtk" << std::endl;
138         vtkFile << "ASCII" << std::endl;
139         vtkFile << "DATASET POLYDATA" << std::endl;
140         vtkFile << "POINTS "<< vertices.size() << " FLOAT" << std::endl;
141 
142         for( size_t i=0; i<vertices.size(); ++i )
143         {
144           vtkFile << vertices[ i ];
145           for( int d=GlobalCoordinate::dimension; d<3; ++d )
146             vtkFile << " 0";
147           vtkFile << std::endl;
148         }
149 
150         // lines, #lines, #entries
151         vtkFile << "LINES " << vertices.size()-1 << " " << (vertices.size()-1)*3 << std::endl;
152 
153         for( size_t i=0; i<vertices.size()-1; ++i )
154           vtkFile << "2 " << i << " " << i+1 << std::endl;
155 
156         vtkFile.close();
157       }
158     }
159   }
160 
161 } // end namespace ALUGridSFC
162 
163 namespace Dune {
164 
165   template <class Coordinate>
166   class SpaceFillingCurveOrdering
167   {
168     typedef ::ALUGridSFC::ZCurve< long int, Coordinate::dimension>      ZCurveOrderingType;
169     typedef ::ALUGridSFC::ZoltanSpaceFillingCurveOrdering< Coordinate > HilbertOrderingType;
170 
171     // type of communicator
172     typedef Dune :: CollectiveCommunication< typename MPIHelper :: MPICommunicator >
173         CollectiveCommunication ;
174   public:
175     enum CurveType { ZCurve, Hilbert, None };
176 
177 #if HAVE_ZOLTAN
178     static const CurveType DefaultCurve = Hilbert ;
179 #else
180     static const CurveType DefaultCurve = ZCurve ;
181 #endif
182 
183   protected:
184     ZCurveOrderingType   zCurve_;
185     HilbertOrderingType  hilbert_;
186 
187     const CurveType curveType_;
188 
189   public:
SpaceFillingCurveOrdering(const CurveType & curveType,const Coordinate & lower,const Coordinate & upper,const CollectiveCommunication & comm=CollectiveCommunication (Dune::MPIHelper::getCommunicator ()))190     SpaceFillingCurveOrdering( const CurveType& curveType,
191                                const Coordinate& lower,
192                                const Coordinate& upper,
193                                const CollectiveCommunication& comm =
194                                CollectiveCommunication( Dune::MPIHelper::getCommunicator() ) )
195       : zCurve_ ( lower, upper, comm )
196       , hilbert_( lower, upper, comm )
197       , curveType_( curveType )
198     {
199     }
200 
201     template <class OtherComm>
SpaceFillingCurveOrdering(const CurveType & curveType,const Coordinate & lower,const Coordinate & upper,const OtherComm & otherComm)202     SpaceFillingCurveOrdering( const CurveType& curveType,
203                                const Coordinate& lower,
204                                const Coordinate& upper,
205                                const OtherComm& otherComm )
206       : zCurve_ ( lower, upper,
207                   otherComm.size() > 1 ? CollectiveCommunication( Dune::MPIHelper::getCommunicator() ) :
208                                          CollectiveCommunication( Dune::MPIHelper::getLocalCommunicator() ) )
209       , hilbert_( lower, upper,
210                   otherComm.size() > 1 ? CollectiveCommunication( Dune::MPIHelper::getCommunicator() ) :
211                                          CollectiveCommunication( Dune::MPIHelper::getLocalCommunicator() ) )
212       , curveType_( curveType )
213     {
214     }
215 
216     // return unique hilbert index in interval [0,1] given an element's center
index(const Coordinate & point) const217     double index( const Coordinate& point ) const
218     {
219       if( curveType_ == ZCurve )
220       {
221         return double( zCurve_.index( point ) );
222       }
223       else if ( curveType_ == Hilbert )
224       {
225         return double( hilbert_.index( point ) );
226       }
227       else
228       {
229         DUNE_THROW(NotImplemented,"Wrong space filling curve ordering selected");
230         return 0.0;
231       }
232     }
233   };
234 
235 } // end namespace Dune
236 
237 #endif // #ifndef DUNE_ALU3DGRID_HSFC_HH
238