Mentions légales du service

Skip to content
Snippets Groups Projects
generateDistributions.cpp 11.38 KiB
/*
 * genarateDistributions.cpp
 *
 *  Created on: 23 mars 2014
 *      Author: Olivier Coulaud
 */


#include <iostream>
#include <fstream>
#include <sstream>
#include <string>
//
#include "Utils/FGlobal.hpp"
#include "Utils/FMath.hpp"
#include "Utils/FPoint.hpp"
#include "Files/FGenerateDistribution.hpp"
#include "Files/FFmaGenericLoader.hpp"
#include "Files/FExportWriter.hpp"

#include "Utils/FParameterNames.hpp"

//
/// \file  generateDistributions.cpp
//!
//! \brief generateDistributions: Driver to generate N points (non)uniformly distributed on a given geometry
//!
//! The goal of this driver is to generate uniform or non uniform points on the following geometries
//!
//!   Uniform : cube, cuboid, sphere, prolate,
//!
//!   Non uniform : ellipsoid, prolate
//!
//!  You can set two kind of physical values depending of your problem. By default all values are between 0 and 1.
//!   If you select the argument -charge (see bellow) the values are between -1 and 1.
//!  The arguments available are
//!
//!  <b> General arguments:</b>
//!     \param   -help (-h)      to see the parameters available in this driver
//!     \param  -N     The number of points in the distribution (default 20000)
//!     \param   -filename name: generic name for files (with extension) and save data
//!                  with following format in name.fma or name.bfma in -bin is set"
//!      \param  -visufmt format for the visu file (vtk, vtp, cvs or cosmo). vtp is the default
//!      \param -extraLength   value    extra length to add to the boxWidth (default 0.0)
//!
//!  <b> Geometry arguments:</b>
//!      \param  -unitCube uniform distribution on unit cube
//!      \param  -cube uniform distribution on a cube
//!          \arg         -length  R - default value for R is 2.0
//!      \param  -unitSphere uniform distribution on unit sphere
//!      \param  -sphere  uniform distribution on  sphere of radius given by
//!          \arg         -radius  R - default value for R is 2.0
//!        \param   -ellipsoid non uniform distribution on  an ellipsoid of aspect ratio given by
//!              \arg          -ar a:b:c   with a, b and c > 0
//!         \param  -prolate ellipsoid with aspect ratio a:a:c  given by
//!                \arg             -ar a:a:c   with  c > a > 0
//!          \param   -plummer (Highly non uniform) plummer distribution (astrophysics)
//!                   \arg         -radius  R - default value 10.0"
//!
//!
//!  <b> Physical values argument:</b>
//!         \param -charge generate physical values between -1 and 1 otherwise generate between 0 and 1
//!         \param -zeromean  the average of the physical values is zero
//!
//!
//! \b examples
//!
//!   generateDistributions -prolate -ar 2:2:4   -N 20000 -filename prolate
//!
//! or
//!
//!   generateDistributions  -cuboid 2:2:4 -filename cuboid  -visufmt vtp -charge  -zeromean
//!



