30
31
32
33#include "implicit_f.inc"
34
35
36
37#include "mvsiz_p.inc"
38
39
40
41
42
43
44 INTEGER NEL
45 INTEGER NC1(*), NC2(*), NC3(*), NC4(*),
46 . NC5(*), NC6(*), NC7(*), NC8(*),
47 . IXS(NIXS,*)
48
49 DOUBLE PRECISION
50 . QF(NEL,24),KS(24,24,NEL),
51 . V1(MVSIZ,3,3),V2(MVSIZ,3,3),V3(MVSIZ,3,3),V4(MVSIZ,3,3),
52 . V5(MVSIZ,3,3),V6(MVSIZ,3,3),V7(MVSIZ,3,3),V8(MVSIZ,3,3),
53 . X1(MVSIZ), Y1(), Z1(MVSIZ),
54 . X2(MVSIZ), Y2(MVSIZ), Z2(MVSIZ),
55 . X3(MVSIZ), Y3(MVSIZ), Z3(MVSIZ),
56 . X4(MVSIZ), Y4(MVSIZ), Z4(MVSIZ),
57 . X5(MVSIZ), Y5(MVSIZ), Z5(MVSIZ),
58 . X6(MVSIZ), Y6(MVSIZ), Z6(MVSIZ),
59 . X7(MVSIZ), Y7(MVSIZ), Z7(MVSIZ),
60 . X8(MVSIZ), Y8(MVSIZ), Z8(MVSIZ)
62 . x(3,*)
63
64
65
66 INTEGER I,II,INOD,J,J2,J3,J4,J5,J6,J7,J8,K
67 double precision
68 . q1,q2,q3, a(3,3), b(3,24),
69 . qx1,qy1,qz1,qx2,qy2,qz2,qx3,qy3,qz3,qx4,qy4,qz4,
70 . qx5,qy5,qz5,qx6,qy6,qz6,qx7,qy7,qz7,qx8,qy8,qz8,
71 . xi1, yi1, zi1, xi2, yi2, zi2, xi3, yi3, zi3,
72 . xi4, yi4, zi4, xi5, yi5, zi5, xi6, yi6, zi6,
73 . xi7, yi7, zi7, xi8, yi8, zi8
74
75
76 ks(1:24,1:24,1:nel) = zero
77
78 DO inod=1,8
79 ii = 3*(inod-1)
80 DO j=1,3
81 DO i=1,nel
82 q1 = qf(i,ii+1)
83 q2 = qf(i,ii+2)
84 q3 = qf(i,ii+3)
85 ks(ii+1,j,i)= q3*v1(i,2,j)-q2*v1(i,3,j)
86 ks(ii+2,j,i)=-q3*v1(i,1,j)+q1*v1(i,3,j)
87 ks(ii+3,j,i)= q2*v1(i,1,j)-q1*v1(i,2,j)
88 j2 = j+3
89 ks(ii+1,j2,i)= q3*v2(i,2,j)-q2*v2(i,3,j)
90 ks(ii+2,j2,i)=-q3*v2(i,1,j)+q1*v2(i,3,j)
91 ks(ii+3,j2,i)= q2*v2(i,1,j)-q1*v2(i,2,j)
92 j3 = j+6
93 ks(ii+1,j3,i)= q3*v3(i,2,j)-q2*v3(i,3,j)
94 ks(ii+2,j3,i)=-q3*v3(i,1,j)+q1*v3(i,3,j)
95 ks(ii+3,j3,i)= q2*v3(i,1,j)-q1*v3(i,2,j)
96 j4 = j+9
97 ks(ii+1,j4,i)= q3*v4(i,2,j)-q2*v4(i,3,j)
98 ks(ii+2,j4,i)=-q3*v4(i,1,j)+q1*v4(i,3,j)
99 ks(ii+3,j4,i)= q2*v4(i,1,j)-q1*v4(i,2,j)
100 j5 = j+12
101 ks(ii+1,j5,i)= q3*v5(i,2,j)-q2*v5(i,3,j)
102 ks(ii+2,j5,i)=-q3*v5(i,1,j)+q1*v5(i,3,j)
103 ks(ii+3,j5,i)= q2*v5(i,1,j)-q1*v5(i,2,j)
104 j6 = j+15
105 ks(ii+1,j6,i)= q3*v6(i,2,j)-q2*v6(i,3,j)
106 ks(ii+2,j6,i)=-q3*v6(i,1,j)+q1*v6(i,3,j)
107 ks(ii+3,j6,i)= q2*v6(i,1,j)-q1*v6(i,2,j)
108 j7 = j+18
109 ks(ii+1,j7,i)= q3*v7(i,2,j)-q2*v7(i,3,j)
110 ks(ii+2,j7,i)=-q3*v7(i,1,j)+q1*v7(i,3,j)
111 ks(ii+3,j7,i)= q2*v7(i,1,j)-q1*v7(i,2,j)
112 j8 = j+21
113 ks(ii+1,j8,i)= q3*v8(i,2,j)-q2*v8(i,3,j)
114 ks(ii+2,j8,i)=-q3*v8(i,1,j)+q1*v8(i,3,j)
115 ks(ii+3,j8,i)= q2*v8(i,1,j)-q1*v8(i,2,j)
116 ENDDO
117 ENDDO
118 ENDDO
119
120 DO inod=1,8
121 ii = 3*(inod-1)
122 DO j=1,3
123 DO i=1,nel
124 q1 = qf(i,ii+1)
125 q2 = qf(i,ii+2)
126 q3 = qf(i,ii+3)
127 ks(j,ii+1,i)=ks(j,ii+1,i)+q3*v1(i,2,j)-q2*v1(i,3,j)
128 ks(j,ii+2,i)=ks(j,ii+2,i)-q3*v1(i,1,j)+q1*v1(i,3,j)
129 ks(j,ii+3,i)=ks(j,ii+3,i)+q2*v1(i,1,j)-q1*v1(i,2,j)
130 j2=j+3
131 ks(j2,ii+1,i)=ks(j2,ii+1,i)+q3*v2(i,2,j)-q2*v2(i,3,j)
132 ks(j2,ii+2,i)=ks(j2,ii+2,i)-q3*v2(i,1,j)+q1*v2(i,3,j)
133 ks(j2,ii+3,i)=ks(j2,ii+3,i)+q2*v2(i,1,j)-q1*v2(i,2,j)
134 j3=j+6
135 ks(j3,ii+1,i)=ks(j3,ii+1,i)+q3*v3(i,2,j)-q2*v3(i,3,j)
136 ks(j3,ii+2,i)=ks(j3,ii+2,i)-q3*v3(i,1,j)+q1*v3(i,3,j)
137 ks(j3,ii+3,i)=ks(j3,ii+3,i)+q2*v3(i,1,j)-q1*v3(i,2,j)
138 j4=j+9
139 ks(j4,ii+1,i)=ks(j4,ii+1,i)+q3*v4(i,2,j)-q2*v4(i,3,j)
140 ks(j4,ii+2,i)=ks(j4,ii+2,i)-q3*v4(i,1,j)+q1*v4(i,3,j)
141 ks(j4,ii+3,i)=ks(j4,ii+3,i)+q2*v4(i,1,j)-q1*v4(i,2,j)
142 j5=j+12
143 ks(j5,ii+1,i)=ks(j5,ii+1,i)+q3*v5(i,2,j)-q2*v5(i,3,j)
144 ks(j5,ii+2,i)=ks(j5,ii+2,i)-q3*v5(i,1,j)+q1*v5(i,3,j)
145 ks(j5,ii+3,i)=ks(j5,ii+3,i)+q2*v5(i,1,j)-q1*v5(i,2,j)
146 j6=j+15
147 ks(j6,ii+1,i)=ks(j6,ii+1,i)+q3*v6(i,2,j)-q2*v6(i,3,j)
148 ks(j6,ii+2,i)=ks(j6,ii+2,i)-q3*v6(i,1,j)+q1*v6(i,3,j)
149 ks(j6,ii+3,i)=ks(j6,ii+3,i)+q2*v6(i,1,j)-q1*v6(i,2,j)
150 j7=j+18
151 ks(j7,ii+1,i)=ks(j7,ii+1,i)+q3*v7(i,2,j)-q2*v7(i,3,j)
152 ks(j7,ii+2,i)=ks(j7,ii+2,i)-q3*v7(i,1,j)+q1*v7(i,3,j)
153 ks(j7,ii+3,i)=ks(j7,ii+3,i)+q2*v7(i,1,j)-q1*v7(i,2,j)
154 j8=j+21
155 ks(j8,ii+1,i)=ks(j8,ii+1,i)+q3*v8(i,2,j)-q2*v8(i,3,j)
156 ks(j8,ii+2,i)=ks(j8,ii+2,i)-q3*v8(i,1,j)+q1*v8(i,3,j)
157 ks(j8,ii+3,i)=ks(j8,ii+3,i)+q2*v8(i,1,j)-q1*v8(i,2,j)
158 ENDDO
159 ENDDO
160 ENDDO
161
162 DO i=1,nel
163 nc1(i)=ixs(2,i)
164 nc2(i)=ixs(3,i)
165 nc3(i)=ixs(4,i)
166 nc4(i)=ixs(5,i)
167 nc5(i)=ixs(6,i)
168 nc6(i)=ixs(7,i)
169 nc7(i)=ixs(8,i)
170 nc8(i)=ixs(9,i)
171 ENDDO
172 DO i=1,nel
173 x1(i)=x(1,nc1(i))
174 y1(i)=x(2,nc1(i))
175 z1(i)=x(3,nc1(i))
176 x2(i)=x(1,nc2(i))
177 y2(i)=x(2,nc2(i))
178 z2(i)=x(3,nc2(i))
179 x3(i)=x(1,nc3(i))
180 y3(i)=x(2,nc3(i))
181 z3(i)=x(3,nc3(i))
182 x4(i)=x(1,nc4(i))
183 y4(i)=x(2,nc4(i))
184 z4(i)=x(3,nc4(i))
185 x5(i)=x(1,nc5(i))
186 y5(i)=x(2,nc5(i))
187 z5(i)=x(3,nc5(i))
188 x6(i)=x(1,nc6(i))
189 y6(i)=x(2,nc6(i))
190 z6(i)=x(3,nc6(i))
191 x7(i)=x(1,nc7(i))
192 y7(i)=x(2,nc7(i))
193 z7(i)=x(3,nc7(i))
194 x8(i)=x(1,nc8(i))
195 y8(i)=x(2,nc8(i))
196 z8(i)=x(3,nc8(i))
197 ENDDO
198 DO i=1,nel
199 xi1 = zero
200 yi1 = zero
201 zi1 = zero
202 xi2 = x2(i)-x1(i)
203 yi2 = y2(i)-y1(i)
204 zi2 = z2(i)-z1(i)
205 xi3 = x3(i)-x1(i)
206 yi3 = y3(i)-y1(i)
207 zi3 = z3(i)-z1(i)
208 xi4 = x4(i)-x1(i)
209 yi4 = y4(i)-y1(i)
210 zi4 = z4(i)-z1(i)
211 xi5 = x5(i)-x1(i)
212 yi5 = y5(i)-y1(i)
213 zi5 = z5(i)-z1(i)
214 xi6 = x6(i)-x1(i)
215 yi6 = y6(i)-y1(i)
216 zi6 = z6(i)-z1(i)
217 xi7 = x7(i)-x1(i)
218 yi7 = y7(i)-y1(i)
219 zi7 = z7(i)-z1(i)
220 xi8 = x8(i)-x1(i)
221 yi8 = y8(i)-y1(i)
222 zi8 = z8(i)-z1(i)
223
224 qx1 = qf(i,1)
225 qy1 = qf(i,2)
226 qz1 = qf(i,3)
227 qx2 = qf(i,4)
228 qy2 = qf(i,5)
229 qz2 = qf(i,6)
230 qx3 = qf(i,7)
231 qy3 = qf(i,8)
232 qz3 = qf(i,9)
233 qx4 = qf(i,10)
234 qy4 = qf(i,11)
235 qz4 = qf(i,12)
236 qx5 = qf(i,13)
237 qy5 = qf(i,14)
238 qz5 = qf(i,15)
239 qx6 = qf(i,16)
240 qy6 = qf(i,17)
241 qz6 = qf(i,18)
242 qx7 = qf(i,19)
243 qy7 = qf(i,20)
244 qz7 = qf(i,21)
245 qx8 = qf(i,22)
246 qy8 = qf(i,23)
247 qz8 = qf(i,24)
248
249 a(1,1) = -zi1*qz1+yi1*qy1-zi2*qz2+yi2*qy2
250 . -zi3*qz3+yi3*qy3-zi4*qz4+yi4*qy4
251 . -zi5*qz5+yi5*qy5-zi6*qz6+yi6*qy6
252 . -zi7*qz7+yi7*qy7-zi8*qz8+yi8*qy8
253 a(2,1) = -xi1*qy1-xi2*qy2-xi3*qy3-xi4*qy4
254 . -xi5*qy5-xi6*qy6-xi7*qy7-xi8*qy8
255 a(3,1) = +xi1*qz1+xi2*qz2+xi3*qz3+xi4*qz4
256 . +xi5*qz5+xi6*qz6+xi7*qz7+xi8*qz8
257 a(1,2) = yi1*qx1+yi2*qx2+yi3*qx3+yi4*qx4
258 . +yi5*qx5+yi6*qx6+yi7*qx7+yi8*qx8
259 a(2,2) = -zi1*qz1-xi1*qx1-zi2*qz2-xi2*qx2
260 . -zi3*qz3-xi3*qx3-zi4*qz4-xi4*qx4
261 . -zi5*qz5-xi5*qx5-zi6*qz6-xi6*qx6
262 . -zi7*qz7-xi7*qx7-zi8*qz8-xi8*qx8
263 a(3,2) = yi1*qz1+yi2*qz2+yi3*qz3+yi4*qz4
264 . +yi5*qz5+yi6*qz6+yi7*qz7+yi8*qz8
265 a(1,3) = zi1*qx1+zi2*qx2+zi3*qx3+zi4*qx4
266 . +zi5*qx5+zi6*qx6+zi7*qx7+zi8*qx8
267 a(2,3) = zi1*qy1+zi2*qy2+zi3*qy3+zi4*qy4
268 . +zi5*qy5+zi6*qy6+zi7*qy7+zi8*qy8
269 a(3,3) = -yi1*qy1-xi1*qx1-yi2*qy2-xi2*qx2
270 . -yi3*qy3-xi3*qx3-yi4*qy4-xi4*qx4
271 . -yi5*qy5-xi5*qx5-yi6*qy6-xi6*qx6
272 . -yi7*qy7-xi7*qx7-yi8*qy8-xi8*qx8
273
274 a(1,2) = half*(a(1,2)+a(2,1))
275 a(1,3) = half*(a(1,3)+a(3,1))
276 a(2,3) = half*(a(2,3)+a(3,2))
277 a(2,1) = a(1,2)
278 a(3,1) = a(1,3)
279 a(3,2) = a(2,3)
280
281 DO j=1,3
282 b(1,j) = a(1,1)*v1(i,1,j)+a(1,2)*v1(i,2,j)+a(1,3)*v1(i,3,j)
283 b(2,j) = a(2,1)*v1(i,1,j)+a(2,2)*v1(i,2,j)+a(2,3)*v1(i,3,j)
284 b(3,j) = a(3,1)*v1(i,1,j)+a(3,2)*v1(i,2,j)+a(3,3)*v1(i,3,j)
285 j2 = j+3
286 b(1,j2) = a(1,1)*v2(i,1,j)+a(1,2)*v2(i,2,j)+a(1,3)*v2(i,3,j)
287 b(2,j2) = a(2,1)*v2(i,1,j)+a(2,2)*v2(i,2,j)+a(2,3)*v2(i,3,j)
288 b(3,j2) = a(3,1)*v2(i,1,j)+a(3,2)*v2(i,2,j)+a(3,3)*v2(i,3,j)
289 j3 = j+6
290 b(1,j3) = a(1,1)*v3(i,1,j)+a(1,2)*v3(i,2,j)+a(1,3)*v3(i,3,j)
291 b(2,j3) = a(2,1)*v3(i,1,j)+a(2,2)*v3(i,2,j)+a(2,3)*v3(i,3,j)
292 b(3,j3) = a(3,1)*v3(i,1,j)+a(3,2)*v3(i,2,j)+a(3,3)*v3(i,3,j)
293 j4 = j+9
294 b(1,j4) = a(1,1)*v4(i,1,j)+a(1,2)*v4(i,2,j)+a(1,3)*v4(i,3,j)
295 b(2,j4) = a(2,1)*v4(i,1,j)+a(2,2)*v4(i,2,j)+a(2,3)*v4(i,3,j)
296 b(3,j4) = a(3,1)*v4(i,1,j)+a(3,2)*v4(i,2,j)+a(3,3)*v4(i,3,j)
297 j5 = j+12
298 b(1,j5) = a(1,1)*v5(i,1,j)+a(1,2)*v5(i,2,j)+a(1,3)*v5(i,3,j)
299 b(2,j5) = a(2,1)*v5(i,1,j)+a(2,2)*v5(i,2,j)+a(2,3)*v5(i,3,j)
300 b(3,j5) = a(3,1)*v5(i,1,j)+a(3,2)*v5(i,2,j)+a(3,3)*v5(i,3,j)
301 j6 = j+15
302 b(1,j6) = a(1,1)*v6(i,1,j)+a(1,2)*v6(i,2,j)+a(1,3)*v6(i,3,j)
303 b(2,j6) = a(2,1)*v6(i,1,j)+a(2,2)*v6(i,2,j)+a(2,3)*v6(i,3,j)
304 b(3,j6) = a(3,1)*v6(i,1,j)+a(3,2)*v6(i,2,j)+a(3,3)*v6(i,3,j)
305 j7 = j+18
306 b(1,j7) = a(1,1)*v7(i,1,j)+a(1,2)*v7(i,2,j)+a(1,3)*v7(i,3,j)
307 b(2,j7) = a(2,1)*v7(i,1,j)+a(2,2)*v7(i,2,j)+a(2,3)*v7(i,3,j)
308 b(3,j7) = a(3,1)*v7(i,1,j)+a(3,2)*v7(i,2,j)+a(3,3)*v7(i,3,j)
309 j8 = j+21
310 b(1,j8) = a(1,1)*v8(i,1,j)+a(1,2)*v8(i,2,j)+a(1,3)*v8(i,3,j)
311 b(2,j8) = a(2,1)*v8(i,1,j)+a(2,2)*v8(i,2,j)+a(2,3)*v8(i,3,j)
312 b(3,j8) = a(3,1)*v8(i,1,j)+a(3,2)*v8(i,2,j)+a(3,3)*v8(i,3,j)
313 ENDDO
314
315 DO k=1,24
316 ks(1,k,i) = ks(1,k,i)+v1(i,1,1)*b(1,k)+v1(i,2,1)*b(2,k)+v1(i,3,1)*b(3,k)
317 ks(2,k,i) = ks(2,k,i)+v1(i,1,2)*b(1,k)+v1(i,2,2)*b(2,k)+v1(i,3,2)*b(3,k
318 ks(3,k,i) = ks
319 ks(4,k,i) = ks(4,k,i)+v2(i,1,1)*b(1,k)+v2(i,2,1)*b(2,k)+v2(i,3,1)*b(3,k)
320 ks(5,k,i) = ks(5,k,i)+v2(i,1,2)*b(1,k)+v2(i,2,2)*b(2,k)+v2(i,3,2)*b(3,k)
321 ks(6,k,i) = ks(6,k,i)+v2(i,1,3)*b(1,k)+v2(i,2,3)*b(2,k)+v2(i,3,3)*b(3,k)
322 ks(7,k,i) = ks(7,k,i
323 ks(8,k,i) = ks(8,k,i)+v3(i,1,2)*b(1,k)+v3(i,2,2)*b(2,k)+v3(i,3,2)*b(3,k)
324 ks(9,k,i) = ks(9,k,i)+v3(i,1,3)*b(1,k)+v3(i,2,3)*b(2,k)+v3(i,3,3)*b(3,k)
325 ks(10,k,i) = ks(10,k,i)+v4(i,1,1)*b(1,k)+v4(i,2,1)*b(2,k)+v4(i,3,1)*b(3,k)
326 ks(11,k,i) = ks(11,k,i)+v4(i,1,2)*b(1,k)+v4(i,2,2)*b(2,k)+v4(i,3,2)*b(3,k)
327 ks(12,k,i) = ks(12,k,i)+v4(i,1,3)*b(1,k)+v4(i,2,3)*b(2,k)+v4(i,3,3)*b(3,k)
328 ks(13,k,i) = ks(13,k,i)+v5(i,1,1)*b(1,k)+v5(i,2,1)*b(2,k)+v5(i,3,1)*b(3,k)
329 ks(14,k,i) = ks(14,k,i)+v5(i,1,2)*b(1,k)+v5(i,2,2)*b(2,k)+v5(i,3,2)*b(3,k)
330 ks(15,k,i) = ks(15,k,i)+v5(i,1,3)*b(1,k)+v5(i,2,3)*b(2,k)+v5(i,3,3)*b(3,k)
331 ks(16,k,i) = ks(16,k,i)+v6(i,1,1)*b(1,k)+v6(i,2,1)*b(2,k)+v6(i,3,1)*b(3,k)
332 ks(17,k,i) = ks(17,k,i)+v6(i,1,2)*b(1,k)+v6(i,2,2)*b(2,k)+v6(i,3,2)*b(3,k)
333 ks(18,k,i) = ks(18,k,i)+v6(i,1,3)*b(1,k)+v6(i,2,3)*b(2,k)+v6(i,3,3)*b(3,k)
334 ks(19,k,i) = ks(19,k,i)+v7(i,1,1)*b(1,k)+v7(i,2,1)*b(2,k)+v7(i,3,1)*b(3,k)
335 ks(20,k,i) = ks(20,k,i)+v7(i,1,2)*b(1,k)+v7(i,2,2)*b(2,k)+v7(i,3,2)*b(3,k)
336 ks(21,k,i) = ks(21,k,i)+v7(i,1,3)*b(1,k)+v7(i,2,3)*b(2,k)+v7(i,3,3)*b(3,k)
337 ks(22,k,i) = ks(22,k,i)+v8(i,1,1)*b(1,k)+v8(i,2,1)*b(2,k)+v8(i,3,1)*b(3,k)
338 ks(23,k,i) = ks(23,k,i)+v8(i,1,2)*b(1,k)+v8(i,2,2)*b(2,k)+v8(i,3,2)*b(3,k)
339 ks(24,k,i) = ks(24,k,i)+v8(i,1,3)*b(1,k)+v8(i,2,3)*b(2,k)+v8(i,3,3)*b(3,k)
340 ENDDO
341 ENDDO
342
343 RETURN