Skip to content
Snippets Groups Projects
Commit 9d1d7e4d authored by Giovanni La Mura's avatar Giovanni La Mura
Browse files

Implement first random particle generator

parent 3f7c55aa
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
src/scripts/config_gen.yml deleted 100644 → 0
+ 0
− 93
View file @ 3f7c55aa
system_settings:
max_host_ram : 0
max_gpu_ram : 0
input_settings:
input_folder : "."
spheres_file : "DEDFB"
geometry_file: "DCLU"
output_settings:
output_folder: "."
output_name : "c_OCLU"
formats : [ "LEGACY", "HDF5" ]
jwtm : 1
particle_settings:
application : "CLUSTER"
n_spheres : 8
n_types : 2
sph_types : [ 1, 1, 1, 1, 2, 2, 2, 2 ]
n_layers : [ 1, 1 ]
radii : [ 6.000e-08, 4.000e-8 ]
rad_frac : [ [ 1.0 ], [ 1.0 ] ]
dielec_id : [ [ 1 ], [ 1 ] ]
material_settings:
diel_flag : 0
extern_diel : 1.00e+00
dielec_path : "../../ref_data"
dielec_file : [ "eps_ashok_C.csv" ]
dielec_fmt : [ "CSV" ]
match_mode : "GRID"
diel_const : [ ]
radiation_settings:
polarization: "LINEAR"
scale_start : 1.00e-06
scale_end : 2.00e-05
scale_step : 5.00e-09
wp : 2.99792e+08
xip : 1.00e+00
step_flag : 0
scale_name : "WAVELENGTH"
geometry_settings:
li : 4
le : 8
npnt : 149
npntts : 300
iavm : 0
isam : 0
in_th_start : 0.0
in_th_step : 0.0
in_th_end : 0.0
in_ph_start : 0.0
in_ph_step : 0.0
in_ph_end : 0.0
sc_th_start : 0.0
sc_th_step : 0.0
sc_th_end : 0.0
sc_ph_start : 0.0
sc_ph_step : 0.0
sc_ph_end : 0.0
x_coords : [
8.00e-08,
-8.00e-08,
0.00e+00,
0.00e+00,
0.00e+00,
0.00e+00,
0.00e+00,
0.00e+00
]
y_coords : [
0.00e+00,
0.00e+00,
8.00e-08,
-8.00e-08,
0.00e+00,
0.00e+00,
0.00e+00,
0.00e+00
]
z_coords : [
0.00e+00,
0.00e+00,
0.00e+00,
0.00e+00,
8.00e-08,
-8.00e-08,
16.00e-08,
-16.00e-08
]
src/scripts/generator.py deleted 100644 → 0
+ 0
− 105
View file @ 3f7c55aa
import math
import random
import yaml
#import pdb
from pathlib import PurePath
from sys import argv
seed = int(argv[1])
random.seed(seed)
config = {}
with open('config_rand.yml', 'r') as stream:
config = yaml.safe_load(stream)
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)]
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
attempts = 0
max_attempts = 100
for i in range(1, nsph):
sph_type_index = config['particle_settings']['sph_types'][i] - 1
vec_rads[i] = config['particle_settings']['radii'][sph_type_index]
is_placed = False
#breakpoint()
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
print(vec_spheres)
print(sorted(vec_spheres, key=lambda item: item['itype']))
src/scripts/model_maker.py
+ 101
− 19
View file @ 9d1d7e4d
#!/bin/python3
#!/usr/bin/env python3
# Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari
#
......@@ -24,13 +24,15 @@
# The script requires python3.
import math
import os
import pdb
import random
import yaml
allow_3d = True
try:
import pyvista as pv
except ModuleNotFound as ex:
except ModuleNotFoundError as ex:
print("WARNING: pyvista not found!")
allow_3d = False
......@@ -140,6 +142,7 @@ def load_model(model_file):
except FileNotFoundError:
print("ERROR: " + model_file + " was not found!")
if model is not None:
max_rad = 0.0
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.")
......@@ -363,8 +366,8 @@ def load_model(model_file):
# 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)
max_rad = float(model['particle_settings']['max_rad'])
random_aggregate(sconf, gconf, rnd_seed, max_rad)
else:
if (len(model['geometry_settings']['x_coords']) != gconf['nsph']):
print("ERROR: coordinate vectors do not match the number of spheres!")
......@@ -373,6 +376,11 @@ def load_model(model_file):
gconf['vec_sph_x'][si] = float(model['geometry_settings']['x_coords'][si])
gconf['vec_sph_y'][si] = float(model['geometry_settings']['y_coords'][si])
gconf['vec_sph_z'][si] = float(model['geometry_settings']['z_coords'][si])
#
if (model['system_settings']['make_3D'] != "0" and allow_3d):
if (max_rad == 0.0):
max_rad = 20.0 * max(sconf['ros'])
write_obj(sconf, gconf, max_rad)
else: # model is None
print("ERROR: could not parse " + model_file + "!")
return (sconf, gconf)
......@@ -500,8 +508,21 @@ 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):
## \brief Generate a random cluster aggregate from YAML configuration options.
#
# This function generates a random aggregate of spheres using radial ejection
# in random directions of new spheres until they become tangent to the
# outermost sphere existing in that direction. 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
# \param max_attempts: `int` Maximum number of attempts to place a particle in any direction
def random_aggregate(scatterer, geometry, seed, max_rad, max_attempts=100):
random.seed(seed)
nsph = scatterer['nsph']
vec_thetas = [0.0 for i in range(nsph)]
......@@ -540,13 +561,14 @@ def random_aggregate(scatterer, geometry, seed, max_rad, make3d=False, max_attem
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):
if (dist2 < 0.9999 * 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.
# breakpoint()
sinthi = math.sin(vec_thetas[i])
sinthj = math.sin(vec_thetas[j])
costhi = math.cos(vec_thetas[i])
......@@ -560,7 +582,17 @@ def random_aggregate(scatterer, geometry, seed, max_rad, make3d=False, max_attem
+ sinthi * sinphi * sinthj * sinphj
+ costhi * costhj
)
rho += 2.0 * vec_rads[j] * cosalpha
D12 = math.sqrt(
vec_spheres[j]['x'] * vec_spheres[j]['x']
+ vec_spheres[j]['y'] * vec_spheres[j]['y']
+ vec_spheres[j]['z'] * vec_spheres[j]['z']
)
O1K = D12 * cosalpha
sinalpha = math.sqrt(1.0 - cosalpha * cosalpha)
sinbetaprime = D12 / (vec_rads[i] + vec_rads[j]) * sinalpha
cosbetaprime = math.sqrt(1.0 - sinbetaprime * sinbetaprime)
Op3K = (vec_rads[i] + vec_rads[j]) * cosbetaprime
rho = O1K + Op3K
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])
......@@ -586,17 +618,6 @@ def random_aggregate(scatterer, geometry, seed, max_rad, make3d=False, max_attem
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.
......@@ -749,6 +770,67 @@ def write_legacy_sconf(conf):
output.close()
return result
## \brief Export the model to a set of OBJ files for 3D visualization.
#
# This function exports the model as a set of OBJ files (one for every
# spherical unit, plus a single scene file) to allow for model visualization
# with 3D software tools.
#
# The spherical units are saved as `sphere_XXXX.obj` files. The material
# information is collected in a `model.mtl` file and the whole scene is written
# to a `model.obj` file.
#
# \param scatterer: `dict` Scatterer configuration dictionary (gets modified)
# \param geometry: `dict` Geometry configuration dictionary (gets modified)
# \param max_rad: `float` Maximum allowed radial extension of the aggregate
def write_obj(scatterer, geometry, max_rad):
color_strings = [
"1.0 1.0 1.0\n", # white
"1.0 0.0 0.0\n", # red
"0.0 0.0 1.0\n", # blue
"0.0 1.0 0.0\n", # green
]
color_names = [
"white", "red", "blue", "green"
]
mtl_file = open("model.mtl", "w")
for mi in range(len(color_strings)):
mtl_line = "newmtl mtl{0:d}\n".format(mi)
mtl_file.write(mtl_line)
color_line = color_strings[mi]
mtl_file.write(" Ka " + color_line)
mtl_file.write(" Ks " + color_line)
mtl_file.write(" Kd " + color_line)
mtl_file.write(" Ns 100.0\n")
mtl_file.write(" Tr 0.0\n")
mtl_file.write(" illum 2\n\n")
mtl_file.close()
pl = pv.Plotter()
for si in range(scatterer['nsph']):
sph_type_index = scatterer['vec_types'][si]
color_by_name = color_names[sph_type_index]
radius = scatterer['ros'][sph_type_index - 1] / max_rad
x = geometry['vec_sph_x'][si] / max_rad
y = geometry['vec_sph_y'][si] / max_rad
z = geometry['vec_sph_z'][si] / max_rad
mesh = pv.Sphere(radius, (x, y, z))
mesh.save("tmp_mesh.obj")
pl.add_mesh(mesh) #, color=color_by_name)
mesh_name = "sphere_{0:04d}.obj".format(si)
in_obj_file = open("tmp_mesh.obj", "r")
out_obj_file = open(mesh_name, "w")
in_line = in_obj_file.readline()
out_obj_file.write(in_line)
out_obj_file.write("mtllib model.mtl\n")
out_obj_file.write("usemtl mtl{0:d}\n".format(sph_type_index))
while (in_line != ""):
in_line = in_obj_file.readline()
out_obj_file.write(in_line)
in_obj_file.close()
out_obj_file.close()
pl.export_obj("model.obj")
os.remove("tmp_mesh.obj")
## \brief Exit code (0 for success)
exit_code = main()
exit(exit_code)
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment