Divide numerical periodic microstructure box into smaller grided cubes or boxes

Hello,

> I am new to Yade

welcome :-)

> I want to divide this sample into cubes and play with each cube by removing certain amount of particle sizes.

please be (much) more specific.
There are too many ambiguous terms in your request.

Ideally provide an example:
- simple version of the packing after compaction (2D would be sufficient for illustration?)
- what the process and result of "divide this sample into cubes" should look like
- what "removing certain amount of particle sizes" should look like
- ...

Cheers
Jan

Hello Jan,

Thank you for your quick reply,

My script used for generating my sample is as follows:

for test in range(1, num_tests):
    O.reset()
    O.periodic=True
    O.cell.hSize=Matrix3(Hbox,0,0, 0,Lbox,0, 0,0,Hbox)
    # Define the material for balls
    O.materials.append(FrictMat(young=young,poisson=poisson,frictionAngle=radians(CompFricDegree),density=dens,label='spheres'))
    #-------------Create a periodic packing of spheres
    mn,mx=Vector3(0,0,0),Vector3(Hbox,Lbox,Lbox) # corners of the initial packing
    sp=pack.SpherePack()
    sp.makeCloud(minCorner=mn,maxCorner=mx,num=num_spheres,psdSizes=D_spheres,psdCumm=CumFrac,distributeMass=True,periodic=True)
    #add sphere packing to simulation
    sp.toSimulation(material='spheres')
    #--------------- Define engines for the numerical simulations
    O.engines=[
        ForceResetter(), #0
        InsertionSortCollider([Bo1_Sphere_Aabb()]), #1
        InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom()], # contact geometry
        [Ip2_FrictMat_FrictMat_FrictPhys()], # contact properties
        [Law2_ScGeom_FrictPhys_CundallStrack()] # contact forces
        ), #2
        GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8), #3, determine step size
        NewtonIntegrator(damping=damp) #4
    ]
    newton = O.engines[4]
    # Then, we reduce the box dimensions in such a way that the stress state reach the target stress tensor
    # First stade of compaction. We compact rapidly the sample until a normal stress exceed 90% of the target confining stress

    # Then, we reduce the box dimensions in such a way that the stress state reach the target stress tensor
    # First stade of compaction. We compact rapidly the sample until a normal stress exceed 90% of the target confining stress
    print("------------ First stade of compaction")
    sm.ServoCompaction1(damp,servo_g1,ConfiningStress,0.9,UnbThreshold)
    # Second stade of compaction. We compact slowly the sample until all the three normal stresses reach 90% of the targer confining stress
    print("------------ Second stade of compaction")
    sm.ServoCompaction1(damp,servo_g2,ConfiningStress,1.0,UnbThreshold)

    newton.damping = damp # reset the numerical damping
    # Reset the friction angle
    O.materials[0].frictionAngle=radians(ShearFricDegree) # radians
    # for existing contacts, set contact friction directly
    for i in O.interactions:
        i.phys.tangensOfFrictionAngle=tan(radians(ShearFricDegree))

    # Final compaction
    sm.FinalCompaction1(damp,servo_g=servo_g2,sigTarget=ConfiningStress,SigThreshold=StressThreshold)

    print('final porosity ', utils.porosity())

    SaveFileName = TagName + str(test)+".yade.gz"
    os.chdir(SavePath)
    O.save(SaveFileName)

What I want to develop:
After the final compaction of my sample, First, I want to divide or discritize my domain (3D box size) into equal subdomains (example 100 smaller boxes) (As if a 3D grid)
Then, I want to choose half of these subdomains randomly and impose a criteria to eliminate particle sizes below a specific imposed diameter.
Thus, I will have at the end a new sample of same domain (size) but less number of particles and less porosity...
I need please help regarding this matter by guides or links or documentations.
Thank you in advance Jan, appreciated...

> I need please help regarding this matter by guides or links or documentations.

I don't think you will find this specific problem as is in the documentation or there is an existing script.
You can however build your script from "building blocks" found in the documentation.
As a start, below is a 2D example.
Please let us know if this is in principle what you are looking for.

Cheers
Jan

###
import random

tileSize = (2, 3)
tiling = (7, 6)

# dummy sample
########################################
colors = {
    0.5: (0,1,1),
    0.2: (1,0,0),
}
tileParticles = (
    ((0.0, 0.0), 0.5),
    ((1.0, 0.4), 0.5),
    ((0.5, 1.3), 0.5),
    ((1.5, 1.6), 0.5),
    ((0.9, 2.4), 0.5),
    ((0.3, 0.7), 0.2),
    ((1.7, 0.6), 0.2),
    ((1.8, 1.0), 0.2),
    ((1.6, 2.4), 0.2),
    ((0.1, 2.4), 0.2),
    ((0.0, 2.0), 0.2),
)
for ix in range(tiling[0]):
    x0 = ix * tileSize[0]
    for iy in range(tiling[1]):
        y0 = iy * tileSize[1]
        for center, radius in tileParticles:
            x = x0 + center[0]
            y = y0 + center[1]
            O.bodies.append(sphere((x,y,0), radius, color=colors[radius]))
########################################

def isInTile(tile,xy):
    x, y = xy
    sizeX, sizeY = tileSize
    x0 = tile[0] * sizeX
    y0 = tile[1] * sizeY
    x1 = x0 + sizeX
    y1 = y0 + sizeY
    return x >= x0 and x < x1 and y >= y0 and y < y1

def eliminate():
    nTiles = tiling[0] * tiling[1]
    tiles = set()
    while len(tiles) < 0.5 * nTiles:
        x = random.randint(0,tiling[0]-1)
        y = random.randint(0,tiling[1]-1)
        tiles.add((x,y))
    for b in O.bodies:
        if b.shape.radius > 0.4: # do nothing for "large" particles
            continue
        x,y,z = b.state.pos
        if any(isInTile(tile,(x,y)) for tile in tiles): # small particles in selected tiles
            O.bodies.erase(b.id)

# see with or without eliminating
#eliminate()

#yade.qt.View()
###

Hello Jan,

Thank you for your help, it was really helpful!

I have one additional question please:
I want to add a random window (of cubic or spherical shape) on my sample so that, this random window will act as a mapping device or radar domain, so that I can delete particles inside this window.
Thus, at the end I will have the original sample with square voids or pores.
To be more precise, I want a function to help me generate a random window shape (spherical shape is preferable).

Thank you again for your help!

Sincerely,

Hello,

deleting particles from certain region is not difficult.
As a starting point, see a MWE below.

cheers
Jan

###
from yade import pack
pred = pack.inAlignedBox((0, 0, 0), (30, 30, 8))
sphs = pack.regularHexa(pred, radius=1, gap=0)
O.bodies.append(sphs)

pores = [ # list of (center,radius) values of the spherical pores
    (Vector3(20, 20, 1), 7),
    (Vector3(6, 6, 6), 5),
]

def isInPore(body, pore):
    bpos, br = body.state.pos, body.shape.radius
    ppos,pr = pore
    dc = (bpos-ppos).norm() # center distance
    return dc + br < pr # including body radius (entire body is inside pore), could be only center-based or mixed

def makePores():
    for b in O.bodies:
        if any(isInPore(b,pore) for pore in pores):
            O.bodies.erase(b.id)

makePores()
###