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