I/O — ADCIRC fort.14 and Gmsh

Round-trip read and write for the ADCIRC v55 fort.14 mesh format and the Gmsh 2.2 ASCII .msh format.

ADCIRC fort.14

read_fort14

admesh.read_fort14

read_fort14(path: 'str | os.PathLike[str] | TextIO') -> Mesh

Parse an ADCIRC v55 fort.14 file into a :class:Mesh.

Applies 1-based → 0-based index conversion and depth → elevation sign flip. Unmapped IBTYPE codes are preserved as plain int values in :attr:BoundarySegment.bc_type.

Raises:

Type	Description
Fort14ParseError	On malformed input. The exception carries line_no, expected, and actual.

Source code in src/admesh/fort14.py

def read_fort14(path: "str | os.PathLike[str] | TextIO") -> Mesh:
    """Parse an ADCIRC v55 fort.14 file into a :class:`Mesh`.

    Applies 1-based → 0-based index conversion and depth → elevation sign
    flip. Unmapped IBTYPE codes are preserved as plain ``int`` values in
    :attr:`BoundarySegment.bc_type`.

    Raises
    ------
    Fort14ParseError
        On malformed input. The exception carries ``line_no``,
        ``expected``, and ``actual``.
    """
    lines, handle = _open_text(path)
    try:
        cursor = _Cursor(lines)
        title = cursor.next_line("AGRID title line").strip()

        header = cursor.next_line("'NE NN' counts line").split()
        if len(header) < 2:
            raise cursor.fail("two integers (NE NN)", " ".join(header))
        n_elements = _parse_int(header[0], cursor, "integer NE")
        n_nodes = _parse_int(header[1], cursor, "integer NN")
        if n_elements < 0 or n_nodes < 0:
            raise cursor.fail(
                "non-negative element/node counts", " ".join(header)
            )

        # Node block: NN lines of "id x y depth"
        nodes = np.empty((n_nodes, 2), dtype=np.float64)
        bathymetry = np.empty(n_nodes, dtype=np.float64)
        for i in range(n_nodes):
            line = cursor.next_line(f"node {i + 1} of {n_nodes}")
            tokens = line.split()
            if len(tokens) < 4:
                raise cursor.fail(
                    "node line 'id x y depth' (4 tokens)", line
                )
            node_id = _parse_int(tokens[0], cursor, "integer node id")
            if node_id != i + 1:
                raise cursor.fail(
                    f"monotonic 1-based node id {i + 1}", tokens[0]
                )
            nodes[i, 0] = _parse_float(tokens[1], cursor, "float x")
            nodes[i, 1] = _parse_float(tokens[2], cursor, "float y")
            # ADCIRC stores depth (positive-down); convert to elevation.
            bathymetry[i] = -_parse_float(tokens[3], cursor, "float depth")

        # Element block: NE lines of "id 3 n1 n2 n3" — 1-based.
        elements = np.empty((n_elements, 3), dtype=np.int64)
        for i in range(n_elements):
            line = cursor.next_line(f"element {i + 1} of {n_elements}")
            tokens = line.split()
            if len(tokens) < 5:
                raise cursor.fail(
                    "element line 'id 3 n1 n2 n3' (5 tokens)", line
                )
            elem_id = _parse_int(tokens[0], cursor, "integer element id")
            if elem_id != i + 1:
                raise cursor.fail(
                    f"monotonic 1-based element id {i + 1}", tokens[0]
                )
            nodes_per_elem = _parse_int(
                tokens[1], cursor, "integer node-count = 3"
            )
            if nodes_per_elem != 3:
                raise cursor.fail(
                    "node-count = 3 (triangle)", tokens[1]
                )
            for j in range(3):
                one_based = _parse_int(
                    tokens[2 + j], cursor, f"integer vertex {j + 1}"
                )
                if not (1 <= one_based <= n_nodes):
                    raise cursor.fail(
                        f"vertex id in [1, {n_nodes}]", tokens[2 + j]
                    )
                elements[i, j] = one_based - 1  # → 0-based

        # Open boundary block.
        n_open_segs_line = cursor.next_line("'NOPE' open-segment count")
        n_open_segs = _parse_int(
            n_open_segs_line.split()[0] if n_open_segs_line.strip() else "",
            cursor,
            "integer NOPE",
        )
        # Total open boundary nodes (informational; we don't validate).
        cursor.next_line("'NETA' total open-boundary node count")

        boundaries: list[BoundarySegment] = []
        for _ in range(n_open_segs):
            seg_line = cursor.next_line("open-segment 'NVDLL [IBTYPE]' line")
            seg_tokens = seg_line.split()
            if not seg_tokens:
                raise cursor.fail(
                    "open-segment node-count line", seg_line
                )
            n_seg_nodes = _parse_int(
                seg_tokens[0], cursor, "integer open-segment node count"
            )
            # IBTYPE optional in open-segment header — defaults to 0 (OPEN).
            # Some fixtures annotate the header with an inline "= Number of …"
            # comment; treat the second token as IBTYPE only if it parses as int.
            bc_code = 0
            if len(seg_tokens) >= 2:
                try:
                    bc_code = int(seg_tokens[1])
                except (TypeError, ValueError):
                    bc_code = 0
            ids = np.empty(n_seg_nodes, dtype=np.int64)
            for j in range(n_seg_nodes):
                tok = cursor.next_line(
                    f"open-segment node {j + 1} of {n_seg_nodes}"
                ).split()
                if not tok:
                    raise cursor.fail("open-segment node id", "")
                one_based = _parse_int(
                    tok[0], cursor, "integer node id"
                )
                if not (1 <= one_based <= n_nodes):
                    raise cursor.fail(
                        f"node id in [1, {n_nodes}]", tok[0]
                    )
                ids[j] = one_based - 1
            boundaries.append(
                BoundarySegment(
                    node_ids=ids,
                    bc_type=_coerce_bc(bc_code),
                    is_open=True,
                )
            )

        # Land boundary block.
        n_land_segs_line = cursor.next_line("'NBOU' land-segment count")
        n_land_segs = _parse_int(
            n_land_segs_line.split()[0] if n_land_segs_line.strip() else "",
            cursor,
            "integer NBOU",
        )
        cursor.next_line("'NVEL' total land-boundary node count")

        for _ in range(n_land_segs):
            seg_line = cursor.next_line("land-segment 'NVELL IBTYPE' line")
            seg_tokens = seg_line.split()
            if len(seg_tokens) < 2:
                raise cursor.fail(
                    "land-segment 'NVELL IBTYPE' (2 ints)", seg_line
                )
            n_seg_nodes = _parse_int(
                seg_tokens[0], cursor, "integer land-segment node count"
            )
            bc_code = _parse_int(seg_tokens[1], cursor, "integer IBTYPE")
            ids = np.empty(n_seg_nodes, dtype=np.int64)
            for j in range(n_seg_nodes):
                tok = cursor.next_line(
                    f"land-segment node {j + 1} of {n_seg_nodes}"
                ).split()
                if not tok:
                    raise cursor.fail("land-segment node id", "")
                one_based = _parse_int(
                    tok[0], cursor, "integer node id"
                )
                if not (1 <= one_based <= n_nodes):
                    raise cursor.fail(
                        f"node id in [1, {n_nodes}]", tok[0]
                    )
                ids[j] = one_based - 1
            boundaries.append(
                BoundarySegment(
                    node_ids=ids,
                    bc_type=_coerce_bc(bc_code),
                    is_open=False,
                )
            )

        # Bathymetry: only return when at least one node carries non-zero
        # depth — pure-zero columns are common in synthetic test meshes
        # and we don't want to round-trip a meaningless column.
        bathy_out = bathymetry if np.any(bathymetry != 0.0) else None

        return Mesh(
            nodes=nodes,
            elements=elements,
            boundaries=tuple(boundaries),
            bathymetry=bathy_out,
            quality=None,
            title=title,
        )
    finally:
        if handle is not None:
            handle.close()

