yade-users team mailing list archive

[Henry <question247021@xxxxxxxxxxxxxxxxxxxxx>]

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

Henry posted a new comment:
Hi Jan, 
    Many thanks for your help, and the script of my file, please check it. Thanks a lot!

# -*- coding: utf-8

from yade import utils,qt,export
from yade import plot
from yade.pack import *

mesh='Circle5cm' #name of gts mesh
surface=gts.read(open(mesh+'.gts'))

execfile('GenSPheresInSurface.py') #import a series of sphere in
"surface"

O.materials.append(FrictMat(young=5.0e10,poisson=.1,density=3000000,frictionAngle=.0))
    
for b in O.bodies:

    if isinstance(b.shape,Sphere): 
        b.state.blockedDOFs='zXY' # blocked movement along Z
        b.state.vel[1]=0.0
    if isinstance(b.shape,Clump):
        Clump_Num=Clump_Num+1
        b.state.blockedDOFs='yXYZ' # blocked movement along Y
        b.state.vel[1]=0.0

O.bodies.append(pack.gtsSurface2Facets(surface,color=(0.8,0.8,0.8),wire=True))

O.engines=[
	ForceResetter(),
	InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()],label='collider'),
	InteractionLoop(
		[Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
		[Ip2_FrictMat_FrictMat_FrictPhys()],
		[Law2_ScGeom_FrictPhys_CundallStrack()],
	),
	NewtonIntegrator(damping=0.9,gravity=(0,0.0,0.0)),
]

O.dt=1.0*utils.PWaveTimeStep()
O.timingEnabled=True;
O.run(2000)

Permit_Min_Force=500.0
Permit_Max_Force=2.0*Permit_Min_Force

def Cal_Interface():
    Rf=0.0
    Contact_Num=0
    for i in O.interactions:
            if not i.isReal:continue
            F0=i.phys.normalForce[0]
            F1=i.phys.normalForce[1]
            f=sqrt(F0*F0+F1*F1)#i.phys.normalForce().norm()
            Rf=Rf+f
            Contact_Num=Contact_Num+1
    return Rf, Contact_Num

def Gen_S():
    while 1:
        Sign=0
        O.run(1000)   
        (Rf, Contact_Num)=Cal_Interface()
        print("Rf=", Rf,"Contact_Num=", Contact_Num)    
        if Contact_Num>0:
            Rf=Rf/Contact_Num
        else:
            Rf=0.0
        if Rf>Permit_Min_Force:
            Mult=0.99
        if Rf<Permit_Max_Force:
            Mult=1.01
        if Rf>Permit_Min_Force and Rf<Permit_Max_Force:
            Sign=1
            Mult=1
        for o in O.bodies: 
            if isinstance(o.shape,Sphere):
                R=o.shape.radius
                o.shape.radius=R*Mult
        print("Contact_Num=", Contact_Num,"Sign=", Sign)
        if Sign==1: 
            break
            
            
Gen_S()

export.text(mesh+'.spheres')
print("Export the Spheres model successfuly!")

Henry

-- 
You received this question notification because you are a member of
yade-users, which is an answer contact for Yade.