WedgeYADE.py 15.9 KB
Newer Older
1
from __future__ import print_function
booncw's avatar
booncw committed
2 3
#CW BOON 2018
# Use the following algorithms:
4
# CW Boon, GT Houlsby, S Utili (2012).  A new algorithm for contact detection between convex polygonal and polyhedral particles in the discrete element method.  Computers and Geotechnics 44, 73-82. 
booncw's avatar
booncw committed
5 6 7
# CW Boon, GT Houlsby, S Utili (2015).  A new rock slicing method based on linear programming.  Computers and Geotechnics 65, 12-29. 

#The code runs some steps to stabilise.  After that a vertical cut is carried out.  And the simulations shows the failure mechanism of the slope.
8
#Display is saved to a vtk file in the "vtk" folder and the user is required to load it using ParaView.  Control the frequency of printing a vtk file using vtkRecorder.iterPeriod in python
booncw's avatar
booncw committed
9 10 11 12 13 14 15 16 17 18 19

#Disclaimer: This script is just for illustration to demonstrate the function of Block Gen and PotentialBlock Code, and not for other purpose outside for this intended use

#To use this script:
#Compile with
#ENABLE_POTENTIAL_BLOCKS=ON, and add  sudo apt-get install
#coinor-clp, 
#coinor-libclp-dev, 
#coinor-libclp1, 
#coinor-libosi1v5

20 21 22 23 24 25 26 27
import os
import errno
try:
   os.mkdir('./vtk/')
except OSError as exc:
   if exc.errno != errno.EEXIST:
      raise
   pass
booncw's avatar
booncw committed
28

29
name ='BlockGeneration'
booncw's avatar
booncw committed
30 31 32

p=BlockGen()
p.maxRatio = 10000.0;
33
p.minSize = 7.0;
booncw's avatar
booncw committed
34 35 36
p.density = 2700 # 23.6/9.81*1000.0 	#kg/m3
p.dampingMomentum = 0.8
p.viscousDamping = 0.0
37 38
p.Kn = 2.0e8  	#
p.Ks = 0.1e8 	#
booncw's avatar
booncw committed
39 40 41
p.frictionDeg = 50.0 #degrees
p.traceEnergy = False
p.defaultDt = 1e-4
42 43
p.rForPP = 0.2 #0.04
p.kForPP = 0.0
booncw's avatar
booncw committed
44 45
p.RForPP = 1800.0
p.gravity = [0,0,9.81]
46
#p.inertiaFactor = 1.0
booncw's avatar
booncw committed
47
p.initialOverlap = 1e-6 
48
p.exactRotation = True
49
#p.shrinkFactor = 1.0
booncw's avatar
booncw committed
50 51 52 53 54 55 56 57 58 59 60
p.boundarySizeXmin = 20.0; #South
p.boundarySizeXmax = 20.0; #North
p.boundarySizeYmin = 20.0; #West
p.boundarySizeYmax = 20; #East 4100
p.boundarySizeZmin = 0.0; #Up
p.boundarySizeZmax = 40.0; #Down
p.persistentPlanes = False
p.jointProbabilistic = True
p.opening = False
p.boundaries = True
p.slopeFace = False
61
p.calContactArea=True
booncw's avatar
booncw committed
62 63 64
p.twoDimension = False
p.unitWidth2D = 9.0
p.intactRockDegradation = True
65
#p.useFaceProperties = False
booncw's avatar
booncw committed
66 67 68 69 70 71 72 73 74 75 76
p.neverErase = False # Must be used when tension is on
p.sliceBoundaries = True
p.filenamePersistentPlanes = ''
p.filenameProbabilistic = './joints/jointC.csv'  #  PLEASE CHEK THIS  
p.filenameOpening = '' #'./joints/opening.csv'
p.filenameBoundaries = ''
p.filenameSlopeFace =''
p.filenameSliceBoundaries = ''
p.directionA=Vector3(1,0,0)
p.directionB=Vector3(0,1,0)
p.directionC=Vector3(0,0,1)
77 78 79 80