write_fort14

admesh.write_fort14

write_fort14(mesh: Mesh, path: 'str | os.PathLike[str] | TextIO', *, precision: int = 6) -> None

Serialize mesh to ADCIRC v55 fort.14 format.

Applies 0-based → 1-based index conversion and elevation → depth sign flip. Coordinates are emitted with precision decimal places.

Parameters:

Name	Type	Description	Default
mesh	Mesh	Triangulation to serialize. Node coordinates, element connectivity, optional bathymetry, and boundary segments are all written.	required
path	str | PathLike[str] | TextIO	Destination file path, or an already-open text stream (written to in place and left open by the caller).	required
precision	int	Number of decimal places for emitted coordinates (default 6). Must be >= 1.	6

Returns:

Type	Description
None	The mesh is written to path for its side effect.

Raises:

Type	Description
ValueError	If precision < 1.

Source code in src/admesh/fort14.py

def write_fort14(
    mesh: Mesh,
    path: "str | os.PathLike[str] | TextIO",
    *,
    precision: int = 6,
) -> None:
    """Serialize ``mesh`` to ADCIRC v55 fort.14 format.

    Applies 0-based → 1-based index conversion and elevation → depth sign
    flip. Coordinates are emitted with ``precision`` decimal places.

    Parameters
    ----------
    mesh : Mesh
        Triangulation to serialize. Node coordinates, element connectivity,
        optional bathymetry, and boundary segments are all written.
    path : str | os.PathLike[str] | TextIO
        Destination file path, or an already-open text stream (written to in
        place and left open by the caller).
    precision : int, optional
        Number of decimal places for emitted coordinates (default ``6``).
        Must be ``>= 1``.

    Returns
    -------
    None
        The mesh is written to ``path`` for its side effect.

    Raises
    ------
    ValueError
        If ``precision < 1``.
    """
    if precision < 1:
        raise ValueError(f"precision must be ≥ 1, got {precision}")

    open_segs = [s for s in mesh.boundaries if s.is_open]
    land_segs = [s for s in mesh.boundaries if not s.is_open]
    bathy = (
        mesh.bathymetry
        if mesh.bathymetry is not None
        else np.zeros(mesh.n_nodes, dtype=np.float64)
    )

    coord_fmt = f"{{:d}} {{:.{precision}f}} {{:.{precision}f}} {{:.{precision}f}}"
    # Open-mode writer
    if hasattr(path, "write"):
        out = path
        close = False
    else:
        out = open(os.fspath(path), "w", encoding="utf-8")
        close = True

    try:
        out.write(f"{mesh.title}\n")
        out.write(f"{mesh.n_elements} {mesh.n_nodes}\n")
        for i in range(mesh.n_nodes):
            out.write(
                coord_fmt.format(
                    i + 1,
                    float(mesh.nodes[i, 0]),
                    float(mesh.nodes[i, 1]),
                    -float(bathy[i]),  # elevation → depth
                )
                + "\n"
            )
        for i in range(mesh.n_elements):
            n0, n1, n2 = mesh.elements[i]
            out.write(
                f"{i + 1} 3 {int(n0) + 1} {int(n1) + 1} {int(n2) + 1}\n"
            )

        n_open_nodes = sum(int(s.node_ids.size) for s in open_segs)
        out.write(f"{len(open_segs)}\n")
        out.write(f"{n_open_nodes}\n")
        for seg in open_segs:
            bc_code = int(seg.bc_type)
            out.write(f"{seg.node_ids.size} {bc_code}\n")
            for nid in seg.node_ids:
                out.write(f"{int(nid) + 1}\n")

        n_land_nodes = sum(int(s.node_ids.size) for s in land_segs)
        out.write(f"{len(land_segs)}\n")
        out.write(f"{n_land_nodes}\n")
        for seg in land_segs:
            bc_code = int(seg.bc_type)
            out.write(f"{seg.node_ids.size} {bc_code}\n")
            for nid in seg.node_ids:
                out.write(f"{int(nid) + 1}\n")
    finally:
        if close:
            out.close()

