Merge pull request #1303 from nschloe/nastran-fix

Nastran fix for small format

Author: nschloe

setup.cfg

@@ -1,6 +1,6 @@
 [metadata]
 name = meshio
-version = 5.3.3
+version = 5.3.4
 author = Nico Schlömer et al.
 author_email = [email protected]
 description = I/O for many mesh formats

src/meshio/_helpers.py

@@ -176,7 +176,8 @@ def write(filename, mesh: Mesh, file_format: str | None = None, **kwargs):
         if key in num_nodes_per_cell:
             if value.shape[1] != num_nodes_per_cell[key]:
                 raise WriteError(
-                    f"Unexpected cells array shape {value.shape} for {key} cells."
+                    f"Unexpected cells array shape {value.shape} for {key} cells. "
+                    + f"Expected shape [:, {num_nodes_per_cell[key]}]."
                 )
         else:
             # we allow custom keys <https://github.com/nschloe/meshio/issues/501> and

src/meshio/nastran/_nastran.py

@@ -1,6 +1,8 @@
 """
 I/O for Nastran bulk data.
 """
+from __future__ import annotations
+
 import numpy as np
 
 from ..__about__ import __version__
@@ -63,18 +65,22 @@ def read_buffer(f):
     cells = []
     cells_id = []
     cell = None
-    cell_type = None
     point_refs = []
     cell_refs = []
     cell_ref = None
 
-    def add_cell(nastran_type, cell, cell_type, cell_ref):
-        cell_type = nastran_to_meshio_type[keyword]
+    def add_cell(nastran_type, cell, cell_ref):
+        cell_type = nastran_to_meshio_type[nastran_type]
         cell = list(map(int, cell))
 
         # Treat 2nd order CTETRA, CPYRA, CPENTA, CHEXA elements
         if len(cell) > num_nodes_per_cell[cell_type]:
-            assert cell_type in ["tetra", "pyramid", "wedge", "hexahedron"]
+            assert cell_type in [
+                "tetra",
+                "pyramid",
+                "wedge",
+                "hexahedron",
+            ], f"Illegal cell type {cell_type}"
             if cell_type == "tetra":
                 cell_type = "tetra10"
                 nastran_type = "CTETRA_"
@@ -90,6 +96,7 @@ def add_cell(nastran_type, cell, cell_type, cell_ref):
 
         cell = _convert_to_vtk_ordering(cell, nastran_type)
 
+        # decide if we should append cell or start a new cell block
         if len(cells) > 0 and cells[-1][0] == cell_type:
             cells[-1][1].append(cell)
             cells_id[-1].append(cell_id)
@@ -114,48 +121,106 @@ def add_cell(nastran_type, cell, cell_type, cell_ref):
         if next_line.startswith("ENDDATA"):
             break
 
-        # read line and merge with all continuation lines (starting with `+`)
-        chunks = _chunk_line(next_line)
+        # read line and merge with all continuation lines (starting with `+` or
+        # `*` or automatic continuation lines in fixed format)
+        chunks = []
+        c, _ = _chunk_line(next_line)
+        chunks.append(c)
         while True:
             next_line = f.readline()
+
             if not next_line:
                 raise ReadError("Premature EOF")
-            next_line = next_line.rstrip()
+
             # Blank lines or comments
             if len(next_line) < 4 or next_line.startswith(("$", "//", "#")):
                 continue
-            elif next_line[0] == "+":
-                # skip the continuation chunk
-                chunks += _chunk_line(next_line)[1:]
+
+            elif next_line[0] in ["+", "*"]:
+                # From
+                # <https://docs.plm.automation.siemens.com/data_services/resources/nxnastran/10/help/en_US/tdocExt/pdf/User.pdf>:
+                # You can manually specify a continuation by using a
+                # continuation identifier. A continuation identifier is a
+                # special character (+ or *) that indicates that the data
+                # continues on another line.
+                assert len(chunks[-1]) <= 10
+                if len(chunks[-1]) == 10:
+                    # This is a continuation line, so the 10th chunk of the
+                    # previous line must also be a continuation indicator.
+                    # Sometimes its first character is a `+`, but it's not
+                    # always present. Anyway, cut it off.
+                    chunks[-1][-1] = None
+                c, _ = _chunk_line(next_line)
+                c[0] = None
+                chunks.append(c)
+
+            elif len(chunks[-1]) == 10 and chunks[-1][-1] == "        ":
+                # automatic continuation: last chunk of previous line and first
+                # chunk of current line are spaces
+                c, _ = _chunk_line(next_line)
+                if c[0] == "        ":
+                    chunks[-1][9] = None
+                    c[0] = None
+                    chunks.append(c)
+                else:
+                    # not a continuation
+                    break
             else:
                 break
 
