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
18extern "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
24namespace 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
162namespace 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
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
Definition alu3dinclude.hh:80
Definition hsfc.hh:24
void printSpaceFillingCurve(const GridView &view, std::string name="sfc", const bool vtk=false)
Definition hsfc.hh:95
ZoltanSpaceFillingCurveOrdering(const Coordinate &lower, const Coordinate &upper, const CollectiveCommunication &comm=CollectiveCommunication(Dune::MPIHelper::getCommunicator()))
Definition hsfc.hh:51
double index(const Coordinate &point) const
Definition hsfc.hh:88
double hilbertIndex(const Coordinate &point) const
Definition hsfc.hh:69
Definition hsfc.hh:166
const CurveType curveType_
Definition hsfc.hh:186
SpaceFillingCurveOrdering(const CurveType &curveType, const Coordinate &lower, const Coordinate &upper, const OtherComm &otherComm)
Definition hsfc.hh:201
ZCurveOrderingType zCurve_
Definition hsfc.hh:183
double index(const Coordinate &point) const
Definition hsfc.hh:216
static const CurveType DefaultCurve
Definition hsfc.hh:179
SpaceFillingCurveOrdering(const CurveType &curveType, const Coordinate &lower, const Coordinate &upper, const CollectiveCommunication &comm=CollectiveCommunication(Dune::MPIHelper::getCommunicator()))
Definition hsfc.hh:189
HilbertOrderingType hilbert_
Definition hsfc.hh:184
CurveType
Definition hsfc.hh:174
@ Hilbert
Definition hsfc.hh:174
@ ZCurve
Definition hsfc.hh:174
@ None
Definition hsfc.hh:174