diff --git a/src/scripts/model_maker.py b/src/scripts/model_maker.py
index 2354daff70947df4471de7ffd84d7638f0846f54..7776fdeb668f3dd2969d62f28b2577f297ad5c3e 100755
--- a/src/scripts/model_maker.py
+++ b/src/scripts/model_maker.py
@@ -24,6 +24,7 @@
# The script requires python3.
import math
+import numpy as np
import os
import pdb
import random
@@ -365,9 +366,13 @@ def load_model(model_file):
if (len_vec_x == 0):
# Generate random cluster
rnd_seed = int(model['system_settings']['rnd_seed'])
- # random_aggregate() checks internally whether application is INCLUSION
max_rad = float(model['particle_settings']['max_rad'])
- random_aggregate(sconf, gconf, rnd_seed, max_rad)
+ # random_aggregate() checks internally whether application is INCLUSION
+ #random_aggregate(sconf, gconf, rnd_seed, max_rad)
+ check = random_compact(sconf, gconf, rnd_seed, max_rad)
+ if (check != 0):
+ print("INFO: stopping with exit code %d."%check)
+ exit(check)
else:
if (len(model['geometry_settings']['x_coords']) != gconf['nsph']):
print("ERROR: coordinate vectors do not match the number of spheres!")
@@ -522,7 +527,10 @@ def print_help():
# \param seed: `int` Seed for the random sequence generation
# \param max_rad: `float` Maximum allowed radial extension of the aggregate
# \param max_attempts: `int` Maximum number of attempts to place a particle in any direction
+# \return result: `int` Function exit code (0 for success, otherwise number of
+# spheres that could not be placed)
def random_aggregate(scatterer, geometry, seed, max_rad, max_attempts=100):
+ result = 0
random.seed(seed)
nsph = scatterer['nsph']
vec_thetas = [0.0 for i in range(nsph)]
@@ -545,6 +553,7 @@ def random_aggregate(scatterer, geometry, seed, max_rad, max_attempts=100):
is_placed = False
while (not is_placed):
if (attempts > max_attempts):
+ result += 1
print("WARNING: could not place sphere %d in allowed radius!"%i)
break # while(not is_placed)
vec_thetas[i] = math.pi * random.random()
@@ -618,8 +627,123 @@ def random_aggregate(scatterer, geometry, seed, max_rad, max_attempts=100):
geometry['vec_sph_y'][sph_index] = sphere['y']
geometry['vec_sph_z'][sph_index] = sphere['z']
sph_index += 1
-# end random_aggregate()
+ return result
+## \brief Generate a random compact cluster from YAML configuration options.
+#
+# This function generates a random aggregate of spheres using the maximum
+# compactness packaging to fill a spherical volume with given maximum radius,
+# then it proceeds by subtracting random spheres from the outer layers, until
+# the aggregate is reduced to the desired number of spheres. The function
+# can only be used if all sphere types have the same radius. The result of the
+# generated model is directly saved in the parameters of the scatterer and
+# geometry configuration dictionaries.
+#
+# \param scatterer: `dict` Scatterer configuration dictionary (gets modified)
+# \param geometry: `dict` Geometry configuration dictionary (gets modified)
+# \param seed: `int` Seed for the random sequence generation
+# \param max_rad: `float` Maximum allowed radial extension of the aggregate
+# \return result: `int` Function exit code (0 for success, otherwise error code)
+def random_compact(scatterer, geometry, seed, max_rad):
+ result = 0
+ random.seed(seed)
+ nsph = scatterer['nsph']
+ n_types = scatterer['configurations']
+ if (0 in scatterer['vec_types']):
+ tincrement = 1 if scatterer['application'] != "INCLUSION" else 2
+ for ti in range(nsph):
+ itype = tincrement + int(n_types * random.random())
+ scatterer['vec_types'][ti] = itype
+ if (max(scatterer['ros']) != min(scatterer['ros'])):
+ result = 1
+ else:
+ radius = scatterer['ros'][0]
+ x_centers = np.arange(-1.0 * max_rad + radius, max_rad, 2.0 * radius)
+ y_centers = np.arange(-1.0 * max_rad + radius, max_rad, math.sqrt(3.0) * radius)
+ z_centers = np.arange(-1.0 * max_rad + radius, max_rad, math.sqrt(3.0) * radius)
+ x_offset = radius
+ y_offset = radius / math.sqrt(3.0)
+ z_offset = 0.0
+ tmp_spheres = []
+ n_cells = len(x_centers) * len(y_centers) * len(z_centers)
+ print("INFO: the cubic space would contain %d spheres."%n_cells)
+ n_max_spheres = int(max_rad * max_rad * max_rad / (radius * radius * radius) * 0.74)
+ print("INFO: the maximum radius allows for %d spheres."%n_max_spheres)
+ for zi in range(len(z_centers)):
+ if (y_offset == 0.0):
+ y_offset = radius / math.sqrt(3.0)
+ else:
+ y_offset = 0.0
+ for yi in range(len(y_centers)):
+ if (x_offset == 0.0):
+ x_offset = radius
+ else:
+ x_offset = 0.0
+ for xi in range(len(x_centers)):
+ x = x_centers[xi] + x_offset
+ y = y_centers[yi] + y_offset
+ z = z_centers[zi]
+ extent = radius + math.sqrt(x * x + y * y + z * z)
+ if (extent < max_rad):
+ tmp_spheres.append({
+ 'itype': 1,
+ 'x': x,
+ 'y': y,
+ 'z': z
+ })
+ #tmp_spheres = [{'itype': 1, 'x': 0.0, 'y': 0.0, 'z': 0.0}]
+ current_n = len(tmp_spheres)
+ print("INFO: before erosion there are %d spheres in use."%current_n)
+ rho = 2.0 * max_rad
+ while (current_n > nsph):
+ theta = 2.0 * math.pi * random.random()
+ phi = math.pi * random.random()
+ x0 = rho * math.sin(theta) * math.cos(phi)
+ y0 = rho * math.sin(theta) * math.sin(phi)
+ z0 = rho * math.cos(theta)
+ closest_index = 0
+ minimum_distance = 1000.0 * max_rad
+ for di in range(len(tmp_spheres)):
+ x1 = tmp_spheres[di]['x']
+ if (x1 == max_rad):
+ continue
+ y1 = tmp_spheres[di]['y']
+ z1 = tmp_spheres[di]['z']
+ distance = math.sqrt(
+ (x1 - x0) * (x1 - x0)
+ + (y1 - y0) * (y1 - y0)
+ + (z1 - z0) * (z1 - z0)
+ )
+ if (distance < minimum_distance):
+ closest_index = di
+ minimum_distance = distance
+ tmp_spheres[closest_index]['x'] = max_rad
+ current_n -= 1
+ vec_spheres = []
+ sph_index = 0
+ for ti in range(len(tmp_spheres)):
+ sphere = tmp_spheres[ti]
+ if (sphere['x'] < max_rad):
+ sphere['itype'] = scatterer['vec_types'][sph_index]
+ sph_index += 1
+ vec_spheres.append(sphere)
+ #pl = pv.Plotter()
+ #for si in range(len(vec_spheres)):
+ # x = vec_spheres[si]['x'] / max_rad
+ # y = vec_spheres[si]['y'] / max_rad
+ # z = vec_spheres[si]['z'] / max_rad
+ # mesh = pv.Sphere(radius / max_rad, (x, y, z))
+ # pl.add_mesh(mesh)
+ #pl.export_obj("scene.obj")
+ sph_index = 0
+ for sphere in sorted(vec_spheres, key=lambda item: item['itype']):
+ scatterer['vec_types'][sph_index] = sphere['itype']
+ geometry['vec_sph_x'][sph_index] = sphere['x']
+ geometry['vec_sph_y'][sph_index] = sphere['y']
+ geometry['vec_sph_z'][sph_index] = sphere['z']
+ sph_index += 1
+ return result
+
## \brief Write the geometry configuration dictionary to legacy format.
#
# \param conf: `dict` Geometry configuration dictionary.