p.saveBlockGenData=True #if True, data of the BlockGen are written in a file; if False, they are displayed on the terminal
p.outputFile="BlockGen_Output.txt" #if empty, the block generation data are not written at all

booncw's avatar
booncw committed
81
p.load()
82
O.engines[1].avoidSelfInteractionMask=2
booncw's avatar
booncw committed
83 84
O.engines[2].lawDispatcher.functors[0].initialOverlapDistance = p.initialOverlap - 1e-6
O.engines[2].lawDispatcher.functors[0].allowBreakage = False
85 86
O.engines[2].lawDispatcher.functors[0].allowViscousAttraction=True
O.engines[2].lawDispatcher.functors[0].traceEnergy=False
booncw's avatar
booncw committed
87 88

O.engines[2].physDispatcher.functors[0].kn_i = 2e9
89
O.engines[2].physDispatcher.functors[0].ks_i = 0.1e9
booncw's avatar
booncw committed
90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106

from yade import plot

rockFriction = 40.0
boundaryFriction = 40.0
targetFriction = 40.0
waterHeight = 460.0  # PLEASE CHEK THIS 460.0 # 
startCountingBrokenBonds = False
minTimeStep = 1000000.0
westBodyId=10 #PLEASE CHEK THIS  
midBodyId =10 #PLEASE CHEK THIS  
eastBodyId =10 #PLEASE CHEK THIS  
originalPositionW = O.bodies[westBodyId].state.pos
originalPositionE = O.bodies[eastBodyId].state.pos
originalPositionM = O.bodies[midBodyId].state.pos
velocityDependency = False

107 108 109

# ----------------------------------------------------------------------------------------------------------------------------------------------- #
# Engines
booncw's avatar
booncw committed
110 111
O.engines[3].label='integration'
O.dt = 10.0e-5 #10e-4
112
O.engines = O.engines + [PotentialBlockVTKRecorder(fileName='./vtk/Wedge', iterPeriod=2000, twoDimension=False, sampleX=70, sampleY=70, sampleZ=70, maxDimension=100.0, REC_INTERACTION=True, REC_VELOCITY=True, label='vtkRecorder')]
113
O.engines = O.engines+ [PyRunner(iterPeriod=2000,command='calTimeStep()')] 
booncw's avatar
booncw committed
114 115 116 117
O.engines=O.engines+[PyRunner(iterPeriod=200,command='myAddPlotData()')]
O.engines = O.engines+ [PyRunner(iterPeriod=200,command='goToNextStage2()',label='dispChecker')] 


118 119 120
# ----------------------------------------------------------------------------------------------------------------------------------------------- #
# Material
O.materials.append(FrictMat(young=-1,poisson=-1,density=p.density,frictionAngle=radians(0.0), label='frictionless'))
booncw's avatar
booncw committed
121 122


123 124
# ----------------------------------------------------------------------------------------------------------------------------------------------- #
# Boundary Plates
booncw's avatar
booncw committed
125 126 127 128 129
thickness = 0.1*p.boundarySizeZmax
wire=False
highlight=False
kPP = p.kForPP
rPP = p.rForPP
130 131 132 133

# ----------------------------------------------------------------------------------------------------------------------------------------------- #
# Bottom plate
bbb = Body()
134
bbb.mask=3
135
color=[0,0.5,1]
booncw's avatar
booncw committed
136 137 138
aPP = [1,-1,0,0,0,0]
bPP = [0,0,1,-1,0,0]
cPP = [0,0,0,0,1,-1]
139 140 141 142 143 144
dPP = [p.boundarySizeXmax-rPP , p.boundarySizeXmax-rPP, p.boundarySizeYmax-rPP, p.boundarySizeYmax-rPP, thickness-rPP, thickness-rPP]
minmaxAabb = 1.05*Vector3( dPP[0], dPP[2], dPP[4] )
bbb.shape=PotentialBlock(k=kPP, r=rPP, R=0.0, a=aPP, b=bPP, c=cPP, d=dPP, id=len(O.bodies), isBoundary=True, color=color, wire=wire, highlight=highlight, AabbMinMax=True, minAabb=minmaxAabb, maxAabb=minmaxAabb, fixedNormal=True, boundaryNormal=Vector3(0,0,-1))
utils._commonBodySetup(bbb,bbb.shape.volume,bbb.shape.inertia,material='frictionless',pos= [0,0,p.boundarySizeZmax+thickness], fixed=True) 
bbb.state.ori = bbb.shape.orientation
O.bodies.append(bbb)
booncw's avatar
booncw committed
145

