From 29d2200db0af9806cb360b996489065c955a5a5e Mon Sep 17 00:00:00 2001
From: Tom Gustafsson <tom.gustafsson@iki.fi>
Date: Sat, 23 Dec 2023 21:08:28 +0200
Subject: [PATCH] reorient boundary normals to outward normals in ElementGlobal

---
 skfem/element/element_global.py | 14 ++++++++++++++
 tests/test_assembly.py          | 21 +++++++++++++++++++++
 tests/test_examples.py          |  3 ++-
 3 files changed, 37 insertions(+), 1 deletion(-)

diff --git a/skfem/element/element_global.py b/skfem/element/element_global.py
index ad71c034..1ac3c777 100644
--- a/skfem/element/element_global.py
+++ b/skfem/element/element_global.py
@@ -166,6 +166,20 @@ def _eval_dofs(self, mesh, tind=None):
                                         -w['n'][itr, 0, :]])
                 w['n'][itr] /= np.linalg.norm(w['n'][itr], axis=0)
 
+            # swap
+            from skfem.assembly import FacetBasis
+            fb = FacetBasis(mesh, mesh.elem(), intorder=0)
+            ix = np.isin(mesh.t2f, mesh.boundary_facets())
+            sortix = np.argsort(mesh.t2f[ix])
+            norms1 = np.vstack((w['n'][:, 0, :][ix][sortix],
+                                w['n'][:, 1, :][ix][sortix]))
+            norms2 = fb.normals[:, :, 0]
+            signs = np.diag(norms1.T @ norms2)
+            signs1 = signs.copy()
+            signs1[sortix] = signs
+            w['n'][:, 0, :][ix] *= signs1
+            w['n'][:, 1, :][ix] *= signs1
+
         # evaluate dofs, gdof implemented in subclasses
         for itr in range(N):
             for jtr in range(N):
diff --git a/tests/test_assembly.py b/tests/test_assembly.py
index 9ac1fe0b..ca921c3e 100644
--- a/tests/test_assembly.py
+++ b/tests/test_assembly.py
@@ -747,3 +747,24 @@ def test_hcurl_projections_3d():
 
     assert_almost_equal(curly, curlx)
     assert_almost_equal(curly, curla)
+
+
+def test_element_global_boundary_normal():
+
+    mesh = MeshTri().refined(3)
+    basis = Basis(mesh, ElementTriMorley())
+
+    D = basis.get_dofs().all('u_n')
+    x = basis.zeros()
+    x[D] = 1
+
+    @BilinearForm
+    def bilinf(u, v, w):
+        return ddot(u.hess, v.hess)
+
+    A = bilinf.assemble(basis)
+    f = unit_load.assemble(basis)
+
+    y = solve(*condense(A, f, D=D, x=x))
+
+    assert (y[basis.get_dofs(elements=True).all('u')] < 0).all()
diff --git a/tests/test_examples.py b/tests/test_examples.py
index a9fa01d9..fae1453c 100644
--- a/tests/test_examples.py
+++ b/tests/test_examples.py
@@ -16,7 +16,8 @@ class TestEx02(TestCase):
 
     def runTest(self):
         import docs.examples.ex02 as ex02
-        self.assertAlmostEqual(np.max(ex02.x), 0.001217973811129439)
+        self.assertAlmostEqual(np.max(ex02.x[ex02.ib.nodal_dofs[0]]),
+                               0.00033840961095522285)
 
 
 class TestEx03(TestCase):