meshio boundaries from cell_sets M*N memory

[jove1]

Function from_meshio 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)

[kinnala]

I will make an attempt to do something about this in #1117.

[kinnala]

I have added another approach in #1117 which is used if ignore_orientation=True is set. I did not have time to test it yet using a large mesh.

[kinnala]

Now I was able to test it. Tetrahedral mesh with 100 000 vertices and 500 000 elements. Has a surface with roughly 6000 facets. It used 0.5 gigabytes of memory and it took around 5 seconds to load with Mesh.load(..., ignore_orientation=True). Without the flag it kept running and running..

[kinnala]

I have now extended this to oriented boundaries, ignore_orientation=True will still shave off half a second in my test case.