146 147 148 149

# ----------------------------------------------------------------------------------------------------------------------------------------------- #
# Lateral plate A
bA = Body()
150
bA.mask=3
151
color=[0,0.5,1]
booncw's avatar
booncw committed
152 153 154
aPP = [1,-1,0,0,0,0]
bPP = [0,0,1,-1,0,0]
cPP = [0,0,0,0,1,-1]
155 156 157 158 159 160 161 162 163 164 165
dPP = [p.boundarySizeXmax-rPP , p.boundarySizeXmax-rPP, thickness-rPP, thickness-rPP, 0.5*p.boundarySizeZmax-rPP, 0.5*p.boundarySizeZmax-rPP]
minmaxAabb = 1.05*Vector3( dPP[0], dPP[2], dPP[4] )
bA.shape=PotentialBlock(k=kPP, r=rPP, R=0.0, a=aPP, b=bPP, c=cPP, d=dPP, id=len(O.bodies), isBoundary=True, color=color, wire=wire, highlight=highlight, AabbMinMax=True, minAabb=minmaxAabb, maxAabb=minmaxAabb, fixedNormal=True, boundaryNormal=Vector3(0,-1,0))
utils._commonBodySetup(bA,bA.shape.volume,bA.shape.inertia,material='frictionless',pos= [0,p.boundarySizeYmax+thickness,0.5*p.boundarySizeZmax], fixed=True) 
bA.state.ori = bA.shape.orientation
O.bodies.append(bA)


# ----------------------------------------------------------------------------------------------------------------------------------------------- #
# Lateral plate B
bB = Body()
166
bB.mask=3
167
color=[0,0.5,1]
booncw's avatar
booncw committed
168 169 170
aPP = [1,-1,0,0,0,0]
bPP = [0,0,1,-1,0,0]
cPP = [0,0,0,0,1,-1]
171 172 173 174 175 176 177 178 179 180 181
dPP = [p.boundarySizeXmax-rPP, p.boundarySizeXmax-rPP, thickness-rPP, thickness-rPP, 0.5*p.boundarySizeZmax-rPP, 0.5*p.boundarySizeZmax-rPP]
minmaxAabb = 1.05*Vector3( dPP[0], dPP[2], dPP[4] )
bB.shape=PotentialBlock(k=kPP, r=rPP, R=0.0, a=aPP, b=bPP, c=cPP, d=dPP, id=len(O.bodies), isBoundary=True, color=color, wire=wire, highlight=highlight, AabbMinMax=True, minAabb=minmaxAabb, maxAabb=minmaxAabb, fixedNormal=True, boundaryNormal=Vector3(0,1,0))
utils._commonBodySetup(bB,bB.shape.volume,bB.shape.inertia,material='frictionless',pos= [0,-p.boundarySizeYmax-thickness,0.5*p.boundarySizeZmax], fixed=True) 
bB.state.ori = bB.shape.orientation
O.bodies.append(bB)


