# yade-users team mailing list archive

## [Question #695192]: wall.append(O.bodies.append(pfacet) in triaxial test

```New question #695192 on Yade:

Hi,

I am doing a triaxial test simulating fiber-reinforced sand in which cylinder and Pfacet are used as fiber and wall(triax.wall). After setting strain to be -0.9 I find that Pfacet doesn't move along coordinate axis as regular triaxial wall. Pfacet can deform and try to fold each other. I need pfacet to move along coordinate axis. can anyone help me with that? I am using ubuntu18.04 and my yade version is 2018.02b.

###################################################
#from random import random as r
import numpy as np
from numpy import *
import math

#### Parameter ####
mi,ma = (-100e-3,-100e-3,-200e-3),(100e-3,100e-3,200e-3)
young=4.0e6
poisson=3
density=1e3
color=[0.,1.,1.]
stabilityThreshold=0.01
#=============================meterials========================================
O.materials.append(FrictMat(young=4.0e6,poisson=.3,frictionAngle=30,density=2630e0,label='sphereMat'))#for sphere

O.materials.append(FrictMat(young=3e9,poisson=.15,frictionAngle=20,density=910e+0,label='extcylMat'))#for sphere-cylinder
O.materials.append(CohFrictMat(young=3e9,poisson=.15,density=910e0,frictionAngle=20,normalCohesion=1e40,shearCohesion=1e40,momentRotationLaw=True,label='intcylMat'))#for cylinder-cylinder

O.materials.append( CohFrictMat( young=3e6,poisson=0.15,density=910e0,frictionAngle=20,normalCohesion=3e100,shearCohesion=3e100,momentRotationLaw=True,label='gridNodeMat' ) )#for gridNodes
#O.materials.append(CohFrictMat(young=3e9,poisson=.15,density=910e6,frictionAngle=20,normalCohesion=1e40,shearCohesion=1e40,momentRotationLaw=True,label='gridNodeMat'))#for gridNodes
O.materials.append(FrictMat(young=4e6,poisson=0.3,density=1000,frictionAngle=20,label='pFacetMat')) #for pfacet

O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=2630,label='walls'))

triax=TriaxialStressController(

maxMultiplier=1.+2e4/young, # spheres growing factor (fast growth)
finalMaxMultiplier=1.+2e3/young, # spheres growing factor (slow growth)
thickness = 0,
max_vel = 3.7e-7,
wall_front_activated=True,
wall_back_activated=True,
wall_top_activated=True,
wall_bottom_activated=True,
wall_left_activated=True,
wall_right_activated=True,

internalCompaction=False, # If true the confining pressure is generated by growing particles
)

#==============================Engines=========================================
O.engines=[
ForceResetter(),
InsertionSortCollider([
Bo1_Sphere_Aabb(),
Bo1_GridConnection_Aabb(),
Bo1_PFacet_Aabb(),
Bo1_Box_Aabb(),
Bo1_Wall_Aabb()
]),
InteractionLoop([
Ig2_Sphere_Sphere_ScGeom(),
Ig2_Box_Sphere_ScGeom(),
Ig2_GridNode_GridNode_GridNodeGeom6D(),
Ig2_Sphere_GridConnection_ScGridCoGeom(),
Ig2_GridConnection_GridConnection_GridCoGridCoGeom(),
Ig2_GridConnection_PFacet_ScGeom(),
Ig2_PFacet_PFacet_ScGeom(),
Ig2_Sphere_PFacet_ScGridCoGeom(),
Ig2_Wall_PFacet_ScGeom(),
Ig2_Wall_Sphere_ScGeom()
],
[
Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True,setCohesionOnNewContacts=False),	# internal cylinder physics
Ip2_FrictMat_FrictMat_FrictPhys(),	# physics for external interactions, i.e., cylinder-cylinder, sphere-sphere, cylinder-sphere
Ip2_FrictMat_CpmMat_FrictPhys(),
],
[
Law2_ScGeom_FrictPhys_CundallStrack(),	# contact law for sphere-sphere
Law2_ScGridCoGeom_FrictPhys_CundallStrack(),	# contact law for cylinder-sphere
Law2_ScGeom6D_CohFrictPhys_CohesionMoment(),	# contact law for "internal" cylinder forces
Law2_GridCoGridCoGeom_FrictPhys_CundallStrack(),	# contact law for cylinder-cylinder interaction
#Law2_GridCoGridCoGeom_CohFrictPhys_CundallStrack(),
Law2_ScGeom_CpmPhys_Cpm(),
]
),
GlobalStiffnessTimeStepper(timestepSafetyCoefficient=0.1,label='ts'),
triax,
NewtonIntegrator(gravity=(0,0,0),damping=0.5,label='newton'),
]

triax.goal1=triax.goal2=triax.goal3=-0.9

RNode=0.00002
nodes=[]
wall=[]
nodes.append(O.bodies.append( gridNode([-100e-3, 100e-3,-200e-3],RNode,wire=False,fixed=True,material='gridNodeMat',color=color) ))
nodes.append(O.bodies.append( gridNode([ 100e-3, 100e-3,-200e-3],RNode,wire=False,fixed=True,material='gridNodeMat',color=color) ))
nodes.append(O.bodies.append( gridNode([ 100e-3,-100e-3,-200e-3],RNode,wire=False,fixed=True,material='gridNodeMat',color=color) ))
nodes.append(O.bodies.append( gridNode([-100e-3,-100e-3,-200e-3],RNode,wire=False,fixed=True,material='gridNodeMat',color=color) ))

nodes.append(O.bodies.append( gridNode([-100e-3, 100e-3, 200e-3],RNode,wire=False,fixed=True,material='gridNodeMat',color=color) ))
nodes.append(O.bodies.append( gridNode([ 100e-3, 100e-3, 200e-3],RNode,wire=False,fixed=True,material='gridNodeMat',color=color) ))
nodes.append(O.bodies.append( gridNode([ 100e-3,-100e-3, 200e-3],RNode,wire=False,fixed=True,material='gridNodeMat',color=color) ))
nodes.append(O.bodies.append( gridNode([-100e-3,-100e-3, 200e-3],RNode,wire=False,fixed=True,material='gridNodeMat',color=color) ))

O.bodies.append(gridConnection(0,1,RNode,color=color,material='gridNodeMat'))
O.bodies.append(gridConnection(0,2,RNode,color=color,material='gridNodeMat'))
O.bodies.append(gridConnection(0,3,RNode,color=color,material='gridNodeMat'))
O.bodies.append(gridConnection(0,4,RNode,color=color,material='gridNodeMat'))
O.bodies.append(gridConnection(0,5,RNode,color=color,material='gridNodeMat'))
O.bodies.append(gridConnection(0,7,RNode,color=color,material='gridNodeMat'))

O.bodies.append(gridConnection(6,1,RNode,color=color,material='gridNodeMat'))
O.bodies.append(gridConnection(6,2,RNode,color=color,material='gridNodeMat'))
O.bodies.append(gridConnection(6,3,RNode,color=color,material='gridNodeMat'))
O.bodies.append(gridConnection(6,4,RNode,color=color,material='gridNodeMat'))
O.bodies.append(gridConnection(6,5,RNode,color=color,material='gridNodeMat'))
O.bodies.append(gridConnection(6,7,RNode,color=color,material='gridNodeMat'))

O.bodies.append(gridConnection(2,1,RNode,color=color,material='gridNodeMat'))
O.bodies.append(gridConnection(2,3,RNode,color=color,material='gridNodeMat'))
O.bodies.append(gridConnection(2,5,RNode,color=color,material='gridNodeMat'))
O.bodies.append(gridConnection(2,7,RNode,color=color,material='gridNodeMat'))

O.bodies.append(gridConnection(4,1,RNode,color=color,material='gridNodeMat'))
O.bodies.append(gridConnection(4,3,RNode,color=color,material='gridNodeMat'))
O.bodies.append(gridConnection(4,7,RNode,color=color,material='gridNodeMat'))
O.bodies.append(gridConnection(4,5,RNode,color=color,material='gridNodeMat'))

O.bodies.append(gridConnection(3,7,RNode,color=color,material='gridNodeMat'))
O.bodies.append(gridConnection(5,1,RNode,color=color,material='gridNodeMat'))
O.bodies.append(gridConnection(5,7,RNode,color=color,material='gridNodeMat'))
O.bodies.append(gridConnection(3,1,RNode,color=color,material='gridNodeMat'))

wall.append(O.bodies.append(pfacet(nodes[0],nodes[1],nodes[2],wire=True,material='pFacetMat',color=color)))
wall.append(O.bodies.append(pfacet(nodes[0],nodes[2],nodes[3],wire=True,material='pFacetMat',color=color)))

wall.append(O.bodies.append(pfacet(nodes[0],nodes[4],nodes[5],wire=True,material='pFacetMat',color=color)))
wall.append(O.bodies.append(pfacet(nodes[5],nodes[1],nodes[0],wire=True,material='pFacetMat',color=color)))

wall.append(O.bodies.append(pfacet(nodes[0],nodes[4],nodes[7],wire=True,material='pFacetMat',color=color)))
wall.append(O.bodies.append(pfacet(nodes[0],nodes[7],nodes[3],wire=True,material='pFacetMat',color=color)))

wall.append(O.bodies.append(pfacet(nodes[6],nodes[5],nodes[4],wire=True,material='pFacetMat',color=color)))
wall.append(O.bodies.append(pfacet(nodes[6],nodes[4],nodes[7],wire=True,material='pFacetMat',color=color)))

wall.append(O.bodies.append(pfacet(nodes[6],nodes[5],nodes[1],wire=True,material='pFacetMat',color=color)))
wall.append(O.bodies.append(pfacet(nodes[6],nodes[1],nodes[2],wire=True,material='pFacetMat',color=color)))

wall.append(O.bodies.append(pfacet(nodes[6],nodes[7],nodes[3],wire=True,material='pFacetMat',color=color)))
wall.append(O.bodies.append(pfacet(nodes[6],nodes[3],nodes[2],wire=True,material='pFacetMat',color=color)))

#=======================================sphere=================================================
sp.makeCloud(mi,ma,porosity=0.66/1.66,psdSizes=[0.5e-3,0.75e-3,0.8e-3,0.9e-3,1.2e-3,2e-3],psdCumm=[0.,0.2,0.4,0.6,0.8,1.0],num=133)
spheres=sp.toSimulation(color=(0,0.3,0.7),material='sphereMat')

#=======================================fiber=================================================
#below codes just to generate cylinders that don't overlap with each other.
#I tried to simplify it somehow error occurred in terminal.
#I think below codes have nothing to do with my problem. Thanks!
length=12e-3
diameter=0.3e-3
Xw=0.25e-2 # fiber content
fiber_vol=(80e-3*(40e-3*0.5)**2*np.pi)*Xw
Vfiber=length*np.pi*(0.5*diameter)**2
numFiber=fiber_vol/Vfiber
print "the numer of fibres =",int(numFiber)
Ne=int(length/(0.9e-3))
print 'Ne=',Ne
nodesIds=[]
cylIds=[]
numFibre=0
target_num=int(numFiber)
target_num=int(6)
####fiberfilewrite
fiberwr=open('fibers.txt','w')
cylNodes=[]
time=0
HalleluYah=0
for n in range(1000):
if numFibre<target_num:
random.seed(n*5)
z0=random.uniform(-200e-3, 200e-3 )
random.seed( )
sita0=random.uniform(-np.pi,np.pi)
random.seed( n+7)
l_y=random.uniform(-100e-3,100e-3)
x0=l_y*sin(sita0)
y0=l_y*cos(sita0)

random.seed(n+8)
sita=random.uniform(-0.5*np.pi, 0.5*np.pi)
random.seed()
phi=random.uniform(0, 2*np.pi)

hz=length*sin(sita)
hx=length*cos(sita)*cos(phi)
hy=length*cos(sita)*sin(phi)

x1=x0+hx
y1=y0+hy
z1=z0+hz

if 0<abs(x1)<100e-3 and 0<abs(y1)<100e-3 and 0<abs(z1)<200e-3 and 0<abs(x0)<100e-3 and 0<abs(y0)<100e-3 and 0<abs(z0)<200e-3 :
if len(cylNodes)<=0:
fiberwr.write(str(numFibre)+'\t'+str(x0)+'\t'+str(y0)+'\t'+str(z0)+'\t'+str(x1)+'\t'+str(y1)+'\t'+str(z1)+'\n')
cylNodes.append([x0,y0,z0,x1,y1,z1])
numFibre+=1
vertices=[]
for i in range(0, Ne+1):
px=float(i)*hx/float(Ne)+x0; py=float(i)*hy/float(Ne)+y0; pz=float(i)*hz/float(Ne)+z0;
vertices.append([px,py,pz])
cylinderConnection(vertices,0.15e-3,nodesIds,cylIds,color=[0.5,0.5,0],fixed=False,highlight=False,intMaterial='intcylMat',extMaterial='extcylMat')
else :
valid=True
for c in cylNodes:
if valid:
A1=np.array([x0,y0,z0])
B1=np.array([x1,y1,z1])
r1=0.15e-3
A2=np.array([c[0],c[1],c[2]])
B2=np.array([c[3],c[4],c[5]])
r2=0.15e-3
l1=B1-A1
l2=B2-A2
l3=A2-A1

delta1=(np.dot(l1,l2))/(np.linalg.norm(l1)*np.linalg.norm(l2))

if abs(abs(delta1)-1)<0.0000000000000001:
p=np.linalg.norm(np.cross(l2,l3))/np.linalg.norm(l2)
if p<r1+r2 :
print("PARALLEL INTERSECT!")
valid=False
else :
print("PARALLEL DO NOT INTERSECT!")
valid=True
else :
l4=np.cross(l1,l2)
ll1=np.cross(l1,l4)
ll2=np.cross(l2,l4)

A11=np.matrix([ll1,[l2[1],-l2[0],0.0],[0.0,l2[2],-l2[1]]])
x11=np.dot(ll1,A1)
x12=l2[1]*A2[0]-l2[0]*A2[1]
x13=l2[2]*A2[1]-l2[1]*A2[2]
b1=np.matrix([[x11],[x12],[x13]])
C1=np.linalg.solve(A11,b1)

A21=np.matrix([ll2,[l1[1],-l1[0],0.0],[0.0,l1[2],-l1[1]]])
x21=np.dot(ll2,A2)
x22=l1[1]*A1[0]-l1[0]*A1[1]
x23=l1[2]*A1[1]-l1[1]*A1[2]
b2=np.matrix([[x21],[x22],[x23]])
C2=np.linalg.solve(A21,b2)

t1=(C1[0]-A1[0])/(B1[0]-A1[0])
t2=(C2[0]-A2[0])/(B2[0]-A2[0])

if t1<0 :
E=A1
elif t1>=0 and t1<=1 :
E=C1
else :
E=B1

if t2<0 :
F=A2
elif t2>=0 and t2<=1 :
F=C2
else :
F=B2
D=np.linalg.norm(E-F)

if D<(r1 + r2) :
print("INTERSECT!")
valid=False

else :
valid=True

if valid:
print(n,valid)
fiberwr.write(str(numFibre)+'\t'+str(x0)+'\t'+str(y0)+'\t'+str(z0)+'\t'+str(x1)+'\t'+str(y1)+'\t'+str(z1)+'\n')
cylNodes.append([x0,y0,z0,x1,y1,z1])
numFibre+=1
vertices=[]
for i in range(0, Ne+1):
px=float(i)*hx/float(Ne)+x0; py=float(i)*hy/float(Ne)+y0; pz=float(i)*hz/float(Ne)+z0;
vertices.append([px,py,pz])
cylinderConnection(vertices,0.15e-3,nodesIds,cylIds,color=[0.5,0.5,0],fixed=False,highlight=False,intMaterial='intcylMat',extMaterial='extcylMat')
#print(n,"A1",x0,y0,z0, "B1",x1,y1,z1)

fiberwr.close()

Thanks!

--