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.
_______________________________________________
Mailing list: https://launchpad.net/~yade-users
Post to : yade-users@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
