Yade

How to quickly reduce unbalancedForce to approximate 0?

Asked by weijie

Dear all,

I am using the PotentialBlock function to do a regular dodecahedron stacking experiment, but I found that after many time steps unbalancedForce still fluctuates around 1, and the dodecahedron is still shaking. Is there any way to quickly reduce the unbalancedForce or let the polyhedron tend to be stationary and stable.

Thanks in advance.
Jie

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Vasileios Angelidakis
Solved:
Last query:
Last reply:

This question was reopened

Revision history for this message
Vasileios Angelidakis (vsangelidakis) said : 		#1

Hi Jie,

Depends on what stage of your analysis you are.

If you just want to stabilise your sample/particle during initialisation, you can use local damping (in the NewtonIntegrator), which is un-physical but will do the job.

If you want to keep some damping during the main stage of your main simulation, I would advise to not use local damping, but use a small amount of viscousDamping in the KnKsPBLaw instead.

If now you want to stabilise your system during initialisation momentarily, you can invoke "calm()" for some timesteps, to reset the velocity of the particles.

Hope this helps,
Vasileios

Revision history for this message
weijie (amandajoe) said : 		#2

Hi Vasileios, and thank you again.

What is the value range of viscousDamping?

Best regards,
Jie

Revision history for this message
Vasileios Angelidakis (vsangelidakis) said : 		#3

Hi Jie,

The viscousDamping parameter takes values from 0 to 1.0 (it is a percentage) and currently it only applies damping on the normalForce (the shearForce remains undamped for the time being).

Regards,
Vasileios

Revision history for this message
weijie (amandajoe) said : 		#4

Thanks Vasileios Angelidakis, that solved my question.

Revision history for this message
weijie (amandajoe) said : 		#5

Hi Vasileios, and thank you again.

If I use the calm () function to reset the velocity of the particles, the program appears as follows: shearForce: -nan -nan -nan, normalForce: -0.65154 1.59526 -0.0581212, viscousNormal: -0 0 -0, viscousShear: 0 0 0, geom normal: -0.377888 0.925237 -0.0337099, effective_phi: 1718.87, shearIncrement: 0 0 0, cs: 4.30304, incidentVs: 0 0 0, id1: 865, id2: 1562, debugShear: -0 0 0, phys-> mobilizedShear : -nan.

Some particles disappear, and the non-disappearing particles will then pass through the surface of the container. How to solve this problem?

Best regards,
Jie

Revision history for this message
Vasileios Angelidakis (vsangelidakis) said : 		#6

Hi Jie,

This is strange. Calm [1] only resets the translational and angular velocity (and the angular momentum), it does not reset forces.
The warning you get comes from [2], which tells you that you have no contact forces (hence the particles pass through boundaries) but I cannot understand why this happens. "Calm" alone cannot cause this behaviour. In order to be more helpful, I need a MWS to reproduce this.

All the best,
Vasileios

[1] https://gitlab.com/yade-dev/trunk/-/blob/master/pkg/dem/Shop_02.cpp#L921-934
[2] https://gitlab.com/yade-dev/trunk/-/blob/master/pkg/dem/KnKsPBLaw.cpp#L626

Revision history for this message
weijie (amandajoe) said : 		#7

Hi Vasileios, and thank you again.

I shortened the code, leaving the part where the main problem occurred. At 21,000 steps, I used the calm () function and the above problem occurred. If I don't use the calm () function, I find that after many steps, the entire system is unstable, unbalancedForce fluctuates around 2, and there is no tendency to decrease, which is obviously inconsistent with the actual reality. Can you help by the way to see what is causing this?

Best regards,
Jie

Following is my code：
#################

from yade import polyhedra_utils,pack,plot,utils,export,qt,ymport
import numpy as np
import math
import random
import os
import errno
try:
 os.mkdir('./vtk/')
except OSError as exc:
 if exc.errno != errno.EEXIST:
  raise
 pass

