-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathFE_template.py
180 lines (165 loc) · 9.8 KB
/
FE_template.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
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
import numpy as np
import torch
def to_torch(x):
return torch.tensor(x).float()
def computeFluidMatrices(dx,dy,mu0):
t2 = dx**2;
t3 = dy**2;
t4 = 1.0/dx;
t5 = 1.0/dy;
t10 = mu0/3.0;
t11 = mu0/4.0;
t19 = (dx*dy)/4.0;
t20 = dx/9.0;
t21 = dx*(2.0/9.0);
t22 = dx/1.8e+1;
t23 = dx*(5.0/1.8e+1);
t24 = dx/3.6e+1;
t25 = dx*(5.0/3.6e+1);
t26 = dy/9.0;
t27 = dy*(2.0/9.0);
t28 = dy/1.8e+1;
t29 = dy*(5.0/1.8e+1);
t30 = dy/3.6e+1;
t31 = dy*(5.0/3.6e+1);
t34 = mu0/9.0;
t35 = mu0/1.2e+1;
t36 = mu0*(4.0/9.0);
t37 = mu0/3.6e+1;
t56 = (dx*dy)/2.25e+2;
t57 = dx*dy*(2.0/2.25e+2);
t58 = dx*dy*(4.0/2.25e+2);
t59 = (dx*dy)/4.5e+2;
t62 = dx*dy*(8.0/2.25e+2);
t63 = dx*dy*(1.6e+1/2.25e+2);
t65 = (dx*dy)/9.0e+2;
t68 = dx*dy*(6.4e+1/2.25e+2);
t6 = t2*2.0;
t7 = t2*7.0;
t8 = t3*2.0;
t9 = t3*7.0;
t12 = t2*8.0;
t13 = t2*3.2e+1;
t14 = -t3;
t17 = t3*8.0;
t18 = t3*3.2e+1;
t32 = -t10;
t33 = -t11;
t40 = -t20;
t41 = -t21;
t42 = -t22;
t43 = -t23;
t44 = -t24;
t45 = -t25;
t46 = -t26;
t47 = -t27;
t48 = -t28;
t49 = -t29;
t50 = -t30;
t51 = -t31;
t52 = -t34;
t53 = -t36;
t54 = -t35;
t55 = -t37;
t60 = -t56;
t61 = -t58;
t64 = -t59;
t15 = -t8;
t16 = -t9;
t38 = -t17;
t39 = -t18;
t66 = t3+t6;
t67 = t2+t8;
t70 = t3+t12;
t71 = t2+t17;
t73 = t6+t9;
t74 = t7+t8;
t75 = t6+t14;
t78 = t7+t17;
t79 = t9+t12;
t80 = t12+t14;
t108 = mu0*t4*t5*(t3-t6)*(-8.0/4.5e+1);
t110 = mu0*t4*t5*(t9-t12)*(-2.0/4.5e+1);
t113 = mu0*t4*t5*(t3-t6)*(-3.2e+1/4.5e+1);
t115 = mu0*t4*t5*(t9-t13)*(-1.0/4.5e+1);
t116 = mu0*t4*t5*(t3-t6)*(8.0/4.5e+1);
t117 = mu0*t4*t5*(t9-t12)*(2.0/4.5e+1);
t118 = mu0*t4*t5*(t3-t12)*(-1.6e+1/4.5e+1);
t119 = mu0*t4*t5*(t3-t6)*(3.2e+1/4.5e+1);
t120 = (mu0*t4*t5*(t9-t13))/4.5e+1;
t122 = mu0*t4*t5*(t9-t12)*(-1.0/9.0e+1);
t123 = mu0*t4*t5*(t3-t12)*(1.6e+1/4.5e+1);
t69 = t2+t15;
t72 = t2+t38;
t76 = t6+t16;
t77 = t7+t15;
t81 = t7+t38;
t82 = t12+t16;
t83 = t7+t39;
t84 = t13+t16;
t85 = mu0*t4*t5*t66*(2.0/4.5e+1);
t86 = mu0*t4*t5*t67*(2.0/4.5e+1);
t87 = (mu0*t4*t5*t70)/4.5e+1;
t88 = (mu0*t4*t5*t71)/4.5e+1;
t89 = mu0*t4*t5*t66*(8.0/4.5e+1);
t90 = mu0*t4*t5*t67*(8.0/4.5e+1);
t91 = mu0*t4*t5*t66*(1.4e+1/4.5e+1);
t92 = mu0*t4*t5*t67*(1.4e+1/4.5e+1);
t93 = (mu0*t4*t5*t66)/9.0e+1;
t94 = (mu0*t4*t5*t67)/9.0e+1;
t105 = mu0*t4*t5*t73*(1.6e+1/4.5e+1);
t106 = mu0*t4*t5*t74*(1.6e+1/4.5e+1);
t111 = mu0*t4*t5*t78*(8.0/4.5e+1);
t112 = mu0*t4*t5*t79*(8.0/4.5e+1);
t95 = -t89;
t96 = mu0*t4*t5*t69*(8.0/4.5e+1);
t97 = -t90;
t98 = mu0*t4*t5*t69*(3.2e+1/4.5e+1);
t99 = -t93;
t100 = -t94;
t102 = mu0*t4*t5*t72*(1.6e+1/4.5e+1);
t103 = (mu0*t4*t5*t76)/4.5e+1;
t104 = (mu0*t4*t5*t77)/4.5e+1;
t109 = mu0*t4*t5*t81*(2.0/4.5e+1);
t114 = (mu0*t4*t5*t83)/4.5e+1;
t121 = (mu0*t4*t5*t81)/9.0e+1;
t101 = -t96;
t107 = -t104;
t124 = -t121;
mt1 = np.array([t92,t11,t124,t35,t100,t55,t103,t54,t114,t32,t85,t34,t88,t34,t117,t10,t97,t53,t11,t91,t54,t107,t55,t99,t35,t122,t10,t109,t34,t87,t34,t86,t32,t120,t53,t95,t124,t54,t92,t33,t103,t35,t100,t37,t114,t10,t117,t32,t88,t52,t85,t52,t97,t36,t35,t107,t33,t91,t54,t122,t37,t99,t32,t109,t10,t120,t52,t86,t52,t87,t36,t95,t100,t55,t103,t54,t92,t11,t124,t35,t88,t34,t117,t10,t114,t32,t85,t34,t97,t53,t55,t99,t35,t122,t11,t91,t54,t107,t34,t86,t32,t120,t10,t109,t34,t87,t53,t95,t103,t35,t100,t37,t124,t54,t92,t33,t88,t52,t85,t52,t114,t10,t117,t32,t97,t36,t54,t122,t37,t99,t35,t107,t33,t91,t52,t86,t52,t87,t32,t109,t10,t120,t36,t95,t114,t10,t114,t32,t88,t34,t88,t52,t111,0.0,t97,t53,t96,0.0,t97,t36,t119,0.0,t32,t109,t10,t109,t34,t86,t52,t86,0.0,t106,t53,t95,0.0,t108,t36,t95,0.0,t123,t85,t34,t117,t10,t117,t32,t85,t52,t97,t53,t105,0.0,t97,t36,t101,0.0,t102,0.0,t34,t87]);
mt2 = np.array([t32,t120,t10,t120,t52,t87,t53,t95,0.0,t112,t36,t95,0.0,t116,0.0,t98,t88,t34,t88,t52,t114,t10,t114,t32,t96,0.0,t97,t36,t111,0.0,t97,t53,t119,0.0,t34,t86,t52,t86,t32,t109,t10,t109,0.0,t108,t36,t95,0.0,t106,t53,t95,0.0,t123,t117,t32,t85,t52,t85,t34,t117,t10,t97,t36,t101,0.0,t97,t53,t105,0.0,t102,0.0,t10,t120,t52,t87,t34,t87,t32,t120,t36,t95,0.0,t116,t53,t95,0.0,t112,0.0,t98,t97,t53,t97,t36,t97,t53,t97,t36,t119,0.0,t102,0.0,t119,0.0,t102,0.0,mu0*t4*t5*t67*(1.28e+2/4.5e+1),0.0,t53,t95,t36,t95,t53,t95,t36,t95,0.0,t123,0.0,t98,0.0,t123,0.0,t98,0.0,mu0*t4*t5*t66*(1.28e+2/4.5e+1)]);
Amu = np.concatenate([mt1,mt2]).reshape((18,18), order='F');
mt3 = np.array([t58,0.0,t60,0.0,t65,0.0,t60,0.0,t57,0.0,t64,0.0,t64,0.0,t57,0.0,t56,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,t60,0.0,t58,0.0,t60,0.0,t65,0.0,t57,0.0,t57,0.0,t64,0.0,t64,0.0,t56,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,t65,0.0,t60,0.0,t58,0.0,t60,0.0,t64,0.0,t57,0.0,t57,0.0,t64,0.0,t56,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,t60,0.0,t65,0.0,t60,0.0,t58,0.0,t64,0.0,t64,0.0,t57,0.0,t57,0.0,t56,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,t57,0.0,t57,0.0,t64,0.0,t64,0.0,t63,0.0,t56,0.0,t61,0.0,t56,0.0,t62,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,t64,0.0,t57,0.0,t57,0.0,t64,0.0,t56,0.0,t63,0.0,t56,0.0,t61,0.0,t62,0.0,0.0,0.0]);
mt4 = np.array([0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,t64,0.0,t64,0.0,t57,0.0,t57,0.0,t61,0.0,t56,0.0,t63,0.0,t56,0.0,t62,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,t57,0.0,t64,0.0,t64,0.0,t57,0.0,t56,0.0,t61,0.0,t56,0.0,t63,0.0,t62,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,t56,0.0,t56,0.0,t56,0.0,t56,0.0,t62,0.0,t62,0.0,t62,0.0,t62,0.0,t68,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0]);
AalphaX = np.concatenate([mt3,mt4]).reshape((18,18), order = 'F');
mt5 = np.array([0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,t58,0.0,t60,0.0,t65,0.0,t60,0.0,t57,0.0,t64,0.0,t64,0.0,t57,0.0,t56,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,t60,0.0,t58,0.0,t60,0.0,t65,0.0,t57,0.0,t57,0.0,t64,0.0,t64,0.0,t56,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,t65,0.0,t60,0.0,t58,0.0,t60,0.0,t64,0.0,t57,0.0,t57,0.0,t64,0.0,t56,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,t60,0.0,t65,0.0,t60,0.0,t58,0.0,t64,0.0,t64,0.0,t57,0.0,t57,0.0,t56,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,t57,0.0,t57,0.0,t64,0.0,t64,0.0,t63,0.0,t56,0.0,t61,0.0,t56,0.0,t62,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,t64]);
mt6 = np.array([0.0,t57,0.0,t57,0.0,t64,0.0,t56,0.0,t63,0.0,t56,0.0,t61,0.0,t62,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,t64,0.0,t64,0.0,t57,0.0,t57,0.0,t61,0.0,t56,0.0,t63,0.0,t56,0.0,t62,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,t57,0.0,t64,0.0,t64,0.0,t57,0.0,t56,0.0,t61,0.0,t56,0.0,t63,0.0,t62,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,t56,0.0,t56,0.0,t56,0.0,t56,0.0,t62,0.0,t62,0.0,t62,0.0,t62,0.0,t68]);
AalphaY = np.concatenate([mt5,mt6]).reshape((18,18), order = 'F');
mt7 = np.array([0.0,t58,0.0,t60,0.0,t65,0.0,t60,0.0,t57,0.0,t64,0.0,t64,0.0,t57,0.0,t56,t58,0.0,t60,0.0,t65,0.0,t60,0.0,t57,0.0,t64,0.0,t64,0.0,t57,0.0,t56,0.0,0.0,t60,0.0,t58,0.0,t60,0.0,t65,0.0,t57,0.0,t57,0.0,t64,0.0,t64,0.0,t56,t60,0.0,t58,0.0,t60,0.0,t65,0.0,t57,0.0,t57,0.0,t64,0.0,t64,0.0,t56,0.0,0.0,t65,0.0,t60,0.0,t58,0.0,t60,0.0,t64,0.0,t57,0.0,t57,0.0,t64,0.0,t56,t65,0.0,t60,0.0,t58,0.0,t60,0.0,t64,0.0,t57,0.0,t57,0.0,t64,0.0,t56,0.0,0.0,t60,0.0,t65,0.0,t60,0.0,t58,0.0,t64,0.0,t64,0.0,t57,0.0,t57,0.0,t56,t60,0.0,t65,0.0,t60,0.0,t58,0.0,t64,0.0,t64,0.0,t57,0.0,t57,0.0,t56,0.0,0.0,t57,0.0,t57,0.0,t64,0.0,t64,0.0,t63,0.0,t56,0.0,t61,0.0,t56,0.0,t62,t57,0.0,t57,0.0,t64,0.0,t64,0.0,t63,0.0,t56,0.0,t61,0.0,t56,0.0,t62,0.0,0.0,t64,0.0,t57,0.0,t57,0.0,t64,0.0,t56,0.0,t63,0.0,t56,0.0,t61,0.0,t62,t64,0.0]);
mt8 = np.array([t57,0.0,t57,0.0,t64,0.0,t56,0.0,t63,0.0,t56,0.0,t61,0.0,t62,0.0,0.0,t64,0.0,t64,0.0,t57,0.0,t57,0.0,t61,0.0,t56,0.0,t63,0.0,t56,0.0,t62,t64,0.0,t64,0.0,t57,0.0,t57,0.0,t61,0.0,t56,0.0,t63,0.0,t56,0.0,t62,0.0,0.0,t57,0.0,t64,0.0,t64,0.0,t57,0.0,t56,0.0,t61,0.0,t56,0.0,t63,0.0,t62,t57,0.0,t64,0.0,t64,0.0,t57,0.0,t56,0.0,t61,0.0,t56,0.0,t63,0.0,t62,0.0,0.0,t56,0.0,t56,0.0,t56,0.0,t56,0.0,t62,0.0,t62,0.0,t62,0.0,t62,0.0,t68,t56,0.0,t56,0.0,t56,0.0,t56,0.0,t62,0.0,t62,0.0,t62,0.0,t62,0.0,t68,0.0]);
AalphaXY = np.concatenate([mt7,mt8]).reshape((18,18), order = 'F');
z1 = np.zeros((5,5))
z2 = np.zeros((18,5))
Amu = np.vstack((np.hstack((Amu, z2)),\
np.hstack((z2.T,z1))))
AalphaX = np.vstack((np.hstack((AalphaX, z2)),\
np.hstack((z2.T,z1))))
AalphaXY = np.vstack((np.hstack((AalphaXY, z2)),\
np.hstack((z2.T,z1))))
AalphaY = np.vstack((np.hstack((AalphaY, z2)),\
np.hstack((z2.T,z1))))
B = np.array([t51,t45,t30,0.0,0.0,0.0,0.0,t24,t26,t43,t28,0.0,0.0,t22,t49,t20,t27,t21,t50,0.0,t31,t45,0.0,t24,0.0,0.0,t46,t43,t29,t20,0.0,t22,t48,0.0,t47,t21,0.0,0.0,0.0,t44,t31,t25,t50,0.0,0.0,t42,t29,t40,t46,t23,t48,0.0,t47,t41,0.0,t44,0.0,0.0,t30,0.0,t51,t25,0.0,t42,t28,0.0,t26,t23,t49,t40,t27,t41]).reshape((18,4), order = 'F');
z3 = np.zeros((18,18))
z4 = np.zeros((18,1))
z5 = np.zeros((4,5))
v1 = np.hstack((z3, B, z4))
v2 = np.hstack((B.T,z5))
v3 = np.zeros((1,23))
B = np.vstack((v1,v2,v3))
matArea = np.array([t19,t19,t19,t19]).reshape((4,1));
v1 = np.zeros((18,23))
v2 = np.hstack((np.zeros((4,22)),matArea))
v3 = np.hstack((np.zeros((1,18)),matArea.T,np.zeros((1,1))))
matArea = np.vstack((v1,v2,v3))
return to_torch(Amu), to_torch(B), to_torch(matArea), to_torch(AalphaX), to_torch(AalphaY), to_torch(AalphaXY)
#Amu, Aalpha, B, matArea = computeFluidMatrices(dx = 0.5,dy = 0.5,mu0 = 1.)