yade-users team mailing list archive

Thread
Date
[Question #707390]: segmentation fault (PFacet)

To: yade-users@xxxxxxxxxxxxxxxxxxx
From: Ruidong LI <question707390@xxxxxxxxxxxxxxxxxxxxx>
Date: Tue, 25 Jul 2023 21:15:29 -0000
Reply-to: question707390@xxxxxxxxxxxxxxxxxxxxx
Sender: noreply@xxxxxxxxxxxxx
New question #707390 on Yade:
https://answers.launchpad.net/yade/+question/707390

Hi! I am trying to replicate the work using PFacet as the flexible membrane for the cylindrical triaxial test. I don't why this code always crashes (segmentation fault) at the stabilization stage. Can someone help me with this issue? Many thanks.

The MWE is presented below:
# import essential modules 
from __future__ import print_function
from yade import pack, ymport, qt, plot, geom
from yade.gridpfacet import *
import gts, os.path, locale, random, math
locale.setlocale(locale.LC_ALL, 'en_US.UTF-8')

###############################################################
###            1. DEFINE VARIABLES AND MATERIALS            ###
###############################################################

# 1.a). define variables
key = 'Triaxial_Undrained'	# file name to be saved 
young=550e6	# normal contact stiffness
compFricDegree = 1.8 # initial contact friction during the confining phase
finalFricDegree = 43	 # contact friction during the deviatoric loading 
poisson = 0.3	# shear-to-normal stiffness ratio
isoStress = 110000	# confining stress
conStress = 100000	 # confining stress for deviatoric loading stage 
width = 1.4e-1	# sample width
height = 2.8e-1	# target sample height(after consolidation) 
height_0 = 3.2e-1	 # initial sample height 
num_spheres=1000	 # number of spheres
R_p = 0.0084	# mean particle radius
rCoff = 10	# thickness of top and bot sphere cap (based on rParticle) 
rParticle = 0.02e-1	# membrane grid seed size
alpha = 8
rate = 0.1	# loading rate (strain rate)
damp = 0.3	# damping coefficient 
targetPorosity = 0.43		# target porosity 
thresholdvalue = 0.05		 # threshold unbalance force
final_rate = 0.1	# strain rate for deviator loading 
thresholdstrain = 0.06	 # threshold axial strain for terminate 
enlargefactor = 1.00

# A.b). create materials for sand spheres and plates 
Sand = O.materials.append(CohFrictMat(
young=young,poisson=poisson,frictionAngle=radians(compFricDegree), alphaKr=0.25,alphaKtw=0,etaRoll=0.005,etaTwist=0, normalCohesion=5e6,shearCohesion=5e6, momentRotationLaw=True,density=2650,label='spheres'
))

# A.c). create membrane materials
 
GridMat = O.materials.append(CohFrictMat( young=100e6,poisson=0.3,density=2650,frictionAngle=radians(0), alphaKr=0,alphaKtw=0,etaRoll=0,etaTwist=0, normalCohesion=1e9,shearCohesion=1e9, momentRotationLaw=True,label='gridNodeMat'
))
pFacetMat = O.materials.append(FrictMat( young=100e6,poisson=0.3,density=2650,frictionAngle=radians(0),label='pFacetMat'
))

# A.d). create TOP & BOT plate materials 
frictMat = O.materials.append(FrictMat(
young=100e6,poisson=0.3,density=2650,frictionAngle=radians(0),label='frictMat'
))

###############################################################
###                 3. DEFINE GLOBAL ENGINES                ###
###############################################################

#**********************************************************************# 
O.engines=[
	ForceResetter(), 
	InsertionSortCollider([
		Bo1_Sphere_Aabb(), 
		Bo1_Facet_Aabb(), 
		Bo1_PFacet_Aabb(), 
		Bo1_GridConnection_Aabb()
	]),
	InteractionLoop(
	[
		Ig2_Sphere_Sphere_ScGeom6D(), 
		Ig2_GridNode_GridNode_GridNodeGeom6D(), 
		Ig2_GridConnection_GridConnection_GridCoGridCoGeom(), 
		Ig2_Sphere_PFacet_ScGridCoGeom(), 
		Ig2_Facet_Sphere_ScGeom6D()
	], 
	[
		Ip2_FrictMat_FrictMat_FrictPhys(), 
		Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(label="cohesiveIp")
	],
	[
		Law2_ScGeom_FrictPhys_CundallStrack(), 
		Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm=True,always_use_moment_law=False,label='cohesiveLaw'),
		Law2_GridCoGridCoGeom_FrictPhys_CundallStrack(), 
		Law2_ScGridCoGeom_FrictPhys_CundallStrack(),
	],
	label="iloop"
	),
	GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=25,timestepSafetyCoefficient=0.8),
	NewtonIntegrator(gravity=(0,0,0),damping=0.1,label='newton')
]
#**********************************************************************# 
O.engines[2].lawDispatcher.functors[1].always_use_moment_law=False 
O.engines[2].physDispatcher.functors[1].setCohesionOnNewContacts=False

