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