O.engines=[
 ForceResetter(),
 InsertionSortCollider([PotentialBlock2AABB()],verletDist=0.001, avoidSelfInteractionMask=2),
 InteractionLoop(
  [Ig2_PB_PB_ScGeom(twoDimension=False, unitWidth2D=1.0, calContactArea=True)],
  [Ip2_FrictMat_FrictMat_KnKsPBPhys(kn_i=1e8, ks_i=1e8, Knormal=1e8, Kshear=1e8, useFaceProperties=False, viscousDamping=0.2)],
  [Law2_SCG_KnKsPBPhys_KnKsPBLaw(label='law',neverErase=False, allowViscousAttraction=True, traceEnergy=False)]
 ),
 #GlobalStiffnessTimeStepper(),
 NewtonIntegrator(damping=0.0,exactAsphericalRot=True,gravity=[0,-9.81,0]),
 PotentialBlockVTKRecorder(fileName='./vtk/cubePBscaled',iterPeriod=300000,twoDimension=False,sampleX=50,sampleY=50,sampleZ=50,maxDimension=0.2,label='vtkRecorder')
]

powderDensity = 3000
O.materials.append(FrictMat(young=-1,poisson=-1,frictionAngle=30.0,density=powderDensity,label='frictionless'))
#######################

meanSize = 0.005
wallThickness = 0.5*meanSize
lengthOfBase = 0.250
heightOfBase = 0.600

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

sp=pack.SpherePack()
mn,mx=Vector3(-0.5*(lengthOfBase-wallThickness),0.5*meanSize,-0.5*(lengthOfBase-wallThickness)),Vector3(0.5*(lengthOfBase-wallThickness),0.5*heightOfBase,0.5*(lengthOfBase-wallThickness))#
R=sqrt(3.0)*0.009654
sp.makeCloud(mn,mx,R,0,-1,False)

r=0.01*meanSize
aa=[0.5257311121191336, -5.836787854364519e-17, 0.5257311121191336, -5.836787854364519e-17, 0.85065080835204, -0.5257311121191336, 0.85065080835204, -0.5257311121191336, -0.85065080835204, 5.836787854364519e-17, 5.836787854364519e-17, -0.8506508083520399]
bb=[0.85065080835204, -0.5257311121191336, -0.85065080835204, 0.5257311121191336, 5.836787854364519e-17, 0.85065080835204, 5.836787854364519e-17, -0.85065080835204, 5.836787854364519e-17, 0.5257311121191336, -0.5257311121191336, -0.0]
cc=[-5.836787854364519e-17, 0.85065080835204, 5.836787854364519e-17, 0.85065080835204, -0.5257311121191336, 5.836787854364519e-17, 0.5257311121191336, 5.836787854364519e-17, 0.5257311121191336, -0.85065080835204, -0.85065080835204, -0.5257311121191336]
dd=[0.01074978698202965, 0.010749786982029651, 0.01074978698202965, 0.01074978698202965, 0.01074978698202965, 0.01074978698202965, 0.010749786982029651, 0.010749786982029651, 0.01074978698202965, 0.01074978698202965, 0.010749786982029651, 0.01074978698202965]
for s in sp:
 b=Body()
 b.mask=1
 b.aspherical=True
 wire=False
 color=Vector3(random.random(),random.random(),random.random())
 highlight=False
 b.shape=PotentialBlock(k=0.0, r=r, R=0.0, a=aa, b=bb, c=cc, d=dd, isBoundary=False, color=color,
wire=wire, highlight=highlight, minAabb=Vector3(R,R,R), maxAabb=Vector3(R,R,R), AabbMinMax=True, fixedNormal=False, id=len(O.bodies))
 utils._commonBodySetup(b, b.shape.volume, b.shape.inertia, material='frictionless',pos=s[0], fixed=False)
 b.state.pos = s[0] #s[0] stores center
 b.state.ori = Quaternion((random.random(),random.random(),random.random()),random.random()) #s[2]
 O.bodies.append(b)
#######################################################

r=0.1*wallThickness
bbb=Body()
bbb.mask=3
wire=False
color=[0,0.5,1]
highlight=False
bbb.shape=PotentialBlock(k=0.0, r=r, R=0.0, a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1], d=[2.5*lengthOfBase-r,2.5*lengthOfBase-r,0.5*wallThickness-r,0.5*wallThickness-r,lengthOfBase/2.0-r,lengthOfBase/2.0-r], id=len(O.bodies), isBoundary=True, color=color, wire=wire, highlight=highlight, AabbMinMax=True, minAabb=1.05*Vector3(lengthOfBase*2.5,0.5*wallThickness,lengthOfBase/2.0), maxAabb=1.05*Vector3(lengthOfBase*2.5,0.5*wallThickness,lengthOfBase/2.0), fixedNormal=False)
utils._commonBodySetup(bbb, bbb.shape.volume, bbb.shape.inertia, material='frictionless', pos=[0,0,0], fixed=True)
bbb.state.pos = [2*lengthOfBase,0,0]
O.bodies.append(bbb)

