Make sure that no incorrect ori is used on boundary, fix deprecation

kinnala · Aug 9, 2024 · 5d8a0b4 · 5d8a0b4
1 parent 4a0cc8e
commit 5d8a0b4
Show file tree

Hide file tree

Showing 3 changed files with 41 additions and 2 deletions.
diff --git a/skfem/element/discrete_field.py b/skfem/element/discrete_field.py
@@ -55,7 +55,10 @@ def __array_finalize__(self, obj):
     def __getitem__(self, key):
         return np.array(self)[key]
 
-    def __array_wrap__(self, out_arr, context=None):
+    def __array_wrap__(self,
+                       out_arr,
+                       context=None,
+                       return_scalar=False):
         # invalidate attributes after ufuncs
         return np.array(out_arr)
 

diff --git a/skfem/io/meshio.py b/skfem/io/meshio.py
@@ -144,7 +144,7 @@ def from_meshio(m,
                 if not ignore_orientation:
                     try:
                         ori = np.zeros_like(boundaries[k], dtype=np.float64)
-                        t1, _ = mtmp.f2t[:, boundaries[k]]
+                        t1, t2 = mtmp.f2t[:, boundaries[k]]
                         if facets.shape[0] == 2:
                             tangents = (mtmp.p[:, facets[1]]
                                         - mtmp.p[:, facets[0]])
@@ -165,6 +165,10 @@ def from_meshio(m,
                                              - mtmp.p[:, mtmp.t[itr, t1]]),
                                           axis=0)
                         ori = 1 * (ori > 0)
+                        # check that the orientation is not 'illegal'
+                        # this might happen for a hole
+                        # as a consequence, choose the only valid ori
+                        ori[t2 == -1] = 0
                         boundaries[k] = OrientedBoundary(boundaries[k],
                                                          ori)
                     except Exception:

diff --git a/tests/test_mesh.py b/tests/test_mesh.py
@@ -900,3 +900,35 @@ def form(w):
         np.testing.assert_almost_equal(form.assemble(fbasis),
                                        volume,
                                        1e-10)
+
+
+def test_hole_orientation():
+
+    # check that a mesh with a tagged hole is oriented correctly
+
+    m = MeshTri.load(MESH_PATH / 'annulus.msh')
+    mig = MeshTri.load(MESH_PATH / 'annulus.msh',
+                       ignore_orientation=True)
+
+    fbasisint = FacetBasis(m, ElementTriP1(),
+                           facets='inter')
+    fbasisext = FacetBasis(m, ElementTriP1(),
+                           facets='exter')
+
+    fbasisintig = FacetBasis(mig, ElementTriP1(),
+                             facets='inter')
+    fbasisextig = FacetBasis(mig, ElementTriP1(),
+                             facets='exter')
+
+    assert (str(type(m.boundaries['inter']))
+            == "<class 'skfem.generic_utils.OrientedBoundary'>")
+
+    np.testing.assert_almost_equal(
+        fbasisint.normals,
+        fbasisintig.normals,
+    )
+
+    np.testing.assert_almost_equal(
+        fbasisext.normals,
+        fbasisextig.normals,
+    )