47
48
49
50#include "implicit_f.inc"
51
52
53
54#include "mvsiz_p.inc"
55
56
57
58
59
60
61 INTEGER, INTENT(IN) :: NEL
63 . nu,
64 . vx1(*),vx2(*),vx3(*),vx4(*),vx5(*),vx6(*),vx7(*),vx8(*),
65 . vy1(*),vy2(*),vy3(*),vy4(*),vy5(*),vy6(*),vy7(*),vy8(*),
66 . vz1(*),vz2(*),vz3(*),vz4(*),vz5(*),vz6(*),vz7(*),vz8(*),
67 . f11(*),f21(*),f31(*),f12(*),f22(*),f32(*),
68 . f13(*),f23(*),f33(*),f14(*),f24(*),f34(*),
69 . f15(*),f25(*),f35(*),f16(*),f26(*),f36(*),
70 . f17(*),f27(*),f37(*),f18(*),f28(*),f38(*),
71 . px1h1(*), px1h2(*), px1h3(*), px1h4(*),
72 . px2h1(*), px2h2(*), px2h3(*), px2h4(*),
73 . px3h1(*), px3h2(*), px3h3(*), px3h4(*),
74 . px4h1(*), px4h2(*), px4h3(*), px4h4(*),
75 . jr_1(mvsiz),js_1(mvsiz),jt_1(mvsiz),fcl(mvsiz)
76
77
78
79 INTEGER I, MX, J,K,IET, MT
81 . h11(mvsiz), h22(mvsiz), h33(mvsiz),
82 . h12(mvsiz), h13(mvsiz), h23(mvsiz),
83 . hgx1(mvsiz), hgx2(mvsiz), hgx3(mvsiz), hgx4(mvsiz),
84 . hgy1(mvsiz), hgy2(mvsiz), hgy3(mvsiz), hgy4(mvsiz),
85 . hgz1(mvsiz), hgz2(mvsiz), hgz3(mvsiz), hgz4(mvsiz),
86 . vx3478, vx2358, vx1467, vx1256,
87 . vy3478, vy2358, vy1467, vy1256,
88 . vz3478, vz2358, vz1467, vz1256,
89 . vx17, vy17, vz17,
90 . vx28, vy28, vz28,
91 . vx35, vy35, vz35,
92 . vx46, vy46, vz46,
93 . jr0(mvsiz),js0(mvsiz),jt0(mvsiz),nfhour(mvsiz,3,4),
94 . fhourt(3,4),dt05,rho0,etmax,nu1,nu2,nu3,nu4,deint,
95 . e_r,e_s,e_t,fac,fac1,fac2,coefh,hq13p,hq13n,hq24p,hq24n,ff
96
97
98 nu1 =two/(one-nu)
99 nu2 =nu*nu1
100 nu3 =two_third*(one + nu)
101 nu4 =nu
102 DO i=1,nel
103 jr0(i) = one/
max(em20,jr_1(i))
104 js0(i) = one/
max(em20,js_1(i))
105 jt0(i) = one/
max(em20,jt_1(i))
106 h11(i) = js0(i)*jt0(i)*jr_1(i)
107 h22(i) = jr0(i)*jt0(i)*js_1(i)
108 h33(i) = jr0(i)*js0(i)*jt_1(i)
109 h12(i) = jt0(i)
110 h13(i) = js0(i)
111 h23(i) = jr0(i)
112 ENDDO
113 DO i=1,nel
114 vx3478=vx3(i)-vx4(i)-vx7(i)+vx8(i)
115 vx2358=vx2(i)-vx3(i)-vx5(i)+vx8(i)
116 vx1467=vx1(i)-vx4(i)-vx6(i)+vx7(i)
117 vx1256=vx1(i)-vx2(i)-vx5(i)+vx6(i)
118
119 vy3478=vy3(i)-vy4(i)-vy7(i)+vy8(i)
120 vy2358=vy2(i)-vy3(i)-vy5(i)+vy8(i)
121 vy1467=vy1(i)-vy4(i)-vy6(i)+vy7(i)
122 vy1256=vy1(i)-vy2(i)-vy5(i)+vy6(i)
123
124 vz3478=vz3(i)-vz4(i)-vz7(i)+vz8(i)
125 vz2358=vz2(i)-vz3(i)-vz5(i)+vz8(i)
126 vz1467=vz1(i)-vz4(i)-vz6(i)+vz7(i)
127 vz1256=vz1(i)-vz2(i)-vz5(i)+vz6(i)
128!
129 hgx3(i)=(vx1467-vx2358)*one_over_8
130 hgx1(i)=(vx1467+vx2358)*one_over_8
131 hgx2(i)=(vx1256-vx3478)*one_over_8
132 hgx4(i)=-(vx1256+vx3478)*one_over_8
133
134 hgy3(i)=(vy1467-vy2358)*one_over_8
135 hgy1(i)=(vy1467+vy2358)*one_over_8
136 hgy2(i)=(vy1256-vy3478)*one_over_8
137 hgy4(i)=-(vy1256+vy3478)*one_over_8
138
139 hgz3(i)=(vz1467-vz2358)*one_over_8
140 hgz1(i)=(vz1467+vz2358)*one_over_8
141 hgz2(i)=(vz1256-vz3478)*one_over_8
142 hgz4(i)=-(vz1256+vz3478)*one_over_8
143 ENDDO
144 DO i=1,nel
145 vx17=vx1(i)-vx7(i)
146 vx28=vx2(i)-vx8(i)
147 vx35=vx3(i)-vx5(i)
148 vx46=vx4(i)-vx6(i)
149 vy17=vy1(i)-vy7(i)
150 vy28=vy2(i)-vy8(i)
151 vy35=vy3(i)-vy5(i)
152 vy46=vy4(i)-vy6(i)
153 vz17=vz1(i)-vz7(i)
154 vz28=vz2(i)-vz8(i)
155 vz35=vz3(i)-vz5(i)
156 vz46=vz4(i)-vz6(i)
157
158
159 hgx1(i)= hgx1(i)
160 & -(px1h1(i)*vx17+px2h1(i)*vx28
161 & +px3h1(i)*vx35+px4h1(i)*vx46)
162 hgy1(i)= hgy1(i)
163 & -(px1h1(i)*vy17+px2h1(i)*vy28
164 & +px3h1(i)*vy35+px4h1(i)*vy46)
165 hgz1(i)= hgz1(i)
166 & -(px1h1(i)*vz17+px2h1(i)*vz28
167 & +px3h1(i)*vz35+px4h1(i)*vz46)
168
169
170
171 hgx2(i)= hgx2(i)
172 & -(px1h2(i)*vx17+px2h2(i)*vx28
173 & +px3h2(i)*vx35+px4h2(i)*vx46)
174 hgy2(i)= hgy2(i)
175 & -(px1h2(i)*vy17+px2h2(i)*vy28
176 & +px3h2(i)*vy35+px4h2(i)*vy46)
177 hgz2(i)= hgz2(i)
178 & -(px1h2(i)*vz17+px2h2(i)*vz28
179 & +px3h2(i)*vz35+px4h2(i)*vz46)
180
181
182 hgx3(i)= hgx3(i)
183 & -(px1h3(i)*vx17+px2h3(i)*vx28
184 & +px3h3(i)*vx35+px4h3(i)*vx46)
185 hgy3(i)= hgy3(i)
186 & -(px1h3(i)*vy17+px2h3(i)*vy28
187 & +px3h3(i)*vy35+px4h3(i)*vy46)
188 hgz3(i)= hgz3(i)
189 & -(px1h3(i)*vz17+px2h3(i)*vz28
190 & +px3h3(i)*vz35+px4h3(i)*vz46)
191
192
193
194 hgx4(i)= hgx4(i)
195 & -(px1h4(i)*vx17+px2h4(i)*vx28
196 & +px3h4(i)*vx35+px4h4(i)*vx46)
197 hgy4(i)= hgy4(i)
198 & -(px1h4(i)*vy17+px2h4(i)*vy28
199 & +px3h4(i)*vy35+px4h4(i)*vy46)
200 hgz4(i)= hgz4(i)
201 & -(px1h4(i)*vz17+px2h4(i)*vz28
202 & +px3h4(i)*vz35+px4h4(i)*vz46)
203 ENDDO
204
205
206 DO i=1,nel
207 fhourt(1,1) = fcl(i)*hgx1(i)
208 fhourt(1,2) = fcl(i)*hgx2(i)
209 fhourt(1,3) = fcl(i)*hgx3(i)
210 fhourt(1,4) = fcl(i)*hgx4(i)
211 fhourt(2,1) = fcl(i)*hgy1(i)
212 fhourt(2,2) = fcl(i)*hgy2(i)
213 fhourt(2,3) = fcl(i)*hgy3(i)
214 fhourt(2,4) = fcl(i)*hgy4(i)
215 fhourt(3,1) = fcl(i)*hgz1(i)
216 fhourt(3,2) = fcl(i)*hgz2(i)
217 fhourt(3,3) = fcl(i)*hgz3(i)
218 fhourt(3,4) = fcl(i)*hgz4(i)
219
220 nfhour(i,1,1) = (h22(i)+h33(i))*fhourt(1,1)
221 . +h12(i)*fhourt(2,2)+h13(i)*fhourt(3,3)
222 nfhour(i,2,2) = (h11(i)+h33(i))*fhourt(2,2)
223 . +h23(i)*fhourt(3,3)+h12(i)*fhourt(1,1)
224 nfhour(i,3,3) = (h11(i)+h22(i))*fhourt(3,3)
225 . +h13(i)*fhourt(1,1)+h23(i)*fhourt(2,2)
226 nfhour(i,1,2) = nu1*h11(i)*fhourt(1,2)+nu2*h12(i)*fhourt(2,1)
227 nfhour(i,1,3) = nu1*h11(i)*fhourt(1,3)+nu2*h13(i)*fhourt(3,1)
228 nfhour(i,2,1) = nu1*h22(i)*fhourt(2,1)+nu2*h12(i)*fhourt(1,2)
229 nfhour(i,3,1) = nu1*h33(i)*fhourt(3,1)+nu2*h13(i)*fhourt(1,3)
230 nfhour(i,2,3) = nu1*h22(i)*fhourt(2,3)+nu2*h23(i)*fhourt(3,2)
231 nfhour(i,3,2) = nu1*h33(i)*fhourt(3,2)+nu2*h23(i)*fhourt(2,3)
232 nfhour(i,1,4) = nu3*h11(i)*fhourt(1,4)
233 nfhour(i,2,4) = nu3*h22(i)*fhourt(2,4)
234 nfhour(i,3,4) = nu3*h33(i)*fhourt(3,4)
235 ENDDO
236
237 DO i=1,nel
238 hq13p = (nfhour(i,1,1)+nfhour(i,1,3))*one_over_8
239 hq13n = (nfhour(i,1,1)-nfhour(i,1,3))*one_over_8
240 hq24p = (nfhour(i,1,2)+nfhour(i,1,4))*one_over_8
241 hq24n = (nfhour(i,1,2)-nfhour(i,1,4))*one_over_8
242 ff =-px1h1(i)*nfhour(i,1,1)-px1h2(i)*nfhour(i,1,2)
243 . -px1h3(i)*nfhour(i,1,3)-px1h4(i)*nfhour(i,1,4)
244 f11(i) =f11(i)-(hq13p+hq24n+ff)
245 f17(i) =f17(i)-(hq13p+hq24p-ff)
246 ff =-px2h1(i)*nfhour(i,1,1)-px2h2(i)*nfhour(i,1,2)
247 . -px2h3(i)*nfhour(i,1,3)-px2h4(i)*nfhour(i,1,4)
248 f12(i) =f12(i)-(hq13n-hq24n+ff)
249 f18(i) =f18(i)-(hq13n-hq24p-ff)
250 ff =-px3h1(i)*nfhour(i,1,1)-px3h2(i)*nfhour(i,1,2)
251 . -px3h3(i)*nfhour(i,1,3)-px3h4(i)*nfhour(i,1,4)
252 f13(i) =f13(i)-(-hq13n-hq24p+ff)
253 f15(i) =f15(i)-(-hq13n-hq24n-ff)
254 ff =-px4h1(i)*nfhour(i,1,1)-px4h2(i)*nfhour(i,1,2)
255 . -px4h3(i)*nfhour(i,1,3)-px4h4(i)*nfhour(i,1,4)
256 f14(i) =f14(i)-(-hq13p+hq24p+ff)
257 f16(i) =f16(i)-(-hq13p+hq24n-ff)
258 ENDDO
259 DO i=1,nel
260 hq13p = (nfhour(i,2,1)+nfhour(i,2,3))*one_over_8
261 hq13n = (nfhour(i,2,1)-nfhour(i,2,3))*one_over_8
262 hq24p = (nfhour(i,2,2)+nfhour(i,2,4))*one_over_8
263 hq24n = (nfhour(i,2,2)-nfhour(i,2,4))*one_over_8
264 ff =-px1h1(i)*nfhour(i,2,1)-px1h2(i)*nfhour(i,2,2)
265 . -px1h3(i)*nfhour(i,2,3)-px1h4(i)*nfhour(i,2,4)
266 f21(i) =f21(i)-(hq13p+hq24n+ff)
267 f27(i) =f27(i)-(hq13p+hq24p-ff)
268 ff =-px2h1(i)*nfhour(i,2,1)-px2h2(i)*nfhour(i,2,2)
269 . -px2h3(i)*nfhour(i,2,3)-px2h4(i)*nfhour(i,2,4)
270 f22(i) =f22(i)-(hq13n-hq24n+ff)
271 f28(i) =f28(i)-(hq13n-hq24p-ff)
272 ff =-px3h1(i)*nfhour(i,2,1)-px3h2(i)*nfhour(i,2,2)
273 . -px3h3(i)*nfhour(i,2,3)-px3h4(i)*nfhour(i,2,4)
274 f23(i) =f23(i)-(-hq13n-hq24p+ff)
275 f25(i) =f25(i)-(-hq13n-hq24n-ff)
276 ff =-px4h1(i)*nfhour(i,2,1)-px4h2(i)*nfhour(i,2,2)
277 . -px4h3(i)*nfhour(i,2,3)-px4h4(i)*nfhour(i,2,4)
278 f24(i) =f24(i)-(-hq13p+hq24p+ff)
279 f26(i) =f26(i)-(-hq13p+hq24n-ff)
280 ENDDO
281 DO i=1,nel
282 hq13p = (nfhour(i,3,1)+nfhour(i,3,3))*one_over_8
283 hq13n = (nfhour(i,3,1)-nfhour(i,3,3))*one_over_8
284 hq24p = (nfhour(i,3,2)+nfhour(i,3,4))*one_over_8
285 hq24n = (nfhour(i,3,2)-nfhour(i,3,4))*one_over_8
286 ff =-px1h1(i)*nfhour(i,3,1)-px1h2(i)*nfhour(i,3,2)
287 . -px1h3(i)*nfhour(i,3,3)-px1h4(i)*nfhour(i,3,4)
288 f31(i) =f31(i)-(hq13p+hq24n+ff)
289 f37(i) =f37(i)-(hq13p+hq24p-ff)
290 ff =-px2h1(i)*nfhour(i,3,1)-px2h2(i)*nfhour(i,3,2)
291 . -px2h3(i)*nfhour(i,3,3)-px2h4(i)*nfhour(i,3,4)
292 f32(i) =f32(i)-(hq13n-hq24n+ff)
293 f38(i) =f38(i)-(hq13n-hq24p-ff)
294 ff =-px3h1(i)*nfhour(i,3,1)-px3h2(i)*nfhour(i,3,2)
295 . -px3h3(i)*nfhour(i,3,3)-px3h4(i)*nfhour(i,3,4)
296 f33(i) =f33(i)-(-hq13n-hq24p+ff)
297 f35(i) =f35(i)-(-hq13n-hq24n-ff)
298 ff =-px4h1(i)*nfhour(i,3,1)-px4h2(i)*nfhour(i,3,2)
299 . -px4h3(i)*nfhour(i,3,3)-px4h4(i)*nfhour(i,3,4)
300 f34(i) =f34(i)-(-hq13p+hq24p+ff)
301 f36(i) =f36(i)-(-hq13p+hq24n-ff)
302 ENDDO
303
304 RETURN