[feat] add meshrefine test, use keywords in meshrefine, 0.3.7

Author: fangq

README.md

@@ -4,7 +4,7 @@
 
 * **Copyright**: (C) Qianqian Fang (2024-2025) <q.fang at neu.edu>
 * **License**: GNU Public License V3 or later
-* **Version**: 0.3.6
+* **Version**: 0.3.7
 * **URL**: [https://pypi.org/project/iso2mesh/](https://pypi.org/project/iso2mesh/)
 * **Homepage**: [https://iso2mesh.sf.net](https://iso2mesh.sf.net)
 * **Github**: [https://github.com/NeuroJSON/pyiso2mesh](https://github.com/NeuroJSON/pyiso2mesh)
@@ -277,7 +277,7 @@ tracked in https://github.com/NeuroJSON/pyiso2mesh/issues/1
 | ✅ `smoothsurf.m` | ✅ tested |
 | ✅ `sortmesh.m` | ⭕️ tested |
 | ✅ `mergemesh.m` | ⭕️ tested |
-| ✅ `meshrefine.m` | ⭕️ tested |
+| ✅ `meshrefine.m` | ✅ tested |
 | ✅ `mergesurf.m` | ⭕️ tested |
 | ✅ `surfboolean.m` | ✅ tested |
 | ✅ `fillsurf.m` | ⭕️ tested |

iso2mesh/__init__.py

@@ -157,7 +157,7 @@
     regpt2surf,
 )
 
-__version__ = "0.3.6"
+__version__ = "0.3.7"
 __all__ = [
     "advancefront",
     "barycentricgrid",

iso2mesh/core.py

@@ -1238,7 +1238,7 @@ def remeshsurf(node, face, opt):
     return newno, newfc
 
 
-def meshrefine(node, elem, *args):
+def meshrefine(node, elem, *args, **kwargs):
     """
     meshrefine - refine a tetrahedral mesh by adding new nodes or constraints
 
@@ -1311,13 +1311,11 @@ def meshrefine(node, elem, *args):
             sizefield = np.array(args[1])
         else:
             newpt = np.array(args[1])
-    else:
-        raise ValueError("meshrefine requires at least 3 inputs")
 
-    if isinstance(opt, dict) and "newnode" in opt:
-        newpt = np.array(opt["newnode"])
-    if isinstance(opt, dict) and "sizefield" in opt:
-        sizefield = np.array(opt["sizefield"])
+    if newpt is None and "newnode" in kwargs:
+        newpt = np.array(kwargs["newnode"])
+    if sizefield is None and "sizefield" in kwargs:
+        sizefield = np.array(kwargs["sizefield"])
 
     exesuff = fallbackexeext(getexeext(), "tetgen")
 
@@ -1327,14 +1325,14 @@ def meshrefine(node, elem, *args):
     moreopt = ""
     setquality = False
 
-    if isinstance(opt, dict) and "reratio" in opt:
-        moreopt += f" -q {opt['reratio']:.10f} "
+    if "reratio" in kwargs:
+        moreopt += f" -q {kwargs['reratio']:.10f} "
         setquality = True
-    if isinstance(opt, dict) and "maxvol" in opt:
-        moreopt += f" -a{opt['maxvol']:.10f} "
+    if "maxvol" in kwargs:
+        moreopt += f" -a{kwargs['maxvol']:.10f} "
 
     externalpt = np.empty((0, 3))
-    if isinstance(opt, dict) and "extcmdopt" in opt and newpt is not None:
+    if "extcmdopt" in kwargs and newpt is not None:
         from scipy.spatial import Delaunay
 
         try:
@@ -1439,7 +1437,7 @@ def meshrefine(node, elem, *args):
         holelist = surfseeds(newnode, face[:, :3])
 
         # mesh the extended space
-        ISO2MESH_TETGENOPT = jsonopt("extcmdopt", "-Y", opt)
+        ISO2MESH_TETGENOPT = kwargs.get("extcmdopt", "-Y")
         try:
             if bothsides.shape[0] >= inface.shape[0]:
                 no, el, _ = surf2mesh(
@@ -1514,7 +1512,7 @@ def meshrefine(node, elem, *args):
             map[map == 0] = np.arange(oldsize, allnode.shape[0]) + 1
 
         # merge the external space with the original mesh
-        el2 = map[el[:, :4] - 1] + 1
+        el2 = map[el[:, :4] - 1]
 
         # label all new elements with -1
         if newelem.shape[1] == 5:

iso2mesh/plot.py

@@ -113,7 +113,7 @@ def plotsurf(node, face, *args, **kwargs):
     if "colormap" in locals() and len(colormap) > 0 and not "facecolors" in kwargs:
         kwargs["facecolors"] = colormap
 
-    if "cmap" in kwargs and not "facecolors" in kwargs:
+    if "cmap" in kwargs and not "facecolors" in kwargs and face:
         node_values = node[:, 3] if node.shape[1] > 3 else node[:, 2]
         face_values = np.array([np.mean(node_values[f]) for f in face[:, :3] - 1])
         norm = Normalize(vmin=face_values.min(), vmax=face_values.max())
@@ -337,6 +337,7 @@ def plotmesh(node, *args, **kwargs):
     opt = []
     face = None
     elem = None
+    node = np.array(node)
 
     # Parse inputs: detect selector strings, face/elem arrays, opts
     for i, a in enumerate(args):
@@ -436,11 +437,15 @@ def plotmesh(node, *args, **kwargs):
 
 def _autoscale_3d(ax, points):
     x, y, z = points[:, 0], points[:, 1], points[:, 2]
-    ax.set_xlim([x.min(), x.max()])
-    ax.set_ylim([y.min(), y.max()])
-    ax.set_zlim([z.min(), z.max()])
-    boxas = [x.max() - x.min(), y.max() - y.min(), z.max() - z.min()]
-    ax.set_box_aspect(boxas)
+    boxas = np.array([x.max() - x.min(), y.max() - y.min(), z.max() - z.min()])
+    if boxas[0] > 0:
+        ax.set_xlim([x.min(), x.max()])
+    if boxas[1] > 0:
+        ax.set_ylim([y.min(), y.max()])
+    if boxas[2] > 0:
+        ax.set_zlim([z.min(), z.max()])
+    if np.all(boxas > 0):
+        ax.set_box_aspect(boxas)
 
 
 def _createaxis(*args, **kwargs):

setup.py

@@ -6,7 +6,7 @@
 setup(
     name="iso2mesh",
     packages=["iso2mesh"],
-    version="0.3.6",
+    version="0.3.7",
     license="GPLv3+",
     description="One-liner 3D Surface and Tetrahedral Mesh Generation Toolbox",
     long_description=readme,

test/run_test.py

@@ -941,6 +941,42 @@ def test_extractloops(self):
         bcutloop.pop()
         self.assertEqual(bcutloop, expected_fc)
 
+    def test_meshrefine_addnode(self):
+        no1, el1, fc1 = meshrefine(self.nbox, self.ebox, [[0.2, 1, 1], [0.2, 1, 2]])
+        self.assertEqual(no1.shape[0] - self.nbox.shape[0], 2)
+        self.assertEqual(el1.shape[0] - self.ebox.shape[0], 7)
+        self.assertEqual(no1[-2:, :].tolist(), [[0.2, 1, 1], [0.2, 1, 2]])
+        no1, el1, fc1 = meshrefine(
+            self.nbox, self.ebox, newnode=[[0.2, 1, 1], [0.2, 1, 2]]
+        )
+        self.assertEqual(no1[-2:, :].tolist(), [[0.2, 1, 1], [0.2, 1, 2]])
+
+    def test_meshrefine_maxvol(self):
+        no1, el1, fc1 = meshrefine(self.nbox, self.ebox, maxvol=0.02)
+        self.assertEqual(np.sum(elemvolume(no1[:, :3], el1[:, :4])), 2)
+        self.assertEqual(np.sum(elemvolume(no1[:, :3], fc1[:, :3])), 10)
+        self.assertTrue(np.max(elemvolume(no1[:, :3], el1[:, :4])) < 0.02)
+        no1, el1, fc1 = meshrefine(self.nbox, self.ebox, maxvol=0.01)
+        self.assertTrue(np.max(elemvolume(no1[:, :3], el1[:, :4])) < 0.01)
+
+    def test_meshrefine_sizefield(self):
+        no1, el1, fc1 = meshrefine(
+            self.nbox, self.ebox, sizefield=(self.nbox[:, 0] + 0.1) * 0.3
+        )
+        self.assertEqual(np.sum(elemvolume(no1[:, :3], el1[:, :4])), 2)
+        self.assertAlmostEqual(
+            np.mean(elemvolume(no1[:, :3], el1[:, :4])), 0.00019615535504119262, 8
+        )
+
+    def test_meshrefine_externalnode(self):
+        node, face, elem = meshasphere([0, 0, 0], 24, 5, 100)
+        extnodes = [[-5, -5, 25], [-5, 5, 25], [5, 5, 25], [5, -5, 25]]
+        no1, el1, fc1 = meshrefine(node, elem, newnode=extnodes, extcmdopt="-Y")
+        self.assertEqual(no1.shape[0] - node.shape[0], 4)
+        self.assertAlmostEqual(
+            np.sum(elemvolume(no1[:, :3], el1[:, :4])), 56628.62907381002, 2
+        )
+
 
 class Test_surfboolean(unittest.TestCase):
     def __init__(self, *args, **kwargs):