#!/usr/bin/env python3 """Write user data to EnSight Gold from Python. This module ships with IST Vector Postprocessor so that any solver, script or notebook can hand it a mesh and some fields without depending on a third-party mesh library. It has **no dependencies** (NumPy is used when present, purely as a fast path) and exposes one function:: from ist_ensight_gold import write_ensight_gold write_ensight_gold( "out/beam.case", nodes=[(0, 0, 0), (1, 0, 0), (1, 1, 0), (0, 1, 0)], elements={"quad4": [(0, 1, 2, 3)]}, node_data={"temperature": [20.0, 35.0, 50.0, 35.0]}, ) That writes ``out/beam.case``, ``out/beam.geo`` and one file per field. Open the ``.case`` file in IST Vector Postprocessor (or EnSight / ParaView). What the writer accepts ----------------------- ``nodes`` ``(nnodes, 3)`` coordinates -- rows, a flat sequence, or a NumPy array. Two-column input is accepted for planar meshes and padded with ``z = 0``. ``elements`` A mapping ``{element_type: connectivity}``, or a sequence of ``(element_type, connectivity)`` pairs when the block order matters. Connectivity is **0-based** by default (see ``node_index_base``) and may be rows or a flat sequence. Element types are the EnSight Gold keywords listed in :data:`ELEMENT_TYPES`. ``node_data`` / ``element_data`` ``{name: values}``. The component count is taken from the shape: =========================== ========================== ================= values written as components =========================== ========================== ================= ``(n,)`` or ``n`` scalars ``scalar per node/element`` 1 ``(n, 3)`` or ``3n`` ``vector per node/element`` 3 ``(n, 6)`` or ``6n`` ``tensor symm ...`` 6 ``(n, 9)`` or ``9n`` ``tensor asym ...`` 9 =========================== ========================== ================= Symmetric tensors are in EnSight's order :data:`SYMMETRIC_TENSOR_ORDER` = ``(11, 22, 33, 12, 13, 23)``; asymmetric tensors are row-major, :data:`ASYMMETRIC_TENSOR_ORDER`. Getting this order wrong mislabels the shear components, so check it against your solver. ``steps`` Transient output: an iterable (a generator is fine -- steps are written as they arrive) of dicts, one per time step:: steps=[{"time": 0.0, "node_data": {"u": u0}}, {"time": 0.5, "node_data": {"u": u1}}, ...] Each step may carry ``time``, ``number`` (the step number that appears in the file names and in the postprocessor's step counter), ``node_data``, ``element_data`` and ``nodes`` (deformed coordinates -- supplying them turns the geometry itself transient). Every step must define the same field names. ``steps`` and ``node_data``/``element_data`` are mutually exclusive. ``parts`` Multi-part output: a list of ``{"name": ..., "elements": {...}}`` dicts instead of a single ``elements`` mapping. Each part is written with just the nodes its elements reference (renumbered locally, as EnSight Gold requires); node fields are sliced to match. Element fields are indexed over the concatenation of every part's blocks, in declaration order. Conventions and limits ---------------------- * Element data is ordered by the concatenation of the connectivity blocks, in the order they are declared -- ``{"tetra4": t, "hexa8": h}`` means all the tetrahedra first, then the hexahedra. * ``binary=True`` writes EnSight "C Binary" (float32, the format EnSight itself stores). ASCII (the default) is written in exponential notation with ``precision`` digits after the point -- 6 by default, which round-trips float32 exactly; ``precision=5`` gives the ``%12.5e`` field EnSight documents. * Non-finite values (NaN, +-Inf) are written as ``0.0``; EnSight has no representation for them. * Field and part names are sanitised: EnSight cannot parse a name containing whitespace. The output base name may not contain whitespace either. Run ``python ist_ensight_gold.py `` to write a small demo data set. """ from __future__ import annotations import array import math import re import struct import sys from pathlib import Path from typing import Any, Iterable, Mapping, Sequence __all__ = [ "write_ensight_gold", "ELEMENT_TYPES", "SYMMETRIC_TENSOR_ORDER", "ASYMMETRIC_TENSOR_ORDER", "EnSightError", ] __version__ = "1.0.0" #: Nodes per element for every EnSight Gold element type this writer emits. #: Mid-side nodes are written verbatim; the postprocessor draws the corner #: topology. ELEMENT_TYPES: "dict[str, int]" = { "point": 1, "bar2": 2, "bar3": 3, "tria3": 3, "tria6": 6, "quad4": 4, "quad8": 8, "tetra4": 4, "tetra10": 10, "pyramid5": 5, "pyramid13": 13, "penta6": 6, "penta15": 15, "hexa8": 8, "hexa20": 20, } #: Component order of a 6-component ``tensor symm`` field. SYMMETRIC_TENSOR_ORDER = ("11", "22", "33", "12", "13", "23") #: Component order of a 9-component ``tensor asym`` field (row-major). ASYMMETRIC_TENSOR_ORDER = ("11", "12", "13", "21", "22", "23", "31", "32", "33") _VAR_KEYWORD = {1: "scalar", 3: "vector", 6: "tensor symm", 9: "tensor asym"} _VAR_SUFFIX = {1: "scl", 3: "vec", 6: "ten", 9: "ten"} _RECORD = 80 # EnSight binary string-record length _F32_MAX = 3.4028234663852886e38 _STEP_KEYS = frozenset(("time", "number", "nodes", "node_data", "element_data")) class EnSightError(ValueError): """Raised when the data handed to :func:`write_ensight_gold` cannot be written.""" # --------------------------------------------------------------------------- # Public entry point # --------------------------------------------------------------------------- def write_ensight_gold( path: "str | Path", nodes: Any, elements: Any = None, *, node_data: "Mapping[str, Any] | None" = None, element_data: "Mapping[str, Any] | None" = None, steps: "Iterable[Mapping[str, Any]] | None" = None, parts: "Sequence[Mapping[str, Any]] | None" = None, description: "str | Sequence[str] | None" = None, part_name: str = "mesh", node_ids: Any = None, element_ids: Any = None, node_index_base: int = 0, binary: bool = False, precision: int = 6, ) -> Path: """Write ``nodes``, ``elements`` and their fields as an EnSight Gold data set. :param path: Output ``.case`` file. Sidecar files are named after its stem and written next to it; parent directories are created. :param nodes: ``(nnodes, 3)`` coordinates (rows, flat sequence or NumPy array). Two columns are accepted and padded with ``z = 0``. :param elements: ``{element_type: connectivity}`` (or a sequence of ``(element_type, connectivity)`` pairs). Omit when using ``parts``. :param node_data: ``{name: values}`` written per node. Steady state only -- use ``steps`` for transient data. :param element_data: ``{name: values}`` written per element, indexed over the concatenation of the connectivity blocks in declaration order. :param steps: Iterable of per-step dicts for transient output; see the module docstring. Mutually exclusive with ``node_data``/``element_data``. :param parts: Multi-part definition, replacing ``elements``. :param description: One or two free-text header lines. :param part_name: Name of the single part, when ``parts`` is not used. :param node_ids: Optional per-node ids (``node id given``). The postprocessor ignores them; EnSight and ParaView display them. :param element_ids: Optional per-element ids, in global element order. :param node_index_base: First node number used by ``elements`` -- ``0`` for Python/NumPy connectivity (the default), ``1`` for Fortran/EnSight. :param binary: Write EnSight "C Binary" (float32) instead of ASCII. :param precision: Digits after the point in ASCII output. ``5`` is the EnSight-documented ``%12.5e``; ``6`` (the default) keeps float32 exact. :returns: The path of the written ``.case`` file. :raises EnSightError: on any inconsistency in the data (bad element type, connectivity out of range, mismatched field length, ...). """ case_path = Path(path) if case_path.suffix.lower() != ".case": case_path = case_path.with_suffix(".case") base = case_path.stem if not base or re.search(r"\s", base): raise EnSightError( "the output base name may not be empty or contain whitespace " "(EnSight cannot parse a file name with a space in it): %r" % base ) out_dir = case_path.parent if str(out_dir): out_dir.mkdir(parents=True, exist_ok=True) if precision < 1 or precision > 17: raise EnSightError("precision must be between 1 and 17, got %r" % (precision,)) float_fmt = "%%%d.%de" % (precision + 7, precision) if steps is not None and (node_data or element_data): raise EnSightError( "pass either steps=... (transient) or node_data=/element_data= " "(steady state), not both" ) if parts is not None and elements is not None: raise EnSightError("pass either elements=... or parts=..., not both") if parts is None and elements is None: raise EnSightError("no elements: pass elements={type: connectivity} or parts=[...]") node_cols = _node_columns(nodes) nnodes = len(node_cols[0]) if nnodes == 0: raise EnSightError("the mesh has no nodes") part_defs = _build_parts(parts, elements, part_name, nnodes, node_index_base) nelements = sum(blk.count for p in part_defs for blk in p.blocks) desc_lines = _description_lines(description) nid = _ids(node_ids, nnodes, "node_ids") eid = _ids(element_ids, nelements, "element_ids") def geo_writer(coords): return lambda fh: _write_geometry(fh, desc_lines, part_defs, coords, nid, eid, binary, float_fmt) # ---- steady state ----------------------------------------------------- if steps is None: geo_name = base + ".geo" _write_file(out_dir / geo_name, geo_writer(node_cols), binary) variables = [] claimed: "dict[str, str]" = {} for location, data in (("node", node_data), ("element", element_data)): for name, values in (data or {}).items(): spec = _field_spec(name, values, location, nnodes, nelements, base, None) _claim(claimed, spec) _write_file( out_dir / spec.file_name, lambda fh, s=spec: _write_variable(fh, s, part_defs, binary, float_fmt), binary, ) variables.append(spec) _write_case(case_path, geo_name, False, variables, [], [], False) return case_path # ---- transient -------------------------------------------------------- times: "list[float]" = [] numbers: "list[int]" = [] used_numbers: "set[int]" = set() claimed: "dict[str, str]" = {} step_specs: "list[_FieldSpec]" = [] field_order: "list[tuple[str, str]] | None" = None moving_geometry = False geo_pattern = base + "_" + "*" * 4 + ".geo" for index, step in enumerate(steps): if not isinstance(step, Mapping): raise EnSightError( "steps[%d] must be a mapping with keys %s, got %r" % (index, sorted(_STEP_KEYS), type(step).__name__) ) unknown = set(step) - _STEP_KEYS if unknown: raise EnSightError( "steps[%d] has unknown key(s) %s; expected any of %s" % (index, sorted(unknown), sorted(_STEP_KEYS)) ) number = int(step.get("number", index)) if number < 0: raise EnSightError("steps[%d]: step number must be >= 0, got %d" % (index, number)) if number in used_numbers: raise EnSightError("steps[%d]: step number %d is already used" % (index, number)) used_numbers.add(number) numbers.append(number) times.append(float(step.get("time", index))) if step.get("nodes") is not None: if index > 0 and not moving_geometry: raise EnSightError( "steps[%d] supplies deformed nodes but earlier steps did not; " "either give every step its own nodes or none" % index ) moving_geometry = True step_cols = _node_columns(step["nodes"]) if len(step_cols[0]) != nnodes: raise EnSightError( "steps[%d]: %d nodes, but the mesh has %d -- a transient geometry " "must keep its node count" % (index, len(step_cols[0]), nnodes) ) elif moving_geometry: raise EnSightError( "steps[%d] has no nodes but an earlier step did; either give every " "step its own nodes or none" % index ) else: step_cols = node_cols if moving_geometry: geo_name = _substitute(geo_pattern, number) _write_file(out_dir / geo_name, geo_writer(step_cols), binary) present = [] for location, key in (("node", "node_data"), ("element", "element_data")): for name, values in (step.get(key) or {}).items(): present.append((location, name)) spec = _field_spec(name, values, location, nnodes, nelements, base, number) if index == 0: _claim(claimed, spec) _write_file( out_dir / spec.file_name, lambda fh, s=spec: _write_variable(fh, s, part_defs, binary, float_fmt), binary, ) if index == 0: step_specs.append(spec) if field_order is None: field_order = present elif sorted(present) != sorted(field_order): missing = sorted(set(field_order) - set(present)) extra = sorted(set(present) - set(field_order)) raise EnSightError( "steps[%d] does not declare the same fields as step 0 " "(missing %s, unexpected %s); EnSight needs one file per field per step" % (index, missing or "none", extra or "none") ) if not numbers: raise EnSightError("steps=... produced no steps") if not moving_geometry: geo_name = base + ".geo" _write_file(out_dir / geo_name, geo_writer(node_cols), binary) model = geo_name else: model = geo_pattern _write_case(case_path, model, moving_geometry, step_specs, times, numbers, True) return case_path # --------------------------------------------------------------------------- # Mesh normalisation # --------------------------------------------------------------------------- class _Block: """One element-type block inside a part.""" __slots__ = ("etype", "conn", "count", "offset") def __init__(self, etype: str, conn: Any, count: int, offset: int): self.etype = etype self.conn = conn # flat, 1-based, local to the part (list or NumPy array) self.count = count self.offset = offset # first global element index of this block class _Part: __slots__ = ("id", "name", "blocks", "node_index", "node_count") def __init__(self, pid: int, name: str, blocks, node_index, node_count: int): self.id = pid self.name = name self.blocks = blocks self.node_index = node_index # None => the part owns every node, in order self.node_count = node_count def _is_seq(x: Any) -> bool: return hasattr(x, "__len__") and not isinstance(x, (str, bytes, bytearray)) def _is_array(x: Any) -> bool: """True for a NumPy-like array (used only as a fast path).""" return hasattr(x, "shape") and hasattr(x, "reshape") and hasattr(x, "dtype") _NUMPY: Any = False def _numpy(): """NumPy if it is installed, else None. Only ever used to go faster.""" global _NUMPY if _NUMPY is False: try: import numpy as np except ImportError: # pragma: no cover - NumPy is optional np = None _NUMPY = np return _NUMPY def _node_columns(nodes: Any) -> "list[Any]": """Normalise coordinates into three column sequences of equal length.""" if nodes is None: raise EnSightError("nodes=None: coordinates are required") if _is_array(nodes): shape = tuple(nodes.shape) if len(shape) == 2 and shape[1] in (2, 3): cols = [nodes[:, j] for j in range(shape[1])] elif len(shape) == 1 and shape[0] % 3 == 0: m = nodes.reshape(shape[0] // 3, 3) cols = [m[:, j] for j in range(3)] else: raise EnSightError( "nodes must have shape (nnodes, 3), (nnodes, 2) or (3*nnodes,), got %r" % (shape,) ) if len(cols) == 2: cols.append([0.0] * len(cols[0])) return cols rows = list(nodes) if rows and _is_seq(rows[0]): width = len(rows[0]) if width not in (2, 3): raise EnSightError("each node needs 2 or 3 coordinates, got %d" % width) xs, ys, zs = [], [], [] for i, row in enumerate(rows): values = list(row) if len(values) != width: raise EnSightError( "node %d has %d coordinates, expected %d" % (i, len(values), width) ) xs.append(float(values[0])) ys.append(float(values[1])) zs.append(float(values[2]) if width == 3 else 0.0) return [xs, ys, zs] flat = [float(v) for v in rows] if len(flat) % 3: raise EnSightError( "a flat coordinate sequence must hold 3 values per node, got %d" % len(flat) ) return [flat[0::3], flat[1::3], flat[2::3]] def _element_blocks(elements: Any, where: str) -> "list[tuple[str, Any]]": if isinstance(elements, Mapping): items = list(elements.items()) elif _is_seq(elements): items = [] for entry in elements: if not (_is_seq(entry) and len(entry) == 2 and isinstance(entry[0], str)): raise EnSightError( "%s: expected {type: connectivity} or [(type, connectivity), ...]" % where ) items.append((entry[0], entry[1])) else: raise EnSightError("%s: expected {type: connectivity} or [(type, connectivity), ...]" % where) if not items: raise EnSightError("%s: no element blocks" % where) return items def _connectivity(etype: str, conn: Any, nnodes: int, base_index: int, where: str) -> Any: """Flatten a connectivity block to validated, flat, 0-based node indices. Returns a NumPy array when NumPy is installed (so a multi-million-element block never becomes a Python list of boxed ints), else a list. """ key = etype.strip().lower() if key not in ELEMENT_TYPES: raise EnSightError( "%s: unknown element type %r; known types are %s" % (where, etype, ", ".join(sorted(ELEMENT_TYPES))) ) npe = ELEMENT_TYPES[key] np = _numpy() def bad(value): return EnSightError( "%s: %r references node %d, outside the mesh's %d nodes " "(node_index_base=%d)" % (where, key, value, nnodes, base_index) ) def wrong_size(total): return EnSightError( "%s: %r needs %d nodes per element but got %d node references " "(not a multiple of %d)" % (where, key, npe, total, npe) ) if np is not None and _is_array(conn): flat = np.asarray(conn).reshape(-1) if flat.size == 0: raise EnSightError("%s: element block %r is empty" % (where, key)) if flat.size % npe: raise wrong_size(int(flat.size)) flat = flat.astype(np.int64, copy=False) if base_index: flat = flat - base_index low, high = int(flat.min()), int(flat.max()) if low < 0: raise bad(low + base_index) if high >= nnodes: raise bad(high + base_index) return flat if _is_array(conn): flat = [int(v) for v in conn.reshape(-1)] else: flat = [] for entry in conn: if _is_seq(entry): flat.extend(int(v) for v in entry) else: flat.append(int(entry)) if not flat: raise EnSightError("%s: element block %r is empty" % (where, key)) if len(flat) % npe: raise wrong_size(len(flat)) out = [] for v in flat: n = v - base_index if n < 0 or n >= nnodes: raise bad(v) out.append(n) return np.asarray(out, dtype=np.int64) if np is not None else out def _build_parts(parts, elements, part_name, nnodes, base_index) -> "list[_Part]": if parts is None: specs = [{"name": part_name, "elements": elements}] share_all_nodes = True else: if not _is_seq(parts) or isinstance(parts, Mapping): raise EnSightError("parts= expects a list of {'name': ..., 'elements': ...} dicts") specs = list(parts) if not specs: raise EnSightError("parts= is empty") share_all_nodes = False known = frozenset(("name", "id", "elements", "node_indices")) result: "list[_Part]" = [] offset = 0 for i, spec in enumerate(specs): if not isinstance(spec, Mapping): raise EnSightError("parts[%d] must be a mapping, got %r" % (i, type(spec).__name__)) unknown = set(spec) - known if unknown: raise EnSightError( "parts[%d] has unknown key(s) %s; expected any of %s" % (i, sorted(unknown), sorted(known)) ) where = "parts[%d]" % i if parts is not None else "elements" name = _clean_text(str(spec.get("name", "part %d" % (i + 1))), 79) or "part %d" % (i + 1) pid = int(spec.get("id", i + 1)) raw_blocks = _element_blocks(spec.get("elements"), where) flat_blocks = [ (etype.strip().lower(), _connectivity(etype, conn, nnodes, base_index, where)) for etype, conn in raw_blocks ] np = _numpy() if share_all_nodes: node_index = None node_count = nnodes lookup = None else: declared = spec.get("node_indices") if declared is not None: node_index = [int(v) - base_index for v in declared] for n in node_index: if n < 0 or n >= nnodes: raise EnSightError( "%s: node_indices references node %d, outside the mesh's %d nodes" % (where, n + base_index, nnodes) ) if np is not None: node_index = np.asarray(node_index, dtype=np.int64) elif np is not None: node_index = np.unique(np.concatenate([c for _, c in flat_blocks])) else: used = set() for _, conn in flat_blocks: used.update(conn) node_index = sorted(used) node_count = len(node_index) if np is not None: # Scatter table: global node -> its 1-based number inside this part. lookup = np.zeros(nnodes, dtype=np.int64) lookup[node_index] = np.arange(1, node_count + 1, dtype=np.int64) else: lookup = {g: k + 1 for k, g in enumerate(node_index)} blocks = [] for etype, conn in flat_blocks: npe = ELEMENT_TYPES[etype] count = len(conn) // npe if lookup is None: local = conn + 1 if _is_array(conn) else [n + 1 for n in conn] elif _is_array(conn): local = lookup[conn] if not bool(local.all()): missing = int(conn[int((local == 0).argmax())]) + base_index raise EnSightError( "%s: element block %r references node %d, which is not in this " "part's node_indices" % (where, etype, missing) ) else: try: local = [lookup[n] for n in conn] except KeyError as exc: raise EnSightError( "%s: element block %r references node %d, which is not in this " "part's node_indices" % (where, etype, int(exc.args[0]) + base_index) ) from None blocks.append(_Block(etype, local, count, offset)) offset += count result.append(_Part(pid, name, blocks, node_index, node_count)) seen = {} for p in result: if p.id in seen: raise EnSightError("two parts share the id %d (%r and %r)" % (p.id, seen[p.id], p.name)) seen[p.id] = p.name return result def _ids(values: Any, count: int, label: str) -> "list[int] | None": if values is None: return None if _is_array(values): out = [int(v) for v in values.reshape(-1)] else: out = [int(v) for v in values] if len(out) != count: raise EnSightError("%s has %d entries, expected %d" % (label, len(out), count)) return out # --------------------------------------------------------------------------- # Field normalisation # --------------------------------------------------------------------------- class _FieldSpec: __slots__ = ("name", "label", "location", "ncomp", "columns", "file_name", "pattern") def __init__(self, name, label, location, ncomp, columns, file_name, pattern): self.name = name self.label = label self.location = location self.ncomp = ncomp self.columns = columns self.file_name = file_name self.pattern = pattern def _ncomp_from_total(total: int, count: int, label: str) -> int: for ncomp in (1, 3, 6, 9): if total == count * ncomp: return ncomp raise EnSightError( "field %r holds %d values for %d items -- expected %d (scalar), %d (vector), " "%d (tensor symm) or %d (tensor asym)" % (label, total, count, count, 3 * count, 6 * count, 9 * count) ) def _field_columns(values: Any, count: int, label: str) -> "tuple[int, list[Any]]": if values is None: raise EnSightError("field %r is None" % label) if _is_array(values): shape = tuple(values.shape) if len(shape) == 2: if shape[0] != count or shape[1] not in (1, 3, 6, 9): raise EnSightError( "field %r has shape %r; expected (%d, c) with c in (1, 3, 6, 9)" % (label, shape, count) ) ncomp = shape[1] matrix = values elif len(shape) == 1: ncomp = _ncomp_from_total(shape[0], count, label) matrix = values.reshape(count, ncomp) else: raise EnSightError("field %r has shape %r; expected 1-D or 2-D" % (label, shape)) return ncomp, [matrix[:, j] for j in range(ncomp)] rows = list(values) if len(rows) == count and rows and _is_seq(rows[0]): ncomp = len(rows[0]) if ncomp not in (1, 3, 6, 9): raise EnSightError( "field %r has %d components per item; expected 1, 3, 6 or 9" % (label, ncomp) ) cols = [[0.0] * count for _ in range(ncomp)] for i, row in enumerate(rows): entry = list(row) if len(entry) != ncomp: raise EnSightError( "field %r: item %d has %d components, expected %d" % (label, i, len(entry), ncomp) ) for j in range(ncomp): cols[j][i] = float(entry[j]) return ncomp, cols flat: "list[float]" = [] for entry in rows: if _is_seq(entry): flat.extend(float(v) for v in entry) else: flat.append(float(entry)) ncomp = _ncomp_from_total(len(flat), count, label) if ncomp == 1: return 1, [flat] return ncomp, [flat[j::ncomp] for j in range(ncomp)] def _field_spec(name, values, location, nnodes, nelements, base, number) -> _FieldSpec: label = _clean_name(name) count = nnodes if location == "node" else nelements if count == 0: raise EnSightError("field %r is per %s but the mesh has none" % (name, location)) ncomp, columns = _field_columns(values, count, label) suffix = _VAR_SUFFIX[ncomp] if location == "node" else "e" + _VAR_SUFFIX[ncomp] stem = "%s.%s" % (base, _clean_file_name(label)) if number is None: pattern = "%s.%s" % (stem, suffix) file_name = pattern else: pattern = "%s_%s.%s" % (stem, "*" * 4, suffix) file_name = _substitute(pattern, number) return _FieldSpec(name, label, location, ncomp, columns, file_name, pattern) def _claim(claimed: "dict[str, str]", spec: _FieldSpec) -> None: """Reserve a field's output file name, so a clash is caught before writing.""" if spec.pattern in claimed: raise EnSightError( "fields %r and %r would both be written to %r -- rename one of them" % (claimed[spec.pattern], spec.name, spec.pattern) ) claimed[spec.pattern] = spec.name def _clean_name(name: Any) -> str: text = _clean_text(str(name), 79) text = re.sub(r"\s+", "_", text) if not text: raise EnSightError("a field name may not be empty") return text def _clean_file_name(label: str) -> str: return re.sub(r"[^A-Za-z0-9_+-]", "_", label) or "field" def _clean_text(text: str, limit: int) -> str: return re.sub(r"[\r\n\t]+", " ", text).strip()[:limit] def _description_lines(description) -> "tuple[str, str]": default = ("EnSight Gold", "written by ist_ensight_gold.py %s" % __version__) if description is None: return default if isinstance(description, str): lines = [description] else: lines = [str(line) for line in description] lines = [_clean_text(line, 79) for line in lines[:2]] while len(lines) < 2: lines.append(default[len(lines)]) return lines[0] or default[0], lines[1] or default[1] def _substitute(pattern: str, number: int) -> str: width = len(re.search(r"\*+", pattern).group(0)) return re.sub(r"\*+", "%0*d" % (width, number), pattern, count=1) # --------------------------------------------------------------------------- # Low-level output # --------------------------------------------------------------------------- def _safe(value: Any) -> float: x = float(value) return x if math.isfinite(x) else 0.0 def _safe32(value: Any) -> float: x = float(value) if not math.isfinite(x): return 0.0 if x > _F32_MAX: return _F32_MAX if x < -_F32_MAX: return -_F32_MAX return x def _record(text: str) -> bytes: raw = text.encode("ascii", "replace")[:_RECORD] return raw + b"\x00" * (_RECORD - len(raw)) def _int32(value: int) -> bytes: return struct.pack(" bytes: if _is_array(column): try: import numpy as np values = np.nan_to_num( np.asarray(column, dtype=np.float64), nan=0.0, posinf=0.0, neginf=0.0 ) np.clip(values, -_F32_MAX, _F32_MAX, out=values) return values.astype(" bytes: if _is_array(values): np = _numpy() if np is not None: return np.asarray(values, dtype=" None: mode = "wb" if binary else "w" kwargs = {} if binary else {"encoding": "ascii", "errors": "replace", "newline": "\n"} with open(path, mode, **kwargs) as fh: body(fh) def _ascii_floats(fh, column, fmt: str) -> None: fh.writelines("%s\n" % (fmt % _safe(v)) for v in column) def _ascii_ints(fh, values) -> None: fh.writelines("%10d\n" % int(v) for v in values) def _write_geometry(fh, desc, part_defs, node_cols, node_ids, element_ids, binary, fmt) -> None: nid_mode = "given" if node_ids is not None else "assign" eid_mode = "given" if element_ids is not None else "assign" x_all, y_all, z_all = node_cols if binary: fh.write(_record("C Binary")) fh.write(_record(desc[0])) fh.write(_record(desc[1])) fh.write(_record("node id " + nid_mode)) fh.write(_record("element id " + eid_mode)) else: fh.write("%s\n%s\n" % desc) fh.write("node id %s\n" % nid_mode) fh.write("element id %s\n" % eid_mode) for part in part_defs: index = part.node_index xs = x_all if index is None else _gather(x_all, index) ys = y_all if index is None else _gather(y_all, index) zs = z_all if index is None else _gather(z_all, index) part_nids = None if node_ids is not None: part_nids = node_ids if index is None else [node_ids[i] for i in index] if binary: fh.write(_record("part")) fh.write(_int32(part.id)) fh.write(_record(part.name)) fh.write(_record("coordinates")) fh.write(_int32(part.node_count)) if part_nids is not None: fh.write(_pack_ints(part_nids)) for column in (xs, ys, zs): fh.write(_pack_floats(column)) else: fh.write("part\n%10d\n%s\n" % (part.id, part.name)) fh.write("coordinates\n%10d\n" % part.node_count) if part_nids is not None: _ascii_ints(fh, part_nids) for column in (xs, ys, zs): _ascii_floats(fh, column, fmt) for block in part.blocks: npe = ELEMENT_TYPES[block.etype] block_eids = ( None if element_ids is None else element_ids[block.offset : block.offset + block.count] ) if binary: fh.write(_record(block.etype)) fh.write(_int32(block.count)) if block_eids is not None: fh.write(_pack_ints(block_eids)) fh.write(_pack_ints(block.conn)) else: fh.write("%s\n%10d\n" % (block.etype, block.count)) if block_eids is not None: _ascii_ints(fh, block_eids) conn = block.conn row_fmt = "%10d" * npe + "\n" if _is_array(conn): rows = conn.reshape(block.count, npe) fh.writelines(row_fmt % tuple(row.tolist()) for row in rows) else: fh.writelines( row_fmt % tuple(conn[start : start + npe]) for start in range(0, len(conn), npe) ) def _write_variable(fh, spec: _FieldSpec, part_defs, binary, fmt) -> None: if binary: fh.write(_record(spec.label)) else: fh.write("%s\n" % spec.label) for part in part_defs: if binary: fh.write(_record("part")) fh.write(_int32(part.id)) else: fh.write("part\n%10d\n" % part.id) if spec.location == "node": index = part.node_index if binary: fh.write(_record("coordinates")) else: fh.write("coordinates\n") for column in spec.columns: values = column if index is None else _gather(column, index) if binary: fh.write(_pack_floats(values)) else: _ascii_floats(fh, values, fmt) else: for block in part.blocks: if binary: fh.write(_record(block.etype)) else: fh.write("%s\n" % block.etype) for column in spec.columns: values = column[block.offset : block.offset + block.count] if binary: fh.write(_pack_floats(values)) else: _ascii_floats(fh, values, fmt) def _gather(column, index): if _is_array(column): return column[index] return [column[i] for i in index] def _write_case(case_path, model, moving_geometry, variables, times, numbers, transient) -> None: lines = ["FORMAT", "type: ensight gold", ""] lines.append("GEOMETRY") if moving_geometry: lines.append("model: 1 %s change_coords_only" % model) else: lines.append("model: %s" % model) if variables: lines.extend(("", "VARIABLE")) for spec in variables: lines.append( "%s per %s: %s%s %s" % ( _VAR_KEYWORD[spec.ncomp], spec.location, "1 " if transient else "", spec.label, spec.pattern, ) ) if transient: lines.extend(("", "TIME")) lines.append("time set: 1") lines.append("number of steps: %d" % len(numbers)) # A constant positive stride is declared as start/increment (what every # EnSight writer emits); anything else needs the explicit number list. increment = numbers[1] - numbers[0] if len(numbers) > 1 else 1 uniform = increment > 0 and all( numbers[i + 1] - numbers[i] == increment for i in range(len(numbers) - 1) ) if uniform: lines.append("filename start number: %d" % numbers[0]) lines.append("filename increment: %d" % increment) else: lines.append("filename numbers:") lines.extend("%d" % n for n in numbers) lines.append("time values:") lines.extend("%.8e" % t for t in times) with open(case_path, "w", encoding="ascii", errors="replace", newline="\n") as fh: fh.write("\n".join(lines)) fh.write("\n") # --------------------------------------------------------------------------- # Demo # --------------------------------------------------------------------------- def _demo(out_dir: Path) -> Path: """Write a small transient two-part data set -- a worked example.""" nodes = [ (0.0, 0.0, 0.0), (1.0, 0.0, 0.0), (1.0, 1.0, 0.0), (0.0, 1.0, 0.0), (0.0, 0.0, 1.0), (1.0, 0.0, 1.0), (1.0, 1.0, 1.0), (0.0, 1.0, 1.0), (2.0, 0.0, 0.0), (2.0, 1.0, 0.0), (2.0, 0.0, 1.0), (2.0, 1.0, 1.0), ] parts = [ {"name": "cube", "id": 1, "elements": {"hexa8": [(0, 1, 2, 3, 4, 5, 6, 7)]}}, {"name": "cap", "id": 2, "elements": {"quad4": [(1, 8, 9, 2), (5, 10, 11, 6)]}}, ] steps = [] for k in range(5): t = 0.25 * k steps.append( { "time": t, "nodes": [(x * (1.0 + 0.05 * k), y, z) for x, y, z in nodes], "node_data": { "temperature": [20.0 + 10.0 * t * (x + y) for x, y, _ in nodes], "displacement": [(0.05 * k * x, 0.0, 0.0) for x, _, _ in nodes], }, "element_data": { # tensor symm: 11, 22, 33, 12, 13, 23 "stress": [(1.0e6 * t, 0.0, 0.0, 0.5e6 * t, 0.0, 0.0)] * 3, "material": [1.0, 2.0, 2.0], }, } ) return write_ensight_gold( out_dir / "demo.case", nodes, parts=parts, steps=steps, description="ist_ensight_gold demo", ) def main(argv: "Sequence[str]") -> int: out_dir = Path(argv[1]) if len(argv) > 1 else Path.cwd() case = _demo(out_dir) print("wrote %s" % case) for sidecar in sorted(case.parent.glob(case.stem + "*")): if sidecar != case: print(" %s" % sidecar.name) return 0 if __name__ == "__main__": raise SystemExit(main(sys.argv))