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 yade.gridpfacet import *
from yade import pack, plot
#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,
 stressMask = 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 = yade.pack.SpherePack()
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:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Jérôme Duriez
Solved:
2021-02-10
Last query:
2021-02-10
Last reply:
2021-02-02
Hanying Zhang (xxxe) said : #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!!

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 : #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?
Thanks for your help!

> -----原始邮件-----
> 发件人: "Bruno Chareyre" <email address hidden>
> 发送时间: 2021-01-27 02:41:01 (星期三)
> 收件人: <email address hidden>
> 抄送:
> 主题: Re: [Question #695192]: wall.append(O.bodies.append(pfacet) in triaxial test
>
> Your question #695192 on Yade changed:
> https://answers.launchpad.net/yade/+question/695192
>
> Status: Open => Answered
>
> 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
>
> --
> If this answers your question, please go to the following page to let us
> know that it is solved:
> https://answers.launchpad.net/yade/+question/695192/+confirm?answer_id=1
>
> If you still need help, you can reply to this email or go to the
> following page to enter your feedback:
> https://answers.launchpad.net/yade/+question/695192
>
> You received this question notification because you asked the question.

Hanying Zhang (xxxe) said : #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.

Thanks for your help!

> 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 : #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 yade.gridpfacet import *
from yade import pack, plot
#from random import random as r
import numpy as np
from numpy import *
import math
from yade import utils, qt
import random
from yade import polyhedra_utils, qt

#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,

 stressMask = 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 = yade.pack.SpherePack()
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 : #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'))
ad= O.bodies.append(gridConnection(aa,dd,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 : #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'))
  File "/usr/lib/x86_64-linux-gnu/yade/py/yade/gridpfacet.py", line 109, in gridConnection
    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!

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

[1] https://www.yade-dem.org/wiki/Howtoask

xxxea (xxxea) said : #10

Thanks for your patience! I have read the article about Pfacet and I think it is what I want for wall as to act as soft membrane.And I tried to clump nodes and pfacs as you suggest but it didn't work. Could you help me with it?(I have commented on this question earlier) Thanks! zhanghanying 邮箱:<email address hidden> 签名由 网易邮箱大师 定制 On 01/28/2021 23:55, Bruno Chareyre wrote: Your question #695192 on Yade changed: https://answers.launchpad.net/yade/+question/695192    Status: Open => Needs information Bruno Chareyre requested more information: 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 [1] https://www.yade-dem.org/wiki/Howtoask -- To answer this request for more information, you can either reply to this email or enter your reply at the following page: https://answers.launchpad.net/yade/+question/695192 You received this question notification because you asked the question.

Best Jérôme Duriez (jduriez) said : #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 : #12

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