generateDistributions.cpp 10.7 KB
Newer Older
1 2 3 4 5 6 7
/*
 * genarateDistributions.cpp
 *
 *  Created on: 23 mars 2014
 *      Author: Olivier Coulaud
 */

8
#include <algorithm>
9 10 11 12
#include <iostream>
#include <fstream>
#include <sstream>
#include <string>
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 "Files/FGenerateDistribution.hpp"
18
#include "Files/FFmaGenericLoader.hpp"
19
#include "Files/FExportWriter.hpp"
20

21
#include "Utils/FParameterNames.hpp"
22

23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77
/**
 * \file
 *
 * \brief Generates 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 -fout name: generic name for files (with extension) and save data with
 *                   following format in name.fma or name.bfma in -bin is set"
 * \param -fvisuout Filename 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 in unit cube
 * \param -cube uniform distribution in a cube
 *     \arg -size LX:LY:LZ - default value for R is 1.0:1.0: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 -ball uniform distribution in ball 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 -size a:b:c with a, b and c > 0
 * \param -prolate ellipsoid with aspect ratio a:a:c given by
 *     \arg -size 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</b>
 *
 *    generateDistributions -prolate -size 2:2:4   -N 20000 -fout prolate
 *
 *  or
 *
 *    generateDistributions -cuboid 2:2:4 -N 100000 -fout cuboid.bfma  -fvisuout cuboid.vtp -charge  -zeromean
 *
 */
78

79

80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99
namespace Param {
    const FParameterNames Ellipsoid
    = {{"-ellipsoid"}, "non uniform distribution on an ellipsoid of aspect ratio given by -size a:b:c with a, b and c > 0"};
    const FParameterNames UnitCube
    = {{"-unitCube"}, "uniform distribution on unit cube"};
    const FParameterNames Cube
    = {{"-cuboid"}, "uniform distribution on rectangular cuboid of size -size a:b:c - default values are 1.0:1.0:2.0 "};
    const FParameterNames UnitSphere
    = {{"-unitSphere"}, "uniform distribution on unit sphere"};
    const FParameterNames Ball
    = {{"-ball"}, "uniform distribution in a ball of radius given by -radius R - default value for R is 2.0"};
    const FParameterNames Sphere
    = {{"-sphere"}, "uniform distribution on sphere of radius given by -radius R - default value for R is 2.0"};
    const FParameterNames Prolate
    = {{"-prolate"}, "ellipsoid with aspect ratio a:a:c given by -size a:a:c with c > a > 0"};
    const FParameterNames Plummer
    = {{"-plummer"}, "(Highly non uniform) plummer distribution (astrophysics) -radius R - default value 10.0"};
    const FParameterNames Size
    = {{"-size"}, "Size of the geometry a:b:c - default values are 1.0:1.0:2.0"};
    const FParameterNames Radius
100
    = {{"-radius"}, "used to specified the radius of the sphere and the plummer distribution or R - default value for R is 2.0"};
101 102 103 104 105 106 107 108 109 110 111 112 113 114 115
    const FParameterNames Charge
    = {{"-charge"}, "generate physical values between -1 and 1 otherwise generate between 0 and 1"};
    const FParameterNames ZM
    = {{"-zeromean"}, "the average of the physical values is zero"};
    const FParameterNames EL
    = {{"-extraLength"}, "-extraLength value extra length to add to the boxWidth"};
}

#define getParamV(name, default)                                \
    FParameters::getValue(argc,argv,(name).options,(default))

#define getParamS(name, default)                                \
    FParameters::getStr(argc,argv,(name).options,(default))

