OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
s10len3.F
Go to the documentation of this file.
1Copyright> OpenRadioss
2Copyright> Copyright (C) 1986-2025 Altair Engineering Inc.
3Copyright>
4Copyright> This program is free software: you can redistribute it and/or modify
5Copyright> it under the terms of the GNU Affero General Public License as published by
6Copyright> the Free Software Foundation, either version 3 of the License, or
7Copyright> (at your option) any later version.
8Copyright>
9Copyright> This program is distributed in the hope that it will be useful,
10Copyright> but WITHOUT ANY WARRANTY; without even the implied warranty of
11Copyright> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12Copyright> GNU Affero General Public License for more details.
13Copyright>
14Copyright> You should have received a copy of the GNU Affero General Public License
15Copyright> along with this program. If not, see <https://www.gnu.org/licenses/>.
16Copyright>
17Copyright>
18Copyright> Commercial Alternative: Altair Radioss Software
19Copyright>
20Copyright> As an alternative to this open-source version, Altair also offers Altair Radioss
21Copyright> software under a commercial license. Contact Altair to discuss further if the
22Copyright> commercial version may interest you: https://www.altair.com/radioss/.
23!||====================================================================
24!|| s10len3 ../engine/source/elements/solid/solide10/s10len3.f
25!||--- called by ------------------------------------------------------
26!|| s10forc3 ../engine/source/elements/solid/solide10/s10forc3.F
27!||====================================================================
28 SUBROUTINE s10len3(
29 1 VOL, NGL, DELTAX, DELTAX2,
30 2 PX, PY, PZ, VOLG,
31 3 VOLGO, RX, RY, RZ,
32 4 SX, SY, SZ, TX,
33 5 TY, TZ, NC, NEL,
34 6 MXT, PM, TAGELEM_SMS,V_PITER,
35 7 NPT, IINT, ISROT, IFORMDT)
36C-----------------------------------------------
37C I m p l i c i t T y p e s
38C-----------------------------------------------
39#include "implicit_f.inc"
40#include "comlock.inc"
41C-----------------------------------------------
42C G l o b a l P a r a m e t e r s
43C-----------------------------------------------
44#include "mvsiz_p.inc"
45C-----------------------------------------------
46C C o m m o n B l o c k s
47C-----------------------------------------------
48#include "com01_c.inc"
49#include "param_c.inc"
50#include "scr17_c.inc"
51#include "sms_c.inc"
52C-----------------------------------------------
53C D u m m y A r g u m e n t s
54C-----------------------------------------------
55 INTEGER, INTENT(IN) :: NPT
56 INTEGER, INTENT(IN) :: IINT
57 INTEGER, INTENT(IN) :: ISROT
58 INTEGER, INTENT(IN) :: IFORMDT
59 INTEGER NGL(*), NC(MVSIZ,10), NEL, MXT(MVSIZ), TAGELEM_SMS(*)
60C REAL
61 my_real
62 . vol(mvsiz,5),deltax(*),deltax2(*),volg(*),
63 . rx(*),ry(*),rz(*), sx(*),sy(*),sz(*), tx(*),ty(*),tz(*),
64 . px(mvsiz,10,5),py(mvsiz,10,5),pz(mvsiz,10,5),
65 . volgo(nel),
66 . pm(npropm,*),v_piter(nel,3,10)
67C-----------------------------------------------
68C L o c a l V a r i a b l e s
69C-----------------------------------------------
70 INTEGER I,IP,N,M,J
71 INTEGER ITER,NITER,INIT_PITER
72 my_real
73 . UL(MVSIZ,10),
74 . PXX(MVSIZ),PYY(MVSIZ),PZZ(MVSIZ),PXY(MVSIZ),PYZ(MVSIZ),PXZ(MVSIZ),FACEIGV(MVSIZ),
75 . QX(MVSIZ,10),QY(MVSIZ,10),QZ(MVSIZ,10)
76 my_real
77 . d(mvsiz),gfac,bfac,ld,p,mass,mref,fac,
78 . aa,bb,a1,a2,a3,a4,
79 . a1x,a2x,a3x,a4x,a1y,a2y,a3y,a4y,a1z,a2z,a3z,a4z,aa0
80 my_real
81 . ux(mvsiz,10), uy(mvsiz,10), uz(mvsiz,10), vx(mvsiz,10), vy(mvsiz,10), vz(mvsiz,10),
82 . nv(mvsiz), ll, bu(mvsiz,6), ninv(mvsiz)
83 my_real
84 . b1(mvsiz), b2(mvsiz), b3(mvsiz), bb4(mvsiz),
85 . c1(mvsiz), c2(mvsiz), c3(mvsiz), c4(mvsiz),
86 . d1(mvsiz), d2(mvsiz), d3(mvsiz), d4(mvsiz),
87 . p4x1(mvsiz), p4x2(mvsiz), p4x3(mvsiz), p4x4(mvsiz),
88 . p4y1(mvsiz), p4y2(mvsiz), p4y3(mvsiz), p4y4(mvsiz),
89 . p4z1(mvsiz), p4z2(mvsiz), p4z3(mvsiz), p4z4(mvsiz),
90 . det6,dd
91 INTEGER ILEAT
92 my_real
93 . ALEAT, NN, WX(10), WY(10), WZ(10)
94C-----------------------------------------------
95 FAC = one/(nine+third) ! cf FACIROT == NINE+THIRD
96 IF(idt1tet10/=0 .AND. isrot==0)THEN
97
98 IF(idtmins==0)THEN
99 DO i=1,nel
100 faceigv(i)=trhee_over_14 ! computing eigenvalues of Me-1 Ke
101 END DO
102 ELSE
103C-------- AMS (possible mix of elements inside the group)
104 DO i=1,nel
105 IF(tagelem_sms(i)==0)THEN
106 faceigv(i)=trhee_over_14 ! computing eigenvalues of Me-1 Ke, TRHEE_OVER_14 = (1/32) / (7/48) rapport des masses
107 ELSE
108 faceigv(i)=one ! compute here maximum eigenvalue K of Ke
109 ! => Retrieve K in mqviscb.F
110 ! Non diagonal mass is computed in mqviscb.F, considering Maximum eigenvalue of Me-1 Ke >= K / Minimum eigenvalue of Me
111 END IF
112 END DO
113 ENDIF
114C
115 niter = idt1tet10-1
116 IF(niter==0)THEN
117C--------------------------------------------------------------------------------------
118C /DT1TET10 Conservative Time Step
119C--------------------------------------------------------------------------------------
120 IF(iformdt==0)THEN
121
122 DO i=1,nel
123 pxx(i)=zero
124 pyy(i)=zero
125 pzz(i)=zero
126C PXY(I)=ZERO
127C PXZ(I)=ZERO
128C PYZ(I)=ZERO
129 END DO
130 DO ip=1,npt
131 DO n=1,4
132 DO i=1,nel
133 pxx(i) =pxx(i) + vol(i,ip)*px(i,n,ip)*px(i,n,ip)
134 pyy(i) =pyy(i) + vol(i,ip)*py(i,n,ip)*py(i,n,ip)
135 pzz(i) =pzz(i) + vol(i,ip)*pz(i,n,ip)*pz(i,n,ip)
136C PXY(I) =PXY(I) + VOL(I,IP)*PX(I,N,IP)*PY(I,N,IP)
137C PXZ(I) =PXZ(I) + VOL(I,IP)*PX(I,N,IP)*PZ(I,N,IP)
138C PYZ(I) =PYZ(I) + VOL(I,IP)*PY(I,N,IP)*PZ(I,N,IP)
139 ENDDO
140 ENDDO
141 DO n=5,10
142 DO i=1,nel
143 pxx(i) =pxx(i) + faceigv(i)*vol(i,ip)*px(i,n,ip)*px(i,n,ip)
144 pyy(i) =pyy(i) + faceigv(i)*vol(i,ip)*py(i,n,ip)*py(i,n,ip)
145 pzz(i) =pzz(i) + faceigv(i)*vol(i,ip)*pz(i,n,ip)*pz(i,n,ip)
146C PXY(I) =PXY(I) + FACEIGV(I)*VOL(I,IP)*PX(I,N,IP)*PY(I,N,IP)
147C PXZ(I) =PXZ(I) + FACEIGV(I)*VOL(I,IP)*PX(I,N,IP)*PZ(I,N,IP)
148C PYZ(I) =PYZ(I) + FACEIGV(I)*VOL(I,IP)*PY(I,N,IP)*PZ(I,N,IP)
149 ENDDO
150 ENDDO
151 ENDDO
152 DO i=1,nel
153 d(i) = pxx(i)+pyy(i)+pzz(i)
154 ENDDO
155
156 ELSEIF(iformdt==1)THEN
157
158 d(1:nel)=zero
159
160 gfac=pm(105,mxt(1))
161 ld =two*sqrt(max(one-gfac,zero))+one
162
163 DO ip=1,npt
164
165 DO i=1,nel
166 pxx(i)=px(i,1,ip)*px(i,1,ip)+px(i,2,ip)*px(i,2,ip)+px(i,3,ip)*px(i,3,ip)+px(i,4,ip)*px(i,4,ip)+
167 . faceigv(i)*(px(i,5,ip)*px(i,5,ip)+px(i,6,ip)*px(i,6,ip)+px(i,7,ip)*px(i,7,ip)+
168 . px(i,8,ip)*px(i,8,ip)+px(i,9,ip)*px(i,9,ip)+px(i,10,ip)*px(i,10,ip))
169 pyy(i)=py(i,1,ip)*py(i,1,ip)+py(i,2,ip)*py(i,2,ip)+py(i,3,ip)*py(i,3,ip)+py(i,4,ip)*py(i,4,ip)+
170 . faceigv(i)*(py(i,5,ip)*py(i,5,ip) +py(i,6,ip)*py(i,6,ip)+py(i,7,ip)*py(i,7,ip)+
171 . py(i,8,ip)*py(i,8,ip)+py(i,9,ip)*py(i,9,ip)+py(i,10,ip)*py(i,10,ip))
172 pzz(i)=pz(i,1,ip)*pz(i,1,ip)+pz(i,2,ip)*pz(i,2,ip)+pz(i,3,ip)*pz(i,3,ip)+pz(i,4,ip)*pz(i,4,ip)+
173 . faceigv(i)*(pz(i,5,ip)*pz(i,5,ip)+pz(i,6,ip)*pz(i,6,ip)+pz(i,7,ip)*pz(i,7,ip)+
174 . pz(i,8,ip)*pz(i,8,ip)+pz(i,9,ip)*pz(i,9,ip)+pz(i,10,ip)*pz(i,10,ip))
175 pxy(i)=px(i,1,ip)*py(i,1,ip)+px(i,2,ip)*py(i,2,ip)+px(i,3,ip)*py(i,3,ip)+px(i,4,ip)*py(i,4,ip)+
176 . faceigv(i)*(px(i,5,ip)*py(i,5,ip)+px(i,6,ip)*py(i,6,ip)+px(i,7,ip)*py(i,7,ip)+
177 . px(i,8,ip)*py(i,8,ip)+px(i,9,ip)*py(i,9,ip)+px(i,10,ip)*py(i,10,ip))
178 pxz(i)=px(i,1,ip)*pz(i,1,ip)+px(i,2,ip)*pz(i,2,ip)+px(i,3,ip)*pz(i,3,ip)+px(i,4,ip)*pz(i,4,ip)+
179 . faceigv(i)*(px(i,5,ip)*pz(i,5,ip)+px(i,6,ip)*pz(i,6,ip)+px(i,7,ip)*pz(i,7,ip)+
180 . px(i,8,ip)*pz(i,8,ip)+px(i,9,ip)*pz(i,9,ip)+px(i,10,ip)*pz(i,10,ip))
181 pyz(i)=py(i,1,ip)*pz(i,1,ip)+py(i,2,ip)*pz(i,2,ip)+py(i,3,ip)*pz(i,3,ip)+py(i,4,ip)*pz(i,4,ip)+
182 . faceigv(i)*(py(i,5,ip)*pz(i,5,ip)+py(i,6,ip)*pz(i,6,ip)+py(i,7,ip)*pz(i,7,ip)+
183 . py(i,8,ip)*pz(i,8,ip)+py(i,9,ip)*pz(i,9,ip)+py(i,10,ip)*pz(i,10,ip))
184 ENDDO
185
186 DO i=1,nel
187 aa = -(pxx(i)+pyy(i)+pzz(i))
188 bb = (pxx(i)*pyy(i)+pxx(i)*pzz(i)+pyy(i)*pzz(i)-pxy(i)**2-pxz(i)**2-pyz(i)**2)
189 p = bb-third*aa*aa
190 d(i) = d(i)+ ld * vol(i,ip) * (two*sqrt(third*max(-p,zero))-third*aa)
191 ENDDO
192 ENDDO
193
194 ELSEIF(iformdt==2)THEN
195
196 d(1:nel)=zero
197
198 gfac=pm(105,mxt(1))
199
200 DO ip=1,npt
201
202 DO i=1,nel
203 pxx(i)=px(i,1,ip)*px(i,1,ip)+px(i,2,ip)*px(i,2,ip)+px(i,3,ip)*px(i,3,ip)+px(i,4,ip)*px(i,4,ip)+
204 . faceigv(i)*(px(i,5,ip)*px(i,5,ip)+px(i,6,ip)*px(i,6,ip)+px(i,7,ip)*px(i,7,ip)+
205 . px(i,8,ip)*px(i,8,ip)+px(i,9,ip)*px(i,9,ip)+px(i,10,ip)*px(i,10,ip))
206 pyy(i)=py(i,1,ip)*py(i,1,ip)+py(i,2,ip)*py(i,2,ip)+py(i,3,ip)*py(i,3,ip)+py(i,4,ip)*py(i,4,ip)+
207 . faceigv(i)*(py(i,5,ip)*py(i,5,ip) +py(i,6,ip)*py(i,6,ip)+py(i,7,ip)*py(i,7,ip)+
208 . py(i,8,ip)*py(i,8,ip)+py(i,9,ip)*py(i,9,ip)+py(i,10,ip)*py(i,10,ip))
209 pzz(i)=pz(i,1,ip)*pz(i,1,ip)+pz(i,2,ip)*pz(i,2,ip)+pz(i,3,ip)*pz(i,3,ip)+pz(i,4,ip)*pz(i,4,ip)+
210 . faceigv(i)*(pz(i,5,ip)*pz(i,5,ip)+pz(i,6,ip)*pz(i,6,ip)+pz(i,7,ip)*pz(i,7,ip)+
211 . pz(i,8,ip)*pz(i,8,ip)+pz(i,9,ip)*pz(i,9,ip)+pz(i,10,ip)*pz(i,10,ip))
212 pxy(i)=px(i,1,ip)*py(i,1,ip)+px(i,2,ip)*py(i,2,ip)+px(i,3,ip)*py(i,3,ip)+px(i,4,ip)*py(i,4,ip)+
213 . faceigv(i)*(px(i,5,ip)*py(i,5,ip)+px(i,6,ip)*py(i,6,ip)+px(i,7,ip)*py(i,7,ip)+
214 . px(i,8,ip)*py(i,8,ip)+px(i,9,ip)*py(i,9,ip)+px(i,10,ip)*py(i,10,ip))
215 pxz(i)=px(i,1,ip)*pz(i,1,ip)+px(i,2,ip)*pz(i,2,ip)+px(i,3,ip)*pz(i,3,ip)+px(i,4,ip)*pz(i,4,ip)+
216 . faceigv(i)*(px(i,5,ip)*pz(i,5,ip)+px(i,6,ip)*pz(i,6,ip)+px(i,7,ip)*pz(i,7,ip)+
217 . px(i,8,ip)*pz(i,8,ip)+px(i,9,ip)*pz(i,9,ip)+px(i,10,ip)*pz(i,10,ip))
218 pyz(i)=py(i,1,ip)*pz(i,1,ip)+py(i,2,ip)*pz(i,2,ip)+py(i,3,ip)*pz(i,3,ip)+py(i,4,ip)*pz(i,4,ip)+
219 . faceigv(i)*(py(i,5,ip)*pz(i,5,ip)+py(i,6,ip)*pz(i,6,ip)+py(i,7,ip)*pz(i,7,ip)+
220 . py(i,8,ip)*pz(i,8,ip)+py(i,9,ip)*pz(i,9,ip)+py(i,10,ip)*pz(i,10,ip))
221 ENDDO
222
223 DO i=1,nel
224 aa = -(pxx(i)+pyy(i)+pzz(i))
225 bb = gfac*(pxx(i)*pyy(i)+pxx(i)*pzz(i)+pyy(i)*pzz(i)-pxy(i)**2-pxz(i)**2-pyz(i)**2)
226 p = bb-third*aa*aa
227
228 d(i) = d(i)+vol(i,ip)*(two*sqrt(third*max(-p,zero))-third*aa)
229 ENDDO
230
231 ENDDO
232
233 END IF ! IF(IFORMDT == ...)
234
235 ELSE ! IF(NITER==0)THEN
236C--------------------------------------------------------------------------------------
237C /DT1TET10/NITER
238C Iterative Power (Could be adapted vs IFORM taking into account Hook's Matrix)
239C--------------------------------------------------------------------------------------
240 init_piter=0
241 IF(ncycle==0)THEN
242 init_piter=1
243 DO n=1,10
244 IF(v_piter(1,1,n)/=zero.OR.v_piter(1,2,n)/=zero.OR.v_piter(1,3,n)/=zero)init_piter=0
245 END DO
246 END IF
247
248 IF(init_piter==1)THEN
249C
250C Initialize at 1st time DT1TET10 has been encountered... Note : U normalized
251 ileat=0
252 DO n=1,10
253 ileat=mod(25173*ileat+13849,65536)
254 aleat=(ileat-32768.)/32768.
255 wx(n)=em03*aleat
256 ileat=mod(25173*ileat+13849,65536)
257 aleat=(ileat-32768.)/32768.
258 wy(n)=em03*aleat
259 ileat=mod(25173*ileat+13849,65536)
260 aleat=(ileat-32768.)/32768.
261 wz(n)=em03*aleat
262 END DO
263
264 nn = zero
265 DO n=1,10
266 nn = nn + wx(n)*wx(n)+wy(n)*wy(n)+wz(n)*wz(n)
267 END DO
268 nn=one/max(em20,nn)
269
270 DO n=1,10
271 wx(n)=nn * wx(n)
272 wy(n)=nn * wy(n)
273 wz(n)=nn * wz(n)
274 END DO
275
276 DO n=1,10
277 DO i=1,nel
278 ux(i,n)=wx(n)
279 uy(i,n)=wy(n)
280 uz(i,n)=wz(n)
281 END DO
282 END DO
283 ELSE
284C
285C Get vector from previous cycle / run
286 ux(1:nel,1:10)=v_piter(1:nel,1,1:10)
287 uy(1:nel,1:10)=v_piter(1:nel,2,1:10)
288 uz(1:nel,1:10)=v_piter(1:nel,3,1:10)
289 END IF
290
291 IF(iformdt==2)THEN
292C
293C It is fine so far as both ratio 3B/rho0 c**2 and G / rho0 c**2
294C are computed using Gmax and c is also computed in the material using Gmax.
295 gfac=half*pm(105,mxt(1)) ! G/(B+4/3G)
296 bfac=three-four*gfac ! 3B/(B+4/3G)
297 ELSE ! conservative formulation in all other cases
298 gfac=half ! over-estimation of G/(B+4/3G)
299 bfac=three !
300 END IF
301
302 DO iter=1,niter
303
304
305 vx(1:nel,1:10)=zero
306 vy(1:nel,1:10)=zero
307 vz(1:nel,1:10)=zero
308
309 DO ip=1,npt
310C
311C Bip.U
312 bu(1:nel,1:6)=zero
313 DO n=1,10
314 DO i=1,nel
315 bu(i,1)=bu(i,1)+px(i,n,ip)*ux(i,n)
316 bu(i,2)=bu(i,2)+py(i,n,ip)*uy(i,n)
317 bu(i,3)=bu(i,3)+pz(i,n,ip)*uz(i,n)
318 bu(i,4)=bu(i,4)+py(i,n,ip)*ux(i,n)+px(i,n,ip)*uy(i,n)
319 bu(i,5)=bu(i,5)+pz(i,n,ip)*uy(i,n)+py(i,n,ip)*uz(i,n)
320 bu(i,6)=bu(i,6)+pz(i,n,ip)*ux(i,n)+px(i,n,ip)*uz(i,n)
321 END DO
322 END DO
323C
324 DO j=1,3
325 DO i=1,nel
326 bu(i,j)=bfac*bu(i,j)
327 END DO
328 END DO
329 DO j=4,6
330 DO i=1,nel
331 bu(i,j)=gfac*bu(i,j)
332 END DO
333 END DO
334C
335C Vol_ip * Transpose(Bip).Bip.U
336 DO n=1,10
337 DO i=1,nel
338 vx(i,n)=vx(i,n)+vol(i,ip)*(px(i,n,ip)*bu(i,1)+py(i,n,ip)*bu(i,4)+pz(i,n,ip)*bu(i,6))
339 vy(i,n)=vy(i,n)+vol(i,ip)*(py(i,n,ip)*bu(i,2)+px(i,n,ip)*bu(i,4)+pz(i,n,ip)*bu(i,5))
340 vz(i,n)=vz(i,n)+vol(i,ip)*(pz(i,n,ip)*bu(i,3)+py(i,n,ip)*bu(i,5)+px(i,n,ip)*bu(i,6))
341 END DO
342 END DO
343
344 END DO
345
346 nv(1:nel) = zero
347 DO n=1,4
348 DO i=1,nel
349 nv(i) = nv(i) + vx(i,n)*vx(i,n)+vy(i,n)*vy(i,n)+vz(i,n)*vz(i,n)
350 END DO
351 END DO
352 DO n=5,10
353 DO i=1,nel
354 vx(i,n)=faceigv(i)*vx(i,n)
355 vy(i,n)=faceigv(i)*vy(i,n)
356 vz(i,n)=faceigv(i)*vz(i,n)
357 nv(i) = nv(i) + vx(i,n)*vx(i,n)+vy(i,n)*vy(i,n)+vz(i,n)*vz(i,n)
358 END DO
359 END DO
360 DO i=1,nel
361 nv(i) =sqrt(nv(i))
362 ninv(i)=one/max(em20,nv(i))
363 END DO
364 DO n=1,10
365 DO i=1,nel
366 vx(i,n)=ninv(i) * vx(i,n)
367 vy(i,n)=ninv(i) * vy(i,n)
368 vz(i,n)=ninv(i) * vz(i,n)
369 END DO
370 END DO
371 ux(1:nel,1:10)=vx(1:nel,1:10)
372 uy(1:nel,1:10)=vy(1:nel,1:10)
373 uz(1:nel,1:10)=vz(1:nel,1:10)
374
375 END DO ! DO ITER=1,NITER
376C
377 DO i=1,nel
378 d(i)=nv(i)
379 END DO
380 v_piter(1:nel,1,1:10)=ux(1:nel,1:10)
381 v_piter(1:nel,2,1:10)=uy(1:nel,1:10)
382 v_piter(1:nel,3,1:10)=uz(1:nel,1:10)
383C
384C 2 * Max diagonal as an (under) estimation of Lambda
385C Because Iterative Power may be too much under-estimated before convergence.
386 DO n=1,10
387 DO i=1,nel
388 ul(i,n) = zero
389 ENDDO
390 ENDDO
391C
392 DO ip=1,npt
393 DO n=1,10
394 DO i=1,nel
395 ul(i,n) =ul(i,n) + vol(i,ip) *
396 + ( px(i,n,ip)*px(i,n,ip) + py(i,n,ip)*py(i,n,ip)
397 + + pz(i,n,ip)*pz(i,n,ip) )
398 ENDDO
399 ENDDO
400 ENDDO
401C
402 DO i=1,nel
403 aa = max(ul(i,1),ul(i,2),ul(i,3),ul(i,4))
404 bb = faceigv(i)*max(ul(i,5),ul(i,6),ul(i,7),ul(i,8),ul(i,9),ul(i,10))
405 d(i) = max(d(i),two*max(aa,bb))
406 ENDDO
407C
408 END IF ! IF(NITER==0)THEN
409
410 DO i=1,nel
411C
412C DELTAX2 is not used w/IDT1TET10
413 deltax2(i) = one
414C
415 mass = volgo(i) ! (multiplied by RHOG0)
416
417 mref = one/thirty2 * mass
418 deltax(i) = two*sqrt(mref/d(i))
419 ENDDO
420
421 ELSEIF(idt1tet10/=0 .AND. isrot==2)THEN
422
423 niter = idt1tet10-1
424 IF(niter==0)THEN
425C--------------------------------------------------------------------------------------
426C /DT1TET10 Conservative Time Step
427C--------------------------------------------------------------------------------------
428 IF(iformdt==0)THEN
429
430 DO i=1,nel
431 pxx(i)=zero
432 pyy(i)=zero
433 pzz(i)=zero
434C PXY(I)=ZERO
435C PXZ(I)=ZERO
436C PYZ(I)=ZERO
437 END DO
438
439 DO ip=1,npt
440c-----------------------------------------------------------------------------
441C Spectral Radius(M-1K)
442C-----------------------------------------------------------------------------
443C
444C Q = M-1 Transpose(B)
445 DO i=1,nel
446 qx(i,1)=px(i,1,ip)+half*(px(i,5,ip)+px(i,7,ip)+px(i,8,ip))
447 qx(i,2)=px(i,2,ip)+half*(px(i,5,ip)+px(i,6,ip)+px(i,9,ip))
448 qx(i,3)=px(i,3,ip)+half*(px(i,6,ip)+px(i,7,ip)+px(i,10,ip))
449 qx(i,4)=px(i,4,ip)+half*(px(i,8,ip)+px(i,9,ip)+px(i,10,ip))
450C QX(I,5)=HALF*(PX(I,1,IP)+PX(I,2,IP))+(HALF+ONE/FACIROT)*PX(I,5,IP)
451C +FOURTH(PX(I,6,IP)+PX(I,7,IP)+PX(I,8,IP)+PX(I,9,IP)
452 qx(i,5) =half*(qx(i,1)+qx(i,2))+fac*px(i,5,ip)
453 qx(i,6) =half*(qx(i,2)+qx(i,3))+fac*px(i,6,ip)
454 qx(i,7) =half*(qx(i,1)+qx(i,3))+fac*px(i,7,ip)
455 qx(i,8) =half*(qx(i,1)+qx(i,4))+fac*px(i,8,ip)
456 qx(i,9) =half*(qx(i,2)+qx(i,4))+fac*px(i,9,ip)
457 qx(i,10)=half*(qx(i,3)+qx(i,4))+fac*px(i,10,ip)
458
459 qy(i,1)=py(i,1,ip)+half*(py(i,5,ip)+py(i,7,ip)+py(i,8,ip))
460 qy(i,2)=py(i,2,ip)+half*(py(i,5,ip)+py(i,6,ip)+py(i,9,ip))
461 qy(i,3)=py(i,3,ip)+half*(py(i,6,ip)+py(i,7,ip)+py(i,10,ip))
462 qy(i,4)=py(i,4,ip)+half*(py(i,8,ip)+py(i,9,ip)+py(i,10,ip))
463 qy(i,5) =half*(qy(i,1)+qy(i,2))+fac*py(i,5,ip)
464 qy(i,6) =half*(qy(i,2)+qy(i,3))+fac*py(i,6,ip)
465 qy(i,7) =half*(qy(i,1)+qy(i,3))+fac*py(i,7,ip)
466 qy(i,8) =half*(qy(i,1)+qy(i,4))+fac*py(i,8,ip)
467 qy(i,9) =half*(qy(i,2)+qy(i,4))+fac*py(i,9,ip)
468 qy(i,10)=half*(qy(i,3)+qy(i,4))+fac*py(i,10,ip)
469
470 qz(i,1)=pz(i,1,ip)+half*(pz(i,5,ip)+pz(i,7,ip)+pz(i,8,ip))
471 qz(i,2)=pz(i,2,ip)+half*(pz(i,5,ip)+pz(i,6,ip)+pz(i,9,ip))
472 qz(i,3)=pz(i,3,ip)+half*(pz(i,6,ip)+pz(i,7,ip)+pz(i,10,ip))
473 qz(i,4)=pz(i,4,ip)+half*(pz(i,8,ip)+pz(i,9,ip)+pz(i,10,ip))
474 qz(i,5) =half*(qz(i,1)+qz(i,2))+fac*pz(i,5,ip)
475 qz(i,6) =half*(qz(i,2)+qz(i,3))+fac*pz(i,6,ip)
476 qz(i,7) =half*(qz(i,1)+qz(i,3))+fac*pz(i,7,ip)
477 qz(i,8) =half*(qz(i,1)+qz(i,4))+fac*pz(i,8,ip)
478 qz(i,9) =half*(qz(i,2)+qz(i,4))+fac*pz(i,9,ip)
479 qz(i,10)=half*(qz(i,3)+qz(i,4))+fac*pz(i,10,ip)
480
481 END DO
482C
483C B M-1 Transpose(B)
484 DO i=1,nel
485 pxx(i)=pxx(i)+
486 . px(i,1,ip)*qx(i,1)+px(i,2,ip)*qx(i,2)+px(i,3,ip)*qx(i,3)+px(i,4,ip)*qx(i,4)+
487 . px(i,5,ip)*qx(i,5)+px(i,6,ip)*qx(i,6)+px(i,7,ip)*qx(i,7)+
488 . px(i,8,ip)*qx(i,8)+px(i,9,ip)*qx(i,9)+px(i,10,ip)*qx(i,10)
489 pyy(i)=pyy(i)+
490 . py(i,1,ip)*qy(i,1)+py(i,2,ip)*qy(i,2)+py(i,3,ip)*qy(i,3)+py(i,4,ip)*qy(i,4)+
491 . py(i,5,ip)*qy(i,5) +py(i,6,ip)*qy(i,6)+py(i,7,ip)*qy(i,7)+
492 . py(i,8,ip)*qy(i,8)+py(i,9,ip)*qy(i,9)+py(i,10,ip)*qy(i,10)
493 pzz(i)=pzz(i)+
494 . pz(i,1,ip)*qz(i,1)+pz(i,2,ip)*qz(i,2)+pz(i,3,ip)*qz(i,3)+pz(i,4,ip)*qz(i,4)+
495 . pz(i,5,ip)*qz(i,5)+pz(i,6,ip)*qz(i,6)+pz(i,7,ip)*qz(i,7)+
496 . pz(i,8,ip)*qz(i,8)+pz(i,9,ip)*qz(i,9)+pz(i,10,ip)*qz(i,10)
497C PXY(I)=PX(I,1,IP)*QY(I,1)+PX(I,2,IP)*QY(I,2)+PX(I,3,IP)*QY(I,3)+PX(I,4,IP)*QY(I,4)+
498C . PX(I,5,IP)*QY(I,5)+PX(I,6,IP)*QY(I,6)+PX(I,7,IP)*QY(I,7)+
499C . PX(I,8,IP)*QY(I,8)+PX(I,9,IP)*QY(I,9)+PX(I,10,IP)*QY(I,10)
500C PXZ(I)=PX(I,1,IP)*QZ(I,1)+PX(I,2,IP)*QZ(I,2)+PX(I,3,IP)*QZ(I,3)+PX(I,4,IP)*QZ(I,4)+
501C . PX(I,5,IP)*QZ(I,5)+PX(I,6,IP)*QZ(I,6)+PX(I,7,IP)*QZ(I,7)+
502C . PX(I,8,IP)*QZ(I,8)+PX(I,9,IP)*QZ(I,9)+PX(I,10,IP)*QZ(I,10)
503C PYZ(I)=PY(I,1,IP)*QZ(I,1)+PY(I,2,IP)*QZ(I,2)+PY(I,3,IP)*QZ(I,3)+PY(I,4,IP)*QZ(I,4)+
504C . PY(I,5,IP)*QZ(I,5)+PY(I,6,IP)*QZ(I,6)+PY(I,7,IP)*QZ(I,7)+
505C . PY(I,8,IP)*QZ(I,8)+PY(I,9,IP)*QZ(I,9)+PY(I,10,IP)*QZ(I,10)
506 ENDDO
507
508 ENDDO
509
510 DO i=1,nel
511 d(i) = pxx(i)+pyy(i)+pzz(i)
512 ENDDO
513
514 ELSEIF(iformdt==1)THEN
515
516 d(1:nel)=zero
517
518 gfac=pm(105,mxt(1))
519 ld =two*sqrt(max(one-gfac,zero))+one
520
521 DO ip=1,npt
522c-----------------------------------------------------------------------------
523C Spectral Radius(M-1K)
524C-----------------------------------------------------------------------------
525C
526C Q = M-1 Transpose(B)
527 DO i=1,nel
528 qx(i,1)=px(i,1,ip)+half*(px(i,5,ip)+px(i,7,ip)+px(i,8,ip))
529 qx(i,2)=px(i,2,ip)+half*(px(i,5,ip)+px(i,6,ip)+px(i,9,ip))
530 qx(i,3)=px(i,3,ip)+half*(px(i,6,ip)+px(i,7,ip)+px(i,10,ip))
531 qx(i,4)=px(i,4,ip)+half*(px(i,8,ip)+px(i,9,ip)+px(i,10,ip))
532C QX(I,5)=HALF*(PX(I,1,IP)+PX(I,2,IP))+(HALF+ONE/FACIROT)*PX(I,5,IP)
533C +FOURTH(PX(I,6,IP)+PX(I,7,IP)+PX(I,8,IP)+PX(I,9,IP)
534 qx(i,5) =half*(qx(i,1)+qx(i,2))+fac*px(i,5,ip)
535 qx(i,6) =half*(qx(i,2)+qx(i,3))+fac*px(i,6,ip)
536 qx(i,7) =half*(qx(i,1)+qx(i,3))+fac*px(i,7,ip)
537 qx(i,8) =half*(qx(i,1)+qx(i,4))+fac*px(i,8,ip)
538 qx(i,9) =half*(qx(i,2)+qx(i,4))+fac*px(i,9,ip)
539 qx(i,10)=half*(qx(i,3)+qx(i,4))+fac*px(i,10,ip)
540
541 qy(i,1)=py(i,1,ip)+half*(py(i,5,ip)+py(i,7,ip)+py(i,8,ip))
542 qy(i,2)=py(i,2,ip)+half*(py(i,5,ip)+py(i,6,ip)+py(i,9,ip))
543 qy(i,3)=py(i,3,ip)+half*(py(i,6,ip)+py(i,7,ip)+py(i,10,ip))
544 qy(i,4)=py(i,4,ip)+half*(py(i,8,ip)+py(i,9,ip)+py(i,10,ip))
545 qy(i,5) =half*(qy(i,1)+qy(i,2))+fac*py(i,5,ip)
546 qy(i,6) =half*(qy(i,2)+qy(i,3))+fac*py(i,6,ip)
547 qy(i,7) =half*(qy(i,1)+qy(i,3))+fac*py(i,7,ip)
548 qy(i,8) =half*(qy(i,1)+qy(i,4))+fac*py(i,8,ip)
549 qy(i,9) =half*(qy(i,2)+qy(i,4))+fac*py(i,9,ip)
550 qy(i,10)=half*(qy(i,3)+qy(i,4))+fac*py(i,10,ip)
551
552 qz(i,1)=pz(i,1,ip)+half*(pz(i,5,ip)+pz(i,7,ip)+pz(i,8,ip))
553 qz(i,2)=pz(i,2,ip)+half*(pz(i,5,ip)+pz(i,6,ip)+pz(i,9,ip))
554 qz(i,3)=pz(i,3,ip)+half*(pz(i,6,ip)+pz(i,7,ip)+pz(i,10,ip))
555 qz(i,4)=pz(i,4,ip)+half*(pz(i,8,ip)+pz(i,9,ip)+pz(i,10,ip))
556 qz(i,5) =half*(qz(i,1)+qz(i,2))+fac*pz(i,5,ip)
557 qz(i,6) =half*(qz(i,2)+qz(i,3))+fac*pz(i,6,ip)
558 qz(i,7) =half*(qz(i,1)+qz(i,3))+fac*pz(i,7,ip)
559 qz(i,8) =half*(qz(i,1)+qz(i,4))+fac*pz(i,8,ip)
560 qz(i,9) =half*(qz(i,2)+qz(i,4))+fac*pz(i,9,ip)
561 qz(i,10)=half*(qz(i,3)+qz(i,4))+fac*pz(i,10,ip)
562
563 END DO
564C
565C B M-1 Transpose(B)
566 DO i=1,nel
567 pxx(i)=px(i,1,ip)*qx(i,1)+px(i,2,ip)*qx(i,2)+px(i,3,ip)*qx(i,3)+px(i,4,ip)*qx(i,4)+
568 . px(i,5,ip)*qx(i,5)+px(i,6,ip)*qx(i,6)+px(i,7,ip)*qx(i,7)+
569 . px(i,8,ip)*qx(i,8)+px(i,9,ip)*qx(i,9)+px(i,10,ip)*qx(i,10)
570 pyy(i)=py(i,1,ip)*qy(i,1)+py(i,2,ip)*qy(i,2)+py(i,3,ip)*qy(i,3)+py(i,4,ip)*qy(i,4)+
571 . py(i,5,ip)*qy(i,5) +py(i,6,ip)*qy(i,6)+py(i,7,ip)*qy(i,7)+
572 . py(i,8,ip)*qy(i,8)+py(i,9,ip)*qy(i,9)+py(i,10,ip)*qy(i,10)
573 pzz(i)=pz(i,1,ip)*qz(i,1)+pz(i,2,ip)*qz(i,2)+pz(i,3,ip)*qz(i,3)+pz(i,4,ip)*qz(i,4)+
574 . pz(i,5,ip)*qz(i,5)+pz(i,6,ip)*qz(i,6)+pz(i,7,ip)*qz(i,7)+
575 . pz(i,8,ip)*qz(i,8)+pz(i,9,ip)*qz(i,9)+pz(i,10,ip)*qz(i,10)
576 pxy(i)=px(i,1,ip)*qy(i,1)+px(i,2,ip)*qy(i,2)+px(i,3,ip)*qy(i,3)+px(i,4,ip)*qy(i,4)+
577 . px(i,5,ip)*qy(i,5)+px(i,6,ip)*qy(i,6)+px(i,7,ip)*qy(i,7)+
578 . px(i,8,ip)*qy(i,8)+px(i,9,ip)*qy(i,9)+px(i,10,ip)*qy(i,10)
579 pxz(i)=px(i,1,ip)*qz(i,1)+px(i,2,ip)*qz(i,2)+px(i,3,ip)*qz(i,3)+px(i,4,ip)*qz(i,4)+
580 . px(i,5,ip)*qz(i,5)+px(i,6,ip)*qz(i,6)+px(i,7,ip)*qz(i,7)+
581 . px(i,8,ip)*qz(i,8)+px(i,9,ip)*qz(i,9)+px(i,10,ip)*qz(i,10)
582 pyz(i)=py(i,1,ip)*qz(i,1)+py(i,2,ip)*qz(i,2)+py(i,3,ip)*qz(i,3)+py(i,4,ip)*qz(i,4)+
583 . py(i,5,ip)*qz(i,5)+py(i,6,ip)*qz(i,6)+py(i,7,ip)*qz(i,7)+
584 . py(i,8,ip)*qz(i,8)+py(i,9,ip)*qz(i,9)+py(i,10,ip)*qz(i,10)
585 ENDDO
586
587 DO i=1,nel
588 aa = -(pxx(i)+pyy(i)+pzz(i))
589 bb = (pxx(i)*pyy(i)+pxx(i)*pzz(i)+pyy(i)*pzz(i)-pxy(i)**2-pxz(i)**2-pyz(i)**2)
590 p = bb-third*aa*aa
591 d(i) = d(i)+ ld * vol(i,ip) * (two*sqrt(third*max(-p,zero))-third*aa)
592 ENDDO
593 ENDDO
594
595 ELSEIF(iformdt==2)THEN
596
597 d(1:nel)=zero
598
599 gfac=pm(105,mxt(1))
600
601 DO ip=1,npt
602c-----------------------------------------------------------------------------
603C Spectral Radius(M-1K)
604C-----------------------------------------------------------------------------
605C
606C Q = M-1 Transpose(B)
607 DO i=1,nel
608 qx(i,1)=px(i,1,ip)+half*(px(i,5,ip)+px(i,7,ip)+px(i,8,ip))
609 qx(i,2)=px(i,2,ip)+half*(px(i,5,ip)+px(i,6,ip)+px(i,9,ip))
610 qx(i,3)=px(i,3,ip)+half*(px(i,6,ip)+px(i,7,ip)+px(i,10,ip))
611 qx(i,4)=px(i,4,ip)+half*(px(i,8,ip)+px(i,9,ip)+px(i,10,ip))
612C QX(I,5)=HALF*(PX(I,1,IP)+PX(I,2,IP))+(HALF+ONE/FACIROT)*PX(I,5,IP)
613C +FOURTH(PX(I,6,IP)+PX(I,7,IP)+PX(I,8,IP)+PX(I,9,IP)
614 qx(i,5) =half*(qx(i,1)+qx(i,2))+fac*px(i,5,ip)
615 qx(i,6) =half*(qx(i,2)+qx(i,3))+fac*px(i,6,ip)
616 qx(i,7) =half*(qx(i,1)+qx(i,3))+fac*px(i,7,ip)
617 qx(i,8) =half*(qx(i,1)+qx(i,4))+fac*px(i,8,ip)
618 qx(i,9) =half*(qx(i,2)+qx(i,4))+fac*px(i,9,ip)
619 qx(i,10)=half*(qx(i,3)+qx(i,4))+fac*px(i,10,ip)
620
621 qy(i,1)=py(i,1,ip)+half*(py(i,5,ip)+py(i,7,ip)+py(i,8,ip))
622 qy(i,2)=py(i,2,ip)+half*(py(i,5,ip)+py(i,6,ip)+py(i,9,ip))
623 qy(i,3)=py(i,3,ip)+half*(py(i,6,ip)+py(i,7,ip)+py(i,10,ip))
624 qy(i,4)=py(i,4,ip)+half*(py(i,8,ip)+py(i,9,ip)+py(i,10,ip))
625 qy(i,5) =half*(qy(i,1)+qy(i,2))+fac*py(i,5,ip)
626 qy(i,6) =half*(qy(i,2)+qy(i,3))+fac*py(i,6,ip)
627 qy(i,7) =half*(qy(i,1)+qy(i,3))+fac*py(i,7,ip)
628 qy(i,8) =half*(qy(i,1)+qy(i,4))+fac*py(i,8,ip)
629 qy(i,9) =half*(qy(i,2)+qy(i,4))+fac*py(i,9,ip)
630 qy(i,10)=half*(qy(i,3)+qy(i,4))+fac*py(i,10,ip)
631
632 qz(i,1)=pz(i,1,ip)+half*(pz(i,5,ip)+pz(i,7,ip)+pz(i,8,ip))
633 qz(i,2)=pz(i,2,ip)+half*(pz(i,5,ip)+pz(i,6,ip)+pz(i,9,ip))
634 qz(i,3)=pz(i,3,ip)+half*(pz(i,6,ip)+pz(i,7,ip)+pz(i,10,ip))
635 qz(i,4)=pz(i,4,ip)+half*(pz(i,8,ip)+pz(i,9,ip)+pz(i,10,ip))
636 qz(i,5) =half*(qz(i,1)+qz(i,2))+fac*pz(i,5,ip)
637 qz(i,6) =half*(qz(i,2)+qz(i,3))+fac*pz(i,6,ip)
638 qz(i,7) =half*(qz(i,1)+qz(i,3))+fac*pz(i,7,ip)
639 qz(i,8) =half*(qz(i,1)+qz(i,4))+fac*pz(i,8,ip)
640 qz(i,9) =half*(qz(i,2)+qz(i,4))+fac*pz(i,9,ip)
641 qz(i,10)=half*(qz(i,3)+qz(i,4))+fac*pz(i,10,ip)
642
643 END DO
644C
645C B M-1 Transpose(B)
646 DO i=1,nel
647 pxx(i)=px(i,1,ip)*qx(i,1)+px(i,2,ip)*qx(i,2)+px(i,3,ip)*qx(i,3)+px(i,4,ip)*qx(i,4)+
648 . px(i,5,ip)*qx(i,5)+px(i,6,ip)*qx(i,6)+px(i,7,ip)*qx(i,7)+
649 . px(i,8,ip)*qx(i,8)+px(i,9,ip)*qx(i,9)+px(i,10,ip)*qx(i,10)
650 pyy(i)=py(i,1,ip)*qy(i,1)+py(i,2,ip)*qy(i,2)+py(i,3,ip)*qy(i,3)+py(i,4,ip)*qy(i,4)+
651 . py(i,5,ip)*qy(i,5) +py(i,6,ip)*qy(i,6)+py(i,7,ip)*qy(i,7)+
652 . py(i,8,ip)*qy(i,8)+py(i,9,ip)*qy(i,9)+py(i,10,ip)*qy(i,10)
653 pzz(i)=pz(i,1,ip)*qz(i,1)+pz(i,2,ip)*qz(i,2)+pz(i,3,ip)*qz(i,3)+pz(i,4,ip)*qz(i,4)+
654 . pz(i,5,ip)*qz(i,5)+pz(i,6,ip)*qz(i,6)+pz(i,7,ip)*qz(i,7)+
655 . pz(i,8,ip)*qz(i,8)+pz(i,9,ip)*qz(i,9)+pz(i,10,ip)*qz(i,10)
656 pxy(i)=px(i,1,ip)*qy(i,1)+px(i,2,ip)*qy(i,2)+px(i,3,ip)*qy(i,3)+px(i,4,ip)*qy(i,4)+
657 . px(i,5,ip)*qy(i,5)+px(i,6,ip)*qy(i,6)+px(i,7,ip)*qy(i,7)+
658 . px(i,8,ip)*qy(i,8)+px(i,9,ip)*qy(i,9)+px(i,10,ip)*qy(i,10)
659 pxz(i)=px(i,1,ip)*qz(i,1)+px(i,2,ip)*qz(i,2)+px(i,3,ip)*qz(i,3)+px(i,4,ip)*qz(i,4)+
660 . px(i,5,ip)*qz(i,5)+px(i,6,ip)*qz(i,6)+px(i,7,ip)*qz(i,7)+
661 . px(i,8,ip)*qz(i,8)+px(i,9,ip)*qz(i,9)+px(i,10,ip)*qz(i,10)
662 pyz(i)=py(i,1,ip)*qz(i,1)+py(i,2,ip)*qz(i,2)+py(i,3,ip)*qz(i,3)+py(i,4,ip)*qz(i,4)+
663 . py(i,5,ip)*qz(i,5)+py(i,6,ip)*qz(i,6)+py(i,7,ip)*qz(i,7)+
664 . py(i,8,ip)*qz(i,8)+py(i,9,ip)*qz(i,9)+py(i,10,ip)*qz(i,10)
665 ENDDO
666C-------------------------------------------------------------------------------
667 DO i=1,nel
668 aa = -(pxx(i)+pyy(i)+pzz(i))
669 bb = gfac*(pxx(i)*pyy(i)+pxx(i)*pzz(i)+pyy(i)*pzz(i)-pxy(i)**2-pxz(i)**2-pyz(i)**2)
670 p = bb-third*aa*aa
671
672 d(i) = d(i)+vol(i,ip)*(two*sqrt(third*max(-p,zero))-third*aa)
673 ENDDO
674
675 ENDDO
676
677 END if! IF(IFORMDT == ...)
678
679 ELSE ! IF(NITER==0)THEN
680C--------------------------------------------------------------------------------------
681C /DT1TET10/NITER
682C Iterative Power (Could be adapted vs IFORM taking into account Hook's Matrix)
683C--------------------------------------------------------------------------------------
684 init_piter=0
685 IF(ncycle==0)THEN
686 init_piter=1
687 DO n=1,10
688 IF(v_piter(1,1,n)/=zero.OR.v_piter(1,2,n)/=zero.OR.v_piter(1,3,n)/=zero)init_piter=0
689 END DO
690
691 END IF
692C
693
694 IF(init_piter==1)THEN
695C
696C Initialize at 1st time DT1TET10 has been encountered... Note : U normalized
697 ileat=0
698 DO n=1,10
699 ileat=mod(25173*ileat+13849,65536)
700 aleat=(ileat-32768.)/32768.
701 wx(n)=em03*aleat
702 ileat=mod(25173*ileat+13849,65536)
703 aleat=(ileat-32768.)/32768.
704 wy(n)=em03*aleat
705 ileat=mod(25173*ileat+13849,65536)
706 aleat=(ileat-32768.)/32768.
707 wz(n)=em03*aleat
708 END DO
709
710 nn = zero
711 DO n=1,10
712 nn = nn + wx(n)*wx(n)+wy(n)*wy(n)+wz(n)*wz(n)
713 END DO
714 nn=one/max(em20,nn)
715
716 DO n=1,10
717 wx(n)=nn * wx(n)
718 wy(n)=nn * wy(n)
719 wz(n)=nn * wz(n)
720 END DO
721
722 DO n=1,10
723 DO i=1,nel
724 ux(i,n)=wx(n)
725 uy(i,n)=wy(n)
726 uz(i,n)=wz(n)
727 END DO
728 END DO
729 ELSE
730C
731C Get vector from previous cycle / run
732 ux(1:nel,1:10)=v_piter(1:nel,1,1:10)
733 uy(1:nel,1:10)=v_piter(1:nel,2,1:10)
734 uz(1:nel,1:10)=v_piter(1:nel,3,1:10)
735 END IF
736
737 IF(iformdt==2)THEN
738C
739C It is fine so far as both ratio 3B/rho0 c**2 and G / rho0 c**2
740C are computed using Gmax and c is also computed in the material using Gmax.
741 gfac=half*pm(105,mxt(1)) ! G/(B+4/3G)
742 bfac=three-four*gfac ! 3B/(B+4/3G)
743 ELSE ! conservative formulation in all other cases
744 gfac=half ! over-estimation of G/(B+4/3G)
745 bfac=three !
746 END IF
747
748 DO iter=1,niter
749
750 vx(1:nel,1:10)=zero
751 vy(1:nel,1:10)=zero
752 vz(1:nel,1:10)=zero
753
754 DO ip=1,npt
755C
756C Bip.U
757 bu(1:nel,1:6)=zero
758 DO n=1,10
759 DO i=1,nel
760 bu(i,1)=bu(i,1)+px(i,n,ip)*ux(i,n)
761 bu(i,2)=bu(i,2)+py(i,n,ip)*uy(i,n)
762 bu(i,3)=bu(i,3)+pz(i,n,ip)*uz(i,n)
763 bu(i,4)=bu(i,4)+py(i,n,ip)*ux(i,n)+px(i,n,ip)*uy(i,n)
764 bu(i,5)=bu(i,5)+pz(i,n,ip)*uy(i,n)+py(i,n,ip)*uz(i,n)
765 bu(i,6)=bu(i,6)+pz(i,n,ip)*ux(i,n)+px(i,n,ip)*uz(i,n)
766 END DO
767 END DO
768C
769 DO j=1,3
770 DO i=1,nel
771 bu(i,j)=bfac*bu(i,j)
772 END DO
773 END DO
774 DO j=4,6
775 DO i=1,nel
776 bu(i,j)=gfac*bu(i,j)
777 END DO
778 END DO
779C
780C Vol_ip * Transpose(Bip).Bip.U
781 DO n=1,10
782 DO i=1,nel
783 vx(i,n)=vx(i,n)+vol(i,ip)*(px(i,n,ip)*bu(i,1)+py(i,n,ip)*bu(i,4)+pz(i,n,ip)*bu(i,6))
784 vy(i,n)=vy(i,n)+vol(i,ip)*(py(i,n,ip)*bu(i,2)+px(i,n,ip)*bu(i,4)+pz(i,n,ip)*bu(i,5))
785 vz(i,n)=vz(i,n)+vol(i,ip)*(pz(i,n,ip)*bu(i,3)+py(i,n,ip)*bu(i,5)+px(i,n,ip)*bu(i,6))
786 END DO
787 END DO
788
789 END DO
790
791C M-1 K . U
792 DO i=1,nel
793 vx(i,1)=vx(i,1)+half*(vx(i,5)+vx(i,7)+vx(i,8))
794 vx(i,2)=vx(i,2)+half*(vx(i,5)+vx(i,6)+vx(i,9))
795 vx(i,3)=vx(i,3)+half*(vx(i,6)+vx(i,7)+vx(i,10))
796 vx(i,4)=vx(i,4)+half*(vx(i,8)+vx(i,9)+vx(i,10))
797C VX(I,5)=HALF*(VX(I,1)+VX(I,2))+(HALF+ONE/FACIROT)*VX(I,5)
798C +FOURTH(VX(I,6)+VX(I,7)+VX(I,8)+VX(I,9)
799 vx(i,5) =half*(vx(i,1)+vx(i,2))+fac*vx(i,5)
800 vx(i,6) =half*(vx(i,2)+vx(i,3))+fac*vx(i,6)
801 vx(i,7) =half*(vx(i,1)+vx(i,3))+fac*vx(i,7)
802 vx(i,8) =half*(vx(i,1)+vx(i,4))+fac*vx(i,8)
803 vx(i,9) =half*(vx(i,2)+vx(i,4))+fac*vx(i,9)
804 vx(i,10)=half*(vx(i,3)+vx(i,4))+fac*vx(i,10)
805
806 vy(i,1)=vy(i,1)+half*(vy(i,5)+vy(i,7)+vy(i,8))
807 vy(i,2)=vy(i,2)+half*(vy(i,5)+vy(i,6)+vy(i,9))
808 vy(i,3)=vy(i,3)+half*(vy(i,6)+vy(i,7)+vy(i,10))
809 vy(i,4)=vy(i,4)+half*(vy(i,8)+vy(i,9)+vy(i,10))
810 vy(i,5) =half*(vy(i,1)+vy(i,2))+fac*vy(i,5)
811 vy(i,6) =half*(vy(i,2)+vy(i,3))+fac*vy(i,6)
812 vy(i,7) =half*(vy(i,1)+vy(i,3))+fac*vy(i,7)
813 vy(i,8) =half*(vy(i,1)+vy(i,4))+fac*vy(i,8)
814 vy(i,9) =half*(vy(i,2)+vy(i,4))+fac*vy(i,9)
815 vy(i,10)=half*(vy(i,3)+vy(i,4))+fac*vy(i,10)
816
817 vz(i,1)=vz(i,1)+half*(vz(i,5)+vz(i,7)+vz(i,8))
818 vz(i,2)=vz(i,2)+half*(vz(i,5)+vz(i,6)+vz(i,9))
819 vz(i,3)=vz(i,3)+half*(vz(i,6)+vz(i,7)+vz(i,10))
820 vz(i,4)=vz(i,4)+half*(vz(i,8)+vz(i,9)+vz(i,10))
821 vz(i,5) =half*(vz(i,1)+vz(i,2))+fac*vz(i,5)
822 vz(i,6) =half*(vz(i,2)+vz(i,3))+fac*vz(i,6)
823 vz(i,7) =half*(vz(i,1)+vz(i,3))+fac*vz(i,7)
824 vz(i,8) =half*(vz(i,1)+vz(i,4))+fac*vz(i,8)
825 vz(i,9) =half*(vz(i,2)+vz(i,4))+fac*vz(i,9)
826 vz(i,10)=half*(vz(i,3)+vz(i,4))+fac*vz(i,10)
827
828 END DO
829
830 nv(1:nel) = zero
831 DO n=1,10
832 DO i=1,nel
833 nv(i) = nv(i) + vx(i,n)*vx(i,n)+vy(i,n)*vy(i,n)+vz(i,n)*vz(i,n)
834 END DO
835 END DO
836 DO i=1,nel
837 nv(i) =sqrt(nv(i))
838 ninv(i)=one/max(em20,nv(i))
839 END DO
840 DO n=1,10
841 DO i=1,nel
842 vx(i,n)=ninv(i) * vx(i,n)
843 vy(i,n)=ninv(i) * vy(i,n)
844 vz(i,n)=ninv(i) * vz(i,n)
845 END DO
846 END DO
847 ux(1:nel,1:10)=vx(1:nel,1:10)
848 uy(1:nel,1:10)=vy(1:nel,1:10)
849 uz(1:nel,1:10)=vz(1:nel,1:10)
850
851 END DO ! DO ITER=1,NITER
852C
853 DO i=1,nel
854 d(i)=nv(i)
855 END DO
856 v_piter(1:nel,1,1:10)=ux(1:nel,1:10)
857 v_piter(1:nel,2,1:10)=uy(1:nel,1:10)
858 v_piter(1:nel,3,1:10)=uz(1:nel,1:10)
859C
860C 2 * Max diagonal as an (under) estimation of Lambda
861C Because Iterative Power may be too much under-estimated before convergence.
862 DO n=1,10
863 DO i=1,nel
864 ul(i,n) = zero
865 ENDDO
866 ENDDO
867C
868 DO ip=1,npt
869C
870C Q = M-1 Transpose(B)
871 DO i=1,nel
872 qx(i,1)=px(i,1,ip)+half*(px(i,5,ip)+px(i,7,ip)+px(i,8,ip))
873 qx(i,2)=px(i,2,ip)+half*(px(i,5,ip)+px(i,6,ip)+px(i,9,ip))
874 qx(i,3)=px(i,3,ip)+half*(px(i,6,ip)+px(i,7,ip)+px(i,10,ip))
875 qx(i,4)=px(i,4,ip)+half*(px(i,8,ip)+px(i,9,ip)+px(i,10,ip))
876C QX(I,5)=HALF*(PX(I,1,IP)+PX(I,2,IP))+(HALF+ONE/FACIROT)*PX(I,5,IP)
877C +FOURTH(PX(I,6,IP)+PX(I,7,IP)+PX(I,8,IP)+PX(I,9,IP)
878 qx(i,5) =half*(qx(i,1)+qx(i,2))+fac*px(i,5,ip)
879 qx(i,6) =half*(qx(i,2)+qx(i,3))+fac*px(i,6,ip)
880 qx(i,7) =half*(qx(i,1)+qx(i,3))+fac*px(i,7,ip)
881 qx(i,8) =half*(qx(i,1)+qx(i,4))+fac*px(i,8,ip)
882 qx(i,9) =half*(qx(i,2)+qx(i,4))+fac*px(i,9,ip)
883 qx(i,10)=half*(qx(i,3)+qx(i,4))+fac*px(i,10,ip)
884
885 qy(i,1)=py(i,1,ip)+half*(py(i,5,ip)+py(i,7,ip)+py(i,8,ip))
886 qy(i,2)=py(i,2,ip)+half*(py(i,5,ip)+py(i,6,ip)+py(i,9,ip))
887 qy(i,3)=py(i,3,ip)+half*(py(i,6,ip)+py(i,7,ip)+py(i,10,ip))
888 qy(i,4)=py(i,4,ip)+half*(py(i,8,ip)+py(i,9,ip)+py(i,10,ip))
889 qy(i,5) =half*(qy(i,1)+qy(i,2))+fac*py(i,5,ip)
890 qy(i,6) =half*(qy(i,2)+qy(i,3))+fac*py(i,6,ip)
891 qy(i,7) =half*(qy(i,1)+qy(i,3))+fac*py(i,7,ip)
892 qy(i,8) =half*(qy(i,1)+qy(i,4))+fac*py(i,8,ip)
893 qy(i,9) =half*(qy(i,2)+qy(i,4))+fac*py(i,9,ip)
894 qy(i,10)=half*(qy(i,3)+qy(i,4))+fac*py(i,10,ip)
895
896 qz(i,1)=pz(i,1,ip)+half*(pz(i,5,ip)+pz(i,7,ip)+pz(i,8,ip))
897 qz(i,2)=pz(i,2,ip)+half*(pz(i,5,ip)+pz(i,6,ip)+pz(i,9,ip))
898 qz(i,3)=pz(i,3,ip)+half*(pz(i,6,ip)+pz(i,7,ip)+pz(i,10,ip))
899 qz(i,4)=pz(i,4,ip)+half*(pz(i,8,ip)+pz(i,9,ip)+pz(i,10,ip))
900 qz(i,5) =half*(qz(i,1)+qz(i,2))+fac*pz(i,5,ip)
901 qz(i,6) =half*(qz(i,2)+qz(i,3))+fac*pz(i,6,ip)
902 qz(i,7) =half*(qz(i,1)+qz(i,3))+fac*pz(i,7,ip)
903 qz(i,8) =half*(qz(i,1)+qz(i,4))+fac*pz(i,8,ip)
904 qz(i,9) =half*(qz(i,2)+qz(i,4))+fac*pz(i,9,ip)
905 qz(i,10)=half*(qz(i,3)+qz(i,4))+fac*pz(i,10,ip)
906
907 END DO
908C
909C Diagonal of M-1 Transpose(B).B
910 DO n=1,10
911 DO i=1,nel
912 ul(i,n) =ul(i,n) + vol(i,ip) *
913 + ( qx(i,n)*px(i,n,ip) + qy(i,n)*py(i,n,ip) + qz(i,n)*pz(i,n,ip) )
914 ENDDO
915 ENDDO
916 ENDDO
917C
918 DO i=1,nel
919 aa = max(ul(i,1),ul(i,2),ul(i,3),ul(i,4))
920 bb = max(ul(i,5),ul(i,6),ul(i,7),ul(i,8),ul(i,9),ul(i,10))
921 d(i) = max(d(i),two*max(aa,bb))
922 ENDDO
923C
924 END IF ! IF(NITER==0)THEN
925
926 DO i=1,nel
927C
928C DELTAX2 is not used w/IDT1TET10
929 deltax2(i) = one
930C
931 mass = volgo(i) ! (multiplied by RHOG0)
932 mref = fourth * mass
933 deltax(i) = two*sqrt(mref/d(i)) !
934
935 ENDDO
936
937 ELSE ! IDT1TET10==0 .OR. ISROT==1
938C
939 IF(isrot == 0)THEN
940C
941 DO n=1,10
942 DO i=1,nel
943 ul(i,n) = zero
944 ENDDO
945 ENDDO
946C
947 DO ip=1,npt
948C
949 DO n=1,10
950 DO i=1,nel
951 ul(i,n) =ul(i,n) + vol(i,ip) *
952 + ( px(i,n,ip)*px(i,n,ip) + py(i,n,ip)*py(i,n,ip)
953 + + pz(i,n,ip)*pz(i,n,ip) )
954 ENDDO
955 ENDDO
956 ENDDO
957C
958 DO i=1,nel
959 aa = max(ul(i,1),ul(i,2),ul(i,3),ul(i,4))
960 bb = max(ul(i,5),ul(i,6),ul(i,7),ul(i,8),ul(i,9),ul(i,10))
961 deltax2(i) = aa/max(aa,bb)
962 aa = aa*thirty2
963 bb = bb*fourty8/seven
964 deltax(i) = sqrt(two*volg(i)/max(aa,bb))
965 ENDDO
966
967 ELSE
968
969 IF(isrot==2.AND.(iint==0.OR.idt1sol>0))THEN
970 DO i=1,nel
971 b1(i) = ty(i)*sz(i) - sy(i)*tz(i)
972 b2(i) = ry(i)*tz(i) - ty(i)*rz(i)
973 b3(i) = sy(i)*rz(i) - ry(i)*sz(i)
974 bb4(i) = -(b1(i) + b2(i) + b3(i))
975C
976 c1(i) = tz(i)*sx(i) - sz(i)*tx(i)
977 c2(i) = rz(i)*tx(i) - tz(i)*rx(i)
978 c3(i) = sz(i)*rx(i) - rz(i)*sx(i)
979 c4(i) = -(c1(i) + c2(i) + c3(i))
980C
981 d1(i) = tx(i)*sy(i) - sx(i)*ty(i)
982 d2(i) = rx(i)*ty(i) - tx(i)*ry(i)
983 d3(i) = sx(i)*ry(i) - rx(i)*sy(i)
984 d4(i) = -(d1(i) + d2(i) + d3(i))
985 det6 = rx(i)*b1(i)+ry(i)*c1(i)+rz(i)*d1(i)
986 dd = one/det6
987 p4x1(i)= b1(i)*dd
988 p4y1(i)= c1(i)*dd
989 p4z1(i)= d1(i)*dd
990 p4x2(i)= b2(i)*dd
991 p4y2(i)= c2(i)*dd
992 p4z2(i)= d2(i)*dd
993 p4x3(i)= b3(i)*dd
994 p4y3(i)= c3(i)*dd
995 p4z3(i)= d3(i)*dd
996 p4x4(i)= bb4(i)*dd
997 p4y4(i)= c4(i)*dd
998 p4z4(i)= d4(i)*dd
999c DELTAX2(I) = ONE
1000 END DO
1001 DO i=1,nel
1002 aa = max(rx(i)*rx(i)+ry(i)*ry(i)+rz(i)*rz(i),
1003 . sx(i)*sx(i)+sy(i)*sy(i)+sz(i)*sz(i),
1004 . tx(i)*tx(i)+ty(i)*ty(i)+tz(i)*tz(i),
1005 . (rx(i)-sx(i))*(rx(i)-sx(i))
1006 + + (ry(i)-sy(i))*(ry(i)-sy(i))
1007 + + (rz(i)-sz(i))*(rz(i)-sz(i)),
1008 . (sx(i)-tx(i))*(sx(i)-tx(i))
1009 + + (sy(i)-ty(i))*(sy(i)-ty(i))
1010 + + (sz(i)-tz(i))*(sz(i)-tz(i)),
1011 . (tx(i)-rx(i))*(tx(i)-rx(i))
1012 + + (ty(i)-ry(i))*(ty(i)-ry(i))
1013 + + (tz(i)-rz(i))*(tz(i)-rz(i)))
1014 deltax2(i) = aa
1015 ENDDO
1016 IF(iformdt==0)THEN
1017 DO i=1,nel
1018 dd = p4x1(i)*p4x1(i)+p4y1(i)*p4y1(i)+p4z1(i)*p4z1(i)
1019 . + p4x2(i)*p4x2(i)+p4y2(i)*p4y2(i)+p4z2(i)*p4z2(i)
1020 . + p4x3(i)*p4x3(i)+p4y3(i)*p4y3(i)+p4z3(i)*p4z3(i)
1021 . + p4x4(i)*p4x4(i)+p4y4(i)*p4y4(i)+p4z4(i)*p4z4(i)
1022 deltax(i) = one / sqrt(dd)
1023 END DO
1024
1025 ELSEIF(iformdt==1)THEN
1026
1027 gfac=pm(105,mxt(1))
1028 ld =two*sqrt(max(one-gfac,zero))+one
1029 DO i=1,nel
1030 pxx(i)=p4x1(i)*p4x1(i)+p4x2(i)*p4x2(i)+p4x3(i)*p4x3(i)+p4x4(i)*p4x4(i)
1031 pyy(i)=p4y1(i)*p4y1(i)+p4y2(i)*p4y2(i)+p4y3(i)*p4y3(i)+p4y4(i)*p4y4(i)
1032 pzz(i)=p4z1(i)*p4z1(i)+p4z2(i)*p4z2(i)+p4z3(i)*p4z3(i)+p4z4(i)*p4z4(i)
1033 pxy(i)=p4x1(i)*p4y1(i)+p4x2(i)*p4y2(i)+p4x3(i)*p4y3(i)+p4x4(i)*p4y4(i)
1034 pxz(i)=p4x1(i)*p4z1(i)+p4x2(i)*p4z2(i)+p4x3(i)*p4z3(i)+p4x4(i)*p4z4(i)
1035 pyz(i)=p4y1(i)*p4z1(i)+p4y2(i)*p4z2(i)+p4y3(i)*p4z3(i)+p4y4(i)*p4z4(i)
1036C
1037 aa = -(pxx(i)+pyy(i)+pzz(i))
1038 bb = (pxx(i)*pyy(i)+pxx(i)*pzz(i)+pyy(i)*pzz(i)-pxy(i)**2-pxz(i)**2-pyz(i)**2)
1039 p = bb-third*aa*aa
1040 dd = two*sqrt(third*max(-p,zero))-third*aa
1041C
1042 dd=ld*dd
1043C
1044 deltax(i) = one / sqrt(dd)
1045 END DO
1046
1047 ELSEIF(iformdt==2)THEN
1048
1049 gfac=pm(105,mxt(1))
1050 DO i=1,nel
1051 pxx(i)=p4x1(i)*p4x1(i)+p4x2(i)*p4x2(i)+p4x3(i)*p4x3(i)+p4x4(i)*p4x4(i)
1052 pyy(i)=p4y1(i)*p4y1(i)+p4y2(i)*p4y2(i)+p4y3(i)*p4y3(i)+p4y4(i)*p4y4(i)
1053 pzz(i)=p4z1(i)*p4z1(i)+p4z2(i)*p4z2(i)+p4z3(i)*p4z3(i)+p4z4(i)*p4z4(i)
1054 pxy(i)=p4x1(i)*p4y1(i)+p4x2(i)*p4y2(i)+p4x3(i)*p4y3(i)+p4x4(i)*p4y4(i)
1055 pxz(i)=p4x1(i)*p4z1(i)+p4x2(i)*p4z2(i)+p4x3(i)*p4z3(i)+p4x4(i)*p4z4(i)
1056 pyz(i)=p4y1(i)*p4z1(i)+p4y2(i)*p4z2(i)+p4y3(i)*p4z3(i)+p4y4(i)*p4z4(i)
1057C
1058 aa = -(pxx(i)+pyy(i)+pzz(i))
1059 bb = gfac*(pxx(i)*pyy(i)+pxx(i)*pzz(i)+pyy(i)*pzz(i)-pxy(i)**2-pxz(i)**2-pyz(i)**2)
1060 p = bb-third*aa*aa
1061 dd = two*sqrt(third*max(-p,zero))-third*aa
1062C
1063 deltax(i) = one / sqrt(dd)
1064 END DO
1065 END if!(IFORMDT==0)THEN
1066 ELSE
1067C
1068 DO i=1,nel
1069
1070 a1x = ry(i)*sz(i)-rz(i)*sy(i)
1071 a1y = rz(i)*sx(i)-rx(i)*sz(i)
1072 a1z = rx(i)*sy(i)-ry(i)*sx(i)
1073 a1 = a1x*a1x+a1y*a1y+a1z*a1z
1074
1075 a2x = sy(i)*tz(i)-sz(i)*ty(i)
1076 a2y = sz(i)*tx(i)-sx(i)*tz(i)
1077 a2z = sx(i)*ty(i)-sy(i)*tx(i)
1078 a2 = a2x*a2x+a2y*a2y+a2z*a2z
1079
1080 a3x = ty(i)*rz(i)-tz(i)*ry(i)
1081 a3y = tz(i)*rx(i)-tx(i)*rz(i)
1082 a3z = tx(i)*ry(i)-ty(i)*rx(i)
1083 a3 = a3x*a3x+a3y*a3y+a3z*a3z
1084
1085 a4x = a1x+a2x+a3x
1086 a4y = a1y+a2y+a3y
1087 a4z = a1z+a2z+a3z
1088 a4 = a4x*a4x+a4y*a4y+a4z*a4z
1089
1090 deltax(i) = six*volg(i)/sqrt(max(a1,a2,a3,a4))
1091cc minoration deltax car inertie constante avec marge
1092cc d'un facteur sqrt(8) sur le cot AA0
1093cc approche a partir du volume initial
1094c AA0 = ((SIX*SQR2*VOLO(I))**TWO_THIRD) * EIGHT
1095c
1096c inertie major e d'un facteur 9.6 (ratio des Lc^2)
1097c sur +grand cot initial
1098c et non plus le cot moyen (SIX*SQR2*VOLO(I))**TWO_THIRD
1099c => il faudrait recalculer ou stocker AA0 avec la m me
1100c formule que AA
1101 aa = max(rx(i)*rx(i)+ry(i)*ry(i)+rz(i)*rz(i),
1102 . sx(i)*sx(i)+sy(i)*sy(i)+sz(i)*sz(i),
1103 . tx(i)*tx(i)+ty(i)*ty(i)+tz(i)*tz(i),
1104 . (rx(i)-sx(i))*(rx(i)-sx(i))
1105 + + (ry(i)-sy(i))*(ry(i)-sy(i))
1106 + + (rz(i)-sz(i))*(rz(i)-sz(i)),
1107 . (sx(i)-tx(i))*(sx(i)-tx(i))
1108 + + (sy(i)-ty(i))*(sy(i)-ty(i))
1109 + + (sz(i)-tz(i))*(sz(i)-tz(i)),
1110 . (tx(i)-rx(i))*(tx(i)-rx(i))
1111 + + (ty(i)-ry(i))*(ty(i)-ry(i))
1112 + + (tz(i)-rz(i))*(tz(i)-rz(i)))
1113c DELTAX(I) = DELTAX(I)*MIN(ONE,AA0/AA)
1114 deltax2(i) = aa
1115 ENDDO
1116 END IF !(IDT1SOL/=0 )THEN
1117 END IF
1118
1119 ENDIF
1120C-----------
1121 1000 FORMAT(/' ZERO OR NEGATIVE VOLUME : 10 NODES TETRAHEDRON NB ',i10,
1122 . ' INTEGRATION POINT NB ',i1/)
1123 1100 FORMAT(/' zero or negative volume : 4 nodes tetrahedron nb ',I10,
1124 . ' integration point nb ',I1/)
1125 2000 FORMAT(/' zero or negative volume : delete 3d-element nb',I10/)
1126 3000 FORMAT(/' zero or negative volume : 3d-element nb:',I10/,
1127 + ' element is switched to small strain option'/)
1128 4000 FORMAT(/' zero or negative volume : 3d-element nb:',I10/)
1129C-----------
1130 RETURN
1131 END
if
if(complex_arithmetic) id
s10len3
subroutine s10len3(vol, ngl, deltax, deltax2, px, py, pz, volg, volgo, rx, ry, rz, sx, sy, sz, tx, ty, tz, nc, nel, mxt, pm, tagelem_sms, v_piter, npt, iint, isrot, iformdt)
Definition s10len3.F:36
max
#define max(a, b)
Definition macros.h:21