generateDistributions.cpp 10.8 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12
/*
 * genarateDistributions.cpp
 *
 *  Created on: 23 mars 2014
 *      Author: Olivier Coulaud
 */


#include <iostream>
#include <fstream>
#include <sstream>
#include <string>
COULAUD Olivier's avatar
COULAUD Olivier committed
13
//
14
#include "Utils/FGlobal.hpp"
15
#include "Utils/FMath.hpp"
COULAUD Olivier's avatar
COULAUD Olivier committed
16
#include "Utils/FPoint.hpp"
17
#include "Utils/FGenerateDistribution.hpp"
COULAUD Olivier's avatar
COULAUD Olivier committed
18 19
#include "Files/FFmaGenericLoader.hpp"

20 21 22 23 24 25 26
//
/// \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
//!
COULAUD Olivier's avatar
COULAUD Olivier committed
27
//!   Uniform : cube, cuboid, sphere, prolate,
28 29 30 31 32 33 34 35
//!
//!   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>
COULAUD Olivier's avatar
COULAUD Olivier committed
36
//!     \param   -help (-h)      to see the parameters available in this driver
COULAUD Olivier's avatar
COULAUD Olivier committed
37
//!     \param  -N     The number of points in the distribution (default 20000)
COULAUD Olivier's avatar
COULAUD Olivier committed
38
//!     \param   -filename name: generic name for files (with extension) and save data
39
//!                  with following format in name.fma or name.bfma in -bin is set"
COULAUD Olivier's avatar
COULAUD Olivier committed
40
//!      \param  -visufmt format for the visu file (vtk, vtp, cvs or cosmo). vtp is the default
COULAUD Olivier's avatar
COULAUD Olivier committed
41
//!      \param -extraLength   value    extra length to add to the boxWidth (default 0.0)
42 43
//!
//!  <b> Geometry arguments:</b>
COULAUD Olivier's avatar
COULAUD Olivier committed
44 45 46
//!      \param  -unitCube uniform distribution on unit cube
//!      \param  -cube uniform distribution on a cube
//!          \arg         -length  R - default value for R is 2.0
47 48
//!      \param  -unitSphere uniform distribution on unit sphere
//!      \param  -sphere  uniform distribution on  sphere of radius given by
49
//!          \arg         -radius  R - default value for R is 2.0
50 51 52 53 54 55 56 57 58 59
//!        \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
COULAUD Olivier's avatar
COULAUD Olivier committed
60 61
//!         \param -zeromean  the average of the physical values is zero
//!
62
//!
63
//! \b examples
64
//!
COULAUD Olivier's avatar
COULAUD Olivier committed
65 66 67
//!   generateDistributions -prolate -ar 2:2:4   -N 20000 -filename prolate
//!
//! or
68
//!
COULAUD Olivier's avatar
COULAUD Olivier committed
69
//!   generateDistributions  -cuboid 2:2:4 -filename cuboid  -visufmt vtp -charge  -zeromean
70
//!
COULAUD Olivier's avatar
COULAUD Olivier committed
71
//
72 73 74

//
//
75 76 77 78 79 80
void genDistusage() {
	std::cout << "Driver to generate N points (non)uniformly distributed on a given geometry"
			<< std::endl;
	std::cout <<	 "Options  "<< std::endl
			<<     "   -help       to see the parameters    " << std::endl
			<<     "   -N       The number of points in the distribution    " << std::endl
COULAUD Olivier's avatar
COULAUD Olivier committed
81
			<<     "   -extraLength   value    extra length to add to the boxWidth"<< std::endl
82 83 84
			<<    std::endl
			<<     "    Distributions   " << std::endl
			<<     "        Uniform on    " << std::endl
COULAUD Olivier's avatar
COULAUD Olivier committed
85 86 87
			<<     "             -unitCube  uniform distribution on unit cube" <<std::endl
			<<     "             -cuboid  uniform distribution on rectangular cuboid of size  a:a:c" <<std::endl
			<<     "                     -lengths   a:a:c - default value for R is 1.0:1.0:2.0" <<std::endl
88 89 90 91 92 93 94 95 96 97 98 99
			<<     "             -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
COULAUD Olivier's avatar
COULAUD Olivier committed
100
			<<     "             -zeromean  the average of the physical values is zero		" <<std::endl<<std::endl
101 102
			<<     "     Output " << std::endl
			<<     "             -filename name: generic name for files (without extension) and save data" <<std::endl
103
			<<     "                     with following format in name.fma or name.bfma in -bin is set" <<std::endl
104 105
			<<     "             -visufmt  vtk, vtp, cosmo or cvs format " <<std::endl;
}
106

