examples: Add utility to export aniso lines

examples/scripts/comma_tools.py

@@ -316,3 +316,51 @@ def prepare_meshio_agglomeration_data(
         # but much faster
         final_data %= modulo_renumbering
     return prepare_meshio_celldata(final_data, cellblocks)
+
+
+def assign_anisotropic_line_data_to_cells(
+    idxs, cells, fc2cc, default_value=-1, modulo_renumbering=None, shuffle=False
+):
+    """Given the CSR definition of the anisotropic lines returned by `CoMMA`, give to
+    each (fine) cell the ID of anistotropic line to which it belongs, otherwise a
+    default values. Since the line definition returned by CoMMA is is terms of coarse
+    cells, knowledge to relate the fine cells to the coarse ones is needed, hence
+    `fc2cc`. It returns a (`numpy`) vector of the same length as `fc2cc` where the value
+    at position `i` is the ID of the anisotropic line to which cell `i` belongs. If cell
+    `i` does not belong to any line, the given default value is used. The line IDs:
+    * are positive and 0-based if the default value is negative
+    * are positive and 1-based if the default value is 0
+    * are negative if the default value is positive
+
+    idxs: array like
+        Indices for lines.
+    cells: array like
+        Flat array containing the all the cells of belonging to the anisotropic lines.
+    fc2cc: array like
+        Result of the agglomeration, ID of the coarse cell to which a fine cell belongs.
+    default_value: int, default = -1
+        Value to assign to cells which do not belong to any line
+    modulo_renumbering: optional, None or int
+        Value of the modulo operation applied to the data.
+    shuffle: bool, default: False
+        Whether to shuffle the ID
+    """
+    n_lines = len(idxs) - 1
+    lines = np.arange(n_lines)
+    f2c = np.asarray(fc2cc)
+    ret = np.full_like(f2c, default_value, dtype=int)
+    if shuffle:
+        np.random.default_rng().shuffle(lines)
+    if modulo_renumbering is not None and modulo_renumbering > 1:
+        lines %= modulo_renumbering
+    if default_value == 0:
+        lines += 1
+    elif default_value > 0:
+        lines = -lines
+    for ln in range(n_lines):
+        # With the following commented definitions, the actual code would me more
+        # readable:
+        # coarse_cells_in_line = cells[idxs[ln] : idxs[ln + 1]]
+        # mask_fine_cells_in_line = np.isin(f2c, coarse_cells_in_line)
+        ret[np.isin(f2c, cells[idxs[ln] : idxs[ln + 1]])] = lines[ln]
+    return ret

examples/scripts/ex_aniso_lines.py

@@ -22,7 +22,11 @@
 from dualGPy.Graph import Graph2D
 from dualGPy.Mesh import Mesh2D
 
-from comma_tools import prepare_meshio_agglomeration_data
+from comma_tools import (
+    assign_anisotropic_line_data_to_cells,
+    prepare_meshio_agglomeration_data,
+    prepare_meshio_celldata,
+)
 
 neigh_type_types = ["Extended", "Pure front advancing"]
 
@@ -182,6 +186,12 @@
     modulo_renumbering=renumber_coarse,
     shuffle=shuffle_coarse,
 )
+out_lines = prepare_meshio_celldata(
+    assign_anisotropic_line_data_to_cells(
+        agglomerationLines_Idx, agglomerationLines, fc_to_cc
+    ),
+    m.mesh.cells,
+)
 
 f2c = np.asarray(fc_to_cc_res, dtype=int)
 # Building a dictionary 'aniso line ID' : [fine cells in line]
@@ -207,6 +217,8 @@
 m.draw_aggl_lines(line_draw, dic, -0.02, 4.02, -0.02, 2.02, show_axes=False)
 
 print(f"Writing in {outname}")
-meshio.Mesh(m.mesh.points, m.mesh.cells, cell_data={"agglomerate": agglo}).write(
-    outname
-)
+meshio.Mesh(
+    m.mesh.points,
+    m.mesh.cells,
+    cell_data={"agglomerate": agglo, "anisotropic_lines": out_lines},
+).write(outname)

examples/scripts/ex_rae.py

@@ -25,6 +25,7 @@
 from dualGPy.Mesh import Mesh2D
 
 from comma_tools import (
+    assign_anisotropic_line_data_to_cells,
     build_coarse_graph,
     compute_neighbourhood_wall_distance,
     prepare_meshio_agglomeration_data,
@@ -59,7 +60,7 @@ def limit_line_length(idxs, cells, max_cells_in_line):
 #################
 input_mesh, input_format = "../meshes/raebis_ansys.msh", "ansys"
 # Number of cells in the aniso lines in the last agglomeration step
-n_cells_in_aniso_line = 2
+# n_cells_in_aniso_line = 2
 ##
 anisotropic = True
 ## CoMMA parameters
@@ -79,7 +80,7 @@ def limit_line_length(idxs, cells, max_cells_in_line):
 # Number of iterations for iterative fine-cell research algorithm
 fc_iter = 1
 # Number of iteration steps to perform
-agglomeration_levels = 1
+agglomeration_levels = 3
 
 # Output-related parameters, they should help with visualization. One can try only one
 # or both at the same time.
@@ -169,7 +170,7 @@ def limit_line_length(idxs, cells, max_cells_in_line):
 # Choosing the whole BL as compliant
 arrayOfFineAnisotropicCompliantCells = squares
 # Weights inversely proportional to the distance so that cells closer to the
-# boundary have higher priority
+# boundary have higher priority (+1 to avoid division by 0)
 weights = np.reciprocal(1.0 + foil_distance.astype(float))
 
 fc_to_cc = np.empty(nb_fc, dtype=CoMMAIndex)
@@ -279,6 +280,17 @@ def limit_line_length(idxs, cells, max_cells_in_line):
         modulo_renumbering=renumber_coarse,
         shuffle=shuffle_coarse,
     )
+    if level == 0:
+        out_data["anisotropic_lines"] = prepare_meshio_celldata(
+            assign_anisotropic_line_data_to_cells(
+                agglomerationLines_Idx,
+                agglomerationLines,
+                fc_to_cc,
+                modulo_renumbering=renumber_coarse,
+                shuffle=shuffle_coarse,
+            ),
+            mesh.mesh.cells,
+        )
     print("OK")
     print()
     if max(fc_to_cc) == 0 and level < agglomeration_levels - 1: