Repository: kinnala/scikit-fem

Issue:
meshio boundaries from cell_sets M*N memory
Function [`from_meshio`](https://github.com/kinnala/scikit-fem/blob/a0d679ed47a07637159431aa91d498905734436e/skfem/io/meshio.py#L138) uses M*N memory when matching facets from cell_set to facets of a mesh.
This leads to out-of-memory error on my real world mesh (148494 facets, 16316 on boundary, 10694 in one of the boundary regions).
Problem is alleviated by using ignore_interior_facets=True,  (I have enough memory for (16316, 10694, 3) but not for (148494, 10694, 3)), but it is still slow and I think it is better to use a hash map (dict):

```
import numpy as np
import meshio
import skfem as sf
import time

# this gets out-of-memory killed
# m = sf.MeshTet.load("model-Cut.msh", ignore_orientation=True) 

# this takes long
m = sf.MeshTet.load("model-Cut.msh", ignore_orientation=True, ignore_interior_facets=True)
print(m)

# try to replicate step by step
m = meshio.read("model-Cut.msh")
print(m)

p = np.ascontiguousarray(m.points.T)
t = np.ascontiguousarray(m.cells_dict["tetra"].T)

mtmp = sf.MeshTet(p, t)

sorted_mesh_cells = np.sort(m.cells_dict['triangle']) # by default sorts last axis

b = mtmp.facets.T
#b = mtmp.facets[:, mtmp.boundary_facets()].T
facet_map = { tuple(f):i for i,f in enumerate(b) } # create lookup table

boundaries = {}
print()
for k, v in m.cell_sets_dict.items():
    if 'triangle' in v and k.split(":")[0] != "gmsh":
        a = sorted_mesh_cells[v['triangle']]

        print(k, a.shape, b.shape)

        # this needs M*N memory !!!
        #print( (a == b[:,None]).shape)

        # use lookup table
        boundaries[k] = np.array([facet_map[tuple(f)] for f in a])
print()

m = sf.MeshTet(p, t, boundaries)
print(m)
```

Discussion:

Patch:
diff --git a/README.md b/README.md
index 66ee8394..a8407ed8 100644
--- a/README.md
+++ b/README.md
@@ -212,8 +212,11 @@ with respect to documented and/or tested features.
 ### Unreleased
 
 - Fixed: Tests now pass with `numpy==2.0rc1`
-- Fixed: `MappingAffine` uses lazy evaluation also for element mappings
-- Fixed: `MeshTet.init_tensor` uses less memory for large meshes
+- Fixed: `MappingAffine` now uses lazy evaluation also for element
+  mappings, in addition to boundary mappings
+- Fixed: `MeshTet.init_tensor` uses significantly less memory for
+  large meshes
+- Fixed: `Mesh.load` uses less memory when loading and matching tags
 - Added: `Basis` has new optional `disable_doflocs` kwarg which
   can be set to `True` to avoid computing `Basis.doflocs`, to reduce memory
   usage
diff --git a/docs/examples/meshes/tagged_gmsh4.msh b/docs/examples/meshes/tagged_gmsh4.msh
new file mode 100644
index 00000000..b7ee83a5
--- /dev/null
+++ b/docs/examples/meshes/tagged_gmsh4.msh
@@ -0,0 +1,240 @@
+$MeshFormat
+4.1 0 8
+$EndMeshFormat
+$PhysicalNames
+3
+1 6 "tagged"
+1 7 "test"
+2 8 "all"
+$EndPhysicalNames
+$Entities
+5 5 1 0
+1 -0.5 -0.5 0 0 
+2 0.5 -0.5 0 0 
+3 0.5 -0.3 0 0 
+4 0 -0.3 0 0 
+5 0 1.3 0 0 
+1 0.5 -0.5 0 0.5 -0.3 0 0 2 2 -3 
+2 0 -0.3 0 0.5 -0.3 0 0 2 3 -4 
+3 0 -0.3 0 0 1.3 0 2 6 7 2 4 -5 
+4 -0.5 -0.5 0 0 1.3 0 0 2 5 -1 
+5 -0.5 -0.5 0 0.5 -0.5 0 0 2 1 -2 
+1 -0.5 -0.5 0 0.5 1.3 0 1 8 5 4 5 1 2 3 
+$EndEntities
+$Nodes
+11 55 1 55
+0 1 0 1
+10
+-0.5 -0.5 0
+0 2 0 1
+11
+0.5 -0.5 0
+0 3 0 1
+12
+0.5 -0.3 0
+0 4 0 1
+1
+0 -0.3 0
+0 5 0 1
+2
+0 1.3 0
+1 1 0 3
+13
+14
+15
+0.5 -0.4500000000000726 0
+0.5 -0.4000000000002213 0
+0.5 -0.3500000000002064 0
+1 2 0 3
+16
+17
+18
+0.3750000000002064 -0.3 0
+0.2500000000006652 -0.3 0
+0.1250000000003326 -0.3 0
+1 3 0 7
+3
+4
+5
+6
+7
+8
+9
+0 -0.1000000000003512 0
+0 0.09999999999904252 0
+0 0.2999999999984534 0
+0 0.4999999999978643 0
+0 0.6999999999982724 0
+0 0.8999999999988308 0
+0 1.099999999999313 0
+1 4 0 7
+19
+20
+21
+22
+23
+24
+25
+-0.06249999999991335 1.075000000000312 0
+-0.1249999999997748 0.8500000000008107 0
+-0.1874999999995735 0.6250000000015353 0
+-0.2499999999993722 0.40000000000226 0
+-0.3124999999996163 0.1750000000013814 0
+-0.3749999999996823 -0.04999999999885629 0
+-0.4374999999998545 -0.2749999999994761 0
+1 5 0 3
+26
+27
+28
+-0.2500000000004945 -0.5 0
+-1.312838726619248e-12 -0.5 0
+0.2499999999993436 -0.5 0
+2 1 0 27
+29
+30
+31
+32
+33
+34
+35
+36
+37
+38
+39
+40
+41
+42
+43
+44
+45
+46
+47
+48
+49
+50
+51
+52
+53
+54
+55
+0.25 -0.3999999999999999 0
+-0.25 -0.3999999999999999 0
+-0.1249999999996861 0.4500000000000621 0
+-0.1249999999996861 0.05000000000112997 0
+0.3750000000003326 -0.3500000000001105 0
+0.2500000000003326 -0.3499999999999999 0
+0.375 -0.4000000000001105 0
+0.125 -0.3499999999999999 0
+0.375 -0.4499999999999999 0
+0.1249999999993436 -0.4499999999999999 0
+0 -0.3999999999999999 0
+-0.1250000000006564 -0.4499999999999999 0
+-0.125 -0.3499999999999999 0
+-0.375 -0.4499999999999999 0
+-0.1874999999995292 0.4250000000011609 0
+-0.1249999999997305 0.6500000000004362 0
+-0.06249999999984306 0.6749999999994464 0
+-0.06249999999988741 0.8749999999998205 0
+-0.06249999999984306 0.4749999999989631 0
+-0.06249999999984306 0.07500000000008623 0
+-0.06249999999984306 -0.124999999999435 0
+-0.06249999999984306 0.2749999999995523 0
+-0.1249999999996861 0.250000000000596 0
+-0.1874999999995292 0.2250000000016949 0
+-0.3124999999998411 -0.2249999999994281 0
+-0.1874999999998431 -0.1749999999994349 0
+-0.2499999999996842 1.136840621640544e-12 0
+$EndNodes
+$Elements
+2 88 1 113
+1 3 1 8
+1 1 3 
+2 3 4 
+3 4 5 
+4 5 6 
+5 6 7 
+6 7 8 
+7 8 9 
+8 9 2 
+2 1 2 80
+34 12 16 15 
+35 16 33 15 
+36 16 17 33 
+37 15 33 14 
+38 17 34 33 
+39 34 35 33 
+40 34 29 35 
+41 33 35 14 
+42 17 18 34 
+43 18 36 34 
+44 18 1 36 
+45 34 36 29 
+46 14 35 13 
+47 35 37 13 
+48 35 29 37 
+49 13 37 11 
+50 11 37 28 
+51 37 38 28 
+52 37 29 38 
+53 28 38 27 
+54 29 39 38 
+55 39 40 38 
+56 39 30 40 
+57 38 40 27 
+58 29 36 39 
+59 36 41 39 
+60 36 1 41 
+61 39 41 30 
+62 27 40 26 
+63 40 42 26 
+64 40 30 42 
+65 26 42 10 
+66 22 43 21 
+67 43 44 21 
+68 43 31 44 
+69 21 44 20 
+70 31 45 44 
+71 45 46 44 
+72 45 8 46 
+73 44 46 20 
+74 31 47 45 
+75 47 7 45 
+76 47 6 7 
+77 45 7 8 
+78 20 46 19 
+79 46 9 19 
+80 46 8 9 
+81 19 9 2 
+82 1 3 49 
+83 3 48 49 
+84 3 4 48 
+85 49 48 32 
+86 4 50 48 
+87 50 51 48 
+88 50 31 51 
+89 48 51 32 
+90 4 5 50 
+91 5 47 50 
+92 5 6 47 
+93 50 47 31 
+94 32 51 52 
+95 51 43 52 
+96 51 31 43 
+97 52 43 22 
+98 10 42 25 
+99 42 53 25 
+100 42 30 53 
+101 25 53 24 
+102 30 54 53 
+103 54 55 53 
+104 54 32 55 
+105 53 55 24 
+106 30 41 54 
+107 41 49 54 
+108 41 1 49 
+109 54 49 32 
+110 24 55 23 
+111 55 52 23 
+112 55 32 52 
+113 23 52 22 
+$EndElements
diff --git a/skfem/io/meshio.py b/skfem/io/meshio.py
index f52e9031..930d2aa4 100644
--- a/skfem/io/meshio.py
+++ b/skfem/io/meshio.py
@@ -1,5 +1,6 @@
 """Import/export any formats supported by meshio."""
 
+import logging
 from pathlib import Path
 
 import meshio
@@ -10,6 +11,9 @@ from skfem.mesh import (MeshTet1, MeshTet2, MeshHex1, MeshHex2, MeshWedge1,
 from skfem.generic_utils import OrientedBoundary
 
 
+logger = logging.getLogger(__name__)
+
+
 MESH_TYPE_MAPPING = {
     'tetra': MeshTet1,
     'tetra10': MeshTet2,
@@ -54,6 +58,9 @@ def from_meshio(m,
                 ignore_orientation=False,
                 ignore_interior_facets=False):
 
+    if ignore_interior_facets:
+        logger.warning("kwarg ignore_interior_facets is unused.")
+
     cells = m.cells_dict
     meshio_type = None
 
@@ -121,56 +128,42 @@ def from_meshio(m,
 
     # parse boundaries from cell_sets
     if m.cell_sets and bnd_type in m.cells_dict:
-        oriented_facets = {
-            k: [tuple(f) for f in m.cells_dict[bnd_type][v[bnd_type]]]
-            for k, v in m.cell_sets_dict.items()
-            if bnd_type in v and k.split(":")[0] != "gmsh"
-        }
-        sorted_facets = {k: [tuple(np.sort(f)) for f in v]
-                         for k, v in oriented_facets.items()}
-        for k, v in oriented_facets.items():
-            if ignore_orientation or ignore_interior_facets:
-                a = np.array(sorted_facets[k])
-                if ignore_interior_facets:
-                    b = mtmp.facets[:, mtmp.boundary_facets()].T
-                else:
-                    b = mtmp.facets.T
-                boundaries[k] = np.nonzero((a == b[:, None])
-                                           .all(axis=2)
-                                           .any(axis=1))[0]
-            else:
-                indices = []
-                oris = []
-                for i, f in enumerate(map(tuple, mtmp.facets.T)):
-                    if f in sorted_facets[k]:
-                        indices.append(i)
-                        ix = sorted_facets[k].index(f)
-                        facet = v[ix]
-                        t1, t2 = mtmp.f2t[:, i]
-                        if t2 == -1:
-                            # orientation on boundary is 0
-                            oris.append(0)
-                            continue
-                        if len(f) == 2:
-                            # rotate tangent to find normal
-                            tangent = (mtmp.p[:, facet[1]]
-                                       - mtmp.p[:, facet[0]])
-                            normal = np.array([-tangent[1], tangent[0]])
-                        elif len(f) == 3:
-                            # cross product to find normal
-                            tangent1 = (mtmp.p[:, facet[1]]
-                                        - mtmp.p[:, facet[0]])
-                            tangent2 = (mtmp.p[:, facet[2]]
-                                        - mtmp.p[:, facet[0]])
-                            normal = -np.cross(tangent1, tangent2)
+        p2f = mtmp.p2f
+        for k, v in m.cell_sets_dict.items():
+            if bnd_type in v and k.split(":")[0] != "gmsh":
+                facets = m.cells_dict[bnd_type][v[bnd_type]].T
+                sorted_facets = np.sort(facets, axis=0)
+                ind = p2f[:, sorted_facets[0]]
+                for itr in range(sorted_facets.shape[0] - 1):
+                    ind = ind.multiply(p2f[:, sorted_facets[itr + 1]])
+                boundaries[k] = np.nonzero(ind)[0]
+
+                if not ignore_orientation:
+                    try:
+                        ori = np.zeros_like(boundaries[k], dtype=np.float64)
+                        t1, _ = mtmp.f2t[:, boundaries[k]]
+                        if facets.shape[0] == 2:
+                            tangents = (mtmp.p[:, facets[1]]
+                                        - mtmp.p[:, facets[0]])
+                            normals = np.array([-tangents[1], tangents[0]])
+                        elif facets.shape[0] == 3:
+                            tangents1 = (mtmp.p[:, facets[1]]
+                                         - mtmp.p[:, facets[0]])
+                            tangents2 = (mtmp.p[:, facets[2]]
+                                         - mtmp.p[:, facets[0]])
+                            normals = -np.cross(tangents1.T, tangents2.T).T
                         else:
                             raise NotImplementedError
-                        # find another vector using node outside of boundary
-                        third = np.setdiff1d(mtmp.t[:, t1],
-                                             np.array(f))[0]
-                        outplane = mtmp.p[:, f[0]] - mtmp.p[:, third]
-                        oris.append(1 if np.dot(normal, outplane) > 0 else 0)
-                boundaries[k] = OrientedBoundary(indices, oris)
+                        for itr in range(mtmp.t.shape[0]):
+                            ori += np.sum(normals
+                                          * (mtmp.p[:, facets[0]]
+                                             - mtmp.p[:, mtmp.t[itr, t1]]),
+                                          axis=0)
+                        ori = 1 * (ori > 0)
+                        boundaries[k] = OrientedBoundary(boundaries[k],
+                                                         ori)
+                    except Exception:
+                        logger.warning("Failure to orient a boundary.")
 
     # MSH 2.2 tag parsing
     if len(boundaries) == 0 and m.cell_data and m.field_data:
@@ -212,7 +205,7 @@ def from_meshio(m,
                 boundaries[find_tagname(tag)] = index[tagindex, 1]
 
         except Exception:
-            pass
+            logger.warning("Failure to parse tags from meshio.")
 
     # attempt parsing skfem tags
     if m.cell_data:
diff --git a/skfem/mesh/mesh.py b/skfem/mesh/mesh.py
index 4b1d6e0f..97d891ec 100644
--- a/skfem/mesh/mesh.py
+++ b/skfem/mesh/mesh.py
@@ -133,6 +133,51 @@ class Mesh:
             self._init_edges()
         return self._t2e
 
+    @property
+    def p2f(self):
+        """Incidence matrix between facets and vertices.
+
+        Examples
+        --------
+        To find facets with point 5:
+
+        >>> import numpy as np
+        >>> from skfem import MeshTet
+        >>> mesh = MeshTet().refined()
+        >>> np.nonzero(mesh.p2f[:, 5])[0]
+        array([33, 34, 35], dtype=int32)
+
+        """
+        from scipy.sparse import coo_matrix
+        facets = self.facets.flatten('C')
+        return coo_matrix(
+            (np.ones(len(facets), dtype=np.int32),
+             (np.concatenate((np.arange(self.nfacets),)
+                             * self.facets.shape[0]), facets)),
+            shape=(self.nfacets, self.nvertices),
+            dtype=np.int32,
+        ).tocsc()
+
+    @property
+    def p2t(self):
+        """Incidence matrix between elements and vertices."""
+        from scipy.sparse import coo_matrix
+        t = self.t.flatten('C')
+        return coo_matrix(
+            (np.ones(len(t), dtype=np.int32),
+             (np.concatenate((np.arange(self.nelements),)
+                             * self.nnodes), t)),
+            shape=(self.nelements, self.nvertices),
+            dtype=np.int32,
+        ).tocsc()
+
+    @property
+    def p2e(self):
+        """Incidence matrix between edges and vertices."""
+        p2t = self.p2t
+        edges = self.edges
+        return p2t[:, edges[0]].multiply(p2t[:, edges[1]])
+
     def dim(self):
         return self.elem.refdom.dim()
 


Test patch:
diff --git a/tests/test_mesh.py b/tests/test_mesh.py
index 95cdc79b..b7db4129 100644
--- a/tests/test_mesh.py
+++ b/tests/test_mesh.py
@@ -512,6 +512,21 @@ def test_meshio_cycle_boundaries(boundaries_only, m):
                            m.boundaries[key].sort())
 
 
+@pytest.mark.parametrize(
+    "mtype, path, ignore_orientation",
+    [
+        (MeshTet, MESH_PATH / 'box.msh', False),
+        (MeshTet, MESH_PATH / 'box.msh', True),
+        (MeshTri, MESH_PATH / 'tagged_gmsh4.msh', False),
+        (MeshTri, MESH_PATH / 'tagged_gmsh4.msh', True),
+    ]
+)
+def test_load_file(mtype, path, ignore_orientation):
+
+    m = mtype.load(path, ignore_orientation=ignore_orientation)
+    assert len(m.boundaries) > 0
+
+
 @pytest.mark.parametrize(
     "m",
     [
@@ -555,7 +570,7 @@ def test_saveload_cycle_vtk(m):
 @pytest.mark.parametrize(
     "fmt, kwargs",
     [
-        ('.msh', {}),
+        ('.msh', {'file_format': 'gmsh'}),
         ('.msh', {'file_format': 'gmsh22'}),
         ('.vtk', {}),
         #('.xdmf', {}),
@@ -757,3 +772,23 @@ def test_remove_duplicates():
 
     assert not m1.is_valid()
     assert m1.remove_duplicate_nodes().is_valid()
+
+
+@pytest.mark.parametrize(
+    "mesh",
+    [
+        MeshTri().refined(3),
+        MeshTet().refined(3),
+        MeshQuad().refined(3),
+        MeshHex().refined(3),
+    ]
+)
+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)
+
+    p2f = mesh.p2f
+    for itr in range(0, 50, 3):
+       assert np.sum((mesh.facets == itr).any(axis=0)) == len(p2f[:, itr].data)