Delete particles dropping by periodicity (System is blocked)
# Script to generate a new fabrication framework or a different microstructure (Heterogeneous sample)
from yade import pack, plot, utils, Matrix3
import sys
import os
sys.path.append(".")
import numpy as np
import math
import yade.qt
from yade.wrapper import *
import matplotlib.pyplot as plt
from yade import export
import time
from collections import defaultdict
import My_Functions as Isg
from scipy.special import expit # Import the sigmoid function
import matplotlib.pyplot as plt
start_time = time.time()
initial_time = O.time
damp = 0.3
g = -9.81
num_spheres = 6000
frequency_record_data = 1000 # save date every this number of steps
Exported_datafile = "0Step3_M1E7_fc20_ID50_N6000.txt"
SavePath = '/home/elshamieh/SHAMIEH edited documents/pores3d_pbc/build_pore_network/New_sample_generation_DEM_Micro_1/Damping_Parametric_study'
LoadPath = '/home/elshamieh/SHAMIEH edited documents/pores3d_pbc/build_pore_network/New_sample_generation_DEM_Micro_1/Damping_Parametric_study'
sampleName = "SAVEStep4_01E7_fc20_ID50_N6000"
AFTER_sampleName = "Step4_01E7_fc20_ID50_N6000"
# Load the numerical sample saved at the end of the isotropic compaction process
os.chdir(LoadPath)
O.load(sampleName + ".yade.gz")
os.chdir(SavePath)
t0 = O.time # initial time
#Set the cohesion in the engine to zero
material = O.materials[0]
material.normalCohesion = 0.0
material.shearCohesion = 0.0
O.engines = [
ForceResetter(), # 0
InsertionSortCollider([Bo1_Sphere_Aabb()]), # 1
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom6D()], # contact geometry
[Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()], # contact properties
[Law2_ScGeom6D_CohFrictPhys_CohesionMoment()] # contact forces
), # 2
GlobalStiffnessTimeStepper(active=1, timeStepUpdateInterval=100, timestepSafetyCoefficient=0.8), # 3, determine step size
NewtonIntegrator(gravity=(0, g, 0), damping=damp), # 4
PyRunner(command='addPlotData()', iterPeriod=frequency_record_data), # engine to add data to plot later at specified steps
]
initial_particles = set() # Set to store initial negative position particles
Box_sizes = Isg.Box_size()
Ly = Box_sizes[1]
cumulative_particles_crossed = 0 # Variable to store the cumulative result
#This function aims to count the crossing particles by periodicity and then deleting/erasing them.
def Track_particles():
global cumulative_particles_crossed # Use the global variable
fines_crossed_filter_this_step = 0
for b in O.bodies:
position_vector = b.state.pos
pos_y = position_vector[1] / Ly
# Check if pos_y is between 0 and 1
if -1 <= pos_y <= 1:
pass
elif -2 <= pos_y <= -1:
print("Particle of ID:", b.id, ",SECOND cross")
elif -3 <= pos_y <= -2:
print("Particle of ID:", b.id, ",THIRD cross")
else:
print("Particle of ID:", b.id, ",4th cross and more")
# Iterate over intervals from -1 to -5000
for threshold in range(0, -50001, -1):
if pos_y <= threshold and b.id not in initial_particles:
fines_crossed_filter_this_step += 1
# Erase the particle
O.bodies.erase(b.id)
break # Exit the loop after erasing the particle
cumulative_particles_crossed += fines_crossed_filter_this_step # Update the cumulative result
print('Number of particles crossing:', fines_crossed_filter_this_step)
print('Cumulative particles crossed:', cumulative_particles_crossed) # Print cumulative result
number_of_bodies = len(O.bodies) - cumulative_particles_crossed
print("Number of bodies:", number_of_bodies) #to show number of particles remaining
return fines_crossed_filter_this_step
# Function adds current values to the history of data, under the names specified
stress = utils.getStress()
def addPlotData():
P = Isg.linear_momentum(dgap) #P is for linear momentum track
nr = Isg.ninteractions(num_spheres) #nr is for floating particles (or no contact particles track)
T = Track_particles() #T is to track crossing particles
plot.addData(i=O.iter, t_ms=1000*(O.time - initial_time), tr=time.time() - start_time, unbalanced=utils.unbalancedForce(),
totEnergy=O.energy.total(), Px=P[0], Py=P[1], Pz=P[2], sxx=-utils.getStress()[0, 0], syy=-utils.getStress()[1, 1], szz=-utils.getStress()[2, 2], floaing_particles_percentage=nr, counting=T)
#Plotting all parameters
plot.plots = {
'tr': 'unbalanced'
}
# Show the plot on the visualizing screen
plot.plot()
plot.plots = {
'tr': (('sxx', 'xr:'), ('syy', 'Pb:'), ('szz', 'xg:'))
}
# Show the plot on the visualizing screen
plot.plot()
plot.plots = {
'tr': ('floaing_particles_percentage')
}
plot.plot()
plot.plots = {
'tr': (('Px', 'xr:'), ('Py', 'Pb:'), ('Pz', 'xg:'))
}
# Show the plot on the visualizing screen
plot.plot()
plot.plots = {
't_ms': 'Gtrack'
}
# Show the plot on the visualizing screen
plot.plot()
plot.plots = {
't_ms': 'counting'
}
# Show the plot on the visualizing screen
plot.plot()
O.save(AFTER_sampleName + ".yade.gz")
Edited by Ezzedine Sousan