reorient boundary normals to outward normals in ElementGlobal

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):

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()

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):