yade-users team mailing list archive

[Gael Lorieul <question680569@xxxxxxxxxxxxxxxxxxxxx>]

Question #680569 on Yade changed:
https://answers.launchpad.net/yade/+question/680569

Gael Lorieul posted a new comment:
Hi again,

It seemed to me that the Voronoi tessellation implemented in Yade is
certainly very good to calculate shear and stresses on homogeneous
packings, but its purpose was not to analyse the geometrical properties
of a packing (including packing factor), nor to export the Voronoi
tessellation, which is what I wanted to do. So instead I am now using
the Voro++ (C++) library which was designed for extracting such
geometrical information and can handle inhomogeneous packings.

I've made a solution involving:
 1. A Yade script exporting the particles' positions and radius
 2. A small C++ code that parses it and calls the Voro++ library to generate the tessellation, which is exported to a file.
 3. A gnuplot script that can visualize that file.

I don't think it is the ideal solution, but it works well and could be
used as a basis for a better version of it. I copy this solution below.

Possibly a better solution could be to include an interface or wrapper
to Voro++ from Yade?

Enjoy!

Gaël


File `01_generatePacking.py`:

    import csv
    from yade import pack

    def saveParticles(particles, fileName):
        with open(fileName, 'w') as csvFile:
            fieldnames = ['id', 'positionX', 'positionY','positionZ','radius']
            csvWriter = csv.DictWriter(csvFile, delimiter=',', quotechar='"',
            fieldnames=fieldnames)
            csvWriter.writeheader()
            for ptcl in particles:
                pos = ptcl.state.pos
                radius = ptcl.shape.radius
                csvWriter.writerow({'id':ptcl.id, 'positionX': pos[0],
                                    'positionY': pos[1], 'positionZ': pos[2],
                                    'radius':radius})


    # GENERATE PARTICLES
    sp = pack.SpherePack()
    sp.makeCloud((-1,-1,-1),(1,1,1),rMean=.2, rRelFuzz=0.8, seed=1)
    sp.toSimulation()

    # SAVE PARTICLES
    particles = O.bodies # Assumes all bodies are particles
    saveParticles(particles, 'particles.csv')

    exit(0) # Exit explicitly otherwise Yade displays the Yade command
prompt


File `02_generateTesselation.cpp` (to compile as [1], the "csv.h" dependencies can be found here [2])

    #include <array>
    #include <iostream>
    #include <string>
    #include <sstream>
    #include <vector>
    #include "voro++.hh"
    #include "dependencies/fast-cpp-csv-parser/csv.h"

    struct Particle
    {
    public:
        Particle() = default;
        Particle(int id, std::array<double,3> position, double radius)
            :id(id), position(position), radius(radius)
        { /*Do nothing*/ }
   
        std::string toString()
        {
            std::ostringstream os;
            os << "\tpos = (" << position[0] << "," << position[1] << "," << position[2] << "), ";
            os << "radius = " << radius << std::endl;
            return os.str();
        }
   
        int id;
        std::array<double,3> position;
        double radius;
    };

    int main()
    {
        // READ PARTICLES' POSITION AND SIZES
        io::CSVReader<5> in("particles.csv");
        in.read_header(io::ignore_extra_column, "id", "positionX", "positionY", "positionZ", "radius");
        int id; double posX, posY, posZ, radius;
        std::vector<Particle> particles = {};
        while(in.read_row(id, posX, posY, posZ, radius))
            { particles.push_back( Particle(id, {posX,posY,posZ}, radius) ); }

        // GENERATE TESSELATION
        const std::vector<double> SWCorner = {-1.0,-1.0,-1.0};
        const std::vector<double> NECorner = {+1.0,+1.0,+1.0};
        const std::vector<int>    nbBlocks = {6,6,6};
        constexpr bool isPeriodicX = false, isPeriodicY = false, isPeriodicZ = false;
        const int nbPtclPerBlock = 8; // Is an estimation (initial allocation size)
        voro::container_poly con(SWCorner[0], NECorner[0], SWCorner[1], NECorner[1],
                                 SWCorner[2], NECorner[2], nbBlocks[0], nbBlocks[1], nbBlocks[2],
                                 isPeriodicX, isPeriodicY, isPeriodicZ, nbPtclPerBlock);
        for(auto ptcl : particles)
            { con.put(ptcl.id, ptcl.position[0], ptcl.position[1], ptcl.position[2], ptcl.radius); }
        con.draw_cells_gnuplot("tesselation.gp");
    }

