Commit 19996b55 authored by COULAUD Olivier's avatar COULAUD Olivier

Add tools to generate uniform and non uniform point distributions

parent 9f03db26
# check if compiling into source directories
STRING(COMPARE EQUAL "${CMAKE_SOURCE_DIR}" "${CMAKE_BINARY_DIR}" insource)
if(insource)
MESSAGE(FATAL_ERROR "${PROJECT_NAME} requires an out of source build. Goto scalfmm/Build and tapes cmake ../")
endif(insource)
project(Examples_scalfmm CXX)
set(EXECUTABLE_OUTPUT_PATH ${CMAKE_BUILD_TYPE})
ADD_DEFINITIONS( ${ScaLFMM_CXX_FLAGS})
# Find all code files
file(
GLOB_RECURSE
source_tests_files
./*.cpp
)
# Adding the project sources dir as an include dir
INCLUDE_DIRECTORIES(
${CMAKE_BINARY_DIR}/Src
${CMAKE_SOURCE_DIR}/Src
)
# Add execs - 1 cpp = 1 exec
foreach(exec ${source_tests_files})
get_filename_component(
execname ${exec}
NAME_WE
)
set(compile_exec "TRUE")
# Test Blas dependency
file(STRINGS "${exec}" lines_blas REGEX "@FUSE_BLAS")
if(lines_blas)
if( NOT ScalFMM_USE_BLAS )
MESSAGE( STATUS "This needs BLAS = ${exec}" )
set(compile_exec "FALSE")
endif()
endif()
# Test FFT dependency
file(STRINGS "${exec}" lines_fft REGEX "@FUSE_FFT")
if(lines_fft)
if( NOT ScalFMM_USE_FFT )
MESSAGE( STATUS "This needs FFT = ${exec}" )
set(compile_exec "FALSE")
endif()
endif()
# Test MPI dependency
file(STRINGS "${exec}" lines_mpi REGEX "@FUSE_MPI")
if(lines_mpi)
if( NOT ScalFMM_USE_MPI )
MESSAGE( STATUS "This needs MPI = ${exec}" )
set(compile_exec "FALSE")
endif()
endif()
# Dependency are OK
if( compile_exec )
add_executable(
${execname}
${exec}
)
SET_TARGET_PROPERTIES(${execname} PROPERTIES ENABLE_EXPORTS TRUE)
target_link_libraries(
${execname}
${scalfmm_lib}
${SCALFMM_LIBRARIES}
)
endif()
endforeach(exec)
/*
* genarateDistributions.cpp
*
* Created on: 23 mars 2014
* Author: Olivier Coulaud
*/
#include <iostream>
#include <fstream>
#include <sstream>
#include <string>
#include "Utils/FMath.hpp"
#include "Utils/FGenerateDistribution.hpp"
/// \file
//! \example generateDistributions.cpp
//! \brief generateDistributions: Driver to generate N points (non)uniformly distributed on a geometry
//! Options
//! -h to see the parameters available in this driver
//! -N The number of points in the distribution
//! Geometry
//! \arg \b -unitSphere uniform distribution on unit sphere
//! \arg -sphere uniform distribution on sphere of radius given by
//! \arg -radius R - default value for R is 2.0
//! \arg -ellipsoid non uniform distribution on an ellipsoid of aspect ratio given by
//! -ar a:b:c with a, b and c > 0
//!
//! -prolate ellipsoid with aspect ratio a:a:c
//! -ar a:a:c with c > a > 0
//! -plummer (Highly non unuiform) plummer distrinution (astrophysics)"
//! -radius R - default value 10.0"
//!
//! -filename name: generic name for files (without extension) and save data
//! with following format in name.xxx or name.bin in -bin is set
//! -visu save output in name.txt
//!
//! \b Physical values
//! -charge generate physical values between -1 and 1 otherwise generate between 0 and 1
void genDistusage() {
std::cout << "Driver to generate N points (non)uniformly distributed on a geometry"
<< std::endl;
std::cout << "Options "<< std::endl
<< " -h to see the parameters " << std::endl
<< " -N The number of points in the distribution " << std::endl
<< std::endl
<< " Distribution " << std::endl
<< " Uniform on " << std::endl
<< " -unitSphere uniform distribution on unit sphere" <<std::endl
<< " -sphere uniform distribution on sphere of radius given by" <<std::endl
<< " -radius R - default value for R is 2.0" <<std::endl
<< " -prolate ellipsoid with aspect ratio a:a:c" <<std::endl
<< " -ar a:a:c with c > a > 0" <<std::endl<<std::endl
<< " Non Uniform on " << std::endl
<< " -ellipsoid non uniform distribution on an ellipsoid of aspect ratio given by" <<std::endl
<< " -ar a:b:c with a, b and c > 0" <<std::endl
<< " -plummer (Highly non unuiform) plummer distrinution (astrophysics)"<<std::endl
<< " -radius R - default value 10.0" <<std::endl
<< " Physical values" <<std::endl
<< " -charge generate physical values between -1 and 1 otherwise generate between 0 and 1 " <<std::endl<<std::endl
<< " Output " << std::endl
<< " -filename name: generic name for files (without extension) and save data" <<std::endl
<< " with following format in name.xxx or name.bin in -bin is set" <<std::endl
<< " -visu save output in name.txt" <<std::endl;
}
int main(int argc, char ** argv){
//
if(FParameters::existParameter(argc, argv, "-h") ||FParameters::existParameter(argc, argv, "-help") ){
genDistusage() ;
exit(-1);
}
const FReal eps = 0.001 ;
const int NbPoints = FParameters::getValue(argc,argv,"-N", 20);
const std::string genericFileName(FParameters::getStr(argc,argv,"-filename", "unifPointDist"));
FReal BoxWith ;
//
// Allocation
//
FReal * particles ;
particles = new FReal[4*NbPoints] ;
memset(particles,0,4*NbPoints*sizeof(FReal));
//
// Generate physical values
//
FReal phyVal, sum,a,b ;
if(FParameters::existParameter(argc, argv, "-charge")){
a= 2.0 ; b = -1.0 ;
}
else {
a= 1.0 ; b = 0.0 ;
}
sum = 0.0 ;
int j = 3 ;
for(int i = 0 ; i< NbPoints; ++i, j+=4){
phyVal = a*getRandom() +b ;
sum += phyVal ;
particles[j] = phyVal ;
}
std::cout << "Sum physical value "<< sum << " Mean Value " << sum/NbPoints<<std::endl ;
//
// Point generation
//
if(FParameters::existParameter(argc, argv, "-unitSphere")){
unifRandonPointsOnUnitSphere(NbPoints, particles) ;
BoxWith = 2.0 ;
}
else if(FParameters::existParameter(argc, argv, "-sphere")){
const FReal Radius = FParameters::getValue(argc,argv,"-radius", 2.0);
unifRandonPointsOnSphere(NbPoints, Radius,particles) ;
BoxWith = 2.0*Radius ;
}
else if(FParameters::existParameter(argc, argv, "-prolate")){
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 ;
if(A != B){
std::cerr << " A /= B in prolate sllipsoide A =B. Your aspect ratio: "<< aspectRatio<<std::endl;
}
std::cout << "A: "<<A<<" B "<< B << " C: " << C<<std::endl;
unifRandonPointsOnProlate(NbPoints,A,C,particles);
BoxWith = C;
}
else if(FParameters::existParameter(argc, argv, "-ellipsoid")){
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;
nonunifRandonPointsOnElipsoid(NbPoints,A,B,C,particles);
BoxWith = FMath::Max( A,FMath::Max( B,C)) ;
}
else if(FParameters::existParameter(argc, argv, "-plummer")){
const FReal Radius = FParameters::getValue(argc,argv,"-radius", 10.0);
unifRandonPlummer(NbPoints, Radius, sum, particles) ;
BoxWith = 2.0*Radius ;
}
else {
std::cout << "Bad geometry option"<< std::endl;
exit(-1) ;
}
if(FParameters::existParameter(argc, argv, "-visu")){
std::ofstream file( genericFileName + ".txt", std::ofstream::out);
if(!file) {
std::cout << "Cannot open file."<< std::endl;
exit(-1) ;
} //
//
// Export data in cvs format
//
std::cout << "Writes in CVS format (visualization) in file "<< genericFileName + ".txt" <<std::endl ;
exportCVS( file, NbPoints, particles) ;
//
// Export data in vtk format
//
}
//
// Generate file for ScalFMM Loader
//
std::ofstream outfile( genericFileName + ".fma", std::ofstream::out);
if(!outfile) {
std::cout << "Cannot open file."<< std::endl;
exit(-1) ;
}
BoxWith += eps ;
std::cout << "Writes in FMA format in file "<< genericFileName + ".fma" <<std::endl ;
std::cout << " Points are in a cube of size "<< BoxWith << " Centered in the Origin"<<std::endl;
//
outfile << NbPoints << " " << BoxWith << " 0.0 0.0 0.0 " << std::endl;
j=0;
for(int i = 0 ; i< NbPoints; ++i, j+=4){
outfile << particles[j] << " " << particles[j+1] << " " << particles[j+2] << " " << particles[j+3] <<std::endl;
}
//
delete particles ;
//
return 1;
}
#ifndef FGENERATEDISTRIBUTION_HPP
#define FGENERATEDISTRIBUTION_HPP
#include <cstdlib>
#include <ctime>
#include <iostream>
#include <fstream>
#include <algorithm>
//
#include "Utils/FMath.hpp"
#include "Utils/FParameters.hpp"
/** return a random number between 0 and 1 */
//double getRandom() {
// return drand48();
//} ;
void initRandom() {
srand( static_cast<unsigned int>(time(0))) ;
} ;
FReal getRandom() {
// return static_cast<FReal>(rdand48());
return static_cast<FReal>(rand()/FReal(RAND_MAX));
} ;
//! \fn unifRandonPointsOnUnitSphere(const int N , FReal * points)
//! \brief Generate N points uniformly distributed on the unit sphere
//!
//! \param N the number of points uniformly randomly sample on the unit sphere
//! \param points array of size 4*N and stores data as follow x,y,z,0,x,y,z,0....
//! \example generateDistributions.cpp
void unifRandonPointsOnUnitSphere(const int N , FReal * points) {
FReal u, v, theta, phi, sinPhi ;
//
initRandom() ;
int j = 0 ;
for (int i = 0 ; i< N ; ++i, j+=4) {
//
u = getRandom() ; v = getRandom() ;
theta = FMath::FTwoPi*u ;
phi = FMath::ACos(2*v-1);
sinPhi = FMath::Sin(phi);
//
points[j] = FMath::Cos(theta)*sinPhi ;
points[j+1] = FMath::Sin(theta)*sinPhi ;
points[j+2] = 2*v-1 ;
//
}
};
//! \fn nonunifRandonPointsOnElipsoid(const int N , const FReal &a, const FReal &b, const FReal &c, FReal * points)
//! \brief Generate N points non uniformly distributed on the ellipsoid 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....
//!
void nonunifRandonPointsOnElipsoid(const int N , const FReal &a, const FReal &b, const FReal &c, FReal * points) {
//
FReal u, v , cosu ;
int j =0 ;
for (int i = 0 ; i< N ; ++i, j+=4) {
u = getRandom() ; v = getRandom() ;
u = FMath::FPi*u - FMath::FPiDiv2; v = FMath::FTwoPi*v - FMath::FPi;
cosu = FMath::Cos(u) ;
points[j] = a*cosu*FMath::Cos(v) ;
points[j+1] = b*cosu*FMath::Sin(v) ;
points[j+2] = c*FMath::Sin(u) ;
}
};
//! \fn nonunifRandonPointsOnElipsoid(const int N , const FReal &a, const FReal &c, FReal * points)
//! \brief Generate N points uniformly distributed on the ellipsoid of aspect ratio a:a:c
//!
//! \param N the number of points
//! \param a the x 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....
//!
void unifRandonPointsOnProlate(const int N , const FReal &a, const FReal &c, FReal * points){
//
FReal u, w,v ,ksi ;
FReal e = (a*a*a*a)/(c*c*c*c) ;
bool isgood = false;
int j =0 , cpt =0 ;
//
for (int i = 0 ; i< N ; ++i, j+=4) {
// Select a random point on the prolate
do {
cpt++ ;
u = getRandom() ; v = getRandom() ;
u = 2.0*u - 1.0; v = FMath::FTwoPi*v;
w =FMath::Sqrt(1-u*u) ;
points[j] = a*w*FMath::Cos(v) ;
points[j+1] = a*w*FMath::Sin(v) ;
points[j+2] = c*u ;
// Accept the position ?
ksi = a*getRandom() ;
// std::cout << "Gradf "<< points[j]*points[j] + points[j+1] *points[j+1] +e*points[j+2] *points[j+2] << std::endl;
isgood = (points[j]*points[j] + points[j+1] *points[j+1] +e*points[j+2] *points[j+2] < ksi*ksi );
} while (isgood);
}
std::cout.precision(4);
std::cout << "Total tested points: "<< cpt << " % of rejected points: "<<100*static_cast<FReal>(cpt-N)/cpt << " %" <<std::endl;
} ;
//! \fn unifRandonPointsOnSphere(const int N , const FReal R, FReal * points)
//! \brief Generate N points uniformly distributed on the sphere of radius R
//!
//! \param N the number of points uniformly randomly sample on the sphere
//! \param R the radius of the sphere
//! \param points array of size 4*N and stores data as follow x,y,z,0,x,y,z,0....
//!
void unifRandonPointsOnSphere(const int N , const FReal R, FReal * points) {
//
unifRandonPointsOnUnitSphere(N , points) ;
int j =0 ;
for (int i = 0 ; i< N ; ++i, j+=4) {
points[j] *= R ;
points[j+1] *= R ;
points[j+2] *= R ;
}
};
//! \fn FReal plummerDist(int & cpt, const FReal &R)
//! \brief Radial Plummer distribution
//!
//! \param cpt : counter to know how many random selections we need to obtain a radius less than R
//! \param R : Radius of the sphere that contains the particles
//! @return Return the radius according to the Plummer distribution either double type or float type
//!
FReal plummerDist(int & cpt, const FReal &R) {
//
FReal radius ,u ;
do {
//radius = getRandom() ;
u = FMath::pow (getRandom() , 2.0/3.0) ;
radius = FMath::Sqrt (u/(1.0-u));
cpt++;
if(radius <=R){
// std::cout << radius << " " <<std::endl;
return static_cast<FReal>(radius);
}
} while (true);
}
//! \fn void unifRandonPlummer(const int N , const FReal R, const FReal M, FReal * points)
//! \brief Build N points following the Plummer distribution
//! First we construct N points uniformly distributed on the unit sphere. Then the radius in construct according to the Plummr distribution.
//!
//! \param N the number of points following the Plummer distribution
//! \param R the radius of the sphere that contains all the points
//! \param M the total mass of all the particles inside the Sphere or radius R
//! \param points array of size 4*N and stores data as follow x,y,z,0,x,y,z,0....
//!
void unifRandonPlummer(const int N , const FReal R, const FReal M, FReal * points) {
//
unifRandonPointsOnUnitSphere(N , points) ;
//
FReal r , rm= 0.0;
// FReal Coeff = 3.0*M/(4.0*FMath::FPi*R*R*R) ;
//am1 = 0 ;//1/FMath::pow(1+R*R,2.5);
int cpt = 0 ;
for (int i = 0,j=0 ; i< N ; ++i, j+=4) {
// u \in []
r = plummerDist(cpt,R) ;
rm = std::max(rm, r);
points[j] *= r ;
points[j+1] *= r ;
points[j+2] *= r ;
}
std::cout << "Total tested points: "<< cpt << " % of rejected points: "
<<100*static_cast<FReal>(cpt-N)/cpt << " %" <<std::endl;
} ;
//! \fn void exportCVS(std::ofstream& file, const int N, const FReal * particles )
//! \brief Export particles in CVS Format
//!
//! Export particles in CVS Format as follow
//! x , y , z , physicalValue
//! It is useful to plot the distribution with paraView
//!
void exportCVS(std::ofstream& file, const int N, const FReal * particles ){
int j = 0;
file << " x , y , z, q " <<std::endl;
for(int i = 0 ; i< N; ++i, j+=4){
file << particles[j] << " , " << particles[j+1] << " , " << particles[j+2] << " , " << particles[j+3] <<std::endl;
}
}
void exportVTK(std::ofstream& file, const int N, const FReal * particles ){
int j = 0;
file << " x , y , z, q " <<std::endl;
for(int i = 0 ; i< N; ++i, j+=4){
file << particles[j] << " , " << particles[j+1] << " , " << particles[j+2] << " , " << particles[j+3] <<std::endl;
}
}
#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