Repository: kinnala/scikit-fem

Issue:
Incorrect edges returned by `Mesh3D.boundary_edges`
I've been working on applying the block-diagonally preconditioned MINRES+AMG Stokes solver of #342 & ex32 to other MeshTet geometries and found some odd errors.  It seems to have nothing to do with the operators but something to do with `.find_dofs` when there are several Gmsh elementary entities in a single Gmsh `Physical Surface`.  The situation here is three Physical Surfaces: inlet, outlet, and wall, with uniform positive pressure applied at the inlet, zero pressure at the outlet, and zero velocity on the wall.  Typically the inlet and outlet are single elementary surfaces and the walls include everything else.  Example: a cube with inlet and outlet on an opposing pair of faces and the other four being walls.

I generate a tetrahedral mesh with marked boundaries in Gmsh, marking the boundaries in two different ways: first with all four walls combined, second with a separate Physical Surface for each.
```
SetFactory("OpenCASCADE");

Box(1) = {0, 0, 0, 1, 1, 1};

Physical Volume("fluid", 13) = {1};

Physical Surface("inlet") = {1};
Physical Surface("outlet") = {2};

Physical Surface("wall") = {3, 4, 5, 6};

Physical Surface("wall3") = {3};
Physical Surface("wall4") = {4};
Physical Surface("wall5") = {5};
Physical Surface("wall6") = {6};
```

I load this with 
```
import numpy as np
import skfem

basis = skfem.InteriorBasis(skfem.MeshTet.load("duct.msh"), skfem.ElementTetP2())

walls = {"separate": np.concatenate([b.all() for k, b in basis.find_dofs().items()
                                     if k.startswith("wall") and k[len("wall") :]]),
         "combined": np.concatenate([b.all() for k, b in basis.find_dofs().items()
                                     if k.endswith("wall")])}
```

Now I was thinking that `walls["separate"]` and `walls["combined"]` would be equal arrays of degrees of freedom on the walls but apparently not.

```python
walls["diff"] = np.setdiff1d(walls["combined"], walls["separate"])
print('In "separate" but not "combined"',  np.setdiff1d(walls["separate"], walls["combined"]))  # empty
print('In "combined" but not "separate"', walls["diff"])  # several
print('In "combined", not "separate", and not edge', np.setdiff1d(walls["diff"], basis.edge_dofs))  # empty
```

