yade-users team mailing list archive

[weijie <question690546@xxxxxxxxxxxxxxxxxxxxx>]

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

    Status: Answered => Open

weijie is still having a problem:
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()

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