fuseDistributions.cpp 10.3 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
/**
 * \file
 * \brief Fuses FMA files to create a new distribution
 *
 * \author Quentin Khan
 * \copyright ScalFmm 2016 INRIA
 * \copyright [CeCILL-C licence](http://www.cecill.info)
 *
 *
 */

#include <algorithm>
#include <fstream>
#include <sstream>
#include <string>
#include <vector>


#include "Files/FFmaGenericLoader.hpp"

21 22 23 24 25 26 27
/**
 * \brief Print programe usage
 *
 * \param progname The name of the program.
 *
 * \note Use `argv[0]` as progname.
 */
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
void usage(const std::string& progname) {
    std::size_t start = progname.find_last_of('/');
    std::string name = progname.substr(start+1);
    std::cout <<
        "usage: " << name <<
        " --file [[-s scale] [-c cx:cy:cz] [-g gx:gy:gz]] filename"
        " -fout output_file"
        " [--extra-length length]"
        "\n"
        "\n"
        "Fuses multiple particle distributions into a bigger one."
        "\n"
        "\n"
        "Options:\n"
        "  -fout output_file\n"
        "    The output file name, must hase .bfma or .fma extension\n"
        "\n"
        "  --file [opts] filename [opts]\n"
        "    Add a .fma or .bfma distibution file. Multiple files may be specified by\n"
        "    adding more --file options. 'opts' is a combination of:\n"
        "      -s scale\n"
        "        Scale the distribution by 'scale' factor.\n"
        "      -c cx:cy:cz\n"
        "        Center the distribution at given coordinates. cx, cy and cz are\n"
        "        floating point numbers.\n"
        "      -g gx:gy:gz\n"
        "        Duplicate the distribution inside a grid of gx by gy by gz dimensions.\n"
        "        gx, gy and gz are integers. The grid center is governed by the -c \n"
        "        option.\n"
57 58 59
        "      -r rx:ry:rz\n"
        "        Rotate the distribution around its x, y and z axes. The rotation \n"
        "        center is the distribution center. rx, ry and rz are in radians.\n"
60 61 62 63 64 65 66 67 68 69 70
        "\n"
        "  --extra-length length\n"
        "    Length to be added to the final box width.\n"
        "\n"
        "  --help\n"
        "    Print this message."
        "\n"
        ;
}


71 72
using FReal = double;

73 74 75
/**
 * \brief Particle object
 */
76
struct Particle {
77 78
    FPoint<FReal> pos; ///< Particle position in space
    FReal val;         ///< Particle physical value
79 80
};

81 82 83 84 85
/**
 * \brief Distribution options
 *
 * Holds the tranformations that must be applied to a particle distribution.
 */
86 87 88 89 90 91 92 93 94 95 96
struct distribution {
    /// Distribution filename
    std::string filename = "";
    /// Distribution offset from center
    FPoint<FReal> offset = {0,0,0};
    /// Distribution rotation around its center
    FPoint<FReal> rot = {0,0,0};
    /// Distribution scale factor
    FReal scale = 1;
};

97 98 99
/**
 * \brief Main program parameters
 */
100
struct parameters {
101
    /// Output distribution file
102
    std::string output_filename;
103
    /// Particles distributions to combine
104
    std::vector<distribution> distributions;
105
    /// Length to add to the output distribution box
106 107 108
    FReal extra_length;
};

109 110 111 112 113 114 115 116
/**
 * \brief Helper method to parse the arguments that apply to one distribution.
 *
 * \param args A vector holding the command line arguments
 * \param i    Index of the current argument in args
 *
 * \return A vector of distribution objects
 */