int main(int argc, char ** argv){
    const FParameterNames LocalOptionEllipsoid = {
        {"-ellipsoid"} ,
        " non uniform distribution on  an ellipsoid of aspect ratio given by \n  -ar a:b:c   with a, b and c > 0\n"
    };
    FHelpDescribeAndExit(argc, argv,
                         ">> Driver to generate N points (non)uniformly distributed on a given geometry.\n"
                         "Options  \n"
                         "   -help       to see the parameters    \n"
                         "   -N       The number of points in the distribution    \n"
                         "   -extraLength   value    extra length to add to the boxWidth\n"
                         "    Distributions   \n"
                         "        Uniform on    \n"
                         "             -unitCube  uniform distribution on unit cube\n"
                         "             -cuboid  uniform distribution on rectangular cuboid of size  a:b:c\n"
                         "                     -lengths   a:b:c - default values are 1.0:1.0:2.0\n"
                         "             -unitSphere  uniform distribution on unit sphere\n"
                         "             -sphere  uniform distribution on  sphere of radius given by\n"
                         "                     -radius  R - default value for R is 2.0\n"
                         "             -prolate ellipsoid with aspect ratio a:a:c\n"
                         "                     -ar a:a:c   with  c > a > 0\n"
                         "        Non Uniform on    \n"
                         "             -ellipsoid non uniform distribution on  an ellipsoid of aspect ratio given by\n"
                         "                     -ar a:b:c   with a, b and c > 0\n"
                         "             -plummer (Highly non unuiform) plummer distrinution (astrophysics)\n"
                         "                     -radius  R - default value 10.0\n"
                         "    Physical values\n"
                         "             -charge generate physical values between -1 and 1 otherwise generate between 0 and 1		\n"
                         "             -zeromean  the average of the physical values is zero		\n",
//                         "     Output \n"
 //                        "             -filename name: generic name for files (without extension) and save data\n"
  //                       "                     with following format in name.fma or name.bfma in -bin is set\n"
   //                      "             -visufmt  vtk, vtp, cosmo or cvs format.",
                          FParameterDefinitions::OutputFile,
                         FParameterDefinitions::NbParticles,FParameterDefinitions::OutputVisuFile,LocalOptionEllipsoid);


    
    typedef double FReal;
    FReal       extraRadius = 0.000 ;

    const FSize NbPoints  = FParameters::getValue(argc,argv,FParameterDefinitions::NbParticles.options,   FSize(20000));
    const std::string genericFileName(FParameters::getStr(argc,argv,FParameterDefinitions::OutputFile.options,   "unifPointDist"));

    FReal BoxWith = 0.0;
    FPoint<FReal> Centre(0.0, 0.0,0.0);
	//
	// Allocation
	//
	FReal * particles ;
	particles = new FReal[4*NbPoints] ;
	memset(particles,0,4*NbPoints*sizeof(FReal));
    FmaRWParticle<FReal, 4,4> *ppart = (FmaRWParticle<FReal, 4,4>*)(&particles[0]);

	//
	// 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<FReal>() +b  ;
		sum              += phyVal ;
		particles[j]       = phyVal ;
	}
	if(FParameters::existParameter(argc, argv, "-zeromean")){
        FReal  rm = FReal(sum)/FReal(NbPoints) ; sum = 0.0 ;
		j = 3 ;
		for(int i = 0 ; i< NbPoints; ++i, j+=4){
			particles[j]    -= rm ;
			sum              += particles[j]   ;
		}
	}
    std::cout << "Sum physical value "<< sum << "   Mean Value " << sum/FReal(NbPoints)<<std::endl ;
	//
	// Point  generation
	//
	if(FParameters::existParameter(argc, argv, "-unitCube")){
		unifRandonPointsOnUnitCube(NbPoints, particles) ;
		Centre.setPosition(0.5,0.5,0.5);
		BoxWith = 1.0 ;
	}
	else if(FParameters::existParameter(argc, argv, "-cuboid")){
		std::string  dd(":"),aspectRatio  = FParameters::getStr(argc,argv,"-lengths",  "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 ;
		unifRandonPointsOnCube(NbPoints, A,B,C,particles) ;
		BoxWith = FMath::Max(A,FMath::Max(B,C) );
		FReal halfBW = BoxWith*0.5;
		Centre.setPosition(halfBW,halfBW,halfBW);
		std::cout << "Cuboid "<< A << ":"<< B<<":"<<C<<std::endl;
	}
	else 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 ellipsoide A =B. Your aspect ratio: "<< aspectRatio<<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));
    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");
//		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 =  2.0*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) ;
	}
    /////////////////////////////////////////////////////////////////////////
	//                                           Save data
    /////////////////////////////////////////////////////////////////////////
	//
    //  Generate FMA file for FFmaGenericLoader<FReal> Loader
	//
	if(FParameters::existParameter(argc, argv, "-extraLength")){
		extraRadius  = FParameters::getValue(argc,argv,"-extraLength",  0.0);
		BoxWith += 2*extraRadius ;
	}
	std::string name(genericFileName);
	std::cout << "Write "<< NbPoints <<" Particles in file " << name << std::endl;
    FFmaGenericWriter<FReal>  writer(name) ;
	writer.writeHeader(Centre,BoxWith, NbPoints, *ppart) ;
	writer.writeArrayOfParticles(ppart, NbPoints);
	std::cout << "    End of writing "<<std::endl;

	//
	//  Generate  file for visualization
	//
    if(FParameters::existParameter(argc, argv, FParameterDefinitions::OutputVisuFile.options)){
        std::string visufile(FParameters::getStr(argc,argv,FParameterDefinitions::OutputVisuFile.options,   "output.vtp"));
         driverExportData(visufile, particles , NbPoints);
	}
	//
	delete [] particles ;

	//
	return 1;
}