Fort14ParseError

admesh.Fort14ParseError

Bases: ValueError

Raised by :func:read_fort14 on malformed input.

Attributes:

Name	Type	Description
line_no	int	1-based line number where the error was detected.
expected	str	Short human-readable description of what was expected.
actual	str	The offending line content (truncated to 120 chars).

Source code in src/admesh/fort14.py

class Fort14ParseError(ValueError):
    """Raised by :func:`read_fort14` on malformed input.

    Attributes
    ----------
    line_no : int
        1-based line number where the error was detected.
    expected : str
        Short human-readable description of what was expected.
    actual : str
        The offending line content (truncated to 120 chars).
    """

    def __init__(self, line_no: int, expected: str, actual: str) -> None:
        self.line_no = line_no
        self.expected = expected
        self.actual = (actual or "")[:120]
        super().__init__(
            f"fort.14 parse error at line {line_no}: expected {expected}; "
            f"got {self.actual!r}"
        )

Gmsh 2.2

read_msh

admesh.read_msh

read_msh(path: 'str | os.PathLike[str] | TextIO') -> Mesh

Read a Gmsh ASCII v2.2 .msh file into a :class:Mesh.

Triangles become Mesh.elements (0-based); dim-1 physical groups become ordered :class:BoundarySegment records; node z becomes bathymetry when any value is non-zero.

