Implement compact random generator

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.