# ----------------------------------------------------------------------------------------------------------------------------------------------- #
# Lateral plate C
bC = Body()
182
bC.mask=3
183
color=[0,0.5,1]
booncw's avatar
booncw committed
184 185 186
aPP = [1,-1,0,0,0,0]
bPP = [0,0,1,-1,0,0]
cPP = [0,0,0,0,1,-1]
187 188 189 190 191 192 193 194 195 196 197
dPP = [thickness-rPP, thickness-rPP, p.boundarySizeYmax-rPP, p.boundarySizeYmax-rPP, 0.5*p.boundarySizeZmax-rPP, 0.5*p.boundarySizeZmax-rPP]
minmaxAabb = 1.05*Vector3( dPP[0], dPP[2], dPP[4] )
bC.shape=PotentialBlock(k=kPP, r=rPP, R=0.0, a=aPP, b=bPP, c=cPP, d=dPP, id=len(O.bodies), isBoundary=True, color=color, wire=wire, highlight=highlight, AabbMinMax=True, minAabb=minmaxAabb, maxAabb=minmaxAabb, fixedNormal=True, boundaryNormal=Vector3(1,0,0))
utils._commonBodySetup(bC,bC.shape.volume,bC.shape.inertia,material='frictionless',pos= [-p.boundarySizeXmax-thickness,0,0.5*p.boundarySizeZmax], fixed=True) 
bC.state.ori = bC.shape.orientation
O.bodies.append(bC)


# ----------------------------------------------------------------------------------------------------------------------------------------------- #
# Lateral plate D
bD = Body()
198
bD.mask=3
199
color=[0,0.5,1]
booncw's avatar
booncw committed
200 201 202
aPP = [1,-1,0,0,0,0]
bPP = [0,0,1,-1,0,0]
cPP = [0,0,0,0,1,-1]
203 204 205 206 207 208
dPP = [thickness-rPP, thickness-rPP, p.boundarySizeYmax-rPP, p.boundarySizeYmax-rPP, 0.5*p.boundarySizeZmax-rPP, 0.5*p.boundarySizeZmax-rPP]
minmaxAabb = 1.05*Vector3( dPP[0], dPP[2], dPP[4] )
bD.shape=PotentialBlock(k=kPP, r=rPP, R=0.0, a=aPP, b=bPP, c=cPP, d=dPP, id=len(O.bodies), isBoundary=True, color=color, wire=wire, highlight=highlight, AabbMinMax=True, minAabb=minmaxAabb, maxAabb=minmaxAabb, fixedNormal=True, boundaryNormal=Vector3(-1,0,0))
utils._commonBodySetup(bD,bD.shape.volume,bD.shape.inertia,material='frictionless',pos= [p.boundarySizeXmax+thickness,0,0.5*p.boundarySizeZmax], fixed=True) 
bD.state.ori = bD.shape.orientation
O.bodies.append(bD)
booncw's avatar
booncw committed
209 210


211
# ----------------------------------------------------------------------------------------------------------------------------------------------- #
booncw's avatar
booncw committed
212 213 214 215 216 217 218
def calTimeStep():
	global minTimeStep
	mkratio = 99999999.9
	maxK = 0.0
	minMass = 1.0e15
	for i in O.interactions:
		if i.isReal==True:
219 220
			dt1 = O.bodies[i.id1].state.mass/i.phys.kn
			dt2 = O.bodies[i.id2].state.mass/i.phys.kn
booncw's avatar
booncw committed
221 222 223 224 225 226 227 228 229 230
			if dt1 < dt2:
				presentDt = 0.15*sqrt(dt1)
				if minTimeStep > presentDt:
					minTimeStep = presentDt
					O.dt = minTimeStep
			else:
				presentDt = 0.15*sqrt(dt2)
				if minTimeStep > presentDt:
					minTimeStep = presentDt
					O.dt = minTimeStep
231 232


233
# ----------------------------------------------------------------------------------------------------------------------------------------------- #
booncw's avatar
booncw committed
234 235 236 237 238 239 240 241
def excavate():
	for b in O.bodies:
		if NorthWall(b.state.pos[0],b.state.pos[1],b.state.pos[2]) <0.001:
			if b.isClumpMember == True and b.shape.isBoundary==False:
				O.bodies.deleteClumpMember(O.bodies[b.clumpId], b)
			elif b.isClump == False and b.shape.isBoundary==False:
				O.bodies.erase(b.id)
				continue
