Implement random generator and include it in model_maker (3f7c55aa) · Commits · Giacomo Mulas / NP_TMcode

src/scripts/generator.py

+23 −13

Original line number	Diff line number	Diff line
		@@ -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']))

src/scripts/model_maker.py

+137 −13

Original line number	Diff line number	Diff line
		@@ -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):
		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'] = [int(str_typ) for str_typ in sph_types]
		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.