From 3f7c55aab7eb30c82ce6220b04003304faa351c6 Mon Sep 17 00:00:00 2001
From: Giovanni La Mura <giovanni.lamura@inaf.it>
Date: Tue, 8 Apr 2025 18:24:10 +0200
Subject: [PATCH] Implement random generator and include it in model_maker
---
src/scripts/generator.py | 36 +++++----
src/scripts/model_maker.py | 150 +++++++++++++++++++++++++++++++++----
2 files changed, 160 insertions(+), 26 deletions(-)
diff --git a/src/scripts/generator.py b/src/scripts/generator.py
index 97491a5f..901f58dd 100644
--- a/src/scripts/generator.py
+++ b/src/scripts/generator.py
@@ -12,7 +12,7 @@ seed = int(argv[1])
random.seed(seed)
config = {}
-with open('config.yml', 'r') as stream:
+with open('config_rand.yml', 'r') as stream:
config = yaml.safe_load(stream)
nsph = int(config['particle_settings']['n_spheres'])
@@ -20,11 +20,19 @@ nsph = int(config['particle_settings']['n_spheres'])
vec_thetas = [0.0 for i in range(nsph)]
vec_phis = [0.0 for i in range(nsph)]
vec_rads = [0.0 for i in range(nsph)]
-vec_sph_x = [0.0 for i in range(nsph)]
-vec_sph_y = [0.0 for i in range(nsph)]
-vec_sph_z = [0.0 for i in range(nsph)]
+n_types = config['particle_settings']['n_types']
+if (len(config['particle_settings']['sph_types']) > 0):
+ if (len(config['particle_settings']['sph_types']) != nsph):
+ print("ERROR: declared number of types does not match the number of spheres!")
+ exit(1)
+else:
+ # Random type generation
+ for ti in range(nsph):
+ itype = 1 + int(n_types * random.random())
+ config['particle_settings']['sph_types'].append(itype)
sph_type_index = config['particle_settings']['sph_types'][0] - 1
+vec_spheres = [{'itype': sph_type_index + 1, 'x': 0.0, 'y': 0.0, 'z': 0.0}]
vec_rads[0] = config['particle_settings']['radii'][sph_type_index]
max_rad = 20.0 * vec_rads[0]
placed_spheres = 1
@@ -48,9 +56,9 @@ for i in range(1, nsph):
j = 0
while (j < i - 1):
j += 1
- dx2 = (x - vec_sph_x[j]) * (x - vec_sph_x[j])
- dy2 = (y - vec_sph_y[j]) * (y - vec_sph_y[j])
- dz2 = (z - vec_sph_z[j]) * (z - vec_sph_z[j])
+ dx2 = (x - vec_spheres[j]['x']) * (x - vec_spheres[j]['x'])
+ dy2 = (y - vec_spheres[j]['y']) * (y - vec_spheres[j]['y'])
+ dz2 = (z - vec_spheres[j]['z']) * (z - vec_spheres[j]['z'])
dist2 = dx2 + dy2 + dz2
rr2 = (vec_rads[i] + vec_rads[j]) * (vec_rads[i] + vec_rads[j])
if (dist2 < rr2):
@@ -83,13 +91,15 @@ for i in range(1, nsph):
# The current direction is filled. Try another one.
attempts += 1
continue # while(not is_placed)
- vec_sph_x[i] = x
- vec_sph_y[i] = y
- vec_sph_z[i] = z
+ vec_spheres.append({
+ 'itype': sph_type_index + 1,
+ 'x': x,
+ 'y': y,
+ 'z': z
+ })
is_placed = True
placed_spheres += 1
attempts = 0
-print(vec_sph_x)
-print(vec_sph_y)
-print(vec_sph_z)
+print(vec_spheres)
+print(sorted(vec_spheres, key=lambda item: item['itype']))
diff --git a/src/scripts/model_maker.py b/src/scripts/model_maker.py
index 21a2eb1b..70b4df23 100755
--- a/src/scripts/model_maker.py
+++ b/src/scripts/model_maker.py
@@ -23,8 +23,17 @@
#
# The script requires python3.
+import math
+import random
import yaml
+allow_3d = True
+try:
+ import pyvista as pv
+except ModuleNotFound as ex:
+ print("WARNING: pyvista not found!")
+ allow_3d = False
+
from pathlib import PurePath
from sys import argv
@@ -131,6 +140,10 @@ def load_model(model_file):
except FileNotFoundError:
print("ERROR: " + model_file + " was not found!")
if model is not None:
+ make_3d = False if model['system_settings']['make_3D'] == "0" else True
+ if (make_3d and not allow_3d):
+ print("WARNING: 3D visualization of models is not available. Disabling.")
+ make_3d = False
# Create the sconf dict
sconf = {
'out_file': PurePath(
@@ -252,11 +265,21 @@ def load_model(model_file):
print("ERROR: %s is not a supported scaling mode!"%(model['radiation_settings']['scale_name']))
return (None, None)
sph_types = model['particle_settings']['sph_types']
- if (len(sph_types) != sconf['nsph'] and len(sph_types) != 0):
- print("ERROR: vector of sphere types does not match the declared number of spheres!")
- return (None, None)
- else:
+ if (len(sph_types) == sconf['nsph']):
sconf['vec_types'] = [int(str_typ) for str_typ in sph_types]
+ else:
+ if (len(sph_types) != 0):
+ print("ERROR: vector of sphere types does not match the declared number of spheres!")
+ return (None, None)
+ else: # len(sph_types) = 0
+ len_vec_x = len(model['geometry_settings']['x_coords'])
+ len_vec_y = len(model['geometry_settings']['y_coords'])
+ len_vec_z = len(model['geometry_settings']['z_coords'])
+ if (len_vec_x != 0 or len_vec_y != 0 or len_vec_z != 0):
+ print("ERROR: cannot assign random types with explicit sphere positions!")
+ return (None, None)
+ else:
+ sconf['vec_types'] = [0 for ti in range(sconf['nsph'])]
if (len(model['particle_settings']['radii']) != sconf['configurations']):
print("ERROR: Declared number of radii does not match number of types!")
return (None, None)
@@ -324,22 +347,24 @@ def load_model(model_file):
gconf['vec_sph_y'] = [0.0 for i in range(gconf['nsph'])]
gconf['vec_sph_z'] = [0.0 for i in range(gconf['nsph'])]
if (gconf['application'] != "SPHERE" or gconf['nsph'] != 1):
- if (len(model['geometry_settings']['x_coords']) != len(model['geometry_settings']['y_coords'])):
+ len_vec_x = len(model['geometry_settings']['x_coords'])
+ len_vec_y = len(model['geometry_settings']['y_coords'])
+ len_vec_z = len(model['geometry_settings']['z_coords'])
+ if (len_vec_x != len_vec_y):
print("ERROR: X and Y coordinate vectors have different lengths!")
return (None, None)
- if (len(model['geometry_settings']['x_coords']) != len(model['geometry_settings']['z_coords'])):
+ if (len_vec_x != len_vec_z):
print("ERROR: X and Z coordinate vectors have different lengths!")
return (None, None)
- if (len(model['geometry_settings']['y_coords']) != len(model['geometry_settings']['z_coords'])):
+ if (len_vec_y != len_vec_z):
print("ERROR: Y and Z coordinate vectors have different lengths!")
return (None, None)
- if (len(model['particle_settings']['sph_types']) == 0 and len(model['geometry_settings']['x_coords']) != 0):
- print("ERROR: random type assignment is not allowed with assigned coordinates!")
- return (None, None)
- # TODO: put here a test to allow for empty vectors in case of random generation
- if (len(model['geometry_settings']['x_coords']) == 0):
+ if (len_vec_x == 0):
# Generate random cluster
- print("DEBUG: would generate random cluster.")
+ rnd_seed = int(model['system_settings']['rnd_seed'])
+ # random_aggregate() checks internally whether application is INCLUSION
+ print("DEBUG: make_3d is %s."%str(make_3d))
+ random_aggregate(sconf, gconf, rnd_seed, 20.0 * max(sconf['ros']), make_3d)
else:
if (len(model['geometry_settings']['x_coords']) != gconf['nsph']):
print("ERROR: coordinate vectors do not match the number of spheres!")
@@ -475,6 +500,105 @@ def print_help():
print("--help Print this help and exit.")
print(" ")
+# random_aggregate()
+def random_aggregate(scatterer, geometry, seed, max_rad, make3d=False, max_attempts=100):
+ random.seed(seed)
+ nsph = scatterer['nsph']
+ vec_thetas = [0.0 for i in range(nsph)]
+ vec_phis = [0.0 for i in range(nsph)]
+ vec_rads = [0.0 for i in range(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
+ sph_type_index = scatterer['vec_types'][0] - 1
+ vec_spheres = [{'itype': sph_type_index + 1, 'x': 0.0, 'y': 0.0, 'z': 0.0}]
+ vec_rads[0] = scatterer['ros'][sph_type_index]
+ placed_spheres = 1
+ attempts = 0
+ for i in range(1, nsph):
+ sph_type_index = scatterer['vec_types'][i] - 1
+ vec_rads[i] = scatterer['ros'][sph_type_index]
+ is_placed = False
+ while (not is_placed):
+ if (attempts > max_attempts):
+ print("WARNING: could not place sphere %d in allowed radius!"%i)
+ break # while(not is_placed)
+ vec_thetas[i] = math.pi * random.random()
+ vec_phis[i] = 2.0 * math.pi * random.random()
+ rho = vec_rads[0] + vec_rads[i]
+ z = rho * math.cos(vec_thetas[i])
+ y = rho * math.sin(vec_thetas[i]) * math.sin(vec_phis[i])
+ x = rho * math.sin(vec_thetas[i]) * math.cos(vec_phis[i])
+ j = 0
+ while (j < i - 1):
+ j += 1
+ dx2 = (x - vec_spheres[j]['x']) * (x - vec_spheres[j]['x'])
+ dy2 = (y - vec_spheres[j]['y']) * (y - vec_spheres[j]['y'])
+ dz2 = (z - vec_spheres[j]['z']) * (z - vec_spheres[j]['z'])
+ dist2 = dx2 + dy2 + dz2
+ rr2 = (vec_rads[i] + vec_rads[j]) * (vec_rads[i] + vec_rads[j])
+ if (dist2 < rr2):
+ # Spheres i and j are compenetrating.
+ # Sphere i is moved out radially until it becomes externally
+ # tangent to sphere j. Then the check is repeated, to verify
+ # that no other sphere was penetrated. The process is iterated
+ # until sphere i is placed or the maximum allowed radius is
+ # reached.
+ sinthi = math.sin(vec_thetas[i])
+ sinthj = math.sin(vec_thetas[j])
+ costhi = math.cos(vec_thetas[i])
+ costhj = math.cos(vec_thetas[j])
+ sinphi = math.sin(vec_phis[i])
+ sinphj = math.sin(vec_phis[j])
+ cosphi = math.cos(vec_phis[i])
+ cosphj = math.cos(vec_phis[j])
+ cosalpha = (
+ sinthi * cosphi * sinthj * cosphj
+ + sinthi * sinphi * sinthj * sinphj
+ + costhi * costhj
+ )
+ rho += 2.0 * vec_rads[j] * cosalpha
+ z = rho * math.cos(vec_thetas[i])
+ y = rho * math.sin(vec_thetas[i]) * math.sin(vec_phis[i])
+ x = rho * math.sin(vec_thetas[i]) * math.cos(vec_phis[i])
+ j = 0
+ continue # while(j < i - 1)
+ if (rho + vec_rads[i] > max_rad):
+ # The current direction is filled. Try another one.
+ attempts += 1
+ continue # while(not is_placed)
+ vec_spheres.append({
+ 'itype': sph_type_index + 1,
+ 'x': x,
+ 'y': y,
+ 'z': z
+ })
+ is_placed = True
+ placed_spheres += 1
+ attempts = 0
+ 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
+ print("DEBUG: make3d is %s"%str(make3d))
+ if make3d:
+ # Do something
+ for vi in range(len(vec_rads)):
+ radius = vec_rads[vi] / max_rad
+ x = vec_spheres[vi]['x'] / max_rad
+ y = vec_spheres[vi]['y'] / max_rad
+ z = vec_spheres[vi]['z'] / max_rad
+ mesh = pv.Sphere(radius, (x, y, z))
+ mesh_name = "sphere_{0:04d}.obj".format(vi)
+ mesh.save(mesh_name)
+# end random_aggregate()
+
## \brief Write the geometry configuration dictionary to legacy format.
#
# \param conf: `dict` Geometry configuration dictionary.
--
GitLab