242 243
	O.bodies.erase(bA.id)
	O.bodies.erase(bC.id)
booncw's avatar
booncw committed
244 245 246 247 248 249 250 251 252 253 254 255 256 257 258


prevDistance = O.bodies[westBodyId].state.pos[0]
tolDistance = 0.003 #0.1
tolDistance2 = 0.05
SRinProgress = False
SRcounter = 0
checkIter = 0
prevKE = 0.0
currentKE = 0.0
tolKE =10e10
initBondedContacts = 0
initDispRate = -1.0
prevDispRate = 0

259

260
# ----------------------------------------------------------------------------------------------------------------------------------------------- #
booncw's avatar
booncw committed
261 262 263 264
def changeKE(newKE):
	global tolKE
	tolKE = newKE

265

266
# ----------------------------------------------------------------------------------------------------------------------------------------------- #
booncw's avatar
booncw committed
267 268 269 270 271
def changeTolDist(newTol):
	global tolDistance
	tolDistance = newTol


272
# ----------------------------------------------------------------------------------------------------------------------------------------------- #
booncw's avatar
booncw committed
273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297
removeDampingBool = False
def goToNextStage2():
	global startCountingBrokenBonds
	global velocityDependency
	global waterHeight
	global boundaryFriction
	global rockFriction
	global targetFriction
	global prevDistance
	global originalPositionW
	global tolDistance
	global tolDistance2
	global checkIter
	global SRinProgress
	global SRcounter
	global prevKE
	global currentKE
	global tolKE
	global initDispRate
	global prevDispRate
	global removeDampingBool
	prevKE = currentKE
	KE = utils.kineticEnergy()
	currentKE = KE
	uf = utils.unbalancedForce()
298
	if O.iter>500 and removeDampingBool == False:
booncw's avatar
booncw committed
299
		removeDamping()
300
		removeDampingBool=True
301
	if O.iter>2000 and SRcounter == 0: # and uf<0.005:
302
		print(O.iter)
booncw's avatar
booncw committed
303 304 305 306 307 308 309 310
		O.pause()
		vtkRecorder.iterPeriod=1
		for b in O.bodies:
			b.state.vel = Vector3(0.0,0.0,0.0)
			b.state.angVel = Vector3(0.0,0.0,0.0)
		calTimeStep()
		excavate()
		dispChecker.iterPeriod=1
311 312 313 314
		SRcounter = 1
		O.step()
		vtkRecorder.iterPeriod=500
		O.run(20000)
booncw's avatar
booncw committed
315 316
		return

317

318
# ----------------------------------------------------------------------------------------------------------------------------------------------- #
booncw's avatar
booncw committed
319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336
def SouthWall(x,y,z):
	Xcentre1 = 0.0 	
	Ycentre1 = 0.0
	Zcentre1 = -0.0
	dip = 90.0
	dipdir = 315.0
	dipdirN = 0.0
	dipN = 90.0-dip
	if dipdir > 180.0:
		dipdirN = dipdir - 180.0
	else:
		dipdirN = dipdir + 180.0
	dipRad = dipN/180.0*pi
	dipdirRad = dipdirN/180.0*pi
	a = cos(dipdirRad)*cos(dipRad)
	b = sin(dipdirRad)*cos(dipRad)
	c = sin(dipRad)
	l = sqrt(a*a + b*b +c*c)
337 338 339
	a = a/l
	b = b/l
	c = c/l
booncw's avatar
booncw committed
340 341 342 343 344 345
	d = a*Xcentre1 + b*Ycentre1 + c*Zcentre1
	plane = a*x+ b*y +c*z - d
	if plane < 0.0:
		plane = 0.0
	return plane

346

