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.
_______________________________________________
Mailing list: https://launchpad.net/~yade-users
Post to : yade-users@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
