Cross-section image

Asked by Mithushan Soundaranathan on 2021-05-05

Hi,

How can I take a 2D cross -section image of my compacts. Similar to a X-ray computed microtomography (XμCT).

Best,
MIthu

Here is my code:
#!/usr/bin/env python
#encoding: ascii

# Testing of the Deformation Enginge with Luding Contact Law
# Modified Oedometric Test
# The reference paper [Haustein2017]

from __future__ import print_function
from yade import utils, plot, timing
from yade import pack
import pandas as pd

o = Omega()

# Physical parameters
fr = 0.54
rho = 1050
Diameter = 0.000079
D=Diameter
r1 = Diameter
r2 = Diameter
k1 = 126000
kp = 12.0*k1
kc = k1 * 0.1
ks = k1 * 0.1
DeltaPMax = Diameter/3.0
Chi1 = 0.34

o.dt = 1.0e-10
particleMass = 4.0/3.0*math.pi*r1*r1*r1*rho

Vi1 = math.sqrt(k1/particleMass)*DeltaPMax*Chi1
PhiF1=0.999
#PhiF1 = DeltaPMax*(kp-k1)*(r1+r2)/(kp*2*r1*r2)

Tab_rad=0.0017
Tab_height=0.0075
Comp_press=3e8

cross_area=math.pi*(Tab_rad**2)
Comp_force=Comp_press*cross_area

#*************************************

# Add material
mat1 = O.materials.append(LudingMat(frictionAngle=fr, density=rho, k1=k1, kp=kp, ks=ks, kc=kc, PhiF=PhiF1, G0 = 0.0))

# Spheres for compression

sp=pack.SpherePack()
sp.makeCloud((-15.0*Diameter,-15.0*Diameter,-40*Diameter),(15.0*Diameter,15.0*Diameter,40.0*Diameter), rMean=Diameter/2.0,rRelFuzz=.18)
#cyl = pack.inCylinder((0,0,0),(0,0,40*Diameter),40*Diameter-0.006)
#sp = pack.filterSpherePack(cyl,sp,True,material=mat1)
sp.toSimulation(material=mat1)

##No of particles
count = 0
for b in O.bodies:
    if isinstance(b.shape, Sphere):
        count +=1
##Mass of particles
m_p=count

######################################################################
#O.bodies.append(geom.facetBox((0,0,0), (4.0*Diameter,4.0*Diameter,4.0*Diameter), wallMask=63-32, material=mat1))
walls=O.bodies.append(yade.geom.facetCylinder((0,0,0),radius=Tab_rad,height=Tab_height,segmentsNumber=20,wallMask=6,material=mat1))

df = pd.DataFrame(columns=['time','porosity_mid_d','porosity_mid_d3'])
df.to_csv('PH101_rev.csv')
from csv import writer
def append_list_as_row(file_name, list_of_elem):
    # Open file in append mode
    with open(file_name, 'a+', newline='') as write_obj:
        # Create a writer object from csv module
        csv_writer = writer(write_obj)
        # Add contents of list as last row in the csv file
        csv_writer.writerow(list_of_elem)

df = pd.DataFrame(columns=['Start_point','End_point','Volume','Porosity'])
df.to_csv('PH101_rev.csv')
from csv import writer
def append_list_as_row(file_name, list_of_elem):
    # Open file in append mode
    with open(file_name, 'a+', newline='') as write_obj:
        # Create a writer object from csv module
        csv_writer = writer(write_obj)
        # Add contents of list as last row in the csv file
        csv_writer.writerow(list_of_elem)

# Add engines
o.engines = [
  ForceResetter(),
  InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1.05),
                         Bo1_Wall_Aabb(),
                         Bo1_Facet_Aabb()
                         ]),
  InteractionLoop(
    [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=1.05),
     Ig2_Facet_Sphere_ScGeom(),
     Ig2_Wall_Sphere_ScGeom()],
    [Ip2_LudingMat_LudingMat_LudingPhys()],
    [Law2_ScGeom_LudingPhys_Basic()]
  ),
  NewtonIntegrator(damping=0.1, gravity=[0, 0, -9.81]),
  #VTKRecorder(fileName='vtk-',recorders=['all'],iterPeriod=10000),
  PyRunner(command='checkForce()', realPeriod=1, label="fCheck"),
  DeformControl(label="DefControl")
]