347
# ----------------------------------------------------------------------------------------------------------------------------------------------- #
booncw's avatar
booncw committed
348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365
def NorthWall(x,y,z):
	Xcentre1 = 0.0 	
	Ycentre1 = 0.0
	Zcentre1 = -0.0
	dip = 90.0
	dipdir = 315.0
	dipdirN = 0.0
	dipN = 90.0-dip
	if dipdir > 180.0:
		dipdirN = dipdir - 180.0
	else:
		dipdirN = dipdir + 180.0
	dipRad = dipN/180.0*pi
	dipdirRad = dipdirN/180.0*pi
	a = cos(dipdirRad)*cos(dipRad)
	b = sin(dipdirRad)*cos(dipRad)
	c = sin(dipRad)
	l = sqrt(a*a + b*b +c*c)
366 367 368
	a = -a/l
	b = -b/l
	c = -c/l
booncw's avatar
booncw committed
369 370 371 372 373 374 375
	d = a*Xcentre1 + b*Ycentre1 + c*Zcentre1
	plane = a*x+ b*y +c*z - d
	if plane < 0.0:
		plane = 0.0
	return plane


376
# ----------------------------------------------------------------------------------------------------------------------------------------------- #
booncw's avatar
booncw committed
377 378 379 380 381 382 383 384
def removeDamping():
	for i in O.interactions:
		i.phys.viscousDamping = 0.5
		i.phys.cumulative_us = 0.0
	O.engines[2].physDispatcher.functors[0].viscousDamping = 0.5
	integration.damping= 0.0


385
# ----------------------------------------------------------------------------------------------------------------------------------------------- #
booncw's avatar
booncw committed
386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401
def myAddPlotData():
	global westBodyId
	global midBodyId
	global eastBodyId
	global originalPositionW
	global originalPositionM
	global originalPositionE
	global boundaryFriction
	KE = utils.kineticEnergy()
	uf = utils.unbalancedForce()
	displacementWx = O.bodies[westBodyId].state.pos[0] - originalPositionW[0]
	displacementW = (O.bodies[westBodyId].state.pos - originalPositionW).norm()
	displacementMx = O.bodies[midBodyId].state.pos[0] - originalPositionM[0]
	displacementM = (O.bodies[midBodyId].state.pos - originalPositionM).norm()
	displacementEx = O.bodies[eastBodyId].state.pos[0] - originalPositionE[0]
	displacementE = (O.bodies[eastBodyId].state.pos - originalPositionE).norm()
402
	plot.addData(timeStep=O.iter,timeStep1=O.iter,timeStep2=O.iter,timeStep3=O.iter,timeStep4=O.iter,timeStep5=O.iter,kineticEn=KE,unbalancedForce=uf, waterLevel=waterHeight,boundary_phi=boundaryFriction,displacement=displacementW,displacementWest=displacementW,dispWx=displacementWx,displacementMid=displacementM, dispMx=displacementMx,displacementEast=displacementE,dispEx=displacementEx)
booncw's avatar
booncw committed
403 404

plot.plots={'timeStep2':('kineticEn'),'timeStep3':(('displacementWest','ro-'),('dispWx','go-')),'timeStep1':(('displacementMid','ro-'),('dispMx','go-')),'timeStep5':(('displacementEast','ro-'),('dispEx','go-')),'timeStep4':('unbalancedForce')} #PLEASE CHECK
405 406
plot.plot() #Uncomment to view plots

booncw's avatar
booncw committed
407

408 409
# ----------------------------------------------------------------------------------------------------------------------------------------------- #
from yade import qt
410 411 412 413 414
try: 
	v=qt.View()
	vaxes=True
	v.viewDir=Vector3(0,-1,0)
	v.eyePosition=Vector3(0,200,0)
booncw's avatar
booncw committed
415

416 417 418 419
	v.eyePosition=Vector3(-77.42657409847706,84.2619166834641,-17.41745783023423)
	v.upVector=Vector3(0.1913254208509842,-0.25732151742252396,-0.9471959776136951)
	v.viewDir=Vector3(0.6412359985709529,-0.697845344639707,0.3191054200439123)
except: pass
booncw's avatar
booncw committed
420 421 422 423 424 425 426 427



O.step()
O.step()
#excavate()
O.step()
calTimeStep()
428
#O.run(20000)
booncw's avatar
booncw committed
429