File 03_visualizeTesselation.py:

    import csv
    import numpy as np
    import PyGnuplot as pg


    def parseParticleFile(fileName):
        particles = []
        with open(fileName, 'r') as csvParticles:
            csvReader = csv.DictReader(csvParticles)
            for entry in csvReader:
                newParticle = {'id': entry['id'],
                               'center': (entry['positionX'], entry['positionY'], entry['positionZ']),
                               'radius' : entry['radius']}
                particles.append(newParticle)
        return particles

    def getCmd_newSphere(particle):
        ctr = particle["center"]
        rad = particle["radius"]
        u, v = np.mgrid[0:2*np.pi:4j, 0:np.pi:3j] # 10, 5
        cmd = ( str(ctr[0]) + " + " + str(rad) + "*cos(u)*cos(v), "
              + str(ctr[1]) + " + " + str(rad) + "*sin(u)*cos(v), "
              + str(ctr[2]) + " + " + str(rad) + "*sin(v) with pm3d" )
        return cmd


    # MAIN SCRIPT

    particlesFile = 'particles.csv'
    tesselationFile = 'tesselation.gp'
    particles = parseParticleFile(particlesFile)

    cmdParticles = "splot [-pi:pi][-pi/2:pi/2] "
    for particle in particles:
        cmdParticles += getCmd_newSphere(particle) + ", "
    cmdParticles = cmdParticles[:-2] # Remove trailing ", "

    cmdTesselation = "'" + tesselationFile + "' with lines lw 3 lc rgb '#000000'"
#print("cmdTesselation = ", cmdTesselation)

    pg.c("set term wxt persist") # So that the window doesn't close directly after finished plotting
    pg.c("unset colorbox")       # Do not display colorbar
    pg.c("set nokey")            # Do not display a legend
    pg.c("set parametric")
    pg.c("set view 60")          # Angle of camera position (0°=from top, 90°=from side, 180°=from bottom)
    pg.c("set view equal xyz")   # Axes have the same scaling
    pg.c("set xrange[-2 : 2]")
    pg.c("set yrange[-2 : 2]")
    pg.c("set zrange[-2 : 2]")
    pg.c("set pm3d depthorder")  # To avoid z-order glitches within pm3d objects
    pg.c("set hidden3d front")   # To draw lines in a way that respects z-order
    pg.c("set pm3d lighting primary 0.4 specular 0.0") # So that particles have some shadding
    pg.c("set palette defined (1 'red', 2 'red')") # Define color in colorbar: in that case paint everything in red
    pg.c(cmdParticles + ", " + cmdTesselation)

Dependencies:
  - Dependencies of `01_generatePacking.py`:
      - Yade (of course!)
  - Dependencies of `02_generateTesselation.cpp`:
      - Voro++ http://math.lbl.gov/voro++/
      - fast-cpp-csv-parser : https://github.com/ben-strasser/fast-cpp-csv-parser/
  - Dependencies of `03_visualizeTesselation.py`:
      - PyGnuplot (get it from pip)
      - gnuplot

Then to run:
    yade 01_generatePacking.py
    ./02_generateTesselation.run
    python3 03_visualizeTesselation.py

---End of instructions---

I place all code within this e-mail under CC0 license.
https://creativecommons.org/choose/zero/

[1]
g++ --std=c++14 -I /usr/local/include/voro++/ 02_generateTesselation.cpp -o 02_generateTesselation.run -lpthread -L /usr/local/lib/ -lvoro++

(Adapt `-I` and `-L` options depending on your installation)


[2]
https://github.com/ben-strasser/fast-cpp-csv-parser/blob/master/csv.h

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.