diff --git a/README.md b/README.md index 66ee83945..a8407ed8f 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 000000000..b7ee83a56 --- /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 f52e90313..930d2aa45 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.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: diff --git a/skfem/mesh/mesh.py b/skfem/mesh/mesh.py index 4b1d6e0f3..97d891eca 100644 --- a/skfem/mesh/mesh.py +++ b/skfem/mesh/mesh.py @@ -133,6 +133,51 @@ def t2e(self): 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() diff --git a/tests/test_mesh.py b/tests/test_mesh.py index 95cdc79b3..b7db4129f 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)