+        # merge chunks according to large field format
+        # large field format: 8 + 16 + 16 + 16 + 16 + 8
+        if chunks[0][0].startswith("GRID*"):
+            new_chunks = []
+            for c in chunks:
+                d = [c[0]]
+
+                if len(c) > 1:
+                    d.append(c[1])
+                if len(c) > 2:
+                    d[-1] += c[2]
+
+                if len(c) > 3:
+                    d.append(c[3])
+                if len(c) > 4:
+                    d[-1] += c[4]
+
+                if len(c) > 5:
+                    d.append(c[5])
+                if len(c) > 6:
+                    d[-1] += c[6]
+
+                if len(c) > 7:
+                    d.append(c[7])
+                if len(c) > 8:
+                    d[-1] += c[8]
+
+                if len(c) > 9:
+                    d.append(c[9])
+
+                new_chunks.append(d)
+
+            chunks = new_chunks
+
+        # flatten
+        chunks = [item for sublist in chunks for item in sublist]
+
+        # remove None (continuation blocks)
+        chunks = [chunk for chunk in chunks if chunk is not None]
+
+        # strip chunks
         chunks = [chunk.strip() for chunk in chunks]
 
         keyword = chunks[0]
 
         # Points
-        if keyword == "GRID":
+        if keyword in ["GRID", "GRID*"]:
             point_id = int(chunks[1])
             pref = chunks[2].strip()
             if len(pref) > 0:
                 point_refs.append(int(pref))
             points_id.append(point_id)
             points.append([_nastran_string_to_float(i) for i in chunks[3:6]])
-        elif keyword == "GRID*":  # large field format: 8 + 16*4 + 8
-            point_id = int(chunks[1] + chunks[2])
-            pref = (chunks[3] + chunks[4]).strip()
-            if len(pref) > 0:
-                point_refs.append(int(pref))
-            points_id.append(point_id)
-            chunks2 = _chunk_line(next_line)
-            next_line = f.readline()
-            points.append(
-                [
-                    _nastran_string_to_float(i + j)
-                    for i, j in [chunks[5:7], chunks[7:9], chunks2[1:3]]
-                ]
-            )
 
         # CellBlock
         elif keyword in nastran_to_meshio_type:
@@ -180,7 +245,7 @@ def add_cell(nastran_type, cell, cell_type, cell_ref):
             cell = [item for item in cell if item != ""]
 
             if cell is not None:
-                add_cell(keyword, cell, cell_type, cell_ref)
+                add_cell(keyword, cell, cell_ref)
 
     # Convert to numpy arrays
     points = np.array(points)
@@ -209,10 +274,20 @@ def add_cell(nastran_type, cell, cell_type, cell_ref):
 
 # There are two basic categories of input data formats in NX Nastran:
 #
-#     "Free" format data, in which the data fields are simply separated by commas. This type of data is known as free field data.
-#     "Fixed" format data, in which your data must be aligned in columns of specific width. There are two subcategories of fixed format data that differ based on the size of the fixed column width:
-#         Small field format, in which a single line of data is divided into 10 fields that can contain eight characters each.
-#         Large field format, in which a single line of input is expanded into two lines The first and last fields on each line are eight columns wide, while the intermediate fields are sixteen columns wide. The large field format is useful when you need greater numerical accuracy.
+# - "Free" format data, in which the data fields are simply separated by
+#   commas. This type of data is known as free field data.
+#
+# - "Fixed" format data, in which your data must be aligned in columns of
+#   specific width. There are two subcategories of fixed format data that differ
+#   based on the size of the fixed column width:
+#
+#     - Small field format, in which a single line of data is divided into 10
+#       fields that can contain eight characters each.
+#
+#     - Large field format, in which a single line of input is expanded into
+#       two lines The first and last fields on each line are eight columns wide,
+#       while the intermediate fields are sixteen columns wide. The large field
+#       format is useful when you need greater numerical accuracy.
 #
 # See: https://docs.plm.automation.siemens.com/data_services/resources/nxnastran/10/help/en_US/tdocExt/pdf/User.pdf
 
@@ -262,7 +337,8 @@ def write(filename, mesh, point_format="fixed-large", cell_format="fixed-small")
         for point_id, x in enumerate(points):
             fx = [float_fmt(k) for k in x]
             pref = str(point_refs[point_id]) if point_refs is not None else ""
-            f.write(grid_fmt.format(point_id + 1, pref, fx[0], fx[1], fx[2]))
+            string = grid_fmt.format(point_id + 1, pref, fx[0], fx[1], fx[2])
+            f.write(string)
 
         # CellBlock
         cell_id = 0
@@ -285,6 +361,7 @@ def write(filename, mesh, point_format="fixed-large", cell_format="fixed-small")
                 cell1 = cell + 1
                 cell1 = _convert_to_nastran_ordering(cell1, nastran_type)
                 conn = sjoin.join(int_fmt.format(nid) for nid in cell1[:nipl1])
