diff --git a/DCNumPro_1.0.py b/DCNumPro_1.0.py
deleted file mode 100644
index 9c9917ea33d44dab7e46df67c46ccb0cea659596..0000000000000000000000000000000000000000
--- a/DCNumPro_1.0.py
+++ /dev/null
@@ -1,4132 +0,0 @@
-""" 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