Implement compact and loose random engines and OBJ exporting capabilities (33ce0d93) · Commits · Giacomo Mulas / NP_TMcode

src/scripts/model_maker.py

+66 −22

Original line number	Diff line number	Diff line
		@@ -24,6 +24,7 @@
		# The script requires python3.

		import math
		import multiprocessing
		import numpy as np
		import os
		import pdb
		@@ -369,10 +370,21 @@ def load_model(model_file):
		max_rad = float(model['particle_settings']['max_rad'])
		# random_aggregate() checks internally whether application is INCLUSION
		#random_aggregate(sconf, gconf, rnd_seed, max_rad)
		rnd_engine = "COMPACT"
		try:
		rnd_engine = model['system_settings']['rnd_engine']
		except KeyError:
		# use compact generator, if no specification is given
		rnd_engine = "COMPACT"
		if (rnd_engine == "COMPACT"):
		check = random_compact(sconf, gconf, rnd_seed, max_rad)
		elif (rnd_engine == "LOOSE"):
		check = random_aggregate(sconf, gconf, rnd_seed, max_rad)
		else:
		print("ERROR: unrecognized random generator engine.")
		return (None, None)
		if (check != 0):
		print("INFO: stopping with exit code %d."%check)
		exit(check)
		print("WARNING: %d sphere(s) could not be placed."%check)
		else:
		if (len(model['geometry_settings']['x_coords']) != gconf['nsph']):
		print("ERROR: coordinate vectors do not match the number of spheres!")
		@@ -386,6 +398,36 @@ def load_model(model_file):
		if (max_rad == 0.0):
		max_rad = 20.0 * max(sconf['ros'])
		write_obj(sconf, gconf, max_rad)
		try:
		max_gpu_ram = int(model['system_settings']['max_gpu_ram'])
		if (max_gpu_ram > 0):
		max_gpu_ram_bytes = max_gpu_ram * 1024 * 1024 * 1024
		matrix_dim = 2 * gconf['nsph'] * gconf['li'] * (gconf['li'] + 2)
		matrix_size_bytes = 16 * matrix_dim * matrix_dim
		if (matrix_size_bytes < max_gpu_ram_bytes):
		max_gpu_processes = int(max_gpu_ram_bytes / matrix_size_bytes)
		print("INFO: system supports up to %d simultaneous processes on GPU."%max_gpu_processes)
		else:
		print("WARNING: estimated matrix size is larger than available GPU memory!")
		else:
		print("INFO: no GPU RAM declared.")
		max_host_ram = int(model['system_settings']['max_host_ram'])
		if (max_host_ram > 0):
		max_host_ram_bytes = max_host_ram * 1024 * 1024 * 1024
		matrix_dim = 2 * gconf['nsph'] * gconf['li'] * (gconf['li'] + 2)
		matrix_size_bytes = 16 * matrix_dim * matrix_dim
		if (matrix_size_bytes < max_host_ram_bytes):
		max_host_processes = int(max_host_ram_bytes / matrix_size_bytes / 2)
		print("INFO: system supports up to %d simultaneous processes."%max_host_processes)
		else:
		print("WARNING: estimated matrix size is larger than available host memory!")
		else:
		print("WARNING: no host RAM declared!")
		except KeyError as ex:
		print(ex)
		print("WARNING: missing system description! Cannot estimate recommended execution.")
		cpu_count = multiprocessing.cpu_count()
		print("INFO: the number of detected CPUs is %d."%cpu_count)
		else: # model is None
		print("ERROR: could not parse " + model_file + "!")
		return (sconf, gconf)
		@@ -554,7 +596,6 @@ def random_aggregate(scatterer, geometry, seed, max_rad, max_attempts=100):
		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()
		vec_phis[i] = 2.0 * math.pi * random.random()
		@@ -662,16 +703,21 @@ def random_compact(scatterer, geometry, seed, max_rad):
		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
		y_offset = radius
		x_layer_offset = radius
		y_layer_offset = radius / math.sqrt(3.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)
		n_max_spheres = int((max_rad / radius) * (max_rad / radius) * (max_rad / radius) * 0.74)
		print("INFO: the maximum radius allows for %d spheres."%n_max_spheres)
		for zi in range(len(z_centers)):
		if (x_layer_offset == 0.0):
		x_layer_offset = radius
		else:
		x_layer_offset = 0.0
		if (y_offset == 0.0):
		y_offset = radius / math.sqrt(3.0)
		y_offset = radius
		else:
		y_offset = 0.0
		for yi in range(len(y_centers)):
		@@ -680,7 +726,7 @@ def random_compact(scatterer, geometry, seed, max_rad):
		else:
		x_offset = 0.0
		for xi in range(len(x_centers)):
		x = x_centers[xi] + x_offset
		x = x_centers[xi] + x_offset + x_layer_offset
		y = y_centers[yi] + y_offset
		z = z_centers[zi]
		extent = radius + math.sqrt(x * x + y * y + z * z)
		@@ -694,10 +740,11 @@ def random_compact(scatterer, geometry, seed, max_rad):
		#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
		rho = 10.0 * max_rad
		discard_rad = 100.0 * max_rad
		while (current_n > nsph):
		theta = 2.0 * math.pi * random.random()
		phi = math.pi * random.random()
		theta = math.pi * random.random()
		phi = 2.0 * 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)
		@@ -705,7 +752,7 @@ def random_compact(scatterer, geometry, seed, max_rad):
		minimum_distance = 1000.0 * max_rad
		for di in range(len(tmp_spheres)):
		x1 = tmp_spheres[di]['x']
		if (x1 == max_rad):
		if (x1 == discard_rad):
		continue
		y1 = tmp_spheres[di]['y']
		z1 = tmp_spheres[di]['z']
		@@ -717,7 +764,7 @@ def random_compact(scatterer, geometry, seed, max_rad):
		if (distance < minimum_distance):
		closest_index = di
		minimum_distance = distance
		tmp_spheres[closest_index]['x'] = max_rad
		tmp_spheres[closest_index]['x'] = discard_rad
		current_n -= 1
		vec_spheres = []
		sph_index = 0
		@@ -900,17 +947,13 @@ def write_legacy_sconf(conf):
		# 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):
		out_dir = scatterer['out_file'].absolute().parent
		out_model_path = Path(out_dir, "model.obj")
		out_material_path = Path(out_dir, "model.mtl")
		out_model_path = Path(str(geometry['out_file']) + ".obj")
		out_material_path = Path(str(geometry['out_file']) + ".mtl")
		color_strings = [
		"1.0 1.0 1.0\n", # white
		"1.0 0.0 0.0\n", # red
		@@ -945,14 +988,15 @@ def write_obj(scatterer, geometry, max_rad):
		pl.add_mesh(mesh, color=None)
		pl.export_obj(str(Path(str(out_dir), "TMP_MODEL.obj")))
		tmp_model_file = open(str(Path(str(out_dir), "TMP_MODEL.obj")), "r")
		out_model_file = open(str(Path(str(out_dir), "model.obj")), "w")
		out_model_file = open(str(out_model_path), "w")
		mtl_line = "mtllib {0:s}\n".format(out_material_path.name)
		sph_index = 0
		sph_type_index = 0
		old_sph_type_index = 0
		str_line = tmp_model_file.readline()
		while (str_line != ""):
		if (str_line.startswith("mtllib")):
		str_line = "mtllib model.mtl\n"
		str_line = mtl_line
		elif (str_line.startswith("g ")):
		sph_index += 1
		sph_type_index = scatterer['vec_types'][sph_index - 1]