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

Go to the source code of this file.

Functions/Subroutines

subroutine spstabw (itask, iparg, ngrounc, igrounc, kxsp, ispcond, ispsym, waspact, sph2sol, wa, wasigsm, war, stab, ixsp, nod2sp, spbuf, x, ipart, ipartsp, xspsym)
subroutine spstabs (itask, iparg, ngrounc, igrounc, kxsp, ispcond, ispsym, waspact, sph2sol, wa, wasigsm, war, stab, ixsp, nod2sp, spbuf, x)
subroutine sph_valpvec (nindx, index, size, sig, val, vec)

Function/Subroutine Documentation

◆ sph_valpvec()

subroutine sph_valpvec ( integer nindx,
integer, dimension(*) index,
integer size,
sig,
val,
vec )

Definition at line 505 of file spstab.F.

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
my_real
#define my_real
Definition cppsort.cpp:32
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
floatmin
void floatmin(int *a, int *b, float *flm)
Definition precision.c:71

◆ spstabs()

subroutine spstabs ( integer itask,
integer, dimension(nparg,*) iparg,
integer ngrounc,
integer, dimension(*) igrounc,
integer, dimension(nisp,*) kxsp,
integer, dimension(nispcond,*) ispcond,
integer, dimension(nspcond,*) ispsym,
integer, dimension(*) waspact,
integer, dimension(*) sph2sol,
wa,
wasigsm,
war,
stab,
integer, dimension(kvoisph,*) ixsp,
integer, dimension(*) nod2sp,
spbuf,
x )

Definition at line 145 of file spstab.F.

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
startimeg
subroutine startimeg(ng)
Definition timer.F:1487
stoptimeg
subroutine stoptimeg(ng)
Definition timer.F:1535
sphbox::ispsymr
integer, dimension(:,:), allocatable ispsymr
Definition sphbox.F:93
sphbox::nsphr
integer nsphr
Definition sphbox.F:83
sph_valpvec
subroutine sph_valpvec(nindx, index, size, sig, val, vec)
Definition spstab.F:506
my_barrier
subroutine my_barrier
Definition machine.F:31

◆ spstabw()

subroutine spstabw ( integer itask,
integer, dimension(nparg,*) iparg,
integer ngrounc,
integer, dimension(*) igrounc,
integer, dimension(nisp,*) kxsp,
integer, dimension(nispcond,*) ispcond,
integer, dimension(nspcond,*) ispsym,
integer, dimension(*) waspact,
integer, dimension(*) sph2sol,
wa,
wasigsm,
war,
stab,
integer, dimension(kvoisph,*) ixsp,
integer, dimension(*) nod2sp,
spbuf,
x,
integer, dimension(lipart1,*) ipart,
integer, dimension(*) ipartsp,
xspsym )

Definition at line 32 of file spstab.F.

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
weight0
subroutine weight0(xi, yi, zi, xj, yj, zj, h, w)
Definition weight.F:34