The "combined" approach has more degrees of freedom on the walls, with the extra ones all being in `basis.edge_dofs`.  (I haven't located them more precisely yet.)  The result is typically a singular matrix.  The "separated" approach seems to lead to a correct matrix

I haven't ascertained whether the discrepancy is in Gmsh, meshio, or `skfem.MeshTet.load` (or in my understanding of `basis.find_dofs`).

Discussion:

Patch:
diff --git a/README.md b/README.md
index 441c9c98..6af9da4b 100644
--- a/README.md
+++ b/README.md
@@ -109,6 +109,13 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### [Unreleased]
 
+#### Added
+- `ElementTetMini`, MINI-element for tetrahedral mesh
+
+#### Fixed
+- `Mesh3d.boundary_edges` incorrectly returned all edges where both nodes are on
+the boundary
+
 ### [1.0.0] - 2020-04-22
 
 #### Added
diff --git a/docs/examples/ex04.py b/docs/examples/ex04.py
index 97850549..5bcafc13 100644
--- a/docs/examples/ex04.py
+++ b/docs/examples/ex04.py
@@ -1,4 +1,4 @@
-"""Contact problem.
+r"""Contact problem.
 
 Mortar methods allow setting interface conditions on non-matching meshes.
 They are useful also when solving variational inequalities such as
diff --git a/docs/examples/ex05.py b/docs/examples/ex05.py
index 0b4026ad..f7b0dd33 100644
--- a/docs/examples/ex05.py
+++ b/docs/examples/ex05.py
@@ -1,4 +1,4 @@
-"""Integral condition.
+r"""Integral condition.
 
 This short example demonstrates the implementation of an integral boundary
  condition
diff --git a/docs/examples/ex06.py b/docs/examples/ex06.py
index 5dc20b69..3ec46d36 100644
--- a/docs/examples/ex06.py
+++ b/docs/examples/ex06.py
@@ -1,4 +1,4 @@
-"""High-order plotting.
+r"""High-order plotting.
 
 This simple example demonstrates the usage and visualisation of biquadratic
 finite element basis. Many plotting tools, including matplotlib, provide tools
diff --git a/skfem/mesh/mesh2d/mesh_tri.py b/skfem/mesh/mesh2d/mesh_tri.py
index ef835f4c..c0d5b383 100644
--- a/skfem/mesh/mesh2d/mesh_tri.py
+++ b/skfem/mesh/mesh2d/mesh_tri.py
@@ -273,7 +273,7 @@ class MeshTri(Mesh2D):
         self.t2f = ixb.reshape((3, self.t.shape[1]))
 
         # build facet-to-triangle mapping: 2 (triangles) x Nedges
-        e_tmp = np.hstack((self.t2f[0, :], self.t2f[1, :], self.t2f[2, :]))
+        e_tmp = np.hstack((self.t2f[0], self.t2f[1], self.t2f[2]))
         t_tmp = np.tile(np.arange(self.t.shape[1]), (1, 3))[0]
 
         e_first, ix_first = np.unique(e_tmp, return_index=True)
@@ -286,7 +286,7 @@ class MeshTri(Mesh2D):
         self.f2t[1, e_last] = t_tmp[ix_last]
 
         # second row to zero if repeated (i.e., on boundary)
-        self.f2t[1, np.nonzero(self.f2t[0, :] == self.f2t[1, :])[0]] = -1
+        self.f2t[1, np.nonzero(self.f2t[0] == self.f2t[1])[0]] = -1
 
     def _uniform_refine(self):
         """Perform a single mesh refine."""
@@ -298,23 +298,23 @@ class MeshTri(Mesh2D):
         t2f = self.t2f + sz
 
         # new vertices are the midpoints of edges
-        new_p = 0.5 * np.vstack((p[0, e[0, :]] + p[0, e[1, :]],
-                                 p[1, e[0, :]] + p[1, e[1, :]]))
+        new_p = 0.5 * np.vstack((p[0, e[0]] + p[0, e[1]],
+                                 p[1, e[0]] + p[1, e[1]]))
         self.p = np.hstack((p, new_p))
 
         # build new triangle definitions
         self.t = np.hstack((
-            np.vstack((t[0, :], t2f[0, :], t2f[2, :])),
-            np.vstack((t[1, :], t2f[0, :], t2f[1, :])),
-            np.vstack((t[2, :], t2f[2, :], t2f[1, :])),
-            np.vstack((t2f[0, :], t2f[1, :], t2f[2, :])),
+            np.vstack((t[0], t2f[0], t2f[2])),
+            np.vstack((t[1], t2f[0], t2f[1])),
+            np.vstack((t[2], t2f[2], t2f[1])),
+            np.vstack((t2f[0], t2f[1], t2f[2])),
         ))
 
         # mapping of indices between old and new facets
         new_facets = np.zeros((2, e.shape[1]), dtype=np.int64)
         ix0 = np.arange(t.shape[1], dtype=np.int64)
         ix1 = ix0 + t.shape[1]
-        ix2 = ix0 + 2*t.shape[1]
+        ix2 = ix0 + 2 * t.shape[1]
 
         # rebuild mappings
         self._build_mappings()
@@ -334,12 +334,12 @@ class MeshTri(Mesh2D):
 
         def sort_mesh(p, t):
             """Make (0, 2) the longest edge in t."""
-            l01 = np.sqrt(np.sum((p[:, t[0, :]] - p[:, t[1, :]])**2, axis=0))
-            l12 = np.sqrt(np.sum((p[:, t[1, :]] - p[:, t[2, :]])**2, axis=0))
-            l02 = np.sqrt(np.sum((p[:, t[0, :]] - p[:, t[2, :]])**2, axis=0))
+            l01 = np.sqrt(np.sum((p[:, t[0]] - p[:, t[1]]) ** 2, axis=0))
+            l12 = np.sqrt(np.sum((p[:, t[1]] - p[:, t[2]]) ** 2, axis=0))
+            l02 = np.sqrt(np.sum((p[:, t[0]] - p[:, t[2]]) ** 2, axis=0))
 
-            ix01 = (l01 > l02)*(l01 > l12)
-            ix12 = (l12 > l01)*(l12 > l02)
+            ix01 = (l01 > l02) * (l01 > l12)
+            ix12 = (l12 > l01) * (l12 > l02)
 
             # row swaps
             tmp = t[2, ix01]
@@ -372,11 +372,11 @@ class MeshTri(Mesh2D):
             ix[facets == 1] = np.arange(np.count_nonzero(facets)) + m.p.shape[1]
             ix = ix[m.t2f] # (0, 1) (1, 2) (0, 2)
 
-            red =   (ix[0, :] >= 0) * (ix[1, :] >= 0) * (ix[2, :] >= 0)
-            blue1 = (ix[0, :] ==-1) * (ix[1, :] >= 0) * (ix[2, :] >= 0)
-            blue2 = (ix[0, :] >= 0) * (ix[1, :] ==-1) * (ix[2, :] >= 0)
-            green = (ix[0, :] ==-1) * (ix[1, :] ==-1) * (ix[2, :] >= 0)
-            rest =  (ix[0, :] ==-1) * (ix[1, :] ==-1) * (ix[2, :] ==-1)
+            red =   (ix[0] >= 0) * (ix[1] >= 0) * (ix[2] >= 0)
+            blue1 = (ix[0] ==-1) * (ix[1] >= 0) * (ix[2] >= 0)
+            blue2 = (ix[0] >= 0) * (ix[1] ==-1) * (ix[2] >= 0)
+            green = (ix[0] ==-1) * (ix[1] ==-1) * (ix[2] >= 0)
+            rest =  (ix[0] ==-1) * (ix[1] ==-1) * (ix[2] ==-1)
 
             # new red elements
             t_red = np.hstack((
@@ -425,6 +425,6 @@ class MeshTri(Mesh2D):
     def element_finder(self):
         from matplotlib.tri import Triangulation
 
-        return Triangulation(self.p[0, :],
-                             self.p[1, :],
+        return Triangulation(self.p[0],
+                             self.p[1],
                              self.t.T).get_trifinder()
diff --git a/skfem/mesh/mesh3d/mesh3d.py b/skfem/mesh/mesh3d/mesh3d.py
index 2b3baf8c..c20afa3c 100644
--- a/skfem/mesh/mesh3d/mesh3d.py
+++ b/skfem/mesh/mesh3d/mesh3d.py
@@ -32,8 +32,14 @@ class Mesh3D(Mesh):
 
     def boundary_edges(self) -> ndarray:
         """Return an array of boundary edge indices."""
-        return np.nonzero(np.isin(self.edges,
-                                  self.boundary_nodes()).all(axis=0))[0]
+        facets = self.boundary_facets()
+        boundary_edges = np.hstack((
+            self.facets[0:2, facets],
+            self.facets[1:3, facets],
+            self.facets[0:3:2, facets],
+        )).T
+        return np.nonzero((self.edges.T[:, None] == boundary_edges)\
+                          .all(-1).any(-1))[0]
 
     def interior_edges(self) -> ndarray:
         """Return an array of interior edge indices."""


Test patch:
diff --git a/tests/test_basis.py b/tests/test_basis.py
index 8d61a27c..93520994 100644
--- a/tests/test_basis.py
+++ b/tests/test_basis.py
@@ -63,12 +63,19 @@ class TestFacetExpansion(TestCase):
     def runTest(self):
 
         m = MeshTet()
+        m.refine(2)
 
         basis = InteriorBasis(m, ElementTetP2())
 
-        arr1 = basis.find_dofs({
-            'kek':m.facets_satisfying(lambda x: x[0] == 0)
-        })['kek'].edge['u']
-        arr2 = basis.edge_dofs[:, m.edges_satisfying(lambda x: x[0] == 0)]
-
-        assert_allclose(arr1, arr2.flatten())
+        for fun in [lambda x: x[0] == 0,
+                    lambda x: x[0] == 1,
+                    lambda x: x[1] == 0,
+                    lambda x: x[1] == 1,
+                    lambda x: x[2] == 0,
+                    lambda x: x[2] == 1]:
+            arr1 = basis.find_dofs({
+                'kek': m.facets_satisfying(fun)
+            })['kek'].edge['u']
+            arr2 = basis.edge_dofs[:, m.edges_satisfying(fun)]
+
+            assert_allclose(arr1, arr2.flatten())
diff --git a/tests/test_convergence.py b/tests/test_convergence.py
index cc2483c3..7a662592 100644
--- a/tests/test_convergence.py
+++ b/tests/test_convergence.py
@@ -224,9 +224,9 @@ class ConvergenceTetP1(ConvergenceQ1):
 
 
 class ConvergenceTetP2(ConvergenceTetP1):
-    rateL2 = 3.0
-    rateH1 = 2.0
-    eps = 0.132
+    rateL2 = 3.32
+    rateH1 = 2.16
+    eps = 0.01
 
     def create_basis(self, m):
         e = ElementTetP2()
diff --git a/tests/test_examples.py b/tests/test_examples.py
index a48ab7a0..6c27a5d9 100644
--- a/tests/test_examples.py
+++ b/tests/test_examples.py
@@ -228,7 +228,7 @@ class TestEx33(TestCase):
 
     def runTest(self):
         from docs.examples.ex33 import x
-        self.assertAlmostEqual(np.max(x), 0.1240571191753624, delta=1e-8)
+        self.assertAlmostEqual(np.max(x), 0.12412747286589934, delta=1e-8)
 
 
 class TestEx34(TestCase):
diff --git a/tests/test_mesh.py b/tests/test_mesh.py
index 030178d6..6646b1cb 100644
--- a/tests/test_mesh.py
+++ b/tests/test_mesh.py
@@ -141,5 +141,18 @@ class SerializeUnserializeCycle(unittest.TestCase):
                 self.assertTrue((m.subdomains[k] == M.subdomains[k]).all())
 
 
+class TestBoundaryEdges(unittest.TestCase):
+
+    def runTest(self):
+        m = MeshTet()
+        # default mesh has all edges on the boundary
+        self.assertTrue(len(m.boundary_edges()) == m.edges.shape[1])
+        m.refine()
+        # check that there is a correct amount of boundary edges:
+        # 12 (cube edges) * 2 (per cube edge) + 6 (cube faces) * 8 (per cube face)
+        # = 72 edges
+        self.assertTrue(len(m.boundary_edges()) == 72)
+
+
 if __name__ == '__main__':
     unittest.main()