+
                 if len(cell1) > nipl1:
                     if cell_format == "free":
                         cflag1 = cflag3 = ""
@@ -312,40 +389,62 @@ def _float_rstrip(x, n=8):
 
 def _float_to_nastran_string(value, length=16):
     """
-    Return a value in NASTRAN scientific notation.
+    From
+    <https://docs.plm.automation.siemens.com/data_services/resources/nxnastran/10/help/en_US/tdocExt/pdf/User.pdf>:
+
+    Real numbers, including zero, must contain a decimal point. You can enter
+    real numbers in a variety of formats. For example, the following are all
+    acceptable versions of the real number, seven:
+    ```
+    7.0   .7E1  0.7+1
+    .70+1 7.E+0 70.-1
+    ```
+
+    This methods converts a float value into the corresponding string. Choose
+    the variant with `E` to make the file less ambigious when edited by a
+    human. (`5.-1` looks like 4.0, not 5.0e-1 = 0.5.)
+
     Examples:
-        1234.56789 --> "1.23456789+3"
-        -0.1234 --> "-1.234-1"
-        3.1415926535897932 --> "3.14159265359+0"
+        1234.56789 --> "1.23456789E+3"
+        -0.1234 --> "-1.234E-1"
+        3.1415926535897932 --> "3.14159265359E+0"
     """
-    aux = length - 2
-    # sfmt = "{" + f":{length}s" + "}"
-    sfmt = "{" + ":s" + "}"
-    pv_fmt = "{" + f":{length}.{aux}e" + "}"
+    out = np.format_float_scientific(value, exp_digits=1, precision=11).replace(
+        "e", "E"
+    )
+    assert len(out) <= 16
+    return out
+    # The following is the manual float conversion. Keep it around for a while in case
+    # we still need it.
+
+    # aux = length - 2
+    # # sfmt = "{" + f":{length}s" + "}"
+    # sfmt = "{" + ":s" + "}"
+    # pv_fmt = "{" + f":{length}.{aux}e" + "}"
 
-    if value == 0.0:
-        return sfmt.format("0.")
+    # if value == 0.0:
+    #     return sfmt.format("0.")
 
-    python_value = pv_fmt.format(value)  # -1.e-2
-    svalue, sexponent = python_value.strip().split("e")
-    exponent = int(sexponent)  # removes 0s
+    # python_value = pv_fmt.format(value)  # -1.e-2
+    # svalue, sexponent = python_value.strip().split("e")
+    # exponent = int(sexponent)  # removes 0s
 
-    sign = "-" if abs(value) < 1.0 else "+"
+    # sign = "-" if abs(value) < 1.0 else "+"
 
-    # the exponent will be added later...
-    sexp2 = str(exponent).strip("-+")
-    value2 = float(svalue)
+    # # the exponent will be added later...
+    # sexp2 = str(exponent).strip("-+")
+    # value2 = float(svalue)
 
-    # the plus 1 is for the sign
-    len_sexp = len(sexp2) + 1
-    leftover = length - len_sexp
-    leftover = leftover - 3 if value < 0 else leftover - 2
-    fmt = "{" + f":1.{leftover:d}f" + "}"
+    # # the plus 1 is for the sign
+    # len_sexp = len(sexp2) + 1
+    # leftover = length - len_sexp
+    # leftover = leftover - 3 if value < 0 else leftover - 2
+    # fmt = "{" + f":1.{leftover:d}f" + "}"
 
-    svalue3 = fmt.format(value2)
-    svalue4 = svalue3.strip("0")
-    field = sfmt.format(svalue4 + sign + sexp2)
-    return field
+    # svalue3 = fmt.format(value2)
+    # svalue4 = svalue3.strip("0")
+    # field = sfmt.format(svalue4 + sign + sexp2)
+    # return field
 
 
 def _nastran_string_to_float(string):
@@ -356,14 +455,17 @@ def _nastran_string_to_float(string):
         return float(string[0] + string[1:].replace("+", "e+").replace("-", "e-"))
 
 
-def _chunk_line(line):
-    if "," in line:  # free format
-        chunks = line.split(",")
-    else:  # fixed format
-        CHUNK_SIZE = 8
-        chunks = [line[i : CHUNK_SIZE + i] for i in range(0, 72, CHUNK_SIZE)]
-    # everything after the 9th chunk is ignored
-    return chunks[:9]
+def _chunk_line(line: str) -> tuple[list[str], bool]:
+    # remove terminal newline
+    assert line[-1] == "\n"
+    line = line[:-1]
+    if "," in line:
+        # free format
+        return line.split(","), True
+    # fixed format
+    CHUNK_SIZE = 8
+    chunks = [line[i : CHUNK_SIZE + i] for i in range(0, len(line), CHUNK_SIZE)]
+    return chunks, False
 
 
 def _convert_to_vtk_ordering(cell, nastran_type):