Module DEMIC

#-------------------------------------------------------------------------------
#Librairies
#-------------------------------------------------------------------------------

from yade import pack, plot, export
import numpy as np
import matplotlib.pyplot as plt
import math, random, pickle

#-------------------------------------------------------------------------------
#User
#-------------------------------------------------------------------------------

# load data 
with open('data/dict_ic', 'rb') as handle:
    dict_ic = pickle.load(handle)

# PSD
n_grains = dict_ic['n_grain']*4
L_dr = []

# Box
Dx = 1
Dy = 1
Dz = 1

# time step
factor_dt_crit = 0.6

# steady-state detection
unbalancedForce_criteria = 0.1

#-------------------------------------------------------------------------------
#Initialisation
#-------------------------------------------------------------------------------

# define wall material (no friction)
O.materials.append(CohFrictMat(young=1, poisson=0.25, frictionAngle=0, density=2650, isCohesive=False, momentRotationLaw=False))
# create box and grains
O.bodies.append(aabbWalls([Vector3(0,0,0),Vector3(Dx,Dy,Dz)], thickness=0., oversizeFactor=1))
# a list of 6 boxes Bodies enclosing the packing, in the order minX, maxX, minY, maxY, minZ, maxZ

# define grain material
O.materials.append(CohFrictMat(young=1, poisson=0.25, frictionAngle=atan(0.05), density=2650,\
                               isCohesive=False, momentRotationLaw=False))
# generate grain
for i in range(n_grains):  
    # define kinetics
    dradius = random.uniform(math.sqrt((0.5*Dx*Dy)/(n_grains*math.pi))/400, math.sqrt((0.5*Dx*Dy)/(n_grains*math.pi))/100)
    center_x = random.uniform(0+dradius*30, Dx-dradius*30)
    center_y = random.uniform(0+dradius*30, Dy-dradius*30)
    center_z = 0
    O.bodies.append(sphere(center=[center_x, center_y, center_z], radius=dradius*30))
    O.bodies[-1].state.blockedDOFs = 'zXY'
    L_dr.append(dradius)

# yade algorithm
O.engines = [
        ForceResetter(),
        # sphere, wall
        InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Box_Aabb()]),
        InteractionLoop(
                # need to handle sphere+sphere and sphere+wall
                # Ig : compute contact point. Ig2_Sphere (3DOF) or Ig2_Sphere6D (6DOF)
                # Ip : compute parameters needed
                # Law : compute contact law with parameters from Ip
                [Ig2_Sphere_Sphere_ScGeom6D(), Ig2_Box_Sphere_ScGeom6D()],
                [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
                [Law2_ScGeom6D_CohFrictPhys_CohesionMoment(always_use_moment_law=True)]
        ),
        NewtonIntegrator(gravity=(0, 0, 0), damping=0.1, label = 'Newton'),
        PyRunner(command='checkUnbalanced_ir_ic()', iterPeriod = 200, label='checker')
]
# time step
O.dt = factor_dt_crit * PWaveTimeStep()

#-------------------------------------------------------------------------------

def compute_ratio_overlap_radius():
    '''
    Compute the mean ratio of the overlap and radius.
    '''
    overlap_average = 0
    radius_average = 0
    n_average = 0
    # iterate on interactions
    for i in O.interactions:
        # only grain-grain contact can be cemented
        if isinstance(O.bodies[i.id1].shape, Sphere) and isinstance(O.bodies[i.id2].shape, Sphere) :
            # compute overlap of this interaction
            b1_x = O.bodies[i.id1].state.pos[0]
            b1_y = O.bodies[i.id1].state.pos[1]
            b2_x = O.bodies[i.id2].state.pos[0]
            b2_y = O.bodies[i.id2].state.pos[1]
            dist = math.sqrt((b1_x-b2_x)**2+(b1_y-b2_y)**2)
            overlap = O.bodies[i.id1].shape.radius + O.bodies[i.id2].shape.radius - dist
            # compute mean value
            overlap_average = overlap_average + overlap
            radius_average = radius_average + (O.bodies[i.id1].shape.radius+O.bodies[i.id2].shape.radius)/2
            n_average = n_average + 1
    if n_average > 0:
        # plot and return the average
        overlap_average = overlap_average/n_average
        radius_average = radius_average/n_average
        return overlap_average/radius_average
    else :
        return 1
    
#-------------------------------------------------------------------------------

def checkUnbalanced_ir_ic():
    '''
    Increase particle radius until a steady-state is found.
    '''
    # the rest will be run only if unbalanced is < .1 (stabilized packing)
    # Compute the ratio of mean summary force on bodies and mean force magnitude on interactions.
    if (2<=avgNumInteractions()) and (avgNumInteractions()<=4) and (unbalancedForce() > unbalancedForce_criteria) :
        return

    if (avgNumInteractions()<4):   
        # increase the radius of particles
        i_L_r = 0
        for b in O.bodies :
            if isinstance(b.shape, Sphere):
                growParticle(b.id, (b.shape.radius+L_dr[i_L_r])/b.shape.radius)
                i_L_r = i_L_r + 1
        # update the dt as the radii change
        O.dt = factor_dt_crit * PWaveTimeStep()
        return
        
    # plot the psd
    binsSizes, binsProc, binsSumCum = psd(bins=10)
    plotPSD(binsSizes, binsProc)
    
    # close the simulation
    stopLoad()
    
#-------------------------------------------------------------------------------

def saveData():
    """
    Save data in .txt file during the ic.
    """
    # pp data
    L_pos = []
    L_radius = []
    for b in O.bodies :
        if isinstance(b.shape, Sphere):
            L_pos.append([b.state.pos[0], b.state.pos[1]])
            L_radius.append(b.shape.radius)

    # create dict
    dict_ic = {'L_pos': L_pos, 'L_radius': L_radius}

    # save
    with open('data/dict_ic', 'wb') as handle:
        pickle.dump(dict_ic, handle, protocol=pickle.HIGHEST_PROTOCOL)

#-------------------------------------------------------------------------------

def stopLoad():
    """
    Close simulation.
    """
    # save at the converged iteration
    saveData()
    # close yade
    O.pause()

    # give order to the user
    print("\nyou can type 'quit()' in the terminal")

#-------------------------------------------------------------------------------

def plotPSD(binsSizes, binsProc):
    """
    This function can be called to plot the evolution of the psd.
    """
    plt.figure(1, figsize=(16,9))
    plt.plot(binsSizes, binsProc)
    plt.title('Particle Size Distribution')
    plt.savefig('output/PSD.png')
    plt.close()

#-------------------------------------------------------------------------------
# start simulation
#-------------------------------------------------------------------------------

O.run()