Use HPC lattice for compact random cluster generation

src/scripts/model_maker.py

@@ -385,17 +385,25 @@ def load_model(model_file):
                     rnd_engine = "COMPACT"
                 if (rnd_engine == "COMPACT"):
                     check = random_compact(sconf, gconf, rnd_seed, max_rad)
-                    if (check == 1):
+                    if (check == -1):
                         print("ERROR: compact random generator works only when all sphere types have the same radius.")
                         return (None, None)
+                    elif (check == -2):
+                        print("ERROR: sub-particle radius larger than particle radius.")
+                        return (None, None)
+                    elif (check == -3):
+                        print("ERROR: requested number of spheres cannot fit in allowed volume.")
+                        return (None, None)
                 elif (rnd_engine == "LOOSE"):
                     # random_aggregate() checks internally whether application is INCLUSION
                     check = random_aggregate(sconf, gconf, rnd_seed, max_rad)
                 else:
                     print("ERROR: unrecognized random generator engine.")
                     return (None, None)
-                if (check != 0):
-                    print("WARNING: %d sphere(s) could not be placed."%check)
+                if (check != sconf['nsph']):
+                    print("WARNING: placed only %d out of %d requested spheres."%(check, sconf['nsph']))
+                    sconf['nsph'] = check
+                    gconf['nsph'] = check
             else:
                 if (len(model['geometry_settings']['x_coords']) != gconf['nsph']):
                     print("ERROR: coordinate vectors do not match the number of spheres!")
@@ -588,15 +596,19 @@ def random_aggregate(scatterer, geometry, seed, max_rad, max_attempts=100):
     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_types = []
     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 (scatterer['application'] == "INCLUSION"):
+            scatterer['vec_types'][0] = 1
     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]
+    vec_types.append(sph_type_index + 1)
     placed_spheres = 1
     attempts = 0
     for i in range(1, nsph):
@@ -670,7 +682,9 @@ def random_aggregate(scatterer, geometry, seed, max_rad, max_attempts=100):
             })
             is_placed = True
             placed_spheres += 1
+            vec_types.append(sph_type_index + 1)
             attempts = 0
+    scatterer['vec_types'] = vec_types
     sph_index = 0
     for sphere in sorted(vec_spheres, key=lambda item: item['itype']):
         scatterer['vec_types'][sph_index] = sphere['itype']
@@ -678,6 +692,7 @@ 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
+    result = placed_spheres
     return result
 
 ## \brief Generate a random compact cluster from YAML configuration options.
@@ -700,45 +715,30 @@ def random_compact(scatterer, geometry, seed, max_rad):
     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
+    radius = scatterer['ros'][0]
+    # Return an error code if types have different radii
     if (max(scatterer['ros']) != min(scatterer['ros'])):
-        result = 1
+        result = -1
+    elif (radius > max_rad):
+        # Requested spheres are larger than the maximum allowed volume.
+        # End function with error code -2.
+        result = -2
     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
-        x_layer_offset = radius
-        y_layer_offset = radius / math.sqrt(3.0)
+        x_centers = np.arange(-1.0 * max_rad + 2.0 * radius, max_rad, 2.0 * radius)
+        x_size = len(x_centers)
+        y_size = int(2.0 * max_rad / ((1.0 + math.sqrt(3.0) / 3.0) * radius))
+        z_size = int(2.0 * max_rad / ((1.0 + 2.0 * math.sqrt(6.0) / 3.0) * radius))
         tmp_spheres = []
-        n_cells = len(x_centers) * len(y_centers) * len(z_centers)
+        n_cells = x_size * y_size * z_size
         print("INFO: the cubic space would contain %d spheres."%n_cells)
-        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
-            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 + x_layer_offset
-                    y = y_centers[yi] + y_offset
-                    z = z_centers[zi]
+        k = 0
+        z = -max_rad + radius
+        while (z < max_rad - radius):
+            j = 0
+            y = -max_rad + radius
+            while (y < max_rad - radius):
+                for i in range(len(x_centers)):
+                    x = (2 * (i + 1) + (j + k) % 2) * radius - max_rad
                     extent = radius + math.sqrt(x * x + y * y + z * z)
                     if (extent < max_rad):
                         tmp_spheres.append({
@@ -747,6 +747,11 @@ def random_compact(scatterer, geometry, seed, max_rad):
                             'y': y,
                             'z': z
                         })
+                #
+                j += 1
+                y = math.sqrt(3.0) * (j + (k % 2) / 3.0) * radius - max_rad + radius
+            k += 1
+            z = 2.0 / 3.0 * math.sqrt(6.0) * k * radius - max_rad + radius
         #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)
@@ -778,20 +783,20 @@ def random_compact(scatterer, geometry, seed, max_rad):
             current_n -= 1
         vec_spheres = []
         sph_index = 0
+        # Generate a vector of types if none is given
+        if (0 in scatterer['vec_types']):
+            tincrement = 1 if scatterer['application'] != "INCLUSION" else 2
+            for ti in range(current_n):
+                itype = tincrement + int(n_types * random.random())
+                scatterer['vec_types'][ti] = itype
+            if (scatterer['application'] == "INCLUSION"):
+                scatterer['vec_types'][0] = 1
         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']
@@ -799,7 +804,7 @@ def random_compact(scatterer, geometry, seed, max_rad):
             geometry['vec_sph_y'][sph_index] = sphere['y']
             geometry['vec_sph_z'][sph_index] = sphere['z']
             sph_index += 1
-    return result
+    return current_n
     
 ## \brief Write the geometry configuration dictionary to legacy format.
 #