Commit 7d6f6083 authored by BLANCHARD Pierre's avatar BLANCHARD Pierre
Browse files

Provided new distribution of particles (hyperbolic paraboloid Z=C*(AX*AX-BY*BY)).

parent a7feed21
......@@ -189,12 +189,24 @@ int main(int argc, char ** argv){
pos = aspectRatio.find(":"); aspectRatio.replace(pos,1," ");
std::stringstream ss(aspectRatio); ss >>A >> B >> C ;
if(A != B){
std::cerr << " A /= B in prolate sllipsoide A =B. Your aspect ratio: "<< aspectRatio<<std::endl;
std::cerr << " A /= B in prolate ellipsoide A =B. Your aspect ratio: "<< aspectRatio<<std::endl;
}
std::cout << "A: "<<A<<" B "<< B << " C: " << C<<std::endl;
std::cout << "A: "<<A<<" B: "<< B << " C: " << C<<std::endl;
unifRandonPointsOnProlate(NbPoints,A,C,particles);
BoxWith = 2.0*C;
} //const FSize NbPoints = FParameters::getValue(argc,argv,FParameterDefinitions::NbParticles.options, FSize(20000));
} //const FSize NbPoints = FParameters::getValue(argc,argv,FParameterDefinitions::NbParticles.options, FSize(20000));
else if(FParameters::existParameter(argc, argv, "-hyperpara")){
std::string dd(":"),aspectRatio = FParameters::getStr(argc,argv,"-ar", "1:1:2");
FReal A,B,C ;
size_t pos = aspectRatio.find(":"); aspectRatio.replace(pos,1," ");
pos = aspectRatio.find(":"); aspectRatio.replace(pos,1," ");
std::stringstream ss(aspectRatio); ss >>A >> B >> C ;
std::cout << "A: "<<A<<" B: "<< B << " C: " << C<<std::endl;
unifRandonPointsOnHyperPara(NbPoints,A,B,C,particles);
BoxWith = 2.0*FMath::Max( A,FMath::Max( B,C)) ;
std::cout << "BoxWith: " << BoxWith<<std::endl;
}
else if(FParameters::existParameter(argc, argv, "-ellipsoid")){
// else if(FParameters::existParameter(argc, argv, "-ellipsoid")){
std::string dd(":"),aspectRatio = FParameters::getStr(argc,argv,"-ar", "1:1:2");
......@@ -203,7 +215,7 @@ int main(int argc, char ** argv){
size_t pos = aspectRatio.find(":"); aspectRatio.replace(pos,1," ");
pos = aspectRatio.find(":"); aspectRatio.replace(pos,1," ");
std::stringstream ss(aspectRatio); ss >>A >> B >> C ;
std::cout << "A: "<<A<<" B "<< B << " C: " << C<<std::endl;
std::cout << "A: "<<A<<" B: "<< B << " C: " << C<<std::endl;
nonunifRandonPointsOnElipsoid(NbPoints,A,B,C,particles);
BoxWith = 2.0*FMath::Max( A,FMath::Max( B,C)) ;
}
......
......@@ -171,6 +171,31 @@ void unifRandonPointsOnProlate(const FSize N , const FReal &a, const FReal &c, F
} ;
//! \fn unifRandonPointsOnHyperPara(const int N , const FReal &a, const FReal &b, const FReal &c, FReal * points)
//! \brief Generate N points uniformly distributed on the hyperbolic paraboloid of aspect ratio a:b:c
//!
//! \param N the number of points
//! \param a the x semi-axe length
//! \param b the y semi-axe length
//! \param c the z semi-axe length
//! \param points array of size 4*N and stores data as follow x,y,z,0,x,y,z,0....
//!
template <class FReal>
void unifRandonPointsOnHyperPara(const FSize N , const FReal &a, const FReal &b, const FReal &c, FReal * points) {
//
FReal u, v ;
FSize j =0 ;
for (FSize i = 0 ; i< N ; ++i, j+=4) {
u = 2.0*getRandom<FReal>() - 1.0 ; v = 2.0*getRandom<FReal>() - 1.0 ;
points[j] = a*u ;
points[j+1] = b*v ;
points[j+2] = c*(u*u - v*v) ;
}
};
//! \fn unifRandonPointsOnSphere(const int N , const FReal R, FReal * points)
//! \brief Generate N points uniformly distributed on the sphere of radius R
......@@ -247,245 +272,5 @@ void unifRandonPlummer(const FSize N , const FReal R, const FReal M, FReal * poi
<<100*static_cast<FReal>(cpt-N)/static_cast<FReal>(cpt) << " %" <<std::endl;
} ;
#ifdef TOREMOVEOC
//! \fn void exportCVS(std::ofstream& file, const FReal * particles , const int N, const int nbDataPerParticle=4)
//! \brief Export particles in CVS Format
//!
//! Export particles in CVS Format as follow
//! x , y , z , physicalValue, P, FX, FY, FY
//! It is useful to plot the distribution with paraView
//!
//! @param file stream to save the data
//! @param N number of particles
//! @param particles array of particles of type FReal (float or double) Its size is N*nbDataPerParticle
//! @param nbDataPerParticle number of values per particles (default value 4)
//!
template <class FReal>
void exportCVS(std::ofstream& file, const FReal * particles , const FSize N, const int nbDataPerParticle=4){
int j = 0;
if (nbDataPerParticle==4){
file << " x , y , z, q " <<std::endl;
}
else {
file << " x , y , z, q P FX FY FZ" <<std::endl;
}
for(int i = 0 ; i< N; ++i, j+=nbDataPerParticle){
file << particles[j] ;
for (int k = 1 ; k< nbDataPerParticle ; ++k) {
file << " , " << particles[j+k] ;
}
file << std::endl;
}
}
//
//! \fn void exportCOSMOS(std::ofstream& file, const FReal * particles, const int N )
//! \brief Export particles in CVS Format
//!
//! Export particles in CVS Format as follow
//! x , y , z , 0.0, 0.0, 0.0, physicalValue
//!
//! @param file stream to save the data
//! @param particles array of particles of type FReal (float or double) Its size is 4*N (X,Y,Z,M)
//! @param N number of particles
template <class FReal>
void exportCOSMOS(std::ofstream& file, const FReal * particles , const FSize N){
FSize j = 0;
file << " x , y , z, q " <<std::endl;
for(FSize i = 0 ; i< N; ++i, j+=4){
file << particles[j] << " " << particles[j+1] << " " << particles[j+2] << " 0.0 0.0 0.0 " << particles[j+3] <<" " << i << std::endl;
}
}
//
//
//! \fn void exportVTK(std::ofstream& file, const FReal * particles, const int N )
//! \brief Export particles in CVS Format
//!
//! Export particles in old polydata Format.
//! A particle is composed of 4 fields x , y , z , physicalValue
//! It is useful to plot the distribution with paraView
//!
//! @param file stream to save the data
//! @param particles array of particles of type FReal (float or double) Its size is 4*N (X,Y,Z,M)
//! @param N number of particles
template <class FReal>
void exportVTK(std::ofstream& VTKfile, const FReal * particles, const FSize N, const int nbDataPerParticle=4 ){
int j = 0;
//---------------------------
// print generic information
//---------------------------
VTKfile << "# vtk DataFile Version 3.0" << "\n";
VTKfile << "# Generated bt exportVTK" << "\n";
VTKfile << "ASCII" << "\n";
VTKfile << "DATASET POLYDATA" << "\n";
//
//---------------------------------
// print nodes ordered by their TAG
//---------------------------------
VTKfile << "POINTS " << N << " float" << "\n";
//
for(FSize i = 0 ; i< N; ++i, j+=nbDataPerParticle){
VTKfile << particles[j] << " " << particles[j+1] << " " << particles[j+2] <<std::endl;
}
// ------------------------------------------
VTKfile << "\n";
VTKfile << "VERTICES " << N << " " << 2*N << "\n";
for(FSize i = 0 ; i< N; ++i){
VTKfile << " 1 " << " " <<i<<std::endl;
}
VTKfile << "POINT_DATA " << N << "\n";
VTKfile << "SCALARS PhysicalValue float 1" << "\n"
<< "LOOKUP_TABLE default" << "\n" ;
j = 0 ;
for(FSize i = 0 ; i< N; ++i, j+=nbDataPerParticle){
VTKfile << particles[j+3] << " " <<std::endl;
}
VTKfile << "\n";
};
//
//
//! \fn void exportVTKxml(std::ofstream& file, const FReal * particles, const int N )
//! \brief Export particles in xml polydata VTK Format
//!
//! Export particles in the xml polydata VTK Format.
//! A particle is composed of 4 fields x , y , z , physicalValue
//! It is useful to plot the distribution with paraView
//!
//! @param file stream to save the data
//! @param particles array of particles of type FReal (float or double) Its size is 4*N (X,Y,Z,M)
//! @param N number of particles
template <class FReal>
void exportVTKxml(std::ofstream& VTKfile, const FReal * particles, const FSize N ){
FSize j = 0;
VTKfile << "<?xml version=\"1.0\"?>" <<std::endl
<< "<VTKFile type=\"PolyData\" version=\"0.1\" byte_order=\"LittleEndian\"> "<<std::endl
<< "<PolyData>"<<std::endl
<< "<Piece NumberOfPoints=\" " << N << " \" NumberOfVerts=\" "<<N <<" \" NumberOfLines=\" 0\" NumberOfStrips=\"0\" NumberOfPolys=\"0\">"<<std::endl
<< "<Points>"<<std::endl
<< "<DataArray type=\"Float64\" NumberOfComponents=\"3\" format=\"ascii\"> "<<std::endl ;
j = 0 ;
for(FSize i = 0 ; i< N; ++i, j+=4){
VTKfile << particles[j] << " " << particles[j+1] << " " << particles[j+2] << " " ;
}
VTKfile <<std::endl<< "</DataArray> "<<std::endl
<< "</Points> "<<std::endl
<< "<PointData Scalars=\"PhysicalValue\" > "<<std::endl
<< "<DataArray type=\"Float64\" Name=\"PhysicalValue\" format=\"ascii\">"<<std::endl ;
j = 0 ;
for(FSize i = 0 ; i< N; ++i, j+=4){
VTKfile << particles[j+3] << " " ;
}
VTKfile <<std::endl << "</DataArray>"<<std::endl
<< " </PointData>"<<std::endl
<< " <CellData>"<<" </CellData>"<<std::endl
<< " <Verts>"<<std::endl
<< " <DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">"<<std::endl ;
for(FSize i = 0 ; i< N; ++i){
VTKfile << i << " " ;
}
VTKfile<<std::endl << "</DataArray>" <<std::endl
<< "<DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">"<<std::endl ;
for(FSize i = 1 ; i< N+1; ++i){
VTKfile << i << " " ;
}
VTKfile<<std::endl << "</DataArray>"<<std::endl
<< " </Verts>"<<std::endl
<< "<Lines></Lines>"<<std::endl
<< "<Strips></Strips>"<<std::endl
<< "<Polys></Polys>"<<std::endl
<< "</Piece>"<<std::endl
<< "</PolyData>"<<std::endl
<< "</VTKFile>"<<std::endl;
} ;
//
//
//
//! \fn void exportVTKxml(std::ofstream& file, const FReal * particles, const int N, const int nbDataPerParticle )
//! \brief Export particles in CVS Format
//!
//! Export particles in the new PolyData Format.
//! A particle is composed of 4 fields x , y , z , physicalValue
//! It is useful to plot the distribution with paraView
//!
//! @param file stream to save the data
//! @param particles array of particles of type FReal (float or double) Its size is nbDataPerParticle*N
//! @param N number of particles
//! @param nbDataPerParticle number of values per particles (default value 4)
template <class FReal>
void exportVTKxml(std::ofstream& VTKfile, const FReal * particles, const FSize N, const int nbDataPerParticle ){
FSize j = 0;
VTKfile << "<?xml version=\"1.0\"?>" <<std::endl
<< "<VTKFile type=\"PolyData\" version=\"0.1\" byte_order=\"LittleEndian\"> "<<std::endl
<< "<PolyData>"<<std::endl
<< "<Piece NumberOfPoints=\" " << N << " \" NumberOfVerts=\" "<<N <<" \" NumberOfLines=\" 0\" NumberOfStrips=\"0\" NumberOfPolys=\"0\">"<<std::endl
<< "<Points>"<<std::endl
<< "<DataArray type=\"Float64\" NumberOfComponents=\"3\" format=\"ascii\"> "<<std::endl ;
j = 0 ;
for(FSize i = 0 ; i< N; ++i, j+=nbDataPerParticle){
VTKfile << particles[j] << " " << particles[j+1] << " " << particles[j+2] << " " ;
}
VTKfile <<std::endl<< "</DataArray> "<<std::endl
<< "</Points> "<<std::endl ;
if (nbDataPerParticle==8 ) {
VTKfile<< "<PointData Scalars=\"PhysicalValue\" > "<<std::endl
<< "<DataArray type=\"Float64\" Name=\"PhysicalValue\" format=\"ascii\">"<<std::endl ;
j = 0 ;
for(int i = 0 ; i< N; ++i, j+=nbDataPerParticle){
VTKfile << particles[j+3] << " " ;
}
VTKfile <<std::endl << "</DataArray>"<<std::endl ;
VTKfile << "<DataArray type=\"Float64\" Name=\"Potential\" format=\"ascii\">"<<std::endl ;
j = 0 ;
for(int i = 0 ; i< N; ++i, j+=nbDataPerParticle){
VTKfile << particles[j+4] << " " ;
}
VTKfile <<std::endl << "</DataArray>"<<std::endl ;
VTKfile<< "<DataArray type=\"Float64\" Name=\"Force\" NumberOfComponents=\"3\" format=\"ascii\"> "<<std::endl ;
j = 0 ;
for(int i = 0 ; i< N; ++i, j+=nbDataPerParticle){
VTKfile << particles[j+5] << " " << particles[j+6] << " " << particles[j+7] << " " ;
}
VTKfile <<std::endl<< "</DataArray> "<<std::endl;
}
else {
VTKfile << "<PointData Scalars=\"PhysicalValue\" > "<<std::endl
<< "<DataArray type=\"Float64\" Name=\"PhysicalValue\" format=\"ascii\">"<<std::endl ;
j = 0 ;
for(int i = 0 ; i< N; ++i, j+=nbDataPerParticle){
VTKfile << particles[j+3] << " " ;
}
VTKfile <<std::endl << "</DataArray>"<<std::endl ;
}
VTKfile << " </PointData>"<<std::endl
<< " <CellData>"<<" </CellData>"<<std::endl
<< " <Verts>"<<std::endl
<< " <DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">"<<std::endl ;
for(FSize i = 0 ; i< N; ++i){
VTKfile << i << " " ;
}
VTKfile<<std::endl << "</DataArray>" <<std::endl
<< "<DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">"<<std::endl ;
for(FSize i = 1 ; i< N+1; ++i){
VTKfile << i << " " ;
}
VTKfile<<std::endl << "</DataArray>"<<std::endl
<< " </Verts>"<<std::endl
<< "<Lines></Lines>"<<std::endl
<< "<Strips></Strips>"<<std::endl
<< "<Polys></Polys>"<<std::endl
<< "</Piece>"<<std::endl
<< "</PolyData>"<<std::endl
<< "</VTKFile>"<<std::endl;
} ;
#endif
//
#endif
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment