-
Notifications
You must be signed in to change notification settings - Fork 0
/
utilities.py
101 lines (80 loc) · 3.04 KB
/
utilities.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
94
95
96
97
98
99
100
101
import numpy as np
import scipy
import sys, os
sys.path.insert(0, os.getcwd()+"/../libigl/python/")
import pyigl as igl
import cPickle as pickle
def dist(p1, p2):
return np.linalg.norm(p1-p2)
def check_mode_symmetry(mode_vec, lev):
minNodeNum = Node.number - len(lev.nodes)
u_f = [[0 for x in range(len(X))] for y in range(len(Y))]
for n in lev.nodes:
p1 = np.array([ mode_vec[2*(n.id-minNodeNum)], mode_vec[2*(n.id-minNodeNum)+1], 0 ])
u_f[n.point[0]][n.point[1]] = np.linalg.norm(p1)
U = np.matrix(u_f)
print(np.allclose(U, U.T, atol=1e-2))
def is_invertible(H):
return H.shape[0]==H.shape[1] and np.linalg.matrix_rank(H) == H.shape[0]
def is_sym_pos_def(x):
# print(sorted(np.linalg.eigvals(x)))
pd = np.all(np.linalg.eigvals(x) >= -1e-8)
sym = np.allclose(x, x.T, atol=1e-10)
return pd and sym
def volume_of_tet(n1, n2, n3, n4):
x = np.array([n1[0] - n4[0], n2[0] - n4[0], n3[0] - n4[0]])
y = np.array([n1[1] - n4[1], n2[1] - n4[1], n3[1] - n4[1]])
z = np.array([n1[2] - n4[2], n2[2] - n4[2], n3[2] - n4[2]])
# print(n1, n2, n3, n4)
Dm = np.matrix([x, y, z])
v = (1.0/6)*np.linalg.det(Dm)
# print(np.linalg.det(Dm))
return v
def PointInTriangle(pt, v1, v2, v3):
#Code from stack overflow here:
#https://stackoverflow.com/questions/20248076/how-do-i-check-if-a-point-is-inside-a-triangle-on-the-line-is-ok-too
def sign(p1, p2, p3):
# print(p1, p2, p3)
return (p1[0] - p3[0]) * (p2[1] - p3[1]) - (p2[0] - p3[0]) * (p1[1] - p3[1])
def PointInAABB(pt, c1, c2):
return c2[0] <= pt[0] <= c1[0] and \
c2[1] <= pt[1] <= c1[1]
b1 = sign(pt, v1, v2) <= 0
b2 = sign(pt, v2, v3) <= 0
b3 = sign(pt, v3, v1) <= 0
return ((b1 == b2) and (b2 == b3)) and \
PointInAABB(pt, map(max, v1, v2, v3), map(min, v1, v2, v3))
def general_eig_solve(A, B):
#pass in A = K matrix, and B = M matrix
# eigvals, eigvecs = scipy.sparse.linalg.eigs(np.matrix(A),k = 4, M = np.matrix(B))
eigvals, eigvecs = scipy.linalg.eigh(np.matrix(A), np.matrix(B), overwrite_a = False, overwrite_b = False)
#eigvals, eigvecs = np.linalg.eig(np.linalg.inv(B)*A)
# print("Check general eigen problem solution")
# print(B)
# print((B.dot(eigvecs)).T.dot(eigvecs))
# print(eigvals)
return eigvals, eigvecs
def mesh_to_obj(name, V, F):
V3 = np.zeros((V.shape[0], V.shape[1]+1))
V3[:,:-1] = V
igl.writeOBJ(name, igl.eigen.MatrixXd(V3), igl.eigen.MatrixXi(F))
def serialize_layers(name, hMesh):
pfile = open(name, 'wb')
pickle.dump(hMesh, pfile, protocol=pickle.HIGHEST_PROTOCOL)
pfile.close()
return 1
def unserialize_layers(name):
pfile = open(name, 'rb')
layers = pickle.load(pfile)
pfile.close()
return layers
def serialize_mesh(name, mesh):
pfile = open(name, 'wb')
pickle.dump(mesh, pfile, protocol=pickle.HIGHEST_PROTOCOL)
pfile.close()
return 1
def unserialize_mesh(name):
pfile = open(name, 'rb')
mesh = pickle.load(pfile)
pfile.close()
return mesh