Ft mfem 3d single-patch IO

splinepy/io/mfem.py

@@ -3,6 +3,8 @@
 Currently hardcoded for 2D-single-patch-splines.
 """
 
+from ast import literal_eval
+
 import gustaf as _gus
 import numpy as _np
 
@@ -31,7 +33,7 @@
 def load(fname):
     """Reads mfem spline and returns a spline.
 
-    Again, only supports 2D single patch.
+    Supports 2D/3D single patch.
 
     Parameters
     -----------
@@ -86,30 +88,45 @@ def load(fname):
 
     # parse nurbs_dict
     # kvs
-    knot_vectors = nurbs_dict["knotvectors"][1:]
-    knot_vectors[0] = eval("[" + knot_vectors[0].replace(" ", ",") + "]")
-    knot_vectors[1] = eval("[" + knot_vectors[1].replace(" ", ",") + "]")
-    # pop some values
-    # ds
-    degrees = [knot_vectors[0].pop(0), knot_vectors[1].pop(0)]
-    ncps = knot_vectors[0].pop(0) * knot_vectors[1].pop(0)
+    knot_vectors_sec = nurbs_dict["knotvectors"][1:]
+    knot_vectors = []
+    for kvs in knot_vectors_sec:
+        knot_vectors.append(literal_eval("[" + kvs.replace(" ", ",") + "]"))
+
+    # parse degrees and get number of control points for consistency checks
+    degrees = []
+    ncps = 1
+    for kv in knot_vectors:
+        # first value in kv is degrees
+        degrees.append(kv.pop(0))
+        # second value is number of control points
+        ncps *= kv.pop(0)
+
     # ws
     weights = nurbs_dict["weights"]
-    weights = [eval(w) for w in weights]  # hopefully not too slow
+    weights = [literal_eval(w) for w in weights]  # hopefully not too slow
     # cps
     control_points = nurbs_dict["Ordering"]
     control_points = [
-        eval(f"[{cp.replace(' ', ',')}]") for cp in control_points
+        literal_eval(f"[{cp.replace(' ', ',')}]") for cp in control_points
     ]
 
     # double check
     # maybe separate them
     if (
         ncps != len(control_points)
         or ncps != len(weights)
-        or int(nurbs_dict["knotvectors"][0]) != 2
+        or int(nurbs_dict["knotvectors"][0]) != len(degrees)
+        or len(degrees) != len(knot_vectors)
     ):