#
bA=Body()
bA.mask=3
wire=False
color=[0,0.5,1]
highlight=False
bA.shape=PotentialBlock(k=0.0, r=r, R=0.0, a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1], d=[0.5*wallThickness-r,0.5*wallThickness-r,0.5*heightOfBase-r,0.5*heightOfBase-r,0.5*lengthOfBase-r,0.5*lengthOfBase-r], id=len(O.bodies), isBoundary=True, color=color, wire=wire, highlight=highlight, AabbMinMax=True, minAabb=1.05*Vector3(0.4*wallThickness,0.5*heightOfBase,0.5*lengthOfBase), maxAabb=1.05*Vector3(0.4*wallThickness,0.5*heightOfBase,0.5*lengthOfBase), fixedNormal=False)
utils._commonBodySetup(bA, bA.shape.volume, bA.shape.inertia, material='frictionless', pos=[0,0,0], fixed=True)
bA.state.pos = [0.5*lengthOfBase,0.5*heightOfBase,0]
O.bodies.append(bA)

bB=Body()
bB.mask=3
wire=False
color=[0,0.5,1]
highlight=False
bB.shape=PotentialBlock(k=0.0, r=r, R=0.0, a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1], d=[0.5*wallThickness-r,0.5*wallThickness-r,0.5*heightOfBase-r,0.5*heightOfBase-r,0.5*lengthOfBase-r,0.5*lengthOfBase-r], id=len(O.bodies), isBoundary=True, color=color, wire=wire, highlight=highlight, AabbMinMax=True, minAabb=1.05*Vector3(0.4*wallThickness,0.5*heightOfBase,0.5*lengthOfBase), maxAabb=1.05*Vector3(0.4*wallThickness,0.5*heightOfBase,0.5*lengthOfBase), fixedNormal=False)
utils._commonBodySetup(bB, bB.shape.volume, bB.shape.inertia, material='frictionless', pos=[0,0,0], fixed=True)
bB.state.pos = [-0.5*lengthOfBase,0.5*heightOfBase,0]
O.bodies.append(bB)

bC=Body()
bC.mask=3
wire=True
color=[0,0.5,1]
highlight=False
bC.shape=PotentialBlock(k=0.0, r=r, R=0.0, a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1], d=[2.5*lengthOfBase-r,2.5*lengthOfBase-r,0.5*heightOfBase-r,0.5*heightOfBase-r,0.5*wallThickness-r,0.5*wallThickness-r], id=len(O.bodies), isBoundary=True, color=color, wire=wire, highlight=highlight, AabbMinMax=True, minAabb=1.05*Vector3(2.5*lengthOfBase,0.5*heightOfBase,0.4*wallThickness), maxAabb=1.05*Vector3(2.5*lengthOfBase,0.5*heightOfBase,0.4*wallThickness), fixedNormal=False)
utils._commonBodySetup(bC, bC.shape.volume, bC.shape.inertia, material='frictionless', pos=[0,0,0], fixed=True)
bC.state.pos = [2*lengthOfBase,0.5*heightOfBase,0.5*lengthOfBase]
O.bodies.append(bC)

bD=Body()
bD.mask=3
wire=False
color=[0,0.5,1]
highlight=False
bD.shape=PotentialBlock(k=0.0, r=r, R=0.0, a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1], d=[2.5*lengthOfBase-r,2.5*lengthOfBase-r,0.5*heightOfBase-r,0.5*heightOfBase-r,0.5*wallThickness-r,0.5*wallThickness-r], id=len(O.bodies), isBoundary=True, color=color, wire=wire, highlight=highlight, AabbMinMax=True, minAabb=1.05*Vector3(2.5*lengthOfBase,0.5*heightOfBase,0.4*wallThickness), maxAabb=1.05*Vector3(2.5*lengthOfBase,0.5*heightOfBase,0.4*wallThickness), fixedNormal=False)
utils._commonBodySetup(bD, bD.shape.volume, bD.shape.inertia, material='frictionless', pos=[0,0,0], fixed=True)
bD.state.pos = [2*lengthOfBase,0.5*heightOfBase,-0.5*lengthOfBase]
O.bodies.append(bD)
##########################################

