-
Notifications
You must be signed in to change notification settings - Fork 1
/
checkMesh.py
93 lines (78 loc) · 2.19 KB
/
checkMesh.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
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
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
import numpy as np
import sys
import meshio
def check(V, E):
"""Check the orientation of each simplex in the mesh.
Parameters
----------
V : ndarray
n x 2 or n x 3 list of coordinates
E : ndarray
n x 3 or n x 4 list of vertices
"""
m = E.shape[1] # dim = m - 1
G = np.zeros((m-1, m-1)) # shape matrix
D = np.zeros(E.shape[0]) # determinants
print(f'Checking {m-1}D mesh')
for j, el in enumerate(E):
G[:,:] = V[el[1:],:m-1] - V[el[0],:m-1]
D[j] = np.linalg.det(G)
Ipos = np.where(D > 0)[0]
Ineg = np.where(D < 0)[0]
Izer = np.where(np.abs(D) < 1e-14)[0]
return Ipos, Ineg, Izer
def test_check():
# 3D (pos)
V = np.array([[1, 0, 0],
[0, 1, 0],
[0, 0, 0],
[0, 0, 1]])
E = np.array([[0, 1, 2, 3]])
Ipos, Ineg, Izer = check(V, E)
assert len(Ipos) == 1
assert len(Ineg) == 0
assert len(Izer) == 0
# (neg)
E = np.array([[1, 0, 2, 3]])
Ipos, Ineg, Izer = check(V, E)
assert len(Ipos) == 0
assert len(Ineg) == 1
assert len(Izer) == 0
# 2D
# CCW (positive)
V = np.array([[1, 0, 0],
[0, 1, 0],
[0, 0, 0]])
E = np.array([[0, 1, 2]])
Ipos, Ineg, Izer = check(V, E)
assert len(Ipos) == 1
assert len(Ineg) == 0
assert len(Izer) == 0
# CW (negative)
E = np.array([[0, 2, 1]])
Ipos, Ineg, Izer = check(V, E)
assert len(Ipos) == 0
assert len(Ineg) == 1
assert len(Izer) == 0
mesh = meshio.read('testmesh.msh')
V = mesh.points
E = mesh.cells_dict['triangle']
Ipos, Ineg, Izer = check(V, E)
assert len(Ipos) == 1
assert len(Ineg) == 1
assert len(Izer) == 0
if __name__ == '__main__':
"""
python check.py mymesh.msh
"""
mname = sys.argv[1]
mesh = meshio.read(mname)
V = mesh.points
if 'tetra' in mesh.cells_dict:
E = mesh.cells_dict['tetra']
else:
E = mesh.cells_dict['triangle']
Ipos, Ineg, Izer = check(V, E)
print(f'found {len(Ipos)} positive elements')
print(f'found {len(Ineg)} negative elements')
print(f'found {len(Izer)} co-planar elements')