###############################################################
###                2. IMPORT CT-BASED PACKING               ###
###############################################################
# C.a). generate random dense sphere pack
pred = pack.inCylinder((0,0,0),(0,0,height_0),.5*width) 
sp = pack.randomDensePack(pred,spheresInCell=num_spheres,radius=R_p,rRelFuzz=0.3, returnSpherePack=True,memoDbg=True,memoizeDb='/tmp/loosePackings11.sqlite')
sand=sp.toSimulation(color=(0,1,1),material=Sand)

# C.b). define different sections of sphere pack
bot = [O.bodies[s] for s in sand if O.bodies[s].state.pos[2]<rParticle*rCoff]
top = [O.bodies[s] for s in sand if O.bodies[s].state.pos[2]>height_0-rParticle*rCoff] 
tot = [O.bodies[s] for s in sand if O.bodies[s].state.pos[2]<=height_0]

# C.c). detect the position of particles in top & bot layer 
top_limit = 0
top_id = 0 
for s in top:
	if s.state.pos[2]>=top_limit: 
		top_limit = s.state.pos[2] 
		top_id = s.id
bot_limit = height_0
bot_id = 0
for s in bot:
	if s.state.pos[2]<=bot_limit: 
		bot_limit = s.state.pos[2] 
		bot_id = s.id

# C.d). create facet wall around particle packing 
facets = []
nw = 45
nh = 15
rCyl2 = 0.5*width / cos(pi/float(nw)) 
for r in range(nw):
	for h in range(nh):
		v1 = Vector3( rCyl2*cos(2*pi*(r+0)/float(nw)), rCyl2*sin(2*pi*(r+0)/float(nw)), height_0*(h+0)/float(nh) )
		v2 = Vector3( rCyl2*cos(2*pi*(r+1)/float(nw)), rCyl2*sin(2*pi*(r+1)/float(nw)), height_0*(h+0)/float(nh) )
		v3 = Vector3( rCyl2*cos(2*pi*(r+1)/float(nw)), rCyl2*sin(2*pi*(r+1)/float(nw)), height_0*(h+1)/float(nh) )
		v4 = Vector3( rCyl2*cos(2*pi*(r+0)/float(nw)), rCyl2*sin(2*pi*(r+0)/float(nw)), height_0*(h+1)/float(nh) )
		f1 = facet((v1,v2,v3),color=(0,0,1),material=frictMat) 
		f2 = facet((v1,v3,v4),color=(0,0,1),material=frictMat) 
		facets.extend((f1,f2))
wall = O.bodies.append(facets)

# C.e). define different sections of facet wall 
for b in wall:
	O.bodies[b].state.blockedDOFs = 'xyzXYZ' 
	O.bodies[b].state.vel = (0,0,0)

# C.f). create top & bot facet plate 
facets3 = []
nw=45
rCyl2 = (.6*width+2*rParticle) / cos(pi/float(nw)) 
for r in range(nw):
	if r%2==0:
		v1 = Vector3( rCyl2*cos(2*pi*(r+0)/float(nw)), rCyl2*sin(2*pi*(r+0)/float(nw)), height_0 )
		v2 = Vector3( rCyl2*cos(2*pi*(r+1)/float(nw)), rCyl2*sin(2*pi*(r+1)/float(nw)), height_0 )
		v3 = Vector3( rCyl2*cos(2*pi*(r+2)/float(nw)), rCyl2*sin(2*pi*(r+2)/float(nw)), height_0 )
		v4 = Vector3( 0, 0, height_0 )
		f1 = facet((v1,v2,v4),color=(0,0,0),material=frictMat)
		f2 = facet((v2,v3,v4),color=(0,0,0),material=frictMat) 
		facets3.extend((f1,f2))
topcap = O.bodies.append(facets3) 
facets3 = []
for r in range(nw):
	if r%2==0:
		v1 = Vector3( rCyl2*cos(2*pi*(r+0)/float(nw)), rCyl2*sin(2*pi*(r+0)/float(nw)), 0 ) 
		v2 = Vector3( rCyl2*cos(2*pi*(r+1)/float(nw)), rCyl2*sin(2*pi*(r+1)/float(nw)), 0 )
		v3 = Vector3( rCyl2*cos(2*pi*(r+2)/float(nw)), rCyl2*sin(2*pi*(r+2)/float(nw)), 0 )
		v4 = Vector3( 0, 0, 0 )
		f1 = facet((v1,v2,v4),color=(0,0,0),material=frictMat) 
		f2 = facet((v2,v3,v4),color=(0,0,0),material=frictMat) 
		facets3.extend((f1,f2))
botcap = O.bodies.append(facets3)

# C.g). define top & bot wall id 
for s in topcap:
	top_id = s
bot_id = 0
for s in botcap:
	bot_id = s

