Repository: kinnala/scikit-fem

Issue:
FacetBasis outward normal vectors are false
Hello
I am new user, and I began with 8.1 version, I encountered the following issue.
FacetBasis didn't work well since upgrading to 10.0.0-3, the outward normal vectors at boundary are sometimes false.
Here an example in order to show this issue, using an annular domain, there are 2 boundaries 'Int' and 'Ext'.
In G_Vh  normal vectors are OK on two boundary 'Int' and 'Ext'
In Gm_Vh normal vectors are OK on  boundary 'Ext'
In Gu_Vh normal vectors are not OK on  boundary 'Int', this could be shown on the figure only two vectors are given on the circular boundary (red points) and in  Gu_Vh.normals.value there are only two normal vectors.
array([[[ 0.89926879,  0.89926879],
        [-0.89926879, -0.89926879],
        ........
        [ 0.89926879,  0.89926879],
        [-0.89926879, -0.89926879]],

       [[-0.43739644, -0.43739644],
        [ 0.43739644,  0.43739644],
        .......
        [-0.43739644, -0.43739644],
        [ 0.43739644,  0.43739644]]])

Here attached the mesh file created with gmsh
[Anneau01.zip](https://github.com/user-attachments/files/16552877/Anneau01.zip)

import os
print(os.environ['PATH'])
from pathlib import Path
import numpy as np
from skfem import *
from skfem.helpers import dot, grad
from skfem.io import from_meshio
import matplotlib.pyplot as plt

mesh = Mesh.load('Anneau01.msh')
element = ElementTriP1()
Vh = Basis(mesh, element)
G_Vh = FacetBasis(mesh, element)
Gm_Vh = FacetBasis(mesh, element, facets=mesh.boundaries['Ext'])
Gu_Vh = FacetBasis(mesh, element, facets=mesh.boundaries['Int'])

plt.figure(1)
plt.plot(G_Vh.normals.value[0],G_Vh.normals.value[1],'*',color='blue')
plt.plot(Gm_Vh.normals.value[0],Gm_Vh.normals.value[1],'*',color='green')
plt.plot(Gu_Vh.normals.value[0],Gu_Vh.normals.value[1],'*',color='red')
plt.show()


Discussion:

Patch:
diff --git a/README.md b/README.md
index 9794ba34..8332a2b7 100644
--- a/README.md
+++ b/README.md
@@ -209,6 +209,10 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html)
 with respect to documented and/or tested features.
 
+### Unreleased
+
+- Fixed: `Mesh.load` returned incorrect orientation for tagged holes
+
 ### [10.0.1] - 2024-08-06
 
 - Fixed: `Mesh.load` returned incorrect orientation for some Gmsh
