OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
spstab.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!|| spstabw ../engine/source/elements/sph/spstab.F
25!||--- called by ------------------------------------------------------
26!|| forintp ../engine/source/elements/forintp.F
27!||--- calls -----------------------------------------------------
28!|| weight0 ../engine/source/elements/sph/weight.F
29!||--- uses -----------------------------------------------------
30!|| sphbox ../engine/share/modules/sphbox.F
31!||====================================================================
32 SUBROUTINE spstabw(
33 1 ITASK ,IPARG ,NGROUNC ,IGROUNC ,KXSP ,
34 2 ISPCOND ,ISPSYM ,WASPACT ,SPH2SOL ,WA ,
35 3 WASIGSM ,WAR ,STAB ,IXSP ,NOD2SP ,
36 4 SPBUF ,X ,IPART ,IPARTSP ,XSPSYM )
37C-----------------------------------------------
38C M o d u l e s
39C-----------------------------------------------
40 USE sphbox
41C-----------------------------------------------
42C I m p l i c i t T y p e s
43C-----------------------------------------------
44#include "implicit_f.inc"
45C-----------------------------------------------
46C C o m m o n B l o c k s
47C-----------------------------------------------
48#include "sphcom.inc"
49#include "param_c.inc"
50#include "scr17_c.inc"
51#include "task_c.inc"
52C-----------------------------------------------
53C D u m m y A r g u m e n t s
54C-----------------------------------------------
55 INTEGER ITASK, IPARG(NPARG,*), NGROUNC, IGROUNC(*), KXSP(NISP,*),
56 . ISPCOND(NISPCOND,*), ISPSYM(NSPCOND,*), WASPACT(*),
57 . SPH2SOL(*), IXSP(KVOISPH,*), NOD2SP(*), IPART(LIPART1,*),
58 . IPARTSP(*)
59 my_real
60 . wa(kwasph,*), wasigsm(6,*), war(10,*), stab(7,*),
61 . spbuf(nspbuf,*), x(3,*), xspsym(3,*)
62C-----------------------------------------------
63 my_real
64 . get_u_geo
65 EXTERNAL get_u_geo
66C-----------------------------------------------
67C L o c a l V a r i a b l e s
68C-----------------------------------------------
69 INTEGER I, J, N, NS, INOD, JNOD, M, NVOIS, NN,
70 . NSPACTF, NSPACTL, IPRT, IGEO, ID, IS,
71 . NVOISS, JS, SM, NC, IACT
72 my_real
73 . xi, yi, zi, xj, yj, zj, di, dmin, dd, wght,
74 . xs, ys, zs, zstab,nul,uno
75C-----------------------------------------------
76 nul=zero
77 uno=one
78 nspactf=1+itask*nsphact/nthread
79 nspactl=(itask+1)*nsphact/nthread
80 DO ns=nspactf,nspactl
81 n=waspact(ns)
82 inod =kxsp(3,n)
83 xi=x(1,inod)
84 yi=x(2,inod)
85 zi=x(3,inod)
86 di =spbuf(1,n)
87 nvois =kxsp(4,n)
88 iprt =ipartsp(n)
89 igeo =ipart(2,iprt)
90 zstab =get_u_geo(7,igeo)
91C
92C indique si il faut calculer SI =-MAX(ZERO,SI)
93C <=> ZSTAB/=0 et il existe voisin a prendre en compte dans spforcp.F
94 stab(7,n)=zero
95 IF(zstab/=zero.AND.nvois/=0)THEN
96 iact=0
97 DO j=1,nvois
98 jnod=ixsp(j,n)
99 IF(jnod>0)THEN
100 m=nod2sp(jnod)
101C
102C Solids to SPH, no interaction if both particles are inactive
103 IF(kxsp(2,n)>0.OR.kxsp(2,m)>0)THEN
104 iact=1
105 EXIT
106 END IF
107 ELSE ! cellule remote
108 nn = -jnod
109C
110C Solidds to SPH, no interaction if both particles are inactive
111 IF(kxsp(2,n)>0.OR.xsphr(13,nn)>0)THEN
112 iact=1
113 EXIT
114 END IF
115 END IF
116 END DO
117 IF(iact==1) THEN
118C-- W(DP)-- DP is the distance of the PPV in initial configuration
119C-- DD = DP/H
120 IF (spbuf(15,n)>em30) THEN
121 dd = spbuf(15,n)
122 ELSE
123 dd=two/three
124 ENDIF
125 CALL weight0(nul,nul,nul,dd,nul,nul,uno,wght)
126 stab(7,n)=zstab/max(em30,wght*wght*wght*wght)
127 ENDIF
128 END IF
129 END DO
130C----------------------------------
131 RETURN
132 END
133!||====================================================================
134!|| spstabs ../engine/source/elements/sph/spstab.F
135!||--- called by ------------------------------------------------------
136!|| forintp ../engine/source/elements/forintp.F
137!||--- calls -----------------------------------------------------
138!|| my_barrier ../engine/source/system/machine.F
139!|| sph_valpvec ../engine/source/elements/sph/spstab.F
140!|| startimeg ../engine/source/system/timer.F
141!|| stoptimeg ../engine/source/system/timer.F
142!||--- uses -----------------------------------------------------
143!|| sphbox ../engine/share/modules/sphbox.F
144!||====================================================================
145 SUBROUTINE spstabs(
146 1 ITASK ,IPARG ,NGROUNC ,IGROUNC ,KXSP ,
147 2 ISPCOND ,ISPSYM ,WASPACT ,SPH2SOL ,WA ,
148 3 WASIGSM ,WAR ,STAB ,IXSP ,NOD2SP ,
149 4 SPBUF ,X )
150C-----------------------------------------------
151C M o d u l e s
152C-----------------------------------------------
153 USE sphbox
154C-----------------------------------------------
155C I m p l i c i t T y p e s
156C-----------------------------------------------
157#include "implicit_f.inc"
158C-----------------------------------------------
159C G l o b a l P a r a m e t e r s
160C-----------------------------------------------
161#include "mvsiz_p.inc"
162C-----------------------------------------------
163C C o m m o n B l o c k s
164C-----------------------------------------------
165#include "vect01_c.inc"
166#include "sphcom.inc"
167#include "param_c.inc"
168#include "task_c.inc"
169C-----------------------------------------------
170C D u m m y A r g u m e n t s
171C-----------------------------------------------
172 INTEGER ITASK, IPARG(NPARG,*), NGROUNC, IGROUNC(*), KXSP(NISP,*),
173 . ISPCOND(NISPCOND,*), ISPSYM(NSPCOND,*), WASPACT(*),
174 . SPH2SOL(*), IXSP(KVOISPH,*), NOD2SP(*)
175C REAL
176 my_real
177 . WA(KWASPH,*), WASIGSM(6,*), WAR(10,*), STAB(7,*),
178 . SPBUF(NSPBUF,*), X(3,*)
179C-----------------------------------------------
180C L o c a l V a r i a b l e s
181C-----------------------------------------------
182 INTEGER I, J, N, NINDX, INDEX(MVSIZ), KS, JS, SM, SS, KR, NR,
183 . IG, NELEM, NG, OFFSET, NEL,
184 . nc, ic, is, islide, nsphrft, nsphrlt, l, SIZE,
185 . ns, inod, jnod, m, nvois,
186 . nspactf, nspactl
187 my_real
188 . sigprv(3,mvsiz), dirprv(3,3,mvsiz),
189 . r1, r2, r3, xi, yi, zi, xj, yj, zj, di, dmin, dd, wght,
190 . xs, ys, zs, sig(6,mvsiz)
191C-----------------------------------------------
192C Boucle parallele dynamique SMP
193C
194!$OMP DO SCHEDULE(DYNAMIC,1)
195c------------------------
196 DO ig = 1, ngrounc
197 ng = igrounc(ig)
198C--------
199C Solids to SPH, particles must be computed when cloud active
200 IF(iparg(8,ng)==1)GOTO 350
201 IF (iddw>0) CALL startimeg(ng)
202 DO nelem = 1,iparg(2,ng),nvsiz
203 offset = nelem - 1
204 nel =iparg(2,ng)
205 nft =iparg(3,ng) + offset
206 iad =iparg(4,ng)
207 ity =iparg(5,ng)
208 isph2sol=iparg(69,ng)
209 ipartsph=0
210 lft=1
211 llt=min(nvsiz,nel-nelem+1)
212 IF(ity==51) THEN
213C-----------
214 nindx=0
215 DO i=lft,llt
216 n =nft+i
217 IF(stab(7,n)/=zero)THEN
218 nindx=nindx+1
219 index(nindx)=i
220 END IF
221 END DO
222C----------------------------------
223C Eigenvalues needed to be calculated in double precision
224C for a simple precision executing
225 size=kwasph
226c IF (IRESP==1) THEN
227c CALL VALPVECDP(AV,SIGPRV,DIRPRV,NEL)
228c ELSE
229 CALL sph_valpvec(nindx,index,SIZE,wa(1,nft+1),sigprv,dirprv)
230c ENDIF
231 DO j=1,nindx
232 i=index(j)
233 n =nft+i
234 r1=-max(zero,sigprv(1,i))
235 r2=-max(zero,sigprv(2,i))
236 r3=-max(zero,sigprv(3,i))
237 stab(1,n) = dirprv(1,1,i)*dirprv(1,1,i)*r1
238 & + dirprv(1,2,i)*dirprv(1,2,i)*r2
239 & + dirprv(1,3,i)*dirprv(1,3,i)*r3
240 stab(2,n) = dirprv(2,2,i)*dirprv(2,2,i)*r2
241 & + dirprv(2,3,i)*dirprv(2,3,i)*r3
242 & + dirprv(2,1,i)*dirprv(2,1,i)*r1
243 stab(3,n) = dirprv(3,3,i)*dirprv(3,3,i)*r3
244 & + dirprv(3,1,i)*dirprv(3,1,i)*r1
245 & + dirprv(3,2,i)*dirprv(3,2,i)*r2
246 stab(4,n) = dirprv(1,1,i)*dirprv(2,1,i)*r1
247 & + dirprv(1,2,i)*dirprv(2,2,i)*r2
248 & + dirprv(1,3,i)*dirprv(2,3,i)*r3
249 stab(5,n) = dirprv(2,2,i)*dirprv(3,2,i)*r2
250 & + dirprv(2,3,i)*dirprv(3,3,i)*r3
251 & + dirprv(2,1,i)*dirprv(3,1,i)*r1
252 stab(6,n) = dirprv(3,3,i)*dirprv(1,3,i)*r3
253 & + dirprv(3,1,i)*dirprv(1,1,i)*r1
254 & + dirprv(3,2,i)*dirprv(1,2,i)*r2
255c STAB(1,N)=-MAX(ZERO,THIRD*(WA(1,N)+WA(2,N)+WA(3,N)))
256c STAB(2,N)=-MAX(ZERO,THIRD*(WA(1,N)+WA(2,N)+WA(3,N)))
257c STAB(3,N)=-MAX(ZERO,THIRD*(WA(1,N)+WA(2,N)+WA(3,N)))
258c STAB(4,N)=ZERO
259c STAB(5,N)=ZERO
260c STAB(6,N)=ZERO
261 END DO
262 ENDIF
263 IF (iddw>0) CALL stoptimeg(ng)
264 END DO
265C--------
266 350 CONTINUE
267 END DO
268!$OMP END DO
269C----------------------------------
270 nsphrft=1+itask*nsphr/nthread
271 nsphrlt=(itask+1)*nsphr/nthread
272 DO n=nsphrft,nsphrlt,mvsiz
273 nindx=0
274 DO l=1,min(nsphrlt-n+1,mvsiz)
275 i=n+l-1
276 IF(stab(7,numsph+i)/=zero)THEN
277 nindx=nindx+1
278 index(nindx)=l
279 END IF
280 ENDDO
281C----------------------------------
282C Eigenvalues needed to be calculated in double precision
283C for a simple precision executing
284 size=10
285c IF (IRESP==1) THEN
286c CALL VALPVECDP(AV,SIGPRV,DIRPRV,NEL)
287c ELSE
288 CALL sph_valpvec(nindx,index,SIZE,war(1,n),sigprv,dirprv)
289c ENDIF
290 DO j=1,nindx
291 i = index(j)
292 nr = n + i - 1
293 kr = numsph + nr
294 r1=-max(zero,sigprv(1,i))
295 r2=-max(zero,sigprv(2,i))
296 r3=-max(zero,sigprv(3,i))
297 stab(1,kr) = dirprv(1,1,i)*dirprv(1,1,i)*r1
298 & + dirprv(1,2,i)*dirprv(1,2,i)*r2
299 & + dirprv(1,3,i)*dirprv(1,3,i)*r3
300 stab(2,kr) = dirprv(2,2,i)*dirprv(2,2,i)*r2
301 & + dirprv(2,3,i)*dirprv(2,3,i)*r3
302 & + dirprv(2,1,i)*dirprv(2,1,i)*r1
303 stab(3,kr) = dirprv(3,3,i)*dirprv(3,3,i)*r3
304 & + dirprv(3,1,i)*dirprv(3,1,i)*r1
305 & + dirprv(3,2,i)*dirprv(3,2,i)*r2
306 stab(4,kr) = dirprv(1,1,i)*dirprv(2,1,i)*r1
307 & + dirprv(1,2,i)*dirprv(2,2,i)*r2
308 & + dirprv(1,3,i)*dirprv(2,3,i)*r3
309 stab(5,kr) = dirprv(2,2,i)*dirprv(3,2,i)*r2
310 & + dirprv(2,3,i)*dirprv(3,3,i)*r3
311 & + dirprv(2,1,i)*dirprv(3,1,i)*r1
312 stab(6,kr) = dirprv(3,3,i)*dirprv(1,3,i)*r3
313 & + dirprv(3,1,i)*dirprv(1,1,i)*r1
314 & + dirprv(3,2,i)*dirprv(1,2,i)*r2
315
316c STAB(1,KR)=-MAX(ZERO,THIRD*(WAR(1,NR)+WAR(2,NR)+WAR(3,NR)))
317c STAB(2,KR)=-MAX(ZERO,THIRD*(WAR(1,NR)+WAR(2,NR)+WAR(3,NR)))
318c STAB(3,KR)=-MAX(ZERO,THIRD*(WAR(1,NR)+WAR(2,NR)+WAR(3,NR)))
319c STAB(4,KR)=ZERO
320c STAB(5,KR)=ZERO
321c STAB(6,KR)=ZERO
322 END DO
323 END DO
324C----------------------------------
325C
326C Particules symetriques
327C
328C /---------------/
329 CALL my_barrier
330C /---------------/
331 nspactf=1+itask*nsphact/nthread
332 nspactl=(itask+1)*nsphact/nthread
333 DO nc=1,nspcond
334 ic=ispcond(2,nc)
335 is=ispcond(3,nc)
336 islide=ispcond(5,nc)
337 IF(islide==0)THEN
338 DO ss=nspactf,nspactl
339 sm=waspact(ss)
340 js=ispsym(nc,sm)
341 ks=numsph+nsphr+js
342 IF(js>0.AND.stab(7,sm)/=zero)THEN
343C nothing changes
344 stab(1,ks)=stab(1,sm)
345 stab(2,ks)=stab(2,sm)
346 stab(3,ks)=stab(3,sm)
347 stab(4,ks)=stab(4,sm)
348 stab(5,ks)=stab(5,sm)
349 stab(6,ks)=stab(6,sm)
350 END IF
351 END DO
352 ELSE !IF(ISLIDE==0)THEN
353 DO is=nspactf,nspactl,mvsiz
354 nindx=0
355 DO l=1,min(nspactl-is+1,mvsiz)
356 ss=is+l-1
357 sm=waspact(ss)
358 js=ispsym(nc,sm)
359 ks=numsph+nsphr+js
360 IF(js>0.AND.stab(7,sm)/=zero)THEN
361 nindx=nindx+1
362 index(nindx)=nindx
363 sig(1,nindx)=wasigsm(1,js)
364 sig(2,nindx)=wasigsm(2,js)
365 sig(3,nindx)=wasigsm(3,js)
366 sig(4,nindx)=wasigsm(4,js)
367 sig(5,nindx)=wasigsm(5,js)
368 sig(6,nindx)=wasigsm(6,js)
369 END IF
370 END DO
371 size=6
372 CALL sph_valpvec(nindx,index,SIZE,sig,sigprv,dirprv)
373 nindx=0
374 DO l=1,min(nspactl-is+1,mvsiz)
375 ss=is+l-1
376 sm=waspact(ss)
377 js=ispsym(nc,sm)
378 ks=numsph+nsphr+js
379 IF(js>0.AND.stab(7,sm)/=zero)THEN
380 nindx=nindx+1
381 r1=-max(zero,sigprv(1,nindx))
382 r2=-max(zero,sigprv(2,nindx))
383 r3=-max(zero,sigprv(3,nindx))
384 stab(1,ks) = dirprv(1,1,nindx)*dirprv(1,1,nindx)*r1
385 & + dirprv(1,2,nindx)*dirprv(1,2,nindx)*r2
386 & + dirprv(1,3,nindx)*dirprv(1,3,nindx)*r3
387 stab(2,ks) = dirprv(2,2,nindx)*dirprv(2,2,nindx)*r2
388 & + dirprv(2,3,nindx)*dirprv(2,3,nindx)*r3
389 & + dirprv(2,1,nindx)*dirprv(2,1,nindx)*r1
390 stab(3,ks) = dirprv(3,3,nindx)*dirprv(3,3,nindx)*r3
391 & + dirprv(3,1,nindx)*dirprv(3,1,nindx)*r1
392 & + dirprv(3,2,nindx)*dirprv(3,2,nindx)*r2
393 stab(4,ks) = dirprv(1,1,nindx)*dirprv(2,1,nindx)*r1
394 & + dirprv(1,2,nindx)*dirprv(2,2,nindx)*r2
395 & + dirprv(1,3,nindx)*dirprv(2,3,nindx)*r3
396 stab(5,ks) = dirprv(2,2,nindx)*dirprv(3,2,nindx)*r2
397 & + dirprv(2,3,nindx)*dirprv(3,3,nindx)*r3
398 & + dirprv(2,1,nindx)*dirprv(3,1,nindx)*r1
399 stab(6,ks) = dirprv(3,3,nindx)*dirprv(1,3,nindx)*r3
400 & + dirprv(3,1,nindx)*dirprv(1,1,nindx)*r1
401 & + dirprv(3,2,nindx)*dirprv(1,2,nindx)*r2
402c STAB(1,KS)=-MAX(ZERO,THIRD*(SIG(1,NINDX)+SIG(2,NINDX)+SIG(3,NINDX)))
403c STAB(2,KS)=-MAX(ZERO,THIRD*(SIG(1,NINDX)+SIG(2,NINDX)+SIG(3,NINDX)))
404c STAB(3,KS)=-MAX(ZERO,THIRD*(SIG(1,NINDX)+SIG(2,NINDX)+SIG(3,NINDX)))
405c STAB(4,KS)=ZERO
406c STAB(5,KS)=ZERO
407c STAB(6,KS)=ZERO
408 END IF
409 END DO
410 END DO
411 END IF
412 END DO
413C-----
414C
415C Particules symetriques de particules remotes
416C
417 DO nc=1,nspcond
418 ic=ispcond(2,nc)
419 is=ispcond(3,nc)
420 islide=ispcond(5,nc)
421 IF(islide==0)THEN
422 DO ss=nsphrft,nsphrlt
423 js=ispsymr(nc,ss)
424 ks=numsph+nsphr+js
425 IF(js>0.AND.stab(7,numsph+ss)/=zero)THEN
426C nothing changes
427 stab(1,ks)=stab(1,numsph+ss)
428 stab(2,ks)=stab(2,numsph+ss)
429 stab(3,ks)=stab(3,numsph+ss)
430 stab(4,ks)=stab(4,numsph+ss)
431 stab(5,ks)=stab(5,numsph+ss)
432 stab(6,ks)=stab(6,numsph+ss)
433 END IF
434 END DO
435 ELSE !IF(ISLIDE==0)THEN
436 DO is=nsphrft,nsphrlt,mvsiz
437 nindx=0
438 DO l=1,min(nsphrlt-is+1,mvsiz)
439 ss=is+l-1
440 js=ispsymr(nc,ss)
441 ks=numsph+nsphr+js
442 IF(js>0.AND.stab(7,numsph+ss)/=zero)THEN
443 nindx=nindx+1
444 index(nindx)=nindx
445 sig(1,nindx)=wasigsm(1,js)
446 sig(2,nindx)=wasigsm(2,js)
447 sig(3,nindx)=wasigsm(3,js)
448 sig(4,nindx)=wasigsm(4,js)
449 sig(5,nindx)=wasigsm(5,js)
450 sig(6,nindx)=wasigsm(6,js)
451 END IF
452 END DO
453 size=6
454 CALL sph_valpvec(nindx,index,SIZE,sig,sigprv,dirprv)
455 nindx=0
456 DO l=1,min(nsphrlt-is+1,mvsiz)
457 ss=is+l-1
458 js=ispsymr(nc,ss)
459 ks=numsph+nsphr+js
460 IF(js>0.AND.stab(7,numsph+ss)/=zero)THEN
461 nindx=nindx+1
462 r1=-max(zero,sigprv(1,nindx))
463 r2=-max(zero,sigprv(2,nindx))
464 r3=-max(zero,sigprv(3,nindx))
465 stab(1,ks) = dirprv(1,1,nindx)*dirprv(1,1,nindx)*r1
466 & + dirprv(1,2,nindx)*dirprv(1,2,nindx)*r2
467 & + dirprv(1,3,nindx)*dirprv(1,3,nindx)*r3
468 stab(2,ks) = dirprv(2,2,nindx)*dirprv(2,2,nindx)*r2
469 & + dirprv(2,3,nindx)*dirprv(2,3,nindx)*r3
470 & + dirprv(2,1,nindx)*dirprv(2,1,nindx)*r1
471 stab(3,ks) = dirprv(3,3,nindx)*dirprv(3,3,nindx)*r3
472 & + dirprv(3,1,nindx)*dirprv(3,1,nindx)*r1
473 & + dirprv(3,2,nindx)*dirprv(3,2,nindx)*r2
474 stab(4,ks) = dirprv(1,1,nindx)*dirprv(2,1,nindx)*r1
475 & + dirprv(1,2,nindx)*dirprv(2,2,nindx)*r2
476 & + dirprv(1,3,nindx)*dirprv(2,3,nindx)*r3
477 stab(5,ks) = dirprv(2,2,nindx)*dirprv(3,2,nindx)*r2
478 & + dirprv(2,3,nindx)*dirprv(3,3,nindx)*r3
479 & + dirprv(2,1,nindx)*dirprv(3,1,nindx)*r1
480 stab(6,ks) = dirprv(3,3,nindx)*dirprv(1,3,nindx)*r3
481 & + dirprv(3,1,nindx)*dirprv(1,1,nindx)*r1
482 & + dirprv(3,2,nindx)*dirprv(1,2,nindx)*r2
483c STAB(1,KS)=-MAX(ZERO,THIRD*(SIG(1,NINDX)+SIG(2,NINDX)+SIG(3,NINDX)))
484c STAB(2,KS)=-MAX(ZERO,THIRD*(SIG(1,NINDX)+SIG(2,NINDX)+SIG(3,NINDX)))
485c STAB(3,KS)=-MAX(ZERO,THIRD*(SIG(1,NINDX)+SIG(2,NINDX)+SIG(3,NINDX)))
486c STAB(4,KS)=ZERO
487c STAB(5,KS)=ZERO
488c STAB(6,KS)=ZERO
489 END IF
490 END DO
491 END DO
492 END IF
493 END DO
494C----------------------------------
495 RETURN
496 END
497C
498!||====================================================================
499!|| sph_valpvec ../engine/source/elements/sph/spstab.F
500!||--- called by ------------------------------------------------------
501!|| spstabs ../engine/source/elements/sph/spstab.F
502!||--- calls -----------------------------------------------------
503!|| floatmin ../common_source/tools/math/precision.c
504!||====================================================================
505 SUBROUTINE sph_valpvec(NINDX,INDEX,SIZE,SIG,VAL,VEC)
506C-----------------------------------------------
507C I m p l i c i t T y p e s
508C-----------------------------------------------
509#include "implicit_f.inc"
510C-----------------------------------------------
511C G l o b a l P a r a m e t e r s
512C-----------------------------------------------
513#include "mvsiz_p.inc"
514C-----------------------------------------------
515C D u m m y A r g u m e n t s
516C-----------------------------------------------
517 INTEGER NINDX, INDEX(*), SIZE
518 my_real
519 . SIG(SIZE,*), VAL(3,*), VEC(9,*)
520C-----------------------------------------------
521C L o c a l V a r i a b l e s
522C-----------------------------------------------
523 INTEGER I, L, N, II, NN, LMAX, LMAXV(MVSIZ),
524 . NINDEX1, NINDEX2, NINDEX3,
525 . index1(mvsiz), index2(mvsiz), index3(mvsiz)
526 my_real
527 . cs(6), str(3,mvsiz), a(3,3,mvsiz), v(3,3,mvsiz), b(3,3,mvsiz),
528 . xmag(3,mvsiz), pr, aa, bb, aaa(mvsiz),
529 . cc, angp, dd, ftpi, ttpi, strmax,
530 . tol1(mvsiz), tol2(mvsiz), xmaxv(mvsiz), norminf(mvsiz),
531 . vmag, s11,
532 . s21, s31, s12, s22, s32, s13, s23, s33, a11, a12, a13, a21,
533 . a22, a23, a31, a32, a33,
534 . mdemi, xmaxinv, flm
535 REAL FLMIN
536C-----------------------------------------------
537C DATA FTPI,TTPI / 4.188790205, 2.094395102 /
538C FTPI=(4/3)*PI, TTPI=(2/3)*PI
539C
540C DEVIATEUR PRINCIPAL DE CONTRAINTE
541C . . . . . . . . . . . . . . . . . . .
542 MDEMI = -half
543 ttpi = acos(mdemi)
544 ftpi = two*ttpi
545C precision minimum dependant du type REAL ou DOUBLE
546 CALL floatmin(cs(1),cs(2),flmin)
547 flm = two*sqrt(flmin)
548 nindex3=0
549 DO ii = 1, nindx
550 nn=index(ii)
551 cs(1) = sig(1,nn)
552 cs(2) = sig(2,nn)
553 cs(3) = sig(3,nn)
554 cs(4) = sig(4,nn)
555 cs(5) = sig(5,nn)
556 cs(6) = sig(6,nn)
557 pr = -(cs(1)+cs(2)+cs(3))* third
558 cs(1) = cs(1) + pr
559 cs(2) = cs(2) + pr
560 cs(3) = cs(3) + pr
561 aaa(nn)=cs(4)**2 + cs(5)**2 + cs(6)**2 - cs(1)*cs(2)
562 & - cs(2)*cs(3) - cs(1)*cs(3)
563 norminf(nn) = max(abs(cs(1)),abs(cs(2)),abs(cs(3)),
564 & abs(cs(4)),abs(cs(5)),abs(cs(6)))
565 norminf(nn) = em10*norminf(nn)
566C cas racine triple
567c AA = MAX(AAA(NN),NORMINF(NN),EM20)
568 aa = max(aaa(nn),norminf(nn))
569C
570 bb=cs(1)*cs(5)**2 + cs(2)*cs(6)**2
571 & + cs(3)*cs(4)**2 - cs(1)*cs(2)*cs(3)
572 & - two*cs(4)*cs(5)*cs(6)
573C
574 cc=-sqrt(twenty7/max(em20,aa))*bb*half/max(em20,aa)
575 cc= min(cc,one)
576 cc= max(cc,-one)
577 angp=acos(cc) * third
578 dd=two*sqrt(aa * third)
579 str(1,nn)=dd*cos(angp)
580 str(2,nn)=dd*cos(angp+ftpi)
581 str(3,nn)=dd*cos(angp+ttpi)
582C
583 val(1,nn) = str(1,nn)-pr
584 val(2,nn) = str(2,nn)-pr
585 val(3,nn) = str(3,nn)-pr
586C renforcement de precision en compression----
587 IF(abs(str(3,nn))>abs(str(1,nn))
588 & .AND.aaa(nn)>norminf(nn)) THEN
589 aa=str(1,nn)
590 str(1,nn)=str(3,nn)
591 str(3,nn)=aa
592 nindex3 = nindex3+1
593 index3(nindex3) = nn
594 ENDIF
595C . . . . . . . . . . .
596C VECTEURS PROPRES
597C . . . . . . . . . . .
598 strmax= max(abs(str(1,nn)),abs(str(3,nn)))
599 tol1(nn)= max(em20,flm*strmax**2)
600 tol2(nn)=flm*strmax/3
601 a(1,1,nn)=cs(1)-str(1,nn)
602 a(2,2,nn)=cs(2)-str(1,nn)
603 a(3,3,nn)=cs(3)-str(1,nn)
604 a(1,2,nn)=cs(4)
605 a(2,1,nn)=cs(4)
606 a(2,3,nn)=cs(5)
607 a(3,2,nn)=cs(5)
608 a(1,3,nn)=cs(6)
609 a(3,1,nn)=cs(6)
610C
611 b(1,1,nn)=a(2,1,nn)*a(3,2,nn)-a(3,1,nn)
612 . *a(2,2,nn)
613 b(1,2,nn)=a(2,2,nn)*a(3,3,nn)-a(3,2,nn)
614 . *a(2,3,nn)
615 b(1,3,nn)=a(2,3,nn)*a(3,1,nn)-a(3,3,nn)
616 . *a(2,1,nn)
617 b(2,1,nn)=a(3,1,nn)*a(1,2,nn)-a(1,1,nn)
618 . *a(3,2,nn)
619 b(2,2,nn)=a(3,2,nn)*a(1,3,nn)-a(1,2,nn)
620 . *a(3,3,nn)
621 b(2,3,nn)=a(3,3,nn)*a(1,1,nn)-a(1,3,nn)
622 . *a(3,1,nn)
623 b(3,1,nn)=a(1,1,nn)*a(2,2,nn)-a(2,1,nn)
624 . *a(1,2,nn)
625 b(3,2,nn)=a(1,2,nn)*a(2,3,nn)-a(2,2,nn)
626 . *a(1,3,nn)
627 b(3,3,nn)=a(1,3,nn)*a(2,1,nn)-a(2,3,nn)
628 . *a(1,1,nn)
629 xmag(1,nn)=sqrt(b(1,1,nn)**2+b(2,1,nn)**2+b(3,1,nn)**2)
630 xmag(2,nn)=sqrt(b(1,2,nn)**2+b(2,2,nn)**2+b(3,2,nn)**2)
631 xmag(3,nn)=sqrt(b(1,3,nn)**2+b(2,3,nn)**2+b(3,3,nn)**2)
632
633 ENDDO
634C
635 nindex1 = 0
636 nindex2 = 0
637 DO ii = 1, nindx
638 nn=index(ii)
639 xmaxv(nn)=max(xmag(1,nn),xmag(2,nn),xmag(3,nn))
640 IF(xmag(1,nn)==xmaxv(nn)) THEN
641 lmaxv(nn) = 1
642 ELSEIF(xmag(2,nn)==xmaxv(nn)) THEN
643 lmaxv(nn) = 2
644 ELSE
645 lmaxv(nn) = 3
646 ENDIF
647 IF(aaa(nn)<norminf(nn)) THEN
648 val(1,nn) = sig(1,nn)
649 val(2,nn) = sig(2,nn)
650 val(3,nn) = sig(3,nn)
651 v(1,1,nn) = one
652 v(2,1,nn) = zero
653 v(3,1,nn) = zero
654 v(1,2,nn) = zero
655 v(2,2,nn) = one
656 v(3,2,nn) = zero
657
658 ELSEIF(xmaxv(nn)>tol1(nn)) THEN
659 nindex1 = nindex1 + 1
660 index1(nindex1) = nn
661 ELSE
662 nindex2 = nindex2 + 1
663 index2(nindex2) = nn
664 ENDIF
665 ENDDO
666C
667#include "vectorize.inc"
668 DO n = 1, nindex1
669 nn = index1(n)
670 lmax = lmaxv(nn)
671 xmaxinv = one/xmaxv(nn)
672 v(1,1,nn)=b(1,lmax,nn)*xmaxinv
673 v(2,1,nn)=b(2,lmax,nn)*xmaxinv
674 v(3,1,nn)=b(3,lmax,nn)*xmaxinv
675 a(1,1,nn)=a(1,1,nn)+str(1,nn)-str(3,nn)
676 a(2,2,nn)=a(2,2,nn)+str(1,nn)-str(3,nn)
677 a(3,3,nn)=a(3,3,nn)+str(1,nn)-str(3,nn)
678C
679 b(1,1,nn)=a(2,1,nn)*v(3,1,nn)-a(3,1,nn)*v(2,1,nn)
680 b(1,2,nn)=a(2,2,nn)*v(3,1,nn)-a(3,2,nn)*v(2,1,nn)
681 b(1,3,nn)=a(2,3,nn)*v(3,1,nn)-a(3,3,nn)*v(2,1,nn)
682 b(2,1,nn)=a(3,1,nn)*v(1,1,nn)-a(1,1,nn)*v(3,1,nn)
683 b(2,2,nn)=a(3,2,nn)*v(1,1,nn)-a(1,2,nn)*v(3,1,nn)
684 b(2,3,nn)=a(3,3,nn)*v(1,1,nn)-a(1,3,nn)*v(3,1,nn)
685 b(3,1,nn)=a(1,1,nn)*v(2,1,nn)-a(2,1,nn)*v(1,1,nn)
686 b(3,2,nn)=a(1,2,nn)*v(2,1,nn)-a(2,2,nn)*v(1,1,nn)
687 b(3,3,nn)=a(1,3,nn)*v(2,1,nn)-a(2,3,nn)*v(1,1,nn)
688 xmag(1,nn)=sqrt(b(1,1,nn)**2+b(2,1,nn)**2+b(3,1,nn)**2)
689 xmag(2,nn)=sqrt(b(1,2,nn)**2+b(2,2,nn)**2+b(3,2,nn)**2)
690 xmag(3,nn)=sqrt(b(1,3,nn)**2+b(2,3,nn)**2+b(3,3,nn)**2)
691C
692 xmaxv(nn)=max(xmag(1,nn),xmag(2,nn),xmag(3,nn))
693 ENDDO
694C
695#include "vectorize.inc"
696 DO n = 1, nindex1
697 nn = index1(n)
698 IF(xmag(1,nn)==xmaxv(nn)) THEN
699 lmaxv(nn) = 1
700 ELSEIF(xmag(2,nn)==xmaxv(nn)) THEN
701 lmaxv(nn) = 2
702 ELSE
703 lmaxv(nn) = 3
704 ENDIF
705C
706 vmag=sqrt(v(1,1,nn)**2+v(2,1,nn)**2)
707 lmax = lmaxv(nn)
708 IF(xmaxv(nn)>tol2(nn))THEN
709 xmaxinv = one/xmaxv(nn)
710 v(1,3,nn)=b(1,lmax,nn)*xmaxinv
711 v(2,3,nn)=b(2,lmax,nn)*xmaxinv
712 v(3,3,nn)=b(3,lmax,nn)*xmaxinv
713 v(1,2,nn)=v(2,3,nn)*v(3,1,nn)-v(2,1,nn)*v(3,3,nn)
714 v(2,2,nn)=v(3,3,nn)*v(1,1,nn)-v(3,1,nn)*v(1,3,nn)
715 v(3,2,nn)=v(1,3,nn)*v(2,1,nn)-v(1,1,nn)*v(2,3,nn)
716 vmag=one/sqrt(v(1,2,nn)**2+v(2,2,nn)**2+v(3,2,nn)**2)
717 v(1,2,nn)=v(1,2,nn)*vmag
718 v(2,2,nn)=v(2,2,nn)*vmag
719 v(3,2,nn)=v(3,2,nn)*vmag
720 ELSEIF(vmag>tol2(nn))THEN
721 v(1,2,nn)=-v(2,1,nn)/vmag
722 v(2,2,nn)=v(1,1,nn)/vmag
723 v(3,2,nn)=zero
724 ELSE
725 v(1,2,nn)=one
726 v(2,2,nn)=zero
727 v(3,2,nn)=zero
728 ENDIF
729 ENDDO
730C . . . . . . . . . . . . .
731C SOLUTION DOUBLE
732C . . . . . . . . . . . . .
733 DO n = 1, nindex2
734 nn = index2(n)
735 xmag(1,nn)=sqrt(a(1,1,nn)**2+a(2,1,nn)**2)
736 xmag(2,nn)=sqrt(a(1,2,nn)**2+a(2,2,nn)**2)
737 xmag(3,nn)=sqrt(a(1,3,nn)**2+a(2,3,nn)**2)
738C
739 xmaxv(nn)=max(xmag(1,nn),xmag(2,nn),xmag(3,nn))
740 ENDDO
741C
742#include "vectorize.inc"
743 DO n = 1, nindex2
744 nn = index2(n)
745 IF(xmag(1,nn)==xmaxv(nn)) THEN
746 lmaxv(nn) = 1
747 ELSEIF(xmag(2,nn)==xmaxv(nn)) THEN
748 lmaxv(nn) = 2
749 ELSE
750 lmaxv(nn) = 3
751 ENDIF
752C
753 lmax = lmaxv(nn)
754 IF(max(abs(a(3,1,nn)),abs(a(3,2,nn)),abs(a(3,3,nn)))
755 & <tol2(nn))THEN
756 xmaxinv = one/xmaxv(nn)
757 v(1,1,nn)= zero
758 v(2,1,nn)= zero
759 v(3,1,nn)= one
760 v(1,2,nn)=-a(2,lmax,nn)*xmaxinv
761 v(2,2,nn)= a(1,lmax,nn)*xmaxinv
762 v(3,2,nn)= zero
763C
764 ELSEIF(xmaxv(nn)>tol2(nn))THEN
765 xmaxinv = one/xmaxv(nn)
766 v(1,1,nn)=-a(2,lmax,nn)*xmaxinv
767 v(2,1,nn)= a(1,lmax,nn)*xmaxinv
768 v(3,1,nn)= zero
769 v(1,2,nn)=-a(3,lmax,nn)*v(2,1,nn)
770 v(2,2,nn)= a(3,lmax,nn)*v(1,1,nn)
771 v(3,2,nn)= a(1,lmax,nn)*v(2,1,nn)-a(2,lmax,nn)*v(1,1,nn)
772 vmag=one/sqrt(v(1,2,nn)**2+v(2,2,nn)**2+v(3,2,nn)**2)
773 v(1,2,nn)=v(1,2,nn)*vmag
774 v(2,2,nn)=v(2,2,nn)*vmag
775 v(3,2,nn)=v(3,2,nn)*vmag
776 ELSE
777 v(1,1,nn) = one
778 v(2,1,nn) = zero
779 v(3,1,nn) = zero
780 v(1,2,nn) = zero
781 v(2,2,nn) = one
782 v(3,2,nn) = zero
783 ENDIF
784 ENDDO
785C
786 DO ii = 1, nindx
787 nn=index(ii)
788 vec(1,nn)=v(1,1,nn)
789 vec(2,nn)=v(2,1,nn)
790 vec(3,nn)=v(3,1,nn)
791 vec(4,nn)=v(1,2,nn)
792 vec(5,nn)=v(2,2,nn)
793 vec(6,nn)=v(3,2,nn)
794 vec(7,nn)=vec(2,nn)*vec(6,nn)-vec(3,nn)*vec(5,nn)
795 vec(8,nn)=vec(3,nn)*vec(4,nn)-vec(1,nn)*vec(6,nn)
796 vec(9,nn)=vec(1,nn)*vec(5,nn)-vec(2,nn)*vec(4,nn)
797 ENDDO
798C reecriture pour contourner probleme sur itanium2 comp 9. + latency=16
799 DO n = 1, nindex3
800 nn = index3(n)
801C str utilise com tableau temporaire au lieu de scalaires temporaires
802 str(1,nn)=vec(7,nn)
803 str(2,nn)=vec(8,nn)
804 str(3,nn)=vec(9,nn)
805 ENDDO
806 DO n = 1, nindex3
807 nn = index3(n)
808 vec(7,nn)=vec(1,nn)
809 vec(8,nn)=vec(2,nn)
810 vec(9,nn)=vec(3,nn)
811 vec(1,nn)=-str(1,nn)
812 vec(2,nn)=-str(2,nn)
813 vec(3,nn)=-str(3,nn)
814 ENDDO
815C
816 RETURN
817 END
818
my_real
#define my_real
Definition cppsort.cpp:32
startimeg
subroutine startimeg(ng)
Definition timer.F:1487
stoptimeg
subroutine stoptimeg(ng)
Definition timer.F:1535
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
sphbox::ispsymr
integer, dimension(:,:), allocatable ispsymr
Definition sphbox.F:93
sphbox::nsphr
integer nsphr
Definition sphbox.F:83
floatmin
void floatmin(int *a, int *b, float *flm)
Definition precision.c:71
spstabw
subroutine spstabw(itask, iparg, ngrounc, igrounc, kxsp, ispcond, ispsym, waspact, sph2sol, wa, wasigsm, war, stab, ixsp, nod2sp, spbuf, x, ipart, ipartsp, xspsym)
Definition spstab.F:37
spstabs
subroutine spstabs(itask, iparg, ngrounc, igrounc, kxsp, ispcond, ispsym, waspact, sph2sol, wa, wasigsm, war, stab, ixsp, nod2sp, spbuf, x)
Definition spstab.F:150
sph_valpvec
subroutine sph_valpvec(nindx, index, size, sig, val, vec)
Definition spstab.F:506
my_barrier
subroutine my_barrier
Definition machine.F:31
weight0
subroutine weight0(xi, yi, zi, xj, yj, zj, h, w)
Definition weight.F:34