dothick.py

#!/usr/bin/env python3
"""
Create grid_thick.txt from CA DOTABLE grid output using a simple
SCALEPPM-style neighborhood thickening:

    If a cell is solid, or lies within R=1 of a solid cell,
    write #liqCells = 0.0 in the output file.

Only the #liqCells column is changed.  Composition, phase, temperature,
orientation, and all other columns are preserved.

Default:
    input  = grid.txt
    output = grid_thick.txt
"""

from __future__ import annotations

import argparse
from pathlib import Path
import numpy as np


COL_X = 0
COL_Y = 1
COL_LCELLS = 2
SOLID_THRESHOLD = 0.5
RTHICK = 1


def parse_args() -> argparse.Namespace:
    p = argparse.ArgumentParser(
        description="Write grid_thick.txt using simple R=1 display thickening."
    )
    p.add_argument("input", nargs="?", default="grid.txt",
                   help="Input grid file [default: grid.txt]")
    p.add_argument("-o", "--output", default="grid_thick.txt",
                   help="Output grid file [default: grid_thick.txt]")
    p.add_argument("--threshold", type=float, default=SOLID_THRESHOLD,
                   help="Solid threshold for #liqCells [default: 0.5]")
    return p.parse_args()


def load_grid(path: Path) -> np.ndarray:
    rows = []

    with path.open("r", encoding="utf-8", errors="replace") as f:
        for line in f:
            s = line.strip()
            if not s:
                continue
            if s.upper().startswith("X Y"):
                continue

            parts = s.split()
            if len(parts) >= 8:
                try:
                    rows.append([
                        int(parts[0]), int(parts[1]),
                        float(parts[2]), float(parts[3]),
                        float(parts[4]), float(parts[5]),
                        float(parts[6]), float(parts[7]),
                    ])
                    continue
                except ValueError:
                    pass

            # Legacy fixed-width fallback:
            # (2I4,3F10.3,2F10.2,F10.4)
            try:
                rows.append([
                    int(line[0:4]), int(line[4:8]),
                    float(line[8:18]), float(line[18:28]),
                    float(line[28:38]), float(line[38:48]),
                    float(line[48:58]), float(line[58:68]),
                ])
            except Exception:
                continue

    if not rows:
        raise ValueError(f"No valid grid data found in {path}")

    return np.asarray(rows, dtype=float)


def simple_r1_thicken(data: np.ndarray, threshold: float) -> np.ndarray:
    out = data.copy()

    x = data[:, COL_X].astype(int)
    y = data[:, COL_Y].astype(int)
    lc = data[:, COL_LCELLS]

    # Accept either 0-based or 1-based grid coordinates.
    x0 = int(x.min())
    y0 = int(y.min())
    ix = x - x0
    iy = y - y0

    nx = int(ix.max()) + 1
    ny = int(iy.max()) + 1

    lc_grid = np.full((ny, nx), np.nan, dtype=float)
    lc_grid[iy, ix] = lc

    solid = np.isfinite(lc_grid) & (lc_grid < threshold)
    display_solid = solid.copy()

    ys, xs = np.where(solid)
    for dy in (-1, 0, 1):
        for dx in (-1, 0, 1):
            yy = ys + dy
            xx = xs + dx
            keep = (xx >= 0) & (xx < nx) & (yy >= 0) & (yy < ny)
            display_solid[yy[keep], xx[keep]] = True

    lc_thick_grid = lc_grid.copy()
    lc_thick_grid[display_solid] = 0.0

    out[:, COL_LCELLS] = lc_thick_grid[iy, ix]
    return out


def save_grid(path: Path, data: np.ndarray) -> None:
    np.savetxt(
        path,
        data,
        fmt=[
            "%4d", "%4d", "%10.3f", "%10.3f",
            "%10.3f", "%10.2f", "%10.2f", "%10.4f",
        ],
        header="X Y #liqCells Lcomp Scomp Phase Temp Orient",
        comments="",
    )


def main() -> int:
    args = parse_args()

    in_path = Path(args.input)
    out_path = Path(args.output)

    data = load_grid(in_path)
    out = simple_r1_thicken(data, args.threshold)
    save_grid(out_path, out)

    changed = int(np.sum(data[:, COL_LCELLS] != out[:, COL_LCELLS]))

    print(f"Input   : {in_path}")
    print(f"Output  : {out_path}")
    print(f"Changed : {changed} #liqCells entries")
    print("Method  : SCALEPPM-style R=1 neighborhood thickening")
    return 0


if __name__ == "__main__":
    raise SystemExit(main())