def checkForce():
    # at the very start, unbalanced force can be low as there is only few
    # contacts, but it does not mean the packing is stable
    if O.iter < 1000:
        return
    # the rest will be run only if unbalanced is < .1 (stabilized packing)
    timing.reset()
    #if unbalancedForce() > 0.2:
    # return
    # add plate at upper box side

    highSphere = 0.0
    for b in O.bodies:
        if highSphere < b.state.pos[2] and isinstance(b.shape, Sphere):
            highSphere = b.state.pos[2]
        else:
            pass

    O.bodies.append(wall(highSphere+0.5*Diameter, axis=2, sense=-1, material=mat1))
    # without this line, the plate variable would only exist inside this
    # function
    global plate
    plate = O.bodies[-1] # the last particles is the plate
    # Wall objects are "fixed" by default, i.e. not subject to forces
    # prescribing a velocity will therefore make it move at constant velocity
    # (downwards)
    plate.state.vel = (0, 0, -30)
    # start plotting the data now, it was not interesting before
    O.engines = O.engines + [PyRunner(command='addPlotData()', iterPeriod=1000)]
    # next time, do not call this function anymore, but the next one
    # (unloadPlate) instead
    fCheck.command = 'unloadPlate()'

def unloadPlate():
    # if the force on plate exceeds maximum load, start unloading
    # if abs(O.forces.f(plate.id)[2]) > 5e-2:
    if abs(O.forces.f(plate.id)[2]) > Comp_force:
        plate.state.vel *= -1
        # next time, do not call this function anymore, but the next one
        # (stopUnloading) instead
        fCheck.command = 'stopUnloading()'

def stopUnloading():
    if abs(O.forces.f(plate.id)[2]) < Comp_force:
        # O.tags can be used to retrieve unique identifiers of the simulation
        # if running in batch, subsequent simulation would overwrite each other's output files otherwise
        # d (or description) is simulation description (composed of parameter values)
        # while the id is composed of time and process number
        # plot.saveDataTxt(O.tags['d.id'] + '.txt')
        plot.saveDataTxt('data'+ O.tags['id'] +'.txt')
        print(timing.stats())
        O.pause()
        #unsat=TwoPhaseFlowEngine()
        #unsat.initialization()
        #por=unsat.getCellPorosity(10)

def addPlotData():
    if not isinstance(O.bodies[-1].shape, Wall):
        plot.addData()
        return
    Fz = O.forces.f(plate.id)[2]
    plot.addData(
        Fz=Fz,
        w=plate.state.pos[2] - (-4*Diameter),
        unbalanced=unbalancedForce(),
        i=O.iter
    )

def defVisualizer():
   with open("data.dat","a") as f:
       for b in O.bodies:
           if isinstance(b.shape, Sphere):
               rData = "{x},{y},{z},{r},{w}\t".format(x = b.state.pos[0],
                                                y = b.state.pos[1],
                                                z = b.state.pos[2],
                                                r = b.shape.radius + b.state.dR,
                                                w = plate.state.pos[2]
                                               )
               f.write(rData)
       f.write("\n")

O.timingEnabled=True
O.run(1, True)
plot.plots={'w':('Fz', None)}
plot.plot()

Question information

Language:
English Edit question
Status:
Answered
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
Last reply:
Revision history for this message
Jan Stránský (honzik) said :
#1

Hello,

please be more specific. Like what shapes will you like for the image, what should be the output (just raster image, some more details about the section, 3D image, just image or some data format etc. etc.)

Also, the code should be a MWE [1], M = minimal, like a few predefined bodies is OK to illustrate the problem.
Materials, running, checkForce, unloadPlate, stopUnloading, addPlotData .... are all irrelevant tot he problem.

Have a look at [2] if it is what you are looking for.

Cheers
Jan

[1] https://www.yade-dem.org/wiki/Howtoask
[2] https://gitlab.com/yade-dev/trunk/-/tree/master/examples/test/paraview-spheres-solid-section

Revision history for this message
Mithushan Soundaranathan (mithushan93) said :
#2

Hi Jan,

I want 2D raster images of the compact in a rectangle shape at different places. Than stack them on each other to create the 3D structure of the compacts. Hope that was more clear.

Best,
Mithu

Revision history for this message
Karol Brzezinski (kbrzezinski) said :
#3

Hi,

please see this topic [1]. I believe that similar issue was raised here more than once, so you can also use search option with keyword "cross section" etc.

Best wishes,
Karol

[1] https://answers.launchpad.net/yade/+question/674789

Revision history for this message
Karol Brzezinski (kbrzezinski) said :
#4

Hi,

I have just realized that the example topic that I provided contains broken link. Sorry for that. There are some potential solutions in other topics [1], but I will try to make a script today or tomorrow if this topic will be still open.

Karol

[1] https://answers.launchpad.net/yade/+question/240738

Revision history for this message
Karol Brzezinski (kbrzezinski) said :
#5

Hi,

please find my script below. It takes about 7 minutes on my computer to take 10 slices with 0.6 Mpx resolution (pixel size 0.000005 m), but obviously you need more slices for your reconstruction.

