yade-users team mailing list archive

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

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

Jan Stránský posted a new comment:
Generating pores "home-made"
Cheers
Jan

###
import scipy.spatial

# https://www2.mps.mpg.de/homes/daly/CSDS/t4h/tetra.htm
# https://math.stackexchange.com/questions/2820212/circumradius-of-a-tetrahedron
# filter out "degenerate" tetrahedrons
def delaunay2poreTetras(delaunay,limit=0.1):
    vertices = delaunay.points
    tetras = delaunay.simplices
    def quality(tetra):
        vs = v0,v1,v2,v3 = [Vector3(vertices[i]) for i in tetra]
        d1 = v1 - v0
        d2 = v2 - v0
        d3 = v3 - v0
        v = d1.dot(d2.cross(d3)) / 6
        v = abs(v) # volume
        a = d1.norm()
        b = d2.norm()
        c = d3.norm()
        A = (v2-v3).norm()
        B = (v1-v3).norm()
        C = (v1-v2).norm()
        aa,bb,cc = a*A, b*B, c*C
        den = (aa+bb+cc)*(aa+bb-cc)*(aa-bb+cc)*(-aa+bb+cc)
        den = abs(den)
        r = sqrt(den)/(24*v) # radius of circum sphere
        q = (9*sqrt(3)/8*v)**(1/3.)/r # quality in range [0,1]
        return q
    return [t for t in tetras if quality(t) > limit]

# some home-made inefficient (but working :-) iterative approach
def tetra2pore(tetra):
    bodies = [O.bodies[int(i)] for i in tetra]
    vertices = [b.state.pos for b in bodies]
    radii = [b.shape.radius for b in bodies]
    center = 0.25*sum(vertices,Vector3.Zero) # initial guess, center of mass of tetrahedron
    for i in range(100): # TODO stop condition
        ds = d1,d2,d3,d4 = [(center-v).norm()-r for v,r in zip(vertices,radii)] # distances from "touching"
        dMin,dMax = min(ds),max(ds)
        iMin,iMax = ds.index(dMin), ds.index(dMax)
        vMin,vMax = vertices[iMin],vertices[iMax]
        dirMin = (center - vMin).normalized()
        dirMax = (center - vMax).normalized()
        # based on min and max distance "push" the center a bit
        dd = (dMax - dMin) * 0.1
        center += dirMin * dd - dirMax * dd
    radius = (center - vertices[0]).norm() - radii[0]
    return center,radius


# testing packing, colored red
O.bodies.append((
    sphere((1.0,2.0,3.0),1.0,color=(1,0,0),mask=0b001),
    sphere((4.0,2.0,3.0),1.6,color=(1,0,0),mask=0b001),
    sphere((2.5,4.0,3.0),1.2,color=(1,0,0),mask=0b001),
    sphere((2.6,3.0,5.0),1.1,color=(1,0,0),mask=0b001),
    sphere((5.0,4.0,5.0),1.4,color=(1,0,0),mask=0b001),
    sphere((5.0,5.0,2.5),1.3,color=(1,0,0),mask=0b001),
    sphere((6.5,3.0,3.0),1.2,color=(1,0,0),mask=0b001),
    sphere((7.4,4.8,4.0),1.2,color=(1,0,0),mask=0b001),
))

centers = [b.state.pos for b in O.bodies]
delaunay = scipy.spatial.Delaunay(centers) # tetrahedralization
tetras = delaunay2poreTetras(delaunay) # only "good" tetrahedrons
pores = [tetra2pore(t) for t in tetras] # pores out of tetrahedrons

# pores, colored cyan
for c,r in pores:
    O.bodies.append(sphere(c,r,color=(0,1,1),mask=0b010))

# tetrahedrons as facets, for visualization
for ii in tetras:
    vs = [delaunay.points[i] for i in ii]
    v1,v2,v3,v4 = [Vector3(list(v)) for v in vs]
    O.bodies.append((
        facet((v1,v2,v3),mask=0b100),
        facet((v1,v2,v4),mask=0b100),
        facet((v1,v3,v4),mask=0b100),
        facet((v2,v3,v4),mask=0b100),
    ))

# improve sphere display quality
yade.qt.Gl1_Sphere.quality = 2
###

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