How to export the edges of a Voronoi tessellation? (no fluid case)

Asked by Gael Lorieul on 2019-04-30

Hi all,

I've seen a picture of a Voronoi Tesselation obtained with Yade on the wiki [1, Fig.2] where the edges of the Voronoi cells are visible. I want to do that (visualizing the Voronoi cells), but I don't really see how to obtain such a result with Yade…

What I've seen:
  - The example code [3] for `TesselationWrapper` does not present explicit calculation and exporting of the tesselation (tesselation is calculated implicitly through call to `computeVolume()`, exporting this information is not mentioned at all).
  - This answered question [4] seems to be very similar to what I am looking for, except that it concerns a Fluid-Particles interaction, whereas I only work with particles (no fluid). The solution of [4] uses the `saveVtk()` function of class `TwoPhaseFlowEngine`. The implementation of the logic seems to be in member function `saveVtk()` of class `DFNBoundingSphere` [7]. In other words, the logic is not whithin the `Tesselation` or `TesselationWrapper` classes but within a module for Fluid-particles interaction. Also note that `TwoPhaseFlowEngine::saveVtk()` has now been commented out in the source code [5] although it can still be found in the documentation [6].
  - Python class `TesselationWrapper` [2] does not seem to propose a function for explicit calculation of Voronoi tesselation.
  - C++ class `TesselationWrapper` contains `computeTesselation()` member function which does perform an explicit calculation of Voronoi tesselation, but it is not part of the `YADE_CLASS_BASE_DOC_*()` macro of the class `TesselationWrapper`, so it is not accessible from the Python interface.
  - The only way to compute the tesselation from Python that I have seen is to call the `TesselationWrapper.computeVolumes()` member function, but it does not output the shape of the cells, only their volume (from what I understand). Moreover (and still from what I understand) volumes can only be obtained one-by-one through successive calls to `TesselationWrapper.volume(sphereId)` for each sphere of the simulation.

Is the feature of exporting the geometry of Voronoi Tesselation cells implemented/supported by Yade?
If not, what would be the best way for me to achieve the result that I want?

Thanks for your answers ;)

Gaël

PS: what is the mysterious "Alpha" in the `getAlpha*()` functions? (for instance `getAlphaFaces()`)

[1] https://www.yade-dem.org/wiki/Triangulation
https://www.yade-dem.org/wiki/File:Cell_volume.png

[2]
https://yade-dev.gitlab.io/trunk/yade.wrapper.html#yade.wrapper.TesselationWrapper

[3]
https://yade-dev.gitlab.io/trunk/tutorial-more-examples-fast.html#tesselationwrapper

[4]
https://answers.launchpad.net/yade/+question/656046

[5]
./pkg/pfv/TwoPhaseFlowEngine.hpp:438
./pkg/pfv/TwoPhaseFlowEngine.hpp:613

[6]
https://yade-dev.gitlab.io/trunk/yade.wrapper.html#yade.wrapper.FlowEngine.saveVtk

[7]
./pkg/pfv/DFNFlow.cpp

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
2019-05-30
Last reply:
2019-05-16
Launchpad Janitor (janitor) said : #1

This question was expired because it remained in the 'Open' state without activity for the last 15 days.

Gael Lorieul (glorieul64) said : #2

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