dune-alugrid  2.6-git
hsfc.hh
Go to the documentation of this file.
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/collectivecommunication.hh>
13 #include <dune/common/parallel/mpicollectivecommunication.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>
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:
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
70  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
89  double index( const Coordinate& point ) const
90  {
91  return hilbertIndex( point );
92  }
93  };
94 
95  template< class GridView >
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>
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:
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_;
186 
188 
189  public:
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>
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
217  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
Definition: alu3dinclude.hh:80
Definition: hsfc.hh:25
void printSpaceFillingCurve(const GridView &view, std::string name="sfc", const bool vtk=false)
Definition: hsfc.hh:96
ZoltanSpaceFillingCurveOrdering(const Coordinate &lower, const Coordinate &upper, const CollectiveCommunication &comm=CollectiveCommunication(Dune::MPIHelper::getCommunicator()))
Definition: hsfc.hh:52
double index(const Coordinate &point) const
Definition: hsfc.hh:89
double hilbertIndex(const Coordinate &point) const
Definition: hsfc.hh:70
Definition: hsfc.hh:167
const CurveType curveType_
Definition: hsfc.hh:187
SpaceFillingCurveOrdering(const CurveType &curveType, const Coordinate &lower, const Coordinate &upper, const OtherComm &otherComm)
Definition: hsfc.hh:202
ZCurveOrderingType zCurve_
Definition: hsfc.hh:184
double index(const Coordinate &point) const
Definition: hsfc.hh:217
static const CurveType DefaultCurve
Definition: hsfc.hh:180
SpaceFillingCurveOrdering(const CurveType &curveType, const Coordinate &lower, const Coordinate &upper, const CollectiveCommunication &comm=CollectiveCommunication(Dune::MPIHelper::getCommunicator()))
Definition: hsfc.hh:190
HilbertOrderingType hilbert_
Definition: hsfc.hh:185
CurveType
Definition: hsfc.hh:175
@ Hilbert
Definition: hsfc.hh:175
@ ZCurve
Definition: hsfc.hh:175
@ None
Definition: hsfc.hh:175