yade-users team mailing list archive

[Xin Wu <question241615@xxxxxxxxxxxxxxxxxxxxx>]

New question #241615 on Yade:
https://answers.launchpad.net/yade/+question/241615

Dear Sir,
I want to write a Simulation and a new Engine with C++ Code. The Simulation see like this:
 
from yade import utils
from numpy import linspace
from numpy import arange
import gts
import itertools
from yade import pack
###########################################################################################
# The components of the batch are:
# 1. table with parameters, one set of parameters per line (ccc.table)
# 2. utils.readParamsFromTable which reads respective line from the parameter file
# 3. the simulation muse be run using yade-batch, not yade
#
# $ yade-batch --job-threads=1 03-oedometric-test.table 03-oedometric-test.py
#
rMean=.0096
rRelFuzz=.0016
maxLoad=4500
#minLoad=500
frictionAngleSt=radians(35)
frictionAngleBo=radians(23.5)
#tc = 0.001
#en = 0.3
#es = 0.3
a=0.05
spheremat=O.materials.append(FrictMat(density=2650,young=175e6,poisson=0.3,frictionAngle=frictionAngleBo))
steel=O.materials.append(FrictMat(young=210e9, poisson=.25,frictionAngle=frictionAngleSt,density=8000))
# load parameters from file if run in batch
# default values are used if not run from batch
#utils.readParamsFromTable(rMean=.0096,rRelFuzz=.016,maxLoad=1e6,minLoad=1e4)
# make rMean, rRelFuzz, maxLoad accessible directly as variables later
#from yade.params.table import *

# create box with free top, and ceate loose packing inside the box
from yade import pack, plot,qt
fIDSI=O.bodies.append(utils.geom.facetBox((.15,.15,.135),(.15,.15,.045),wallMask=15,material=steel))
fIDSII=O.bodies.append(utils.geom.facetBox((.15,.15,.045),(.15,.15,.045),wallMask=31,material=steel))
sp=pack.SpherePack()
sp.makeCloud((0,0,0),(0.3,0.3,0.3250),rMean=rMean,rRelFuzz=rRelFuzz)
sp.toSimulation(material=spheremat)

O.engines=[
   ForceResetter(),
   # sphere, facet, wall
   InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Wall_Aabb()]),
   InteractionLoop(
      # the loading plate is a wall, we need to handle sphere+sphere, sphere+facet, sphere+wall
      [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
      [Ip2_FrictMat_FrictMat_FrictPhys()],
      [Law2_ScGeom_FrictPhys_CundallStrack()]
   ),
   GravityEngine(gravity=(0,0,-9.81)),
   NewtonIntegrator(damping=0.5),
   # the label creates an automatic variable referring to this engine
   # we use it below to change its attributes from the functions called
   qt.SnapshotEngine(iterPeriod=60000,fileBase='/home/wuxin/Desktop/Ra-RW-235/78Fr07-',label='snapshooter'),
   PyRunner(command='checkUnbalanced()',realPeriod=1,label='checker'),
   TranslationEngine(translationAxis=[1,0,0],velocity=0,ids=fIDSI,label='Transl'),
]
O.dt=.25*utils.PWaveTimeStep()
#Transl.velocity=0.01
#qt.Controller()
#qt.View()

#### show how to use makeClumpTemplate():


#dyad:
relRadList1 = [.007,.0036]
relPosList1 = [[1,0,0],[0.001,0,0]]

#peanut:
relRadList2 = [.0036,0.007,.0036]
relPosList2 = [[0.0060,0,0],[0,0,0],[-0.0060,0,0]]

templates= []
templates.append(clumpTemplate(relRadii=relRadList1,relPositions=relPosList1))
templates.append(clumpTemplate(relRadii=relRadList2,relPositions=relPosList2))
templates= []
templates.append(clumpTemplate(relRadii=relRadList1,relPositions=relPosList1))
templates.append(clumpTemplate(relRadii=relRadList2,relPositions=relPosList2))

