Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Tools for matching topological entities #1117

Merged
merged 8 commits into from
Apr 13, 2024
Merged

Tools for matching topological entities #1117

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
0c22d8f
Add Mesh.f2p
kinnala Apr 9, 2024
2239115
Rename to p2f and add p2t, tests
kinnala Apr 10, 2024
71748d4
add p2e
kinnala Apr 12, 2024
3f86ddd
add another test and tagged gmsh4 file
kinnala Apr 12, 2024
ee34da2
use p2f to match boundaries in from_meshio
kinnala Apr 12, 2024
e9f365c
optimize also loading of oriented boundaries
kinnala Apr 13, 2024
33ac627
Update changelog
kinnala Apr 13, 2024
35fca29
fix flake8
kinnala Apr 13, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 5 additions & 2 deletions README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
240 changes: 240 additions & 0 deletions docs/examples/meshes/tagged_gmsh4.msh
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
91 changes: 42 additions & 49 deletions skfem/io/meshio.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
"""Import/export any formats supported by meshio."""

import logging
from pathlib import Path

import meshio
Expand All @@ -10,6 +11,9 @@
from skfem.generic_utils import OrientedBoundary


logger = logging.getLogger(__name__)


MESH_TYPE_MAPPING = {
'tetra': MeshTet1,
'tetra10': MeshTet2,
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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:
Expand Down
Loading
Loading