def myAddPlotData():
 global wallThickness
 global meanSize
 uf=utils.unbalancedForce()
 if isnan(uf):
  uf = 1.0
 KE = utils.kineticEnergy()

 for b in O.bodies:
  if b.state.pos[1] < -5.0*meanSize and len(b.state.blockedDOFs)==0: #i.e. fixed==False
   O.bodies.erase(b.id)

 plot.addData(timeStep1=O.iter,timeStep2=O.iter,timeStep3=O.iter,timeStep4=O.iter,time=O.time,unbalancedForce=uf,kineticEn=KE)

from yade import plot
plot.plots={'timeStep1':('unbalancedForce'),'timeStep2':('kineticEn'),}
plot.plot() #Uncomment to view plots
O.engines=O.engines+[PyRunner(iterPeriod=20,command='myAddPlotData()')]
O.engines=O.engines+[PyRunner(iterPeriod=20,command='still()')]
O.dt = 0.2*sqrt(0.3*O.bodies[0].state.mass/3.0e5)

def still():
 if O.iter>21000 and O.iter<22000:
  for b in O.bodies:
   calm()

from yade import qt
v=qt.View()
v.sceneRadius=3.0

O.saveTmp()

Revision history for this message
weijie (amandajoe) said : 		#8

You can change the step size of the calm () function in my code to make it appear earlier, and the same problem will occur.

Revision history for this message
Best Vasileios Angelidakis (vsangelidakis) said : 		#9

Hi Jie,

I had a look in your script.

There are several reasons causing instability:
- Your particles pass through the boundaries because the particles forming the boundaries are very thin. A value of wallThickness=0.05 worked well for me.
- frictionAngle must be given in radians, not degrees
- Some particles were erased because you had a statement to erase them inside the myAddPlotData() function

Some further advice:
- Use larger values for "r". If your dodecahedra have an inradius d=0.107, then I would advice using r=0.053 (around half the min(d) value). Please recheck the size of your dodecahedra, because I used d=array(dd)-r in the shape class and I am not sure if you had subtracted "r" from the dd values already.
- calm() resets the velocities of all bodies; no need to put it in a loop

Last, I replaced for you the PotentialBlockVTKRecorder with export.exportPotentialBlocks, so you don't have to worry about using minAabb, maxAabb. Also, the "id" parameter is no longer used in the PotentialBlock shape class.

You can compare the updated script below with your previous script using meld [1]

All the best,
Vasileios

[1] https://meldmerge.org/

from yade import polyhedra_utils,pack,plot,utils,export,qt
import numpy as np
import math
import random
import os
import errno
try:
 os.mkdir('./vtk/')
except OSError as exc:
 if exc.errno != errno.EEXIST:
  raise
  pass

#-------------------------------------------
# Engines
Kn=1e8
Ks=Kn/10

O.engines=[
 ForceResetter(),
 InsertionSortCollider([PotentialBlock2AABB()],verletDist=0.005, avoidSelfInteractionMask=2),
 InteractionLoop(
  [Ig2_PB_PB_ScGeom(calContactArea=True)],
  [Ip2_FrictMat_FrictMat_KnKsPBPhys(kn_i=Kn, ks_i=Ks, Knormal=Kn, Kshear=Ks, viscousDamping=0.2)],
  [Law2_SCG_KnKsPBPhys_KnKsPBLaw(allowViscousAttraction=True, traceEnergy=False)]
 ),
 NewtonIntegrator(damping=0.0,exactAsphericalRot=True,gravity=[0,-9.81,0]),
]

O.materials.append(FrictMat(young=-1,poisson=-1,frictionAngle=radians(30.0),density=3000,label='frictionless'))

#-------------------------------------------
# Dimensions
meanSize = 0.05
wallThickness = 0.5*meanSize
lengthOfBase = 0.250
heightOfBase = 0.600