117 118 119 120 121 122 123 124 125 126 127
std::vector<distribution> subparse_file(const std::vector<std::string>& args, std::size_t& i) {
    std::stringstream sstr;
    // Grid layout
    unsigned int gx = 1, gy = 1, gz = 1;
    // Final distributions, one per grid part
    std::vector<distribution> distributions;
    // Distributions options
    distribution dist;

    while(i < args.size() && args[i] != "--file") {
        sstr.clear();
128
        if(false) {// false for alignment of the other cases
129
        } else if(args[i] == "-s") {
130
            // Parse scale argument
131 132 133 134
            ++i;
            sstr.str(args.at(i));
            sstr >> dist.scale;
        } else if (args[i] == "-c") {
135
            // Parse offset from center argument
136
            ++i;
137
            char c; // Used to discard the ':' from argument format
138 139 140
            sstr.str(args.at(i));
            sstr >> dist.offset[0] >> c >> dist.offset[1] >> c >> dist.offset[2];
        } else if(args[i] == "-g") {
141
            // Parse grid layout
142
            ++i;
143
            char c; // Used to discard the ':' from argument format
144 145
            sstr.str(args.at(i));
            sstr >> gx >> c >> gy >> c >> gz;
146
        } else if(args[i] == "-r") {
147
            // Parse rotation angle
148 149 150 151
            ++i;
            char c; // Used to discard the ':' from argument format
            sstr.str(args.at(i));
            sstr >> dist.rot[0] >> c >> dist.rot[1] >> c >> dist.rot[2];
152 153 154 155 156 157 158 159 160 161
        } else {
            if(dist.filename != "") {
                --i;
                break;
            }
            dist.filename = args[i];
        }
        ++i;
    }

162 163 164 165
    if(i < args.size() && args[i] == "--file") {
        --i;
    }

166 167
    // If a grid layout was specified, duplicate the distribution for each cell
    // in the grid
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
    if(gx > 1 || gy > 1 || gz > 1) {
        // Compute offset of lowest left grid offset
        FFmaGenericLoader<FReal> loader(dist.filename);
        FReal box_width = loader.getBoxWidth() * dist.scale;
        dist.offset[0] -= (gx-1) * box_width / 2;
        dist.offset[1] -= (gy-1) * box_width / 2;
        dist.offset[2] -= (gz-1) * box_width / 2;

        // Create one distribution for each part of the grid layout
        for(unsigned int x = 0; x < gx; ++x) {
            for(unsigned int y = 0; y < gy; ++y) {
                for(unsigned int z = 0; z < gz; ++z) {
                    distribution tmp_dist = dist;
                    tmp_dist.offset[0] += x * box_width;
                    tmp_dist.offset[1] += y * box_width;
                    tmp_dist.offset[2] += z * box_width;
                    distributions.push_back(tmp_dist);
                }
            }
        }
    } else {
        distributions.push_back(dist);
    }

    return distributions;
}

195 196 197 198 199
/**
 * \brief Main argument parse function
 *
 * \param args A vector holding the command line arguments
 */
200 201
parameters parse(const std::vector<std::string>& args) {
    parameters params;
202
    std::stringstream sstr; // Use to convert strings to numbers
203
    for(std::size_t i = 1; i < args.size(); ++i) {
204 205 206
        if(args[i] == "--help") {
            usage(args[0]);
            exit(0);
207
        } else if(args[i] == "--file") {
208
            // Parse a new distribution specification
209 210 211 212 213
            ++i;
            auto ds = subparse_file(args, i);
            params.distributions.insert(params.distributions.end(),
                                        ds.begin(), ds.end());
        } else if(args[i] == "--extra-length") {
214
            // Parse extra length parameter
215 216 217 218
            ++i;
            sstr.str(args.at(i));
            sstr >> params.extra_length;
        } else if(args[i] == "-fout") {
219
            // Parse output file name
220 221 222 223 224 225 226 227 228
            ++i;
            params.output_filename = args.at(i);
        } else {
            std::cerr << "Unknown or misplaced parameters: " << args[i] << '\n';
        }
    }
    return params;
}

