OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
s10len3.F File Reference
#include "implicit_f.inc"
#include "mvsiz_p.inc"
#include "vect01_c.inc"
#include "param_c.inc"
#include "scr17_c.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine s10len3 (vol, ngl, deltax, deltax2, px, py, pz, volu, voln, volg, rx, ry, rz, sx, sy, sz, tx, ty, tz, nel, mxt, pm, v_piter, iint)

Function/Subroutine Documentation

◆ s10len3()

subroutine s10len3 ( vol,
integer, dimension(*) ngl,
deltax,
deltax2,
px,
py,
pz,
volu,
voln,
volg,
rx,
ry,
rz,
sx,
sy,
sz,
tx,
ty,
tz,
integer nel,
integer, dimension(mvsiz) mxt,
pm,
v_piter,
integer iint )

Definition at line 29 of file s10len3.F.

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