# D.h). calculate porosity 
V_sand = 0
num_sand = 0 
for b in sand:
	r = O.bodies[b].shape.radius
	V_sand += 1.3333333*3.1416*r*r*r 
	num_sand +=1
porosity = 1-V_sand/(.25*width*width*3.1416*height_0)
print('v_sand= ',V_sand,' number of sand: ',num_sand,'porosity is: ',porosity)
O.pause()
###############################################################
###                   4. DEFINE FUNCTIONS                   ###
###############################################################

####################################################################### #	D. DEFINING ADD-ON FUNCTIONS	# #######################################################################
# D.a). a function for saving variables 
def plotAddData():
	f1 = sum(O.forces.f(b)[2] for b in topcap) 
	f2 = sum(O.forces.f(b)[2] for b in botcap) 
	f11 = sum(O.forces.f(b.id)[2] for b in top) 
	f22 = sum(O.forces.f(b.id)[2] for b in bot) 
	fa = abs(.5*(f2-f1))
	fa1 = abs(.5*(f22-f11))
	e = (top[0].state.displ()[2] - bot[0].state.displ()[2])/height_0-e_ini 
	f4 = 0
	r_cum = 0
	count = 0
	Area = 0

	for f in membrane_grid:
		f.shape.color = Vector3(0,0,0) 
		x = f.state.pos[0]
		y = f.state.pos[1]
		z = f.state.pos[2]
		dist = math.sqrt(x*x+y*y) 
		n = Vector3(x/dist,y/dist,0) 
		a = 2*3.1416*dist/(360/alpha)*6*rParticle*(1+e)
		if z<=O.bodies[top_id].state.pos[2] and z>=O.bodies[bot_id].state.pos[2]: 
			count += 1
			r_local = dist 
			r_cum += r_local
	r_avg = r_cum/count-rParticle
	Area = r_avg*2*3.1416*(O.bodies[top_id].state.pos[2] - O.bodies[bot_id].state.pos[2])
	area_avg = Area/count
	s = fa/(3.1416*r_avg*r_avg) 
	s1 = fa1/(3.1416*r_avg*r_avg) 
	for b in membrane_grid:
		x = b.state.pos[0] 
		y = b.state.pos[1]
		z = b.state.pos[2]
		dist = math.sqrt(x*x+y*y) 
		n = Vector3(x/dist,y/dist,0)
		a = 2*3.1416*dist/(360/alpha)*6*rParticle*(1+e) 	
		if z<=O.bodies[top_id].state.pos[2] and z>=O.bodies[bot_id].state.pos[2]: 
			f_local = O.forces.f(b.id)
			length = math.sqrt(f_local[0]*f_local[0]+f_local[1]*f_local[1]+f_local[2]*f_local[2])
			cos_theta = (n[0]*f_local[0]+n[1]*f_local[1]+n[2]*f_local[2])/length
			p_normal = (length*cos_theta/a) 
			f4 += (p_normal)
	p = abs(f4/count/1000)
	h = O.bodies[top_id].state.pos[2] - O.bodies[bot_id].state.pos[2]
	VV = h*r_avg*r_avg*3.1416 
	dV = VV-V_ini
	ev = -((O.bodies[top_id].state.pos[2] - O.bodies[bot_id].state.pos[2])*r_avg*r_avg*3.1416-V_ini)/V_ini
	er = (r_avg-R_avg)/R_avg
	if (abs(e*100)>thresholdstrain*100): 
		O.pause()
	plot.addData(
		i = O.iter,
		q = (abs(s)-p*1000)/1000,
		q1 = (abs(s1)-conStress)/1000, 
		p = p,
		u= conStress/1000-p, 
		e = -e*100,
		ev = ev*100,
	)
	for b in tot:
		b.shape.color=scalarOnColorScale(b.state.rot().norm(),0,pi/2.) 
	return (dV,e)

