yade-users team mailing list archive

[Mithushan Soundaranathan <question701015@xxxxxxxxxxxxxxxxxxxxx>]

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

    Status: Answered => Open

Mithushan Soundaranathan is still having a problem:
Hi Robert,

Thank you for your help, I implemented your suggestion and defined as
PoreSize in the MWE. However, I get an empty vertices dataset.

I want track the pore size change over time during swelling (radius
increase) of the particle. Do I need include a flow in there over time,
or can still use it by define dead=1 in the engine? How do change the
name of the vertices.txt file so it differ for each iteration?

Best regards,
Mithu

Sorry for the code, here is the new MWE:
from __future__ import print_function
from yade import utils, plot, timing
from yade import pack

o = Omega()
o.dt = 1.0e-8

# Add material
mat1 = O.materials.append(LudingMat(frictionAngle=fr, density=rho, k1=k1, kp=kp, ks=ks, kc=kc, PhiF=PhiF1, G0 = 0.0))

# Spheres for compression and walls
sp=pack.SpherePack()
sp.makeCloud((-7.5*Diameter,-7.5*Diameter,-30*Diameter),(7.5*Diameter,7.5*Diameter,5.0*Diameter), rMean=Diameter/2.0,rRelFuzz=0.18,num=round(no_p))
sp.toSimulation(material=mat1)
walls=O.bodies.append(yade.geom.facetCylinder((0,0,0),radius=Tab_rad,height=Cyl_height,segmentsNumber=20,wallMask=6,material=mat1))


# Add engines
o.engines = [
  ForceResetter(),
  InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1.05),
                         Bo1_Wall_Aabb(),
                         Bo1_Facet_Aabb()
                         ]),
  InteractionLoop(
    [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=1.05),
     Ig2_Facet_Sphere_ScGeom(),
     Ig2_Wall_Sphere_ScGeom()],
    [Ip2_LudingMat_LudingMat_LudingPhys()],
    [Law2_ScGeom_LudingPhys_Basic()]
  ),
  NewtonIntegrator(damping=0.1, gravity=[0, 0, g]),
  FlowEngine(dead=1,label="flow"),
  PyRunner(command='checkForce()', realPeriod=1, label="fCheck"),
  #DeformControl(label="DefControl")
]


def checkForce():
    if O.iter < 3000000:
        return
    if unbalancedForce() > 1:
        return
    global upper_punch
    upper_punch=O.bodies.append(geom.facetCylinder((0,0,((-Cyl_height/2)+0.0001)+utils.aabbDim()[2]),Tab_rad-.00001,0,segmentsNumber=50,wallMask=1))
    for i in upper_punch:
        body= O.bodies[i]
        body.state.vel = (0,0,-0.02)
    global lower_punch
    lower_punch= O.bodies.append(geom.facetCylinder((0,0,(-Cyl_height/2)-0.0001),Tab_rad-.00001,0,segmentsNumber=50,wallMask=1))
    for n in lower_punch:
        body= O.bodies[n]
        body.state.vel = (0,0,0.02)
    O.engines = O.engines + [PyRunner(command='storeData()', iterPeriod=25000)]+ [PyRunner(command='saveData()', iterPeriod=100000)]
    fCheck.command = 'unloadPlate()'

def checkForce():
    if O.iter < 3000000:
        return
    if unbalancedForce() > 1:
        return
    fCheck.command = 'unloadPlate()'

def unloadPlate():
    if ((force_up > Comp_force_up) and (force_lp > Comp_force_lp)):
        for i in upper_punch:
            body= O.bodies[i]
            body.state.vel = (0,0,0.04)
        for n in lower_punch:
            body= O.bodies[n]
            body.state.vel = (0,0,-0.04)
        fCheck.command = 'stopUnloading()'

def stopUnloading():
    if pos_up[2]> pos_lp[2]+utils.aabbDim()[2]+0.0002: 
        for j in upper_punch: O.bodies.erase(j)
        fCheck.command = 'Savecheck()'
    
    
def Savecheck():    
    if save==0:
        o.engines = o.engines+[PyRunner(command='PoreSize()', iterPeriod=100000)]
        fCheck.dead = True # (!!!)
        storeData.dead=True
        saveData.dead=True
   

def PoreSize():
    flow.printVertices()


if save==0:
    read_filename='PH101_'+str(tab_porosity)+'_'+str(tab_height)+'mm.xml'
    o.load(read_filename)
    utils.loadVars('lower_punch')
    from yade.params.lower_punch import *
    for b in O.bodies:
        if isinstance(b.shape, Sphere):
            #b.state.blockedDOFs = 'xyXY' 

            r=b.shape.radius
            oldm=b.state.mass
            oldI=b.state.inertia

            m=oldm*3./4./r
            b.state.mass=m

            b.state.inertia[0] = 15./16./r*oldI[0]	#inertia with respect to x and y axes are not used and the computation here is wrong
            b.state.inertia[1] = 15./16./r*oldI[1]  #inertia with respect to x and y axes are not used and the computation here is wrong
            b.state.inertia[2] = 15./16./r*oldI[2]  #only inertia with respect to z axis is usefull
    o.dt=1e-6

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