#-------------------------------------------
# Make cloud
sp=pack.SpherePack()
mn,mx=Vector3(-0.5*(lengthOfBase-wallThickness),0.5*meanSize,-0.5*(lengthOfBase-wallThickness)),Vector3(0.5*(lengthOfBase-wallThickness),0.5*heightOfBase,0.5*(lengthOfBase-wallThickness))#
R=sqrt(3.0)*0.009654 #You can find exact formulae for the circumradius of dodecahedra
sp.makeCloud(mn,mx,R,0,-1,False)

#r=0.01*meanSize
aa=[0.5257311121191336, -5.836787854364519e-17, 0.5257311121191336, -5.836787854364519e-17, 0.85065080835204, -0.5257311121191336, 0.85065080835204, -0.5257311121191336, -0.85065080835204, 5.836787854364519e-17, 5.836787854364519e-17, -0.8506508083520399]
bb=[0.85065080835204, -0.5257311121191336, -0.85065080835204, 0.5257311121191336, 5.836787854364519e-17, 0.85065080835204, 5.836787854364519e-17, -0.85065080835204, 5.836787854364519e-17, 0.5257311121191336, -0.5257311121191336, -0.0]
cc=[-5.836787854364519e-17, 0.85065080835204, 5.836787854364519e-17, 0.85065080835204, -0.5257311121191336, 5.836787854364519e-17, 0.5257311121191336, 5.836787854364519e-17, 0.5257311121191336, -0.85065080835204, -0.85065080835204, -0.5257311121191336]
dd=[0.01074978698202965, 0.010749786982029651, 0.01074978698202965, 0.01074978698202965, 0.01074978698202965, 0.01074978698202965, 0.010749786982029651, 0.010749786982029651, 0.01074978698202965, 0.01074978698202965, 0.010749786982029651, 0.01074978698202965]

r=min(dd)/2 #Suggested value

for s in sp:
 b=Body()
 b.mask=1
 b.aspherical=True
 color=Vector3(random.random(),random.random(),random.random())
 b.shape=PotentialBlock(k=0.0, r=r, R=0.0, a=aa, b=bb, c=cc, d=np.array(dd)-r, isBoundary=False, color=color, wire=False, AabbMinMax=True, fixedNormal=False)
 utils._commonBodySetup(b, b.shape.volume, b.shape.inertia, material='frictionless',pos=s[0], fixed=False)
 b.state.pos = s[0] #s[0] stores center
 b.state.ori = Quaternion((random.random(),random.random(),random.random()),random.random()) #s[2]
 O.bodies.append(b)

#-------------------------------------------
# Bottom plate
color=[0,0.5,1]
r=0.15*wallThickness

bbb=Body()
bbb.mask=3
wire=False
bbb.shape=PotentialBlock(k=0.0, r=r, R=0.0, a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1], d=[2.5*lengthOfBase-r,2.5*lengthOfBase-r,0.5*wallThickness-r,0.5*wallThickness-r,lengthOfBase/2.0-r,lengthOfBase/2.0-r], isBoundary=False, color=color, wire=wire, AabbMinMax=True, fixedNormal=False)
utils._commonBodySetup(bbb, bbb.shape.volume, bbb.shape.inertia, material='frictionless', pos=[0,0,0], fixed=True)
bbb.state.pos = [2*lengthOfBase,0,0]
O.bodies.append(bbb)

bA=Body()
bA.mask=3
bA.shape=PotentialBlock(k=0.0, r=r, R=0.0, a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1], d=[0.5*wallThickness-r,0.5*wallThickness-r,0.5*heightOfBase-r,0.5*heightOfBase-r,0.5*lengthOfBase-r,0.5*lengthOfBase-r], isBoundary=False, color=color, wire=False, AabbMinMax=True, fixedNormal=False)
utils._commonBodySetup(bA, bA.shape.volume, bA.shape.inertia, material='frictionless', pos=[0,0,0], fixed=True)
bA.state.pos = [0.5*lengthOfBase,0.5*heightOfBase,0]
O.bodies.append(bA)