# D.b). a function for adding force (servo-controlled of lateral wall) 
def addforce():
	h_sample = O.bodies[top_id].state.pos[2] - O.bodies[bot_id].state.pos[2] #print 'height is ',h_sample
	r_cum = 0
	count = 0
	f4 = 0
	for b in membrane_grid:
		x = b.state.pos[0] 
		y = b.state.pos[1]
		z = b.state.pos[2]
		dist = math.sqrt(x*x+y*y) 
		n = Vector3(x/dist,y/dist,0)
		a = 2*3.1416*dist/(360/alpha)*6*rParticle*(1+e) 
		if z<=O.bodies[top_id].state.pos[2] and z>=O.bodies[bot_id].state.pos[2]: 
			f_local = O.forces.f(b.id)
			length = math.sqrt(f_local[0]*f_local[0]+f_local[1]*f_local[1]+f_local[2]*f_local[2])
			cos_theta = (n[0]*f_local[0]+n[1]*f_local[1]+n[2]*f_local[2])/length
			p_normal = (length*cos_theta/a)
			f4 += (p_normal)
			count += 1 
	p = abs(f4/count/1000) 
	for f in membrane_grid:
		f.shape.color = Vector3(0,0,0) 
		x = f.state.pos[0]
		y = f.state.pos[1]
		z = f.state.pos[2]
		dist = math.sqrt(x*x+y*y)
		n = Vector3(x/dist,y/dist,0) 
		a = 2*3.1416*dist/(360/alpha)*6*rParticle*(1+e)
		if z<=O.bodies[top_id].state.pos[2] and z>=O.bodies[bot_id].state.pos[2]: 
			r_local = dist
			r_cum += r_local 
	r_avg = r_cum/count-rParticle
	Volume = r_avg*r_avg*3.1416*h_sample 
	delV = Volume - V_ini
	for b in topcap:
		O.bodies[b].state.blockedDOFs = 'xyzXYZ' 
		O.bodies[b].state.vel = (0,0,-vel_a)
	for b in botcap:
		O.bodies[b].state.blockedDOFs = 'xyzXYZ' 
		O.bodies[b].state.vel = (0,0,vel_a)
	for f in membrane_grid:
		f.shape.color = Vector3(0,0,0) 
		x = f.state.pos[0]
		y = f.state.pos[1]
		z = f.state.pos[2]
		dist = math.sqrt(x*x+y*y) 
		n = Vector3(x/dist,y/dist,0) 
		a = 2*3.1416*dist/(360/alpha)*6*rParticle#*(1+e)
		if z<=O.bodies[top_id].state.pos[2] and z>=O.bodies[bot_id].state.pos[2]: 
			f.state.blocked = 'z'
			f.state.vel = (dist/h_sample)*vel_a*n 
			if delV>0:
				f.state.vel = -3.0*(dist/h_sample)*vel_a*n
			else:
				f.state.vel = 3.0*(dist/h_sample)*vel_a*n
		else:
			f.state.vel = 0*n

# D.c). a function for recording data 
def checkrecord():
	plot.saveDataTxt('results_'+key)

# D.d). a function used for consolidation 
def confining():
	e_ini = (top[0].state.displ()[2] - bot[0].state.displ()[2])/height_0 
	f1 = sum(O.forces.f(b)[2] for b in topcap)
	f2 = sum(O.forces.f(b)[2] for b in botcap) 
	f4 = 0
	r_cum = 0
	count = 0
	a = 2*3.1416*(.5*width+rParticle)/(360/alpha)*6*rParticle 
	for f in membrane_grid:
		f.shape.color = Vector3(0,0,0)
		x = f.state.pos[0] 
		y = f.state.pos[1]
		z = f.state.pos[2]
		dist = math.sqrt(x*x+y*y) 
		n = Vector3(x/dist,y/dist,0)
		a = 2*3.1416*dist/(360/alpha)*6*rParticle 
		f.state.vel = -vel_ini_r*n 
	for b in topcap:
		O.bodies[b].state.blockedDOFs = 'xyzXYZ' 
		O.bodies[b].state.vel = (0,0,-vel_ini_a)
	for b in botcap:
		O.bodies[b].state.blockedDOFs = 'xyzXYZ' 
		O.bodies[b].state.vel = (0,0,vel_ini_a)
	for b in membrane_grid:
		x = b.state.pos[0] 
		y = b.state.pos[1]
		z = b.state.pos[2]
		dist = math.sqrt(x*x+y*y) 
		n = Vector3(x/dist,y/dist,0)
		a = 2*3.1416*dist/(360/alpha)*6*rParticle 
		if z<=O.bodies[top_id].state.pos[2] and z>=O.bodies[bot_id].state.pos[2]: 
			f_local = O.forces.f(b.id)
			length = math.sqrt(f_local[0]*f_local[0]+f_local[1]*f_local[1]+f_local[2]*f_local[2])
			cos_theta = (n[0]*f_local[0]+n[1]*f_local[1]+n[2]*f_local[2])/length
			p_normal = (length*cos_theta/a) 
			f4 += (p_normal)
			r_cum += dist 
			count += 1
	r_avg = r_cum/count-rParticle 
	fa = abs(.5*(f2-f1))
	p_r = f4/count/1000
	p_a = fa/(3.1416*0.25*width*width)/1000
	e_ini2 = ((O.bodies[top_id].state.pos[2] - O.bodies[bot_id].state.pos[2])- height_0)/height_0
	V_ini = (O.bodies[top_id].state.pos[2] - O.bodies[bot_id].state.pos[2])*r_avg*r_avg*3.1416
	R_avg = r_avg
	H_ini = O.bodies[top_id].state.pos[2] - O.bodies[bot_id].state.pos[2]
	porosity = 1-V_sand/((R_avg)*(R_avg)*3.1416*H_ini)
	return (p_r,p_a,e_ini,V_ini,R_avg,H_ini,porosity)

