Commit 65e208d3 by Quentin Khan

### Merge branch 'master' into adaptive

parents 15453e60 16823939
This diff is collapsed.
 ... ... @@ -16,7 +16,12 @@ #ifndef FGENERATEDISTRIBUTION_HPP #define FGENERATEDISTRIBUTION_HPP // @author O. Coulaud /** * \file * \brief Distribution generation implementations * \author O. Coulaud */ #include #include ... ... @@ -27,250 +32,247 @@ #include "Utils/FMath.hpp" #include "Utils/FParameters.hpp" /** return a random number between 0 and 1 */ /** * \brief Seed the random number generator using current time */ void initRandom() { srand48( static_cast(time(nullptr))) ; } ; template FReal getRandom() { return static_cast(drand48()); //return static_cast(rand()/FReal(RAND_MAX)); } ; //! \fn unifRandonPointsOnUnitCube(const int N , FReal * points) //! \brief Generate N points uniformly distributed on the unit cube srand48(static_cast(time(nullptr))); } //! //! \param N the number of points uniformly randomly sample on the unit cube //! \param points array of size 4*N and stores data as follow x,y,z,0,x,y,z,0.... //! \example generateDistributions.cpp /** * \brief Generate a random number * \tparam FReal Floating point type * \return A random number in [0,1] */ template void unifRandonPointsOnUnitCube(const FSize N , FReal * points) { // initRandom() ; int j = 0; for (FSize i = 0 ; i< N ; ++i, j+=4) { // points[j] = getRandom() ; points[j+1] = getRandom() ; points[j+2] = getRandom() ; // } }; //! \fn unifRandonPointsOnCube(const int N , FReal * points) //! \brief Generate N points uniformly distributed on the cube of length R FReal getRandom() { return static_cast(drand48()); } //! //! \param N the number of points uniformly randomly sample on the unit cube //! \param Lx the the X-length of the cube //! \param Ly the the Y-length of the cube //! \param Lz the the Z-length of the cube //! \param points array of size 4*N and stores data as follow x,y,z,0,x,y,z,0.... //! \example generateDistributions.cpp /** * \brief Generate points uniformly inside a cuboid * * \tparam FReal Floating point type * * \param N the number of points uniformly randomly sample on the unit cube * \param Lx the the X-length of the cuboid * \param Ly the the Y-length of the cuboid * \param Lz the the Z-length of the cuboid * \param points array of size 4*N and stores data as follow x,y,z,0,x,y,z,0... */ template void unifRandonPointsOnCube(const FSize N , const FReal& Lx, const FReal &Ly, const FReal& Lz, FReal * points) { // unifRandonPointsOnUnitCube(N , points) ; FSize j =0 ; for (FSize i = 0 ; i< N ; ++i, j+=4) { points[j] *= Lx ; points[j+1] *= Ly ; points[j+2] *= Lz ; } }; //! \fn unifRandonPointsOnUnitSphere(const int N , FReal * points) void unifRandomPointsInCube(const FSize N, const FReal& Lx, const FReal& Ly, const FReal& Lz, FReal* points) { initRandom(); for(FSize i = 0, j = 0 ; i< N ; ++i, j+=4) { points[j] = getRandom() * Lx; points[j+1] = getRandom() * Ly; points[j+2] = getRandom() * Lz; } } //! \brief Generate N points uniformly distributed on the unit sphere /** * \brief Generate points uniformly inside a ball * * \tparam FReal Floating point type * * \param R the ball radius * \param points array of size 4*N and stores data as follow x,y,z,0,x,y,z,0... */ template void unifRandomPointsInBall(const FSize N, const FReal R, FReal* points) { initRandom(); //! //! \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 template void unifRandonPointsOnUnitSphere(const FSize N , FReal * points) { FReal u, v, theta, phi, sinPhi ; // initRandom() ; FSize j = 0 ; for (FSize 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) auto is_in_sphere = [&R](FReal* p) { return p[0]*p[0] + p[1]*p[1] + p[2]*p[2] < R*R; }; //! \brief Generate N points non uniformly distributed on the ellipsoid of aspect ratio a:b:c for(FSize i = 0, j = 0 ; i< N ; ++i, j+=4) { do { points[j] = (getRandom() - 0.5) * 2 * R; points[j+1] = (getRandom() - 0.5) * 2 * R; points[j+2] = (getRandom() - 0.5) * 2 * R; } while(! is_in_sphere(points + j)); } } //! //! \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.... //! /** * \brief Generate N points non uniformly distributed on the ellipsoid of aspect ratio a:b:c * * \tparam FReal Floating point type * * \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 void nonunifRandonPointsOnElipsoid(const FSize N , const FReal &a, const FReal &b, const FReal &c, FReal * points) { // FReal u, v , cosu ; FSize j =0 ; for (FSize 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) void nonunifRandomPointsOnElipsoid(const FSize N, const FReal& a, const FReal& b, const FReal& c, FReal* points) { FReal u, v, cosu; for (FSize i = 0, j = 0 ; i< N ; ++i, j+=4) { u = FMath::FPi() * getRandom() - FMath::FPiDiv2(); v = FMath::FTwoPi() * getRandom() - 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); } } //! \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.... //! /** * \brief Generate N points uniformly distributed on the ellipsoid of aspect ratio a:a:c * * \tparam FReal Floating point type * * \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.... */ template void unifRandonPointsOnProlate(const FSize 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; FSize j =0 , cpt =0 ; // for (FSize 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(cpt-N)/static_cast(cpt) << " %" <() - 1.0; v = FMath::FTwoPi() * getRandom(); 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(); 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(cpt-N) / static_cast(cpt) << " %" << std::endl; } //! \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.... //! /** * \brief Generate N points uniformly distributed on the hyperbolic paraboloid of aspect ratio a:b:c * * \tparam FReal Floating point type * * \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 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() - 1.0 ; v = 2.0*getRandom() - 1.0 ; points[j] = a*u ; points[j+1] = b*v ; points[j+2] = c*(u*u - v*v) ; void unifRandomPointsOnHyperPara(const FSize N, const FReal &a, const FReal &b, const FReal &c, FReal * points) { FReal u, v; for (FSize i = 0, j = 0 ; i< N ; ++i, j+=4) { u = 2.0 * getRandom() - 1.0; v = 2.0 * getRandom() - 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 //! //! \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.... //! /** * \brief Generate N points uniformly distributed on the sphere of radius R * * \tparam FReal Floating point type * * \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... */ template void unifRandonPointsOnSphere(const FSize N , const FReal R, FReal * points) { // unifRandonPointsOnUnitSphere(N , points) ; FSize j =0 ; for (FSize i = 0 ; i< N ; ++i, j+=4) { points[j] *= R ; points[j+1] *= R ; points[j+2] *= R ; } void unifRandomPointsOnSphere(const FSize N, const FReal R, FReal* points) { initRandom(); FReal u, v, theta, phi, sinPhi; for (FSize i = 0, j = 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 * R; points[j+1] = FMath::Sin(theta) * sinPhi * R; points[j+2] = (2*v - 1) * R; } }; //! \fn void 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 //! /** * \brief Radial Plummer distribution * * \tparam FReal Floating point type * * \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 The radius according to the Plummer distribution */ template FReal plummerDist(FSize cpt, const FReal &R) { // FReal radius ,u ; do { // u = FMath::pow (getRandom() , 2.0/3.0) ; radius = FMath::Sqrt (u/(1.0-u)); cpt++; if(radius <=R){ // std::cout << radius << " " <(radius); } } while (true); FReal plummerDist(FSize& cpt, const FReal &R) { FReal radius, u; while(true) { u = FMath::pow(getRandom(), 2.0/3.0); radius = FMath::Sqrt(u/(1.0-u)); cpt++; if(radius <= R) { return static_cast(radius); } } } //! \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.... /** * \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. * * \tparam FReal Floating point type * * \param N the number of points following the Plummer distribution * \param R the radius of the sphere that contains all the points * \param points array of size 4*N and stores data as follow x,y,z,0,x,y,z,0.... */ template void unifRandonPlummer(const FSize 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); FSize cpt = 0 ; for (FSize 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(cpt-N)/static_cast(cpt) << " %" <(N, 1, points); FReal r, rm = 0.0; FSize cpt = 0; for (FSize i = 0, j = 0 ; i< N ; ++i, j+=4) { 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(cpt-N) / static_cast(cpt) << " %" << std::endl; } // #endif
 ... ... @@ -17,6 +17,7 @@ // ==== CMAKE ==== // Keep in private GIT // @SCALFMM_PRIVATE // @FUSE_BLAS #include ... ...
 ... ... @@ -16,7 +16,7 @@ // Keep in private GIT // @SCALFMM_PRIVATE // @FUSE_BLAS #define DEBUG_SPHERICAL_M2L #define BLAS_SPHERICAL_COMPRESS #define BLAS_M2L_P ... ...
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