107

108 109
int main(int argc, char ** argv){
	//
COULAUD Olivier's avatar
COULAUD Olivier committed
110
	if(FParameters::existParameter(argc, argv, "-h")||FParameters::existParameter(argc, argv, "-help")|| (argc < 3 )){
111 112 113
		genDistusage() ;
		exit(-1);
	}
COULAUD Olivier's avatar
COULAUD Olivier committed
114
	 FReal       extraRadius = 0.000 ;
COULAUD Olivier's avatar
COULAUD Olivier committed
115
	const int NbPoints  = FParameters::getValue(argc,argv,"-N",   20000);
116
	const std::string genericFileName(FParameters::getStr(argc,argv,"-filename",   "unifPointDist"));
COULAUD Olivier's avatar
COULAUD Olivier committed
117 118
	FReal BoxWith = 0.0;
	FPoint Centre(0.0, 0.0,0.0);
119 120 121 122 123 124
	//
	// Allocation
	//
	FReal * particles ;
	particles = new FReal[4*NbPoints] ;
	memset(particles,0,4*NbPoints*sizeof(FReal));
COULAUD Olivier's avatar
COULAUD Olivier committed
125 126
	FmaBasicParticle *ppart = (FmaBasicParticle*)(&particles[0]);

127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144
	//
	// 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 ;
	}
COULAUD Olivier's avatar
COULAUD Olivier committed
145 146 147 148 149 150 151 152
	if(FParameters::existParameter(argc, argv, "-zeromean")){
		FReal  rm = sum/NbPoints ; sum = 0.0 ;
		j = 3 ;
		for(int i = 0 ; i< NbPoints; ++i, j+=4){
			particles[j]    -= rm ;
			sum              += particles[j]   ;
		}
	}
153 154 155 156
	std::cout << "Sum physical value "<< sum << "   Mean Value " << sum/NbPoints<<std::endl ;
	//
	// Point  generation
	//
COULAUD Olivier's avatar
COULAUD Olivier committed
157 158 159
	if(FParameters::existParameter(argc, argv, "-unitCube")){
		unifRandonPointsOnUnitCube(NbPoints, particles) ;
		Centre.setPosition(0.5,0.5,0.5);
COULAUD Olivier's avatar
COULAUD Olivier committed
160
		BoxWith = 1.0 ;
COULAUD Olivier's avatar
COULAUD Olivier committed
161 162 163 164 165 166 167 168
	}
	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) ;
COULAUD Olivier's avatar
COULAUD Olivier committed
169
		BoxWith = FMath::Max(A,FMath::Max(B,C) );
COULAUD Olivier's avatar
COULAUD Olivier committed
170
		FReal halfBW = BoxWith*0.5;
COULAUD Olivier's avatar
COULAUD Olivier committed
171 172 173
		Centre.setPosition(halfBW,halfBW,halfBW);
	}
	else if(FParameters::existParameter(argc, argv, "-unitSphere")){
174
		unifRandonPointsOnUnitSphere(NbPoints, particles) ;
COULAUD Olivier's avatar
COULAUD Olivier committed
175
		BoxWith = 2.0 ;
176 177 178 179 180 181 182 183 184 185 186 187 188 189 190
	}
	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;
		}
191
		std::cout << "A: "<<A<<" B "<< B << " C: " << C<<std::endl;
192
		unifRandonPointsOnProlate(NbPoints,A,C,particles);
COULAUD Olivier's avatar
COULAUD Olivier committed
193
		BoxWith =  2.0*C;
194 195 196 197 198 199 200 201 202
	}
	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);