229 230 231 232 233 234
/**
 * \brief Rotate a particle around the origin
 *
 * \param p    particle to rotate
 * \param dist distribution transformation specification
 */
235 236 237 238 239 240 241 242 243 244 245 246 247 248 249
void rotate(Particle& p, const distribution& dist) {
    // Rotate around x axis
    if(dist.rot[0] > 1e-5 || dist.rot[0] < -1e-5) {
        FReal alpha = dist.rot[0];
        p.pos[1] = p.pos[1] * cos(alpha) - p.pos[2] * sin(alpha);
        p.pos[2] = p.pos[1] * sin(alpha) + p.pos[2] * cos(alpha);
    }
    // Rotate around y axis
    if(dist.rot[1] > 1e-5 || dist.rot[1] < -1e-5) {
        FReal alpha = dist.rot[1];
        p.pos[0] =  p.pos[0] * cos(alpha) + p.pos[2] * sin(alpha);
        p.pos[2] = -p.pos[0] * sin(alpha) + p.pos[2] * cos(alpha);
    }
    // Rotate around z axis
    if(dist.rot[2] > 1e-5 || dist.rot[2] < -1e-5) {
250
        FReal alpha = dist.rot[2];
251 252 253 254 255
        p.pos[0] = p.pos[0] * cos(alpha) - p.pos[1] * sin(alpha);
        p.pos[1] = p.pos[0] * sin(alpha) + p.pos[1] * cos(alpha);
    }
}

256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276


int main(int argc, char** argv) {
    auto params = parse({argv,argv+argc});

    // Fail early if output file raises an error
    FFmaGenericWriter<FReal> writer(params.output_filename);

    // Fuse particle distributions
    std::vector<Particle> particles;
    FReal axis_max = 0;

    for(distribution& dist : params.distributions) {
        // Load particles into array
        FFmaGenericLoader<FReal> loader(dist.filename);
        const std::size_t count = loader.getParticleCount();
        // Particle array: x1, y1, z1, val1, x2, y2...
        particles.reserve(particles.size() + count);

        FPoint<FReal> center = loader.getBoxCenter();

277
        // Load each particle and apply the distribution transformations
278 279 280 281 282 283 284 285
        Particle p;
        for(std::size_t i = 0; i < count; ++i) {
            loader.fillParticle(&p.pos, &p.val);
            // Move distribution center to origin
            p.pos -= center;
            // Scale distribution
            p.pos *= dist.scale;
            // Rotate distribution
286
            rotate(p, dist);
287 288 289 290 291 292 293 294 295 296 297 298 299
            // Move to new position
            p.pos += dist.offset;
            // Add particle to list
            particles.push_back(p);

            // Save particle x,y,z min/max to compute final box
            axis_max = std::max(std::abs(p.pos[0]), axis_max);
            axis_max = std::max(std::abs(p.pos[1]), axis_max);
            axis_max = std::max(std::abs(p.pos[2]), axis_max);
        }
    }

    // Write final distribution
300

301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330
    FPoint<FReal> center(0,0,0);
    // Compute final box width
    FReal box_width = 2 * (axis_max + params.extra_length);

    // Write header
    writer.writeHeader(center, box_width, particles.size(), 8, 4);

    // Write all particles

    // Buffer avoids duplicating particle vector
    std::vector<FReal> buffer;
    buffer.reserve(4*1024); // Avoid reallocations, size is a multiple of 4

    auto cur = particles.begin();
    auto sentinel = particles.end();

    // Fill and write buffer until we're done
    while(cur != sentinel) {
        buffer.clear();
        while(buffer.size() != buffer.capacity() && cur != sentinel) {
            buffer.push_back(cur->pos[0]);
            buffer.push_back(cur->pos[1]);
            buffer.push_back(cur->pos[2]);
            buffer.push_back(cur->val);
            ++cur;
        }
        writer.writeArrayOfReal(buffer.data(), 4, buffer.size()/4);
    }

}