Exporting a packing geometry made in yade in a stl format

Asked by Rioual on 2020-02-17

Hello,

I would like to export a packing geometry of spheres created in yade to a .stl format in order to open it
 and read it in a second stage with comsol.
How can I do that ??

Thanks,

Vincent

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Jan Stránský
Solved:
2020-02-21
Last query:
2020-02-21
Last reply:
2020-02-17
Jan Stránský (honzik) said : #1

Hello,

If it is possible (I don't know comsol), use the spherical data (center + radius).

STL [1] is a file format for triangulated surfaces. So if you need to use STL:
- triangulate the spherical packing
- save triangles in a .stl file [1]

cheers
Jan

[1] https://en.wikipedia.org/wiki/STL_(file_format)

Rioual (francois-rioual-v) said : #2

Hello Jan,

Thanks for your reply but:

- How do you triangulate the spherical packing ??

- "save triangles": do you mean i have to build by hands
 the stl file as explained in wikipedia i.e. save as a "built up"
 text file ???

Best wishes,

V.

Robert Caulk (rcaulk) said : #3

gmsh is one option [1]

[1]http://gmsh.info/

Jan Stránský (honzik) said : #4

> - How do you triangulate the spherical packing ??

depending on the meaning of a question:

1) triangulate it such that it suits your needs in comsol

2)
def triangulateOneSphere(b):
   ...
   return triangles
triangles = [triangulateOneSphere(b) for b in O.bodies]

triangulateOneSphere: google "sphere triangulation", depends on how fine you want/need the triangulation

> - "save triangles": do you mean i have to build by hands the stl file as explained in wikipedia i.e. save as a "built up" text file ???

yes. There might be some python stl library for easier processing

Before continuing:
- why do you want to use comsol (what is the background)? is there an alternative, like Paraview?
- is there another file format / method you can use more naturally (than triangulating spheres)? Like [1]?

cheers
Jan

https://www.comsol.com/forum/thread/211861/input-geometry-from-coordinates

Rioual (francois-rioual-v) said : #5

Dear Jan,

- My problem is related to the problem of evaporation dynamics in a porous material made with Yade so Comsol is indeed relevant.
- Maybe i will have to find another way to treat the problem from the coordinates and radii of the spheres calculated in Yade if making a stl file is too tough....

Thank you for your input,

V.

Best Jan Stránský (honzik) said : #6

well, it is not tough (MWE below, probably could be simplified by a meshing/triangulation library)
I just found it unnecessary without enough information
cheers
Jan

###
import stl # pip install stl

class Triangle:
   def __init__(self,p1,p2,p3):
      self.points = [Vector3(p) for p in (p1,p2,p3)]
   def scale(self,factor):
      self.points = [p*factor for p in self.points]
      return self
   def shift(self,v):
      self.points = [p+v for p in self.points]
      return self
   def refine(self): # split 1 triangle into 4
      p1,p2,p3 = self.points
      p4 = (.5*(p1+p2)).normalized()
      p5 = (.5*(p2+p3)).normalized()
      p6 = (.5*(p3+p1)).normalized()
      return [
         Triangle(p1,p4,p6),
         Triangle(p2,p5,p4),
         Triangle(p3,p6,p5),
         Triangle(p4,p5,p6),
      ]
   def toStlFacet(self):
      minieigen2stl = lambda v: stl.Vector3d(*tuple(v))
      p1,p2,p3 = self.points
      normal = (p2-p1).cross(p3-p1).normalized()
      normal = minieigen2stl(normal)
      vertices = [minieigen2stl(p) for p in self.points]
      return stl.Facet(normal,vertices)

def triangulateUnitSphere(): # octahedron, you can use whatever instead, e.g. icosahedron
   xm,xp = [(v,0,0) for v in (-1,+1)]
   ym,yp = [(0,v,0) for v in (-1,+1)]
   zm,zp = [(0,0,v) for v in (-1,+1)]
   return [
      Triangle(xp,yp,zp),
      Triangle(xm,yp,zp),
      Triangle(xp,ym,zp),
      Triangle(xm,ym,zp),
      Triangle(xp,yp,zm),
      Triangle(xm,yp,zm),
      Triangle(xp,ym,zm),
      Triangle(xm,ym,zm),
   ]

def triangulateOneSphere(sph,smooth=1):
   ret = triangulateUnitSphere()
   #
   for ismooth in range(smooth):
      ret2 = []
      for t in ret:
         ret2.extend(t.refine())
      ret = ret2
   #
   c = sph.state.pos
   r = sph.shape.radius
   for t in ret:
      t.scale(r).shift(c)
   #
   return ret

def triangulateSpheres(smooth=1):
   ret = []
   for b in O.bodies:
      if not isinstance(b.shape,Sphere): continue
      ret.extend(triangulateOneSphere(b,smooth))
   return ret

def spheres2stl(fName,smooth=1):
   triangles = triangulateSpheres(smooth)
   triangles = [t.toStlFacet() for t in triangles]
   solid = stl.Solid()
   solid.facets = triangles
   with open(fName,"w") as f:
      solid.write_ascii(f)

O.bodies.append((
   sphere((0,0,0),1),
   sphere((5,4,3),2),
   sphere((0,1,4),3),
))

spheres2stl("sphs.stl")
spheres2stl("sphs-smooth.stl",smooth=3)
###

Rioual (francois-rioual-v) said : #7

Thanks Jan Stránský, that solved my question.