From dd1e1ad5fd1ebf2fa5d195806dcb2e3ff10695fe Mon Sep 17 00:00:00 2001 From: Vili Kohonen <43872467+vohonen@users.noreply.github.com> Date: Thu, 13 Jun 2024 14:07:27 +0300 Subject: [PATCH] Fix Mesh.p2e (#1131) Co-authored-by: vohonen --- skfem/mesh/mesh.py | 13 +++++++++++++ tests/test_mesh.py | 22 +++++++++++++++++----- 2 files changed, 30 insertions(+), 5 deletions(-) diff --git a/skfem/mesh/mesh.py b/skfem/mesh/mesh.py index f85a10b1..2f2a5b22 100644 --- a/skfem/mesh/mesh.py +++ b/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]]) diff --git a/tests/test_mesh.py b/tests/test_mesh.py index 906c0eae..de9f74ad 100644 --- a/tests/test_mesh.py +++ b/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():