Source code for gustaf.io.mixd

"""gustaf/gustaf/io/mixd.py.

io functions for mixd.
"""

import os
import struct

import numpy as np

from gustaf.faces import Faces
from gustaf.io.ioutils import abs_fname, check_and_makedirs
from gustaf.utils import log
from gustaf.volumes import Volumes


[docs] def load( simplex=True, volume=False, fname=None, mxyz=None, mien=None, mrng=None ): """mixd load. To avoid reading minf, all the crucial info can be given as params. Default input will try to import `mxyz`, `mien`, `mrng` from current location and assumes mesh is 2D triangle. Parameters ----------- simplex: bool Default is True. Is it triangle based? volume: bool Default is False. Is it 3D? fname: str Default is None. Specify your mixd file names with ".xns" postfix. Ex) "gustaf.xns" will try to load "gustaf.mxyz", "gustaf.mien", and "gustaf.mrng" mxyz: str Default is None. mien: str Default is None. mrng: str Default is None. This is optional. Returns -------- mesh: Faces or Volumes """ # figure out input type specified_input = mxyz is not None # bare minimum input fname_input = (fname is not None) and not specified_input default_input = not (fname_input or specified_input) if default_input: mxyz = "mxyz" mien = "mien" mrng = "mrng" elif fname_input: fname = abs_fname(fname) fbase, ext = os.path.splitext(fname) if os.path.basename(fbase) == "_": fbase = fbase[:-1] else: fbase += "." mxyz = fbase + "mxyz" mien = fbase + "mien" mrng = fbase + "mrng" # vertices vertices = np.fromfile(mxyz, dtype=">d").astype(np.float64) # connec connec = None try: connec = (np.fromfile(mien, dtype=">i") - int(1)).astype(np.int32) except BaseException: log.debug(f"mien file, `{mien}`, does not exist. Skipping.") # boundary conditions bcs = dict() try: bcs_in = np.fromfile(mrng, dtype=">i").astype(np.int32) # flattened uniq_bcs_in = np.unique(bcs_in) uniq_bcs_in = uniq_bcs_in[uniq_bcs_in > 0] # keep only natural nums sub_elem_ids = np.arange(bcs_in.size) for ubci in uniq_bcs_in: bcs.update({str(ubci): sub_elem_ids[bcs_in == ubci]}) except BaseException: log.debug(f"mrng file, `{mrng}`, does not exist. Skipping.") # reshape vertices vertices = vertices.reshape(-1, 3) if volume else vertices.reshape(-1, 2) # reshape connec if connec is not None: ncol = int(3) if simplex and not volume else int(4) ncol = int(8) if ncol == int(4) and volume and not simplex else ncol connec = connec.reshape(-1, ncol) mesh = Volumes(vertices, connec) if volume else Faces(vertices, connec) # bc if len(bcs) != 0: mesh.BC = bcs return mesh
[docs] def export( mesh, fname, space_time=False, ): """Export in mixd format. Supports triangle, quadrilateral, tetrahedron, and hexahedron semi-discrete and (flat) space-time mesh output. Parameters ----------- mesh: Faces or Volumes fname: str space_time : bool Export Mesh as Space-Time Slab for discontinuous space-time Returns -------- None """ # did you give us an acceptable mesh? acceptable_shapes = ["tri", "quad", "tet", "hexa"] whatami = mesh.whatami if whatami not in acceptable_shapes: raise NotImplementedError( f"Sorry, we can't export {whatami}-shape in mixd format." ) # prepare export location fname = abs_fname(fname) check_and_makedirs(fname) # basic infos dim = mesh.vertices.shape[1] big_endian_int = ">i" big_endian_double = ">d" # prep files fbase, ext = os.path.splitext(fname) if ext.startswith(".xns"): # frequently used case in practice. no base export if os.path.basename(fbase) == "_": fdir = os.path.dirname(fbase) vert_file = os.path.join(fdir, "mxyz") connec_file = os.path.join(fdir, "mien") bc_file = os.path.join(fdir, "mrng") info_file = os.path.join(fdir, "minf") else: fbase += "." vert_file = fbase + "mxyz" connec_file = fbase + "mien" bc_file = fbase + "mrng" info_file = fbase + "minf" else: raise NotImplementedError("`mixd` format only supports xns.") # write v with open(vert_file, "wb") as vf: for v in mesh.vertices.ravel(): vf.write(struct.pack(big_endian_double, v)) if space_time: for v in mesh.vertices.ravel(): vf.write(struct.pack(big_endian_double, v)) # write connec with open(connec_file, "wb") as cf: for c in mesh.elements.ravel() + 1: cf.write(struct.pack(big_endian_int, c)) # write bc with open(bc_file, "wb") as bf: boundaries = make_mrng(mesh) for b in boundaries: bf.write(struct.pack(big_endian_int, b)) # write info with open(info_file, "w") as infof: # if and inf... just can't infof.write(f"# dim: {dim}\n") infof.write(f"# mesh type: {whatami}\n\n") st_factor = 2 if space_time else 1 infof.write(f"nn {int(mesh.vertices.shape[0] * st_factor)}\n") infof.write(f"ne {int(mesh.elements.shape[0])}\n") infof.write(f"nsd {dim}\n") infof.write(f"nen {int(mesh.elements.shape[1] * st_factor)}\n") if space_time: infof.write("space-time on\n\n\n") else: infof.write("semi-discrete on\n\n\n") # BC guide infof.write("# boundary name : referenced number.\n") for i, bname in enumerate(mesh.BC.keys()): infof.write(f"# {bname} : {i + 1}\n") # signature infof.write("\n\n\n# MIXD generated using `gustaf`.\n")
[docs] def make_mrng(mesh): """ Builds and return mrng array based on `mesh.BC` Supports `Faces` and `Volumes`. Parameters ----------- mesh: Faces or Volumes Number of participating elements Returns -------- boundaries : ndarray The mrng-array. """ # determine number of sub elements whatami = mesh.whatami nbelem = 3 if whatami.startswith("quad") or whatami.startswith("tet"): nbelem += 1 elif whatami.startswith("hexa"): nbelem += 3 # init boundaries with -1, as it is the value for non-boundary. # alternatively, they could be (-1 * neighbor_elem_id). # But they aren't. boundaries = np.empty(mesh.elements.shape[0] * nbelem, dtype=int) boundaries[:] = -1 for i, belem_ids in enumerate(mesh.BC.values()): boundaries[belem_ids] = i + 1 # bid starts at 1 return boundaries