Parameters:

Name	Type	Description	Default
path	str | PathLike[str] | TextIO	Source .msh file path, or an already-open text stream.	required

Returns:

Type	Description
Mesh	Mesh with 0-based triangle connectivity, ordered boundary segments recovered from dim-1 physical groups, and bathymetry from node z.

Raises:

Type	Description
GmshParseError	If the file is not ASCII Gmsh format 2.x or is otherwise malformed.

Source code in src/admesh/gmsh.py

def read_msh(path: "str | os.PathLike[str] | TextIO") -> Mesh:
    """Read a Gmsh ASCII v2.2 ``.msh`` file into a :class:`Mesh`.

    Triangles become ``Mesh.elements`` (0-based); dim-1 physical groups
    become ordered :class:`BoundarySegment` records; node ``z`` becomes
    ``bathymetry`` when any value is non-zero.

    Parameters
    ----------
    path : str | os.PathLike[str] | TextIO
        Source ``.msh`` file path, or an already-open text stream.

    Returns
    -------
    Mesh
        Mesh with 0-based triangle connectivity, ordered boundary segments
        recovered from dim-1 physical groups, and bathymetry from node ``z``.

    Raises
    ------
    GmshParseError
        If the file is not ASCII Gmsh format 2.x or is otherwise malformed.
    """
    lines, handle = _open_text(path)
    cursor = _Cursor(lines)
    try:
        names: dict[int, str] = {}
        coords: dict[int, tuple[float, float, float]] = {}
        tris: list[tuple[int, int, int]] = []
        # physical tag -> ordered list of (a, b) 0-based line endpoints
        lines_by_tag: dict[int, list[tuple[int, int]]] = {}
        saw_nodes = False

        while True:
            try:
                line = cursor.next_line("$EndElements")
            except GmshParseError:
                break
            line = line.strip()
            if not line:
                continue
            if line == "$MeshFormat":
                hdr = cursor.next_line("version_number file_type data_size").split()
                if not hdr or not hdr[0].startswith("2"):
                    raise cursor.fail("ASCII format 2.x header", " ".join(hdr))
                if len(hdr) >= 2 and hdr[1] != "0":
                    raise cursor.fail("ASCII file_type 0 (binary unsupported)", hdr[1])
                _expect_end(cursor, "$EndMeshFormat")
            elif line == "$PhysicalNames":
                n = _parse_int(cursor.next_line("physical-name count"), cursor)
                for _ in range(n):
                    parts = cursor.next_line("dim tag \"name\"").split(maxsplit=2)
                    if len(parts) < 3:
                        raise cursor.fail("dim tag \"name\"", " ".join(parts))
                    tag = _parse_int(parts[1], cursor)
                    names[tag] = parts[2].strip().strip('"')
                _expect_end(cursor, "$EndPhysicalNames")
            elif line == "$Nodes":
                saw_nodes = True
                n = _parse_int(cursor.next_line("node count"), cursor)
                for _ in range(n):
                    p = cursor.next_line("id x y z").split()
                    if len(p) < 4:
                        raise cursor.fail("id x y z", " ".join(p))
                    nid = _parse_int(p[0], cursor)
                    coords[nid] = (
                        _parse_float(p[1], cursor),
                        _parse_float(p[2], cursor),
                        _parse_float(p[3], cursor),
                    )
                _expect_end(cursor, "$EndNodes")
            elif line == "$Elements":
                n = _parse_int(cursor.next_line("element count"), cursor)
                for _ in range(n):
                    p = cursor.next_line("id type n_tags ...").split()
                    if len(p) < 3:
                        raise cursor.fail("id type n_tags ...", " ".join(p))
                    etype = _parse_int(p[1], cursor)
                    n_tags = _parse_int(p[2], cursor)
                    tags = [_parse_int(t, cursor) for t in p[3 : 3 + n_tags]]
                    conn = [_parse_int(t, cursor) for t in p[3 + n_tags :]]
                    if etype == _GMSH_TRI:
                        if len(conn) != 3:
                            raise cursor.fail("3 triangle node ids", " ".join(p))
                        tris.append((conn[0], conn[1], conn[2]))
                    elif etype == _GMSH_LINE:
                        if len(conn) != 2:
                            raise cursor.fail("2 line node ids", " ".join(p))
                        phys = tags[0] if tags else 0
                        lines_by_tag.setdefault(phys, []).append((conn[0], conn[1]))
                    # Other element types (points, quads, …) ignored for MVP.
                _expect_end(cursor, "$EndElements")
            # Unknown sections are skipped silently (Gmsh forward-compat).

        if not saw_nodes:
            raise cursor.fail("a $Nodes section", "<none>")

        return _assemble_mesh(coords, tris, lines_by_tag, names, cursor)
    finally:
        if handle is not None and hasattr(handle, "close"):
            handle.close()

