1 |
vert0 = {x0, y0, z0} |
2 |
vert1 = {x1, y1, z1} |
3 |
vert2 = {x2, y2, z2} |
4 |
|
5 |
centroid = (vert0 + vert1 + vert2) / 3 |
6 |
a = vert0 - vert1 |
7 |
b = vert0 - vert2 |
8 |
c = vert1 - vert2 |
9 |
|
10 |
u0 = -a |
11 |
v0 = centroid - vert0 |
12 |
s0 = Sqrt[Cross[u0,v0] . Cross[u0,v0]] / 2 |
13 |
|
14 |
u1 = -c |
15 |
v1 = centroid - vert1 |
16 |
s1 = Sqrt[Cross[u1,v1] . Cross[u1,v1]] / 2 |
17 |
|
18 |
u2 = b |
19 |
v2 = centroid - vert2 |
20 |
s2 = Sqrt[Cross[u2,v2] . Cross[u2,v2]] / 2 |
21 |
|
22 |
|
23 |
sc0 = (centroid + vert1 + vert0) / 3 |
24 |
dr0 = centroid - sc0 |
25 |
l0 = 1 / (dr0.dr0) |
26 |
|
27 |
G0 = (IdentityMatrix[3] + Outer[Times, dr0, dr0]*l0)*Sqrt[l0]*s0*8/(Pi*eta) |
28 |
|
29 |
sc1 = (centroid + vert1 + vert2) / 3 |
30 |
dr1 = centroid - sc1 |
31 |
l1 = 1 / (dr1.dr1) |
32 |
|
33 |
G1 = (IdentityMatrix[3] + Outer[Times, dr1, dr1]*l1)*Sqrt[l1]*s1*8/(Pi*eta) |
34 |
|
35 |
sc2 = (centroid + vert2 + vert0) / 3 |
36 |
dr2 = centroid - sc2 |
37 |
l2 = 1 / (dr2.dr2) |
38 |
|
39 |
G2 = (IdentityMatrix[3] + Outer[Times, dr2, dr2]*l2)*Sqrt[l2]*s2*8/(Pi*eta) |
40 |
|
41 |
|
42 |
Hmat = G0 + G1 + G2 |
43 |
|
44 |
Xi = Inverse[Hmat] |