bB=Body()
bB.mask=3
bB.shape=PotentialBlock(k=0.0, r=r, R=0.0, a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1], d=[0.5*wallThickness-r,0.5*wallThickness-r,0.5*heightOfBase-r,0.5*heightOfBase-r,0.5*lengthOfBase-r,0.5*lengthOfBase-r], isBoundary=False, color=color, wire=False, AabbMinMax=True, fixedNormal=False)
utils._commonBodySetup(bB, bB.shape.volume, bB.shape.inertia, material='frictionless', pos=[0,0,0], fixed=True)
bB.state.pos = [-0.5*lengthOfBase,0.5*heightOfBase,0]
O.bodies.append(bB)

bC=Body()
bC.mask=3
bC.shape=PotentialBlock(k=0.0, r=r, R=0.0, a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1], d=[2.5*lengthOfBase-r,2.5*lengthOfBase-r,0.5*heightOfBase-r,0.5*heightOfBase-r,0.5*wallThickness-r,0.5*wallThickness-r], isBoundary=False, color=color, wire=True, AabbMinMax=True, fixedNormal=False)
utils._commonBodySetup(bC, bC.shape.volume, bC.shape.inertia, material='frictionless', pos=[0,0,0], fixed=True)
bC.state.pos = [2*lengthOfBase,0.5*heightOfBase,0.5*lengthOfBase]
O.bodies.append(bC)

bD=Body()
bD.mask=3
bD.shape=PotentialBlock(k=0.0, r=r, R=0.0, a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1], d=[2.5*lengthOfBase-r,2.5*lengthOfBase-r,0.5*heightOfBase-r,0.5*heightOfBase-r,0.5*wallThickness-r,0.5*wallThickness-r], isBoundary=False, color=color, wire=False, AabbMinMax=True, fixedNormal=False)
utils._commonBodySetup(bD, bD.shape.volume, bD.shape.inertia, material='frictionless', pos=[0,0,0], fixed=True)
bD.state.pos = [2*lengthOfBase,0.5*heightOfBase,-0.5*lengthOfBase]
O.bodies.append(bD)

#-------------------------------------------
# Plot results
def myAddPlotData():
 global wallThickness
 global meanSize
 uf=utils.unbalancedForce()
 if isnan(uf):
  uf = 1.0
 KE = utils.kineticEnergy()

# for b in O.bodies:
# if b.state.pos[1] < -5.0*meanSize and len(b.state.blockedDOFs)==0: #i.e. fixed==False
# O.bodies.erase(b.id)

 plot.addData(timeStep1=O.iter,timeStep2=O.iter,timeStep3=O.iter,timeStep4=O.iter,time=O.time,unbalancedForce=uf,kineticEn=KE)

from yade import plot
plot.plots={'timeStep1':('unbalancedForce'),'timeStep2':('kineticEn'),}
plot.plot() #Uncomment to view plots
O.engines=O.engines+[PyRunner(iterPeriod=20,command='myAddPlotData()')]
O.engines=O.engines+[PyRunner(iterPeriod=200,command='still()',firstIterRun=4000)] #enable this engine after 4000 iterations
O.dt = 0.2*sqrt(0.3*O.bodies[0].state.mass/Kn)*100

#-------------------------------------------
# Invoke calm()
def still():
 KE = utils.kineticEnergy()
 if KE<0.01:
  calm()

#-------------------------------------------
# exportPotentialBlocks
from yade import export
vtkExporter_PotentialBlocks = export.VTKExporter('vtk/pb')

def vtkExport():
 vtkExporter_PotentialBlocks.exportPotentialBlocks(what=dict(n='b.id'))

O.engines=O.engines+[PyRunner(iterPeriod=500,command='vtkExport()')]

from yade import qt
v=qt.View()
v.sceneRadius=3.0

O.saveTmp()

Revision history for this message
Vasileios Angelidakis (vsangelidakis) said : 		#10

A correction: I used meanSize = 0.05. wallThickness is half that.

Cheers,
Vasileios

Revision history for this message
weijie (amandajoe) said : 		#11

Thanks Vasileios Angelidakis, that solved my question.

Revision history for this message
Bruno Chareyre (bruno-chareyre) said : 		#12

A short comment for future reference: in my experience tricks like "calm()" are needed only when there is another problem in parameters or boundary conditions (e.g. too dynamic).
When the parameters are not doomed calm() only makes the whole thing slower (which is logical: if a system is unbalanced, freezing it will not make it converge faster to equilibrium). It is thus more efficient to track the initial problem when possible.
Bruno