# wall.append(O.bodies.append(pfacet) in triaxial test

Asked by Hanying Zhang on 2021-01-26

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!

## Question information

Language:
English Edit question
Status:
Solved
For:
Assignee:
No assignee Edit question
Solved by:
Jérôme Duriez
Solved:
2021-02-10
Last query:
2021-02-10
2021-02-02
 Hanying Zhang (xxxe) said on 2021-01-26: #1

I set Pfacet attribute wire=true to check body moving better, but it seems that the triax engine doesn‘t work. So is it workable to use Pfacet simulating triaxial test by set "wall.append(O.bodies.append(pfacet(nodes[0],nodes[1],nodes[2],wire=True,material='pFacetMat',color=color)))"? Thanks!!

 Bruno Chareyre (bruno-chareyre) said on 2021-01-26: #2

Hi,
I would suggest to use boxes for the walls since they are completely flat. It would minimize the number of problems for a start.

Using pfacets should work but not the way you do it.
First problem: triax.wallIds is undefined, so triax engine doesn't know where to find walls .
Second problem: you can't impose the motion of a PFacet, it always follow the motion of the nodes. One exception is if you clump facet+nodes into a clump for each of the walls. Then the clumps' ids should be used to define wallIds.
I hope it helps
Bruno

 xxxea (xxxea) said on 2021-01-27: #3

Hi，Bruno.
> I would suggest to use boxes for the walls since they are completely flat. It would minimize the number of problems for a start.
I have tried boxes for the walls and it worked well except boxes can't interact with cylinder, and that's why I use Pfacet for walls.

> First problem: triax.wallIds is undefined, so triax engine doesn't know where to find walls .
So I should define a wallIDs list then append Pfacet-wall to it?

> Second problem: you can't impose the motion of a PFacet, it always follow the motion of the nodes. One exception is if you clump facet+nodes into a clump for each of the walls. Then the clumps' ids should be used to define wallIds.
One question, there are 8 nodes and I need to generate 12 Pfacets, so every node would be reused if I clump pfacet+nodes into a clump for each of the walls. Is that workable?

Besides, I find only Ig2_GridConnection_PFacet_ScGeom can generate contact between cylinder and other kind of wall-like body, and I have been working on fiber-reinforced sand so cylinder is the kind of body that I must use. Don't really know is there any other body that I can use to simulate wall?

> -----原始邮件-----
> 发件人: "Bruno Chareyre" <email address hidden>
> 发送时间: 2021-01-27 02:41:01 (星期三)
> 抄送:
> 主题: Re: [Question #695192]: wall.append(O.bodies.append(pfacet) in triaxial test
>
>
>
> Bruno Chareyre proposed the following answer:
> Hi,
> I would suggest to use boxes for the walls since they are completely flat. It would minimize the number of problems for a start.
>
> Using pfacets should work but not the way you do it.
> First problem: triax.wallIds is undefined, so triax engine doesn't know where to find walls .
> Second problem: you can't impose the motion of a PFacet, it always follow the motion of the nodes. One exception is if you clump facet+nodes into a clump for each of the walls. Then the clumps' ids should be used to define wallIds.
> I hope it helps
> Bruno
>
> --
> know that it is solved:
>
> If you still need help, you can reply to this email or go to the
> following page to enter your feedback:
>

 Hanying Zhang (xxxe) said on 2021-01-27: #4

Hi，Bruno.
> I would suggest to use boxes for the walls since they are completely flat. It would minimize the number of problems for a start.
I have tried boxes for the walls and it worked well except boxes can't interact with cylinder, and that's why I use Pfacet for walls.

> First problem: triax.wallIds is undefined, so triax engine doesn't know where to find walls .
So I should define a wallIDs list then append Pfacet-wall to it?

> Second problem: you can't impose the motion of a PFacet, it always follow the motion of the nodes. One exception is if you clump facet+nodes into a clump for each of the walls. Then the clumps' ids should be used to define wallIds.
One question, there are 8 nodes and I need to generate 12 Pfacets, so every node would be reused if I clump pfacet+nodes into a clump for each of the walls. Is that workable?

Besides, I find only Ig2_GridConnection_PFacet_ScGeom can generate contact between cylinder and other kind of wall-like body, and I have been working on fiber-reinforced sand so cylinder is the kind of body that I must use. Don't really know is there any other body that I can use to simulate wall?
Or, can I just impose the motion of a Pfacet without triax engine?Is there any attribute like O.area.force(I make up it)?
My idea is that:
(1)Generate sample by Pfacets pushing spheres and cylinders till they achieve some certain positions(in a smaller box)
(2)triaxial test on sample.

 Bruno Chareyre (bruno-chareyre) said on 2021-01-27: #5

> I have tried boxes for the walls and it worked well except boxes can't interact with cylinder

Sure of that? Any example script?

> One question, there are 8 nodes and I need to generate 12 Pfacets, so every node would be reused if I clump pfacet+nodes into a clump for each of the walls.

I would expect 4 nodes / 2 facets for a rectangle, but whatever the numbers you can clump them all together in one body, hence the nodes won't be repeated.

Bruno

 Hanying Zhang (xxxe) said on 2021-01-28: #6

>Sure of that? Any example script?

Just replace PFacet code with:
walls=aabbWalls([mi,ma],thickness=0,material='walls')
wallIds=O.bodies.append(walls)
(try script below)

>I would expect 4 nodes / 2 facets for a rectangle, but whatever the numbers you can clump them all together in one body, hence the nodes won't be repeated.

4 nodes / 2 facets for a rectangle.
a box has 8 vertices and 6 sides(12 Pfacets).
hence 8 nodes would be repeated?

Can I just add force on Pfacets?

Thanks!
#################################
#run and you can see cylinders go through wall like there is not a wall

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

#parameters

#mi,ma = (-100e-3,-100e-3,-200e-3),(100e-3,100e-3,200e-3)
mi,ma = (-10e-3,-10e-3,-20e-3),(10e-3,10e-3,20e-3)

color=[0.,1.,1.]

#### Parameter ####
young=4.0e6
poisson=3
density=1e3
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,20),damping=0.5,label='newton'),
]