O.bodies.replaceByClumps(templates,[.2,.3])
# the following checkUnbalanced, unloadPlate and stopUnloading functions are all called by the 'checker'
# (the last engine) one after another; this sequence defines progression of different stages of the
# simulation, as each of the functions, when the condition is satisfied, updates 'checker' to call
# the next function when it is run from within the simulation next time

def AngVel():
    for b in O.bodies:
        if isinstance(b.shape,Sphere):
           #b.state.blockedDOFs='Y'
           b.state.angVel[1]=0
           #b.state.angVel[1]=0


# check whether the gravity deposition has already finished
# if so, add wall on the top of the packing and start the oedometric test

def checkUnbalanced():
    # at the very start, unbalanced force can be low as there is only few contacts, but it does not mean the packing is stable
    if O.iter<240000: return 
    # the rest will be run only if unbalanced is < .1 (stabilized packing)
    if utils.unbalancedForce()>.05: return 
    # add plate at the position on the top of the packing
    # the maximum finds the z-coordinate of the top of the topmost particle
    #Transl.velocity=0.01
    fIDSIII=O.bodies.append(utils.wall(max([b.state.pos[2]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)]),axis=2,sense=-1,material=spheremat))
    global plate        # without this line, the plate variable would only exist inside this function
    plate=O.bodies[fIDSIII]  # the last particles is the plate
    # Wall objects are "fixed" by default, i.e. not subject to forces
    # prescribing a velocity will therefore make it move at constant velocity (downwards)
    plate.state.vel=(0,0,-.025)
    # start plotting the data now, it was not interesting before
    O.engines=O.engines+[PyRunner(command='addPlotData()',iterPeriod=1000)]
    #O.engines=O.engines+[PyRunner(command='AngVel()',iterPeriod=1,)]
    # next time, do not call this function anymore, but the next one (unloadPlate) instead
    checker.command='unloadPlate()'
  

def unloadPlate():
   # if the force on plate exceeds maximum load, start unloading
   if abs(O.forces.f(plate.id)[2])>=maxLoad:
      O.engines=O.engines+[PyRunner(command='ServoContorl()',iterPeriod=1)]
      O.engines=O.engines+[PyRunner(command='AngVel()',iterPeriod=1)]
      # next time, do not call this function anymore, but the next one (stopUnloading) instead
      checker.command='stopUnloading()'

def stopUnloading():
    if abs(O.forces.f(plate.id)[2])>(maxLoad-1):
       if abs(O.forces.f(plate.id)[2])<(maxLoad+1):
          Transl.velocity=0.0010


def ServoContorl():
    Fzs=O.forces.f(plate.id)[2]
    if Fzs==0 :
       Fzs=0.000000000000001
    global G
    global KN
    KN=0.000000000000000001
    for i in plate.intrs() :
        KN=KN+i.phys.kn
    G=a/(KN*(O.dt))
    plate.state.vel[2]=(G*(Fzs-maxLoad))
 
def addPlotData():
    Fz=O.forces.f(plate.id)[2]
    F = 0
    global i
    for i in fIDSI:
        F += O.forces.f(i)[0]
        PX1=(O.bodies[i].state.pos[0]-O.bodies[i].state.refPos[0])
        SF=((0.30-PX1)*0.30)
    plot.addData(t=O.time,Fx=(-F/SF),PX=PX1,Fz=Fz,w=plate.state.pos[2]-plate.state.refPos[2],i=O.iter)
    if (O.bodies[i].state.pos[0]-O.bodies[i].state.refPos[0])>0.050:
       plot.saveDataTxt('175e6p03RW235.txt',vars=('i','PX','Fx','w','Fz'))
       O.pause()

# besides unbalanced force evolution, also plot the displacement-force diagram
plot.plots={'i':('w','Fz',),'PX':('Fx',)}
plot.plot()

############################################################################################################################
#O.saveTmp()
qt.Controller()
qt.View()
r=qt.Renderer()
#r.lightPos=Vector3(0,0,50)

##############################################################################################################################

O.run()
Can anyone give me some Examples to indicate, how can I write a new Engine .
For Example wiche Classe shall i use in the Head data.
And how can I use the Camke to Compile the Code (with whiche commands).

regards and Happy new Year
Xin

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