yade-users team mailing list archive

[Jan Stránský <question681014@xxxxxxxxxxxxxxxxxxxxx>]

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

    Status: Open => Answered

Jan Stránský proposed the following answer:
Hi Xuesong,

nice task.
Since the packing is porous and in reality there is no "top surface", triangulation of particles would nor help.. 

A few options (depending on your needs and requirements):
- enlarge the particles "a bit". Then you could get "real surface" either analytically or by triangulation
- using some kind of a potential function of the spheres and assume the surface at certain level of this potential (see below)
- ...

### a MWE for the "potential" approach
# method params (you can play with the values)
level = .98 # closer to 1 -> closer to actual surface
factor = 2 # higher -> closer to actual surface
# dimensions of the packing
dim = 10.
# using only [dd,dim-dd] interval
dd = 2.5
# createa grid of the surface
n = 50

# "potential" of a body. =1 inside or on the surface, decreases to 0 for increasing distance from surface
def potential1(x,y,z,body,factor=1.):
   xyz = Vector3(x,y,z)
   p = body.state.pos
   r = body.shape.radius
   d = (xyz-p).norm() - r # distance from surface
   f = factor/r
   if d < 0: # inside
      ret = 1
   else: # outside
      ret = exp(-f*d) # could be other functions, like 1/x or so...
   return ret

# sum of potentials of all bodies. Very inoptimal (!) using all bodies, a few "boudary bodies" would be enough
def potential(x,y,z,factor=1.):
   return sum(potential1(x,y,z,b,factor) for b in O.bodies)

# for given x,y find z such that potential=level
# The closer is level to 1, the closer is the final z to actual surface
# using bisection method with initial boundaries z1,z2
def findZ(level,x,y,z1,z2,factor=1.):
   for i in range(20): # TODO a more sohisticated end condition
      z = .5*(z1+z2)
      l = potential(x,y,z,factor)
      if l < level:
         z2 = z
      else:
         z1 = z
   return z

# create some packing
pred = pack.inAlignedBox((0,0,0),(dim,dim,dim))
sphs = pack.randomDensePack(pred,1,rRelFuzz=.5,spheresInCell=100)
O.bodies.append(sphs)

surf = [[None for y in range(n+1)] for x in range(n+1)]
for ix in range(n+1):
   x = dd + ix*(dim-2*dd)/n
   for iy in range(n+1):
      y = dd + iy*(dim-2*dd)/n
      z = findZ(level,x,y,.5*dim,2*dim,factor)
      surf[ix][iy] = (x,y,z)
# create facets to see the result
for ix in range(n):
   for iy in range(n):
      v1 = surf[ix+0][iy+0]
      v2 = surf[ix+1][iy+0]
      v3 = surf[ix+0][iy+1]
      v4 = surf[ix+1][iy+1]
      f1 = facet((v1,v2,v3))
      f2 = facet((v2,v3,v4))
      O.bodies.append((f1,f2))

# export to vtk for Paraview
from yade import export
vtk = export.VTKExporter("surf")
vtk.exportSpheres()
vtk.exportFacets()
###

cheers
Jan

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