1#ifndef DUNE_ALU3DGRID_HSFC_HH
2#define DUNE_ALU3DGRID_HSFC_HH
9#include <dune/common/exceptions.hh>
11#include <dune/common/parallel/mpihelper.hh>
12#include <dune/common/parallel/communication.hh>
13#include <dune/common/parallel/mpicommunication.hh>
15#include <dune/alugrid/impl/parallel/zcurve.hh>
17#include <dune/alugrid/impl/parallel/aluzoltan.hh>
20 extern double Zoltan_HSFC_InvHilbert3d (Zoltan_Struct *zz,
double *coord);
21 extern double Zoltan_HSFC_InvHilbert2d (Zoltan_Struct *zz,
double *coord);
27 template <
class Coordinate>
31 typedef Dune :: CollectiveCommunication< typename Dune :: MPIHelper :: MPICommunicator >
32 CollectiveCommunication ;
35 typedef void Zoltan_Struct;
36 typedef CollectiveCommunication Zoltan;
40 typedef double zoltan_hsfc_inv_t(Zoltan_Struct *zz,
double *coord);
42 static const int dimension = Coordinate::dimension;
47 zoltan_hsfc_inv_t* hsfcInv_;
53 const Coordinate& upper,
54 const CollectiveCommunication& comm =
55 CollectiveCommunication( Dune::MPIHelper::getCommunicator() ) )
59 hsfcInv_( dimension == 3 ? Zoltan_HSFC_InvHilbert3d : Zoltan_HSFC_InvHilbert2d ),
73 assert( point.size() == (
unsigned int)dimension );
75 Coordinate center( 0 ) ;
77 for(
int d=0; d<dimension; ++d )
78 center[ d ] = (point[ d ] - lower_[ d ]) / length_[ d ];
81 return hsfcInv_( zz_.Get_C_Handle(), ¢er[ 0 ] );
83 DUNE_THROW(Dune::SystemError,
"Zoltan not found, cannot use Zoltan's Hilbert curve");
89 double index(
const Coordinate& point )
const
95 template<
class Gr
idView >
98 typedef typename GridView :: template Codim< 0 > :: Iterator Iterator ;
99 typedef typename Iterator :: Entity :: Geometry :: GlobalCoordinate GlobalCoordinate ;
101 std::vector< GlobalCoordinate > vertices;
102 vertices.reserve( view.indexSet().size( 0 ) );
104 const Iterator endit = view.template end< 0 > ();
105 for(Iterator it = view.template begin< 0 > (); it != endit; ++ it )
107 GlobalCoordinate center = (*it).geometry().center();
108 vertices.push_back( center );
111 std::stringstream gnufilename;
112 gnufilename << name <<
".gnu";
113 if( view.grid().comm().size() > 1 )
114 gnufilename <<
"." << view.grid().comm().rank();
116 std::ofstream gnuFile ( gnufilename.str() );
119 for(
size_t i=0; i<vertices.size(); ++i )
121 gnuFile << vertices[ i ] << std::endl;
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() );
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;
142 for(
size_t i=0; i<vertices.size(); ++i )
144 vtkFile << vertices[ i ];
145 for(
int d=GlobalCoordinate::dimension; d<3; ++d )
147 vtkFile << std::endl;
151 vtkFile <<
"LINES " << vertices.size()-1 <<
" " << (vertices.size()-1)*3 << std::endl;
153 for(
size_t i=0; i<vertices.size()-1; ++i )
154 vtkFile <<
"2 " << i <<
" " << i+1 << std::endl;
165 template <
class Coordinate>
168 typedef ::ALUGridSFC::ZCurve< long int, Coordinate::dimension> ZCurveOrderingType;
172 typedef Dune :: CollectiveCommunication< typename MPIHelper :: MPICommunicator >
173 CollectiveCommunication ;
191 const Coordinate& lower,
192 const Coordinate& upper,
193 const CollectiveCommunication& comm =
194 CollectiveCommunication( Dune::MPIHelper::getCommunicator() ) )
195 :
zCurve_ ( lower, upper, comm )
201 template <
class OtherComm>
203 const Coordinate& lower,
204 const Coordinate& upper,
205 const OtherComm& otherComm )
207 otherComm.size() > 1 ? CollectiveCommunication(
Dune::MPIHelper::getCommunicator() ) :
208 CollectiveCommunication(
Dune::MPIHelper::getLocalCommunicator() ) )
210 otherComm.size() > 1 ? CollectiveCommunication(
Dune::MPIHelper::getCommunicator() ) :
211 CollectiveCommunication(
Dune::MPIHelper::getLocalCommunicator() ) )
217 double index(
const Coordinate& point )
const
221 return double(
zCurve_.index( point ) );
229 DUNE_THROW(NotImplemented,
"Wrong space filling curve ordering selected");
Definition: alu3dinclude.hh:63
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
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