# D.e). a function for stablization 
def stable():
	for f in membrane_grid:
		x = f.state.pos[0] 
		y = f.state.pos[1]
		z = f.state.pos[2]
		dist = math.sqrt(x*x+y*y) 
		n = Vector3(x/dist,y/dist,0)
		a = 2*3.1416*dist/(360/alpha)*6*rParticle 
		f.state.blockedDOFs = 'xyzXYZ' 
		f.state.vel = 0*n

	for b in topcap:
		O.bodies[b].state.blockedDOFs = 'xyzXYZ' 
		O.bodies[b].state.vel = (0,0,0)

	for b in botcap:
		O.bodies[b].state.blockedDOFs = 'xyzXYZ'
		O.bodies[b].state.vel = (0,0,0)
def compress():
	for b in wall:
		O.bodies[b].state.blockedDOFs = 'xyzXYZ' 
		O.bodies[b].state.vel = (0,0,0)

####################################################################### #	E. APPLYING CONFINING STRESS TO FLEXIBLE MEMBRANE	# #######################################################################
# E.a). adding corrosponding python function 
O.engines=O.engines+[
	PyRunner(command='compress()',iterPeriod=1), 
	PyRunner(command='plotAddData',iterPeriod=1)
]

# E.b). compress until target porosity 
vel_ini_a = rate*height_0
vel_ini_r = rate*height_0 
for b in topcap:
	O.bodies[b].state.blockedDOFs = 'xyzXYZ' 
	O.bodies[b].state.vel = (0,0,-vel_ini_a)
for b in botcap:
	O.bodies[b].state.blockedDOFs = 'xyzXYZ' 
	O.bodies[b].state.vel = (0,0,vel_ini_a)
while 1:
	O.run(100,True)
	h = (O.bodies[top_id].state.pos[2]-O.bodies[bot_id].state.pos[2])
	V = h*0.25*width*width*3.1416 
	porosity = 1-V_sand/V
	print('porosity: ',porosity, ' height: ', h)
	if (porosity <= targetPorosity):
		print( 'compression stage finished!')
		break
h_ini = h # record height after compression

# E.c). ADD CONFINING STRESS
# E.c.1). remove facet wall 
for b in wall:
	O.bodies.erase(b)

# E.c.2). create membrane around particle packing by reading mesh file 
def membraneCylinderGenerator(width, height, nw, nh, rNode, posC, materialNodes, materialPfacets):
	'''
	width: diameter of cylinder
	height: height of cylinder
	nw: number of facets along perimeter
	nh: number of facets along height
	rNode: radius of grid node
	posC: center of cylindrical bottom surface-->Vector3(x,y,z)
	materialNodes: material of grid nodes
	materialPfacets: material of PFacets
	'''
	nodesIds = [] # ID list of gridnodes
	cylIds = []   # ID list of cylinders
	pfIds = []    # ID list of pfacets
	colornode = [255/255,102/255,0/255]
	colorcylinder = [0,0,0]
	colorpfacet1 = [248/255,248/255,168/255]
	colorpfacet2 = [156/255,160/255,98/255]
	# create gridnodes
	rCyl2 = 0.5*width / cos(pi/float(nw))
	for r in range(nw):
 		for h in range(nh+1):
    			v = Vector3(rCyl2*cos(2*pi*r/float(nw)),rCyl2*sin(2*pi*r/float(nw)), height*h/float(nh))+posC
    			V = (O.bodies.append(gridNode(v,rNode,wire=False,fixed=False, material=materialNodes,color=colornode)) )
    			nodesIds.append(V)
	# create grid connection
	nodeIdsMin = min(nodesIds)
	for r in range(nw):
		for h in range(nh):
			if r == nw-1:
				V1 = nodeIdsMin+(r+0)*(nh+1)+h+0
				V2 = nodeIdsMin+(r+0)*(nh+1)+h+1
				V3 = nodeIdsMin+(0+0)*(nh+1)+h+1
				V4 = nodeIdsMin+(0+0)*(nh+1)+h+0
			else:
				V1 = nodeIdsMin+(r+0)*(nh+1)+h+0
				V2 = nodeIdsMin+(r+0)*(nh+1)+h+1
				V3 = nodeIdsMin+(r+1)*(nh+1)+h+1
				V4 = nodeIdsMin+(r+1)*(nh+1)+h+0
			#create grid connection 
			cylIds.append(O.bodies.append(gridConnection(V1,V2,rNode,material=materialNodes,color=colorcylinder)))
			cylIds.append(O.bodies.append(gridConnection(V1,V4,rNode,material=materialNodes,color=colorcylinder)))
			cylIds.append(O.bodies.append(gridConnection(V1,V3,rNode,material=materialNodes,color=colorcylinder)))
			if h == nh-1:
				cylIds.append(O.bodies.append(gridConnection(V2,V3,rNode,material=materialNodes,color=colorcylinder)))
	# create PFacet
	for r in range(nw):
		for h in range(nh):
			if r == nw-1:
				V1 = nodeIdsMin+(r+0)*(nh+1)+h+0
				V2 = nodeIdsMin+(r+0)*(nh+1)+h+1
				V3 = nodeIdsMin+(0+0)*(nh+1)+h+1
				V4 = nodeIdsMin+(0+0)*(nh+1)+h+0
			else:
				V1 = nodeIdsMin+(r+0)*(nh+1)+h+0
				V2 = nodeIdsMin+(r+0)*(nh+1)+h+1
				V3 = nodeIdsMin+(r+1)*(nh+1)+h+1
				V4 = nodeIdsMin+(r+1)*(nh+1)+h+0
			#create Pfacets
			pfIds.append(O.bodies.append(pfacet(V1,V3,V2,wire=True,material=materialPfacets,color=colorpfacet1)))
			pfIds.append(O.bodies.append(pfacet(V1,V4,V3,wire=True,material=materialPfacets,color=colorpfacet2)))
	return nodesIds,cylIds,pfIds
