Tools for matching topological entities (#1117)

* Add Mesh.p2f, Mesh.p2e, Mesh.p2t * use p2f to match boundaries in from_meshio * optimize also loading of oriented boundaries
kinnala · Apr 13, 2024 · 9c87116 · 9c87116
1 parent 5ce46d9
commit 9c87116
Show file tree

Hide file tree

Showing 5 changed files with 368 additions and 52 deletions.
diff --git 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
@@ -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
@@ -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.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 find_tagname(tag):
                 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: