diff --git a/DCNumPro_1.0.py b/DCNumPro_1.0.py
new file mode 100644
index 0000000000000000000000000000000000000000..9c9917ea33d44dab7e46df67c46ccb0cea659596
--- /dev/null
+++ b/DCNumPro_1.0.py
@@ -0,0 +1,4132 @@
+""" Configuration start """
+# Perform on cluster? False = no, True = yes
+cluster = False
+# Graphically display particles in 3D? False = no, True = yes
+display_particles = False
+# Final test for overlaps? False = no, True = yes
+final_test = True
+# Display initial distribution vs. distribution of placed particles? False = no, True = yes
+display_distribution = True
+# Display detailed progress: status of particle at each step? False = no, True = yes
+# If False, only overall percentage progress will be displayed
+display_detailed_progress = False
+# Display solver warnings? False = no, True = yes
+display_solver_warnings = False
+# Output Excel file
+name_excel_file = 'Test_300'
+# Number of particles
+n_particles_sim = 500
+# Number of classes
+n_classes = 10
+# Median
+mu = 500
+# Standard deviation
+sigma = 130
+# Estimated hold-up
+holdup_est = 0.6
+# Random seed
+seed = 5
+# Number of cores for multiprocessing. 0 = automated determination of available cores
+n_cores = 6
+# Display particle within its cell? False = no, True = yes
+# Is used for debugging or illustriation
+display_particle_cell = False
+# Specify holdup steps
+holdup_delta = 0.025
+# Stepsize for random loose packing simulation
+stepsize = 1
+""" Configuration end """
+""" Import packages """
+import os
+import time
+import sys
+import math
+import random
+import xlwt
+import warnings
+import multiprocessing
+import numpy as np
+import pandas as pd
+import sympy as sy
+import scipy.stats as stats
+import matplotlib.pyplot as plt
+from scipy.integrate import quad
+from scipy.stats import norm
+from scipy.stats import lognorm
+from cycler import cycler
+from functools import reduce
+from scipy.optimize import fsolve
+from scipy.optimize import minimize
+from tess import Container
+from functools import partial
+from scipy.spatial import Delaunay
+from scipy.spatial import ConvexHull
+from mpl_toolkits.mplot3d.art3d import Poly3DCollection
+from tqdm import tqdm
+from collections import defaultdict
+from collections import Counter
+""" Simulation starts from here """
+# Solver warning messages
+if display_solver_warnings == False:
+ warnings.filterwarnings("ignore")
+# If calculation on cluster overwrite parameters above
+# Specify parameter in the .job file
+if cluster == True:
+ import argparse
+ parser = argparse.ArgumentParser(description='Train Mask R-CNN to detect particles.')
+ parser.add_argument('--name_excel_file', required=True)
+ parser.add_argument('--n_particles_sim', required=True, type=int)
+ parser.add_argument('--n_classes', required=True, type=int)
+ parser.add_argument('--holdup_est', required=True, type=float)
+ parser.add_argument('--seed', required=True, type=int)
+ parser.add_argument('--n_cores', required=True, type=int)
+ parser.add_argument('--holdup_delta', required=True, type=float)
+ parser.add_argument('--stepsize', required=True, type=float)
+ parser.add_argument('--mu', required=True, type=float)
+ parser.add_argument('--sigma', required=True, type=float)
+ args = parser.parse_args()
+ NAME_EXCEL_FILE = (args.name_excel_file + '.xls')
+ N_PARTICLES_SIM = args.n_particles_sim
+ N_CLASSES = args.n_classes
+ HOLDUP_EST = args.holdup_est
+ SEED = args.seed
+ N_CORES = args.n_cores
+ HOLDUP_DELTA = args.holdup_delta
+ STEPSIZE = args.stepsize
+ mu = args.mu
+ sigma = args.sigma
+ ROOT_DIR = os.path.join("/rwthfs/rz/cluster", os.path.abspath("../output"))
+else:
+ NAME_EXCEL_FILE = (name_excel_file + '.xls')
+ N_PARTICLES_SIM = int(n_particles_sim)
+ N_CLASSES = int(n_classes)
+ HOLDUP_EST = float(holdup_est)
+ SEED = int(seed)
+ N_CORES = int(n_cores)
+ HOLDUP_DELTA = float(holdup_delta)
+ STEPSIZE = stepsize
+ ROOT_DIR = os.path.abspath("../output")
+
+""" Functions """
+def calc_distribution(N_PARTICLES_SIM, N_CLASSES):
+ """
+ Generate nuber distribution of particles based on the defined distribution parameters mu and sigma calculated from D_32.
+ Args:
+ D_32 (list of floats): Sauter mean diameter.
+ N_PARTICLES_SIM (float): Number of particles.
+ N_CLASSES (float): Number of classes.
+ Returns:
+ d_particles (list of floats): Particle diameters.
+ """
+ # CDF parameter for log-normal DSD
+ # Mean of log values
+ mu_log = (np.log((mu**2)/(((sigma**2)+(mu**2))**0.5)))
+ # Standard deviation of log values
+ sigma_log = ((np.log(((sigma**2)/(mu**2))+1))**0.5)
+ # DSD boundaries which are d_05 and d_95 corresponding to the 5th and 95th percentile of the DSD
+ d_05, d_95 = calc_DSD_boundaries(mu_log, sigma_log)
+ # Class interval
+ delta_class = (d_95-d_05)/N_CLASSES
+ # Class diameters
+ d_classes = np.linspace((d_05+delta_class/2), (d_95-delta_class/2), N_CLASSES)
+ # Determine minimum and maximum class diameter
+ d_min = np.min(d_classes)
+ d_max = np.max(d_classes)
+ # Class boundaries
+ boundaries = [d_05 + i * delta_class for i in range(N_CLASSES+1)]
+ # Class sizes
+ class_sizes, q0_initial, q3_initial, volume_total = calc_class_sizes(mu_log, sigma_log, boundaries, N_PARTICLES_SIM, d_classes)
+ # Sauter mean diameter
+ D_32 = np.sum(q0_initial * d_classes**3) / np.sum(q0_initial * d_classes**2)
+ # Actual number of particles
+ n_particles_init = int(np.sum(class_sizes))
+ # Create array q0 [class_sizes, d_classes]
+ q0 = np.concatenate((class_sizes.reshape(-1,1), d_classes.reshape(-1,1)), axis=1)
+
+ return q0, n_particles_init, q0_initial, q3_initial, d_classes, d_min, d_max, volume_total, d_05, d_95, delta_class, D_32
+
+def calc_DSD_boundaries(mu, sigma):
+ """
+ Calculate DSD boundaries d_05 and d_95 corresponding to the 5th and 95th percentile of the DSD.
+ Args:
+ mu (float): Mean of the logarithm of the values.
+ sigma (float): Standard deviation of the logarithm of the values.
+ Returns:
+ d_min (float): Minimum diameter.
+ d_max (float): Maximum diameter.
+ """
+ # Find the z-value corresponding to the 5th and 95th percentile of the standard normal distribution
+ d_05_percentile = norm.ppf(0.05)
+ d_95_percentile = norm.ppf(0.95)
+ # Calculate the 5th and 95th percentile of the log-normal distribution
+ d_05 = np.exp(mu + sigma * d_05_percentile)
+ d_95 = np.exp(mu + sigma * d_95_percentile)
+
+ return d_05, d_95
+
+def calc_class_sizes(mu, sigma, boundaries, N_PARTICLES_SIM, d_classes):
+ """
+ Calculate class sizes for a log-normal distribution with parameters mu and sigma.
+ Args:
+ mu (float): Mean of the logarithm of the values.
+ sigma (float): Standard deviation of the logarithm of the values.
+ boundaries (list of floats): Boundaries of the classes.
+ Returns:
+ class_sizes (list of floats): Quantity for each class defined by the boundaries.
+ """
+ # Ensure boundaries are sorted
+ boundaries = sorted(boundaries)
+ # Create a log-normal distribution object
+ # Scale parameter for lognorm is sigma
+ # Scale parameter exp(mu) for shifting the distribution
+ dist = lognorm(s=sigma, scale=np.exp(mu))
+ # Calculate probabilities for each interval
+ q3_initial = []
+ for i in range(len(boundaries) - 1):
+ lower_bound = boundaries[i]
+ upper_bound = boundaries[i + 1]
+ lb = dist.cdf(lower_bound)
+ ub = dist.cdf(upper_bound)
+ # Probability that a value falls within the interval [lower_bound, upper_bound]
+ prob = dist.cdf(upper_bound) - dist.cdf(lower_bound)
+ # Calculate class sizes
+ q3_initial.append(prob)
+ q3_initial = np.array(q3_initial)
+ # Convert class volume probabilities into number probabilities
+ q0_initial = q3_initial / d_classes**3
+ # Normalize number probabilities.
+ # It ignores that the class volume probabilities include only 90vol% of the distribution.
+ # However, since we use number probabilities to distribute a certain number of particles
+ # within these considered 90vol%, this calculation is valid
+ q0_initial /= np.sum(q0_initial)
+ # Distribute the nubmer of particles to the classes considering number probabilities and round
+ class_sizes = N_PARTICLES_SIM * np.array(q0_initial)
+ np.round(class_sizes, decimals=0, out=class_sizes)
+ # Total particle volume
+ volume_total = np.sum(class_sizes * (math.pi * d_classes**3 / 6))
+
+ return class_sizes, q0_initial, q3_initial, volume_total
+
+def calc_particles(q0, n_particles_init):
+ """
+ Fill an array with particles from number distribution q0 and determine their drop coordinates randomly.
+ Args:
+ q0 (array): Number distribution.
+ n_particles_init (float): Actualnumber of particles.
+ Returns:
+ M_particles_initial (array): Array with particles and their drop coordinates
+ """
+ # Array columns:
+ # 0: x-coordinate, 1: y-coordinate, 2: z-coordinate, 3: index, 4: diameter, 5: class number (0 to n)
+ M_particles_initial = np.zeros(shape=(int(n_particles_init), 6))
+ # Assign particle diameters and class numbers
+ j = 0
+ for k in range(len(q0)):
+ n = int(q0[k, 0])
+ for i in range(n):
+ M_particles_initial[j, 4] = q0[k, 1]
+ M_particles_initial[j, 5] = k
+ j = j + 1
+ # Assign indices and randomly generated coordinates
+ for i in range(len(M_particles_initial)):
+ M_particles_initial[i, 0] = random.uniform(-A_0 + (M_particles_initial[i, 4]) / 2, A_0 - (M_particles_initial[i, 4]) / 2)
+ M_particles_initial[i, 1] = random.uniform(-A_0 + (M_particles_initial[i, 4]) / 2, A_0 - (M_particles_initial[i, 4]) / 2)
+ M_particles_initial[i, 2] = None
+ M_particles_initial[i, 3] = i
+ np.random.shuffle(M_particles_initial)
+
+ return M_particles_initial
+
+def particle_info(M_particles, index):
+ # x-coordinate
+ x = M_particles[index, 0]
+ # y-coordinate
+ y = M_particles[index, 1]
+ # z-coordinate
+ z = M_particles[index, 2]
+ # Index
+ p_index = M_particles[index, 3]
+ # Diameter
+ p_diameter = M_particles[index, 4]
+ # Class
+ p_class = int(M_particles[index, 5])
+
+ return x, y, z, p_index, p_diameter, p_class
+
+def find_z_drop(y, x, M_particle_steady, test_range, d_max_steady):
+ """ Determine drop position of particles """
+ # Define testing area
+ y_min = y - test_range
+ y_max = y + test_range
+ x_min = x - test_range
+ x_max = x + test_range
+ # Check particles in steady state for interraction with testing area:
+ # in y direction
+ particle_to_y = np.argwhere(np.logical_and(M_particle_steady[:, 2] > y_min, M_particle_steady[:, 2] < y_max))
+ # in x direction
+ particle_to_x = np.argwhere(np.logical_and(M_particle_steady[:, 3] > x_min, M_particle_steady[:, 3] < x_max))
+ # Combine the results of both tests (arrays with indices of particles which interract with testing area)
+ particle_to_check = reduce(np.intersect1d, (particle_to_y, particle_to_x))
+ # Convert array to list
+ particle_to_check = np.array(particle_to_check).tolist()
+ # Create list with z-coordinates of interracting particles
+ z_list = []
+ for i in particle_to_check:
+ value = M_particle_steady[i, 1]
+ z_list.append(value)
+ # Check list for entries and pick the max value
+ if not z_list:
+ max_z_value = 0
+ else:
+ max_z_value = max(z_list)
+ # Calculate dropping height
+ z = (max_z_value + d_max_steady)*1.2
+
+ return z
+
+def calc_overlap(d_i, z_i, y_i, x_i, M_particle_steady, test_range):
+ """
+ Function that checks overlaps of an unstable particle.
+ It also double check the collision criteria of the calc_sedimentation function.
+ If calc_sedimentation function works properly, only one overlapping should be identified.
+ """
+ # Variable to identifiycollisions with walls
+ walls = 0
+ # Counting variable for number of collisions with steady state particles
+ n_overlaps = 0
+ # Array with indices of particles in steady state which overlap with the considered particle
+ overlaps = []
+ """ Define region of interrest """
+ z_min = z_i - test_range
+ if z_min < 0:
+ z_min = 0
+ z_max = z_i + test_range
+ y_min = y_i - test_range
+ y_max = y_i + test_range
+ x_min = x_i - test_range
+ x_max = x_i + test_range
+ """ Checking particles within the testing range for possible overlappings """
+ # Arrays in x, y and z direction with indices of particles, which could overlap with the considered particle
+ x_array = np.argwhere(np.logical_and(M_particle_steady[:, 3] > x_min, M_particle_steady[:, 3] < x_max))
+ y_array = np.argwhere(np.logical_and(M_particle_steady[:, 2] > y_min, M_particle_steady[:, 2] < y_max))
+ z_array = np.argwhere(np.logical_and(M_particle_steady[:, 1] > z_min, M_particle_steady[:, 1] < z_max))
+ # Combine the arrays
+ particle_to_check = reduce(np.intersect1d, (x_array, y_array, z_array))
+ # Convert array to list
+ particle_to_check = np.array(particle_to_check).tolist()
+ """ Check each particle that may overlap for an actual overlapping """
+ for i in particle_to_check:
+ # Retrieve information about particles that may overlap with the particle under consideration from the steady state array
+ # Coordinates
+ z_j = float(M_particle_steady[i, 1])
+ y_j = float(M_particle_steady[i, 2])
+ x_j = float(M_particle_steady[i, 3])
+ # Diameter
+ d_j = float(M_particle_steady[i, 5])
+ # distance between the particle's centers when they collide
+ dist_required = (d_i + d_j) / 2
+ # Actual distance between the particle's centers
+ distance_particles = math.sqrt((z_i-z_j)**2 + (y_i-y_j)**2 + (x_i-x_j)**2)
+ # Check whether particles overlap and store the overlapping particle in array
+ if dist_required - distance_particles > 1e-9:
+ overlaps.append(i)
+ # Remove duplicates from a list
+ overlaps = list(set(overlaps))
+ # Determine number of collisions
+ n_overlaps = len(overlaps)
+ """ Now checking for overlappings with walls """
+ # Wall in positive x direction
+ if (x_i + d_i / 2) >= A_0:
+ walls = 1
+ # Wall in negative x direction
+ elif (x_i - d_i / 2) <= -A_0:
+ walls = -1
+ # Wall in positive y direction
+ elif (y_i + d_i / 2) >= A_0:
+ walls = 2
+ # Wall in negative y direction
+ elif (y_i - d_i / 2) <= -A_0:
+ walls = -2
+ # Wall in positive x direction and positive y direction
+ elif (x_i + d_i / 2) >= A_0 and (y_i + d_i / 2) >= A_0:
+ walls = 3
+ # Wall in negative x direction and positive y direction
+ elif (x_i - d_i / 2) <= -A_0 and (y_i + d_i / 2) >= A_0:
+ walls = -3
+ # Wall in positive x direction and negative y direction
+ elif (x_i + d_i / 2) >= A_0 and (y_i - d_i / 2) <= -A_0:
+ walls = 4
+ # Wall in negative x direction and negaive y direction
+ elif (x_i - d_i / 2) <= -A_0 and (y_i - d_i / 2) <= -A_0:
+ walls = -4
+
+ return n_overlaps, overlaps, walls
+
+def position_update1(y_i, x_i, d_i, M_particle_steady, overlaps):
+ """ Function that updates the position of the considered particle overlapping with one steady state particle. """
+ # Coordinates of the steady state particle
+ y1 = float(M_particle_steady[overlaps[0], 2])
+ x1 = float(M_particle_steady[overlaps[0], 3])
+ # Diameter of the steady state particle
+ d1 = float(M_particle_steady[overlaps[0], 5])
+ # Required distance between the centers of the considered particle and the steady state particle
+ dist_1i = (d1 + d_i) / 2
+ # Actual distance between the centers of the considered particle and the steady state particle in x and y coordinate
+ delta_x = x_i - x1
+ delta_y = y_i - y1
+ dist_1i_current = np.sqrt(delta_x**2 + delta_y**2)
+ # Calculate the scaling factor
+ scaling_factor = dist_1i / dist_1i_current
+ # Calculate the new position
+ x = x1 + delta_x * scaling_factor
+ y = y1 + delta_y * scaling_factor
+
+ return y, x
+
+def position_update2(z_i, y_i, x_i, d_i, M_particle_steady, overlaps):
+ """
+ Function to update the position of the considered particle if it overlaps with two particles.
+ There are two mathematical solutions of the new position, but only one logical.
+ A solver is used to determine these solutions.
+ First, an intersection circle plane is calculated corresponding to the intersection area of two auxiliary spheres.
+ The intersection plane is used to identify the guesses of the both solutions of the considered particle's new position.
+ Moreover, the plane is used to check wheter the solver solutions are correct.
+ The centers of the auxiliary spheres corresponds to the center coordinates of the both steady state particles.
+ The radii of the auxiliary spheres corresponds to the sums of the steady state particles' radii and the considered particle's radius.
+ To identify the logical solution of the both solutions, the distances between the original position and the new positions are calculated.
+ The new position with the lowest distance to the original position is the logical solution.
+ """
+ # Used in case an error occurs during position update.
+ error_position_update2 = 0
+ # Coordinates of the first particle in steady state
+ z1 = float(M_particle_steady[overlaps[0], 1])
+ y1 = float(M_particle_steady[overlaps[0], 2])
+ x1 = float(M_particle_steady[overlaps[0], 3])
+ # Diameter of the first particle in steady state
+ d1 = float(M_particle_steady[overlaps[0], 5])
+ # Center coordinates of the first auxiliary sphere
+ c1 = np.array([z1, y1, x1])
+ # Radius of the first auxiliary sphere
+ r1 = (d_i + d1) / 2
+ # Coordinates of the second particle in steady state
+ z2 = float(M_particle_steady[overlaps[1], 1])
+ y2 = float(M_particle_steady[overlaps[1], 2])
+ x2 = float(M_particle_steady[overlaps[1], 3])
+ # Diameter of the second particle in steady state
+ d2 = float(M_particle_steady[overlaps[1], 5])
+ # Center coordinates of the second auxiliary sphere
+ c2 = np.array([z2, y2, x2])
+ # Radius of the second auxiliary sphere
+ r2 = (d_i + d2) / 2
+ """ Calculation of the intersection circle plane """
+ # Vector between the centers
+ V_centers_12 = c2 - c1
+ # Distance between the centers
+ dist_centers_12 = np.linalg.norm(V_centers_12)
+ # Normal vector of the intersection circle plane
+ V_unit_centers_12 = V_centers_12 / dist_centers_12
+ # Distance from the first sphere's center to the intersection circle plane
+ dist_c1_plane = (r1**2 - r2**2 + dist_centers_12**2) / (2 * dist_centers_12)
+ # Center of the intersection circle (point on the intersection plane closest to sphere1)
+ c_ic = c1 + dist_c1_plane * V_unit_centers_12
+ """
+ Calculation of the reflected position of the considered particle.
+ First, the vector between the centers of the considered particle and the intersection circle is determined.
+ Then, this vector is reflected such a way that the projection of this vector on the yx plane is rotated by 180 degrees on the yx plane.
+ The coordiantes of the reflected position are used as second guess for the solver to determine the second solver solution.
+ """
+ # Vector between the centers of the considered particle and the intersection circle
+ V_centers_ici = np.array([z_i, y_i, x_i]) - np.array(c_ic)
+ # Reflected vector
+ V_reflected = np.array([V_centers_ici[0], -V_centers_ici[1], -V_centers_ici[2]])
+ # Reflected position
+ P_reflected = np.array(c_ic) + V_reflected
+ # Coordinates for the initial guess
+ x_01 = x_i
+ y_01 = y_i
+ x_02 = P_reflected[2]
+ y_02 = P_reflected[1]
+ # Initial guess values
+ initial_guess1 = [x_01, y_01]
+ initial_guess2 = [x_02, y_02]
+ """
+ Solver to find possible new positions of the considered particle
+ The first guess is the logical one. The second is just to determine the second solver solution.
+ The second solver solution is just to check wheter the first one is the right one
+ """
+ x_s1, y_s1 = fsolve(solver_position2, initial_guess1, args=(c1, c2, r1, r2, z_i))
+ x_s2, y_s2 = fsolve(solver_position2, initial_guess2, args=(c1, c2, r1, r2, z_i))
+ """ Check wheter the solver found two different solutions or the same twice """
+ # Error value is 1 if the solver found the same position twice
+ if x_s1 == x_s2 and y_s1 == y_s2:
+ error_position_update2 = 1
+ # Exclude values behind the simulation walls
+ if x_s1 < -A_0 or x_s1 > A_0:
+ x_s1 = None
+ if x_s2 < -A_0 or x_s2 > A_0:
+ x_s2 = None
+ if y_s1 < -A_0 or y_s1 > A_0:
+ y_s1 = None
+ if y_s2 < -A_0 or y_s2 > A_0:
+ y_s2 = None
+ """
+ Determine the actual new position of the considered particle from the both solver solutions.
+ Must be the closest position to the original position of the considered particle
+ """
+ # Choose the new position based on the validity of values
+ if y_s1 == None or x_s1 == None:
+ if y_s2 != None and x_s2 != None:
+ x = x_s2
+ y = y_s2
+ else:
+ error_position_update2 = 3
+ return y_i, x_i, error_position_update2
+ elif y_s2 == None or x_s2 == None:
+ if y_s1 != None or x_s1 != None:
+ x = x_s1
+ y = y_s1
+ else:
+ error_position_update2 = 3
+ return y_i, x_i, error_position_update2
+ else:
+ # Vectors between the centers of the considered particle's original position and the new positions
+ V_centers_new1 = np.array([z_i, y_i, x_i]) - np.array([z_i, y_s1, x_s1])
+ V_centers_new2 = np.array([z_i, y_i, x_i]) - np.array([z_i, y_s2, x_s2])
+ # Distances between these centers
+ dist_centers_new1 = np.linalg.norm(V_centers_new1)
+ dist_centers_new2 = np.linalg.norm(V_centers_new2)
+ # Choose new position based on the shortest distance
+ if dist_centers_new1 <= dist_centers_new2:
+ x = x_s1
+ y = y_s1
+ else:
+ x = x_s2
+ y = y_s2
+ """
+ Check wheter the new positions' coordinates are located on the intersection plane
+ Just double checking wheter the solver found correct solutions
+ The equation of the intersection plane is Az+By+Cx+D=0, where A, B, C and D are the equation coefficients.
+ """
+ # Determine the plane equation coefficients
+ # The normal vector to the plane is the same as the V_centers_12
+ coeff_A, coeff_B, coeff_C = V_unit_centers_12
+ # Calculate D using the point on the plane
+ coeff_D = -(coeff_A * c_ic[0] + coeff_B * c_ic[1] + coeff_C * c_ic[2])
+ # Check whether the new positions of considered particle is located on the intersection plane
+ on_plane = check_point_on_plane(z_i, y, x, coeff_A, coeff_B, coeff_C, coeff_D)
+ # Error value is 2 if the position of the first solver solution is not located on the intersection plane
+ if on_plane == False:
+ error_position_update2 = 2
+
+ return y, x, error_position_update2
+
+def solver_position2(params, c1, c2, r1, r2, z_i):
+ """ Solver to calculate the new position of the considered particle if it overlaps with two particles """
+ # Solver variables
+ x, y = params
+ # Equation of a circle. The circles coresponds here to the auxiliary spheres
+ eq1 = (x - c1[2])**2 + (y - c1[1])**2 + (z_i - c1[0])**2 - r1**2
+ eq2 = (x - c2[2])**2 + (y - c2[1])**2 + (z_i - c2[0])**2 - r2**2
+ return [eq1, eq2]
+
+def check_point_on_plane(z_i, y, x, coeff_A, coeff_B, coeff_C, coeff_D):
+ """ Function that checks whether a point lies on a plane or not (plane equation in coordinate form) """
+ # Calculate the left-hand side of the plane equation
+ check_zero = z_i*coeff_A + y*coeff_B + x*coeff_C + coeff_D
+ # Check if the left-hand side is approximately zero (account for floating-point precision)
+ on_plane = abs(check_zero) < 1e-9
+ return on_plane
+
+def position_update3(z_i, y_i, x_i, d_i, M_particle_steady, overlaps):
+ """
+ Function to update the position of the considered particle if it overlaps with three particles.
+ There are two mathematical solutions of the new position, but only one logical.
+ A solver is used to determine these solutions.
+ First, a triangle plane is calculated which crosses the centers of the three steady state particles.
+ Moreover, the centroid of the triangle is determined.
+ The triangle plane is used to identify the guesses of the both solutions of the considered particle's new position.
+ It is tested whether the both solutions differ and if they are correct.
+ To identify the logical solution of the both solutions, the z-coordinates of the both soultions and the centroid are compared.
+ The new position with the z-coordinate above the centroid is the logical solution.
+ """
+ # Used in case an error occurs during position update.
+ error_position_update3 = 0
+ # Coordinates of the first particle in steady state
+ z1 = float(M_particle_steady[overlaps[0], 1])
+ y1 = float(M_particle_steady[overlaps[0], 2])
+ x1 = float(M_particle_steady[overlaps[0], 3])
+ # Diameter of the first particle in steady state
+ d1 = float(M_particle_steady[overlaps[0], 5])
+ # Center coordinates of the first auxiliary sphere
+ c1 = np.array([z1, y1, x1])
+ # Radius of the first auxiliary sphere
+ r1 = (d_i + d1) / 2
+ # Coordinates of the second particle in steady state
+ z2 = float(M_particle_steady[overlaps[1], 1])
+ y2 = float(M_particle_steady[overlaps[1], 2])
+ x2 = float(M_particle_steady[overlaps[1], 3])
+ # Diameter of the second particle in steady state
+ d2 = float(M_particle_steady[overlaps[1], 5])
+ # Center coordinates of the second auxiliary sphere
+ c2 = np.array([z2, y2, x2])
+ # Radius of the second auxiliary sphere
+ r2 = (d_i + d2) / 2
+ # Coordinates of the third particle in steady state
+ z3 = float(M_particle_steady[overlaps[2], 1])
+ y3 = float(M_particle_steady[overlaps[2], 2])
+ x3 = float(M_particle_steady[overlaps[2], 3])
+ # Diameter of the third particle in steady state
+ d3 = float(M_particle_steady[overlaps[2], 5])
+ # Center coordinates of the third auxiliary sphere
+ c3 = np.array([z3, y3, x3])
+ # Radius of the third auxiliary sphere
+ r3 = (d_i + d3) / 2
+ """ Calculation of the triangle's plane """
+ # Calculate vectors AB and AC
+ AB = np.array([z2 - z1, y2 - y1, x2 - x1])
+ AC = np.array([z3 - z1, y3 - y1, x3 - x1])
+ # Compute the normal vector to the plane (cross product of AB and AC)
+ V_triangle_norm = np.cross(AB, AC)
+ A, B, C = V_triangle_norm
+ # Calculate D using the plane equation Ax + By + Cz + D = 0
+ # Using point A (x1, y1, z1) to find D
+ D = -(A * z1 + B * y1 + C * x1)
+ """
+ Calculation of the reflected position of the considered particle.
+ First, the distance between the center of the considered particle and the triangle's plane is determined.
+ Then, the position of the considered particle is reflected about the triangle's plane
+ The coordiantes of the reflected position are used as second guess for the solver to determine the second solver solution.
+ """
+ # Calculate the distance from the considered particle's center to the plane
+ distance_pi = (A * z_i + B * y_i + C * x_i + D) / np.linalg.norm(V_triangle_norm)
+ # Calculate the reflection point on the plane
+ P_reflected = np.array([z_i, y_i, x_i]) - 2 * distance_pi * V_triangle_norm / np.linalg.norm(V_triangle_norm)
+ """ Solver to find possible new position of the considered particle """
+ # Coordinates for the initial guess
+ x_01 = x_i
+ y_01 = y_i
+ z_01 = z_i
+ x_02 = P_reflected[2]
+ y_02 = P_reflected[1]
+ z_02 = P_reflected[0]
+ # Initial guess values
+ initial_guess1 = [x_01, y_01, z_01]
+ initial_guess2 = [x_02, y_02, z_02]
+ # Solver
+ x_s1, y_s1, z_s1 = fsolve(solver_position3, initial_guess1, args=(c1, c2, c3, r1, r2, r3))
+ x_s2, y_s2, z_s2 = fsolve(solver_position3, initial_guess2, args=(c1, c2, c3, r1, r2, r3))
+ """ Check wheter the solver found two different solutions or the same twice """
+ # Error value is 1 if the solver found the same position twice
+ if x_s1 == x_s2 and y_s1 == y_s2 and z_s1 == z_s2:
+ error_position_update3 = 1
+ # Exclude values behind the simulation walls
+ if x_s1 < -A_0 or x_s1 > A_0:
+ x_s1 = None
+ if x_s2 < -A_0 or x_s2 > A_0:
+ x_s2 = None
+ if y_s1 < -A_0 or y_s1 > A_0:
+ y_s1 = None
+ if y_s2 < -A_0 or y_s2 > A_0:
+ y_s2 = None
+ """
+ Determine the actual new position of the considered particle from the both solver solutions.
+ Must be the closest position to the original position of the considered particle
+ """
+ # Z-coordinate of polygon's centroid
+ z_cen = (z1 + z2 + z3) / 3
+ # If the center of the considered particle is located above the centroid,
+ # the coordinates of the first solver solution are assigned to the center of the considered particle
+ if x_s1 != None and y_s1 != None and z_s1 > z_cen:
+ z = z_s1
+ y = y_s1
+ x = x_s1
+ # If the center of the considered particle is located above the centroid,
+ # the coordinates of the second solver solution are assigned to the center of the considered particle
+ elif x_s2 != None and y_s2 != None and z_s2 > z_cen:
+ z = z_s2
+ y = y_s2
+ x = x_s2
+ else:
+ error_position_update3 = 2
+ return z_i, y_i, x_i, error_position_update3
+ """
+ Check wheter the distances between the centers of the three steady state particles
+ and the new center position of the considered particle are in the samerange as the sums of their radii.
+ """
+ same_range = check_distances3(r1, r2, r3, c1, c2, c3, z, y, x)
+ # Error value is 2 if the distance between the steady state particles and the position of the solver solution
+ # is not in the same range as the sums of their radii
+ if same_range == False:
+ error_position_update3 = 3
+ return z_i, y_i, x_i, error_position_update3
+
+ return z, y, x, error_position_update3
+
+def solver_position3(params, c1, c2, c3, r1, r2, r3):
+ """ Solver to calculate the new position of the considered particle if it overlaps with three particles """
+ # Solver variables
+ x, y, z = params
+ # Equation of a circle. The circles coresponds here to the auxiliary spheres
+ eq1 = (x - c1[2])**2 + (y - c1[1])**2 + (z - c1[0])**2 - r1**2
+ eq2 = (x - c2[2])**2 + (y - c2[1])**2 + (z - c2[0])**2 - r2**2
+ eq3 = (x - c3[2])**2 + (y - c3[1])**2 + (z - c3[0])**2 - r3**2
+ return [eq1, eq2, eq3]
+
+def check_distances3(r1, r2, r3, c1, c2, c3, z_s, y_s, x_s):
+ """ Function that checks wheter the distances between the centers of the three steady state particles
+ and the new center position of the considered particle are in the same range as the sums of their radii. """
+ # True if distances are in the same range as the sums of radii
+ same_range = True
+ # Vectors between the centers of the steady state particles and the new center positions of the considered particle
+ V_centers_1i = np.array([z_s, y_s, x_s]) - c1
+ V_centers_2i = np.array([z_s, y_s, x_s]) - c2
+ V_centers_3i = np.array([z_s, y_s, x_s]) - c3
+ # Distance between the centers of the steady state particles and the new center position of the considered particle
+ dist_centers_1i = np.linalg.norm(V_centers_1i)
+ dist_centers_2i = np.linalg.norm(V_centers_2i)
+ dist_centers_3i = np.linalg.norm(V_centers_3i)
+ # Check if the left-hand side is approximately zero (account for floating-point precision)
+ if abs(dist_centers_1i-r1) > 1e-9 or abs(dist_centers_2i-r2) > 1e-9 or abs(dist_centers_3i-r3) > 1e-9:
+ same_range = False
+
+ return same_range
+
+def check_steadystate3(z_i, y_i, x_i, M_particle_steady, overlaps):
+ """
+ Function to check the steady state.
+ As soon as considered particle collides with three steady state particles, a steady state check is performed.
+ Considered particle is in steady state if its center lies within the y-x-plane of a triangle located between the
+ steady state particle centers.
+ Barycentric coordinate method is used to determine if the center of the considered particle lies within the polygon.
+ Moreover, barycentric coordinates are used to determine the subdomain of the triangle where the center lies.
+ Based on the subdomain, the particle rolling direction is identified.
+ """
+ # Used in case an error occurs during steadystate check.
+ error_check_steadystate3 = 0
+ # Variable for steady state check condition
+ steady_check = 0
+ # Coordinates of the first particle in steady state
+ z1 = float(M_particle_steady[overlaps[0], 1])
+ y1 = float(M_particle_steady[overlaps[0], 2])
+ x1 = float(M_particle_steady[overlaps[0], 3])
+ # Coordinates of the second particle in steady state
+ z2 = float(M_particle_steady[overlaps[1], 1])
+ y2 = float(M_particle_steady[overlaps[1], 2])
+ x2 = float(M_particle_steady[overlaps[1], 3])
+ # Coordinates of the third particle in steady state
+ z3 = float(M_particle_steady[overlaps[2], 1])
+ y3 = float(M_particle_steady[overlaps[2], 2])
+ x3 = float(M_particle_steady[overlaps[2], 3])
+ # Z-coordinate of triangle's centroid
+ z_cen = (z1 + z2 + z3) / 3
+ # Check whether the center of the considered particle is located above the centroid
+ # Double checking at this point, since already handled in the position_update3 function.
+ # Error value is 1 when the center is located below the centroid
+ if z_i < z_cen:
+ error_check_steadystate3 = 1
+ return steady_check, overlaps, error_check_steadystate3
+ """ Barycentric coordinate technique """
+ # Denominator
+ denominator = ((y2 - y3) * (x1 - x3) + (x3 - x2) * (y1 - y3))
+ # Check whether the triangle points are collinear
+ # Error value is 2 when the triangle points are collinear
+ if denominator == 0:
+ error_check_steadystate3 = 2
+ return steady_check, overlaps, error_check_steadystate3
+ # Barycentric coordinates
+ a = ((y2 - y3) * (x_i - x3) + (x3 - x2) * (y_i - y3)) / denominator
+ b = ((y3 - y1) * (x_i - x3) + (x1 - x3) * (y_i - y3)) / denominator
+ c = 1 - a - b
+ # Determine the subdomain of the triangle and thus wheter it is in steady state or in which direction it is rolling
+ if (0 <= a <= 1) and (0 <= b <= 1) and (0 <= c <= 1):
+ steady_check = 3
+ elif a < 0 and b >= 0 and c >= 0:
+ # Rolling along steady state particle Nr.2 and Nr.3
+ overlaps = np.delete(overlaps, 0)
+ steady_check = 2
+ elif b < 0 and a >= 0 and c >= 0:
+ # Rolling along steady state particle Nr.1 and Nr.3
+ overlaps = np.delete(overlaps, 1)
+ steady_check = 2
+ elif c < 0 and a >= 0 and b >= 0:
+ # Rolling along steady state particle Nr.1 and Nr.2
+ overlaps = np.delete(overlaps, 2)
+ steady_check = 2
+ elif a < 0 and b < 0 and c > 1:
+ # Rolling along steady state particle Nr.3
+ overlaps = np.delete(overlaps, [0, 1])
+ steady_check = 1
+ elif a < 0 and c < 0 and b > 1:
+ # Rolling along steady state particle Nr.2
+ overlaps = np.delete(overlaps, [0, 2])
+ steady_check = 1
+ elif b < 0 and c < 0 and a > 1:
+ # Rolling along steady state particle Nr.1
+ overlaps = np.delete(overlaps, [1, 2])
+ steady_check = 1
+ else:
+ # Error value is 3 when the subdomain of the triangle could not be determined correctly
+ error_check_steadystate3 = 3
+
+ return steady_check, overlaps, error_check_steadystate3
+
+def position_update1w(z_i, y_i, x_i, d_i, M_particle_steady, overlaps, walls):
+ """
+ Function to update the position of the considered particle if it overlaps with one particle and one wall.
+ There are two mathematical solutions of the new position, but only one logical.
+ A solver is used to determine these solutions.
+ First, an intersection circle plane is calculated corresponding to the intersection area of one auxiliary sphere and the wall.
+ The intersection plane is used to identify the guesses of the both solutions of the considered particle's new position.
+ Moreover, the plane is used to check wheter the solver solutions are correct.
+ The center of the auxiliary sphere corresponds to the center coordinates of the steady state particle.
+ The radius of the auxiliary sphere corresponds to the sums of the steady state particle's radius and the considered particle's radius.
+ To identify the logical solution of the both solutions, the distances between the original position and the new positions are calculated.
+ The new position with the lowest distance to the original position is the logical solution.
+ """
+ # Used in case an error occurs during position update.
+ error_position_update1w = 0
+ # Coordinates of the particle in steady state
+ z1 = float(M_particle_steady[overlaps[0], 1])
+ y1 = float(M_particle_steady[overlaps[0], 2])
+ x1 = float(M_particle_steady[overlaps[0], 3])
+ # Diameter of the particle in steady state
+ d1 = float(M_particle_steady[overlaps[0], 5])
+ # Center coordinates of the auxiliary sphere
+ c1 = np.array([z1, y1, x1])
+ # Radius of the auxiliary sphere
+ r1 = (d_i + d1) / 2
+ # Coordinates of the wall
+ z_w1 = z1
+ if walls == 1:
+ x_w1 = A_0
+ y_w1 = y1
+ elif walls == -1:
+ x_w1 = -A_0
+ y_w1 = y1
+ elif walls == 2:
+ x_w1 = x1
+ y_w1 = A_0
+ elif walls == -2:
+ x_w1 = x1
+ y_w1 = -A_0
+ # Center coordinates of the wall
+ c_w1 = np.array([z_w1, y_w1, x_w1])
+ # Distance between wall and center of considered particle if they are in contact
+ dist_wi = d_i / 2
+ """ Calculation of the intersection circle plane """
+ # Vector between the centers
+ V_centers_1w = c_w1 - c1
+ # Distance between the centers
+ dist_centers_1w = np.linalg.norm(V_centers_1w)
+ # Normal vector of the intersection circle plane
+ V_unit_centers_1w = V_centers_1w / dist_centers_1w
+ # Distance from the first sphere's center to the intersection circle plane
+ dist_c1_plane = dist_centers_1w - dist_wi
+ # Center of the intersection circle (point on the intersection plane closest to sphere1)
+ c_ic = c1 + dist_c1_plane * V_unit_centers_1w
+ """
+ Calculation of the reflected position of the considered particle.
+ First, the vector between the centers of the considered particle and the intersection circle is determined.
+ Then, this vector is reflected such a way that the projection of this vector on the yx plane is rotated by 180 degrees on the yx plane.
+ The coordiantes of the reflected position are used as second guess for the solver to determine the second solver solution.
+ """
+ # Vector between the centers of the considered particle and the intersection circle
+ V_centers_ici = np.array([z_i, y_i, x_i]) - np.array(c_ic)
+ # Reflected vector
+ V_reflected = np.array([V_centers_ici[0], -V_centers_ici[1], -V_centers_ici[2]])
+ # Reflected position
+ P_reflected = np.array(c_ic) + V_reflected
+ # Coordinates for the initial guess
+ x_01 = x_i
+ y_01 = y_i
+ x_02 = P_reflected[2]
+ y_02 = P_reflected[1]
+ # If contact is with X-wall
+ if walls == 1 or walls == -1:
+ # Initial guess values
+ initial_guess1 = y_01
+ initial_guess2 = y_02
+ # X-coordinate of the new position
+ if walls == 1:
+ x_s = A_0 - d_i / 2
+ elif walls == -1:
+ x_s = -A_0 + d_i / 2
+ """
+ Solver to find possible new positions of the considered particle
+ The first guess is the logical one. The second is just to determine the second solver solution.
+ The second solver solution is just to check wheter the first one is the right one
+ """
+ y_s1 = fsolve(solver_position1w, initial_guess1, args=(c1, r1, z_i, x_s, walls))
+ y_s2 = fsolve(solver_position1w, initial_guess2, args=(c1, r1, z_i, x_s, walls))
+ y_s1 = y_s1[0]
+ y_s2 = y_s2[0]
+ """ Check wheter the solver found two different solutions or the same twice """
+ # Error value is 1 if the solver found the same position twice
+ if y_s1 == y_s2:
+ error_position_update1w = 1
+ # Exclude values behind the simulation walls
+ if y_s1 < -A_0 or y_s1 > A_0:
+ y_s1 = None
+ if y_s2 < -A_0 or y_s2 > A_0:
+ y_s2 = None
+ """
+ Determine the actual new position of the considered particle from the both solver solutions.
+ Must be the closest position to the original position of the considered particle
+ """
+ # Coose the new position based on the validity of values
+ if y_s1 == None and y_s2 != None:
+ y = y_s2
+ elif y_s2 == None and y_s1 != None:
+ y = y_s1
+ elif y_s1 == None and y_s2 == None:
+ error_position_update1w = 3
+ return y_i, x_i, error_position_update1w
+ else:
+ # Vectors between the centers of the considered particle's original position and the new positions
+ V_centers_new1 = np.array([z_i, y_i, x_i]) - np.array([z_i, y_s1, x_s])
+ V_centers_new2 = np.array([z_i, y_i, x_i]) - np.array([z_i, y_s2, x_s])
+ # Distances between these centers
+ dist_centers_new1 = np.linalg.norm(V_centers_new1)
+ dist_centers_new2 = np.linalg.norm(V_centers_new2)
+ # Choose new position based on the shortest distance
+ if dist_centers_new1 <= dist_centers_new2:
+ y = y_s1
+ else:
+ y = y_s2
+ x = x_s
+ # If contact is with Y-wall
+ elif walls == 2 or walls == -2:
+ # Initial guess values
+ initial_guess1 = x_01
+ initial_guess2 = x_02
+ # Y-coordinate of the new position
+ if walls == 2:
+ y_s = A_0 - d_i / 2
+ elif walls == -2:
+ y_s = -A_0 + d_i / 2
+ """
+ Solver to find possible new positions of the considered particle
+ The first guess is the logical one. The second is just to determine the second solver solution.
+ The second solver solution is just to check wheter the first one is the right one
+ """
+ x_s1 = fsolve(solver_position1w, initial_guess1, args=(c1, r1, z_i, y_s, walls))
+ x_s2 = fsolve(solver_position1w, initial_guess2, args=(c1, r1, z_i, y_s, walls))
+ x_s1 = x_s1[0]
+ x_s2 = x_s2[0]
+ """ Check wheter the solver found two different solutions or the same twice """
+ # Error value is 1 if the solver found the same position twice
+ if x_s1 == x_s2:
+ error_position_update1w = 1
+ # Exclude values behind the simulation walls
+ if x_s1 < -A_0 or x_s1 > A_0:
+ x_s1 = None
+ if x_s2 < -A_0 or x_s2 > A_0:
+ x_s2 = None
+ """
+ Determine the actual new position of the considered particle from the both solver solutions.
+ Must be the closest position to the original position of the considered particle
+ """
+ # Coose the new position based on the validity of values
+ if x_s1 == None and x_s2 != None:
+ x = x_s2
+ elif x_s2 == None and x_s1 != None:
+ x = x_s1
+ elif x_s1 == None and x_s2 == None:
+ error_position_update1w = 3
+ return y_i, x_i, error_position_update1w
+ else:
+ # Vectors between the centers of the considered particle's original position and the new positions
+ V_centers_new1 = np.array([z_i, y_i, x_i]) - np.array([z_i, y_s, x_s1])
+ V_centers_new2 = np.array([z_i, y_i, x_i]) - np.array([z_i, y_s, x_s2])
+ # Distances between these centers
+ dist_centers_new1 = np.linalg.norm(V_centers_new1)
+ dist_centers_new2 = np.linalg.norm(V_centers_new2)
+ # Choose new position based on the shortest distance
+ if dist_centers_new1 <= dist_centers_new2:
+ x = x_s1
+ else:
+ x = x_s2
+ # Assign x-coordinate
+ y = y_s
+ """
+ Check wheter the new positions' coordinates are located on the intersection plane
+ Just double checking wheter the solver found correct solutions
+ The equation of the intersection plane is Az+By+Cx+D=0, where A, B, C and D are the equation coefficients.
+ """
+ # Determine the plane equation coefficients
+ # The normal vector to the plane is the same as the V_centers_12
+ coeff_A, coeff_B, coeff_C = V_unit_centers_1w
+ # Calculate D using the point on the plane
+ coeff_D = -(coeff_A * c_ic[0] + coeff_B * c_ic[1] + coeff_C * c_ic[2])
+ # Check whether the new positions of considered particle is located on the intersection plane
+ on_plane = check_point_on_plane(z_i, y, x, coeff_A, coeff_B, coeff_C, coeff_D)
+ # Error value is 2 if the position of the first solver solution is not located on the intersection plane
+ if on_plane == False:
+ error_position_update1w = 2
+ return y_i, x_i, error_position_update1w
+
+ return y, x, error_position_update1w
+
+def solver_position1w(params, c1, r1, z_i, par, walls):
+ """ Solver to calculate the new position of the considered particle if it overlaps with one particle and one wall """
+ if walls == 1 or walls == -1:
+ # Solver variables
+ y = params
+ # Equation of a circle. The circles coresponds here to the auxiliary spheres
+ eq1 = (par - c1[2])**2 + (y - c1[1])**2 + (z_i - c1[0])**2 - r1**2
+ return eq1
+ elif walls == 2 or walls == -2:
+ # Solver variables
+ x = params
+ # Equation of a circle. The circles coresponds here to the auxiliary spheres
+ eq1 = (x - c1[2])**2 + (par - c1[1])**2 + (z_i - c1[0])**2 - r1**2
+ return eq1
+
+def position_update2w(z_i, y_i, x_i, d_i, M_particle_steady, overlaps, walls):
+ """
+ Function to update the position of the considered particle if it overlaps with two particles and one wall.
+ There are two mathematical solutions of the new position, but only one logical.
+ A solver is used to determine these solutions.
+ First, a triangle plane is calculated which crosses the centers of the two steady state particles and the wall.
+ Moreover, the centroid of the triangle is determined.
+ The triangle plane is used to identify the guesses of the both solutions of the considered particle's new position.
+ It is tested whether the both solutions differ and if they are correct.
+ To identify the logical solution of the both solutions, the z-coordinates of the both soultions and the centroid are compared.
+ The new position with the z-coordinate above the centroid is the logical solution.
+ """
+ # Used in case an error occurs during position update.
+ error_position_update2w = 0
+ # Coordinates of the first particle in steady state
+ z1 = float(M_particle_steady[overlaps[0], 1])
+ y1 = float(M_particle_steady[overlaps[0], 2])
+ x1 = float(M_particle_steady[overlaps[0], 3])
+ # Diameter of the first particle in steady state
+ d1 = float(M_particle_steady[overlaps[0], 5])
+ # Center coordinates of the first auxiliary sphere
+ c1 = np.array([z1, y1, x1])
+ # Radius of the first auxiliary sphere
+ r1 = (d_i + d1) / 2
+ # Coordinates of the second particle in steady state
+ z2 = float(M_particle_steady[overlaps[1], 1])
+ y2 = float(M_particle_steady[overlaps[1], 2])
+ x2 = float(M_particle_steady[overlaps[1], 3])
+ # Diameter of the second particle in steady state
+ d2 = float(M_particle_steady[overlaps[1], 5])
+ # Center coordinates of the second auxiliary sphere
+ c2 = np.array([z2, y2, x2])
+ # Radius of the second auxiliary sphere
+ r2 = (d_i + d2) / 2
+ # Coordinates of the wall
+ z_w1 = (z1+z2)/2
+ if walls == 1:
+ x_w1 = A_0
+ y_w1 = (y1+y2)/2
+ elif walls == -1:
+ x_w1 = -A_0
+ y_w1 = (y1+y2)/2
+ elif walls == 2:
+ x_w1 = (x1+x2)/2
+ y_w1 = A_0
+ elif walls == -2:
+ x_w1 = (x1+x2)/2
+ y_w1 = -A_0
+ """ Calculation of the triangle's plane """
+ # Calculate vectors AB and AC
+ AB = np.array([z2 - z1, y2 - y1, x2 - x1])
+ AC = np.array([z_w1 - z1, y_w1 - y1, x_w1 - x1])
+ # Compute the normal vector to the plane (cross product of AB and AC)
+ V_triangle_norm = np.cross(AB, AC)
+ A, B, C = V_triangle_norm
+ # Calculate D using the plane equation Ax + By + Cz + D = 0
+ # Using point A (x1, y1, z1) to find D
+ D = -(A * z1 + B * y1 + C * x1)
+ """
+ Calculation of the reflected position of the considered particle.
+ First, the distance between the center of the considered particle and the triangle's plane is determined.
+ Then, the position of the considered particle is reflected about the triangle's plane
+ The coordiantes of the reflected position are used as second guess for the solver to determine the second solver solution.
+ """
+ # Calculate the distance from the considered particle's center to the plane
+ distance_pi = (A * z_i + B * y_i + C * x_i + D) / np.linalg.norm(V_triangle_norm)
+ # Calculate the reflection point on the plane
+ P_reflected = np.array([z_i, y_i, x_i]) - 2 * distance_pi * V_triangle_norm / np.linalg.norm(V_triangle_norm)
+ """ Solver to find possible new position of the considered particle """
+ # Coordinates for the initial guess
+ x_01 = x_i
+ y_01 = y_i
+ z_01 = z_i
+ x_02 = P_reflected[2]
+ y_02 = P_reflected[1]
+ z_02 = P_reflected[0]
+ # If contact is with X-wall
+ if walls == 1 or walls == -1:
+ # Initial guess values
+ initial_guess1 = [y_01, z_01]
+ initial_guess2 = [y_02, z_02]
+ # X-coordinate of the new position
+ if walls == 1:
+ x_s1 = A_0 - d_i / 2
+ x_s2 = A_0 - d_i / 2
+ elif walls == -1:
+ x_s1 = -A_0 + d_i / 2
+ x_s2 = -A_0 + d_i / 2
+ # Solver
+ y_s1, z_s1 = fsolve(solver_position2w, initial_guess1, args=(c1, c2, r1, r2, x_s1, walls))
+ y_s2, z_s2 = fsolve(solver_position2w, initial_guess2, args=(c1, c2, r1, r2, x_s2, walls))
+ """ Check wheter the solver found two different solutions or the same twice """
+ # Error value is 1 if the solver found the same position twice
+ if y_s1 == y_s2 and z_s1 == z_s2:
+ error_position_update2w = 1
+ # Exclude values behind the simulation walls
+ if y_s1 < -A_0 or y_s1 > A_0:
+ y_s1 = None
+ if y_s2 < -A_0 or y_s2 > A_0:
+ y_s2 = None
+ # If contact is with Y-wall
+ if walls == 2 or walls == -2:
+ # Initial guess values
+ initial_guess1 = [x_01, z_01]
+ initial_guess2 = [x_02, z_02]
+ # X-coordinate of the new position
+ if walls == 2:
+ y_s1 = A_0 - d_i / 2
+ y_s2 = A_0 - d_i / 2
+ elif walls == -2:
+ y_s1 = -A_0 + d_i / 2
+ y_s2 = -A_0 + d_i / 2
+ # Solver
+ x_s1, z_s1 = fsolve(solver_position2w, initial_guess1, args=(c1, c2, r1, r2, y_s1, walls))
+ x_s2, z_s2 = fsolve(solver_position2w, initial_guess2, args=(c1, c2, r1, r2, y_s2, walls))
+ """ Check wheter the solver found two different solutions or the same twice """
+ # Error value is 1 if the solver found the same position twice
+ if x_s1 == x_s2 and z_s1 == z_s2:
+ error_position_update2w = 1
+ # Exclude values behind the simulation walls
+ if x_s1 < -A_0 or x_s1 > A_0:
+ y_s1 = None
+ if x_s2 < -A_0 or x_s2 > A_0:
+ x_s2 = None
+ """
+ Determine the actual new position of the considered particle from the both solver solutions.
+ Must be the closest position to the original position of the considered particle
+ """
+ # Z-coordinate of triangles's centroid
+ z_cen = (z1 + z2 + z_w1) / 3
+ # If the center of the considered particle is located above the centroid,
+ # the coordinates of the first solver solution are assigned to the center of the considered particle
+ if x_s1 != None and y_s1 != None and z_s1 > z_cen:
+ z = z_s1
+ y = y_s1
+ x = x_s1
+ # If the center of the considered particle is located above the centroid,
+ # the coordinates of the second solver solution are assigned to the center of the considered particle
+ elif x_s2 != None and y_s2 != None and z_s2 > z_cen:
+ z = z_s2
+ y = y_s2
+ x = x_s2
+ else:
+ error_position_update2w = 2
+ return z_i, y_i, x_i, error_position_update2w
+ """
+ Check wheter the distances between the centers of the three steady state particles
+ and the new center position of the considered particle are in the samerange as the sums of their radii.
+ """
+ same_range = check_distances2w(r1, r2, c1, c2, d_i, z, y, x, walls)
+ # Error value is 2 if the distance between the steady state particles and the position of the solver solution
+ # is not in the same range as the sums of their radii
+ if same_range == False:
+ error_position_update2w = 3
+ return z_i, y_i, x_i, error_position_update2w
+
+ return z, y, x, error_position_update2w
+
+def solver_position2w(params, c1, c2, r1, r2, par, walls):
+ """ Solver to calculate the new position of the considered particle if it overlaps with two particles and one wall """
+ if walls == 1 or walls == -1:
+ # Solver variables
+ y, z = params
+ # Equation of a circle. The circles coresponds here to the auxiliary spheres
+ eq1 = (par - c1[2])**2 + (y - c1[1])**2 + (z - c1[0])**2 - r1**2
+ eq2 = (par - c2[2])**2 + (y - c2[1])**2 + (z - c2[0])**2 - r2**2
+ return [eq1, eq2]
+ elif walls == 2 or walls == -2:
+ # Solver variables
+ x, z = params
+ # Equation of a circle. The circles coresponds here to the auxiliary spheres
+ eq1 = (x - c1[2])**2 + (par - c1[1])**2 + (z - c1[0])**2 - r1**2
+ eq2 = (x - c2[2])**2 + (par - c2[1])**2 + (z - c2[0])**2 - r2**2
+ return [eq1, eq2]
+
+def check_distances2w(r1, r2, c1, c2, d_i, z_s, y_s, x_s, walls):
+ """ Function that checks wheter the distances between the centers of the two steady state particles, the wall,
+ and the new center position of the considered particle are in the same range as the sums of their radii. """
+ # True if distances are in the same range as the sums of radii
+ same_range = True
+ # Vectors between the centers of the steady state particles and the new center positions of the considered particle
+ V_centers_1i = np.array([z_s, y_s, x_s]) - c1
+ V_centers_2i = np.array([z_s, y_s, x_s]) - c2
+ # Distance between the centers of the steady state particles and the new center position of the considered particle
+ dist_centers_1i = np.linalg.norm(V_centers_1i)
+ dist_centers_2i = np.linalg.norm(V_centers_2i)
+ # Distance between the center of the considered particle and the wall
+ if walls == 1 or walls == -1:
+ dist_wi = A_0 - abs(x_s)
+ elif walls == 2 or walls == -2:
+ dist_wi = A_0 - abs(y_s)
+ # Check if the left-hand side is approximately zero (account for floating-point precision)
+ if abs(dist_centers_1i-r1) > 1e-9 or abs(dist_centers_2i-r2) > 1e-9 or abs(dist_wi-d_i/2) > 1e-9:
+ same_range = False
+
+ return same_range
+
+def check_steadystate2w(z_i, y_i, x_i, particle_i, M_particle_steady, overlaps, walls):
+ """
+ Function to check the steady state.
+ As soon as considered particle collides with two steady state particles and one wall, a steady state check is performed.
+ Considered particle is in steady state if its center lies within the y-x-plane of a triangle located between the
+ steady state particle centers and the wall.
+ Barycentric coordinate method is used to determine if the center of the considered particle lies within the polygon.
+ Moreover, barycentric coordinates are used to determine the subdomain of the triangle where the center lies.
+ Based on the subdomain, the particle rolling direction is identified.
+ """
+ # Used in case an error occurs during steadystate check.
+ error_check_steadystate3 = 0
+ # Variable for steady state check condition
+ steady_check = 0
+ # Coordinates of the first particle in steady state
+ z1 = float(M_particle_steady[overlaps[0], 1])
+ y1 = float(M_particle_steady[overlaps[0], 2])
+ x1 = float(M_particle_steady[overlaps[0], 3])
+ # Coordinates of the second particle in steady state
+ z2 = float(M_particle_steady[overlaps[1], 1])
+ y2 = float(M_particle_steady[overlaps[1], 2])
+ x2 = float(M_particle_steady[overlaps[1], 3])
+ # Coordinates of the wall
+ z_w1 = z_i
+ if walls == 1:
+ x_w1 = A_0
+ y_w1 = y_i
+ elif walls == -1:
+ x_w1 = -A_0
+ y_w1 = y_i
+ elif walls == 2:
+ x_w1 = x_i
+ y_w1 = A_0
+ elif walls == -2:
+ x_w1 = x_i
+ y_w1 = -A_0
+ # Z-coordinate of triangle's centroid
+ z_cen = (z1 + z2 + z_w1) / 3
+ # Check whether the center of the considered particle is located above the centroid
+ # Double checking at this point, since already handled in the position_update3 function.
+ # Error value is 1 when the center is located below the centroid
+ if z_i < z_cen:
+ error_check_steadystate3 = 1
+ return steady_check, overlaps, error_check_steadystate3
+ """ Barycentric coordinate technique """
+ # Denominator
+ denominator = ((y2 - y_w1) * (x1 - x_w1) + (x_w1 - x2) * (y1 - y_w1))
+ # Check whether the triangle points are collinear
+ # Error value is 2 when the triangle points are collinear
+ if denominator == 0:
+ error_check_steadystate3 = 2
+ return steady_check, overlaps, error_check_steadystate3
+ # Barycentric coordinates
+ a = ((y2 - y_w1) * (x_i - x_w1) + (x_w1 - x2) * (y_i - y_w1)) / denominator
+ b = ((y_w1 - y1) * (x_i - x_w1) + (x1 - x_w1) * (y_i - y_w1)) / denominator
+ c = 1 - a - b
+ # Determine the subdomain of the triangle and thus wheter it is in steady state or in which direction it is rolling
+ if (0 <= a <= 1) and (0 <= b <= 1) and (0 <= c <= 1):
+ steady_check = 3
+ elif a < 0 and b >= 0 and c >= 0:
+ # Rolling along steady state particle Nr.2 and wall
+ overlaps = np.delete(overlaps, 0)
+ steady_check = 4
+ elif b < 0 and a >= 0 and c >= 0:
+ # Rolling along steady state particle Nr.1 and wall
+ overlaps = np.delete(overlaps, 1)
+ steady_check = 4
+ elif c < 0 and a >= 0 and b >= 0:
+ # Rolling along steady state particle Nr.1 and Nr.2
+ steady_check = 2
+ elif a < 0 and c < 0 and b > 1:
+ # Rolling along steady state particle Nr.2
+ overlaps = np.delete(overlaps, 0)
+ steady_check = 1
+ elif b < 0 and c < 0 and a > 1:
+ # Rolling along steady state particle Nr.1
+ overlaps = np.delete(overlaps, 1)
+ steady_check = 1
+ else:
+ # Error value is 3 when the subdomain of the triangle could not be determined correctly
+ error_check_steadystate3 = 3
+
+ return steady_check, overlaps, error_check_steadystate3
+
+def position_update1w2(z_i, y_i, x_i, d_i, M_particle_steady, overlaps, walls):
+ """
+ Function to update the position of the considered particle if it overlaps with one particle and two walls.
+ There are two mathematical solutions of the new position, but only one logical.
+ A solver is used to determine these solutions.
+ First, a triangle plane is calculated which crosses the center of the steady state particle and the walls.
+ Moreover, the centroid of the triangle is determined.
+ The triangle plane is used to identify the guesses of the both solutions of the considered particle's new position.
+ It is tested whether the both solutions differ and if they are correct.
+ To identify the logical solution of the both solutions, the z-coordinates of the both soultions and the centroid are compared.
+ The new position with the z-coordinate above the centroid is the logical solution.
+ """
+ # Used in case an error occurs during position update.
+ error_position_update1w2 = 0
+ # Coordinates of the first particle in steady state
+ z1 = float(M_particle_steady[overlaps[0], 1])
+ y1 = float(M_particle_steady[overlaps[0], 2])
+ x1 = float(M_particle_steady[overlaps[0], 3])
+ # Diameter of the first particle in steady state
+ d1 = float(M_particle_steady[overlaps[0], 5])
+ # Center coordinates of the first auxiliary sphere
+ c1 = np.array([z1, y1, x1])
+ # Radius of the first auxiliary sphere
+ r1 = (d_i + d1) / 2
+ # Coordinates of the wall
+ z_w1 = z1
+ z_w2 = z1
+ if walls == 3:
+ x_w1 = A_0
+ y_w1 = y_i
+ x_w2 = x_i
+ y_w2 = A_0
+ elif walls == -3:
+ x_w1 = -A_0
+ y_w1 = y_i
+ x_w2 = x_i
+ y_w2 = A_0
+ elif walls == 4:
+ x_w1 = A_0
+ y_w1 = y_i
+ x_w2 = x_i
+ y_w2 = -A_0
+ elif walls == -4:
+ x_w1 = -A_0
+ y_w1 = y_i
+ x_w2 = x_i
+ y_w2 = -A_0
+ """ Calculation of the triangle's plane """
+ # Calculate vectors AB and AC
+ AB = np.array([z_w2 - z1, y_w2 - y1, x_w2 - x1])
+ AC = np.array([z_w1 - z1, y_w1 - y1, x_w1 - x1])
+ # Compute the normal vector to the plane (cross product of AB and AC)
+ V_triangle_norm = np.cross(AB, AC)
+ A, B, C = V_triangle_norm
+ # Calculate D using the plane equation Ax + By + Cz + D = 0
+ # Using point A (x1, y1, z1) to find D
+ D = -(A * z1 + B * y1 + C * x1)
+ """
+ Calculation of the reflected position of the considered particle.
+ First, the distance between the center of the considered particle and the triangle's plane is determined.
+ Then, the position of the considered particle is reflected about the triangle's plane
+ The coordiantes of the reflected position are used as second guess for the solver to determine the second solver solution.
+ """
+ # Calculate the distance from the considered particle's center to the plane
+ distance_pi = (A * z_i + B * y_i + C * x_i + D) / np.linalg.norm(V_triangle_norm)
+ # Calculate the reflection point on the plane
+ P_reflected = np.array([z_i, y_i, x_i]) - 2 * distance_pi * V_triangle_norm / np.linalg.norm(V_triangle_norm)
+ """ Solver to find possible new position of the considered particle """
+ # Coordinates for the initial guess
+ z_01 = z_i
+ z_02 = P_reflected[0]
+ # Initial guess values
+ initial_guess1 = z_01
+ initial_guess2 = z_02
+ # X and y-coordinates of the new position
+ if walls == 3:
+ x = A_0 - d_i / 2
+ y = A_0 - d_i / 2
+ elif walls == -3:
+ x = -A_0 + d_i / 2
+ y = A_0 - d_i / 2
+ elif walls == 4:
+ x = A_0 - d_i / 2
+ y = -A_0 + d_i / 2
+ elif walls == -4:
+ x = -A_0 + d_i / 2
+ y = -A_0 + d_i / 2
+ # Solver
+ z_s1 = fsolve(solver_position1w2, initial_guess1, args=(c1, r1, x, y))
+ z_s2 = fsolve(solver_position1w2, initial_guess2, args=(c1, r1, x, y))
+ z_s1 = z_s1[0]
+ z_s2 = z_s2[0]
+ """ Check wheter the solver found two different solutions or the same twice """
+ # Error value is 1 if the solver found the same position twice
+ if z_s1 == z_s2:
+ error_position_update1w2 = 1
+ """
+ Determine the actual new position of the considered particle from the both solver solutions.
+ Must be the closest position to the original position of the considered particle
+ """
+ # Z-coordinate of polygon's centroid
+ z_cen = z1
+ # If the center of the considered particle is located above the centroid,
+ # the coordinates of the first solver solution are assigned to the center of the considered particle
+ if z_s1 > z_cen:
+ z = z_s1
+ # If the center of the considered particle is located above the centroid,
+ # the coordinates of the second solver solution are assigned to the center of the considered particle
+ elif z_s2 > z_cen:
+ z = z_s2
+ else:
+ error_position_update1w2 = 2
+ return z_i, y_i, x_i, error_position_update1w2
+ """
+ Check wheter the distances between the centers of the three steady state particles
+ and the new center position of the considered particle are in the samerange as the sums of their radii.
+ """
+ same_range = check_distances1w2(r1, c1, d_i, z, y, x)
+ # Error value is 2 if the distance between the steady state particles and the position of the solver solution
+ # is not in the same range as the sums of their radii
+ if same_range == False:
+ error_position_update1w2 = 3
+ return z_i, y_i, x_i, error_position_update1w2
+
+ return z, y, x, error_position_update1w2
+
+def solver_position1w2(params, c1, r1, x_s, y_s):
+ """ Solver to calculate the new position of the considered particle if it overlaps one particle and two walls """
+ # Solver variables
+ z = params
+ # Equation of a circle. The circles coresponds here to the auxiliary spheres
+ eq1 = (x_s - c1[2])**2 + (y_s - c1[1])**2 + (z - c1[0])**2 - r1**2
+ return eq1
+
+def check_distances1w2(r1, c1, d_i, z_s, y_s, x_s):
+ """ Function that checks wheter the distances between the center of one steady state particles, the walls,
+ and the new center position of the considered particle are in the same range as the sums of their radii. """
+ # True if distances are in the same range as the sums of radii
+ same_range = True
+ # Vectors between the centers of the steady state particles and the new center positions of the considered particle
+ V_centers_1i = np.array([z_s, y_s, x_s]) - c1
+ # Distance between the centers of the steady state particles and the new center position of the considered particle
+ dist_centers_1i = np.linalg.norm(V_centers_1i)
+ # Distances between the center of the considered particle and the walls
+ dist_w1i = A_0 - abs(x_s)
+ dist_w2i = A_0 - abs(y_s)
+ # Check if the left-hand side is approximately zero (account for floating-point precision)
+ if abs(dist_centers_1i-r1) > 1e-9 or abs(dist_w2i-d_i/2) > 1e-9 or abs(dist_w1i-d_i/2) > 1e-9:
+ same_range = False
+
+ return same_range
+
+def check_steadystate1w2(z_i, y_i, x_i, M_particle_steady, overlaps, walls):
+ """
+ Function to check the steady state.
+ As soon as considered particle collides with one steady state particle and two walls, a steady state check is performed.
+ Considered particle is in steady state if its center lies within the y-x-plane of a triangle located between the
+ steady state particle center and the walls.
+ Barycentric coordinate method is used to determine if the center of the considered particle lies within the polygon.
+ Moreover, barycentric coordinates are used to determine the subdomain of the triangle where the center lies.
+ Based on the subdomain, the particle rolling direction is identified.
+ """
+ # Used in case an error occurs during steadystate check.
+ error_check_steadystate1w2 = 0
+ # Variable for steady state check condition
+ steady_check = 0
+ # Coordinates of the first particle in steady state
+ z1 = float(M_particle_steady[overlaps[0], 1])
+ y1 = float(M_particle_steady[overlaps[0], 2])
+ x1 = float(M_particle_steady[overlaps[0], 3])
+ # Coordinates of the wall
+ if walls == 3:
+ x_w1 = A_0
+ y_w1 = y_i
+ x_w2 = x_i
+ y_w2 = A_0
+ elif walls == -3:
+ x_w1 = -A_0
+ y_w1 = y_i
+ x_w2 = x_i
+ y_w2 = A_0
+ elif walls == 4:
+ x_w1 = A_0
+ y_w1 = y_i
+ x_w2 = x_i
+ y_w2 = -A_0
+ elif walls == -4:
+ x_w1 = -A_0
+ y_w1 = y_i
+ x_w2 = x_i
+ y_w2 = -A_0
+ # Z-coordinate of triangle's centroid
+ z_cen = z1
+ # Check whether the center of the considered particle is located above the centroid
+ # Double checking at this point, since already handled in the position_update3 function.
+ # Error value is 1 when the center is located below the centroid
+ if z_i < z_cen:
+ error_check_steadystate1w2 = 1
+ return steady_check, overlaps, error_check_steadystate1w2
+ """ Barycentric coordinate technique """
+ # Denominator
+ denominator = ((y_w2 - y_w1) * (x1 - x_w1) + (x_w1 - x_w2) * (y1 - y_w1))
+ # Check whether the triangle points are collinear
+ # Error value is 2 when the triangle points are collinear
+ if denominator == 0:
+ error_check_steadystate1w2 = 2
+ return steady_check, overlaps, error_check_steadystate1w2
+ # Barycentric coordinates
+ a = ((y_w2 - y_w1) * (x_i - x_w1) + (x_w1 - x_w2) * (y_i - y_w1)) / denominator
+ b = ((y_w1 - y1) * (x_i - x_w1) + (x1 - x_w1) * (y_i - y_w1)) / denominator
+ c = 1 - a - b
+ # Determine the subdomain of the triangle and thus wheter it is in steady state or in which direction it is rolling
+ if (0 <= a <= 1) and (0 <= b <= 1) and (0 <= c <= 1):
+ steady_check = 3
+ elif b < 0 and a >= 0 and c >= 0:
+ # Rolling along steady state particle Nr.1 and wall Nr.2
+ if walls == 3 or walls == -3:
+ walls = 2
+ elif walls == 4 or walls == -4:
+ walls = -2
+ steady_check = 4
+ elif c < 0 and a >= 0 and b >= 0:
+ # Rolling along steady state particle Nr.1 and wall Nr.1
+ if walls == 3 or walls == 4:
+ walls = 1
+ elif walls == -3 or walls == -4:
+ walls = -1
+ steady_check = 4
+ elif b < 0 and c < 0 and a > 1:
+ # Rolling along steady state particle Nr.1
+ steady_check = 1
+ else:
+ # Error value is 3 when the subdomain of the triangle could not be determined correctly
+ error_check_steadystate1w2 = 3
+
+ return steady_check, overlaps, walls, error_check_steadystate1w2
+
+def position_update2w2(z_i, y_i, x_i, d_i, M_particle_steady, overlaps, walls):
+ """
+ Function to update the position of the considered particle if it overlaps with two particles and two walls.
+ There are two mathematical solutions of the new position, but only one logical.
+ A solver is used to determine these solutions.
+ First, a polygon plane is calculated which crosses the centers of the two steady state particles and the walls.
+ Moreover, the centroid of the polygon is determined.
+ The polygon plane is used to identify the guesses of the both solutions of the considered particle's new position.
+ It is tested whether the both solutions differ and if they are correct.
+ To identify the logical solution of the both solutions, the z-coordinates of the both soultions and the centroid are compared.
+ The new position with the z-coordinate above the centroid is the logical solution.
+ """
+ # Used in case an error occurs during position update.
+ error_position_update2w2 = 0
+ # Coordinates of the first particle in steady state
+ z1 = float(M_particle_steady[overlaps[0], 1])
+ y1 = float(M_particle_steady[overlaps[0], 2])
+ x1 = float(M_particle_steady[overlaps[0], 3])
+ # Diameter of the first particle in steady state
+ d1 = float(M_particle_steady[overlaps[0], 5])
+ # Radius of the first auxiliary sphere
+ r1 = (d_i + d1) / 2
+ # Center coordinates of the first auxiliary sphere
+ c1 = np.array([z1, y1, x1])
+ # Coordinates of the second particle in steady state
+ z2 = float(M_particle_steady[overlaps[1], 1])
+ y2 = float(M_particle_steady[overlaps[1], 2])
+ x2 = float(M_particle_steady[overlaps[1], 3])
+ # Diameter of the second particle in steady state
+ d2 = float(M_particle_steady[overlaps[1], 5])
+ # Radius of the first auxiliary sphere
+ r2 = (d_i + d2) / 2
+ # Center coordinates of the second auxiliary sphere
+ c2 = np.array([z2, y2, x2])
+ # Coordinates of the wall
+ z_w1 = (z1+z2)/2
+ z_w2 = (z1+z2)/2
+ """ Solver to find possible new position of the considered particle """
+ # Array with the center coordinates and radii of the steady state particles
+ particles = [
+ {'center': c1, 'diameter': d1},
+ {'center': c2, 'diameter': d2}]
+ # Coordinates for the initial guess
+ x_01 = x_i
+ y_01 = y_i
+ z_01 = z_i
+ # Initial guess values
+ initial_guess = [x_01, y_01, z_01]
+ # Bounds
+ bnds = [(None, None), (-A_0, A_0), (-A_0, A_0)]
+ # Perform the optimization
+ result = minimize(solver_position2w2, initial_guess, args=(particles, d_i), method='L-BFGS-B', bounds=bnds, tol=1e-10)
+ # Extract the coordinates
+ z_s, y_s, x_s = result.x[:3]
+ """
+ Check wheter the distances between the centers of the four steady state particles
+ and the new center position of the considered particle are in the samerange as the sums of their radii.
+ """
+ # Check if the left-hand side is approximately zero (account for floating-point precision)
+ same_range = check_distances2w2(r1, r2, c1, c2, d_i, z_s, y_s, x_s)
+ # Error value is 1 if the distance between the steay state particles and the position of the first solver solution
+ # is not in the same range as the sums of their radii
+ if same_range == False:
+ error_position_update2w2 = 1
+ return z_i, y_i, x_i, error_position_update2w2
+ """
+ Check whether the new position of the considered particles is located above the centroid of the steady state particles
+ """
+ # Z-coordinate of polygon's centroid
+ z_cen = (z1 + z2 + z_w1 + z_w2) / 4
+ # If the center of the considered particle is located above the centroid,
+ # the new coordinates are assigned to the center of the considered particle
+ if z_s > z_cen:
+ z = z_s
+ y = y_s
+ x = x_s
+ # If the center of the considered particle is located below the centroid,
+ # the new coordinates are not correct
+ else:
+ error_position_update2w2 = 2
+ return z_i, y_i, x_i, error_position_update2w2
+
+ return z, y, x, error_position_update2w2
+
+def solver_position2w2(params, particles, d_i):
+ """ Solver to calculate the new position of the considered particle if it overlaps with two particles and two walls """
+ # Center coordinates of the considered particle are the optimizing parameters
+ c_i = np.array(params[:3])
+ # Mini
+ total_error = 0
+ # Calculate the total error between the actual and required distances of the particle centers
+ for particle in particles:
+ # Center of the steady state particle
+ c_j = particle['center']
+ # Diameter of the steady state particle
+ d_j = particle['diameter']
+ # Actual distance between the centers of the considered particle and the steady state particle
+ distance = np.linalg.norm(c_i - c_j)
+ # Required distance between the centers of the considered particle and the steady state particle
+ required_distance = (d_i + d_j)/2
+ # Minimization function
+ total_error += 100*(distance - required_distance)/required_distance
+ # Actual distance between the center of the considered particle and the wall
+ distance_w1 = A_0 - abs(c_i[2])
+ distance_w2 = A_0 - abs(c_i[1])
+ # Required distance between the center of the considered particle and the steady state particle
+ required_distance_w = d_i/2
+ # Update minimization function
+ total_error += 100*abs(distance_w1 - required_distance_w)/required_distance_w
+ total_error += 100*abs(distance_w2 - required_distance_w)/required_distance_w
+
+ return total_error
+
+def check_distances2w2(r1, r2, c1, c2, d_i, z_s, y_s, x_s):
+ """ Function that checks wheter the distances between the centers of the two steady state particles, the walls,
+ and the new center position of the considered particle are in the same range as the sums of their radii. """
+ # True if distances are in the same range as the sums of radii
+ same_range = True
+ # Vectors between the centers of the steady state particles and the new center positions of the considered particle
+ V_centers_1i = np.array([z_s, y_s, x_s]) - c1
+ V_centers_2i = np.array([z_s, y_s, x_s]) - c2
+ # Distance between the centers of the steady state particles and the new center position of the considered particle
+ dist_centers_1i = np.linalg.norm(V_centers_1i)
+ dist_centers_2i = np.linalg.norm(V_centers_2i)
+ # Distance between the center of the considered particle and the wall
+ dist_w1i = A_0 - abs(x_s)
+ dist_w2i = A_0 - abs(y_s)
+ # Check if the left-hand side is approximately zero (account for floating-point precision)
+ if abs(dist_centers_1i-r1) > 1e-9 or abs(dist_centers_2i-r2) > 1e-9 or abs(dist_w1i-d_i/2) > 1e-9 or abs(dist_w2i-d_i/2) > 1e-9:
+ same_range = False
+
+ return same_range
+
+def check_steadystate2w2(z_i, y_i, x_i, M_particle_steady, overlaps, walls):
+ """
+ Function to check the steady state.
+ As soon as considered particle collides with two steady state particles and two walls, a steady state check is performed.
+ Considered particle is in steady state if its center lies within the y-x-plane of a polygon located between the
+ steady state particle centers and the walls.
+ Barycentric coordinate method is used to determine if the center of the condsidered particle lies within the polygon.
+ For this purpose, the polygon is devided into two triangles.
+ Moreover, barycentric coordinates are used to determine the subdomain of the polygon where the point lies.
+ Based on the subdomain, the particle rolling direction is identified.
+ """
+ # Used in case an error occurs during steadystate check.
+ error_check_steadystate2w2 = 0
+ # Variable for steady state check condition
+ steady_check = 0
+ # Coordinates of the first particle in steady state
+ z1 = float(M_particle_steady[overlaps[0], 1])
+ y1 = float(M_particle_steady[overlaps[0], 2])
+ x1 = float(M_particle_steady[overlaps[0], 3])
+ # Coordinates of the second particle in steady state
+ z2 = float(M_particle_steady[overlaps[1], 1])
+ y2 = float(M_particle_steady[overlaps[1], 2])
+ x2 = float(M_particle_steady[overlaps[1], 3])
+ # Coordinates of the wall
+ z_w1 = (z1+z2)/2
+ z_w2 = (z1+z2)/2
+ if walls == 3:
+ x_w1 = A_0
+ y_w1 = y_i
+ x_w2 = x_i
+ y_w2 = A_0
+ elif walls == -3:
+ x_w1 = -A_0
+ y_w1 = y_i
+ x_w2 = x_i
+ y_w2 = A_0
+ elif walls == 4:
+ x_w1 = A_0
+ y_w1 = y_i
+ x_w2 = x_i
+ y_w2 = -A_0
+ elif walls == -4:
+ x_w1 = -A_0
+ y_w1 = y_i
+ x_w2 = x_i
+ y_w2 = -A_0
+ # Z-coordinate of polygon's centroid
+ z_cen = (z1 + z2 + z_w1 + z_w2) / 4
+ # Check whether the center of the considered particle is located above the centroid
+ # Double checking at this point, since already handled in the position_update4 function.
+ # Error value is 1 when the center is located below the centroid
+ if z_i < z_cen:
+ error_check_steadystate2w2 = 1
+ return steady_check, overlaps, error_check_steadystate2w2
+ """
+ The polygon is devided into two triangles.
+ The first triangle is defined by the centers of the steady state particles 1, 2 and 4. 1=A, 2=B, 4=D
+ The second triangle is defined by the centers of the steady state particles 2, 3 and 4. 2=B, 3=C, 4=D
+ Barycentric coordinate technique is applied on each triangle
+ """
+ # Denominator of the first triangle
+ denominator1 = ((y2 - y_w1) * (x1 - x_w1) + (x_w1 - x2) * (y1 - y_w1))
+ # Denominator of the second triangle
+ denominator2 = ((y_w2 - y_w1) * (x2 - x_w1) + (x_w1 - x_w2) * (y2 - y_w1))
+ # Check whether the triangle points are collinear
+ # Error value is 2 when the triangle points are collinear
+ if denominator1 == 0 or denominator2 == 0:
+ error_check_steadystate2w2 = 2
+ return steady_check, overlaps, error_check_steadystate2w2
+ # Barycentric coordinates of the first triangle
+ a = ((y2 - y_w1) * (x_i - x_w1) + (x_w1 - x2) * (y_i - y_w1)) / denominator1
+ b1 = ((y_w1 - y1) * (x_i - x_w1) + (x1 - x_w1) * (y_i - y_w1)) / denominator1
+ d1 = 1 - a - b1
+ # Barycentric coordinates of the second triangle
+ b2 = ((y_w2 - y_w1) * (x_i - x_w1) + (x_w1 - x_w2) * (y_i - y_w1)) / denominator2
+ c = ((y_w1 - y2) * (x_i - x_w1) + (x2 - x_w1) * (y_i - y_w1)) / denominator2
+ d2 = 1 - b2 - c
+ # Determine the subdomain of the polygon and thus wheter it is in steady state or in which direction it is rolling
+ if (0 <= a <= 1) and (0 <= b1 <= 1) and (0 <= d1 <= 1):
+ steady_check = 3
+ elif (0 <= b2 <= 1) and (0 <= c <= 1) and (0 <= d2 <= 1):
+ steady_check = 3
+ elif b1 < 0 and d1 < 0:
+ # Rolling along steady state particle Nr.1
+ overlaps = np.delete(overlaps, 1)
+ steady_check = 1
+ elif d1 < 0 and d2 < 0:
+ # Rolling along steady state particle Nr.2
+ overlaps = np.delete(overlaps, 0)
+ steady_check = 1
+ elif d1 < 0 and b1 >= 0 and d2 >= 0:
+ # Rolling along steady state particle Nr.1 and Nr.2
+ steady_check = 2
+ elif d2 < 0 and b2 >= 0 and d1 >= 0:
+ # Rolling along steady state particle Nr.2 and wall Nr.1
+ if walls == 3 or walls == 4:
+ walls = 1
+ elif walls == -3 or walls == -4:
+ walls = -1
+ overlaps = np.delete(overlaps, 0)
+ steady_check = 4
+ elif b1 < 0 and d1 >= 0 and b2 >= 0:
+ # Rolling along steady state particle Nr.1 and wall Nr.2
+ if walls == 3 or walls == -3:
+ walls = 2
+ elif walls == 4 or walls == -4:
+ walls = -2
+ overlaps = np.delete(overlaps, 1)
+ steady_check = 4
+ else:
+ # Error value is 3 when the subdomain of the polygon could not be determined correctly
+ error_check_steadystate2w2 = 3
+
+ return steady_check, overlaps, walls, error_check_steadystate2w2
+
+def plot_particles(M_particle_steady):
+ """ Function to graphically display spheres in 3D """
+ fig = plt.figure()
+ ax = fig.add_subplot(111, projection='3d')
+ # Center coordinates of the steady state particles
+ centers = M_particle_steady[:,[1,2,3]]
+ # Radii
+ radii = M_particle_steady[:,5]/2
+ # Iterate through all spheres
+ for center, radius in zip(centers, radii):
+ # Create a grid of points
+ u = np.linspace(0, 2 * np.pi, 100)
+ v = np.linspace(0, np.pi, 100)
+ x = radius * np.outer(np.cos(u), np.sin(v)) + center[2] # x-axis in the (z, y, x) order
+ y = radius * np.outer(np.sin(u), np.sin(v)) + center[1] # y-axis in the (z, y, x) order
+ z = radius * np.outer(np.ones(np.size(u)), np.cos(v)) + center[0] # z-axis in the (z, y, x) order
+ # Plot the surface
+ ax.plot_surface(x, y, z, color='b', alpha=0.3)
+ # Set the same limits for x, y, and z axes
+ ax.set_xlim(-1000, 1000)
+ ax.set_ylim(-1000, 1000)
+ ax.set_zlim(0, 500)
+ # Set the aspect ratio to be equal
+ ax.set_aspect('equal', adjustable='box')
+ # Set labels
+ ax.set_xlabel('X')
+ ax.set_ylabel('Y')
+ ax.set_zlabel('Z')
+ # Show the plot
+ plt.show()
+
+def voro_tess(M_particle_tess, tess_start, tess_end, diameters_0):
+ """ Function to perform Tesselation """
+ cntr = Container(M_particle_tess, limits=(tess_start, tess_end), radii=(diameters_0/2), periodic=False)
+ # Cell volumes
+ volumes_cells = [round(v.volume(), 3) for v in cntr]
+ # Cell neighbors
+ cell_neighbors = [v.neighbors() for v in cntr]
+ # Nested tuple structure with x,y,z coordinates of the cell face normals in form:
+ # [[(x,y,z), (x,y,z), ...], [(x,y,z), (x,y,z), ...], ...]
+ cell_face_normals = [v.normals() for v in cntr]
+ # Nested tuple structure with x,y,z coordinates of the cell verticies in form:
+ # [[(x,y,z), (x,y,z), ...], [(x,y,z), (x,y,z), ...], ...]
+ cell_vertices = [v.vertices() for v in cntr]
+ # A list of the indices of the vertices of each face.
+ vertices_indices = [v.face_vertices() for v in cntr]
+
+ return volumes_cells, cell_neighbors, cell_face_normals, cell_vertices, vertices_indices
+
+def calc_grow(particle_i, volumes_particles, diameters, cell_face_fixed, distances_faces_max, cell_face_normals, vertices_indices, cell_vertices, display_particle_cell, particle_status, cell_status):
+ # Particle volume
+ volume_particle_i = volumes_particles[particle_i]
+ # Particle radius
+ radius_particle_i = diameters[particle_i] / 2
+ # List of fixed faces of the considered particle
+ cell_face_fixed_i = cell_face_fixed[particle_i]
+ # Verticies of the cell corresponding to the particle
+ cell_vertices_i = cell_vertices[particle_i]
+ # List of maximal distances between particle center and cell faces for the considered particle
+ distances_faces_max_i = distances_faces_max[particle_i]
+ # Normal vectors of cell faces corresponding to the particle
+ cell_face_normals_i = cell_face_normals[particle_i]
+ # Indices of the vertices of each cell face corresponding to the particle.
+ vertices_indices_i = vertices_indices[particle_i]
+ # Cell status
+ cell_status_i = cell_status[particle_i]
+ # Initial guess for solver
+ initial_guess = radius_particle_i
+ # Bounds
+ bnds = [(0, np.max(distances_faces_max_i))]
+ # Solver method
+ solver_method = 'SLSQP'
+ if cell_status_i == False:
+ # Perform the optimization
+ result = minimize(solver_distances, initial_guess, args=(cell_face_fixed_i, distances_faces_max_i, volume_particle_i, cell_face_normals_i, vertices_indices_i, cell_vertices_i, display_particle_cell, radius_particle_i), method=solver_method, bounds=bnds, tol=1e-3)
+ # Obtain the value of the minimized function
+ minimized_value = result.fun
+ # If value to high, mark particle as invalid
+ if minimized_value > 1.0:
+ particle_status[particle_i] = False
+ # Extract the distance
+ distance_face_var = result.x[0]
+ # Iterate through each face of the cell to update the list of fixed cell faces
+ for face in range(len(cell_face_normals_i)):
+ if cell_face_fixed[particle_i][face] == True:
+ continue
+ else:
+ # If the distance is less than or equal to the radius of the sphere, the face intersects with the sphere.
+ # Thus it is set to fixed
+ if (np.abs(distance_face_var) - distances_faces_max[particle_i][face]) > 1e-9:
+ cell_face_fixed[particle_i][face] = True
+ else:
+ # Iterate through each face of the cell to update the list of fixed cell faces
+ for face in range(len(cell_face_normals_i)):
+ if cell_face_fixed[particle_i][face] == True:
+ continue
+ else:
+ cell_face_fixed[particle_i][face] = True
+
+ return cell_face_fixed, particle_status
+
+def solver_distances(params, cell_face_fixed_i, distances_faces_max_i, volume_particle_i, cell_face_normals_i, vertices_indices_i, cell_vertices_i, display_particle_cell, radius_particle_i):
+ """ Solver to calculate the distance between the particle center and the variable cell faces """
+ # Variable distance from particle center to the cell face
+ distance_face_var = params
+ # Minimization value
+ total_error = 0
+ # Calculate new verticies
+ cell_vertices_i_new = calc_verticies(cell_face_normals_i, vertices_indices_i, distances_faces_max_i, distance_face_var, cell_face_fixed_i, cell_vertices_i, display_particle_cell, radius_particle_i)
+ # Check arrayfor NaN values
+ if np.isnan(cell_vertices_i_new).any():
+ total_error = 10000
+ else:
+ # Calculate polyhedron volume
+ volume_polyhedron = calc_volume_polyhedron(cell_vertices_i_new)
+ # Minimization function
+ total_error = 100*abs(volume_particle_i - volume_polyhedron) / volume_particle_i
+
+ return total_error
+
+def calc_verticies(cell_face_normals_i, vertices_indices_i, distances_faces_max_i, distance_face_var, cell_face_fixed_i, cell_vertices_i, display_particle_cell, radius_particle_i):
+ """ Function to calculate new verticies coordinates """
+ # Predefine list with aim distances between particle center and cell faces
+ distances_faces = []
+ # Chech each face of polyhedron
+ for face_index in range(len(cell_face_normals_i)):
+ # If face is from beginning fixed set distance to particle radius
+ if cell_face_fixed_i[face_index] == True:
+ distance_face = distances_faces_max_i[face_index]
+ # If the distance of variable face is bigger than the distance between particle center and the cell face (max. distance)
+ # set the distance equal to the distance between particle center and the cell face
+ elif distance_face_var >= distances_faces_max_i[face_index]:
+ distance_face = distances_faces_max_i[face_index]
+ # Else set distance to variable distance value
+ else:
+ distance_face = distance_face_var
+ # Update list
+ distances_faces.append(distance_face)
+ # Calculate the difference between the aim distances and maximal distances
+ delta_distances = [abs(a - b) for a, b in zip(distances_faces_max_i, distances_faces)]
+ # Number of vertices
+ n_verticies = len(cell_vertices_i)
+ # Predefine list for new vertices
+ cell_vertices_i_new = [(None, None, None) for _ in range(n_verticies)]
+ # Consider each vertex of polyhedron
+ for vertex_index in range(len(cell_vertices_i)):
+ # Find every face that is connected with the considered vertex
+ face_index = [i for i, subarray in enumerate(vertices_indices_i) if vertex_index in subarray]
+ # Define move distances of vetex
+ distances_move = [delta_distances[i] for i in face_index]
+ # Define move vectors of vetex
+ vectors_move = np.array([cell_face_normals_i[i] for i in face_index])
+ # Compute the displacement vector
+ displacement = sum(distance * normal for distance, normal in zip(distances_move, vectors_move))
+ # Calculate the new position coordinates of vertex
+ verticies_coordinates = cell_vertices_i[vertex_index] - displacement
+ # Update list
+ cell_vertices_i_new[vertex_index] = verticies_coordinates
+ # Convert the list of arrays to a 2D NumPy array
+ cell_vertices_i_new = np.vstack(cell_vertices_i_new)
+
+ """ Graphically display particles in 3D """
+ if display_particle_cell == True:
+ plot_particle_cell(cell_vertices_i_new, vertices_indices_i, cell_face_fixed_i, radius_particle_i)
+
+ return cell_vertices_i_new
+
+def plot_particle_cell(cell_vertices_i_new, vertices_indices_i, cell_face_fixed_i, radius_particle_i):
+ """ Function to graphically display particle and cell """
+ # Boolean list indicating which faces to mark in red
+ mark_faces = cell_face_fixed_i
+ # Calculate the convex hull
+ #hull = ConvexHull(cell_vertices_i_new)
+ # Extract faces
+ #faces = hull.simplices
+ # Create a 3D plot
+ fig = plt.figure()
+ ax = fig.add_subplot(111, projection='3d')
+ # Function to get coordinates of vertices for each face
+ def get_face_vertices(face_indices):
+ return [cell_vertices_i_new[idx] for idx in face_indices]
+ # Create Poly3DCollection from faces
+ poly3d = Poly3DCollection([get_face_vertices(face) for face in vertices_indices_i], alpha=0.25, edgecolors='r')
+ # Set face colors based on mark_faces list
+ face_colors = ['red' if mark else 'blue' for mark in mark_faces]
+ poly3d.set_facecolor(face_colors)
+ # Add collection to plot
+ ax.add_collection3d(poly3d)
+ # Plot vertices
+ ax.scatter(cell_vertices_i_new[:, 0], cell_vertices_i_new[:, 1], cell_vertices_i_new[:, 2])
+ # Create a grid of points of the particle
+ u = np.linspace(0, 2 * np.pi, 100)
+ v = np.linspace(0, np.pi, 100)
+ x = radius_particle_i * np.outer(np.cos(u), np.sin(v)) + center_particle_i[0] # x-axis in the (z, y, x) order
+ y = radius_particle_i * np.outer(np.sin(u), np.sin(v)) + center_particle_i[1] # y-axis in the (z, y, x) order
+ z = radius_particle_i * np.outer(np.ones(np.size(u)), np.cos(v)) + center_particle_i[2] # z-axis in the (z, y, x) order
+ # Plot the surface
+ ax.plot_surface(x, y, z, color='b', alpha=0.3)
+ # Label vertices
+ for i, v in enumerate(cell_vertices_i_new):
+ ax.text(v[0], v[1], v[2], str(i), color='black', fontsize=10)
+ # Set the aspect ratio to be equal
+ ax.set_aspect('equal', adjustable='box')
+ # Set labels
+ ax.set_xlabel('X')
+ ax.set_ylabel('Y')
+ ax.set_zlabel('Z')
+ # Show plot
+ plt.show()
+
+def calc_volume_polyhedron(cell_vertices_i_new):
+ """ Function to calculates the polyhedron volume based on the Tetrahedral Decomposition method"""
+ # Polyhedron volume
+ volume_polyhedron = 0
+ # Perform Delaunay triangulation to get tetrahedra
+ delaunay = Delaunay(cell_vertices_i_new)
+ # Consider each tetrahedron
+ for simplex in delaunay.simplices:
+ # Get tetrahedron verticies
+ verticies_tetrahedron = cell_vertices_i_new[simplex]
+ # Calculate the tetrahedron volumes
+ volume_tetrahedron = calc_tetrahedron_volume(verticies_tetrahedron)
+ # Calculate the volume of the polyhedron by summing tetrahedron volumes
+ volume_polyhedron += volume_tetrahedron
+
+ return volume_polyhedron
+
+def calc_tetrahedron_volume(tetra):
+ """ Function to calculates the volume of a tetrahedron based on its vertices. """
+ # Tetrahedron vertices
+ a, b, c, d = tetra
+ # Calculate tetrahedron volume
+ volume_tetrahedron = abs(np.dot(a-d, np.cross(b-d, c-d))) / 6
+
+ return volume_tetrahedron
+
+def solver_grow_factor(params, holdup_end, volume_cube, diameters_0, volumes_cells):
+ """ Solver to calculate the grow factor """
+ # Variable grow factor
+ grow_factor = params
+ # Minimization value
+ total_error = 0
+ # New diameters
+ diameters_new = diameters_0 * grow_factor
+ # New volumes
+ volumes_particles_new = (np.pi / 6) * diameters_new **(3)
+ # Correct volumes
+ for i in range(len(volumes_particles_new)):
+ if volumes_particles_new[i] - volumes_cells[i] > 1e-9:
+ volumes_particles_new[i] = volumes_cells[i]
+ # New holdup
+ holdup_new = np.sum(volumes_particles_new) / volume_cube
+ #
+ total_error = 100*abs(holdup_end - holdup_new) / holdup_end
+
+ return total_error
+
+def calc_contact_probabilities(average_contacts, contact_probabilities, cell_face_fixed, diameters_neighbors, diameters_0, grow_step, particle_status):
+ """
+ This function calculates the average contact number of each class and the contact probability between classes.
+ Particles which are in contact with walls at their initial state are excluded.
+ If solver issues occur at growing simulation of particles those particles are excluded.
+ Wall contacts are excluded.
+ Only contacts between fixed faces are included in the calculation.
+ """
+ """ Count contacts """
+ # Initialize a dictionary to hold the contact counts
+ contact_counts = defaultdict(lambda: defaultdict(int))
+ # Iterate over each particle and its neighbors to count the contacts
+ for particle_i in range(len(diameters_0)):
+ particle = diameters_0[particle_i]
+ for neighbor_j in range(len(diameters_neighbors[particle_i])):
+ neighbor = diameters_neighbors[particle_i][neighbor_j]
+ # Exclude particle with status "False", wall contacts (None values) and faces which are not fixed
+ if particle_status[particle_i] != False and neighbor != None and cell_face_fixed[particle_i][neighbor_j] == True:
+ # Increase number of contacts
+ contact_counts[particle][neighbor] += 1
+ """ Calculate total contacts """
+ # Initialize a dictionary to hold the total contacts for each particle class
+ total_contacts = defaultdict(int)
+ # Sum the total contacts for each particle class
+ for class_particle in contact_counts:
+ total_contacts[class_particle] += sum(contact_counts[class_particle].values())
+ # Sort value_counts by keys (values) in ascending order and convert to a list of tuples
+ total_contacts_sorted = sorted(total_contacts.items())
+ # Extract only the counts (second column)
+ total_contacts_values = [count for value, count in total_contacts_sorted]
+ """ Calculate class sizes """
+ # Count occurrences of each diameter
+ class_size = Counter()
+ # Iterate through data_list and include values based on include_list
+ for class_diameter, status in zip(diameters_0, particle_status):
+ if status == True:
+ class_size[class_diameter] += 1
+ # Sort value_counts by keys (values) in ascending order and extract only the counts
+ class_size_values = [class_size[value] for value in sorted(class_size)]
+ """ Calculate average contact nubmers based on class sizes and total contacts """
+ # Perform element-wise division
+ for i in range(len(class_size_values)):
+ average_contacts[grow_step][i] = total_contacts_values[i] / class_size_values[i]
+ """ Calculate contact probability based on the contact counts and total contacts """
+ # Initialize a dictionary to hold the probabilities
+ contact_probabilities_dict = defaultdict(lambda: defaultdict(float))
+ # Iterate over each particle and its neighbors to calculate the probabilities
+ for class_particle, classes_neighbors in contact_counts.items():
+ for class_neighbor, count_contacts in classes_neighbors.items():
+ # Avoid division by zero
+ if total_contacts[class_particle] > 0:
+ # Calculate the contact probabilities
+ contact_probabilities_dict[class_particle][class_neighbor] = 100 * count_contacts / total_contacts[class_particle]
+ # Convert the contact_probabilities dictionary into a sorted list (sorted in ascending order)
+ contact_counts_sorted = sorted(set(contact_counts))
+ # Iterate over each particle and its neighbors to update the list in ascending order
+ for class_particle_i in range(len(contact_counts_sorted)):
+ class_particle = contact_counts_sorted[class_particle_i]
+ row = []
+ for class_neighbor in contact_counts_sorted:
+ row.append(contact_probabilities_dict[class_particle].get(class_neighbor, 0))
+ contact_probabilities[grow_step][class_particle_i] = row
+
+ return average_contacts, contact_probabilities
+
+def calc_q(diameters, n_particle, d_classes):
+ ' Function to calculate the number and volume distributions'
+ # Predefine array for class sizes after random loose packing simulation
+ class_sizes = np.zeros(len(d_classes), dtype=int)
+ # Categorize cell diameters into classes
+ for particle in diameters:
+ for i in range(len(d_classes)):
+ lb_class = d_classes[i] - delta_class/2
+ ub_class = d_classes[i] + delta_class/2
+ if lb_class <= particle < ub_class:
+ class_sizes[i] += 1
+ break
+ # Final number distribution after growing simulation
+ q0 = class_sizes / float(n_particle)
+ # Final volume distribution after growing simulation
+ q3 = 0.9*(q0 * (np.pi/6) * d_classes**3) / np.sum(q0 * (np.pi/6) * d_classes**3)
+
+ return q0, q3
+
+""" MAIN """
+if __name__ == "__main__":
+ # Time stamp
+ start_time_rlp = time.time()
+ # Set certain variables to integer
+ SEED = int(SEED)
+ N_PARTICLES_SIM = int(N_PARTICLES_SIM)
+ print('Starting simulation with seed', SEED)
+ # Set random seed
+ random.seed(SEED)
+ np.random.seed(SEED)
+ sys.stdout.write("Generating particle size distribution...")
+ # Number distribution of particles
+ q0, n_particles_init, q0_initial, q3_initial, d_classes, d_min, d_max, volume_total, d_05, d_95, delta_class, D_32 = calc_distribution(N_PARTICLES_SIM, N_CLASSES)
+ # Volume of the simulation environment
+ volume_sim = volume_total / HOLDUP_EST
+ # Cuboid's edge lenght (µm)
+ L_0 = volume_sim ** (1/3)
+ # Half of cuboid's edge lenght (µm)
+ A_0 = L_0 / 2
+ # Mix particles in array and determine their coordinates randomly --> initial array for particles
+ M_particles_initial = calc_particles(q0, n_particles_init)
+ # Array for particles in steady state.
+ # Each row corresponds to one particle.
+ # Columns: 0: index of this array, 1: z-coordinate, 2: y-coordinate, 3: x-coordinate, 4: steady state status,
+ # 5: particle diameter, 6: index of initial array, 7: particle class
+ # further columns (len(q0)) corresponds to number of classes of the distribution.
+ # Thus, contact particle classes can be assigned to each particle in array.
+ M_particle_steady = np.zeros((1, (8 + len(q0))))
+ sys.stdout.write("completed\n")
+ # Use tqdm to track progress
+ with tqdm(total=n_particles_init, desc='Random loose packing simulation progress', unit='particle') as pbar:
+ # Drop and roll loop:
+ for particle_i in range(len(M_particles_initial)):
+ if display_detailed_progress == True:
+ print('Particle Nr.:', particle_i)
+ """ Take information of considered particles from initial particle array to define its initial state """
+ x, y, z, index_i, d, class_i = particle_info(M_particles_initial, particle_i)
+ """ Define test range: the distance between the centers of the considered particle and the bigest particlein steady state"""
+ # List of all particles in steady state
+ d_steady = M_particle_steady[:, 5]
+ # Diameter of biggest particle in steady state
+ d_max_steady = np.max(d_steady)
+ # Test range
+ test_range = 1.2*((d + d_max_steady) / 2)
+ # The first particle is simply set because there is nothing there before.
+ if particle_i == 0:
+ # Calculate new z-coordinate
+ z = 0 + d / 2
+ """Fill steady state array"""
+ # Particle index
+ M_particle_steady[0, 0] = particle_i
+ # z-coordinate
+ M_particle_steady[0, 1] = z
+ # y-coordinate
+ M_particle_steady[0, 2] = y
+ # x-coordinate
+ M_particle_steady[0, 3] = x
+ # Contact with ground, 1=yes
+ M_particle_steady[0, 4] = 1
+ # Particle diameter
+ M_particle_steady[0, 5] = d
+ # Initial index
+ M_particle_steady[0, 6] = index_i
+ # Particle class
+ M_particle_steady[0, 7] = class_i
+ continue
+ # Calculate dropping height based on the location of the heighest particle
+ z = find_z_drop(y, x, M_particle_steady, test_range, d_max_steady)
+ # Initial dropping position
+ z_drop = z
+ y_drop = y
+ x_drop = x
+ """
+ While loop for particle sedimentation.
+ In this loop it is checked whether a particle sediments or collides with other particles.
+ If the considered particle does not overlap with other particles, it is moved downwards.
+ If the considered particle overlaps with other particles,
+ depending on the number of overlappings, its position is updated by the corresponding function.
+ The position is updated so that the particles touch each other but do not overlap.
+ After the position update the considered particleis tested for overlapping again because it could
+ overlap with other particles due to the update.
+ Depending on the number of overlaps, the position is updated over and over until no overlapping after
+ the position update occurs. Only when no overlappings are present it is moved downwards.
+ When three or four overlappings are present, the particle is checked for steady state.
+ If it is not in steady state, it is calculated at which particle it would roll along in the next loop step
+ and the corresponding position update function is applied. This is neccessary to prevent the particle to get stucked.
+ This loop is repeated agian and again until the considered particle reches its steady state position at the ground,
+ at three or at four particles.
+ """
+ # Counter variables and lists
+ # Number of overplaps
+ n_overlaps = 0
+ # List of overlaps
+ overlaps = []
+ # Variable for overlapping walls
+ walls = 0
+ # Error Variables
+ error_position_update2 = 0
+ error_position_update3 = 0
+ error_position_update4 = 0
+ error_position_update1w = 0
+ error_position_update2w = 0
+ error_position_update3w = 0
+ error_position_update1w2 = 0
+ error_position_update2w2 = 0
+ error_check_steadystate3 = 0
+ error_check_steadystate4 = 0
+ error_check_steadystate2w = 0
+ error_check_steadystate3w = 0
+ error_check_steadystate1w2 = 0
+ error_check_steadystate2w2 = 0
+ # Argument for the while loop
+ steady = 0
+ # Aborting variable for the while loop
+ counter_sedi = 0
+ # Borders of the aborting variable for the while loop
+ counter_sedi_max1 = 10000
+ counter_sedi_max2 = 2 * counter_sedi_max1
+ counter_sedi_max3 = 3 * counter_sedi_max1
+ # Error counter
+ n_error = 0
+ # If fatal error is True, the while loop is leaved
+ fatal_error = False
+ # While loop fpr particle sedimentation
+ while steady == 0:
+ # Increase aborting variable for the while loop
+ counter_sedi = counter_sedi + 1
+ # Try again if it takes too long or an error occurs
+ if counter_sedi == counter_sedi_max1 or n_error == 1:
+ z = z_drop
+ y = y_drop
+ x = y_drop
+ if n_error == 1:
+ n_error = n_error + 1
+ continue
+ # Assign new coordinates if it takes again too long or an error occurs again
+ if counter_sedi == counter_sedi_max2 or n_error == 3:
+ x = random.uniform(-A_0 + d / 2, A_0 - d / 2)
+ y = random.uniform(-A_0 + d / 2, A_0 - d / 2)
+ z = find_z_drop(y, x, M_particle_steady, test_range, d_max_steady)
+ if n_error == 3:
+ n_error = n_error + 1
+ continue
+ # Abort while loop by abort variable if trying again and assigning new coordinates does not work
+ # or an error occurs three times
+ if counter_sedi > counter_sedi_max3 or n_error == 5:
+ n_error = 0
+ if display_detailed_progress == True:
+ print('Timeout in the sedimentation loop or too many errors. Check sedimentation loop.')
+ fatal_error = True
+ break
+ # Calculate overlapping/collision
+ n_overlaps, overlaps, walls = calc_overlap(d, z, y, x, M_particle_steady, test_range)
+ # If the considered particle does not overlap with other particles, it is moved downwards.
+ if n_overlaps == 0:
+ # If z == d/2, the considered particle reached the ground
+ if z == d/2:
+ # If steady = 1, steady state position is reched on the ground
+ steady = 1
+ else:
+ z = z - STEPSIZE
+ # Overlap with one steady state particle in total
+ elif n_overlaps == 1:
+ # index of the first steady state particle, the considered particle collides with
+ index1 = overlaps[0]
+ # No overlappings with walls
+ if walls == 0:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i ,'collides with one steady state particle.')
+ # If the considered particle overlaps with one steady state particle,
+ # its position is updated by position_update1 function.
+ y, x = position_update1(y, x, d, M_particle_steady, overlaps)
+ # Overlap with one wall
+ elif walls == 1 or walls == -1 or walls == 2 or walls == -2:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i ,'collides with one steady state particle and wall.')
+ # If the considered particle overlaps with one steady state particle and one wall,
+ # its position is updated by position_update1w function.
+ y_s, x_s, error_position_update1w = position_update1w(z, y, x, d, M_particle_steady, overlaps, walls)
+ # If no range error occurs, the new coordinates are assignet
+ if error_position_update1w == 0:
+ y = y_s
+ x = x_s
+ # Error value is 1 if the solver found the same new position of the considered particle twice.
+ # The error is only displayed but the calculation is continued, as the position could be correct
+ elif error_position_update1w == 1:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update1w function. Solver found the same new position twice. Calculation proceeds with the found value.')
+ y = y_s
+ x = x_s
+ # Error value is 2 if the distance between the centers of the two steady state particles and the position
+ # of the first solver solution (the logical one) are not in the same range as the sums of their radii
+ elif error_position_update1w == 2:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update1w function. The logical solution of the new position determined by the solver lies not in the correct range. ')
+ z = z + STEPSIZE / 2
+ continue
+ # Error value is 3 if the both solver solutions lies behind the simulation borders
+ elif error_position_update1w == 3:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update1w function. The both solver solutions lies behind the simulation borders. ')
+ n_error = n_error + 1
+ continue
+ # Overlap with two walls
+ elif walls == 3 or walls == -3 or walls == 4 or walls == -4:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i ,'collides with one steady state particle and two walls.')
+ # If the considered particle overlaps with one steady state particle and two walls,
+ # its position is updated by position_update1w function.
+ z_s, y_s, x_s, error_position_update1w2 = position_update1w2(z, y, x, d, M_particle_steady, overlaps, walls)
+ # Error value is 1 if the solver found the same new position of the considered particle twice.
+ # The error is only displayed but the calculation is continued, as the position could be correct
+ if error_position_update1w2 == 1:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update1w2 function. Solver found the same new position twice. Calculation proceeds with the found value.')
+ # Error value is 2 if the position of the considered particle is located below the centroid of the steady state particles
+ elif error_position_update1w2 == 2:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update1w2 function. Position of the considered particle is located below the centroid of the steady state particles.')
+ z = z + STEPSIZE / 2
+ continue
+ # Error value is 3 if the position of the first solver solution (the logical one) is not located on the intersection plane
+ elif error_position_update1w2 == 3:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update1w2 function. The logical new position determnined by the solver is not located on the intersection plane. ')
+ z = z + STEPSIZE / 2
+ continue
+ # If particleis in contact with one steady state particle and two walls, it could be in steady state
+ # Checking for overlaps after position update
+ n_overlaps, overlaps1w2, walls = calc_overlap(d, z_s, y_s, x_s, M_particle_steady, test_range)
+ # If no overlaps occure, the particle is tested for steady state
+ if n_overlaps == 0:
+ # With three contact partners, the considered particle can be in a steady state,
+ # which is tested by the check_steadystate1w2 function
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i ,'is tested for steady state at one particle and two walls.')
+ steady_check, overlaps, walls, error_check_steadystate1w2 = check_steadystate1w2(z_s, y_s, x_s, M_particle_steady, overlaps, walls)
+ # Error value is 1 if particles center is located below the triangle's centroid of the steady state particles
+ if error_check_steadystate1w2 == 1:
+ if display_detailed_progress == True:
+ print('Error in the check_steadystate3 function. The center of the considered particle is located below the triangles centroid of the steady state particles.')
+ n_error = n_error + 1
+ continue
+ # Error value is 2 if the triangle points of the steady state particles are collinear
+ elif error_check_steadystate1w2 == 2:
+ if display_detailed_progress == True:
+ print('Error in the check_steadystate3 function. The triangle points of the steady state particles are collinear.')
+ n_error = n_error + 1
+ continue
+ # Error value is 3 if the subdomain of the triangle could not be determined correctly
+ elif error_check_steadystate1w2 == 3:
+ if display_detailed_progress == True:
+ print('Error in the check_steadystate3 function. The subdomain of the triangle could not be determined correctly.')
+ n_error = n_error + 1
+ continue
+ # If steady_check == 0, the particle is in steady state
+ if steady_check == 3:
+ # Assigning the coordinates of the considered partcile to its steady state coordinates at three particles
+ z = z_s
+ y = y_s
+ x = x_s
+ # If steady = 2, steady state position is reched on one particle and two walls
+ steady = 6
+ # If steady_check == 1, the particle is not in steady state and would roll along one steady state particle
+ # Thus its position is updated by position_update1 function
+ elif steady_check == 1:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i, 'exits the steady state check and is now rolling along one particle')
+ # Function for position update at one contact.
+ y, x = position_update1(y, x, d, M_particle_steady, overlaps)
+ # If steady_check == 2, the particle is not in steady state and would roll along one steady state particle and wall
+ # Thus its position is updated by position_update2 function
+ elif steady_check == 4:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i, 'exits the steady state check and is now rolling along one particle and wall')
+ # Function for position update at two contacts.
+ y_s, x_s, error_position_update1w = position_update1w(z, y, x, d, M_particle_steady, overlaps, walls)
+ # If no range error occurs, the new coordinates are assignet
+ if error_position_update1w == 0:
+ y = y_s
+ x = x_s
+ # Error value is 1 if the solver found the same new position of the considered particle twice.
+ # The error is only displayed but the calculation is continued, as the position could be correct
+ elif error_position_update1w == 1:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update1w function. Solver found the same new position twice. Calculation proceeds with the found value.')
+ y = y_s
+ x = x_s
+ # Error value is 2 if the distance between the centers of the two steady state particles and the position
+ # of the first solver solution (the logical one) are not in the same range as the sums of their radii
+ elif error_position_update1w == 2:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update1w function. The logical solution of the new position determined by the solver lies not in the correct range. ')
+ z = z + STEPSIZE / 2
+ continue
+ # Error value is 3 if the both solver solutions lies behind the simulation borders
+ elif error_position_update1w == 3:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update1w function. The both solver solutions lies behind the simulation borders. ')
+ n_error = n_error + 1
+ continue
+ continue
+ # If overlaps occur, the considered particle is moved 1/10 of the step upwards an handled again in the next loop.
+ else:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i ,'collides with three steady state particles and two walls.')
+ z = z + STEPSIZE / 2
+ continue
+ # Checking for overlaps after position update
+ n_overlaps, overlaps, walls = calc_overlap(d, z, y, x, M_particle_steady, test_range)
+ # Overlap with two steady state particles in total
+ if n_overlaps == 1:
+ # index of the second steady state particle, the considered particle collides with
+ index2 = overlaps[0]
+ # Update overlaps array
+ overlaps = np.array([index1, index2])
+ # No overlappings with walls
+ if walls == 0:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i ,'collides with two steady state particles.')
+ # If the considered particle overlaps with one steady state particle after position update at one particle,
+ # it is in contact with two particles. Thus, its position is updated by position_update2 function.
+ y_s, x_s, error_position_update2 = position_update2(z, y, x, d, M_particle_steady, overlaps)
+ # If no range error occurs, the new coordinates are assignet
+ if error_position_update2 == 0:
+ y = y_s
+ x = x_s
+ # Error value is 1 if the solver found the same new position of the considered particle twice.
+ # The error is only displayed but the calculation is continued, as the position could be correct
+ elif error_position_update2 == 1:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update2 function. Solver found the same new position twice. Calculation proceeds with the found value.')
+ y = y_s
+ x = x_s
+ # Error value is 2 if the distance between the centers of the two steady state particles and the position
+ # of the first solver solution (the logical one) are not in the same range as the sums of their radii
+ elif error_position_update2 == 2:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update2 function. The logical solution of the new position determined by the solver lies not in the correct range. ')
+ z = z + STEPSIZE / 2
+ continue
+ # Error value is 3 if the both solver solutions lies behind the simulation borders
+ elif error_position_update2 == 3:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update1w function. The both solver solutions lies behind the simulation borders. ')
+ n_error = n_error + 1
+ continue
+ # Overlap with one wall
+ elif walls == 1 or walls == -1 or walls == 2 or walls == -2:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i ,'collides with two steady state particles and one wall.')
+ # If the considered particle overlaps with one steady state particle and one wall after position update at one particle,
+ # it is in contact with two particles and one wall. Thus, its position is updated by position_update2w function.
+ z_s, y_s, x_s, error_position_update2w = position_update2w(z, y, x, d, M_particle_steady, overlaps, walls)
+ # Error value is 1 if the solver found the same new position of the considered particle twice.
+ # The error is only displayed but the calculation is continued, as the position could be correct
+ if error_position_update2w == 1:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update2w function. Solver found the same new position twice. Calculation proceeds with the found value.')
+ # Error value is 2 if the distance between the centers of the two steady state particles and the position
+ # of the first solver solution (the logical one) are not in the same range as the sums of their radii
+ elif error_position_update2w == 2:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update2w function. The logical solution of the new position determined by the solver lies not in the correct range. ')
+ z = z + STEPSIZE / 2
+ continue
+ # If particle is in contact with two steady state particle and one wall, it could be in steady state
+ # Checking for overlaps after position update
+ n_overlaps, overlaps2w, walls = calc_overlap(d, z_s, y_s, x_s, M_particle_steady, test_range)
+ # If the particle overlaps with a second wall, its position is updated in the next loop
+ if walls == 3 or walls == -3 or walls == 4 or walls == -4 or walls == 0:
+ continue
+ # If no overlaps occure, the particle is tested for steady state
+ if n_overlaps == 0:
+ # With three contact partners, the considered particle can be in a steady state,
+ # which is tested by the check_steadystate1w2 function
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i ,'is tested for steady state at two particles and one wall.')
+ steady_check, overlaps, error_check_steadystate2w = check_steadystate2w(z_s, y_s, x_s, particle_i, M_particle_steady, overlaps, walls)
+ # Error value is 1 if particles center is located below the triangle's centroid of the steady state particles
+ if error_check_steadystate2w == 1:
+ if display_detailed_progress == True:
+ print('Error in the check_steadystate3 function. The center of the considered particle is located below the triangles centroid of the steady state particles.')
+ n_error = n_error + 1
+ continue
+ # Error value is 2 if the triangle points of the steady state particles are collinear
+ elif error_check_steadystate2w == 2:
+ if display_detailed_progress == True:
+ print('Error in the check_steadystate3 function. The triangle points of the steady state particles are collinear.')
+ n_error = n_error + 1
+ continue
+ # Error value is 3 if the subdomain of the triangle could not be determined correctly
+ elif error_check_steadystate2w == 3:
+ if display_detailed_progress == True:
+ print('Error in the check_steadystate3 function. The subdomain of the triangle could not be determined correctly.')
+ n_error = n_error + 1
+ continue
+ # If steady_check == 0, the particle is in steady state
+ if steady_check == 3:
+ # Assigning the coordinates of the considered partcile to its steady state coordinates at three particles
+ z = z_s
+ y = y_s
+ x = x_s
+ # If steady = 2, steady state position is reched at two particles and one wall
+ steady = 4
+ # If steady_check == 1, the particle is not in steady state and would roll along one steady state particle
+ # Thus its position is updated by position_update1 function
+ elif steady_check == 1:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i, 'exits the steady state check and is now rolling along one particle')
+ # Function for position update at one contact.
+ y, x = position_update1(y, x, d, M_particle_steady, overlaps)
+ # If steady_check == 2, the particle is not in steady state and would roll along two steady state particles
+ # Thus its position is updated by position_update2 function
+ elif steady_check == 2:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i, 'exits the steady state check and is now rolling along two particles')
+ # Function for position update at two contacts.
+ y_s, x_s, error_position_update2 = position_update2(z, y, x, d, M_particle_steady, overlaps)
+ # If no range error occurs, the new coordinates are assignet
+ if error_position_update2 == 0:
+ y = y_s
+ x = x_s
+ # Error value is 1 if the solver found the same new position of the considered particle twice.
+ # The error is only displayed but the calculation is continued, as the position could be correct
+ elif error_position_update2 == 1:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update2 function. Solver found the same new position twice. Calculation proceeds with the found value.')
+ y = y_s
+ x = x_s
+ # Error value is 2 if the distance between the centers of the two steady state particles and the position
+ # of the first solver solution (the logical one) are not in the same range as the sums of their radii
+ elif error_position_update2 == 2:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update2 function. The logical solution of the new position determined by the solver lies not in the correct range. ')
+ z = z + STEPSIZE / 2
+ continue
+ # Error value is 3 if the both solver solutions lies behind the simulation borders
+ elif error_position_update2 == 3:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update1w function. The both solver solutions lies behind the simulation borders. ')
+ n_error = n_error + 1
+ continue
+ # If steady_check == 4, the particle is not in steady state and would roll along one steady state particle and wall
+ # Thus its position is updated by position_update1w function
+ elif steady_check == 4:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i, 'exits the steady state check and is now rolling along one particle and wall')
+ # Function for position update at two contacts.
+ y_s, x_s, error_position_update1w = position_update1w(z, y, x, d, M_particle_steady, overlaps, walls)
+ # If no range error occurs, the new coordinates are assignet
+ if error_position_update1w == 0:
+ y = y_s
+ x = x_s
+ # Error value is 1 if the solver found the same new position of the considered particle twice.
+ # The error is only displayed but the calculation is continued, as the position could be correct
+ elif error_position_update1w == 1:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update1w function. Solver found the same new position twice. Calculation proceeds with the found value.')
+ y = y_s
+ x = x_s
+ # Error value is 2 if the distance between the centers of the two steady state particles and the position
+ # of the first solver solution (the logical one) are not in the same range as the sums of their radii
+ elif error_position_update1w == 2:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update1w function. The logical solution of the new position determined by the solver lies not in the correct range. ')
+ z = z + STEPSIZE / 2
+ continue
+ # Error value is 3 if the both solver solutions lies behind the simulation borders
+ elif error_position_update1w == 3:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update1w function. The both solver solutions lies behind the simulation borders. ')
+ n_error = n_error + 1
+ continue
+ continue
+ # If overlaps occur, the considered particle is moved 1/10 of the step upwards an handled again in the next loop.
+ else:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i ,'collides with three steady state particles and two walls.')
+ z = z + STEPSIZE / 2
+ continue
+ # Overlap with two walls
+ elif walls == 3 or walls == -3 or walls == 4 or walls == -4:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i ,'collides with one steady state particle and two walls.')
+ # If the considered particle overlaps with one steady state particle and two walls after position update at one particle,
+ # it is in contact with two particles and two walls. Thus, its position is updated by position_update2w function.
+ z_s, y_s, x_s, error_position_update2w2 = position_update2w2(z, y, x, d, M_particle_steady, overlaps, walls)
+ # Error value is 1 if the solver found the same new position of the considered particle twice.
+ # The error is only displayed but the calculation is continued, as the position could be correct
+ if error_position_update2w2 == 1:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update2w2 function. Solver found the same new position twice. Calculation proceeds with the found value.')
+ # Error value is 2 if the distance between the centers of the two steady state particles and the position
+ # of the first solver solution (the logical one) are not in the same range as the sums of their radii
+ elif error_position_update2w2 == 2:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update2w2 function. The logical solution of the new position determined by the solver lies not in the correct range. ')
+ z = z + STEPSIZE / 2
+ continue
+ # If particle is in contact with two steady state particle and two walls, it could be in steady state
+ # Checking for overlaps after position update
+ n_overlaps, overlaps2w2, walls = calc_overlap(d, z_s, y_s, x_s, M_particle_steady, test_range)
+ # If the particle overlaps with onlyone wall, its position is updated in the next loop
+ if walls == 1 or walls == -1 or walls == 2 or walls == -2 or walls == 0:
+ continue
+ # If no overlaps occure, the particle is tested for steady state
+ if n_overlaps == 0:
+ # With three contact partners, the considered particle can be in a steady state,
+ # which is tested by the check_steadystate1w2 function
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i ,'is tested for steady state at two particles and two walls.')
+ steady_check, overlaps, walls, error_check_steadystate2w2 = check_steadystate2w2(z_s, y_s, x_s, M_particle_steady, overlaps, walls)
+ # Error value is 1 if particles center is located below the triangle's centroid of the steady state particles
+ if error_check_steadystate2w2 == 1:
+ if display_detailed_progress == True:
+ print('Error in the check_steadystate3 function. The center of the considered particle is located below the triangles centroid of the steady state particles.')
+ n_error = n_error + 1
+ continue
+ # Error value is 2 if the triangle points of the steady state particles are collinear
+ elif error_check_steadystate2w2 == 2:
+ if display_detailed_progress == True:
+ print('Error in the check_steadystate3 function. The triangle points of the steady state particles are collinear.')
+ n_error = n_error + 1
+ continue
+ # Error value is 3 if the subdomain of the triangle could not be determined correctly
+ elif error_check_steadystate2w2 == 3:
+ if display_detailed_progress == True:
+ print('Error in the check_steadystate3 function. The subdomain of the triangle could not be determined correctly.')
+ n_error = n_error + 1
+ continue
+ # If steady_check == 0, the particle is in steady state
+ if steady_check == 3:
+ # Assigning the coordinates of the considered partcile to its steady state coordinates at three particles
+ z = z_s
+ y = y_s
+ x = x_s
+ # If steady = 2, steady state position is reched at two particles and two walls
+ steady = 7
+ # If steady_check == 1, the particle is not in steady state and would roll along one steady state particle
+ # Thus its position is updated by position_update1 function
+ elif steady_check == 1:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i, 'exits the steady state check and is now rolling along one particle')
+ # Function for position update at one contact.
+ y, x = position_update1(y, x, d, M_particle_steady, overlaps)
+ # If steady_check == 2, the particle is not in steady state and would roll along two steady state particles
+ # Thus its position is updated by position_update2 function
+ elif steady_check == 2:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i, 'exits the steady state check and is now rolling along two particles')
+ # Function for position update at two contacts.
+ y_s, x_s, error_position_update2 = position_update2(z, y, x, d, M_particle_steady, overlaps)
+ # If no range error occurs, the new coordinates are assignet
+ if error_position_update2 == 0:
+ y = y_s
+ x = x_s
+ # Error value is 1 if the solver found the same new position of the considered particle twice.
+ # The error is only displayed but the calculation is continued, as the position could be correct
+ elif error_position_update2 == 1:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update2 function. Solver found the same new position twice. Calculation proceeds with the found value.')
+ y = y_s
+ x = x_s
+ # Error value is 2 if the distance between the centers of the two steady state particles and the position
+ # of the first solver solution (the logical one) are not in the same range as the sums of their radii
+ elif error_position_update2 == 2:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update2 function. The logical solution of the new position determined by the solver lies not in the correct range. ')
+ z = z + STEPSIZE / 2
+ continue
+ # Error value is 3 if the both solver solutions lies behind the simulation borders
+ elif error_position_update2 == 3:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update1w function. The both solver solutions lies behind the simulation borders. ')
+ n_error = n_error + 1
+ continue
+ # If steady_check == 4, the particle is not in steady state and would roll along one steady state particle and wall
+ # Thus its position is updated by position_update1w function
+ elif steady_check == 4:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i, 'exits the steady state check and is now rolling along one particle and wall')
+ # Function for position update at two contacts.
+ y_s, x_s, error_position_update1w = position_update1w(z, y, x, d, M_particle_steady, overlaps, walls)
+ # If no range error occurs, the new coordinates are assignet
+ if error_position_update1w == 0:
+ y = y_s
+ x = x_s
+ # Error value is 1 if the solver found the same new position of the considered particle twice.
+ # The error is only displayed but the calculation is continued, as the position could be correct
+ elif error_position_update1w == 1:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update1w function. Solver found the same new position twice. Calculation proceeds with the found value.')
+ y = y_s
+ x = x_s
+ # Error value is 2 if the distance between the centers of the two steady state particles and the position
+ # of the first solver solution (the logical one) are not in the same range as the sums of their radii
+ elif error_position_update1w == 2:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update1w function. The logical solution of the new position determined by the solver lies not in the correct range. ')
+ z = z + STEPSIZE / 2
+ continue
+ # Error value is 3 if the both solver solutions lies behind the simulation borders
+ elif error_position_update1w == 3:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update1w function. The both solver solutions lies behind the simulation borders. ')
+ n_error = n_error + 1
+ continue
+ continue
+ # If overlaps occur, the considered particle is moved 1/10 of the step upwards an handled again in the next loop.
+ else:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i ,'collides with three steady state particles and two walls.')
+ z = z + STEPSIZE / 2
+ continue
+ # Checking for overlaps after position update
+ n_overlaps, overlaps, walls = calc_overlap(d, z, y, x, M_particle_steady, test_range)
+ # Three overlappings in total
+ if n_overlaps == 1:
+ # index of the third steady state particle, the considered particle collides with
+ index3 = overlaps[0]
+ # Update overlaps array
+ overlaps = np.array([index1, index2, index3])
+ # No overlappings with walls
+ if walls == 0:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i ,'collides with three steady state particles.')
+ # If the considered particle overlaps with one steady state particle after position update at two particles,
+ # it is in contact with three particles. Thus, its position is updated by position_update3 function.
+ z_s, y_s, x_s, error_position_update3 = position_update3(z, y, x, d, M_particle_steady, overlaps)
+ # Error value is 1 if the solver found the same new position of the considered particle twice.
+ # The error is only displayed but the calculation is continued, as the position could be correct
+ if error_position_update3 == 1:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update3 function. Solver found the same new position twice. Calculation proceeds with the found value.')
+ # Error value is 2 if the position of the considered particle is located below the centroid of the steady state particles
+ elif error_position_update3 == 2:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update3 function. Position of the considered particle is located below the centroid of the steady state particles.')
+ z = z + STEPSIZE / 2
+ continue
+ # Error value is 3 if the position of the first solver solution (the logical one) is not located on the intersection plane
+ elif error_position_update3 == 3:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update3 function. The logical new position determnined by the solver is not located on the intersection plane. ')
+ z = z + STEPSIZE / 2
+ continue
+ # If the considered particle overlaps with more than one wall after position update at three particles,
+ # it is moved 1/10 of the step upwards an handled again in the next loop.
+ else:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i ,'collides with three steady state particles and walls.')
+ z = z + STEPSIZE / 2
+ continue
+ # Checking for overlaps after position update
+ n_overlaps, overlaps3, error_overlap = calc_overlap(d, z_s, y_s, x_s, M_particle_steady, test_range)
+ # If value is True, an error has occurred in the calc_overlapping function.
+ # The considered particle is stored in the steady state array with NaN values.
+ if error_overlap == True:
+ if display_detailed_progress == True:
+ print('Error in the calc_overlapping function. Check calc_overlapping function.')
+ fatal_error = True
+ break
+ # If no overlaps occure, the particle is tested for steady state
+ if n_overlaps == 0:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i ,'is checked for steady state at three steady state particles.')
+ # With three contact partners, the considered particle can be in a steady state,
+ # which is tested by the check_steadystate3 function
+ steady_check, overlaps, error_check_steadystate3 = check_steadystate3(z_s, y_s, x_s, M_particle_steady, overlaps)
+ # Error value is 1 if particles center is located below the triangle's centroid of the steady state particles
+ if error_check_steadystate3 == 1:
+ if display_detailed_progress == True:
+ print('Error in the check_steadystate3 function. The center of the considered particle is located below the triangles centroid of the steady state particles.')
+ n_error = n_error + 1
+ continue
+ # Error value is 2 if the triangle points of the steady state particles are collinear
+ elif error_check_steadystate3 == 2:
+ if display_detailed_progress == True:
+ print('Error in the check_steadystate3 function. The triangle points of the steady state particles are collinear.')
+ n_error = n_error + 1
+ continue
+ # Error value is 3 if the subdomain of the triangle could not be determined correctly
+ elif error_check_steadystate3 == 3:
+ if display_detailed_progress == True:
+ print('Error in the check_steadystate3 function. The subdomain of the triangle could not be determined correctly.')
+ n_error = n_error + 1
+ continue
+ # If steady_check == 0, the particle is in steady state
+ if steady_check == 3:
+ # Assigning the coordinates of the considered partcile to its steady state coordinates at three particles
+ z = z_s
+ y = y_s
+ x = x_s
+ # If steady = 2, steady state position is reched on three particles
+ steady = 2
+ # If steady_check == 1, the particle is not in steady state and would roll along one steady state particle
+ # Thus its position is updated by position_update1 function
+ elif steady_check == 1:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i, 'exits the steady state check and is now rolling along one particle')
+ # Function for position update at one contact.
+ y, x = position_update1(y, x, d, M_particle_steady, overlaps)
+ # If steady_check == 2, the particle is not in steady state and would roll along two steady state particles
+ # Thus its position is updated by position_update2 function
+ elif steady_check == 2:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i, 'exits the steady state check and is now rolling along two particles')
+ # Function for position update at two contacts.
+ y_s, x_s, error_position_update2 = position_update2(z, y, x, d, M_particle_steady, overlaps)
+ # If no range error occurs, the new coordinates are assignet
+ if error_position_update2 == 0:
+ y = y_s
+ x = x_s
+ # Error value is 1 if the solver found the same new position of the considered particle twice.
+ # The error is only displayed but the calculation is continued, as the position could be correct
+ elif error_position_update2 == 1:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update2 function. Solver found the same new position twice. Calculation proceeds with the found value.')
+ y = y_s
+ x = x_s
+ # Error value is 2 if the distance between the centers of the two steady state particles and the position
+ # of the first solver solution (the logical one) are not in the same range as the sums of their radii
+ elif error_position_update2 == 2:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update2 function. The logical solution of the new position determined by the solver lies not in the correct range. ')
+ z = z + STEPSIZE / 2
+ continue
+ # Error value is 3 if the both solver solutions lies behind the simulation borders
+ elif error_position_update2 == 3:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update1w function. The both solver solutions lies behind the simulation borders. ')
+ n_error = n_error + 1
+ continue
+ # If the considered particle overlaps with more than one steady state particle after position update at three contacts,
+ # it is in contact with more than four particles.
+ # In this case, the particle is moved 1/10 of the step upwards an handled again in the next loop.
+ elif n_overlaps > 0:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i, 'collides with more than three particles')
+ z = z + STEPSIZE / 2
+ continue
+ # If the considered particle overlaps with more than two steady state particle after position update at two contacts,
+ # it is in contact with more than four particles.
+ # In this case, the particle is moved 1/10 of the step upwards an handled again in the next loop.
+ elif n_overlaps > 1:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i, 'collides with more than three particles')
+ z = z + STEPSIZE / 2
+ continue
+ # Overlap with three steady state particles in total
+ elif n_overlaps == 2:
+ # Indices of the second and third steady state particles, the considered particle collides with
+ index2 = overlaps[0]
+ index3 = overlaps[1]
+ # Update overlaps array
+ overlaps = np.array([index1, index2, index3])
+ # No overlappings with walls
+ if walls == 0:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i ,'collides with three steady state particles.')
+ # If the considered particle overlaps with two steady state particle after position update at one contact,
+ # it is in contact with three particles. Thus, its position is updated by position_update3 function.
+ z_s, y_s, x_s, error_position_update3 = position_update3(z, y, x, d, M_particle_steady, overlaps)
+ # Error value is 1 if the solver found the same new position of the considered particle twice.
+ # The error is only displayed but the calculation is continued, as the position could be correct
+ if error_position_update3 == 1:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update3 function. Solver found the same new position twice. Calculation proceeds with the found value.')
+ # Error value is 2 if the position of the considered particle is located below the centroid of the steady state particles
+ elif error_position_update3 == 2:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update3 function. Position of the considered particle is located below the centroid of the steady state particles.')
+ z = z + STEPSIZE / 2
+ continue
+ # Error value is 3 if the position of the first solver solution (the logical one) is not located on the intersection plane
+ elif error_position_update3 == 3:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update3 function. The logical new position determnined by the solver is not located on the intersection plane. ')
+ z = z + STEPSIZE / 2
+ continue
+ # If the considered particle overlaps with more than one wall after position update at three particles,
+ # it is moved 1/10 of the step upwards an handled again in the next loop.
+ else:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i ,'collides with three steady state particles and walls.')
+ z = z + STEPSIZE / 2
+ continue
+ # Checking for overlaps after position update
+ n_overlaps, overlaps3, error_overlap = calc_overlap(d, z_s, y_s, x_s, M_particle_steady, test_range)
+ # If value is True, an error has occurred in the calc_overlapping function.
+ # The considered particle is stored in the steady state array with NaN values.
+ if error_overlap == True:
+ if display_detailed_progress == True:
+ print('Error in the calc_overlapping function. Check calc_overlapping function.')
+ fatal_error = True
+ break
+ # If no overlaps occure, the particle is tested for steady state
+ if n_overlaps == 0:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i ,'is checked for steady state at three steady state particles.')
+ # With three contact partners, the considered particle can be in a steady state,
+ # which is tested by the check_steadystate3 function
+ steady_check, overlaps, error_check_steadystate3 = check_steadystate3(z_s, y_s, x_s, M_particle_steady, overlaps)
+ # Error value is 1 if particles center is located below the triangle's centroid of the steady state particles
+ if error_check_steadystate3 == 1:
+ if display_detailed_progress == True:
+ print('Error in the check_steadystate3 function. The center of the considered particle is located below the triangles centroid of the steady state particles.')
+ n_error = n_error + 1
+ continue
+ # Error value is 2 if the triangle points of the steady state particles are collinear
+ elif error_check_steadystate3 == 2:
+ if display_detailed_progress == True:
+ print('Error in the check_steadystate3 function. The triangle points of the steady state particles are collinear.')
+ n_error = n_error + 1
+ continue
+ # Error value is 3 if the subdomain of the triangle could not be determined correctly
+ elif error_check_steadystate3 == 3:
+ if display_detailed_progress == True:
+ print('Error in the check_steadystate3 function. The subdomain of the triangle could not be determined correctly.')
+ n_error = n_error + 1
+ continue
+ # If steady_check == 0, the particle is in steady state
+ if steady_check == 3:
+ # Assigning the coordinates of the considered partcile to its steady state coordinates at three particles
+ z = z_s
+ y = y_s
+ x = x_s
+ # If steady = 2, steady state position is reched on three particles
+ steady = 2
+ # If steady_check == 1, the particle is not in steady state and would roll along one steady state particle
+ # Thus its position is updated by position_update1 function
+ elif steady_check == 1:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i, 'exits the steady state check and is now rolling along one particle')
+ # Function for position update at one contact.
+ y, x = position_update1(y, x, d, M_particle_steady, overlaps)
+ # If steady_check == 2, the particle is not in steady state and would roll along two steady state particles
+ # Thus its position is updated by position_update2 function
+ elif steady_check == 2:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i, 'exits the steady state check and is now rolling along two particles')
+ # Function for position update at two contacts.
+ y_s, x_s, error_position_update2 = position_update2(z, y, x, d, M_particle_steady, overlaps)
+ # If no range error occurs, the new coordinates are assignet
+ if error_position_update2 == 0:
+ y = y_s
+ x = x_s
+ # Error value is 1 if the solver found the same new position of the considered particle twice.
+ # The error is only displayed but the calculation is continued, as the position could be correct
+ elif error_position_update2 == 1:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update2 function. Solver found the same new position twice. Calculation proceeds with the found value.')
+ y = y_s
+ x = x_s
+ # Error value is 2 if the distance between the centers of the two steady state particles and the position
+ # of the first solver solution (the logical one) are not in the same range as the sums of their radii
+ elif error_position_update2 == 2:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update2 function. The logical solution of the new position determined by the solver lies not in the correct range. ')
+ z = z + STEPSIZE / 2
+ continue
+ # Error value is 3 if the both solver solutions lies behind the simulation borders
+ elif error_position_update2 == 3:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update1w function. The both solver solutions lies behind the simulation borders. ')
+ n_error = n_error + 1
+ continue
+ # If the considered particle overlaps with more than one steady state particle after position update at three contacts,
+ # it is in contact with more than four particles.
+ # In this case, the particle is moved 1/10 of the step upwards an handled again in the next loop.
+ elif n_overlaps > 0:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i, 'collides with more than three particles')
+ z = z + STEPSIZE / 2
+ continue
+ # If the considered particle overlaps with more than three steady state particle after position update at one contact,
+ # it is in contact with more than four particles.
+ # In this case, the particle is moved 1/10 of the step upwards an handled again in the next loop.
+ elif n_overlaps > 2:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i, 'collides with more than three particles')
+ z = z + STEPSIZE / 2
+ continue
+ # Two overlaps with steady state particles in total
+ elif n_overlaps == 2:
+ # Indices of the second and third steady state particles, the considered particle collides with
+ index1 = overlaps[0]
+ index2 = overlaps[1]
+ # No overlappings with walls
+ if walls == 0:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i ,'collides with two steady state particles.')
+ # If the considered particle overlaps with two steady state particles,
+ # its position is updated by position_update2 function.
+ y_s, x_s, error_position_update2 = position_update2(z, y, x, d, M_particle_steady, overlaps)
+ # If no range error occurs, the new coordinates are assignet
+ if error_position_update2 == 0:
+ y = y_s
+ x = x_s
+ # Error value is 1 if the solver found the same new position of the considered particle twice.
+ # The error is only displayed but the calculation is continued, as the position could be correct
+ elif error_position_update2 == 1:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update2 function. Solver found the same new position twice. Calculation proceeds with the found value.')
+ y = y_s
+ x = x_s
+ # Error value is 2 if the distance between the centers of the two steady state particles and the position
+ # of the first solver solution (the logical one) are not in the same range as the sums of their radii
+ elif error_position_update2 == 2:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update2 function. The logical solution of the new position determined by the solver lies not in the correct range. ')
+ z = z + STEPSIZE / 2
+ continue
+ # Error value is 3 if the both solver solutions lies behind the simulation borders
+ elif error_position_update2 == 3:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update1w function. The both solver solutions lies behind the simulation borders. ')
+ n_error = n_error + 1
+ continue
+ # Overlap with one wall
+ elif walls == 1 or walls == -1 or walls == 2 or walls == -2:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i ,'collides with two steady state particles and one wall.')
+ # If the considered particle overlaps with one steady state particle and one wall after position update at one particle,
+ # it is in contact with two particles and one wall. Thus, its position is updated by position_update2w function.
+ z_s, y_s, x_s, error_position_update2w = position_update2w(z, y, x, d, M_particle_steady, overlaps, walls)
+ # Error value is 1 if the solver found the same new position of the considered particle twice.
+ # The error is only displayed but the calculation is continued, as the position could be correct
+ if error_position_update2w == 1:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update2w function. Solver found the same new position twice. Calculation proceeds with the found value.')
+ # Error value is 2 if the distance between the centers of the two steady state particles and the position
+ # of the first solver solution (the logical one) are not in the same range as the sums of their radii
+ elif error_position_update2w == 2:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update2w function. The logical solution of the new position determined by the solver lies not in the correct range. ')
+ z = z + STEPSIZE / 2
+ continue
+ # If particle is in contact with two steady state particle and one wall, it could be in steady state
+ # Checking for overlaps after position update
+ n_overlaps, overlaps2w, walls = calc_overlap(d, z_s, y_s, x_s, M_particle_steady, test_range)
+ # If the particle overlaps with a second wall, its position is updated in the next loop
+ if walls == 3 or walls == -3 or walls == 4 or walls == -4 or walls == 0:
+ continue
+ # If no overlaps occure, the particle is tested for steady state
+ if n_overlaps == 0:
+ # With three contact partners, the considered particle can be in a steady state,
+ # which is tested by the check_steadystate1w2 function
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i ,'is tested for steady state at two particles and one wall.')
+ steady_check, overlaps, error_check_steadystate2w = check_steadystate2w(z_s, y_s, x_s, particle_i, M_particle_steady, overlaps, walls)
+ # Error value is 1 if particles center is located below the triangle's centroid of the steady state particles
+ if error_check_steadystate2w == 1:
+ if display_detailed_progress == True:
+ print('Error in the check_steadystate3 function. The center of the considered particle is located below the triangles centroid of the steady state particles.')
+ n_error = n_error + 1
+ continue
+ # Error value is 2 if the triangle points of the steady state particles are collinear
+ elif error_check_steadystate2w == 2:
+ if display_detailed_progress == True:
+ print('Error in the check_steadystate3 function. The triangle points of the steady state particles are collinear.')
+ n_error = n_error + 1
+ continue
+ # Error value is 3 if the subdomain of the triangle could not be determined correctly
+ elif error_check_steadystate2w == 3:
+ if display_detailed_progress == True:
+ print('Error in the check_steadystate3 function. The subdomain of the triangle could not be determined correctly.')
+ n_error = n_error + 1
+ continue
+ # If steady_check == 0, the particle is in steady state
+ if steady_check == 3:
+ # Assigning the coordinates of the considered partcile to its steady state coordinates at three particles
+ z = z_s
+ y = y_s
+ x = x_s
+ # If steady = 2, steady state position is reched at two particles and one wall
+ steady = 4
+ # If steady_check == 1, the particle is not in steady state and would roll along one steady state particle
+ # Thus its position is updated by position_update1 function
+ elif steady_check == 1:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i, 'exits the steady state check and is now rolling along one particle')
+ # Function for position update at one contact.
+ y, x = position_update1(y, x, d, M_particle_steady, overlaps)
+ # If steady_check == 2, the particle is not in steady state and would roll along two steady state particles
+ # Thus its position is updated by position_update2 function
+ elif steady_check == 2:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i, 'exits the steady state check and is now rolling along two particles')
+ # Function for position update at two contacts.
+ y_s, x_s, error_position_update2 = position_update2(z, y, x, d, M_particle_steady, overlaps)
+ # If no range error occurs, the new coordinates are assignet
+ if error_position_update2 == 0:
+ y = y_s
+ x = x_s
+ # Error value is 1 if the solver found the same new position of the considered particle twice.
+ # The error is only displayed but the calculation is continued, as the position could be correct
+ elif error_position_update2 == 1:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update2 function. Solver found the same new position twice. Calculation proceeds with the found value.')
+ y = y_s
+ x = x_s
+ # Error value is 2 if the distance between the centers of the two steady state particles and the position
+ # of the first solver solution (the logical one) are not in the same range as the sums of their radii
+ elif error_position_update2 == 2:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update2 function. The logical solution of the new position determined by the solver lies not in the correct range. ')
+ z = z + STEPSIZE / 2
+ continue
+ # Error value is 3 if the both solver solutions lies behind the simulation borders
+ elif error_position_update2 == 3:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update1w function. The both solver solutions lies behind the simulation borders. ')
+ n_error = n_error + 1
+ continue
+ # If steady_check == 4, the particle is not in steady state and would roll along one steady state particle and wall
+ # Thus its position is updated by position_update1w function
+ elif steady_check == 4:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i, 'exits the steady state check and is now rolling along one particle and wall')
+ # Function for position update at two contacts.
+ y_s, x_s, error_position_update1w = position_update1w(z, y, x, d, M_particle_steady, overlaps, walls)
+ # If no range error occurs, the new coordinates are assignet
+ if error_position_update1w == 0:
+ y = y_s
+ x = x_s
+ # Error value is 1 if the solver found the same new position of the considered particle twice.
+ # The error is only displayed but the calculation is continued, as the position could be correct
+ elif error_position_update1w == 1:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update1w function. Solver found the same new position twice. Calculation proceeds with the found value.')
+ y = y_s
+ x = x_s
+ # Error value is 2 if the distance between the centers of the two steady state particles and the position
+ # of the first solver solution (the logical one) are not in the same range as the sums of their radii
+ elif error_position_update1w == 2:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update1w function. The logical solution of the new position determined by the solver lies not in the correct range. ')
+ z = z + STEPSIZE / 2
+ continue
+ # Error value is 3 if the both solver solutions lies behind the simulation borders
+ elif error_position_update1w == 3:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update1w function. The both solver solutions lies behind the simulation borders. ')
+ n_error = n_error + 1
+ continue
+ continue
+ # If overlaps occur, the considered particle is moved 1/10 of the step upwards an handled again in the next loop.
+ else:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i ,'collides with three steady state particles and two walls.')
+ z = z + STEPSIZE / 2
+ continue
+ # Overlap with two walls
+ elif walls == 3 or walls == -3 or walls == 4 or walls == -4:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i ,'collides with one steady state particle and two walls.')
+ # If the considered particle overlaps with one steady state particle and two walls after position update at one particle,
+ # it is in contact with two particles and two walls. Thus, its position is updated by position_update2w function.
+ z_s, y_s, x_s, error_position_update2w2 = position_update2w2(z, y, x, d, M_particle_steady, overlaps, walls)
+ # Error value is 1 if the solver found the same new position of the considered particle twice.
+ # The error is only displayed but the calculation is continued, as the position could be correct
+ if error_position_update2w2 == 1:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update2w2 function. Solver found the same new position twice. Calculation proceeds with the found value.')
+ # Error value is 2 if the distance between the centers of the two steady state particles and the position
+ # of the first solver solution (the logical one) are not in the same range as the sums of their radii
+ elif error_position_update2w2 == 2:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update2w2 function. The logical solution of the new position determined by the solver lies not in the correct range. ')
+ z = z + STEPSIZE / 2
+ continue
+ # If particle is in contact with two steady state particle and two walls, it could be in steady state
+ # Checking for overlaps after position update
+ n_overlaps, overlaps2w2, walls = calc_overlap(d, z_s, y_s, x_s, M_particle_steady, test_range)
+ # If the particle overlaps with onlyone wall, its position is updated in the next loop
+ if walls == 1 or walls == -1 or walls == 2 or walls == -2 or walls == 0:
+ continue
+ # If no overlaps occure, the particle is tested for steady state
+ if n_overlaps == 0:
+ # With three contact partners, the considered particle can be in a steady state,
+ # which is tested by the check_steadystate1w2 function
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i ,'is tested for steady state at two particles and two walls.')
+ steady_check, overlaps, walls, error_check_steadystate2w2 = check_steadystate2w2(z_s, y_s, x_s, M_particle_steady, overlaps, walls)
+ # Error value is 1 if particles center is located below the triangle's centroid of the steady state particles
+ if error_check_steadystate2w2 == 1:
+ if display_detailed_progress == True:
+ print('Error in the check_steadystate3 function. The center of the considered particle is located below the triangles centroid of the steady state particles.')
+ n_error = n_error + 1
+ continue
+ # Error value is 2 if the triangle points of the steady state particles are collinear
+ elif error_check_steadystate2w2 == 2:
+ if display_detailed_progress == True:
+ print('Error in the check_steadystate3 function. The triangle points of the steady state particles are collinear.')
+ n_error = n_error + 1
+ continue
+ # Error value is 3 if the subdomain of the triangle could not be determined correctly
+ elif error_check_steadystate2w2 == 3:
+ if display_detailed_progress == True:
+ print('Error in the check_steadystate3 function. The subdomain of the triangle could not be determined correctly.')
+ n_error = n_error + 1
+ continue
+ # If steady_check == 0, the particle is in steady state
+ if steady_check == 3:
+ # Assigning the coordinates of the considered partcile to its steady state coordinates at three particles
+ z = z_s
+ y = y_s
+ x = x_s
+ # If steady = 2, steady state position is reched at two particles and two walls
+ steady = 7
+ # If steady_check == 1, the particle is not in steady state and would roll along one steady state particle
+ # Thus its position is updated by position_update1 function
+ elif steady_check == 1:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i, 'exits the steady state check and is now rolling along one particle')
+ # Function for position update at one contact.
+ y, x = position_update1(y, x, d, M_particle_steady, overlaps)
+ # If steady_check == 2, the particle is not in steady state and would roll along two steady state particles
+ # Thus its position is updated by position_update2 function
+ elif steady_check == 2:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i, 'exits the steady state check and is now rolling along two particles')
+ # Function for position update at two contacts.
+ y_s, x_s, error_position_update2 = position_update2(z, y, x, d, M_particle_steady, overlaps)
+ # If no range error occurs, the new coordinates are assignet
+ if error_position_update2 == 0:
+ y = y_s
+ x = x_s
+ # Error value is 1 if the solver found the same new position of the considered particle twice.
+ # The error is only displayed but the calculation is continued, as the position could be correct
+ elif error_position_update2 == 1:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update2 function. Solver found the same new position twice. Calculation proceeds with the found value.')
+ y = y_s
+ x = x_s
+ # Error value is 2 if the distance between the centers of the two steady state particles and the position
+ # of the first solver solution (the logical one) are not in the same range as the sums of their radii
+ elif error_position_update2 == 2:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update2 function. The logical solution of the new position determined by the solver lies not in the correct range. ')
+ z = z + STEPSIZE / 2
+ continue
+ # Error value is 3 if the both solver solutions lies behind the simulation borders
+ elif error_position_update2 == 3:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update1w function. The both solver solutions lies behind the simulation borders. ')
+ n_error = n_error + 1
+ continue
+ # If steady_check == 4, the particle is not in steady state and would roll along one steady state particle and wall
+ # Thus its position is updated by position_update1w function
+ elif steady_check == 4:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i, 'exits the steady state check and is now rolling along one particle and wall')
+ # Function for position update at two contacts.
+ y_s, x_s, error_position_update1w = position_update1w(z, y, x, d, M_particle_steady, overlaps, walls)
+ # If no range error occurs, the new coordinates are assignet
+ if error_position_update1w == 0:
+ y = y_s
+ x = x_s
+ # Error value is 1 if the solver found the same new position of the considered particle twice.
+ # The error is only displayed but the calculation is continued, as the position could be correct
+ elif error_position_update1w == 1:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update1w function. Solver found the same new position twice. Calculation proceeds with the found value.')
+ y = y_s
+ x = x_s
+ # Error value is 2 if the distance between the centers of the two steady state particles and the position
+ # of the first solver solution (the logical one) are not in the same range as the sums of their radii
+ elif error_position_update1w == 2:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update1w function. The logical solution of the new position determined by the solver lies not in the correct range. ')
+ z = z + STEPSIZE / 2
+ continue
+ # Error value is 3 if the both solver solutions lies behind the simulation borders
+ elif error_position_update1w == 3:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update1w function. The both solver solutions lies behind the simulation borders. ')
+ n_error = n_error + 1
+ continue
+ continue
+ # If overlaps occur, the considered particle is moved 1/10 of the step upwards an handled again in the next loop.
+ else:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i ,'collides with three steady state particles and two walls.')
+ z = z + STEPSIZE / 2
+ continue
+ # Checking for overlaps after position update
+ n_overlaps, overlaps, walls = calc_overlap(d, z, y, x, M_particle_steady, test_range)
+ # Three overlaps with steady state particles in total
+ if n_overlaps == 1:
+ # Index of the second and third steady state particles, the considered particle collides with
+ index3 = overlaps[0]
+ # Update overlaps array
+ overlaps = np.array([index1, index2, index3])
+ # No overlappings with walls
+ if walls == 0:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i ,'collides with three steady state particles.')
+ # If the considered particle overlaps with one steady state particle after position update at two contacts,
+ # it is in contact with three particles. Thus, its position is updated by position_update3 function.
+ z_s, y_s, x_s, error_position_update3 = position_update3(z, y, x, d, M_particle_steady, overlaps)
+ # Error value is 1 if the solver found the same new position of the considered particle twice.
+ # The error is only displayed but the calculation is continued, as the position could be correct
+ if error_position_update3 == 1:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update3 function. Solver found the same new position twice. Calculation proceeds with the found value.')
+ # Error value is 2 if the position of the considered particle is located below the centroid of the steady state particles
+ elif error_position_update3 == 2:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update3 function. Position of the considered particle is located below the centroid of the steady state particles.')
+ z = z + STEPSIZE / 2
+ continue
+ # Error value is 3 if the position of the first solver solution (the logical one) is not located on the intersection plane
+ elif error_position_update3 == 3:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update3 function. The logical new position determnined by the solver is not located on the intersection plane. ')
+ z = z + STEPSIZE / 2
+ continue
+ # If the considered particle overlaps with more than one wall after position update at three particles,
+ # it is moved 1/10 of the step upwards an handled again in the next loop.
+ else:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i ,'collides with three steady state particles and walls.')
+ z = z + STEPSIZE / 2
+ continue
+ # Checking for overlaps after position update
+ n_overlaps, overlaps3, error_overlap = calc_overlap(d, z_s, y_s, x_s, M_particle_steady, test_range)
+ # If value is True, an error has occurred in the calc_overlapping function.
+ # The considered particle is stored in the steady state array with NaN values.
+ if error_overlap == True:
+ if display_detailed_progress == True:
+ print('Error in the calc_overlapping function. Check calc_overlapping function.')
+ fatal_error = True
+ break
+ # If no overlaps occure, the particle is tested for steady state
+ if n_overlaps == 0:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i ,'is checked for steady state at three steady state particles.')
+ # With three contact partners, the considered particle can be in a steady state,
+ # which is tested by the check_steadystate3 function
+ steady_check, overlaps, error_check_steadystate3 = check_steadystate3(z_s, y_s, x_s, M_particle_steady, overlaps)
+ # Error value is 1 if particles center is located below the triangle's centroid of the steady state particles
+ if error_check_steadystate3 == 1:
+ if display_detailed_progress == True:
+ print('Error in the check_steadystate3 function. The center of the considered particle is located below the triangles centroid of the steady state particles.')
+ n_error = n_error + 1
+ continue
+ # Error value is 2 if the triangle points of the steady state particles are collinear
+ elif error_check_steadystate3 == 2:
+ if display_detailed_progress == True:
+ print('Error in the check_steadystate3 function. The triangle points of the steady state particles are collinear.')
+ n_error = n_error + 1
+ continue
+ # Error value is 3 if the subdomain of the triangle could not be determined correctly
+ elif error_check_steadystate3 == 3:
+ if display_detailed_progress == True:
+ print('Error in the check_steadystate3 function. The subdomain of the triangle could not be determined correctly.')
+ n_error = n_error + 1
+ continue
+ # If steady_check == 0, the particle is in steady state
+ if steady_check == 3:
+ # Assigning the coordinates of the considered partcile to its steady state coordinates at three particles
+ z = z_s
+ y = y_s
+ x = x_s
+ # If steady = 2, steady state position is reched on three particles
+ steady = 2
+ # If steady_check == 1, the particle is not in steady state and would roll along one steady state particle
+ # Thus its position is updated by position_update1 function
+ elif steady_check == 1:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i, 'exits the steady state check and is now rolling along one particle')
+ # Function for position update at one contact.
+ y, x = position_update1(y, x, d, M_particle_steady, overlaps)
+ # If steady_check == 2, the particle is not in steady state and would roll along two steady state particles
+ # Thus its position is updated by position_update2 function
+ elif steady_check == 2:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i, 'exits the steady state check and is now rolling along two particles')
+ # Function for position update at two contacts.
+ y_s, x_s, error_position_update2 = position_update2(z, y, x, d, M_particle_steady, overlaps)
+ # If no range error occurs, the new coordinates are assignet
+ if error_position_update2 == 0:
+ y = y_s
+ x = x_s
+ # Error value is 1 if the solver found the same new position of the considered particle twice.
+ # The error is only displayed but the calculation is continued, as the position could be correct
+ elif error_position_update2 == 1:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update2 function. Solver found the same new position twice. Calculation proceeds with the found value.')
+ y = y_s
+ x = x_s
+ # Error value is 2 if the distance between the centers of the two steady state particles and the position
+ # of the first solver solution (the logical one) are not in the same range as the sums of their radii
+ elif error_position_update2 == 2:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update2 function. The logical solution of the new position determined by the solver lies not in the correct range. ')
+ z = z + STEPSIZE / 2
+ continue
+ # Error value is 3 if the both solver solutions lies behind the simulation borders
+ elif error_position_update2 == 3:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update1w function. The both solver solutions lies behind the simulation borders. ')
+ n_error = n_error + 1
+ continue
+ # If the considered particle overlaps with more than one steady state particle after position update at three contacts,
+ # it is in contact with more than four particles.
+ # In this case, the particle is moved 1/10 of the step upwards an handled again in the next loop.
+ elif n_overlaps > 0:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i, 'collides with more than three particles')
+ z = z + STEPSIZE / 2
+ continue
+ # If the considered particle overlaps with more than two steady state particle after position update at two contacts,
+ # it is in contact with more than four particles.
+ # In this case, the particle is moved 1/10 of the step upwards an handled again in the next loop.
+ elif n_overlaps > 1:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i, 'collides with more than four particles')
+ z = z + STEPSIZE / 2
+ continue
+ # Three overlaps with steady state particles in total
+ elif n_overlaps == 3:
+ # Indices of the first, second and third steady state particles, the considered particle collides with
+ index1 = overlaps[0]
+ index2 = overlaps[1]
+ index3 = overlaps[2]
+ # No overlappings with walls
+ if walls == 0:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i ,'collides with three steady state particles.')
+ # If the considered particle overlaps with three steady state particles,
+ # its position is updated by position_update3 function.
+ z_s, y_s, x_s, error_position_update3 = position_update3(z, y, x, d, M_particle_steady, overlaps)
+ # Error value is 1 if the solver found the same new position of the considered particle twice.
+ # The error is only displayed but the calculation is continued, as the position could be correct
+ if error_position_update3 == 1:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update3 function. Solver found the same new position twice. Calculation proceeds with the found value.')
+ # Error value is 2 if the position of the considered particle is located below the centroid of the steady state particles
+ elif error_position_update3 == 2:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update3 function. Position of the considered particle is located below the centroid of the steady state particles.')
+ z = z + STEPSIZE / 2
+ continue
+ # Error value is 3 if the position of the first solver solution (the logical one) is not located on the intersection plane
+ elif error_position_update3 == 3:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update3 function. The logical new position determnined by the solver is not located on the intersection plane. ')
+ z = z + STEPSIZE / 2
+ continue
+ # If the considered particle overlaps with more than one wall after position update at three particles,
+ # it is moved 1/10 of the step upwards an handled again in the next loop.
+ else:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i ,'collides with three steady state particles and walls.')
+ z = z + STEPSIZE / 2
+ continue
+ # Checking for overlaps after position update
+ n_overlaps, overlaps3, error_overlap = calc_overlap(d, z_s, y_s, x_s, M_particle_steady, test_range)
+ # If no overlaps occure, the particle is tested for steady state
+ if n_overlaps == 0:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i ,'is checked for steady state at three steady state particles.')
+ # With three contact partners, the considered particle can be in a steady state,
+ # which is tested by the check_steadystate3 function
+ steady_check, overlaps, error_check_steadystate3 = check_steadystate3(z_s, y_s, x_s, M_particle_steady, overlaps)
+ # Error value is 1 if particles center is located below the triangle's centroid of the steady state particles
+ if error_check_steadystate3 == 1:
+ if display_detailed_progress == True:
+ print('Error in the check_steadystate3 function. The center of the considered particle is located below the triangles centroid of the steady state particles.')
+ n_error = n_error + 1
+ continue
+ # Error value is 2 if the triangle points of the steady state particles are collinear
+ elif error_check_steadystate3 == 2:
+ if display_detailed_progress == True:
+ print('Error in the check_steadystate3 function. The triangle points of the steady state particles are collinear.')
+ n_error = n_error + 1
+ continue
+ # Error value is 3 if the subdomain of the triangle could not be determined correctly
+ elif error_check_steadystate3 == 3:
+ if display_detailed_progress == True:
+ print('Error in the check_steadystate3 function. The subdomain of the triangle could not be determined correctly.')
+ n_error = n_error + 1
+ continue
+ # If steady_check == 0, the particle is in steady state
+ if steady_check == 3:
+ # Assigning the coordinates of the considered partcile to its steady state coordinates at three particles
+ z = z_s
+ y = y_s
+ x = x_s
+ # If steady = 2, steady state position is reched on three particles
+ steady = 2
+ # If steady_check == 1, the particle is not in steady state and would roll along one steady state particle
+ # Thus its position is updated by position_update1 function
+ elif steady_check == 1:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i, 'exits the steady state check and is now rolling along one particle')
+ # Function for position update at one contact.
+ y, x = position_update1(y, x, d, M_particle_steady, overlaps)
+ # If steady_check == 2, the particle is not in steady state and would roll along two steady state particles
+ # Thus its position is updated by position_update2 function
+ elif steady_check == 2:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i, 'exits the steady state check and is now rolling along two particles')
+ # Function for position update at two contacts.
+ y_s, x_s, error_position_update2 = position_update2(z, y, x, d, M_particle_steady, overlaps)
+ # If no range error occurs, the new coordinates are assignet
+ if error_position_update2 == 0:
+ y = y_s
+ x = x_s
+ # Error value is 1 if the solver found the same new position of the considered particle twice.
+ # The error is only displayed but the calculation is continued, as the position could be correct
+ elif error_position_update2 == 1:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update2 function. Solver found the same new position twice. Calculation proceeds with the found value.')
+ y = y_s
+ x = x_s
+ # Error value is 2 if the distance between the centers of the two steady state particles and the position
+ # of the first solver solution (the logical one) are not in the same range as the sums of their radii
+ elif error_position_update2 == 2:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update2 function. The logical solution of the new position determined by the solver lies not in the correct range. ')
+ z = z + STEPSIZE / 2
+ continue
+ # Error value is 3 if the both solver solutions lies behind the simulation borders
+ elif error_position_update2 == 3:
+ if display_detailed_progress == True:
+ print('Error in the error_position_update1w function. The both solver solutions lies behind the simulation borders. ')
+ n_error = n_error + 1
+ continue
+ # If the considered particle overlaps with more than one steady state particle after position update at three contacts,
+ # it is in contact with more than four particles.
+ # In this case, the particle is moved 1/10 of the step upwards an handled again in the next loop.
+ elif n_overlaps > 0:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i, 'collides with more than four particles')
+ z = z + STEPSIZE / 2
+ continue
+ # If the considered particle overlaps with more than four steady state particle,
+ # it is in contact with more than four particles.
+ # In this case, the particle is moved 1/10 of the step upwards an handled again in the next loop.
+ elif n_overlaps > 3:
+ if display_detailed_progress == True:
+ print('Particle Nr.', particle_i, 'collides with more than four particles')
+ z = z + STEPSIZE / 2
+ continue
+ # If particle reaches the ground, it is stored in steady state array.
+ if z <= d / 2:
+ z = d / 2
+ """ Display final particle status and store the steady state particle in the steady state array. """
+ # Display if the considered particle is in steady state
+ if display_detailed_progress == True:
+ if steady == 1:
+ print('Particle Nr.', particle_i, 'reaches the ground.')
+ elif steady == 2:
+ print('Particle Nr.', particle_i, 'has reached steady state on three particles.')
+ elif steady == 3:
+ print('Particle Nr.', particle_i, 'has reached steady state on three particles.')
+ elif steady == 4:
+ print('Particle Nr.', particle_i, 'has reached steady state on two particles and one wall.')
+ elif steady == 5:
+ print('Particle Nr.', particle_i, 'has reached steady state on three particles and one wall.')
+ elif steady == 6:
+ print('Particle Nr.', particle_i, 'has reached steady state on one particles and two walls.')
+ elif steady == 7:
+ print('Particle Nr.', particle_i, 'has reached steady state on two particles and two walls.')
+ # If error=True, timeout in the calc_sedimentation function occured.
+ # The considered particle is stored in the steady state array with NaN values.
+ if fatal_error == True:
+ M_particle_steady = np.vstack([M_particle_steady, np.zeros((1, 8 + len(q0)))])
+ M_particle_steady[len(M_particle_steady) - 1, 0] = particle_i
+ M_particle_steady[len(M_particle_steady) - 1, 1] = np.nan
+ M_particle_steady[len(M_particle_steady) - 1, 2] = np.nan
+ M_particle_steady[len(M_particle_steady) - 1, 3] = np.nan
+ M_particle_steady[len(M_particle_steady) - 1, 4] = np.nan
+ M_particle_steady[len(M_particle_steady) - 1, 5] = d
+ M_particle_steady[len(M_particle_steady) - 1, 6] = index_i
+ M_particle_steady[len(M_particle_steady) - 1, 7] = class_i
+ # The considered particle is stored in the steady state array.
+ else:
+ M_particle_steady = np.vstack([M_particle_steady, np.zeros((1, 8 + len(q0)))])
+ M_particle_steady[len(M_particle_steady) - 1, 0] = particle_i
+ M_particle_steady[len(M_particle_steady) - 1, 1] = z
+ M_particle_steady[len(M_particle_steady) - 1, 2] = y
+ M_particle_steady[len(M_particle_steady) - 1, 3] = x
+ M_particle_steady[len(M_particle_steady) - 1, 4] = steady
+ M_particle_steady[len(M_particle_steady) - 1, 5] = d
+ M_particle_steady[len(M_particle_steady) - 1, 6] = index_i
+ M_particle_steady[len(M_particle_steady) - 1, 7] = class_i
+ # Store contact partners in the steady state array.
+ # Determine number of collisions
+ n_overlaps_final = len(overlaps)
+ if n_overlaps_final == 1:
+ # Particle classes of the particles in contact
+ class1 = int(M_particle_steady[overlaps[0], 7])
+ # The classes of the particles in contact are assigned to the considered particle
+ M_particle_steady[len(M_particle_steady) - 1, 8 + class1] = M_particle_steady[len(M_particle_steady) - 1, 8 + class1] + 1
+ # The class of the considered particle is assigned to the particles in contact
+ M_particle_steady[overlaps[0], 8 + class_i] = M_particle_steady[overlaps[0], 8 + class_i] + 1
+ elif n_overlaps_final == 2:
+ # Particle classes of the particles in contact
+ class1 = int(M_particle_steady[overlaps[0], 7])
+ class2 = int(M_particle_steady[overlaps[1], 7])
+ # The classes of the particles in contact are assigned to the considered particle
+ M_particle_steady[len(M_particle_steady) - 1, 8 + class1] = M_particle_steady[len(M_particle_steady) - 1, 8 + class1] + 1
+ M_particle_steady[len(M_particle_steady) - 1, 8 + class2] = M_particle_steady[len(M_particle_steady) - 1, 8 + class2] + 1
+ # The class of the considered particle is assigned to the particles in contact
+ M_particle_steady[overlaps[0], 8 + class_i] = M_particle_steady[overlaps[0], 8 + class_i] + 1
+ M_particle_steady[overlaps[1], 8 + class_i] = M_particle_steady[overlaps[1], 8 + class_i] + 1
+ elif n_overlaps_final == 3:
+ # Particle classes of the particles in contact
+ class1 = int(M_particle_steady[overlaps[0], 7])
+ class2 = int(M_particle_steady[overlaps[1], 7])
+ class3 = int(M_particle_steady[overlaps[2], 7])
+ # The classes of the particles in contact are assigned to the considered particle
+ M_particle_steady[len(M_particle_steady) - 1, 8 + class1] = M_particle_steady[len(M_particle_steady) - 1, 8 + class1] + 1
+ M_particle_steady[len(M_particle_steady) - 1, 8 + class2] = M_particle_steady[len(M_particle_steady) - 1, 8 + class2] + 1
+ M_particle_steady[len(M_particle_steady) - 1, 8 + class3] = M_particle_steady[len(M_particle_steady) - 1, 8 + class3] + 1
+ # The class of the considered particle is assigned to the particles in contact
+ M_particle_steady[overlaps[0], 8 + class_i] = M_particle_steady[overlaps[0], 8 + class_i] + 1
+ M_particle_steady[overlaps[1], 8 + class_i] = M_particle_steady[overlaps[1], 8 + class_i] + 1
+ M_particle_steady[overlaps[2], 8 + class_i] = M_particle_steady[overlaps[2], 8 + class_i] + 1
+ elif n_overlaps_final == 4:
+ # Particle classes of the particles in contact
+ class1 = int(M_particle_steady[overlaps[0], 7])
+ class2 = int(M_particle_steady[overlaps[1], 7])
+ class3 = int(M_particle_steady[overlaps[2], 7])
+ class4 = int(M_particle_steady[overlaps[3], 7])
+ # The classes of the particles in contact are assigned to the considered particle
+ M_particle_steady[len(M_particle_steady) - 1, 8 + class1] = M_particle_steady[len(M_particle_steady) - 1, 8 + class1] + 1
+ M_particle_steady[len(M_particle_steady) - 1, 8 + class2] = M_particle_steady[len(M_particle_steady) - 1, 8 + class2] + 1
+ M_particle_steady[len(M_particle_steady) - 1, 8 + class3] = M_particle_steady[len(M_particle_steady) - 1, 8 + class3] + 1
+ M_particle_steady[len(M_particle_steady) - 1, 8 + class4] = M_particle_steady[len(M_particle_steady) - 1, 8 + class4] + 1
+ # The class of the considered particle is assigned to the particles in contact
+ M_particle_steady[overlaps[0], 8 + class_i] = M_particle_steady[overlaps[0], 8 + class_i] + 1
+ M_particle_steady[overlaps[1], 8 + class_i] = M_particle_steady[overlaps[1], 8 + class_i] + 1
+ M_particle_steady[overlaps[2], 8 + class_i] = M_particle_steady[overlaps[2], 8 + class_i] + 1
+ M_particle_steady[overlaps[3], 8 + class_i] = M_particle_steady[overlaps[3], 8 + class_i] + 1
+ # Update the progress bar
+ pbar.update(1)
+ # Find NaN values in the steady state array
+ del_list = np.argwhere(np.isnan(M_particle_steady[:, 1]))
+ print('Number of particles that could not be placed:', len(del_list))
+ print('Diameter of particles that could not be placed:', M_particle_steady[del_list, 5])
+ # Delete NaN values from the steady state array
+ M_particle_steady = np.delete(M_particle_steady, del_list, 0)
+ # Number of set particles
+ n_particles_placed = len(M_particle_steady)
+ """ Check the whole particle packing for particles outside the simulation room boarders. """
+ # Mask the values between the bounds in x direction
+ mask_x = (M_particle_steady[:, 3] <= A_0) & (M_particle_steady[:, 3] >= -A_0)
+ # Mask the values between the bounds in y direction
+ mask_y = (M_particle_steady[:, 2] <= A_0) & (M_particle_steady[:, 2] >= -A_0)
+ # Combine the results of both tests
+ mask_xy = [False if not mask_x[i] else mask_y[i] for i in range(len(mask_x))]
+ # Extract particles within the compartiment
+ M_particle_filtered = M_particle_steady[mask_xy]
+ # Number of particles after test
+ n_particles_filtered = len(M_particle_filtered)
+ # Difference between before and after test
+ n_particles_delta = n_particles_placed - n_particles_filtered
+ print('Number of particles that were placed outside simulation boarders (deleted):', n_particles_delta)
+ """ Check the whole particle packing for overlappings. """
+ if final_test == True:
+ # Counting variable for overlaps
+ n_overlaps_total = 0
+ # Consider each particles in the packing
+ for i in range(len(M_particle_filtered)):
+ # Skip NaN values
+ if np.isnan(M_particle_filtered[i, 1]):
+ continue
+ # Retrieve information from the steady state array
+ # Coordinates
+ z1 = float(M_particle_filtered[i, 1])
+ y1 = float(M_particle_filtered[i, 2])
+ x1 = float(M_particle_filtered[i, 3])
+ # Diameter
+ d1 = float(M_particle_filtered[i, 5])
+ # Define particle under consideration
+ particle1 = z1, y1, x1, d1
+ """ Define region of interrest """
+ z_min = z1 - test_range
+ if z_min < 0:
+ z_min = 0
+ z_max = z1 + test_range
+ y_min = y1 - test_range
+ y_max = y1 + test_range
+ x_min = x1 - test_range
+ x_max = x1 + test_range
+ """ Checking particles within the testing range for possible overlappings """
+ # Arrays in x, y and z direction with indices of particles, which could overlap with the considered particle
+ x_array = np.argwhere(np.logical_and(M_particle_filtered[:, 3] > x_min, M_particle_filtered[:, 3] < x_max))
+ y_array = np.argwhere(np.logical_and(M_particle_filtered[:, 2] > y_min, M_particle_filtered[:, 2] < y_max))
+ z_array = np.argwhere(np.logical_and(M_particle_filtered[:, 1] > z_min, M_particle_filtered[:, 1] < z_max))
+ # Combine the arrays --> Particles to check
+ particle_to_check = reduce(np.intersect1d, (x_array, y_array, z_array))
+ # Convert array to list
+ particle_to_check = np.array(particle_to_check).tolist()
+ """ Check each particle that may overlap for an actual overlapping """
+ for j in range(len(particle_to_check)):
+ if i != j:
+ # Retrieve information from the steady state array
+ # Coordinates
+ z2 = float(M_particle_filtered[j, 1])
+ y2 = float(M_particle_filtered[j, 2])
+ x2 = float(M_particle_filtered[j, 3])
+ # Diameter
+ d2 = float(M_particle_filtered[j, 5])
+ # Define particle
+ particle2 = z2, y2, x2, d2
+ # Required min distance between particles
+ dist_required = (d1 + d2) / 2
+ # Actual distance between particles
+ dist_actual = ((x2 - x1) ** 2 + (y2 - y1) ** 2 + (z2 - z1) ** 2) ** 0.5
+ # If actual distanze is below the min distance, they overlap
+ if (dist_required - dist_actual) > 1e-9:
+ # Update counting variable
+ n_overlaps_total = n_overlaps_total + 1
+ # Since every particle is checked with every particle, the overlaps are counted double
+ n_overlaps_total = int(n_overlaps_total / 2)
+ print('Total overlaps after final check:', n_overlaps_total)
+ print('Random loose packing simulation time:')
+ print("--- %s seconds ---" % (time.time() - start_time_rlp))
+ # Time stamp
+ start_time_growing = time.time()
+ """ Delete upper part of particle packing with lower hold-up than average """
+ # z-coordinate of the top particle
+ z_particle_max = np.max(M_particle_filtered[:, 1])
+ # Number of testing compartiments
+ n_compartiments = int(round(z_particle_max / (D_32/2), 0))
+ # Testing compartment height
+ delta_z = z_particle_max / n_compartiments
+ # Volume of testing compartiment
+ volume_compartiment = delta_z * L_0**2
+ # Predefine a list to calculate compartiment holdups
+ holdup_compartiments = np.zeros(n_compartiments)
+ # Iterate through each comartiment and calculate its holdup
+ for i in range(n_compartiments):
+ # Mask the values between the defined compartiment bounds
+ mask_holdup = (M_particle_filtered[:, 1] <= ((n_compartiments-i) * delta_z)) & (M_particle_filtered[:, 1] >= ((n_compartiments-i-1) * delta_z))
+ # Extract particles within the compartiment
+ M_particle_compartiment = M_particle_filtered[mask_holdup]
+ # Volume of particles within the compartiment
+ volume_particles_holdup = np.sum((np.pi / 6) * M_particle_compartiment[:, 5]**3)
+ # Compartiment holdup
+ holdup_compartiment_i = volume_particles_holdup / volume_compartiment
+ # Update list of compartiment holdups
+ holdup_compartiments[i] = holdup_compartiment_i
+ # Upper boandary of average holdup calculation
+ bound_up_holdup = int(round(n_compartiments * 0.3, 0))
+ # Average holdup in packing
+ holdup_average = np.mean(holdup_compartiments[bound_up_holdup:n_compartiments-1])
+ for holdup_boandary in range(n_compartiments):
+ if abs(holdup_average - holdup_compartiments[holdup_boandary]) < 0.05:
+ break
+ # Max. compartiment height
+ z_compartiment_max = (n_compartiments - holdup_boandary) * delta_z
+ # Extract final particles
+ M_particle_final = M_particle_filtered[M_particle_filtered[:, 1] <= z_compartiment_max]
+ # Final height of packing
+ z_packing = max(M_particle_final[:,1] + M_particle_final[:,5] / 2)
+ # Final number of particles in packing
+ n_particle_final_rlp = len(M_particle_final)
+ print('Final number of particles after deleting the top part of the packing with low holdup:', n_particle_final_rlp)
+ """ Graphically display particles in 3D """
+ if display_particles == True:
+ plot_particles(M_particle_final)
+ """
+ Tesselation.
+ From here the coordinates order is x,y,z since the tesselation package uses this order.
+ Perform Tesselation and get information from Tesselation.
+ The cells are polyhedrons.
+ The lists are defined as follows:
+ volumes_cells: cell volumes
+ cell_neighbors: cell neighbors
+ cell_face_normals: normal vectors of cell faces
+ cell_vertices: cell verticies
+ vertices_indices: indices of the vertices of each cell face
+ """
+ sys.stdout.write("Tesselation...")
+ # Starting coordinates for Tesselation, x, y, z
+ tess_start = (-A_0, -A_0, 0)
+ # Final coordinates for Tesselation, x, y, z
+ tess_end = (A_0, A_0, z_packing)
+ # Create an array with coordinates for tesselation
+ M_particle_tess = M_particle_final[:, [3, 2, 1]]
+ # Initial particle diameters
+ diameters_0 = M_particle_final[:,5]
+ # Calculate tesselation
+ volumes_cells, cell_neighbors, cell_face_normals, cell_vertices, vertices_indices = voro_tess(M_particle_tess, tess_start, tess_end, diameters_0)
+ sys.stdout.write("completed\n")
+ # Convert list to array
+ volumes_cells = np.array(volumes_cells)
+ """ Calculate initial/max distances between particle center and polyhedron faces
+ and determine cell faces which are from beginning in contact with other cells due to initial particle contact """
+ sys.stdout.write("Determining initial particle conditions...")
+ # Predefine list for initial/max distances between particle center and polyhedron faces
+ distances_faces_max = [
+ [
+ [None for _ in sub_inner_tuple]
+ for sub_inner_tuple in inner_tuple
+ ]
+ for inner_tuple in cell_face_normals
+ ]
+ # Predefine list for cell faces which are from beginning in contact with other cells due to initial particle contact
+ cell_face_fixed = [
+ [
+ [None for _ in sub_inner_tuple]
+ for sub_inner_tuple in inner_tuple
+ ]
+ for inner_tuple in cell_face_normals
+ ]
+ # Check each particle from particles array
+ for particle_i in range(len(M_particle_tess)):
+ # List of fixed cell faces
+ cell_face_fixed_i = []
+ # List of distances between particle center and faces at each iteration step for volume grow
+ distances_faces_max_i = []
+ # Particle center coordinates
+ center_particle_i = M_particle_tess[particle_i]
+ # Initial particle radius
+ radius_particle_i_0 = diameters_0[particle_i] / 2
+ # Verticies of the cell corresponding to the particle
+ cell_vertices_i = cell_vertices[particle_i]
+ # Indices of the vertices of each cell face corresponding to the particle.
+ vertices_indices_i = vertices_indices[particle_i]
+ # Normal vectors of cell faces corresponding to the particle
+ cell_face_normals_i = cell_face_normals[particle_i]
+ # Iterate through each face of the cell
+ for face in range(len(cell_face_normals_i)):
+ # Indices of the vertices of the considered cell face
+ indices_cell_face = vertices_indices_i[face]
+ # Vector between the particle center and one point on the face area (vertex)
+ vector_center_vertex = center_particle_i - cell_vertices_i[indices_cell_face[0]]
+ # Distance between particle/cell center and the cell face
+ distance_face = np.abs(np.dot(vector_center_vertex, cell_face_normals_i[face]) / np.linalg.norm(cell_face_normals_i[face]))
+ # If the distance is less than or equal to the radius of the sphere, the face intersects with the sphere.
+ # Thus it is set to fixed
+ if np.abs(np.abs(distance_face) - radius_particle_i_0) < 1e-9:
+ face_fixed = True
+ else:
+ face_fixed = False
+ # Update lists
+ distances_faces_max_i.append(distance_face)
+ cell_face_fixed_i.append(face_fixed)
+ # Update lists
+ distances_faces_max[particle_i] = distances_faces_max_i
+ cell_face_fixed[particle_i] = cell_face_fixed_i
+ """ Calculate the contact probability of the initial state """
+ # Create a new list with cell neighbor diameters of each particle
+ # None values are walls
+ diameters_neighbors = [
+ [diameters_0[i] if i >= 0 else None for i in sublist]
+ for sublist in cell_neighbors
+ ]
+ # List to label particles. Label is False if particle is in contact with walls at its initial state
+ # or if the solver found no solution during the growing simulation
+ particle_status_0 = list([True] * len(diameters_0))
+ # Check conditions and update particle_status list accordingly
+ # If None values (walls) of the diameters_neighbors list are at the same position as True values (fixed faces) of the
+ # cell_face_fixed list, the particle are in contact with walls at their initial state. Thus their label is set to False
+ for i in range(len(particle_status_0)):
+ for j in range(len(diameters_neighbors[i])):
+ if diameters_neighbors[i][j] is None and cell_face_fixed[i][j]:
+ particle_status_0[i] = False
+ # Number of particles not in contact with walls
+ n_particle_nowall = int(np.sum(particle_status_0))
+ # Volume of cells not in contact with walls
+ volumes_cells_nowall = volumes_cells[particle_status_0]
+ # Initial diameters of particles not in contact with walls
+ diameters_nowall_0 = diameters_0[particle_status_0]
+ # Initial volumes of particles not in contact with walls
+ volumes_particles_nowall_0 = (np.pi / 6) * diameters_nowall_0 **(3)
+ # Grow factor. This factor is defined by the lowest cell volume to particle volume ratio
+ # since the volume of each particle can be increased max. to this value due to the required constant size distribution
+ # Only particles/cells not in contact with walls are considered
+ # Calculate the grow factor only where the particle_status_0 is True
+ grow_factor_limit = round(np.min((volumes_cells_nowall / volumes_particles_nowall_0) ** (1/3)), 2)
+ grow_factor_max = round(np.max((volumes_cells_nowall / volumes_particles_nowall_0) ** (1/3)), 2)
+ if grow_factor_limit <= 1.0:
+ print('Fatal error: Tesselation failed. Try another seed')
+ sys.exit()
+ # Volume of cube surrounding the packing
+ volume_cube = z_packing * L_0**2
+ # Final particle diameters at constant growing
+ diameters_final = diameters_0 * grow_factor_limit
+ # Final particle volumes at constant growing
+ volumes_particles_final = (np.pi / 6) * diameters_final **(3)
+ # Correct volumes
+ for i in range(len(volumes_particles_final)):
+ if volumes_particles_final[i] - volumes_cells[i] > 1e-9:
+ volumes_particles_final[i] = volumes_cells[i]
+ # Packing final holdup at constant growing
+ holdup_final = np.sum(volumes_particles_final) / volume_cube
+ # Starting holdup for grow simulation
+ holdup_start = 0.525
+ # Grow steps at constant dsd size increase
+ grow_steps_const = int(math.floor((holdup_final - holdup_start) / HOLDUP_DELTA))
+ # Grow steps at non-constant dsd size increase
+ grow_steps_nonconst = int(math.floor((1.0 - holdup_final) / 0.2))
+ # Total grow steps
+ grow_steps = grow_steps_const + grow_steps_nonconst
+ # Predefine list with grow factors
+ grow_factor_list=[]
+ # Grow factor the simulations starts from
+ for i in range(grow_steps+1):
+ if i <= grow_steps_const:
+ holdup_end = holdup_start + HOLDUP_DELTA * i
+ else:
+ holdup_end = 1.0 - (grow_steps + 1 - i) * 0.2
+ # Initial guess for solver
+ initial_guess = grow_factor_limit
+ # Bounds
+ bnds = [(1, grow_factor_max)]
+ # Solver method
+ solver_method = 'SLSQP'
+ # Perform the optimization
+ result = minimize(solver_grow_factor, initial_guess, args=(holdup_end, volume_cube, diameters_0, volumes_cells), method=solver_method, bounds=bnds, tol=1e-5)
+ # Obtain the value of the minimized function
+ minimized_value = result.fun
+ # Obtain grow factor value
+ grow_factor = result.x[0]
+ grow_factor_list.append(grow_factor)
+ # Predefine list of contact probabilities
+ contact_probabilities = [[[[] for _ in range(N_CLASSES)] for _ in range(N_CLASSES)] for _ in range(grow_steps+4)]
+ # Predefine list of average contacts
+ average_contacts = [[None] * N_CLASSES for _ in range(grow_steps+4)]
+ # Determine average contacts and contact probabilities of initial state and update the list
+ average_contacts, contact_probabilities = calc_contact_probabilities(average_contacts, contact_probabilities, cell_face_fixed, diameters_neighbors, diameters_0, 0, particle_status_0)
+ # Set all faces to fixed which is the case in final state
+ cell_face_fixed_final = [[True for _ in row] for row in cell_face_fixed]
+ # Determine average contacts and contact probabilities of final state and update the list
+ average_contacts, contact_probabilities = calc_contact_probabilities(average_contacts, contact_probabilities, cell_face_fixed_final, diameters_neighbors, diameters_0, grow_steps+3, particle_status_0)
+ # Final distribution after random loose packing simulation
+ q0_rlp, q3_rlp = calc_q(diameters_nowall_0, n_particle_nowall, d_classes)
+ # Equivalent particle diameters from cell volumes
+ diameters_cells = ((6/np.pi) * volumes_cells)**(1/3)
+ # Diameters of droplets not in contact with walls
+ diameters_cells_nowall = diameters_cells[particle_status_0]
+ # Max. diameter
+ d_cells_max = np.max(diameters_cells_nowall)
+ # Update distribution classes
+ if d_cells_max > d_95:
+ # Final number of classes
+ n_classes_final = math.ceil((d_cells_max - d_05) / delta_class)
+ # Final class diameters
+ d_classes_final = np.array([d_05+delta_class/2 + i * delta_class for i in range(n_classes_final)])
+ else:
+ d_classes_final = d_classes
+ # Final distribution after growing simulation
+ q0_final, q3_final = calc_q(diameters_cells_nowall, n_particle_nowall, d_classes_final)
+ # Particle initial volumes
+ volumes_particles = (np.pi / 6) * diameters_0 **(3)
+ # Total particle volume
+ volume_particle_total = np.sum(volumes_particles)
+ # Predefine list of packing holdups
+ holdup_packing = np.zeros(grow_steps+4)
+ # Packing initial holdup
+ holdup_packing[0] = volume_particle_total / volume_cube
+ # Packing final holdup
+ holdup_packing[grow_steps+3] = 1.0
+ #
+ sys.stdout.write("completed\n")
+ # List to label particles. Label is False if particle's volume is smaller than its cell volume.
+ cell_status = list([False] * len(diameters_0))
+ # For loop to simulate growing
+ for i_grow in range(grow_steps + 2):
+ print('Particle grow simulation. Step:', i_grow+1, 'of', (grow_steps + 2))
+ # Increase diameters
+ if i_grow <= grow_steps_const:
+ diameters = diameters_0 * grow_factor_list[i_grow]
+ elif i_grow == (grow_steps_const + 1):
+ diameters = diameters_0 * grow_factor_limit
+ else:
+ diameters = diameters_0 * grow_factor_list[i_grow-1]
+ # Update volumes
+ volumes_particles = (np.pi / 6) * diameters **(3)
+ # Check and update cell status
+ for i in range(len(cell_status)):
+ if volumes_particles[i] - volumes_cells[i] > 1e-9:
+ volumes_particles[i] = volumes_cells[i]
+ cell_status[i] = True
+ # Final total particle volume
+ volume_particle_total = np.sum(volumes_particles)
+ # Update packing holdup
+ holdup_packing[i_grow+1] = volume_particle_total / volume_cube
+ # List of numbers to be squared
+ particle_indices = list(range(0, len(M_particle_tess)))
+ sys.stdout.write("Initializing multiprocessing...")
+ # Initialize manager for shared data
+ with multiprocessing.Manager() as manager:
+ # Lists shared with manager
+ cell_face_fixed = manager.list([manager.list(sublist) for sublist in cell_face_fixed])
+ particle_status = manager.list(particle_status_0)
+ # Number of processes to use
+ if N_CORES == 0:
+ num_processes = multiprocessing.cpu_count()
+ else:
+ num_processes = N_CORES
+ # Create a pool of worker processes
+ with multiprocessing.Pool(processes=num_processes) as pool:
+ # Create partial function for parallel processing
+ calc_partial = partial(calc_grow,
+ volumes_particles = volumes_particles,
+ diameters = diameters,
+ cell_face_fixed = cell_face_fixed,
+ distances_faces_max = distances_faces_max,
+ cell_face_normals = cell_face_normals,
+ vertices_indices = vertices_indices,
+ cell_vertices = cell_vertices,
+ display_particle_cell = display_particle_cell,
+ particle_status = particle_status,
+ cell_status = cell_status)
+ # List to keep track of AsyncResult objects
+ async_results = []
+ # Distribute the tasks among the worker processes using apply_async
+ for index in particle_indices:
+ result = pool.apply_async(calc_partial, (index,))
+ async_results.append(result)
+ sys.stdout.write("completed\n")
+ # Use tqdm to track progress
+ with tqdm(total=len(async_results), desc='Particle grow simulation progress', unit='particle') as pbar:
+ results = []
+ for result in async_results:
+ results.append(result.get())
+ # Update progress bar after each task completes
+ pbar.update(1)
+ # No need to close the pool and wait, 'with' statement handles it
+ # Convert shared lists to regular lists for further processing
+ cell_face_fixed = [list(sublist) for sublist in cell_face_fixed]
+ particle_status = list(particle_status)
+ # Determine contact probabilities and update the list
+ average_contacts, contact_probabilities = calc_contact_probabilities(average_contacts, contact_probabilities, cell_face_fixed, diameters_neighbors, diameters_0, i_grow+1, particle_status)
+ # Final number of particles considered for the growing simulation
+ n_particle_final_gs = int(np.sum(particle_status))
+ """ Save results as Excel sheet. """
+ # Convert the arrays into dataframes
+ holdup_packing_df = pd.DataFrame(100*holdup_packing)
+ average_contacts_df = pd.DataFrame(average_contacts)
+ contact_probabilities_df = pd.DataFrame(contact_probabilities)
+ q0_initial_df = pd.DataFrame(100*q0_initial)
+ q0_rlp_df = pd.DataFrame(100*q0_rlp)
+ q0_final_df = pd.DataFrame(100*q0_final)
+ q3_initial_df = pd.DataFrame(100*q3_initial)
+ q3_rlp_df = pd.DataFrame(100*q3_rlp)
+ q3_final_df = pd.DataFrame(100*q3_final)
+ # Round class diameters
+ np.round(d_classes_final, decimals=0, out=d_classes_final)
+ d_classes_final_df = pd.DataFrame(d_classes_final)
+ # Define headers
+ headers_simulation = ['Seed', 'Final number of particles']
+ headers_holdup_packing = ['Grow step', 'Hold-up']
+ headers_d_classes = ['Class', 'Class diameter']
+ header_average_contacts = ['Grow step']
+ headers_dsd = ['q0_initial', 'q0_rlp', 'q0_final', 'q3_initial', 'q3_rlp', 'q3_final']
+ for i in range(N_CLASSES):
+ header_average_contacts.append(f'Class {i+1}')
+ header_contact_probabilities = []
+ for i in range(N_CLASSES):
+ header_contact_probabilities.append(f'Class {i+1}')
+ # Create a new Excel workbook
+ workbook = xlwt.Workbook()
+ # Create a new Excel sheet for the holdup data
+ sheet_main = workbook.add_sheet('Main')
+ # Write simulation settings headers to Excel
+ for row_index, header in enumerate(headers_simulation):
+ sheet_main.write(row_index, 0, header)
+ # Write simulation settings data to Excel
+ sheet_main.write(0, 1, SEED)
+ sheet_main.write(1, 1, n_particle_final_rlp)
+ sheet_main.write(1, 2, n_particle_nowall)
+ sheet_main.write(1, 3, n_particle_final_gs)
+ # Write holdup headers to Excel
+ for row_index, header in enumerate(headers_holdup_packing):
+ sheet_main.write(row_index + 3, 0, header)
+ # Write holdup data to Excel
+ for col_index, value in enumerate(holdup_packing_df.values):
+ sheet_main.write(3, col_index + 1, col_index)
+ sheet_main.write(4, col_index + 1, value[0])
+ # Write classes headers to Excel
+ for row_index, header in enumerate(headers_d_classes):
+ sheet_main.write(row_index + 6, 0, header)
+ # Write classes data to Excel
+ for col_index in range(len(d_classes_final)):
+ sheet_main.write(6, col_index + 1, col_index + 1)
+ for col_index, value in enumerate(d_classes_final_df.values):
+ sheet_main.write(7, col_index + 1, value[0])
+ # Write DSD headers to Excel
+ for row_index, header in enumerate(headers_dsd):
+ sheet_main.write(row_index + 9, 0, header)
+ # Write DSD data to Excel
+ for col_index, value in enumerate(q0_initial_df.values):
+ sheet_main.write(9, col_index + 1, value[0])
+ for col_index, value in enumerate(q0_rlp_df.values):
+ sheet_main.write(10, col_index + 1, value[0])
+ for col_index, value in enumerate(q0_final_df.values):
+ sheet_main.write(11, col_index + 1, value[0])
+ for col_index, value in enumerate(q3_initial_df.values):
+ sheet_main.write(12, col_index + 1, value[0])
+ for col_index, value in enumerate(q3_rlp_df.values):
+ sheet_main.write(13, col_index + 1, value[0])
+ for col_index, value in enumerate(q3_final_df.values):
+ sheet_main.write(14, col_index + 1, value[0])
+ # Create a new Excel sheet for the contacts data
+ sheet_contacts = workbook.add_sheet('Contacts')
+ # Write contacts headers to Excel
+ for col_index, header in enumerate(header_average_contacts):
+ sheet_contacts.write(0, col_index, header)
+ # Write contacts data to Excel
+ for row_index, values in enumerate(average_contacts_df.values):
+ sheet_contacts.write(row_index + 1, 0, row_index)
+ for col_index, value in enumerate(values):
+ sheet_contacts.write(row_index + 1, col_index + 1, value)
+ # Iterate through the nested dataframe of contact probabilities
+ for sheet_index, sheet_data in enumerate(contact_probabilities_df.values):
+ # Create a new Excel sheet for each inner list
+ sheet_probabilities = workbook.add_sheet(f'Probabilities Step {sheet_index}')
+ # Write probabilities headers to Excel
+ for col_index, header in enumerate(header_contact_probabilities):
+ sheet_probabilities.write(0, col_index + 1, header)
+ sheet_probabilities.write(col_index + 1, 0, header)
+ # Iterate through each sublist
+ for row_index, sublist in enumerate(sheet_data):
+ # Write contact probabilities to Excel
+ for col_index, value in enumerate(sublist):
+ sheet_probabilities.write(row_index + 1, col_index + 1, value)
+ # Save the workbook to a file
+ file_name = os.path.join(ROOT_DIR, NAME_EXCEL_FILE)
+ workbook.save(file_name)
+ """ Graphically display particles in 3D """
+ if display_particles == True:
+ plot_particles(M_particle_final)
+ print('Particle growing simulation time:')
+ print("--- %s seconds ---" % (time.time() - start_time_growing))
+ print('Total simulation time:')
+ print("--- %s seconds ---" % (time.time() - start_time_rlp))
+ sss = 1
\ No newline at end of file