# 4.2.2). generate flexible membrane
shiftfactor = O.bodies[bot_id].state.pos[2]-((height-h_ini)/2) 
nodesIds,cylIds,pfIds = membraneCylinderGenerator(width=width, height=O.bodies[top_id].state.pos[2]-O.bodies[bot_id].state.pos[2], nw=40, nh=25, rNode=width/100,
posC=Vector3(0,0,O.bodies[bot_id].state.pos[2]), materialNodes=GridMat, materialPfacets=pFacetMat)

'''
shiftfactor = O.bodies[bot_id].state.pos[2]-((height-h_ini)/2) 
nodesIds,cylIds,pfIds = gtsPFacet(
'Mesh_cylinder.gts',shift=(0,0,shiftfactor),scale=1,radius=rParticle,wire=False,fixed=False,
materialNodes='gridNodeMat',material='pFacetMat',color=Vector3(0.5,1,0.5)
)
'''
# E.c.3). define different sections of membrane 
membrane_grid = [O.bodies[s] for s in nodesIds ] 
membrane_pfacet = [O.bodies[s] for s in pfIds]

# E.c.4). run one interaction 
for f in membrane_grid:
	f.shape.color = Vector3(0,0,0) 
	x = f.state.pos[0]
	y = f.state.pos[1]
	z = f.state.pos[2]
	dist = math.sqrt(x*x+y*y) 
	n = Vector3(x/dist,y/dist,0)
	if z<=O.bodies[top_id].state.pos[2] and z>=O.bodies[bot_id].state.pos[2]: 
		f.state.vel = -vel_ini_r*n
O.engines[2].physDispatcher.functors[1].setCohesionNow=True 
O.engines[5].dead = True
while 1:
	O.run(1,True) 
	break

# E.c.5). redefine engine #**********************************************************************# 
O.engines=[
	ForceResetter(), 
	InsertionSortCollider(
	[
		Bo1_Sphere_Aabb(), 
		Bo1_Facet_Aabb(),
		Bo1_PFacet_Aabb(),
		Bo1_GridConnection_Aabb()
	]
	),
	InteractionLoop( 
	[
		Ig2_Sphere_Sphere_ScGeom6D(), 
		Ig2_GridNode_GridNode_GridNodeGeom6D(), 
		Ig2_GridConnection_GridConnection_GridCoGridCoGeom(), 
		Ig2_Sphere_PFacet_ScGridCoGeom(), 
		Ig2_Facet_Sphere_ScGeom6D()
	], 
	[
		Ip2_FrictMat_FrictMat_FrictPhys(), 
		Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(label="cohesiveIp")
	], 
	[
		Law2_ScGeom_FrictPhys_CundallStrack(), 
		Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm=True,always_use_moment_law=False,label='cohesiveLaw'), 
		Law2_GridCoGridCoGeom_FrictPhys_CundallStrack(), Law2_ScGridCoGeom_FrictPhys_CundallStrack(),
	],
	label="iloop"
	),
	GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=25,timestepSafetyCoefficient=0.8),
	NewtonIntegrator(gravity=(0,0,0),damping=0.3,label='newton'), 
	PyRunner(command='confining()',iterPeriod=1), 
	PyRunner(command='plotAddData',iterPeriod=1)
]
#**********************************************************************#

# set final friction angle, enable cohesion 
setContactFriction(radians(finalFricDegree)) 
O.engines[2].lawDispatcher.functors[1].always_use_moment_law=True 
O.engines[2].physDispatcher.functors[1].setCohesionNow=True 
O.engines[2].physDispatcher.functors[1].setCohesionOnNewContacts=False

# E.c.6). confining
# some initial parameters 
p_a = 0
p_r = 0
e_ini = 0
V_ini = 0
R_avg = 0
H_ini = 0
 