int main(int argc, char ** argv){
116

117 118 119 120 121 122 123 124 125 126 127 128
    FHelpDescribeAndExit(
        argc, argv,
        ">> Driver to generate N points (non)uniformly distributed on a given geometry.\n"
        "Options  \n"
        "   -help       to see the parameters    ",
        FParameterDefinitions::OutputFile, FParameterDefinitions::NbParticles,
        FParameterDefinitions::OutputVisuFile,
        Param::UnitCube, Param::Cube,      Param::UnitSphere, Param::Sphere,
        Param::Radius,   Param::Ellipsoid, Param::Prolate,    Param::Plummer,
        Param::Ball,
        Param::Charge,   Param::ZM,        Param::EL,    Param::Size
        );
129

130
    using FReal = double;
131

132
    const FSize NbPoints = getParamV(FParameterDefinitions::NbParticles, FSize(20000));
133

134
    FReal extraRadius = 0.000;
135
    FReal BoxWith = 0.0;
136
    FPoint<FReal> Centre(0.0, 0.0,0.0);
137

138 139 140 141 142 143 144 145 146 147 148 149 150
    // Allocate particle array
    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 sum = 0;
    FReal a = 1.0;
    FReal b = 0.0;
    if(FParameters::existParameter(argc, argv, "-charge")){
        a = 2.0; b = -1.0;
    }
151

152 153 154 155
    for(int i = 0, j = 3 ; i< NbPoints; ++i, j+=4){
        particles[j] = a * getRandom<FReal>() + b;
        sum += particles[j] ;
    }
156

157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234
    if(FParameters::existParameter(argc, argv, "-zeromean")){
        FReal rm = FReal(sum) / FReal(NbPoints) ;
        sum -= static_cast<FReal>(NbPoints) * rm;
        for(int i = 0, j = 3 ; i< NbPoints; ++i, j+=4){
            particles[j] -= rm ;
        }
    }

    std::cout << "Physical value sum: " << sum
              << " mean: " << sum / FReal(NbPoints)
              << std::endl;

    // Read arguments
    // Radius
    const FReal Radius = getParamV(Param::Radius,  2.0);
    // Aspect ratio
    std::string aspectRatio = getParamS(Param::Size, "1:1:2");
    std::replace(aspectRatio.begin(), aspectRatio.end(), ':', ' ');
    FReal A, B, C;
    std::stringstream(aspectRatio) >> A >> B >> C;

    // Point  generation
    if(FParameters::existParameter(argc, argv, "-unitCube")) {
        unifRandomPointsInCube<FReal>(NbPoints, 1, 1, 1, particles);
        Centre.setPosition(0.5,0.5,0.5);
        BoxWith = 1.0;
        std::cout << "Unit cube "<< std::endl;
    }
    else if(FParameters::existParameter(argc, argv, "-ball")) {
        unifRandomPointsInBall<FReal>(NbPoints, Radius, particles);
        BoxWith = 2.0 * Radius;
        std::cout << "Ball radius: " << Radius << std::endl;
    }
    else if(FParameters::existParameter(argc, argv, "-cuboid")) {
        unifRandomPointsInCube(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")) {
        unifRandomPointsOnSphere<FReal>(NbPoints, 1.0, particles);
        BoxWith = 2.0;
    }
    else if(FParameters::existParameter(argc, argv, "-sphere")) {
        unifRandomPointsOnSphere(NbPoints, Radius, particles);
        BoxWith = 2.0 * Radius;
        std::cout << "Sphere radius: " << Radius << std::endl;
    }
    else if(FParameters::existParameter(argc, argv, "-prolate")) {
        if(A != B){
            std::cerr << " A != B in prolate ellipsoid. Your aspect ratio: "
                      << aspectRatio << std::endl;
        }
        std::cout << "Prolate A: " << A << " B: " << B << " C: " << C << std::endl;
        unifRandomPointsOnProlate(NbPoints, A, C, particles);
        BoxWith = 2.0 * C;
    }
    else if(FParameters::existParameter(argc, argv, "-hyperpara")) {
        unifRandomPointsOnHyperPara(NbPoints, A, B, C, particles);
        BoxWith = 2.0 * FMath::Max(A, FMath::Max(B, C));
        std::cout << "Hyperpara "<< A << ":"<< B<<":"<<C<<std::endl;
        std::cout << "BoxWith: " << BoxWith << std::endl;

    }
    else if(FParameters::existParameter(argc, argv, "-ellipsoid")){
        nonunifRandomPointsOnElipsoid(NbPoints, A, B, C, particles);
        BoxWith =  2.0 * FMath::Max(A, FMath::Max(B, C));
        std::cout << "Ellipsoid " << A << ":" << B << ":" << C << std::endl;
    }
    else if(FParameters::existParameter(argc, argv, "-plummer")){
        unifRandomPlummer(NbPoints, Radius, particles);
        BoxWith = 2.0 * Radius;
        std::cout << "Plummer radius: " << Radius << std::endl;
    }
    else {
        std::cout << "Bad geometry option"<< std::endl;
        exit(-1);
235
    }
236

COULAUD Olivier's avatar
COULAUD Olivier committed
237
    /////////////////////////////////////////////////////////////////////////
238
    //                                           Save data
COULAUD Olivier's avatar
COULAUD Olivier committed
239
    /////////////////////////////////////////////////////////////////////////
240
    //  Generate FMA file for FFmaGenericLoader<FReal> Loader
241 242 243 244 245 246 247 248 249 250
    if(FParameters::existParameter(argc, argv, "-extraLength")){
        extraRadius = FParameters::getValue(argc, argv, "-extraLength",  0.0);
        BoxWith += 2 * extraRadius;
    }
    const std::string name(getParamS(FParameterDefinitions::OutputFile, "unifPointDist"));
    std::cout << "Write "<< NbPoints <<" particles to '" << name << "'" << std::endl;
    FFmaGenericWriter<FReal> writer(name);
    writer.writeHeader(Centre, BoxWith, NbPoints, *ppart);
    writer.writeArrayOfParticles(ppart, NbPoints);
    std::cout << "End of writing" <<std::endl;
251

252
    //  Generate  file for visualization
253 254 255 256
//    if(FParameters::existParameter(argc, argv, FParameterDefinitions::OutputVisuFile.options)){
//        std::string outfilename(FParameters::getStr(argc,argv,FParameterDefinitions::OutputFile.options,   "output.vtp"));
//        driverExportData(outfilename, particles , NbPoints,loader.getNbRecordPerline() );
//    }
257 258 259 260 261 262
    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;
263
}