Rock Joints Generation & TriPatchSet()

the Following is my script:

from gengeo import *
#from read_jointset import *
import sys

# An example python script to generate a bonded rectangular prism with Joints

# particle properties
minRadius = 0.2
maxRadius = 1.0

# Define region extremities using the bounding box
minPoint = Vector3(-10.0,-10.0,-10.0)
maxPoint = Vector3(10.0,10.0,10.0)

# Define the volume to be filled with spheres:

box = BoxWithPlanes3D (
   minPoint = minPoint,
   maxPoint = maxPoint
)

# boundary planes
box.addPlane(Plane(minPoint,Vector3(1.0,0.0,0.0)))
box.addPlane(Plane(minPoint, Vector3(0.0,1.0,0.0)))
box.addPlane(Plane(minPoint, Vector3(0.0,0.0,1.0)))
box.addPlane(Plane(maxPoint, Vector3(-1.0,0.0,0.0)))
box.addPlane(Plane(maxPoint, Vector3(0.0,-1.0,0.0)))
box.addPlane(Plane(maxPoint, Vector3(0.0,0.0,-1.0)))

# Create a multi-group neighbour table to contain the particles:
mntable = MNTable3D (
   minPoint = minPoint,
   maxPoint = maxPoint,
   gridSize = 2.5*maxRadius,
   numGroups = 1
)

# Fill the volume with particles:
packer = InsertGenerator3D (
   minRadius = minRadius,
   maxRadius = maxRadius,
   insertFails = 1000,
   maxIterations = 1000,
   tolerance = 1.0e-6,
   seed=42
)

packer.generatePacking(
   volume = box,
   ntable = mntable,
   groupID = 0,
   tag = 2
)

#set up joint set
joints=TriPatchSet()
joints.addTriangle(
    Point1=Vector3(-8.63342787782,-5.68558653544,11.6066017178),
    Point2=Vector3(10.6334278778,7.68558653544,-9.6066017178),
    Point3=Vector3(1.97317383998,12.6855865354,-9.6066017178),
    tag=12
)
joints.addTriangle(
    Point1=Vector3(0.0268261600231,-10.6855865354,11.6066017178),
    Point2=Vector3(10.6334278778,7.68558653544,-9.6066017178),
    Point3=Vector3(-8.63342787782,-5.68558653544,11.6066017178),
    tag=12
)
print "Bounding Box of Joint Set: " , joints.getMinPoint()

#bonding
tol=1.0e-2
basetag=1 # tag for the "normal" bonds
mntable.generateBondsWithJointSet(joints,0,tol,basetag)

# write a geometry and a VTK file
mntable.write("box_"+str(maxRadius)+".geo",1)
mntable.write("box_"+str(maxRadius)+".vtu",2)

Hi Henry,

When generating a geometry containing a joint set (or TrianglePatchSet), one must use special Volume classes. In your case, change BoxWithPlanes3D to BoxWithJointSet. You will also need to add the jointset to the Volume before generating the packing. For your script, add the following command after the multiple box.addPlane(..) commands.

box.addJoints(joints)

To tag particles near the joints with a different tag to the rest of the particles, add the following command just before writing the geometry file:

mntable.tagParticlesAlongJoints (
   jointset = joints,
   distance = 1.1*maxRadius,
   tag = 3,
   mask = -1,
   groupID = 0
)

Cheers,

Dion.

Dear Dion,
   Thanks for your help. Acrossing your suggestion I modified my script.
   but only the partiles along the line Point1 & Point2 of the joints are taged.

  to simplify the joints and check the results, I changed the jointset.
  and following is my modified script . Thanks a lot!

Henry

from gengeo import *
#from read_jointset import *
import sys

# An example python script to generate a bonded rectangular prism with Joints

# particle properties
minRadius = 0.5
maxRadius = 2.0

# Define region extremities using the bounding box
minPoint = Vector3(-10.0,-10.0,-10.0)
maxPoint = Vector3(10.0,10.0,10.0)

#set up joint set
joints=TriPatchSet()

# Define the volume to be filled with spheres:

box = BoxWithJointSet (
   minPoint = minPoint,
   maxPoint = maxPoint
)

# boundary planes
box.addPlane(Plane(minPoint,Vector3(1.0,0.0,0.0)))
box.addPlane(Plane(minPoint, Vector3(0.0,1.0,0.0)))
box.addPlane(Plane(minPoint, Vector3(0.0,0.0,1.0)))
box.addPlane(Plane(maxPoint, Vector3(-1.0,0.0,0.0)))
box.addPlane(Plane(maxPoint, Vector3(0.0,-1.0,0.0)))
box.addPlane(Plane(maxPoint, Vector3(0.0,0.0,-1.0)))

box.addJoints(joints)

# Create a multi-group neighbour table to contain the particles:
mntable = MNTable3D (
   minPoint = minPoint,
   maxPoint = maxPoint,
   gridSize = 2.5*maxRadius,
   numGroups = 1
)

# Fill the volume with particles:
packer = InsertGenerator3D (
   minRadius = minRadius,
   maxRadius = maxRadius,
   insertFails = 1000,
   maxIterations = 1000,
   tolerance = 1.0e-6,
   seed=42
)