porosity = 0 
# velocity
vel_ini_a = 0.05*rate*height_0 
vel_ini_r = 0.05*rate*height_0
# loops (fast-slow) for reaching target confining stress 
while 1:
	O.run(10,True) 
	(p_r,p_a,e_ini,V_ini,R_avg,H_ini,porosity)=confining() 
	p_r = abs(p_r)
	pressure = max(p_r,p_a) 
	dif = p_r-p_a
	if (p_a > isoStress/1000): 
		vel_ini_a = -abs(vel_ini_a) 
		if (p_r > isoStress/1000):
			vel_ini_r = -abs(vel_ini_r)
		else:
			vel_ini_r = abs(vel_ini_r)
	elif (p_a <= isoStress/1000):
		if (p_r > isoStress/1000): 
			vel_ini_a = abs(vel_ini_a) 
			vel_ini_r = -abs(vel_ini_r)
		else:
 			if (pressure<0.9*isoStress/1000): 
 				if dif > 5:
 					vel_ini_a = 1.05*abs(vel_ini_a)
 				elif dif < -5:
 					vel_ini_r = 1.05*abs(vel_ini_r)
 			if (pressure>=0.85*isoStress/1000 and pressure<=isoStress/1000): 
 				if dif > 1:
 					if dif > 5:
 						vel_ini_a = 1.5*abs(vel_ini_a)
 					else:
 						vel_ini_a = 1.01*abs(vel_ini_a)
 				elif dif < -1:
 					if dif < -5:
 						vel_ini_r = 1.5*abs(vel_ini_r)
 					else:
 						vel_ini_r = 1.01*abs(vel_ini_r)
	mean = (p_r+p_a)/2 
	unb=unbalancedForce()
	print( 'p= ',p_r,' q= ',p_a,' porosity= ',porosity,' unbalanced force: ',unb,' vr:',vel_ini_r,' va:',vel_ini_a)
	if abs(isoStress/1000-mean)/(isoStress/1000)<0.15 and abs(dif) <15:
		print( 'initial strain: ',e_ini) 
		print( 'initial volume: ',V_ini)
		print( 'Confining stage I is finished!')
		break
while 1:
	for f in membrane_grid:
		f.state.blockedDOFs = 'xyzXYZ' 
	O.run(1,True)
	(p_r,p_a,e_ini,V_ini,R_avg,H_ini,porosity)=confining() 
	p_r = abs(p_r)
	pressure = max(p_r,p_a) 
	dif = p_r-p_a
	if (p_a > isoStress/1000): 
		vel_ini_a = -abs(vel_ini_a) 
		if (p_r > isoStress/1000):
			vel_ini_r = -abs(vel_ini_r)
		else:
			vel_ini_r = abs(vel_ini_r)
	elif (p_a <= isoStress/1000):
		if (p_r > isoStress/1000): 
			vel_ini_a = abs(vel_ini_a) 
			vel_ini_r = -abs(vel_ini_r)
		else:
			if (pressure<0.9*isoStress/1000):
				if dif > 5:
					vel_ini_a = 1.05*abs(vel_ini_a) 
				elif dif < -5:
					vel_ini_r = 1.05*abs(vel_ini_r)
			if (pressure>=0.9*isoStress/1000 and pressure<=isoStress/1000): 
				if dif > 1:
					if dif > 5:
						vel_ini_a = 1.1*abs(vel_ini_a)
					else:
						vel_ini_a = 1.01*abs(vel_ini_a)
				elif dif < -1:
					if dif < -5:
						vel_ini_r = 1.1*abs(vel_ini_r)
					else:
						vel_ini_r = 1.01*abs(vel_ini_r)
	mean = (p_r+p_a)/2 
	unb=unbalancedForce()
	print( 'p= ',p_r,' q= ',p_a,' porosity= ',porosity,' unbalanced force: ',unb,' vr:',vel_ini_r,' va:',vel_ini_a )
	if abs(isoStress/1000-mean)/(isoStress/1000)<0.05 and abs(dif) <10:
		print( 'initial strain: ',e_ini )
		print( 'initial volume: ',V_ini )
		print( 'Confining stage II is finished!' )
		break
print( 'V_ini= ',V_ini )
 