write_msh

admesh.write_msh

write_msh(mesh: Mesh, path: 'str | os.PathLike[str] | TextIO', *, precision: int = 6) -> None

Serialize mesh to Gmsh ASCII v2.2.

Emits a dim-2 domain physical group for the triangles and one dim-1 group per boundary segment (<label>_<index>). Node z carries bathymetry/elevation (0 when unset).

Parameters:

Name	Type	Description	Default
mesh	Mesh	Triangulation to serialize. Triangles, boundary segments, and optional bathymetry are written.	required
path	str | PathLike[str] | TextIO	Destination file path, or an already-open text stream (written to in place and left open by the caller).	required
precision	int	Number of decimal places for emitted coordinates (default 6). Must be >= 1.	6

Returns:

Type	Description
None	The mesh is written to path for its side effect.

Raises:

Type	Description
ValueError	If precision < 1.

Source code in src/admesh/gmsh.py

def write_msh(
    mesh: Mesh,
    path: "str | os.PathLike[str] | TextIO",
    *,
    precision: int = 6,
) -> None:
    """Serialize ``mesh`` to Gmsh ASCII v2.2.

    Emits a dim-2 ``domain`` physical group for the triangles and one
    dim-1 group per boundary segment (``<label>_<index>``). Node ``z``
    carries bathymetry/elevation (``0`` when unset).

    Parameters
    ----------
    mesh : Mesh
        Triangulation to serialize. Triangles, boundary segments, and optional
        bathymetry are written.
    path : str | os.PathLike[str] | TextIO
        Destination file path, or an already-open text stream (written to in
        place and left open by the caller).
    precision : int, optional
        Number of decimal places for emitted coordinates (default ``6``).
        Must be ``>= 1``.

    Returns
    -------
    None
        The mesh is written to ``path`` for its side effect.

    Raises
    ------
    ValueError
        If ``precision < 1``.
    """
    if precision < 1:
        raise ValueError(f"precision must be ≥ 1, got {precision}")

    bathy = (
        mesh.bathymetry
        if mesh.bathymetry is not None
        else np.zeros(mesh.n_nodes, dtype=np.float64)
    )

    # Allocate physical tags: 1 = domain surface, 2.. = boundary segments.
    phys: list[tuple[int, int, str]] = [(2, 1, "domain")]  # (dim, tag, name)
    seg_tags: list[int] = []
    label_counts: dict[str, int] = {}
    for seg in mesh.boundaries:
        code = int(seg.bc_type)
        label = _BC_TO_LABEL.get(code, f"bc{code}")
        idx = label_counts.get(label, 0)
        label_counts[label] = idx + 1
        tag = len(phys) + 1
        phys.append((1, tag, f"{label}_{idx}"))
        seg_tags.append(tag)

    coord_fmt = f"{{:d}} {{:.{precision}f}} {{:.{precision}f}} {{:.{precision}f}}"

    if hasattr(path, "write"):
        out, close = path, False
    else:
        out, close = open(os.fspath(path), "w", encoding="utf-8"), True

    try:
        out.write("$MeshFormat\n2.2 0 8\n$EndMeshFormat\n")

        out.write("$PhysicalNames\n")
        out.write(f"{len(phys)}\n")
        for dim, tag, name in phys:
            out.write(f'{dim} {tag} "{name}"\n')
        out.write("$EndPhysicalNames\n")

        out.write("$Nodes\n")
        out.write(f"{mesh.n_nodes}\n")
        for i in range(mesh.n_nodes):
            out.write(
                coord_fmt.format(
                    i + 1,
                    float(mesh.nodes[i, 0]),
                    float(mesh.nodes[i, 1]),
                    float(bathy[i]),
                )
                + "\n"
            )
        out.write("$EndNodes\n")

        # Count elements: triangles + boundary line segments.
        n_line_elems = sum(max(int(s.node_ids.size) - 1, 0) for s in mesh.boundaries)
        out.write("$Elements\n")
        out.write(f"{mesh.n_elements + n_line_elems}\n")
        eid = 0
        for i in range(mesh.n_elements):
            eid += 1
            n0, n1, n2 = mesh.elements[i]
            out.write(
                f"{eid} {_GMSH_TRI} 2 1 1 "
                f"{int(n0) + 1} {int(n1) + 1} {int(n2) + 1}\n"
            )
        for seg, tag in zip(mesh.boundaries, seg_tags):
            ids = seg.node_ids
            for k in range(int(ids.size) - 1):
                eid += 1
                out.write(
                    f"{eid} {_GMSH_LINE} 2 {tag} {tag} "
                    f"{int(ids[k]) + 1} {int(ids[k + 1]) + 1}\n"
                )
        out.write("$EndElements\n")
    finally:
        if close:
            out.close()

GmshParseError

admesh.GmshParseError

Bases: ValueError

Raised by :func:read_msh on malformed or unsupported input.

Attributes:

Name	Type	Description
line_no	int	1-based line number where the error was detected.
expected	str	Short human-readable description of what was expected.
actual	str	The offending line content (truncated to 120 chars).

Source code in src/admesh/gmsh.py

class GmshParseError(ValueError):
    """Raised by :func:`read_msh` on malformed or unsupported input.

    Attributes
    ----------
    line_no : int
        1-based line number where the error was detected.
    expected : str
        Short human-readable description of what was expected.
    actual : str
        The offending line content (truncated to 120 chars).
    """

    def __init__(self, line_no: int, expected: str, actual: str) -> None:
        self.line_no = line_no
        self.expected = expected
        self.actual = (actual or "")[:120]
        super().__init__(
            f".msh parse error at line {line_no}: expected {expected}; "
            f"got {self.actual!r}"
        )