triax.goal1=triax.goal2=triax.goal3=-0.9
walls=aabbWalls([mi,ma],thickness=0,material='walls')
wallIds=O.bodies.append(walls)
#=======================================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=================================================
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(36)
####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(-20e-3, 20e-3 )
random.seed( )
sita0=random.uniform(-np.pi,np.pi)
random.seed( n+7)
l_y=random.uniform(-10e-3,10e-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)<10e-3 and 0<abs(y1)<10e-3 and 0<abs(z1)<20e-3 and 0<abs(x0)<10e-3 and 0<abs(y0)<10e-3 and 0<abs(z0)<20e-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()

 Hanying Zhang (xxxe) said on 2021-01-28: #7

I got problems when clumping nodes and pfacets into a clump. There are several ways to creat Pfacets and I tried two of them, both did't work. Here is one way I tried. Please help me with it. Thanks!

####################################
nodes=[]
#up
aa= nodes.append(O.bodies.append( gridNode([-110e-3, 110e-3,-200e-3],RNode,wire=True,fixed=False,material='gridNodeMat') ))
bb=nodes.append(O.bodies.append( gridNode([ 110e-3, 110e-3,-200e-3],RNode,wire=True,fixed=False,material='gridNodeMat') ))
cc= nodes.append(O.bodies.append( gridNode([ 110e-3,-110e-3,-200e-3],RNode,wire=True,fixed=False,material='gridNodeMat') ))
dd= nodes.append(O.bodies.append( gridNode([-110e-3,-110e-3,-200e-3],RNode,wire=True,fixed=False,material='gridNodeMat') ))

ab= O.bodies.append(gridConnection(aa,bb,RNode,wire=True,material='gridNodeMat'))
ac= O.bodies.append(gridConnection(aa,cc,RNode,wire=True,material='gridNodeMat'))
bc= O.bodies.append(gridConnection(bb,cc,RNode,wire=True,material='gridNodeMat'))
bd= O.bodies.append(gridConnection(bb,dd,RNode,wire=True,material='gridNodeMat'))

abc= O.bodies.append(pfacet(aa,bb,cc,wire=True,material='pFacetMat',color=color))
acd= O.bodies.append(pfacet(aa,cc,dd,wire=True,material='pFacetMat',color=color))

clump1=O.bodies.appendClumped([\
aa, bb, cc, dd, ab, ac, ad, bc, bd, abc, acd
])

###################################################
Thanks!

 Hanying Zhang (xxxe) said on 2021-01-28: #8

And the error says:

File "1-28.py", line 124, in <module>
ab= O.bodies.append(gridConnection(aa,bb,RNode,wire=True,material='gridNodeMat'))
sph1=O.bodies[id1] ; sph2=O.bodies[id2]
ArgumentError: Python argument types in
BodyContainer.__getitem__(BodyContainer, NoneType)
did not match C++ signature:
__getitem__(pyBodyContainer {lvalue}, int)

Thanks!

 Bruno Chareyre (bruno-chareyre) said on 2021-01-28: #9

That's too many questions for me, sorry.
Simple answer: you did not show a problem with Box's, you showed a problem with Wall's (aabbWalls())
Walls are not Boxes. :)
This being said I would have expected this to work with walls too.

Could you please show the pfacet-box or pfacet-wall problem in a MWE [1]?

Bruno

 xxxea (xxxea) said on 2021-01-28: #10

 Jérôme Duriez (jduriez) said on 2021-02-02: #11

In #7 and #8 I think you're kind of misusing the functions gridNode and gridConnection (which I do not know at all :-) )

What the error says is that you gave something else than an integer (a Body id typically) at some point where only an integer (a Body id typically) was expected.

Maybe your aa, bb, .. should be these integers (as it would be with something like aa = O.bodies.append(..))

Instead of the (None) return type of Python list.append function (as it is now with aa = nodes.append(O.bodies.append(..))

 Hanying Zhang (xxxe) said on 2021-02-10: #12

Thanks Jérôme Duriez, that solved my question.