COULAUD Olivier's avatar
COULAUD Olivier committed
203
		BoxWith =  2.0*FMath::Max( A,FMath::Max( B,C)) ;
204
	}
205 206 207 208 209
	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 ;
	}
210

211 212 213 214
	else {
		std::cout << "Bad geometry option"<< std::endl;
		exit(-1) ;
	}
COULAUD Olivier's avatar
COULAUD Olivier committed
215 216 217
    /////////////////////////////////////////////////////////////////////////
	//                                           Save data
    /////////////////////////////////////////////////////////////////////////
218
	//
COULAUD Olivier's avatar
COULAUD Olivier committed
219
	//  Generate FMA file for FFmaGenericLoader Loader
220
	//
COULAUD Olivier's avatar
COULAUD Olivier committed
221 222 223 224
	if(FParameters::existParameter(argc, argv, "-extraLength")){
		extraRadius  = FParameters::getValue(argc,argv,"-extraLength",  0.0);
		BoxWith += 2*extraRadius ;
	}
COULAUD Olivier's avatar
COULAUD Olivier committed
225 226 227 228 229 230
	bool binaryMode = false;
	std::string name(genericFileName);

	if(  FParameters::existParameter(argc, argv, "-bin")){
		binaryMode = true;
		name += ".bfma";
COULAUD Olivier's avatar
COULAUD Olivier committed
231
	}
COULAUD Olivier's avatar
COULAUD Olivier committed
232 233
	else {
		name += ".fma";
234
	}
COULAUD Olivier's avatar
COULAUD Olivier committed
235 236 237
	FFmaGenericWriter  writer(name,binaryMode) ;
	writer.writeHeader(Centre,BoxWith, NbPoints, *ppart) ;
	writer.writeArrayOfParticles(ppart, NbPoints);
238
	//
COULAUD Olivier's avatar
COULAUD Olivier committed
239
	//  Generate  file for visualization
240 241 242 243 244 245 246 247 248 249 250 251 252 253 254
	//
	if(FParameters::existParameter(argc, argv, "-visufmt")){
		std::string visufile(""), fmt(FParameters::getStr(argc,argv,"-visufmt",   "vtp"));
		if( fmt == "vtp" ){
			visufile = genericFileName + ".vtp" ;
		}
		else	if( fmt == "vtk" ){
			visufile = genericFileName + ".vtk" ;
		}
		else if( fmt == "cosmo" ){
			visufile = genericFileName + ".cosmo" ;
		}
		else {
			visufile =   genericFileName + ".csv" ;
		}
255
		std::ofstream file( visufile.c_str(), std::ofstream::out);
256 257 258 259 260
		if(!file) {
			std::cout << "Cannot open file."<< std::endl;
			exit(-1)	;
		}	//
		//
COULAUD Olivier's avatar
COULAUD Olivier committed
261
		// Export data in VTK format
262 263 264
		//
		if( fmt == "vtp" ){
			std::cout << "Writes in XML VTP format  (visualization) in file "<< visufile <<std::endl ;
COULAUD Olivier's avatar
COULAUD Olivier committed
265
			exportVTKxml( file, particles, NbPoints)  ;
266 267 268
		}
		else		if( fmt == "vtk" ){
			std::cout << "Writes in VTK format  (visualization) in file "<< visufile <<std::endl ;
COULAUD Olivier's avatar
COULAUD Olivier committed
269
			exportVTK( file, particles, NbPoints)  ;
270 271 272
		}
		else if( fmt == "cosmo" ){
			std::cout << "Writes in COSMO format  (visualization) in file "<< visufile <<std::endl ;
COULAUD Olivier's avatar
COULAUD Olivier committed
273
			exportCOSMOS( file, particles, NbPoints)  ;
274 275 276
		}
		else {
			std::cout << "Writes in CVS format  (visualization) in file "<<visufile<<std::endl ;
COULAUD Olivier's avatar
COULAUD Olivier committed
277
			exportCVS( file, particles, NbPoints)  ;
278 279 280
		}
	}
	//
281 282 283 284 285
	delete particles ;

	//
	return 1;
}