From cf431d8ad86f79c9c51346b5c6922939a6d8e1c8 Mon Sep 17 00:00:00 2001 From: Tom Gustafsson Date: Thu, 17 Aug 2023 12:30:58 +0300 Subject: [PATCH 1/6] Add ignore_orientation flag --- skfem/io/meshio.py | 67 +++++++++++++++++++++++++--------------------- 1 file changed, 36 insertions(+), 31 deletions(-) diff --git a/skfem/io/meshio.py b/skfem/io/meshio.py index a058cfeb..0163d7fc 100644 --- a/skfem/io/meshio.py +++ b/skfem/io/meshio.py @@ -50,7 +50,8 @@ def from_meshio(m, out=None, int_data_to_sets=False, - force_meshio_type=None): + force_meshio_type=None, + ignore_orientation=False): cells = m.cells_dict meshio_type = None @@ -126,37 +127,41 @@ def from_meshio(m, } 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(): - 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) - 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) + if ignore_orientation: + a = np.array(sorted_facets[k]) + 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) + 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) # MSH 2.2 tag parsing if len(boundaries) == 0 and m.cell_data and m.field_data: From c8320b2589ff35f73c292c3ea01eeb7385262139 Mon Sep 17 00:00:00 2001 From: Tom Gustafsson Date: Wed, 25 Oct 2023 08:43:49 +0300 Subject: [PATCH 2/6] Test ignore_orientation --- skfem/io/meshio.py | 13 +++++++++---- tests/test_mesh.py | 11 +++++++++-- 2 files changed, 18 insertions(+), 6 deletions(-) diff --git a/skfem/io/meshio.py b/skfem/io/meshio.py index 0163d7fc..84877511 100644 --- a/skfem/io/meshio.py +++ b/skfem/io/meshio.py @@ -131,7 +131,9 @@ def from_meshio(m, if ignore_orientation: a = np.array(sorted_facets[k]) b = mtmp.facets.T - boundaries[k] = np.nonzero((a == b[:, None]).all(axis=2).any(axis=1))[0] + boundaries[k] = np.nonzero((a == b[:, None]) + .all(axis=2) + .any(axis=1))[0] else: indices = [] oris = [] @@ -147,12 +149,15 @@ def from_meshio(m, continue if len(f) == 2: # rotate tangent to find normal - tangent = mtmp.p[:, facet[1]] - mtmp.p[:, facet[0]] + 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]] + tangent1 = (mtmp.p[:, facet[1]] + - mtmp.p[:, facet[0]]) + tangent2 = (mtmp.p[:, facet[2]] + - mtmp.p[:, facet[0]]) normal = -np.cross(tangent1, tangent2) else: raise NotImplementedError diff --git a/tests/test_mesh.py b/tests/test_mesh.py index 4542ff1b..ae480620 100644 --- a/tests/test_mesh.py +++ b/tests/test_mesh.py @@ -571,7 +571,14 @@ def test_saveload_cycle_vtk(m): MeshTet(), ] ) -def test_saveload_cycle_tags(fmt, kwargs, m): +@pytest.mark.parametrize( + "ignore_orientation", + [ + True, + False, + ] +) +def test_saveload_cycle_tags(fmt, kwargs, m, ignore_orientation): m = (m .refined(2) @@ -582,7 +589,7 @@ def test_saveload_cycle_tags(fmt, kwargs, m): with NamedTemporaryFile(suffix=fmt) as f: m.save(f.name, point_data={'foo': m.p[0]}, **kwargs) out = ['point_data', 'cells_dict'] - m2 = Mesh.load(f.name, out=out) + m2 = Mesh.load(f.name, out=out, ignore_orientation=ignore_orientation) assert_array_equal(m.p, m2.p) From 836d9d9ac9eb1a7956387cb61432fbdcd8f7275c Mon Sep 17 00:00:00 2001 From: Tom Gustafsson Date: Wed, 25 Oct 2023 08:45:29 +0300 Subject: [PATCH 3/6] Fix syntax error --- skfem/io/meshio.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skfem/io/meshio.py b/skfem/io/meshio.py index 84877511..515f2a0a 100644 --- a/skfem/io/meshio.py +++ b/skfem/io/meshio.py @@ -150,7 +150,7 @@ def from_meshio(m, if len(f) == 2: # rotate tangent to find normal tangent = (mtmp.p[:, facet[1]] - - mtmp.p[:, facet[0]] + - mtmp.p[:, facet[0]]) normal = np.array([-tangent[1], tangent[0]]) elif len(f) == 3: # cross product to find normal From 7fe0f99afd269816b8cc23bfb640d9c9a34d3792 Mon Sep 17 00:00:00 2001 From: Tom Gustafsson Date: Wed, 25 Oct 2023 08:49:28 +0300 Subject: [PATCH 4/6] Ignore types in matplotlib --- skfem/visuals/matplotlib.py | 1 + 1 file changed, 1 insertion(+) diff --git a/skfem/visuals/matplotlib.py b/skfem/visuals/matplotlib.py index 1421a09e..4c7b49d5 100644 --- a/skfem/visuals/matplotlib.py +++ b/skfem/visuals/matplotlib.py @@ -1,3 +1,4 @@ +# type: ignore """Drawing meshes and solutions using matplotlib.""" from functools import singledispatch From 1847358f85679f8290f0b7aa2ad1bc9ef1668d55 Mon Sep 17 00:00:00 2001 From: Tom Gustafsson Date: Wed, 25 Oct 2023 09:02:22 +0300 Subject: [PATCH 5/6] Add another ignore_interior_facets flag --- skfem/io/meshio.py | 10 +++++++--- tests/test_mesh.py | 14 ++++++++++++-- 2 files changed, 19 insertions(+), 5 deletions(-) diff --git a/skfem/io/meshio.py b/skfem/io/meshio.py index 515f2a0a..b3ee2a78 100644 --- a/skfem/io/meshio.py +++ b/skfem/io/meshio.py @@ -51,7 +51,8 @@ def from_meshio(m, out=None, int_data_to_sets=False, force_meshio_type=None, - ignore_orientation=False): + ignore_orientation=False, + ignore_interior_facets=False): cells = m.cells_dict meshio_type = None @@ -128,9 +129,12 @@ def from_meshio(m, 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: + if ignore_orientation or ignore_interior_facets: a = np.array(sorted_facets[k]) - b = mtmp.facets.T + 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] diff --git a/tests/test_mesh.py b/tests/test_mesh.py index 12fa11ae..5f862a9a 100644 --- a/tests/test_mesh.py +++ b/tests/test_mesh.py @@ -578,7 +578,14 @@ def test_saveload_cycle_vtk(m): False, ] ) -def test_saveload_cycle_tags(fmt, kwargs, m, ignore_orientation): +@pytest.mark.parametrize( + "ignore_interior_facets", + [ + True, + False, + ] +) +def test_saveload_cycle_tags(fmt, kwargs, m, ignore_orientation, ignore_interior_facets): m = (m .refined(2) @@ -589,7 +596,10 @@ def test_saveload_cycle_tags(fmt, kwargs, m, ignore_orientation): with NamedTemporaryFile(suffix=fmt) as f: m.save(f.name, point_data={'foo': m.p[0]}, **kwargs) out = ['point_data', 'cells_dict'] - m2 = Mesh.load(f.name, out=out, ignore_orientation=ignore_orientation) + m2 = Mesh.load(f.name, + out=out, + ignore_orientation=ignore_orientation, + ignore_interior_facets=ignore_interior_facets) assert_array_equal(m.p, m2.p) From e7f0722be7e7c9d036e590e26e5784e96416dafe Mon Sep 17 00:00:00 2001 From: Tom Gustafsson Date: Wed, 25 Oct 2023 09:07:48 +0300 Subject: [PATCH 6/6] Add changelog entry --- README.md | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/README.md b/README.md index ad604d7b..05d8e302 100644 --- a/README.md +++ b/README.md @@ -212,6 +212,11 @@ with respect to documented and/or tested features. ### Unreleased - Removed: Python 3.7 support +- Added: `Mesh.load` supports new keyword arguments + `ignore_orientation=True` and `ignore_interior_facets=True` which + will both speed up the loading of larger three-dimensional meshes by + ignoring facet orientation and all tags not on the boundary, + respectively. - Fixed: `MeshTet` uniform refine was reindexing subdomains incorrectly ## [8.1.0] - 2023-06-16