Fix Mesh.p2e (#1131)

Co-authored-by: vohonen <[email protected]>

skfem/mesh/mesh.py

@@ -174,6 +174,19 @@ def p2t(self):
     @property
     def p2e(self):
         """Incidence matrix between edges and vertices."""
+        from scipy.sparse import coo_matrix
+        edges = self.edges.flatten('C')
+        return coo_matrix(
+            (np.ones(len(edges), dtype=np.int32),
+             (np.concatenate((np.arange(self.nedges),)
+                             * self.edges.shape[0]), edges)),
+            shape=(self.nedges, self.nvertices),
+            dtype=np.int32,
+        ).tocsc()
+
+    @property
+    def e2t(self):
+        """Incidence matrix between elements and edges."""
         p2t = self.p2t
         edges = self.edges
         return p2t[:, edges[0]].multiply(p2t[:, edges[1]])

tests/test_mesh.py

@@ -7,8 +7,8 @@
 from numpy.testing import assert_array_equal, assert_almost_equal
 
 from skfem.mesh import (Mesh, MeshHex, MeshLine, MeshQuad, MeshTet, MeshTri,
-                        MeshTri2, MeshQuad2, MeshTet2, MeshHex2, MeshLine1DG,
-                        MeshQuad1DG, MeshHex2, MeshTri1DG)
+                        MeshTet1, MeshHex1, MeshLine1DG, MeshQuad1DG, 
+                        MeshTri1DG, MeshTri2, MeshQuad2, MeshTet2, MeshHex2)
 from skfem.assembly import Basis, LinearForm, Functional, FacetBasis
 from skfem.element import (ElementTetP1, ElementTriP0, ElementQuad0,
                            ElementHex0)
@@ -787,12 +787,24 @@ def test_incidence(mesh):
 
     p2t = mesh.p2t
     for itr in range(0, 50, 3):
-       assert np.sum((mesh.t == itr).any(axis=0)) == len(p2t[:, itr].data)
+        assert np.sum((mesh.t == itr).any(axis=0)) == len(p2t[:, itr].data)
 
     p2f = mesh.p2f
     for itr in range(0, 50, 3):
-       assert np.sum((mesh.facets == itr).any(axis=0)) == len(p2f[:, itr].data)
-
+        assert np.sum((mesh.facets == itr).any(axis=0)) == len(p2f[:, itr].data)
+
+    if isinstance(mesh, (MeshTet1, MeshHex1)):
+        p2e = mesh.p2e
+        for itr in range(0, 50, 3):
+            assert np.sum((mesh.edges == itr).any(axis=0)) == len(p2e[:, itr].data)
+
+        e2t = mesh.e2t
+        for itr in range(0, 50, 3):
+            e = mesh.edges[:, itr]
+            datalen = np.sum(
+                (mesh.t == e[0]).any(axis=0) & (mesh.t == e[1]).any(axis=0)
+            )
+            assert datalen == len(e2t[:, itr].data)
 
 def test_restrict_tags_boundary():