Domain decomposition VTK convert

LAMMPS에서 load balancing을 하게되면, 시뮬레이션 타임병 load balancing된 domain decomposition 정보가 텍스트 파일로 저장된다. 저장되는 정보는 sub domain의 node와 grid정보이다.

이걸 paraview에서 읽어드릴려면, 변환해 줘야 되는데, 아래와 같이 하면 된다.

 

#!/usr/bin/env python3
import sys

input_filename = "./info.balancing"  # or ./info.balancing, etc.

with open(input_filename, 'r') as f:
    lines = f.readlines()

lines = [ln.strip() for ln in lines]

# We'll store data in a dictionary keyed by timestep value:
#   data_dict[timestep] = {
#       "node_lines": [],
#       "cube_lines": [],
#   }
data_dict = {}

current_ts = None
reading_nodes = False
reading_cubes = False

i = 0
n = len(lines)

while i < n:
    line = lines[i]

    if line.startswith("ITEM: TIMESTEP"):
        # The next line holds the numeric timestep
        i += 1
        ts_val = int(lines[i])  # e.g. "10000" or "20000"
        # Ensure we have an entry in data_dict
        if ts_val not in data_dict:
            data_dict[ts_val] = {
                "node_lines": [],
                "cube_lines": []
            }
        current_ts = ts_val
        reading_nodes = False
        reading_cubes = False

    elif line.startswith("ITEM: NODES"):
        # from here on, lines go to node_lines until next ITEM
        reading_nodes = True
        reading_cubes = False

    elif line.startswith("ITEM: CUBES"):
        # from here on, lines go to cube_lines until next ITEM
        reading_nodes = False
        reading_cubes = True

    elif line.startswith("ITEM:"):
        # Some other "ITEM" (like BOX BOUNDS, NUMBER OF NODES, etc.)
        # We can just ignore or handle if needed
        reading_nodes = False
        reading_cubes = False

    else:
        # It's a data line
        if current_ts is not None:
            if reading_nodes and line:
                data_dict[current_ts]["node_lines"].append(line)
            elif reading_cubes and line:
                data_dict[current_ts]["cube_lines"].append(line)

    i += 1


def write_vtk_for_timestep(ts_val, node_lines, cube_lines):
    """Given node lines and cube lines, write domain_tXXXX.vtk."""
    # Parse node lines -> dictionary of coords
    id_to_coords = {}
    for ln in node_lines:
        # Format: ID  Type  X  Y  Z
        toks = ln.split()
        node_id = int(toks[0])
        x, y, z = float(toks[2]), float(toks[3]), float(toks[4])
        id_to_coords[node_id] = (x, y, z)

    sorted_ids = sorted(id_to_coords.keys())
    id_to_new = {}
    coords_list = []
    for i, old_id in enumerate(sorted_ids):
        id_to_new[old_id] = i
        coords_list.append(id_to_coords[old_id])

    num_points = len(coords_list)

    # Parse cube lines -> list of 8 corners each
    cells = []
    for ln in cube_lines:
        toks = ln.split()
        # skip first two (cubeID, cubeType), read n1..n8
        n1, n2, n3, n4, n5, n6, n7, n8 = map(int, toks[2:10])
        c = [
            id_to_new[n1],
            id_to_new[n2],
            id_to_new[n3],
            id_to_new[n4],
            id_to_new[n5],
            id_to_new[n6],
            id_to_new[n7],
            id_to_new[n8],
        ]
        cells.append(c)

    num_cells = len(cells)
    out_name = f"domain_t{ts_val}.vtk"

    with open(out_name, 'w') as out:
        out.write("# vtk DataFile Version 2.0\n")
        out.write("LAMMPS Domain Decomposition Example\n")
        out.write("ASCII\n")
        out.write("DATASET UNSTRUCTURED_GRID\n")

        # Points
        out.write(f"POINTS {num_points} float\n")
        for (x, y, z) in coords_list:
            out.write(f"{x} {y} {z}\n")

        # Cells
        out.write(f"\nCELLS {num_cells} {num_cells*(8+1)}\n")
        for c in cells:
            out.write("8 " + " ".join(map(str, c)) + "\n")

        # Cell types
        out.write(f"\nCELL_TYPES {num_cells}\n")
        for _ in range(num_cells):
            out.write("12\n")  # 12 = VTK_HEXAHEDRON

    print(f"Wrote {out_name}  (Points={num_points}, Cells={num_cells})")


# Finally, process each timestep in ascending order
for ts_val in sorted(data_dict.keys()):
    node_lines = data_dict[ts_val]["node_lines"]
    cube_lines = data_dict[ts_val]["cube_lines"]

    # If either is empty, you probably have partial data
    # but in your final file you do have both for each TS.
    if not node_lines:
        print(f"WARNING: Timestep {ts_val} has no NODES!")
    if not cube_lines:
        print(f"WARNING: Timestep {ts_val} has no CUBES!")

    write_vtk_for_timestep(ts_val, node_lines, cube_lines)

print("Done.")

 

paraview에서는 Surface with Edge가 볼만 하다.