-        raise ValueError("Inconsistent spline info in " + fname)
+        raise ValueError(
+            f"Inconsistent spline info in {fname}: "
+            f"ncps = {ncps}, "
+            f"len(control_points) = {len(control_points)}, "
+            f"len(weights) = {len(weights)}, "
+            f"len(degrees) = {len(degrees)}, "
+            f"len(knot_vectors) = {len(knot_vectors)}"
+        )
 
     mfem_dict_spline = {
         "degrees": degrees,
@@ -152,7 +169,7 @@ def read_solution(
         collect = False
         while line is not None:
             if collect:  # check this first. should be true most time
-                solution.append(eval(f"[{line.replace(' ', ',')}]"))
+                solution.append(literal_eval(f"[{line.replace(' ', ',')}]"))
             elif hotkey in line:
                 collect = True
 
@@ -579,51 +596,68 @@ def export(fname, nurbs, precision=10):
             "4",
             "",
         )
+    elif nurbs.para_dim == 3:
+        elements_sec = _form_lines("elements", "1", "1 5 0 1 2 3 4 5 6 7", "")
+        boundary_sec = _form_lines(
+            "boundary",
+            "6",
+            "1 3 3 2 1 0",
+            "2 3 4 5 6 7",
+            "3 3 0 1 5 4",
+            "4 3 1 2 6 5",
+            "5 3 2 3 7 6",
+            "6 3 3 0 4 7",
+            "",
+        )
+        edges_sec = _form_lines(
+            "edges",
+            "12",
+            "0 0 1",
+            "0 3 2",
+            "0 4 5",
+            "0 7 6",
+            "1 0 3",
+            "1 1 2",
+            "1 4 7",
+            "1 5 6",
+            "2 0 4",
+            "2 1 5",
+            "2 2 6",
+            "2 3 7",
+            "",
+        )
+        vertices_sec = _form_lines(
+            "vertices",
+            "8",
+            "",
+        )
+    else:
+        raise ValueError("mfem output support 2D and 3D splines.")
 
-        # I am not sure if mixed order is allowed, but in case not, let's
-        # match orders
-        max_degree = max(nurbs.degrees)
-        for i, d in enumerate(nurbs.degrees):
-            d_diff = int(max_degree - d)
-            if d_diff > 0:
-                for _ in range(d_diff):
-                    nurbs.elevate_degrees(i)
-
-        cnr = nurbs.control_mesh_resolutions
-
-        # double-check
-        if not (nurbs.degrees == nurbs.degrees[0]).all():
-            raise RuntimeError(
-                "Something went wrong trying to match degrees of nurbs "
-                + "before export."
-            )
+    cnr = nurbs.control_mesh_resolutions
 
-        # This is reusable
-        def kv_sec(spline):
-            kvs = _form_lines(
-                "knotvectors",
-                str(len(spline.knot_vectors)),
+    # This is reusable
+    def kv_sec(spline):
+        kvs = _form_lines(
+            "knotvectors",
+            str(len(spline.knot_vectors)),
+        )
+        kvs2 = ""
+        for i in range(spline.para_dim):
+            kvs2 += _form_lines(
+                str(spline.degrees[i])
+                + " "
+                + str(cnr[i])
+                + " "
+                + str(list(spline.knot_vectors[i]))[1:-1].replace(",", "")
             )
-            kvs2 = ""
-            for i in range(spline.para_dim):
-                kvs2 += _form_lines(
-                    str(spline.degrees[i])
-                    + " "
-                    + str(cnr[i])
-                    + " "
-                    + str(list(spline.knot_vectors[i]))[1:-1].replace(",", "")
-                )
 
-            kvs = kvs + kvs2
+        kvs = kvs + kvs2
 
-            kvs += "\n"  # missing empty line
+        kvs += "\n"  # missing empty line
+        return kvs
 
-            return kvs
-
-        knotvectors_sec = kv_sec(nurbs)
-
-    else:
-        raise NotImplementedError
+    knotvectors_sec = kv_sec(nurbs)
 
     # disregard inverse
     reorder_ids, _ = dof_mapping(nurbs)
@@ -649,7 +683,7 @@ def kv_sec(spline):
 
     fe_space_sec = _form_lines(
         "FiniteElementSpace",
-        "FiniteElementCollection: NURBS" + str(nurbs.degrees[0]),
+        "FiniteElementCollection: NURBS",
         "VDim: " + str(nurbs.dim),
         "Ordering: 1",
         "",

tests/io/test_mfem_export.py

@@ -129,3 +129,39 @@ def test_mfem_export(to_tmpf, are_stripped_lines_same):
             assert are_stripped_lines_same(
                 base_file.readlines(), tmp_read.readlines(), True
             )
+
+
+def test_mfem_single_patch_export(to_tmpf, are_splines_equal, np_rng):
+    """
+    Test MFEM single patch export routine
+    """
+    # create 2d/3d splines of all supporting types
+    box = splinepy.helpme.create.box(1, 2)
+    boxes = []
+    boxes.extend([box, box.rationalbezier, box.bspline, box.nurbs])
+    boxes.extend([b.create.extruded([0, 0, 3]) for b in boxes[:4]])
+    for spline in boxes:
+        dim = spline.dim
+
+        # make it a bit complex
+        for i in range(dim):
+            if not spline.has_knot_vectors:
+                break
+            spline.insert_knots(i, np_rng.random([5]))
+
+        spline.elevate_degrees(list(range(dim)) * 2)
+
+        for i in range(dim):
+            if not spline.has_knot_vectors:
+                break
+            spline.insert_knots(i, np_rng.random([5]))
+
+        # Test Output
+        with tempfile.TemporaryDirectory() as tmpd:
+            tmpf = to_tmpf(tmpd)
+            splinepy.io.mfem.export(tmpf, spline)
+
+            loaded_spline = splinepy.io.mfem.load(tmpf)
+
+            # all mfem export is nurbs
+            assert are_splines_equal(spline.nurbs, loaded_spline)