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(), ¢er[ 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