# E.c.7). stablize #**********************************************************************# 
O.engines=[
	ForceResetter(), 
	InsertionSortCollider([
		Bo1_Sphere_Aabb(), 
		Bo1_Facet_Aabb(),
		Bo1_PFacet_Aabb(),  
		Bo1_GridConnection_Aabb()
	]),
	InteractionLoop(
	[
		Ig2_Sphere_Sphere_ScGeom6D(), 
		Ig2_GridNode_GridNode_GridNodeGeom6D(), 
		Ig2_GridConnection_GridConnection_GridCoGridCoGeom(), 
		Ig2_Sphere_PFacet_ScGridCoGeom(), 
		Ig2_Facet_Sphere_ScGeom6D()
	], 
	[
		Ip2_FrictMat_FrictMat_FrictPhys(), 
		Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(label="cohesiveIp")
	], 
	[
		Law2_ScGeom_FrictPhys_CundallStrack(), 
		Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm=True,always_use_moment_law=False,label='cohesiveLaw'), 
		Law2_GridCoGridCoGeom_FrictPhys_CundallStrack(), 
		Law2_ScGridCoGeom_FrictPhys_CundallStrack(),
	],
	label="iloop"
	),
	GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=25,timestepSafetyCoefficient=0.8),
	NewtonIntegrator(gravity=(0,0,0),damping=0.3,label='newton'), 
	PyRunner(command='stable()',iterPeriod=1), 
	PyRunner(command='plotAddData',iterPeriod=1)
]
#**********************************************************************# 
O.engines[2].lawDispatcher.functors[1].always_use_moment_law=True 
O.engines[2].physDispatcher.functors[1].setCohesionNow=True 
O.engines[2].physDispatcher.functors[1].setCohesionOnNewContacts=False
vel_a = abs(vel_ini_a) 
while 1:
	O.run(100,True) 
	unb=unbalancedForce()
	print( 'unbalanced force: ',unb)
	e_ini = (top[0].state.displ()[2] - bot[0].state.displ()[2])/height_0 
	f1 = sum(O.forces.f(b)[2] for b in topcap)
	f2 = sum(O.forces.f(b)[2] for b in botcap) 
	f4 = 0
	r_cum = 0
	count = 0
	for b in membrane_grid:
		x = b.state.pos[0] 
		y = b.state.pos[1]
		z = b.state.pos[2]
		dist = math.sqrt(x*x+y*y) 
		n = Vector3(x/dist,y/dist,0)
		a = 2*3.1416*dist/(360/alpha)*6*rParticle#*(1+e) 
		if z<=O.bodies[top_id].state.pos[2] and z>=O.bodies[bot_id].state.pos[2]: 
			f_local = O.forces.f(b.id)
			length = math.sqrt(f_local[0]*f_local[0]+f_local[1]*f_local[1]+f_local[2]*f_local[2])
			cos_theta = (n[0]*f_local[0]+n[1]*f_local[1]+n[2]*f_local[2])/length
			p_normal = (length*cos_theta/a) 
			f4 += (p_normal)
			r_cum += dist 
			count += 1
	r_avg = r_cum/count-rParticle 
	fa = abs(.5*(f2-f1))
	p_r = f4/count/1000
	p_a = fa/(3.1416*r_avg*r_avg)/1000 
	print( 'pr=', p_r, ' pa=',p_a )
	if unb <= thresholdvalue: 
		break

####################################################################### #	F. APPLYING DEVIATOR LOADING	# #######################################################################
# F.a). redifine engines 
O.engines=[
	ForceResetter(), 
	InsertionSortCollider(
		[
		Bo1_Sphere_Aabb(aabbEnlargeFactor=enlargefactor), 
		Bo1_Facet_Aabb(), 
		Bo1_PFacet_Aabb(),
		Bo1_GridConnection_Aabb()
		]
	),
	InteractionLoop( 
		[
		Ig2_Sphere_Sphere_ScGeom6D(interactionDetectionFactor=enlargefactor), 
		Ig2_GridNode_GridNode_GridNodeGeom6D(), 
		Ig2_GridConnection_GridConnection_GridCoGridCoGeom(), 
		Ig2_Sphere_PFacet_ScGridCoGeom(), 
		Ig2_Facet_Sphere_ScGeom6D()
		], 
		[
		Ip2_FrictMat_FrictMat_FrictPhys(), 
		Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(label="cohesiveIp")
		], 
		[
		Law2_ScGeom_FrictPhys_CundallStrack(), 
		Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm=True,always_use_moment_law=False,label='cohesiveLaw'), 
		Law2_GridCoGridCoGeom_FrictPhys_CundallStrack(), 
		Law2_ScGridCoGeom_FrictPhys_CundallStrack(),
		],
		label="iloop"
	),
	GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=25,timestepSafetyCoefficient=0.8),
	NewtonIntegrator(gravity=(0,0,0),damping=0.3,label='newton')
]

O.engines[2].lawDispatcher.functors[1].always_use_moment_law=True 
O.engines[2].physDispatcher.functors[1].setCohesionNow=True 
O.engines[2].physDispatcher.functors[1].setCohesionOnNewContacts=False 
O.engines=O.engines+[
	PyRunner(command='addforce()',iterPeriod=1,label='force'), 
	PyRunner(command='plotAddData()',iterPeriod=1,label='recorder'), 
	PyRunner(command='checkrecord()',realPeriod=10,label='checker')
]

# F.b). define the velocity of membrane walls to maintain the volume constant condition 
vel_a = final_rate*abs(vel_ini_a)
vel_r = vel_a*.5*width/height 
Vel_r = vel_r
conStress = p_r*1000 ####################################################################### #	G. UTILITIES	# #######################################################################
# G.a). time step (recommanded by YADE) 
O.dt=0.3*PWaveTimeStep()
t = O.dt

# G.b). funtion for plot
 
plot.plots={'e':('q','p'),'e':('u','ev')} 
plot.plot()
O.saveTmp()
O.timingEnabled=1 #################################################################### #	G. RUN	# #################################################################### 
print ('=========================')
print ("start triaxial simulation")
print ('=========================')
O.run()


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