Question #695192 on Yade changed:
https://answers.launchpad.net/yade/+question/695192
Status: Answered => Open
Hanying Zhang is still having a problem:
>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()
--
You received this question notification because your team yade-users is
an answer contact for Yade.
_______________________________________________
Mailing list: https://launchpad.net/~yade-users
Post to : [email protected]
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