diff --git a/docs/examples/meshes/annulus.msh b/docs/examples/meshes/annulus.msh
new file mode 100644
index 00000000..93f5155e
--- /dev/null
+++ b/docs/examples/meshes/annulus.msh
@@ -0,0 +1,271 @@
+$MeshFormat
+4.1 0 8
+$EndMeshFormat
+$PhysicalNames
+3
+1 7 "exter"
+1 8 "inter"
+2 9 "all"
+$EndPhysicalNames
+$Entities
+2 2 1 0
+2 0.1 0 0 0 
+3 0.5 0 0 0 
+2 -0.1000001 -0.1000001 -1.000000000310317e-07 0.1000001 0.1000001 1.000000000310317e-07 1 8 2 2 -2 
+3 -0.5000000999999999 -0.5000000999999999 -1e-07 0.5000000999999999 0.5000000999999999 1e-07 1 7 2 3 -3 
+1 -0.7757299569548867 -0.7765209529517499 -1e-07 0.7773807226745383 0.7765209529517506 1e-07 1 9 2 3 -2 
+$EndEntities
+$Nodes
+5 60 1 60
+0 2 0 1
+1
+0.1 0 0
+0 3 0 1
+2
+0.5 0 0
+1 2 0 6
+3
+4
+5
+6
+7
+8
+0.06234898018587345 0.07818314824680291 0
+-0.02225209339563107 0.09749279121818245 0
+-0.09009688679024162 0.04338837391175643 0
+-0.09009688679024234 -0.04338837391175492 0
+-0.02225209339563258 -0.09749279121818211 0
+0.06234898018587293 -0.07818314824680334 0
+1 3 0 14
+9
+10
+11
+12
+13
+14
+15
+16
+17
+18
+19
+20
+21
+22
+0.4567727288213005 0.2033683215379 0
+0.3345653031794293 0.371572412738697 0
+0.1545084971874739 0.4755282581475767 0
+-0.05226423163382667 0.4972609476841367 0
+-0.2499999999999997 0.4330127018922195 0
+-0.4045084971874736 0.2938926261462368 0
+-0.4890738003669028 0.1039558454088799 0
+-0.4890738003669028 -0.1039558454088798 0
+-0.4045084971874737 -0.2938926261462365 0
+-0.2500000000000006 -0.433012701892219 0
+-0.05226423163382756 -0.4972609476841366 0
+0.1545084971874732 -0.4755282581475769 0
+0.3345653031794286 -0.3715724127386976 0
+0.4567727288213003 -0.2033683215379005 0
+2 1 0 38
+23
+24
+25
+26
+27
+28
+29
+30
+31
+32
+33
+34
+35
+36
+37
+38
+39
+40
+41
+42
+43
+44
+45
+46
+47
+48
+49
+50
+51
+52
+53
+54
+55
+56
+57
+58
+59
+60
+-0.1662660720751967 -0.0001311801848818683 0
+-0.1045567449378255 -0.1308075055877821 0
+-0.1024111798352759 0.1136488548414939 0
+0.06431348789077392 -0.1636876743640242 0
+0.02487665115318162 0.1500884344384296 0
+0.1629248213409642 0.07753041812638373 0
+0.1488830826225129 -0.07169831376082697 0
+-0.3209390580020073 -0.1451136089557044 0
+-0.3145402434788418 0.1421080805038236 0
+-0.1032529705984628 -0.3445272085041622 0
+-0.1059571516202201 0.3140208936830532 0
+0.1646819614582139 0.2889900680375265 0
+0.1581086517838572 -0.292595328446286 0
+0.3026764156309572 0.07199706676686718 0
+-0.3520064066612191 -0.0006320669365808482 0
+-0.2379744838939296 -0.2716090417956659 0
+-0.2517379717048167 0.2793380604463241 0
+0.03921590899191962 -0.3742254697096736 0
+0.03921590899192019 0.3742254697096739 0
+0.3045000554910642 0.2213068982403297 0
+0.3045362052741258 -0.2212585044490235 0
+0.3576425992375117 -0.08308252471836414 0
+0.2588843611999955 -0.02461302275662827 0
+0.20994004363952 -0.1529119138267931 0
+0.2090337779492169 0.1596474668918353 0
+0.1770227538768362 1.068388044943303e-07 0
+0.1061963146536865 0.169702022396523 0
+-0.03344465612838528 -0.1750662965616162 0
+0.001153272113847848 -0.2616100911278346 0
+-0.06554426470640014 0.1861664489380913 0
+0.02724722694203124 0.2471975697668277 0
+-0.1617759847287665 -0.07832172546951767 0
+-0.1984773500942226 -0.1680957205934023 0
+-0.2508378444878755 -0.04860037888487612 0
+-0.1607132009296628 0.07587910713981455 0
+-0.1844062125741907 0.1687752223551687 0
+-0.1127591563904036 -0.2252850401436459 0
+-0.2439645438503664 0.04451873273149497 0
+$EndNodes
+$Elements
+3 120 1 120
+1 2 1 7
+1 1 3 
+2 3 4 
+3 4 5 
+4 5 6 
+5 6 7 
+6 7 8 
+7 8 1 
+1 3 1 15
+8 2 9 
+9 9 10 
+10 10 11 
+11 11 12 
+12 12 13 
+13 13 14 
+14 14 15 
+15 15 16 
+16 16 17 
+17 17 18 
+18 18 19 
+19 19 20 
+20 20 21 
+21 21 22 
+22 22 2 
+2 1 2 98
+23 28 48 36 
+24 26 46 29 
+25 36 47 28 
+26 35 46 26 
+27 26 51 35 
+28 33 58 52 
+29 34 53 49 
+30 36 48 45 
+31 45 46 44 
+32 29 46 45 
+33 17 30 16 
+34 15 31 14 
+35 34 49 47 
+36 19 32 18 
+37 13 33 12 
+38 44 46 43 
+39 11 34 10 
+40 21 35 20 
+41 9 36 2 
+42 52 53 33 
+43 42 47 36 
+44 39 58 33 
+45 17 38 30 
+46 31 39 14 
+47 30 37 16 
+48 15 37 31 
+49 1 28 3 
+50 8 29 1 
+51 3 27 4 
+52 7 26 8 
+53 6 24 7 
+54 4 25 5 
+55 5 23 6 
+56 32 38 18 
+57 13 39 33 
+58 19 40 32 
+59 33 41 12 
+60 54 56 55 
+61 22 43 21 
+62 21 43 35 
+63 2 44 22 
+64 10 42 9 
+65 20 40 19 
+66 12 41 11 
+67 18 38 17 
+68 14 39 13 
+69 16 37 15 
+70 11 41 34 
+71 35 40 20 
+72 34 42 10 
+73 36 44 2 
+74 58 60 57 
+75 9 42 36 
+76 50 59 51 
+77 35 51 40 
+78 41 53 34 
+79 32 59 38 
+80 37 60 31 
+81 45 48 29 
+82 33 53 41 
+83 40 51 32 
+84 36 45 44 
+85 51 59 32 
+86 30 56 37 
+87 43 46 35 
+88 34 47 42 
+89 55 56 30 
+90 38 55 30 
+91 31 60 58 
+92 31 58 39 
+93 47 49 28 
+94 38 59 55 
+95 25 58 57 
+96 22 44 43 
+97 56 60 37 
+98 54 55 24 
+99 50 51 26 
+100 27 53 52 
+101 7 50 26 
+102 27 52 4 
+103 4 52 25 
+104 24 50 7 
+105 29 48 1 
+106 1 48 28 
+107 28 49 3 
+108 6 54 24 
+109 25 57 5 
+110 3 49 27 
+111 5 57 23 
+112 23 54 6 
+113 23 56 54 
+114 55 59 24 
+115 23 60 56 
+116 24 59 50 
+117 57 60 23 
+118 26 29 8 
+119 52 58 25 
+120 49 53 27 
+$EndElements
diff --git a/skfem/element/discrete_field.py b/skfem/element/discrete_field.py
index d188d195..ab548369 100644
--- a/skfem/element/discrete_field.py
+++ b/skfem/element/discrete_field.py
@@ -55,7 +55,10 @@ class DiscreteField(ndarray):
     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
index d16d1b53..08e65c82 100644
--- 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:


Test patch:
diff --git a/tests/test_mesh.py b/tests/test_mesh.py
index abf42c65..35c26848 100644
--- a/tests/test_mesh.py
+++ b/tests/test_mesh.py
@@ -900,3 +900,35 @@ def test_load_orientation(mesh, volume, etype):
         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,
+    )