Cheers,
KB

#######################
import numpy as np
from PIL import Image

n_slices = 50

x_s = np.array([])
y_s = np.array([])
z_s = np.array([])
radii = np.array([])
for b in O.bodies:
 if isinstance(b.shape, Sphere):
  pos = b.state.pos
  x_s = np.append(x_s,pos[0])
  y_s = np.append(y_s,pos[1])
  z_s = np.append(z_s,pos[2])
  radii = np.append(radii,b.shape.radius)

px_size = 0.000005 # pixel size (m)
(x_min,y_min,z_min),(x_max,y_max,z_max) = aabbExtrema()
real_h = z_max-z_min
real_w = x_max - x_min # I assume here that slices will be perpendicular to y-axis

img_h = ceil(real_h/px_size) # image height in px
img_w = ceil(real_w/px_size) # image height in px

"""
The idea is as follows:
- divide the cross section into squares (representing circles)
- for the center of every square (pixel) check whether any center of the sphere is closer than its radius
- if True color the pixel white, otherwise black

One could probabily optimize this. I only took adance on the vectorisation in numpy.
I have made a "manual" decision for this algorithm: since the number of pixels is greater than number of bodies I loop over spheres and do vectorized operation on the "pixel" array.

"""

x_px = np.tile(np.arange(x_min+px_size/2,x_min+px_size*img_w,px_size),img_h)
z_px = np.repeat(np.arange(z_min+px_size/2,z_min+px_size*img_h,px_size),img_w)

for y_slice in np.linspace(y_min,y_max,n_slices+2)[1:-1]:# take n_slices slices (linspace prepares n_slices+2, but I ommit first and last since it is not interesting)
 y_px = np.repeat(y_slice,img_h*img_w)#
 pixels_state = np.zeros(img_h*img_w).astype(bool)

 count = 0
 for i in range(len(radii)):
  radius = radii[i]
  x = x_s[i]
  y = y_s[i]
  z = z_s[i]
  if abs(y-y_slice)<=radius:# I make use of simplified assumption, I only check pixels for the spheres that are close enough to my plane.
   pixels_state_tmp = np.sqrt((x-x_px)**2+(y-y_px)**2+(z-z_px)**2)<radius
   pixels_state = pixels_state + pixels_state_tmp

 img_array = 256*pixels_state.reshape(img_h,img_w).astype(np.uint8)

 img = Image.fromarray(img_array)
 img.save('slice_y_{:.5f}.gif'.format(y_slice))

Revision history for this message
Mithushan Soundaranathan (mithushan93) said :
#6

Hi Karol,

Thank you very much, it very helpful.
Just a quick question. I am taking cross section of compacted spheres, will it account for that. Since you defined in your code if isinstance(b.shape, Sphere).

Best,
Mithu

Revision history for this message
Karol Brzezinski (kbrzezinski) said :
#7

Hi,

you have only spheres in simulation. After compaction, the are going to be compacted (overlaping) spheres - nothing changes in my opinion.

Karol

PS Maybe you should change the last line if you will "slice denser". Save files using six or seven decimal places

img.save('slice_y_{:.7f}.gif'.format(y_slice))

Revision history for this message
Mithushan Soundaranathan (mithushan93) said :
#8

Hi Karol,

Thank you very much, sorry was not expressing myself clearly. I meant to say, will this method capture the deformed particles?

Best,
Mithu

Revision history for this message
Karol Brzezinski (kbrzezinski) said :
#10

Hi,

my bad. I totally forgot that I had a problem with the Deformation Engine and I commented it out (I assumed it is some sort of servomechanism like Triaxial Engine).

I definitely didn't account for deformation. After a skim-reading the Haustein [1] paper, I think that you can easily account for that. Instead of the initial particle radii, use "extended radii". (If I understand correctly, the change in the shape is actually increasing the radius).

Best wishes,
Karol

[1] http://dx.doi.org/https://doi.org/10.1016/j.softx.2017.05.001

Revision history for this message
Mithushan Soundaranathan (mithushan93) said :
#11

Hi Karol,

The deformation was not visible in the images, thank you for your help.

Best regards,
Mithu

Revision history for this message
Karol Brzezinski (kbrzezinski) said :
#12

Hi,

thank you for the update. Am I correct that the deformed particles in this simulation are still spherical, just bigger? Is it visible during a preview of the simulation? Can you access the information about increased radii? If so, the previous script should solve the issue (if you use radii of deformed particles instead of undeformed).

Best wishes,
Karol

Can you help with this problem?

Provide an answer of your own, or ask Mithushan Soundaranathan for more information if necessary.

To post a message you must log in.