packer.generatePacking(
   volume = box,
   ntable = mntable,
   groupID = 0,
   tag = 2
)
joints.addTriangle(
    Point1=Vector3(0,-10,10),
    Point2=Vector3(0,10,-10),
    Point3=Vector3(0,-10,-10),
    tag=4
)

print "Bounding Box of Joint Set: " , joints.getMinPoint()

#bonding
tol=1.0e-2
basetag=1 # tag for the "normal" bonds
mntable.generateBondsWithJointSet(joints,0,tol,basetag)

#o tag particles near the joints with a different tag to the rest of the particles, add the following command just before writing the geometry file:
mntable.tagParticlesAlongJoints (
   jointset = joints,
   distance = 1.1*maxRadius,
   tag = 3,
   mask = -1,
   groupID = 0
)

# write a geometry and a VTK file
mntable.write("box_"+str(maxRadius)+".geo",1)
mntable.write("box_"+str(maxRadius)+".vtu",2)

Hi Henry,

I can reproduce the bug using your script or variants of it but I cannot figure out what is wrong. If I use my own script (see below) it works correctly, despite being essentially the same as yours. I will have to ask you to try to work out what is wrong in your script and get back to us.

Cheers,

Dion.

from gengeo import *

def makeGeometry(numTrys=1000):
 #An example python script to generate a bonded rectangular block with 2 joints
 xdim=20.0
 ydim=xdim
 zdim=xdim

 # Define region extremities:
 maxRadius = 2.0
 minRadius = 0.5
 minPoint = Vector3(0.0,0.0,0.0)
 maxPoint = Vector3(xdim,ydim,zdim)

 # Define the volume to be filled with spheres:
 # (e.g. a box bounded by planes)
 box = BoxWithJointSet (minPoint,maxPoint)

 box.addPlane(Plane(minPoint,Vector3(1.0,0.0,0.0)))
 box.addPlane(Plane(minPoint, Vector3(0.0,1.0,0.0)))
 box.addPlane(Plane(minPoint, Vector3(0.0,0.0,1.0)))
 box.addPlane(Plane(maxPoint, Vector3(-1.0,0.0,0.0)))
 box.addPlane(Plane(maxPoint, Vector3(0.0,-1.0,0.0)))
 box.addPlane(Plane(maxPoint, Vector3(0.0,0.0,-1.0)))

 # define empty joint set
 Joints=TriPatchSet()

        # add triangles to joint set to define location of pre-existing crack
 Joints.addTriangle(
    Vector3(0.0,ydim/2.0,0.0),
    Vector3(0.0,ydim/2.0,zdim),
    Vector3(xdim/2.0,ydim/2.0,zdim),
    1
 )
 Joints.addTriangle(
    Vector3(0.0,ydim/2.0,0.0),
    Vector3(xdim/2.0,ydim/2.0,0.0),
    Vector3(xdim/2.0,ydim/2.0,zdim),
    1
 )

 # add joints to box
 box.addJoints(Joints)

 # Create a multi-group neighbour table to contain the particles:
 mntable = MNTable3D (
    minPoint = minPoint,
    maxPoint = maxPoint,
    gridSize = 2.5*maxRadius,
    numGroups = 1
 )

 # Fill the volume with particles:
 packer = InsertGenerator3D (
    minRadius = minRadius,
    maxRadius = maxRadius,
    insertFails = numTrys,
    maxIterations = 1000,
    tolerance = 1.0e-6
 )

 packer.generatePacking(
    volume = box,
    ntable = mntable,
    groupID = 0
 )

 # create bonds between neighbouring particles:
 #mntable.generateBonds(
 # groupID = 0,
 # tolerance = 1.0e-5,
 # bondID = 0
 #)
 mntable.generateBondsWithJointSet(Joints,0,1.0e-5,0)

 #tag particles along joints
 mntable.tagParticlesAlongJoints(Joints,0.01,2,2,0)

# mntable.tagParticlesAlongPlane(Plane(Vector3(0,0,0),Vector3(0,1,0)), 2.5*maxRadius, 3, 0)
# mntable.tagParticlesAlongPlane(Plane(Vector3(0,ydim,0),Vector3(0,-1,0)), 2.5*maxRadius, 4, 0)

 # print the porosity:
 volume = xdim*ydim*zdim
 porosity = (volume - mntable.getSumVolume(groupID=0))/volume
 print "Porosity: ", porosity

 # write a geometry file
 mntable.write('geometry.geo', 1)
 mntable.write('geometry.vtu', 2)

if __name__=="__main__":
   makeGeometry()

Dear Dion,
     Thanks a lot.
     I find the errors in my script:
   the addTriangle should be used as:

    joints.addTriangle(
    Vector3(0,-10,10),
    Vector3(0,10,-10),
    Vector3(0,-10,-10),
    4
  )

not :
joints.addTriangle(
    Point1=Vector3(0,-10,10),
    Point2=Vector3(0,10,-10),
    Point3=Vector3(0,-10,-10),
    tag=4
)

Thanks!

Best Regard,
Henry