OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps33.F File Reference
#include "implicit_f.inc"
#include "scr05_c.inc"
#include "mvsiz_p.inc"
#include "com01_c.inc"
#include "vectorize.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine sigeps33 (nel, nuparam, nuvar, nfunc, ifunc, npf, tf, time, timestep, uparam, rho0, rho, volume, eint, epspxx, epspyy, epspzz, epspxy, epspyz, epspzx, depsxx, depsyy, depszz, depsxy, depsyz, depszx, epsxx, epsyy, epszz, epsxy, epsyz, epszx, sigoxx, sigoyy, sigozz, sigoxy, sigoyz, sigozx, signxx, signyy, signzz, signxy, signyz, signzx, sigvxx, sigvyy, sigvzz, sigvxy, sigvyz, sigvzx, soundsp, viscmax, uvar, off)
subroutine jacobiew (a, n, ew, ev, nrot)
subroutine dreh (b, tr, ii, jj, ken)
subroutine valpvec (sig, val, vec, nel)
subroutine valpvecdp (sig, val, vec, nel)
subroutine epsf2u (nel, fxx, fxy, fxz, fyx, fyy, fyz, fzx, fzy, fzz, uxx, uyy, uzz, uxy, uyz, uxz)
subroutine valpvecop (sig, val, vec, nel)
subroutine valpvec_v (sig, val, vec, nel)
subroutine valpvecdp_v (sig, val, vec, nel)
subroutine valpvecop_v (sig, val, vec, nel)

Function/Subroutine Documentation

◆ dreh()

subroutine dreh ( b,
tr,
integer ii,
integer jj,
integer ken )

Definition at line 738 of file sigeps33.F.

739C-------------------------------------------------------KK010587-
740C KEN=0 -> ROTATION OF MATRIX B BY TR = TR(T)*B*TR
741C KEN=1 -> ROTATION OF MATRIX B BY TR(T) = TR*B*TR(T)
742C II AND JJ POINT TO THE SUBMATRIX TO BE ROTATED; DEFAULT II=JJ=1
743C----------------------------------------------------------------
744C-----------------------------------------------
745C I m p l i c i t T y p e s
746C-----------------------------------------------
747#include "implicit_f.inc"
748 my_real
749 . b(3,3),tr(3,3),lk(3,3),x
750 INTEGER KEN,II,JJ
751 INTEGER I,J,K,I1,J1
752C----------------------------------------------------------------
753 IF(ii<=0) ii=1
754 IF(jj<=0) jj=1
755
756 DO 20 i=1,3
757 DO 20 j=1,3
758 j1=jj+j-1
759 x=0.0
760 DO 10 k=1,3
761C IF(TR(K,I)==ZERO) GOTO 10
762 i1=k+ii-1
763C IF(B(I1,J1)==ZERO) GOTO 10
764 IF(ken/=1) x=x+tr(k,i)*b(i1,j1)
765 IF(ken==1) x=x+tr(i,k)*b(i1,j1)
766 10 CONTINUE
767 20 lk(i,j)=x
768
769 DO 40 i=1,3
770 DO 40 j=1,3
771 x=zero
772 DO 30 k=1,3
773C IF(TR(K,J)==ZERO) GOTO 30
774C IF(LK(I,K)==ZERO) GOTO 30
775 IF(ken/=1) x=x+lk(i,k)*tr(k,j)
776 IF(ken==1) x=x+lk(i,k)*tr(j,k)
777 30 CONTINUE
778 40 b(ii+i-1,jj+j-1)=x
779
780 RETURN
781
my_real
#define my_real
Definition cppsort.cpp:32

◆ epsf2u()

subroutine epsf2u ( integer nel,
fxx,
fxy,
fxz,
fyx,
fyy,
fyz,
fzx,
fzy,
fzz,
uxx,
uyy,
uzz,
uxy,
uyz,
uxz )

Definition at line 1443 of file sigeps33.F.

1448C-----------------------------------------------
1449C I M P L I C I T T Y P E S
1450C-----------------------------------------------
1451#include "implicit_f.inc"
1452C-----------------------------------------------
1453C G L O B A L P A R A M E T E R S
1454C-----------------------------------------------
1455#include "mvsiz_p.inc"
1456C-----------------------------------------------
1457C C O M M O N
1458C-----------------------------------------------
1459#include "scr05_c.inc"
1460C----------------------------------------------------------------
1461C I N P U T A R G U M E N T S
1462C----------------------------------------------------------------
1463 INTEGER NEL
1464 my_real
1465 . fxx(nel) , fxy(nel), fxz(nel),
1466 . fyx(nel) , fyy(nel), fyz(nel),
1467 . fzx(nel) , fzy(nel), fzz(nel)
1468C----------------------------------------------------------------
1469C O U T P U T A R G U M E N T S
1470C----------------------------------------------------------------
1471 my_real
1472 . uxx(nel) , uxy(nel), uxz(nel),
1473 . uyy(nel) , uyz(nel), uzz(nel)
1474C----------------------------------------------------------------
1475C I N P U T O U T P U T A R G U M E N T S
1476C----------------------------------------------------------------
1477C----------------------------------------------------------------
1478C L O C A L V A R I B L E S
1479C----------------------------------------------------------------
1480 INTEGER I,J,K
1481 my_real
1482 . ev(3,mvsiz),av(6,mvsiz),evv(3,mvsiz),dirprv(3,3,mvsiz)
1483C----------------------------------------------------------------
1484C--------- [C]=[F]^t[F] strain-----
1485 DO i=1,nel
1486 av(1,i)=fxx(i)*fxx(i)+fyx(i)*fyx(i)+fzx(i)*fzx(i)
1487 av(2,i)=fyy(i)*fyy(i)+fxy(i)*fxy(i)+fzy(i)*fzy(i)
1488 av(3,i)=fzz(i)*fzz(i)+fzx(i)*fzx(i)+fzy(i)*fzy(i)
1489 av(4,i)=fxx(i)*fxy(i)+fyx(i)*fyy(i)+fzx(i)*fzy(i)
1490 av(6,i)=fxx(i)*fxz(i)+fyx(i)*fyz(i)+fzx(i)*fzz(i)
1491 av(5,i)=fxz(i)*fxy(i)+fyz(i)*fyy(i)+fzz(i)*fzy(i)
1492 ENDDO
1493C for a simple precision executing
1494 IF (iresp==1) THEN
1495 CALL valpvecdp(av,evv,dirprv,nel)
1496 ELSE
1497 CALL valpvec(av,evv,dirprv,nel)
1498 ENDIF
1499 DO i=1,nel
1500 ev(1,i)=sqrt(evv(1,i))
1501 ev(2,i)=sqrt(evv(2,i))
1502 ev(3,i)=sqrt(evv(3,i))
1503 ENDDO
1504* TRANSFORM PRINCIPAL U TO GLOBAL DIRECTIONS
1505 DO i=1,nel
1506C
1507 uxx(i) = dirprv(1,1,i)*dirprv(1,1,i)*ev(1,i)
1508 . + dirprv(1,2,i)*dirprv(1,2,i)*ev(2,i)
1509 . + dirprv(1,3,i)*dirprv(1,3,i)*ev(3,i)
1510c
1511 uyy(i) = dirprv(2,2,i)*dirprv(2,2,i)*ev(2,i)
1512 . + dirprv(2,3,i)*dirprv(2,3,i)*ev(3,i)
1513 . + dirprv(2,1,i)*dirprv(2,1,i)*ev(1,i)
1514c
1515 uzz(i) = dirprv(3,3,i)*dirprv(3,3,i)*ev(3,i)
1516 . + dirprv(3,1,i)*dirprv(3,1,i)*ev(1,i)
1517 . + dirprv(3,2,i)*dirprv(3,2,i)*ev(2,i)
1518c
1519 uxy(i) = dirprv(1,1,i)*dirprv(2,1,i)*ev(1,i)
1520 . + dirprv(1,2,i)*dirprv(2,2,i)*ev(2,i)
1521 . + dirprv(1,3,i)*dirprv(2,3,i)*ev(3,i)
1522c
1523 uyz(i) = dirprv(2,2,i)*dirprv(3,2,i)*ev(2,i)
1524 . + dirprv(2,3,i)*dirprv(3,3,i)*ev(3,i)
1525 . + dirprv(2,1,i)*dirprv(3,1,i)*ev(1,i)
1526c
1527 uxz(i) = dirprv(3,3,i)*dirprv(1,3,i)*ev(3,i)
1528 . + dirprv(3,1,i)*dirprv(1,1,i)*ev(1,i)
1529 . + dirprv(3,2,i)*dirprv(1,2,i)*ev(2,i)
1530 ENDDO
1531C
1532 RETURN
valpvecdp
subroutine valpvecdp(sig, val, vec, nel)
Definition sigeps33.F:1112
valpvec
subroutine valpvec(sig, val, vec, nel)
Definition sigeps33.F:794

◆ jacobiew()

subroutine jacobiew ( a,
integer n,
ew,
ev,
integer nrot )

Definition at line 606 of file sigeps33.F.

607C-------------------------------------------------------KK141189-
608C COMPUTATION OF ALL EIGENVALUES AND EIGENVECTORS OF A SYMMETRIC
609C MATRIX A BY THE JACOBI ALGORITHM
610C
611C A(N,N) EIGENWERTPROBLEM
612C N DIMENSION OF A
613C EW(N) EIGENVALUES
614C EV(N,N) EIGENVEKTORS
615C NROT NUMBER OF ROTATIONS
616C MAXA MAXIMUM ELEMENT OF A
617C-----------------------------------------------
618C I m p l i c i t T y p e s
619C-----------------------------------------------
620#include "implicit_f.inc"
621 INTEGER NN
622 parameter(nn=9)
623
624 INTEGER N,NROT
625 my_real
626 . a(n,n), ew(n), ev(n,n)
627 . , b(nn), z(nn)
628 INTEGER IZ,IS,ITER,J
629 my_real
630 . sumrs,eps,g,h,t,c,s,tau,theta
631C----------------------------------------------------------------
632 DO 130 iz=1,n
633 DO 120 is=1,n
634C PRACTICAL FOR RADIOSS (PASS ONLY LOWER DIAGONAL MATRIX)
635 IF(iz>is) a(is,iz) = a(iz,is)
636 ev(iz,is)=zero
637 120 CONTINUE
638 b(iz)=a(iz,iz)
639 ew(iz)=b(iz)
640 z(iz)=0.
641 ev(iz,iz)=one
642 130 CONTINUE
643
644 nrot=0
645
646C START ITERATION
647
648 DO 240 iter = 1,50
649 sumrs = zero
650
651C SUM OF THE OFF DIAGONALS
652 DO 150 iz=1,n-1
653 DO 140 is=iz+1,n
654 sumrs=sumrs+abs(a(iz,is))
655 140 CONTINUE
656 150 CONTINUE
657
658 IF (sumrs ==zero) GOTO 9000
659 IF (iter > 4) THEN
660 eps = zero
661 ELSE
662 eps = one_fifth*sumrs/n**2
663 ENDIF
664
665 DO 220 iz=1,n-1
666 DO 210 is=iz+1,n
667 g = 100. * abs(a(iz,is))
668 IF (iter>4 .AND. abs(ew(iz))+g==abs(ew(iz))
669 & .AND. abs(ew(is))+g==abs(ew(is))) THEN
670 a(iz,is)=zero
671 ELSE IF (abs(a(iz,is)) > eps) THEN
672 h = ew(is)-ew(iz)
673 IF (abs(h)+g==abs(h)) THEN
674 t = a(iz,is)/h
675 ELSE
676 theta = half*h/a(iz,is)
677 t=one/(abs(theta)+sqrt(one+theta**2))
678 IF (theta < zero) t=-t
679 ENDIF
680 c=one/sqrt(one+t**2)
681 s=t*c
682 tau=s/(one+c)
683 h=t*a(iz,is)
684 z(iz)=z(iz)-h
685 z(is)=z(is)+h
686 ew(iz)=ew(iz)-h
687 ew(is)=ew(is)+h
688 a(iz,is)=zero
689 DO 160 j=1,iz-1
690 g=a(j,iz)
691 h=a(j,is)
692 a(j,iz)=g-s*(h+g*tau)
693 a(j,is)=h+s*(g-h*tau)
694 160 CONTINUE
695 DO 170 j=iz+1,is-1
696 g=a(iz,j)
697 h=a(j,is)
698 a(iz,j)=g-s*(h+g*tau)
699 a(j,is)=h+s*(g-h*tau)
700 170 CONTINUE
701 DO 180 j=is+1,n
702 g=a(iz,j)
703 h=a(is,j)
704 a(iz,j)=g-s*(h+g*tau)
705 a(is,j)=h+s*(g-h*tau)
706 180 CONTINUE
707 DO 190 j=1,n
708 g=ev(j,iz)
709 h=ev(j,is)
710 ev(j,iz)=g-s*(h+g*tau)
711 ev(j,is)=h+s*(g-h*tau)
712 190 CONTINUE
713 nrot=nrot+1
714 ENDIF
715 210 CONTINUE
716 220 CONTINUE
717
718 DO 230 iz=1,n
719 b(iz)=b(iz)+z(iz)
720 ew(iz)=b(iz)
721 z(iz)=zero
722 230 CONTINUE
723 240 CONTINUE
724
725 9000 CONTINUE
726
727 RETURN
728

◆ sigeps33()

subroutine sigeps33 ( integer nel,
integer nuparam,
integer nuvar,
integer nfunc,
integer, dimension(nfunc) ifunc,
integer, dimension(*) npf,
tf,
time,
timestep,
uparam,
rho0,
rho,
volume,
eint,
epspxx,
epspyy,
epspzz,
epspxy,
epspyz,
epspzx,
depsxx,
depsyy,
depszz,
depsxy,
depsyz,
depszx,
epsxx,
epsyy,
epszz,
epsxy,
epsyz,
epszx,
sigoxx,
sigoyy,
sigozz,
sigoxy,
sigoyz,
sigozx,
signxx,
signyy,
signzz,
signxy,
signyz,
signzx,
sigvxx,
sigvyy,
sigvzz,
sigvxy,
sigvyz,
sigvzx,
soundsp,
viscmax,
uvar,
off )

Definition at line 35 of file sigeps33.F.

46
47C-----------------------------------------------
48C I m p l i c i t T y p e s
49C-----------------------------------------------
50#include "implicit_f.inc"
51#include "scr05_c.inc"
52C-----------------------------------------------
53C G l o b a l P a r a m e t e r s
54C-----------------------------------------------
55#include "mvsiz_p.inc"
56#include "com01_c.inc"
57C----------------------------------------------------------------
58C I N P U T A R G U M E N T S
59C----------------------------------------------------------------
60 INTEGER NEL, NUPARAM, NUVAR
61 my_real
62 . time , timestep , uparam(nuparam),
63 . rho(nel), rho0(nel), volume(nel), eint(nel),
64 . epspxx(nel), epspyy(nel), epspzz(nel),
65 . epspxy(nel), epspyz(nel), epspzx(nel),
66 . depsxx(nel), depsyy(nel), depszz(nel),
67 . depsxy(nel), depsyz(nel), depszx(nel),
68 . epsxx(nel), epsyy(nel), epszz(nel),
69 . epsxy(nel), epsyz(nel), epszx(nel),
70 . sigoxx(nel), sigoyy(nel), sigozz(nel),
71 . sigoxy(nel), sigoyz(nel), sigozx(nel)
72C----------------------------------------------------------------
73C O U T P U T A R G U M E N T S
74C----------------------------------------------------------------
75 my_real
76 . signxx(nel), signyy(nel), signzz(nel),
77 . signxy(nel), signyz(nel), signzx(nel),
78 . sigvxx(nel), sigvyy(nel), sigvzz(nel),
79 . sigvxy(nel), sigvyz(nel), sigvzx(nel),
80 . soundsp(nel), viscmax(nel)
81C----------------------------------------------------------------
82C I N P U T O U T P U T A R G U M E N T S
83C----------------------------------------------------------------
84 my_real uvar(nel,nuvar), off(nel)
85C----------------------------------------------------------------
86C VARIABLES FOR FUNCTION INTERPOLATION
87C----------------------------------------------------------------
88 INTEGER NPF(*), NFUNC, IFUNC(NFUNC)
89 my_real finter,tf(*)
90 EXTERNAL finter
91C----------------------------------------------------------------
92C L O C A L V A R I B L E S
93C----------------------------------------------------------------
94 INTEGER I,J,K,L,NROT,KEN
95 INTEGER VPVEC,ICASE
96 my_real
97 . ey,c1,c2,et,vmu,vmu0
98 . , a,b,c,d,var
99 . , axx,bxx,cxx,ayy,byy,cyy,azz,bzz,czz
100 . , axy,bxy,cxy,ayz,byz,cyz,azx,bzx,czx
101 . , c1xx,c2xx,etxx,vmuxx, c1yy,c2yy,etyy,vmuyy
102 . , c1xy,c2xy,gtxy,vmuxy, c1yz,c2yz,gtyz,vmuyz
103 . , c1zx,c2zx,gtzx,vmuzx, p0,phi,gama0,fac,fac1
104 . , s(3,3),sigpr(3),dirpr(3,3)
105 . , c1zz,c2zz,etzz,vmuzz,sigt_cutoff,sigc_cutoff
106 my_real
107 . e(mvsiz),edot(mvsiz),dav,e1,e2,e3,e4,e5,e6,epsp
108 . , epet(mvsiz),emet(mvsiz),epets(mvsiz),emets(mvsiz)
109 . , emxx(mvsiz),emyy(mvsiz),emzz(mvsiz)
110 . , gmxy(mvsiz),gmyz(mvsiz),gmzx(mvsiz)
111 . , epetxx(mvsiz),epetyy(mvsiz),epetzz(mvsiz)
112 . , emetxx(mvsiz),emetyy(mvsiz),emetzz(mvsiz)
113 . , gpgtxy(mvsiz),gpgtyz(mvsiz),gpgtzx(mvsiz)
114 . , gmgtxy(mvsiz),gmgtyz(mvsiz),gmgtzx(mvsiz)
115 . , syxx(mvsiz),syyy(mvsiz),syzz(mvsiz)
116 . , syxy(mvsiz),syyz(mvsiz),syzx(mvsiz)
117 . , sigair(mvsiz),gama(mvsiz),syield(mvsiz)
118 . , sigsxx(mvsiz),sigsyy(mvsiz),sigszz(mvsiz)
119 . , sigsxy(mvsiz),sigsyz(mvsiz),sigszx(mvsiz)
120 . , dexx(mvsiz),deyy(mvsiz),dezz(mvsiz)
121 . , dexy(mvsiz),deyz(mvsiz),dezx(mvsiz)
122 . , dedtxx(mvsiz),dedtyy(mvsiz),dedtzz(mvsiz)
123 . , dedtxy(mvsiz),dedtyz(mvsiz),dedtzx(mvsiz)
124 . , dsdtxx(mvsiz),dsdtyy(mvsiz),dsdtzz(mvsiz)
125 . , dsdtxy(mvsiz),dsdtyz(mvsiz),dsdtzx(mvsiz)
126C
127 my_real
128 . sigv(mvsiz,6), sigprv(mvsiz,3), dirprv(mvsiz,3,3)
129C
130
131C----------------------------------------------------------------
132
133C INITIALIZATION
134 IF(isigi == 0) THEN
135 IF(time==zero) THEN
136 DO 1000 j=1,nuvar
137 DO 1000 i=1,nel
138 uvar(i,j)=zero
139 1000 CONTINUE
140 ENDIF
141 ENDIF
142
143 ken = uparam(1)
144 ey = uparam(2)
145 a = uparam(3)
146 b = uparam(4)
147 c = uparam(5)
148 p0 = uparam(6)
149 phi = uparam(7)
150 gama0 = uparam(8)
151 !!
152 vpvec = 0
153C COMPUTE VOLUMETRIC STRAIN-AIR PRESSURE
154 DO 900 i=1,nel
155 gama(i) = rho0(i)/rho(i)-one+gama0
156 var = -(p0*gama(i))/(one+gama(i)-phi+em15)
157 sigair(i) = max(zero,var)
158 900 CONTINUE
159 icase = abs(ken) + 1
160 SELECT CASE (icase)
161C----------------------------------
162 CASE(1)
163C----------------------------------
164 fac = uparam(9)
165 fac1 = uparam(10)
166
167 DO 100 i=1,nel
168 IF(ifunc(1)/=0) THEN
169 syield(i)=fac*finter(ifunc(1),gama(i),npf,tf,b*c)
170 ELSE
171 syield(i) = abs(a+b*(one+c*gama(i)))
172 ENDIF
173 IF(ifunc(2)/=0)THEN
174 dav = (epspxx(i)+epspyy(i)+epspzz(i))*third
175 e1 = epspxx(i) - dav
176 e2 = epspyy(i) - dav
177 e3 = epspzz(i) - dav
178 e4 = half*epspxy(i)
179 e5 = half*epspyz(i)
180 e6 = half*epspzx(i)
181 epsp =half*(e1**2+e2**2+e3**2) +e4**2+e5**2+e6**2
182 epsp = sqrt(three*epsp)/three_half
183 syield(i) = fac1*syield(i)*(finter(ifunc(2),epsp,npf,tf,b*c))
184 ENDIF
185 100 CONTINUE
186
187 DO 110 i=1,nel
188 sigsxx(i)=sigoxx(i)+sigair(i)
189 sigsyy(i)=sigoyy(i)+sigair(i)
190 sigszz(i)=sigozz(i)+sigair(i)
191 sigsxy(i)=sigoxy(i)
192 sigsyz(i)=sigoyz(i)
193 sigszx(i)=sigozx(i)
194 110 CONTINUE
195
196C COMP UTE TRIAL STRESSES
197 DO 120 i=1,nel
198 sigsxx(i)=sigsxx(i)+ey*epspxx(i)*timestep
199 sigsyy(i)=sigsyy(i)+ey*epspyy(i)*timestep
200 sigszz(i)=sigszz(i)+ey*epspzz(i)*timestep
201 sigsxy(i)=sigsxy(i)+ey*epspxy(i)*timestep*half
202 sigsyz(i)=sigsyz(i)+ey*epspyz(i)*timestep*half
203 sigszx(i)=sigszx(i)+ey*epspzx(i)*timestep*half
204 120 CONTINUE
205C si ancienne methode alors ....
206 IF (vpvec==1) THEN
207 DO 130 i=1,nel
208 s(1,1)=sigsxx(i)
209 s(2,1)=sigsxy(i)
210 s(2,2)=sigsyy(i)
211 s(3,1)=sigszx(i)
212 s(3,2)=sigsyz(i)
213 s(3,3)=sigszz(i)
214
215C TRANSFORM GLOBAL STRESSES (S_IJ) TO PRINCIPAL STRESSES (SIGPR_K)
216 CALL jacobiew(s,3,sigpr,dirpr,nrot)
217
218C YIELD CRITERIA - SCALING
219 sigpr(1)=min(syield(i),abs(sigpr(1)))*sign(one,sigpr(1))
220 sigpr(2)=min(syield(i),abs(sigpr(2)))*sign(one,sigpr(2))
221 sigpr(3)=min(syield(i),abs(sigpr(3)))*sign(one,sigpr(3))
222
223C TRANSFORM PRINCIPAL STRESSES (SIGMA_K) TO GLOBAL STRESSES (SIGMA_IJ)
224 DO 131 k=1,3
225 DO 131 l=1,3
226 s(k,l)=zero
227 IF(k==l) s(k,l)= sigpr(k)
228 131 CONTINUE
229 CALL dreh(s,dirpr,1,1,1)
230
231 sigsxx(i)=s(1,1)
232 sigsxy(i)=s(2,1)
233 sigsyy(i)=s(2,2)
234 sigszx(i)=s(3,1)
235 sigsyz(i)=s(3,2)
236 sigszz(i)=s(3,3)
237 130 CONTINUE
238C fin ancienne methode
239C sino n
240 ELSE
241 DO i = 1, nel
242 sigv(i,1) = sigsxx(i)
243 sigv(i,2) = sigsyy(i)
244 sigv(i,3) = sigszz(i)
245 sigv(i,4) = sigsxy(i)
246 sigv(i,5) = sigsyz(i)
247 sigv(i,6) = sigszx(i)
248 ENDDO
249C for a simple precision executing
250 IF (iresp==1) THEN
251 CALL valpvecdp_v(sigv,sigprv,dirprv,nel)
252 ELSE
253 CALL valpvec_v(sigv,sigprv,dirprv,nel)
254 ENDIF
255 DO i = 1, nel
256C YIEL D CRITERIA - SCALING
257 sigprv(i,1)=sign(min(syield(i),abs(sigprv(i,1))),sigprv(i,1))
258 sigprv(i,2)=sign(min(syield(i),abs(sigprv(i,2))),sigprv(i,2))
259 sigprv(i,3)=sign(min(syield(i),abs(sigprv(i,3))),sigprv(i,3))
260C TRANSFORM PRINCIPAL STRESSES (SIGMA_K) TO GLOBAL STRESSES (SIGMA_IJ)
261 sigsxx(i) = dirprv(i,1,1)*dirprv(i,1,1)*sigprv(i,1)
262 . + dirprv(i,1,2)*dirprv(i,1,2)*sigprv(i,2)
263 . + dirprv(i,1,3)*dirprv(i,1,3)*sigprv(i,3)
264 sigsyy(i) = dirprv(i,2,2)*dirprv(i,2,2)*sigprv(i,2)
265 . + dirprv(i,2,3)*dirprv(i,2,3)*sigprv(i,3)
266 . + dirprv(i,2,1)*dirprv(i,2,1)*sigprv(i,1)
267 sigszz(i) = dirprv(i,3,3)*dirprv(i,3,3)*sigprv(i,3)
268 . + dirprv(i,3,1)*dirprv(i,3,1)*sigprv(i,1)
269 . + dirprv(i,3,2)*dirprv(i,3,2)*sigprv(i,2)
270 sigsxy(i) = dirprv(i,1,1)*dirprv(i,2,1)*sigprv(i,1)
271 . + dirprv(i,1,2)*dirprv(i,2,2)*sigprv(i,2)
272 . + dirprv(i,1,3)*dirprv(i,2,3)*sigprv(i,3)
273 sigsyz(i) = dirprv(i,2,2)*dirprv(i,3,2)*sigprv(i,2)
274 . + dirprv(i,2,3)*dirprv(i,3,3)*sigprv(i,3)
275 . + dirprv(i,2,1)*dirprv(i,3,1)*sigprv(i,1)
276 sigszx(i) = dirprv(i,3,3)*dirprv(i,1,3)*sigprv(i,3)
277 . + dirprv(i,3,1)*dirprv(i,1,1)*sigprv(i,1)
278 . + dirprv(i,3,2)*dirprv(i,1,2)*sigprv(i,2)
279 ENDDO
280 ENDIF
281C
282C COMPUTE UPDATED STRESSES AND SOUND SPEED
283 DO 140 i=1,nel
284 signxx(i)=off(i)*sigsxx(i)-sigair(i)
285 signyy(i)=off(i)*sigsyy(i)-sigair(i)
286 signzz(i)=off(i)*sigszz(i)-sigair(i)
287 signxy(i)=off(i)*sigsxy(i)
288 signyz(i)=off(i)*sigsyz(i)
289 signzx(i)=off(i)*sigszx(i)
290
291 soundsp(i)=sqrt(ey/rho0(i))
292 140 CONTINUE
293
294C--------------------------
295 CASE(2)
296C--------------------------
297 c1 = uparam(9)
298 c2 = uparam(10)
299 et = uparam(11)
300 vmu = uparam(12)
301 vmu0 = uparam(13)
302 fac = uparam(14)
303 fac1 = uparam(15)
304
305C COMPUTE YOUNG MODULUS, ETC.
306 DO 200 i=1,nel
307 edot(i)=max(
308 & abs(epspxx(i)),abs(epspyy(i)),abs(epspzz(i)),
309 & abs(epspxy(i)),abs(epspyz(i)),abs(epspzx(i)))
310 e(i)=c1*edot(i)+c2
311 e(i)=max(e(i),ey)
312 epet(i)=(e(i)+et)/vmu
313 emet(i)=(e(i)*et)/vmu
314 epets(i)=(e(i)+et)/vmu0
315 emets(i)=two*(e(i)*et)/vmu0
316 200 CONTINUE
317
318C COMPUTE YIELD STRESS
319 DO 210 i=1,nel
320 IF(ifunc(1)/=0) THEN
321 syield(i)=fac*finter(ifunc(1),gama(i),npf,tf,b*c)
322 ELSE
323 syield(i) = abs(a+b*(one+c*gama(i)))
324 ENDIF
325 IF(ifunc(2)/=0)THEN
326 dav = (epspxx(i)+epspyy(i)+epspzz(i))*third
327 e1 = epspxx(i) - dav
328 e2 = epspyy(i) - dav
329 e3 = epspzz(i) - dav
330 e4 = half*epspxy(i)
331 e5 = half*epspyz(i)
332 e6 = half*epspzx(i)
333 epsp =half*(e1**2+e2**2+e3**2) +e4**2+e5**2+e6**2
334 epsp = sqrt(three*epsp)/three_half
335 syield(i) = fac1*syield(i)*(finter(ifunc(2),epsp,npf,tf,b*c))
336 ENDIF
337 210 CONTINUE
338
339 DO 220 i=1,nel
340 sigsxx(i)=sigoxx(i)+sigair(i)
341 sigsyy(i)=sigoyy(i)+sigair(i)
342 sigszz(i)=sigozz(i)+sigair(i)
343 sigsxy(i)=sigoxy(i)
344 sigsyz(i)=sigoyz(i)
345 sigszx(i)=sigozx(i)
346
347 dexx(i)=epsxx(i)
348 deyy(i)=epsyy(i)
349 dezz(i)=epszz(i)
350 dexy(i)=epsxy(i)* half
351 deyz(i)=epsyz(i)* half
352 dezx(i)=epszx(i)* half
353
354 dedtxx(i)=epspxx(i)
355 dedtyy(i)=epspyy(i)
356 dedtzz(i)=epspzz(i)
357 dedtxy(i)=epspxy(i) * half
358 dedtyz(i)=epspyz(i) * half
359 dedtzx(i)=epspzx(i) * half
360 220 CONTINUE
361C COMPUTE STRESS RATES
362 DO 230 i=1,nel
363 dsdtxx(i)=e(i)*dedtxx(i)-epet(i)*sigsxx(i)+emet(i)*dexx(i)
364 dsdtyy(i)=e(i)*dedtyy(i)-epet(i)*sigsyy(i)+emet(i)*deyy(i)
365 dsdtzz(i)=e(i)*dedtzz(i)-epet(i)*sigszz(i)+emet(i)*dezz(i)
366 dsdtxy(i)=e(i)*dedtxy(i)-epets(i)*sigsxy(i)+emets(i)*dexy(i)
367 dsdtyz(i)=e(i)*dedtyz(i)-epets(i)*sigsyz(i)+emets(i)*deyz(i)
368 dsdtzx(i)=e(i)*dedtzx(i)-epets(i)*sigszx(i)+emets(i)*dezx(i)
369 230 CONTINUE
370
371C COMPUTE TRIAL STRESSES
372 DO 240 i=1,nel
373 sigsxx(i)=sigsxx(i)+dsdtxx(i)*timestep
374 sigsyy(i)=sigsyy(i)+dsdtyy(i)*timestep
375 sigszz(i)=sigszz(i)+dsdtzz(i)*timestep
376 sigsxy(i)=sigsxy(i)+dsdtxy(i)*timestep
377 sigsyz(i)=sigsyz(i)+dsdtyz(i)*timestep
378 sigszx(i)=sigszx(i)+dsdtzx(i)*timestep
379 240 CONTINUE
380 IF(ken<0)GOTO 255
381C
382C si ancienne methode alors ....
383 IF (vpvec==1) THEN
384 DO 250 i=1,nel
385 s(1,1)=sigsxx(i)
386 s(2,1)=sigsxy(i)
387 s(2,2)=sigsyy(i)
388 s(3,1)=sigszx(i)
389 s(3,2)=sigsyz(i)
390 s(3,3)=sigszz(i)
391
392C TRANSFORM GLOBAL STRESSES (S_IJ) TO PRINCIPAL STRESSES (SIGPR_K)
393 CALL jacobiew(s,3,sigpr,dirpr,nrot)
394
395C YIELD CRITERIA - SCALING
396 sigpr(1)=min(syield(i),abs(sigpr(1)))*sign(one,sigpr(1))
397 sigpr(2)=min(syield(i),abs(sigpr(2)))*sign(one,sigpr(2))
398 sigpr(3)=min(syield(i),abs(sigpr(3)))*sign(one,sigpr(3))
399
400C TRANSFORM PRINCIPAL STRESSES (SIGMA_K) TO GLOBAL STRESSES (SIGMA_IJ)
401 DO 251 k=1,3
402 DO 251 l=1,3
403 s(k,l)=zero
404 IF(k==l) s(k,l)= sigpr(k)
405 251 CONTINUE
406 CALL dreh(s,dirpr,1,1,1)
407
408 sigsxx(i)=s(1,1)
409 sigsxy(i)=s(2,1)
410 sigsyy(i)=s(2,2)
411 sigszx(i)=s(3,1)
412 sigsyz(i)=s(3,2)
413 sigszz(i)=s(3,3)
414 250 CONTINUE
415 ELSE
416C fin ancienne methode
417 DO i = 1, nel
418 sigv(i,1) = sigsxx(i)
419 sigv(i,2) = sigsyy(i)
420 sigv(i,3) = sigszz(i)
421 sigv(i,4) = sigsxy(i)
422 sigv(i,5) = sigsyz(i)
423 sigv(i,6) = sigszx(i)
424 ENDDO
425 CALL valpvec_v(sigv,sigprv,dirprv,nel)
426 DO i = 1, nel
427C YIELD CRITERIA - SCALING
428 sigprv(i,1)=sign(min(syield(i),abs(sigprv(i,1))),sigprv(i,1))
429 sigprv(i,2)=sign(min(syield(i),abs(sigprv(i,2))),sigprv(i,2))
430 sigprv(i,3)=sign(min(syield(i),abs(sigprv(i,3))),sigprv(i,3))
431C TRANSFORM PRINCIPAL STRESSES (SIGMA_K) TO GLOBAL STRESSES (SIGMA_IJ)
432 sigsxx(i) = dirprv(i,1,1)*dirprv(i,1,1)*sigprv(i,1)
433 . + dirprv(i,1,2)*dirprv(i,1,2)*sigprv(i,2)
434 . + dirprv(i,1,3)*dirprv(i,1,3)*sigprv(i,3)
435 sigsyy(i) = dirprv(i,2,2)*dirprv(i,2,2)*sigprv(i,2)
436 . + dirprv(i,2,3)*dirprv(i,2,3)*sigprv(i,3)
437 . + dirprv(i,2,1)*dirprv(i,2,1)*sigprv(i,1)
438 sigszz(i) = dirprv(i,3,3)*dirprv(i,3,3)*sigprv(i,3)
439 . + dirprv(i,3,1)*dirprv(i,3,1)*sigprv(i,1)
440 . + dirprv(i,3,2)*dirprv(i,3,2)*sigprv(i,2)
441 sigsxy(i) = dirprv(i,1,1)*dirprv(i,2,1)*sigprv(i,1)
442 . + dirprv(i,1,2)*dirprv(i,2,2)*sigprv(i,2)
443 . + dirprv(i,1,3)*dirprv(i,2,3)*sigprv(i,3)
444 sigsyz(i) = dirprv(i,2,2)*dirprv(i,3,2)*sigprv(i,2)
445 . + dirprv(i,2,3)*dirprv(i,3,3)*sigprv(i,3)
446 . + dirprv(i,2,1)*dirprv(i,3,1)*sigprv(i,1)
447 sigszx(i) = dirprv(i,3,3)*dirprv(i,1,3)*sigprv(i,3)
448 . + dirprv(i,3,1)*dirprv(i,1,1)*sigprv(i,1)
449 . + dirprv(i,3,2)*dirprv(i,1,2)*sigprv(i,2)
450 ENDDO
451 ENDIF
452C fin chgt
453 255 CONTINUE
454C COMPUTE UPDATED STRESSES AND SOUND SPEED
455 DO 260 i=1,nel
456 signxx(i)=off(i)*sigsxx(i)-sigair(i)
457 signyy(i)=off(i)*sigsyy(i)-sigair(i)
458 signzz(i)=off(i)*sigszz(i)-sigair(i)
459 signxy(i)=off(i)*sigsxy(i)
460 signyz(i)=off(i)*sigsyz(i)
461 signzx(i)=off(i)*sigszx(i)
462
463 soundsp(i)=sqrt(e(i)/rho0(i))
464 260 CONTINUE
465C--------------------------
466 CASE(3)
467C--------------------------
468 fac = uparam(9)
469 fac1 = uparam(10)
470 sigt_cutoff = uparam(11)
471 sigc_cutoff = sigt_cutoff
472 DO i=1,nel
473 IF(ifunc(1)/=0) THEN
474 syield(i)=fac*finter(ifunc(1),gama(i),npf,tf,b*c)
475 ELSE
476 syield(i) = abs(a+b*(one+c*gama(i)))
477 ENDIF
478 IF(gama(i) > zero .AND. (gama(i)-uvar(i,1)) < zero) sigc_cutoff = zero
479 IF(ifunc(2)/=0)THEN
480 dav = (epspxx(i)+epspyy(i)+epspzz(i))*third
481 e1 = epspxx(i) - dav
482 e2 = epspyy(i) - dav
483 e3 = epspzz(i) - dav
484 e4 = half*epspxy(i)
485 e5 = half*epspyz(i)
486 e6 = half*epspzx(i)
487 epsp =half*(e1**2+e2**2+e3**2) +e4**2+e5**2+e6**2
488 epsp = sqrt(three*epsp)/three_half
489 syield(i) = fac1*syield(i)*(finter(ifunc(2),epsp,npf,tf,b*c))
490 ENDIF
491 ENDDO
492
493 DO i=1,nel
494 sigsxx(i)=sigoxx(i)+sigair(i)
495 sigsyy(i)=sigoyy(i)+sigair(i)
496 sigszz(i)=sigozz(i)+sigair(i)
497 sigsxy(i)=sigoxy(i)
498 sigsyz(i)=sigoyz(i)
499 sigszx(i)=sigozx(i)
500 ENDDO
501
502C COMPUTE TRIAL STRESSES
503 DO i=1,nel
504 sigsxx(i)=sigsxx(i)+ey*epspxx(i)*timestep
505 sigsyy(i)=sigsyy(i)+ey*epspyy(i)*timestep
506 sigszz(i)=sigszz(i)+ey*epspzz(i)*timestep
507 sigsxy(i)=sigsxy(i)+ey*epspxy(i)*timestep*half
508 sigsyz(i)=sigsyz(i)+ey*epspyz(i)*timestep*half
509 sigszx(i)=sigszx(i)+ey*epspzx(i)*timestep*half
510 ENDDO
511 DO i = 1, nel
512 sigv(i,1) = sigsxx(i)
513 sigv(i,2) = sigsyy(i)
514 sigv(i,3) = sigszz(i)
515 sigv(i,4) = sigsxy(i)
516 sigv(i,5) = sigsyz(i)
517 sigv(i,6) = sigszx(i)
518 ENDDO
519C for a simple precision executing
520 IF (iresp==1) THEN
521 CALL valpvecdp_v(sigv,sigprv,dirprv,nel)
522 ELSE
523 CALL valpvec_v(sigv,sigprv,dirprv,nel)
524 ENDIF
525 DO i =1, nel
526C YIELD CRITERIA - SCALING
527 IF(abs(gama(i)) < em10) THEN
528 sigprv(i,1)=sign(min(abs(sigprv(i,1)), sigt_cutoff),sigprv(i,1))
529 sigprv(i,2)=sign(min(abs(sigprv(i,2)), sigt_cutoff),sigprv(i,2))
530 sigprv(i,3)=sign(min(abs(sigprv(i,3)), sigt_cutoff),sigprv(i,3))
531 ELSEIF(gama(i) < zero )then
532 !!IF(DELTA_GAMA(I) < ZERO)then
533 sigprv(i,1)=sign(min(syield(i),abs(sigprv(i,1))),sigprv(i,1))
534 sigprv(i,2)=sign(min(syield(i),abs(sigprv(i,2))),sigprv(i,2))
535 sigprv(i,3)=sign(min(syield(i),abs(sigprv(i,3))),sigprv(i,3))
536 !! ELSE
537 sigprv(i,1)=min(sigprv(i,1), sigt_cutoff)
538 sigprv(i,2)=min(sigprv(i,2), sigt_cutoff)
539 sigprv(i,3)=min(sigprv(i,3), sigt_cutoff)
540 !! ENDIF
541 ELSE
542 sigprv(i,1)=sign(min(abs(sigprv(i,1)), sigt_cutoff),sigprv(i,1))
543 sigprv(i,2)=sign(min(abs(sigprv(i,2)), sigt_cutoff),sigprv(i,2))
544 sigprv(i,3)=sign(min(abs(sigprv(i,3)), sigt_cutoff),sigprv(i,3))
545 !!SIGPRV(I,1)=MIN(SIGPRV(I,1), SIGT_CUTOFF)
546 !!SIGPRV(I,2)=MIN(SIGPRV(I,2), SIGT_CUTOFF)
547 !!SIGPRV(I,3)=MIN(SIGPRV(I,3), SIGT_CUTOFF)
548 !! IF(DELTA_GAMA(I) < ZERO) THEN
549 sigprv(i,1)=max(sigprv(i,1), zero)
550 sigprv(i,2)=max(sigprv(i,2), zero)
551 sigprv(i,3)=max(sigprv(i,3), zero)
552 !! ENDIF
553 ENDIF
554C TRANSF ORM PRINCIPAL STRESSES (SIGMA_K) TO GLOBAL STRESSES (SIGMA_IJ)
555 sigsxx(i) = dirprv(i,1,1)*dirprv(i,1,1)*sigprv(i,1)
556 . + dirprv(i,1,2)*dirprv(i,1,2)*sigprv(i,2)
557 . + dirprv(i,1,3)*dirprv(i,1,3)*sigprv(i,3)
558!!
559 sigsyy(i) = dirprv(i,2,2)*dirprv(i,2,2)*sigprv(i,2)
560 . + dirprv(i,2,3)*dirprv(i,2,3)*sigprv(i,3)
561 . + dirprv(i,2,1)*dirprv(i,2,1)*sigprv(i,1)
562!!
563 sigszz(i) = dirprv(i,3,3)*dirprv(i,3,3)*sigprv(i,3)
564 . + dirprv(i,3,1)*dirprv(i,3,1)*sigprv(i,1)
565 . + dirprv(i,3,2)*dirprv(i,3,2)*sigprv(i,2)
566!!
567 sigsxy(i) = dirprv(i,1,1)*dirprv(i,2,1)*sigprv(i,1)
568 . + dirprv(i,1,2)*dirprv(i,2,2)*sigprv(i,2)
569 . + dirprv(i,1,3)*dirprv(i,2,3)*sigprv(i,3)
570!!
571 sigsyz(i) = dirprv(i,2,2)*dirprv(i,3,2)*sigprv(i,2)
572 . + dirprv(i,2,3)*dirprv(i,3,3)*sigprv(i,3)
573 . + dirprv(i,2,1)*dirprv(i,3,1)*sigprv(i,1)
574!!
575 sigszx(i) = dirprv(i,3,3)*dirprv(i,1,3)*sigprv(i,3)
576 . + dirprv(i,3,1)*dirprv(i,1,1)*sigprv(i,1)
577 . + dirprv(i,3,2)*dirprv(i,1,2)*sigprv(i,2)
578!!
579 ENDDO
580C fin chgt
581C COMPUTE UPDATED STRESSES AND SOUND SPEED
582C
583 DO i=1,nel
584 signxx(i)=off(i)*sigsxx(i)-sigair(i)
585 signyy(i)=off(i)*sigsyy(i)-sigair(i)
586 signzz(i)=off(i)*sigszz(i)-sigair(i)
587 signxy(i)=off(i)*sigsxy(i)
588 signyz(i)=off(i)*sigsyz(i)
589 signzx(i)=off(i)*sigszx(i)
590 soundsp(i)=sqrt(ey/rho0(i))
591 ENDDO
592 END SELECT
593
594 RETURN
595
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
dreh
subroutine dreh(b, tr, ii, jj, ken)
Definition sigeps33.F:739
jacobiew
subroutine jacobiew(a, n, ew, ev, nrot)
Definition sigeps33.F:607
valpvecdp_v
subroutine valpvecdp_v(sig, val, vec, nel)
Definition sigeps33.F:2231
valpvec_v
subroutine valpvec_v(sig, val, vec, nel)
Definition sigeps33.F:1902

◆ valpvec()

subroutine valpvec ( sig,
val,
vec,
integer nel )

Definition at line 793 of file sigeps33.F.

794C-----------------------------------------------
795C I m p l i c i t T y p e s
796C-----------------------------------------------
797#include "implicit_f.inc"
798C-----------------------------------------------
799C G l o b a l P a r a m e t e r s
800C-----------------------------------------------
801#include "mvsiz_p.inc"
802C-----------------------------------------------
803C D u m m y A r g u m e n t s
804C-----------------------------------------------
805 my_real
806 . sig(6,*), val(3,*), vec(9,*)
807 INTEGER NEL
808C-----------------------------------------------
809C L o c a l V a r i a b l e s
810C-----------------------------------------------
811 INTEGER I, L, N, NN, LMAX, LMAXV(MVSIZ),NINDEX3, INDEX3(MVSIZ),
812 . NINDEX1, NINDEX2, INDEX1(MVSIZ), INDEX2(MVSIZ)
813 my_real
814 . cs(6), str(3,mvsiz), a(3,3,mvsiz), v(3,3,mvsiz), b(3,3,mvsiz),
815 . xmag(3,mvsiz), pr, aa, bb, aaa(mvsiz),
816 . cc, angp, dd, ftpi, ttpi, strmax,
817 . tol1(mvsiz), tol2(mvsiz), xmaxv(mvsiz), norminf(mvsiz),
818 . vmag, s11,
819 . s21, s31, s12, s22, s32, s13, s23, s33, a11, a12, a13, a21,
820 . a22, a23, a31, a32, a33,
821 . mdemi, xmaxinv, flm
822 REAL FLMIN
823C-----------------------------------------------
824C DATA FTPI,TTPI / 4.188790205, 2.094395102 /
825C FTPI=(4/3)*PI, TTPI=(2/3)*PI
826C
827C DEVIATEUR PRINCIPAL DE CONTRAINTE
828C . . . . . . . . . . . . . . . . . . .
829 mdemi = -half
830 ttpi = acos(mdemi)
831 ftpi = two*ttpi
832C precision minimum dependant du type REAL ou DOUBLE
833 CALL floatmin(cs(1),cs(2),flmin)
834 flm = two*sqrt(flmin)
835 nindex3=0
836 DO nn = 1, nel
837 cs(1) = sig(1,nn)
838 cs(2) = sig(2,nn)
839 cs(3) = sig(3,nn)
840 cs(4) = sig(4,nn)
841 cs(5) = sig(5,nn)
842 cs(6) = sig(6,nn)
843 pr = -(cs(1)+cs(2)+cs(3))* third
844 cs(1) = cs(1) + pr
845 cs(2) = cs(2) + pr
846 cs(3) = cs(3) + pr
847 aaa(nn)=cs(4)**2 + cs(5)**2 + cs(6)**2 - cs(1)*cs(2)
848 & - cs(2)*cs(3) - cs(1)*cs(3)
849 norminf(nn) = max(abs(cs(1)),abs(cs(2)),abs(cs(3)),
850 & abs(cs(4)),abs(cs(5)),abs(cs(6)))
851 norminf(nn) = em10*norminf(nn)
852C cas racine triple
853c AA = MAX(AAA(NN),NORMINF(NN),EM20)
854 aa = max(aaa(nn),norminf(nn))
855C
856 bb=cs(1)*cs(5)**2 + cs(2)*cs(6)**2
857 & + cs(3)*cs(4)**2 - cs(1)*cs(2)*cs(3)
858 & - two*cs(4)*cs(5)*cs(6)
859C
860 cc=-sqrt(twenty7/max(em20,aa))*bb*half/max(em20,aa)
861 cc= min(cc,one)
862 cc= max(cc,-one)
863 angp=acos(cc) * third
864 dd=two*sqrt(aa * third)
865 str(1,nn)=dd*cos(angp)
866 str(2,nn)=dd*cos(angp+ftpi)
867 str(3,nn)=dd*cos(angp+ttpi)
868C
869 val(1,nn) = str(1,nn)-pr
870 val(2,nn) = str(2,nn)-pr
871 val(3,nn) = str(3,nn)-pr
872C renforcement de precision en compression----
873 IF(abs(str(3,nn))>abs(str(1,nn))
874 & .AND.aaa(nn)>norminf(nn)) THEN
875 aa=str(1,nn)
876 str(1,nn)=str(3,nn)
877 str(3,nn)=aa
878 nindex3 = nindex3+1
879 index3(nindex3) = nn
880 ENDIF
881C . . . . . . . . . . .
882C VECTEURS PROPRES
883C . . . . . . . . . . .
884 strmax= max(abs(str(1,nn)),abs(str(3,nn)))
885 tol1(nn)= max(em20,flm*strmax**2)
886 tol2(nn)=flm*strmax/3
887 a(1,1,nn)=cs(1)-str(1,nn)
888 a(2,2,nn)=cs(2)-str(1,nn)
889 a(3,3,nn)=cs(3)-str(1,nn)
890 a(1,2,nn)=cs(4)
891 a(2,1,nn)=cs(4)
892 a(2,3,nn)=cs(5)
893 a(3,2,nn)=cs(5)
894 a(1,3,nn)=cs(6)
895 a(3,1,nn)=cs(6)
896C
897 b(1,1,nn)=a(2,1,nn)*a(3,2,nn)-a(3,1,nn)
898 . *a(2,2,nn)
899 b(1,2,nn)=a(2,2,nn)*a(3,3,nn)-a(3,2,nn)
900 . *a(2,3,nn)
901 b(1,3,nn)=a(2,3,nn)*a(3,1,nn)-a(3,3,nn)
902 . *a(2,1,nn)
903 b(2,1,nn)=a(3,1,nn)*a(1,2,nn)-a(1,1,nn)
904 . *a(3,2,nn)
905 b(2,2,nn)=a(3,2,nn)*a(1,3,nn)-a(1,2,nn)
906 . *a(3,3,nn)
907 b(2,3,nn)=a(3,3,nn)*a(1,1,nn)-a(1,3,nn)
908 . *a(3,1,nn)
909 b(3,1,nn)=a(1,1,nn)*a(2,2,nn)-a(2,1,nn)
910 . *a(1,2,nn)
911 b(3,2,nn)=a(1,2,nn)*a(2,3,nn)-a(2,2,nn)
912 . *a(1,3,nn)
913 b(3,3,nn)=a(1,3,nn)*a(2,1,nn)-a(2,3,nn)
914 . *a(1,1,nn)
915 xmag(1,nn)=sqrt(b(1,1,nn)**2+b(2,1,nn)**2+b(3,1,nn)**2)
916 xmag(2,nn)=sqrt(b(1,2,nn)**2+b(2,2,nn)**2+b(3,2,nn)**2)
917 xmag(3,nn)=sqrt(b(1,3,nn)**2+b(2,3,nn)**2+b(3,3,nn)**2)
918
919 ENDDO
920C
921 nindex1 = 0
922 nindex2 = 0
923 DO nn = 1, nel
924 xmaxv(nn)=max(xmag(1,nn),xmag(2,nn),xmag(3,nn))
925 IF(xmag(1,nn)==xmaxv(nn)) THEN
926 lmaxv(nn) = 1
927 ELSEIF(xmag(2,nn)==xmaxv(nn)) THEN
928 lmaxv(nn) = 2
929 ELSE
930 lmaxv(nn) = 3
931 ENDIF
932 IF(aaa(nn)<norminf(nn)) THEN
933 val(1,nn) = sig(1,nn)
934 val(2,nn) = sig(2,nn)
935 val(3,nn) = sig(3,nn)
936 v(1,1,nn) = one
937 v(2,1,nn) = zero
938 v(3,1,nn) = zero
939 v(1,2,nn) = zero
940 v(2,2,nn) = one
941 v(3,2,nn) = zero
942
943 ELSEIF(xmaxv(nn)>tol1(nn)) THEN
944 nindex1 = nindex1 + 1
945 index1(nindex1) = nn
946 ELSE
947 nindex2 = nindex2 + 1
948 index2(nindex2) = nn
949 ENDIF
950 ENDDO
951C
952#include "vectorize.inc"
953 DO n = 1, nindex1
954 nn = index1(n)
955 lmax = lmaxv(nn)
956 xmaxinv = one/xmaxv(nn)
957 v(1,1,nn)=b(1,lmax,nn)*xmaxinv
958 v(2,1,nn)=b(2,lmax,nn)*xmaxinv
959 v(3,1,nn)=b(3,lmax,nn)*xmaxinv
960 a(1,1,nn)=a(1,1,nn)+str(1,nn)-str(3,nn)
961 a(2,2,nn)=a(2,2,nn)+str(1,nn)-str(3,nn)
962 a(3,3,nn)=a(3,3,nn)+str(1,nn)-str(3,nn)
963C
964 b(1,1,nn)=a(2,1,nn)*v(3,1,nn)-a(3,1,nn)*v(2,1,nn)
965 b(1,2,nn)=a(2,2,nn)*v(3,1,nn)-a(3,2,nn)*v(2,1,nn)
966 b(1,3,nn)=a(2,3,nn)*v(3,1,nn)-a(3,3,nn)*v(2,1,nn)
967 b(2,1,nn)=a(3,1,nn)*v(1,1,nn)-a(1,1,nn)*v(3,1,nn)
968 b(2,2,nn)=a(3,2,nn)*v(1,1,nn)-a(1,2,nn)*v(3,1,nn)
969 b(2,3,nn)=a(3,3,nn)*v(1,1,nn)-a(1,3,nn)*v(3,1,nn)
970 b(3,1,nn)=a(1,1,nn)*v(2,1,nn)-a(2,1,nn)*v(1,1,nn)
971 b(3,2,nn)=a(1,2,nn)*v(2,1,nn)-a(2,2,nn)*v(1,1,nn)
972 b(3,3,nn)=a(1,3,nn)*v(2,1,nn)-a(2,3,nn)*v(1,1,nn)
973 xmag(1,nn)=sqrt(b(1,1,nn)**2+b(2,1,nn)**2+b(3,1,nn)**2)
974 xmag(2,nn)=sqrt(b(1,2,nn)**2+b(2,2,nn)**2+b(3,2,nn)**2)
975 xmag(3,nn)=sqrt(b(1,3,nn)**2+b(2,3,nn)**2+b(3,3,nn)**2)
976C
977 xmaxv(nn)=max(xmag(1,nn),xmag(2,nn),xmag(3,nn))
978 ENDDO
979C
980#include "vectorize.inc"
981 DO n = 1, nindex1
982 nn = index1(n)
983 IF(xmag(1,nn)==xmaxv(nn)) THEN
984 lmaxv(nn) = 1
985 ELSEIF(xmag(2,nn)==xmaxv(nn)) THEN
986 lmaxv(nn) = 2
987 ELSE
988 lmaxv(nn) = 3
989 ENDIF
990C
991 vmag=sqrt(v(1,1,nn)**2+v(2,1,nn)**2)
992 lmax = lmaxv(nn)
993 IF(xmaxv(nn)>tol2(nn))THEN
994 xmaxinv = one/xmaxv(nn)
995 v(1,3,nn)=b(1,lmax,nn)*xmaxinv
996 v(2,3,nn)=b(2,lmax,nn)*xmaxinv
997 v(3,3,nn)=b(3,lmax,nn)*xmaxinv
998 v(1,2,nn)=v(2,3,nn)*v(3,1,nn)-v(2,1,nn)*v(3,3,nn)
999 v(2,2,nn)=v(3,3,nn)*v(1,1,nn)-v(3,1,nn)*v(1,3,nn)
1000 v(3,2,nn)=v(1,3,nn)*v(2,1,nn)-v(1,1,nn)*v(2,3,nn)
1001 vmag=one/sqrt(v(1,2,nn)**2+v(2,2,nn)**2+v(3,2,nn)**2)
1002 v(1,2,nn)=v(1,2,nn)*vmag
1003 v(2,2,nn)=v(2,2,nn)*vmag
1004 v(3,2,nn)=v(3,2,nn)*vmag
1005 ELSEIF(vmag>tol2(nn))THEN
1006 v(1,2,nn)=-v(2,1,nn)/vmag
1007 v(2,2,nn)=v(1,1,nn)/vmag
1008 v(3,2,nn)=zero
1009 ELSE
1010 v(1,2,nn)=one
1011 v(2,2,nn)=zero
1012 v(3,2,nn)=zero
1013 ENDIF
1014 ENDDO
1015C . . . . . . . . . . . . .
1016C SOLUTION DOUBLE
1017C . . . . . . . . . . . . .
1018 DO n = 1, nindex2
1019 nn = index2(n)
1020 xmag(1,nn)=sqrt(a(1,1,nn)**2+a(2,1,nn)**2)
1021 xmag(2,nn)=sqrt(a(1,2,nn)**2+a(2,2,nn)**2)
1022 xmag(3,nn)=sqrt(a(1,3,nn)**2+a(2,3,nn)**2)
1023C
1024 xmaxv(nn)=max(xmag(1,nn),xmag(2,nn),xmag(3,nn))
1025 ENDDO
1026C
1027#include "vectorize.inc"
1028 DO n = 1, nindex2
1029 nn = index2(n)
1030 IF(xmag(1,nn)==xmaxv(nn)) THEN
1031 lmaxv(nn) = 1
1032 ELSEIF(xmag(2,nn)==xmaxv(nn)) THEN
1033 lmaxv(nn) = 2
1034 ELSE
1035 lmaxv(nn) = 3
1036 ENDIF
1037C
1038 lmax = lmaxv(nn)
1039 IF(max(abs(a(3,1,nn)),abs(a(3,2,nn)),abs(a(3,3,nn)))
1040 & <tol2(nn))THEN
1041 xmaxinv = one/xmaxv(nn)
1042 v(1,1,nn)= zero
1043 v(2,1,nn)= zero
1044 v(3,1,nn)= one
1045 v(1,2,nn)=-a(2,lmax,nn)*xmaxinv
1046 v(2,2,nn)= a(1,lmax,nn)*xmaxinv
1047 v(3,2,nn)= zero
1048C
1049 ELSEIF(xmaxv(nn)>tol2(nn))THEN
1050 xmaxinv = one/xmaxv(nn)
1051 v(1,1,nn)=-a(2,lmax,nn)*xmaxinv
1052 v(2,1,nn)= a(1,lmax,nn)*xmaxinv
1053 v(3,1,nn)= zero
1054 v(1,2,nn)=-a(3,lmax,nn)*v(2,1,nn)
1055 v(2,2,nn)= a(3,lmax,nn)*v(1,1,nn)
1056 v(3,2,nn)= a(1,lmax,nn)*v(2,1,nn)-a(2,lmax,nn)*v(1,1,nn)
1057 vmag=one/sqrt(v(1,2,nn)**2+v(2,2,nn)**2+v(3,2,nn)**2)
1058 v(1,2,nn)=v(1,2,nn)*vmag
1059 v(2,2,nn)=v(2,2,nn)*vmag
1060 v(3,2,nn)=v(3,2,nn)*vmag
1061 ELSE
1062 v(1,1,nn) = one
1063 v(2,1,nn) = zero
1064 v(3,1,nn) = zero
1065 v(1,2,nn) = zero
1066 v(2,2,nn) = one
1067 v(3,2,nn) = zero
1068 ENDIF
1069 ENDDO
1070C
1071 DO nn = 1, nel
1072 vec(1,nn)=v(1,1,nn)
1073 vec(2,nn)=v(2,1,nn)
1074 vec(3,nn)=v(3,1,nn)
1075 vec(4,nn)=v(1,2,nn)
1076 vec(5,nn)=v(2,2,nn)
1077 vec(6,nn)=v(3,2,nn)
1078 vec(7,nn)=vec(2,nn)*vec(6,nn)-vec(3,nn)*vec(5,nn)
1079 vec(8,nn)=vec(3,nn)*vec(4,nn)-vec(1,nn)*vec(6,nn)
1080 vec(9,nn)=vec(1,nn)*vec(5,nn)-vec(2,nn)*vec(4,nn)
1081 ENDDO
1082C reecriture pour contourner probleme sur itanium2 comp 9. + latency=16
1083 DO n = 1, nindex3
1084 nn = index3(n)
1085C str utilise com tableau temporaire au lieu de scalaires temporaires
1086 str(1,nn)=vec(7,nn)
1087 str(2,nn)=vec(8,nn)
1088 str(3,nn)=vec(9,nn)
1089 ENDDO
1090 DO n = 1, nindex3
1091 nn = index3(n)
1092 vec(7,nn)=vec(1,nn)
1093 vec(8,nn)=vec(2,nn)
1094 vec(9,nn)=vec(3,nn)
1095 vec(1,nn)=-str(1,nn)
1096 vec(2,nn)=-str(2,nn)
1097 vec(3,nn)=-str(3,nn)
1098 ENDDO
1099C
1100 RETURN
floatmin
void floatmin(int *a, int *b, float *flm)
Definition precision.c:71

◆ valpvec_v()

subroutine valpvec_v ( sig,
val,
vec,
integer nel )

Definition at line 1901 of file sigeps33.F.

1902C-----------------------------------------------
1903C I m p l i c i t T y p e s
1904C-----------------------------------------------
1905#include "implicit_f.inc"
1906C-----------------------------------------------
1907C G l o b a l P a r a m e t e r s
1908C-----------------------------------------------
1909#include "mvsiz_p.inc"
1910C-----------------------------------------------
1911C D u m m y A r g u m e n t s
1912C-----------------------------------------------
1913 my_real
1914 . sig(mvsiz,6), val(mvsiz,3), vec(mvsiz,9)
1915 INTEGER NEL
1916C-----------------------------------------------
1917C L o c a l V a r i a b l e s
1918C-----------------------------------------------
1919 INTEGER I, L, N, NN, LMAX, LMAXV(MVSIZ),NINDEX3, INDEX3(MVSIZ),
1920 . NINDEX1, NINDEX2, INDEX1(MVSIZ), INDEX2(MVSIZ)
1921 my_real
1922 . cs(6), str(mvsiz,3), a(mvsiz,3,3), v(mvsiz,3,3), b(mvsiz,3,3),
1923 . xmag(mvsiz,3), pr, aa, bb, aaa(mvsiz),
1924 . cc, angp, dd, ftpi, ttpi, strmax,
1925 . tol1(mvsiz), tol2(mvsiz), xmaxv(mvsiz), norminf(mvsiz),
1926 . vmag, s11,
1927 . s21, s31, s12, s22, s32, s13, s23, s33, a11, a12, a13, a21,
1928 . a22, a23, a31, a32, a33,
1929 . mdemi, xmaxinv, flm
1930 REAL FLMIN
1931C-----------------------------------------------
1932C DATA FTPI,TTPI / 4.188790205, 2.094395102 /
1933C FTPI=(4/3)*PI, TTPI=(2/3)*PI
1934C
1935C DEVIATEUR PRINCIPAL DE CONTRAINTE
1936C . . . . . . . . . . . . . . . . . . .
1937 mdemi = -half
1938 ttpi = acos(mdemi)
1939 ftpi = two*ttpi
1940C precision minimum dependant du type REAL ou DOUBLE
1941 CALL floatmin(cs(1),cs(2),flmin)
1942 flm = two*sqrt(flmin)
1943 nindex3=0
1944 DO nn = 1, nel
1945 cs(1) = sig(nn,1)
1946 cs(2) = sig(nn,2)
1947 cs(3) = sig(nn,3)
1948 cs(4) = sig(nn,4)
1949 cs(5) = sig(nn,5)
1950 cs(6) = sig(nn,6)
1951 pr = -(cs(1)+cs(2)+cs(3))* third
1952 cs(1) = cs(1) + pr
1953 cs(2) = cs(2) + pr
1954 cs(3) = cs(3) + pr
1955 aaa(nn)=cs(4)**2 + cs(5)**2 + cs(6)**2 - cs(1)*cs(2)
1956 & - cs(2)*cs(3) - cs(1)*cs(3)
1957 norminf(nn) = max(abs(cs(1)),abs(cs(2)),abs(cs(3)),
1958 & abs(cs(4)),abs(cs(5)),abs(cs(6)))
1959 norminf(nn) = em10*norminf(nn)
1960C cas racine triple
1961c AA = MAX(AAA(NN),NORMINF(NN),EM20)
1962 aa = max(aaa(nn),norminf(nn))
1963C
1964 bb=cs(1)*cs(5)**2 + cs(2)*cs(6)**2
1965 & + cs(3)*cs(4)**2 - cs(1)*cs(2)*cs(3)
1966 & - two*cs(4)*cs(5)*cs(6)
1967C
1968 cc=-sqrt(twenty7/max(em20,aa))*bb*half/max(em20,aa)
1969 cc= min(cc,one)
1970 cc= max(cc,-one)
1971 angp=acos(cc) * third
1972 dd=two*sqrt(aa * third)
1973 str(nn,1)=dd*cos(angp)
1974 str(nn,2)=dd*cos(angp+ftpi)
1975 str(nn,3)=dd*cos(angp+ttpi)
1976C
1977 val(nn,1) = str(nn,1)-pr
1978 val(nn,2) = str(nn,2)-pr
1979 val(nn,3) = str(nn,3)-pr
1980C renforcement de precision en compression----
1981 IF(abs(str(nn,3))>abs(str(nn,1))
1982 & .AND.aaa(nn)>norminf(nn)) THEN
1983 aa=str(nn,1)
1984 str(nn,1)=str(nn,3)
1985 str(nn,3)=aa
1986 nindex3 = nindex3+1
1987 index3(nindex3) = nn
1988 ENDIF
1989C . . . . . . . . . . .
1990C VECTEURS PROPRES
1991C . . . . . . . . . . .
1992 strmax= max(abs(str(nn,1)),abs(str(nn,3)))
1993 tol1(nn)= max(em20,flm*strmax**2)
1994 tol2(nn)=flm*strmax/3
1995 a(nn,1,1)=cs(1)-str(nn,1)
1996 a(nn,2,2)=cs(2)-str(nn,1)
1997 a(nn,3,3)=cs(3)-str(nn,1)
1998 a(nn,1,2)=cs(4)
1999 a(nn,2,1)=cs(4)
2000 a(nn,2,3)=cs(5)
2001 a(nn,3,2)=cs(5)
2002 a(nn,1,3)=cs(6)
2003 a(nn,3,1)=cs(6)
2004C
2005 b(nn,1,1)=a(nn,2,1)*a(nn,3,2)-a(nn,3,1)
2006 . *a(nn,2,2)
2007 b(nn,1,2)=a(nn,2,2)*a(nn,3,3)-a(nn,3,2)
2008 . *a(nn,2,3)
2009 b(nn,1,3)=a(nn,2,3)*a(nn,3,1)-a(nn,3,3)
2010 . *a(nn,2,1)
2011 b(nn,2,1)=a(nn,3,1)*a(nn,1,2)-a(nn,1,1)
2012 . *a(nn,3,2)
2013 b(nn,2,2)=a(nn,3,2)*a(nn,1,3)-a(nn,1,2)
2014 . *a(nn,3,3)
2015 b(nn,2,3)=a(nn,3,3)*a(nn,1,1)-a(nn,1,3)
2016 . *a(nn,3,1)
2017 b(nn,3,1)=a(nn,1,1)*a(nn,2,2)-a(nn,2,1)
2018 . *a(nn,1,2)
2019 b(nn,3,2)=a(nn,1,2)*a(nn,2,3)-a(nn,2,2)
2020 . *a(nn,1,3)
2021 b(nn,3,3)=a(nn,1,3)*a(nn,2,1)-a(nn,2,3)
2022 . *a(nn,1,1)
2023 xmag(nn,1)=sqrt(b(nn,1,1)**2+b(nn,2,1)**2+b(nn,3,1)**2)
2024 xmag(nn,2)=sqrt(b(nn,1,2)**2+b(nn,2,2)**2+b(nn,3,2)**2)
2025 xmag(nn,3)=sqrt(b(nn,1,3)**2+b(nn,2,3)**2+b(nn,3,3)**2)
2026
2027 ENDDO
2028C
2029 nindex1 = 0
2030 nindex2 = 0
2031 DO nn = 1, nel
2032 xmaxv(nn)=max(xmag(nn,1),xmag(nn,2),xmag(nn,3))
2033 IF(xmag(nn,1)==xmaxv(nn)) THEN
2034 lmaxv(nn) = 1
2035 ELSEIF(xmag(nn,2)==xmaxv(nn)) THEN
2036 lmaxv(nn) = 2
2037 ELSE
2038 lmaxv(nn) = 3
2039 ENDIF
2040 IF(aaa(nn)<norminf(nn)) THEN
2041 val(nn,1) = sig(nn,1)
2042 val(nn,2) = sig(nn,2)
2043 val(nn,3) = sig(nn,3)
2044 v(nn,1,1) = one
2045 v(nn,2,1) = zero
2046 v(nn,3,1) = zero
2047 v(nn,1,2) = zero
2048 v(nn,2,2) = one
2049 v(nn,3,2) = zero
2050
2051 ELSEIF(xmaxv(nn)>tol1(nn)) THEN
2052 nindex1 = nindex1 + 1
2053 index1(nindex1) = nn
2054 ELSE
2055 nindex2 = nindex2 + 1
2056 index2(nindex2) = nn
2057 ENDIF
2058 ENDDO
2059C
2060#include "vectorize.inc"
2061 DO n = 1, nindex1
2062 nn = index1(n)
2063 lmax = lmaxv(nn)
2064 xmaxinv = one/xmaxv(nn)
2065 v(nn,1,1)=b(nn,1,lmax)*xmaxinv
2066 v(nn,2,1)=b(nn,2,lmax)*xmaxinv
2067 v(nn,3,1)=b(nn,3,lmax)*xmaxinv
2068 a(nn,1,1)=a(nn,1,1)+str(nn,1)-str(nn,3)
2069 a(nn,2,2)=a(nn,2,2)+str(nn,1)-str(nn,3)
2070 a(nn,3,3)=a(nn,3,3)+str(nn,1)-str(nn,3)
2071C
2072 b(nn,1,1)=a(nn,2,1)*v(nn,3,1)-a(nn,3,1)*v(nn,2,1)
2073 b(nn,1,2)=a(nn,2,2)*v(nn,3,1)-a(nn,3,2)*v(nn,2,1)
2074 b(nn,1,3)=a(nn,2,3)*v(nn,3,1)-a(nn,3,3)*v(nn,2,1)
2075 b(nn,2,1)=a(nn,3,1)*v(nn,1,1)-a(nn,1,1)*v(nn,3,1)
2076 b(nn,2,2)=a(nn,3,2)*v(nn,1,1)-a(nn,1,2)*v(nn,3,1)
2077 b(nn,2,3)=a(nn,3,3)*v(nn,1,1)-a(nn,1,3)*v(nn,3,1)
2078 b(nn,3,1)=a(nn,1,1)*v(nn,2,1)-a(nn,2,1)*v(nn,1,1)
2079 b(nn,3,2)=a(nn,1,2)*v(nn,2,1)-a(nn,2,2)*v(nn,1,1)
2080 b(nn,3,3)=a(nn,1,3)*v(nn,2,1)-a(nn,2,3)*v(nn,1,1)
2081 xmag(nn,1)=sqrt(b(nn,1,1)**2+b(nn,2,1)**2+b(nn,3,1)**2)
2082 xmag(nn,2)=sqrt(b(nn,1,2)**2+b(nn,2,2)**2+b(nn,3,2)**2)
2083 xmag(nn,3)=sqrt(b(nn,1,3)**2+b(nn,2,3)**2+b(nn,3,3)**2)
2084C
2085 xmaxv(nn)=max(xmag(nn,1),xmag(nn,2),xmag(nn,3))
2086 ENDDO
2087C
2088#include "vectorize.inc"
2089 DO n = 1, nindex1
2090 nn = index1(n)
2091 IF(xmag(nn,1)==xmaxv(nn)) THEN
2092 lmaxv(nn) = 1
2093 ELSEIF(xmag(nn,2)==xmaxv(nn)) THEN
2094 lmaxv(nn) = 2
2095 ELSE
2096 lmaxv(nn) = 3
2097 ENDIF
2098C
2099 vmag=sqrt(v(nn,1,1)**2+v(nn,2,1)**2)
2100 lmax = lmaxv(nn)
2101 IF(xmaxv(nn)>tol2(nn))THEN
2102 xmaxinv = one/xmaxv(nn)
2103 v(nn,1,3)=b(nn,1,lmax)*xmaxinv
2104 v(nn,2,3)=b(nn,2,lmax)*xmaxinv
2105 v(nn,3,3)=b(nn,3,lmax)*xmaxinv
2106 v(nn,1,2)=v(nn,2,3)*v(nn,3,1)-v(nn,2,1)*v(nn,3,3)
2107 v(nn,2,2)=v(nn,3,3)*v(nn,1,1)-v(nn,3,1)*v(nn,1,3)
2108 v(nn,3,2)=v(nn,1,3)*v(nn,2,1)-v(nn,1,1)*v(nn,2,3)
2109 vmag=one/sqrt(v(nn,1,2)**2+v(nn,2,2)**2+v(nn,3,2)**2)
2110 v(nn,1,2)=v(nn,1,2)*vmag
2111 v(nn,2,2)=v(nn,2,2)*vmag
2112 v(nn,3,2)=v(nn,3,2)*vmag
2113 ELSEIF(vmag>tol2(nn))THEN
2114 v(nn,1,2)=-v(nn,2,1)/vmag
2115 v(nn,2,2)=v(nn,1,1)/vmag
2116 v(nn,3,2)=zero
2117 ELSE
2118 v(nn,1,2)=one
2119 v(nn,2,2)=zero
2120 v(nn,3,2)=zero
2121 ENDIF
2122 ENDDO
2123C . . . . . . . . . . . . .
2124C SOLUTION DOUBLE
2125C . . . . . . . . . . . . .
2126 DO n = 1, nindex2
2127 nn = index2(n)
2128 xmag(nn,1)=sqrt(a(nn,1,1)**2+a(nn,2,1)**2)
2129 xmag(nn,2)=sqrt(a(nn,1,2)**2+a(nn,2,2)**2)
2130 xmag(nn,3)=sqrt(a(nn,1,3)**2+a(nn,2,3)**2)
2131C
2132 xmaxv(nn)=max(xmag(nn,1),xmag(nn,2),xmag(nn,3))
2133 ENDDO
2134C
2135#include "vectorize.inc"
2136 DO n = 1, nindex2
2137 nn = index2(n)
2138 IF(xmag(nn,1)==xmaxv(nn)) THEN
2139 lmaxv(nn) = 1
2140 ELSEIF(xmag(nn,2)==xmaxv(nn)) THEN
2141 lmaxv(nn) = 2
2142 ELSE
2143 lmaxv(nn) = 3
2144 ENDIF
2145C
2146 lmax = lmaxv(nn)
2147 IF(max(abs(a(nn,3,1)),abs(a(nn,3,2)),abs(a(nn,3,3)))
2148 & <tol2(nn))THEN
2149 xmaxinv = one/xmaxv(nn)
2150 v(nn,1,1)= zero
2151 v(nn,2,1)= zero
2152 v(nn,3,1)= one
2153 v(nn,1,2)=-a(nn,2,lmax)*xmaxinv
2154 v(nn,2,2)= a(nn,1,lmax)*xmaxinv
2155 v(nn,3,2)= zero
2156C
2157 ELSEIF(xmaxv(nn)>tol2(nn))THEN
2158 xmaxinv = one/xmaxv(nn)
2159 v(nn,1,1)=-a(nn,2,lmax)*xmaxinv
2160 v(nn,2,1)= a(nn,1,lmax)*xmaxinv
2161 v(nn,3,1)= zero
2162 v(nn,1,2)=-a(nn,3,lmax)*v(nn,2,1)
2163 v(nn,2,2)= a(nn,3,lmax)*v(nn,1,1)
2164 v(nn,3,2)= a(nn,1,lmax)*v(nn,2,1)-a(nn,2,lmax)*v(nn,1,1)
2165 vmag=one/sqrt(v(nn,1,2)**2+v(nn,2,2)**2+v(nn,3,2)**2)
2166 v(nn,1,2)=v(nn,1,2)*vmag
2167 v(nn,2,2)=v(nn,2,2)*vmag
2168 v(nn,3,2)=v(nn,3,2)*vmag
2169 ELSE
2170 v(nn,1,1) = one
2171 v(nn,2,1) = zero
2172 v(nn,3,1) = zero
2173 v(nn,1,2) = zero
2174 v(nn,2,2) = one
2175 v(nn,3,2) = zero
2176 ENDIF
2177 ENDDO
2178C
2179 DO nn = 1,nel
2180 vec(nn,1)=v(nn,1,1)
2181 vec(nn,2)=v(nn,2,1)
2182 vec(nn,3)=v(nn,3,1)
2183 vec(nn,4)=v(nn,1,2)
2184 vec(nn,5)=v(nn,2,2)
2185 vec(nn,6)=v(nn,3,2)
2186 vec(nn,7)=vec(nn,2)*vec(nn,6)-vec(nn,3)*vec(nn,5)
2187 vec(nn,8)=vec(nn,3)*vec(nn,4)-vec(nn,1)*vec(nn,6)
2188 vec(nn,9)=vec(nn,1)*vec(nn,5)-vec(nn,2)*vec(nn,4)
2189 ENDDO
2190C reecriture pour contourner probleme sur itanium2 comp 9. + latency=16
2191 DO n = 1, nindex3
2192 nn = index3(n)
2193C str utilise com tableau temporaire au lieu de scalaires temporaires
2194 str(nn,1)=vec(nn,7)
2195 str(nn,2)=vec(nn,8)
2196 str(nn,3)=vec(nn,9)
2197 ENDDO
2198 DO n = 1, nindex3
2199 nn = index3(n)
2200 vec(nn,7)=vec(nn,1)
2201 vec(nn,8)=vec(nn,2)
2202 vec(nn,9)=vec(nn,3)
2203 vec(nn,1)=-str(nn,1)
2204 vec(nn,2)=-str(nn,2)
2205 vec(nn,3)=-str(nn,3)
2206 ENDDO
2207C
2208 RETURN

◆ valpvecdp()

subroutine valpvecdp ( sig,
val,
vec,
integer nel )

Definition at line 1111 of file sigeps33.F.

1112C-----------------------------------------------
1113C I m p l i c i t T y p e s
1114C-----------------------------------------------
1115#include "implicit_f.inc"
1116C-----------------------------------------------
1117C G l o b a l P a r a m e t e r s
1118C-----------------------------------------------
1119#include "mvsiz_p.inc"
1120C-----------------------------------------------
1121C D u m m y A r g u m e n t s
1122C-----------------------------------------------
1123 my_real
1124 . sig(6,*), val(3,*), vec(9,*)
1125 INTEGER NEL
1126C-----------------------------------------------
1127C L o c a l V a r i a b l e s
1128C-----------------------------------------------
1129 INTEGER I, L, N, NN, LMAX, LMAXV(MVSIZ),NINDEX3, INDEX3(MVSIZ),
1130 . NINDEX1, NINDEX2, INDEX1(MVSIZ), INDEX2(MVSIZ)
1131 double precision
1132 . cs(6), str(3,mvsiz), a(3,3,mvsiz), v(3,3,mvsiz), b(3,3,mvsiz),
1133 . xmag(3,mvsiz), pr, aa, bb, aaa(mvsiz),
1134 . cc, angp, dd, ftpi, ttpi, strmax,
1135 . tol1(mvsiz), tol2(mvsiz), xmaxv(mvsiz), norminf(mvsiz),
1136 . vmag, s11,
1137 . s21, s31, s12, s22, s32, s13, s23, s33, a11, a12, a13, a21,
1138 . a22, a23, a31, a32, a33,
1139 . mdemi, xmaxinv, flm,
1140 . valdp(3,mvsiz),vecdp(9,mvsiz)
1141 REAL FLMIN
1142C-----------------------------------------------
1143C DATA FTPI,TTPI / 4.188790205, 2.094395102 /
1144C FTPI=(4/3)*PI, TTPI=(2/3)*PI
1145C
1146C DEVIATEUR PRINCIPAL DE CONTRAINTE
1147C . . . . . . . . . . . . . . . . . . .
1148 mdemi = -half
1149 ttpi = acos(mdemi)
1150 ftpi = two*ttpi
1151C precision minimum dependant du type REAL ou DOUBLE
1152 CALL floatmin(cs(1),cs(2),flmin)
1153 flm = two*sqrt(flmin)
1154 nindex3=0
1155 DO nn = 1, nel
1156 cs(1) = sig(1,nn)
1157 cs(2) = sig(2,nn)
1158 cs(3) = sig(3,nn)
1159 cs(4) = sig(4,nn)
1160 cs(5) = sig(5,nn)
1161 cs(6) = sig(6,nn)
1162 pr = -(cs(1)+cs(2)+cs(3)) * third
1163 cs(1) = cs(1) + pr
1164 cs(2) = cs(2) + pr
1165 cs(3) = cs(3) + pr
1166 aaa(nn)=cs(4)**2 + cs(5)**2 + cs(6)**2 - cs(1)*cs(2)
1167 & - cs(2)*cs(3) - cs(1)*cs(3)
1168 norminf(nn) = max(abs(cs(1)),abs(cs(2)),abs(cs(3)),
1169 & abs(cs(4)),abs(cs(5)),abs(cs(6)))
1170 norminf(nn) = em10*norminf(nn)
1171C cas racine triple
1172c AA = MAX(AAA(NN),NORMINF(NN),EM20)
1173 aa = max(aaa(nn),norminf(nn))
1174C
1175 bb=cs(1)*cs(5)**2 + cs(2)*cs(6)**2
1176 & + cs(3)*cs(4)**2 - cs(1)*cs(2)*cs(3)
1177 & - two*cs(4)*cs(5)*cs(6)
1178C
1179 cc=-sqrt(twenty7/max(em20,aa))*bb*half/max(em20,aa)
1180 cc= min(cc,one)
1181 cc= max(cc,-one)
1182 angp=acos(cc) * third
1183 dd=two*sqrt(aa * third)
1184 str(1,nn)=dd*cos(angp)
1185 str(2,nn)=dd*cos(angp+ftpi)
1186 str(3,nn)=dd*cos(angp+ttpi)
1187C
1188 valdp(1,nn) = str(1,nn)-pr
1189 valdp(2,nn) = str(2,nn)-pr
1190 valdp(3,nn) = str(3,nn)-pr
1191C renforcement de precision en compression simple----
1192 IF(abs(str(3,nn))>abs(str(1,nn))
1193 & .AND.aaa(nn)>norminf(nn)) THEN
1194 aa=str(1,nn)
1195 str(1,nn)=str(3,nn)
1196 str(3,nn)=aa
1197 nindex3 = nindex3+1
1198 index3(nindex3) = nn
1199 ENDIF
1200C . . . . . . . . . . .
1201C VECTEURS PROPRES
1202C . . . . . . . . . . .
1203 strmax= max(abs(str(1,nn)),abs(str(3,nn)))
1204 tol1(nn)= max(em20,flm*strmax**2)
1205 tol2(nn)=flm*strmax * third
1206 a(1,1,nn)=cs(1)-str(1,nn)
1207 a(2,2,nn)=cs(2)-str(1,nn)
1208 a(3,3,nn)=cs(3)-str(1,nn)
1209 a(1,2,nn)=cs(4)
1210 a(2,1,nn)=cs(4)
1211 a(2,3,nn)=cs(5)
1212 a(3,2,nn)=cs(5)
1213 a(1,3,nn)=cs(6)
1214 a(3,1,nn)=cs(6)
1215C
1216 b(1,1,nn)=a(2,1,nn)*a(3,2,nn)-a(3,1,nn)
1217 . *a(2,2,nn)
1218 b(1,2,nn)=a(2,2,nn)*a(3,3,nn)-a(3,2,nn)
1219 . *a(2,3,nn)
1220 b(1,3,nn)=a(2,3,nn)*a(3,1,nn)-a(3,3,nn)
1221 . *a(2,1,nn)
1222 b(2,1,nn)=a(3,1,nn)*a(1,2,nn)-a(1,1,nn)
1223 . *a(3,2,nn)
1224 b(2,2,nn)=a(3,2,nn)*a(1,3,nn)-a(1,2,nn)
1225 . *a(3,3,nn)
1226 b(2,3,nn)=a(3,3,nn)*a(1,1,nn)-a(1,3,nn)
1227 . *a(3,1,nn)
1228 b(3,1,nn)=a(1,1,nn)*a(2,2,nn)-a(2,1,nn)
1229 . *a(1,2,nn)
1230 b(3,2,nn)=a(1,2,nn)*a(2,3,nn)-a(2,2,nn)
1231 . *a(1,3,nn)
1232 b(3,3,nn)=a(1,3,nn)*a(2,1,nn)-a(2,3,nn)
1233 . *a(1,1,nn)
1234 xmag(1,nn)=sqrt(b(1,1,nn)**2+b(2,1,nn)**2+b(3,1,nn)**2)
1235 xmag(2,nn)=sqrt(b(1,2,nn)**2+b(2,2,nn)**2+b(3,2,nn)**2)
1236 xmag(3,nn)=sqrt(b(1,3,nn)**2+b(2,3,nn)**2+b(3,3,nn)**2)
1237
1238 ENDDO
1239C
1240 nindex1 = 0
1241 nindex2 = 0
1242 DO nn = 1, nel
1243 xmaxv(nn)=max(xmag(1,nn),xmag(2,nn),xmag(3,nn))
1244 IF(xmag(1,nn)==xmaxv(nn)) THEN
1245 lmaxv(nn) = 1
1246 ELSEIF(xmag(2,nn)==xmaxv(nn)) THEN
1247 lmaxv(nn) = 2
1248 ELSE
1249 lmaxv(nn) = 3
1250 ENDIF
1251 IF(aaa(nn)<norminf(nn)) THEN
1252 valdp(1,nn) = sig(1,nn)
1253 valdp(2,nn) = sig(2,nn)
1254 valdp(3,nn) = sig(3,nn)
1255 v(1,1,nn) = one
1256 v(2,1,nn) = zero
1257 v(3,1,nn) = zero
1258 v(1,2,nn) = zero
1259 v(2,2,nn) = one
1260 v(3,2,nn) = zero
1261 ELSEIF(xmaxv(nn)>tol1(nn)) THEN
1262 nindex1 = nindex1 + 1
1263 index1(nindex1) = nn
1264 ELSE
1265 nindex2 = nindex2 + 1
1266 index2(nindex2) = nn
1267 ENDIF
1268 ENDDO
1269C
1270#include "vectorize.inc"
1271 DO n = 1, nindex1
1272 nn = index1(n)
1273 lmax = lmaxv(nn)
1274 xmaxinv = one/xmaxv(nn)
1275 v(1,1,nn)=b(1,lmax,nn)*xmaxinv
1276 v(2,1,nn)=b(2,lmax,nn)*xmaxinv
1277 v(3,1,nn)=b(3,lmax,nn)*xmaxinv
1278 a(1,1,nn)=a(1,1,nn)+str(1,nn)-str(3,nn)
1279 a(2,2,nn)=a(2,2,nn)+str(1,nn)-str(3,nn)
1280 a(3,3,nn)=a(3,3,nn)+str(1,nn)-str(3,nn)
1281C
1282 b(1,1,nn)=a(2,1,nn)*v(3,1,nn)-a(3,1,nn)*v(2,1,nn)
1283 b(1,2,nn)=a(2,2,nn)*v(3,1,nn)-a(3,2,nn)*v(2,1,nn)
1284 b(1,3,nn)=a(2,3,nn)*v(3,1,nn)-a(3,3,nn)*v(2,1,nn)
1285 b(2,1,nn)=a(3,1,nn)*v(1,1,nn)-a(1,1,nn)*v(3,1,nn)
1286 b(2,2,nn)=a(3,2,nn)*v(1,1,nn)-a(1,2,nn)*v(3,1,nn)
1287 b(2,3,nn)=a(3,3,nn)*v(1,1,nn)-a(1,3,nn)*v(3,1,nn)
1288 b(3,1,nn)=a(1,1,nn)*v(2,1,nn)-a(2,1,nn)*v(1,1,nn)
1289 b(3,2,nn)=a(1,2,nn)*v(2,1,nn)-a(2,2,nn)*v(1,1,nn)
1290 b(3,3,nn)=a(1,3,nn)*v(2,1,nn)-a(2,3,nn)*v(1,1,nn)
1291 xmag(1,nn)=sqrt(b(1,1,nn)**2+b(2,1,nn)**2+b(3,1,nn)**2)
1292 xmag(2,nn)=sqrt(b(1,2,nn)**2+b(2,2,nn)**2+b(3,2,nn)**2)
1293 xmag(3,nn)=sqrt(b(1,3,nn)**2+b(2,3,nn)**2+b(3,3,nn)**2)
1294C
1295 xmaxv(nn)=max(xmag(1,nn),xmag(2,nn),xmag(3,nn))
1296 ENDDO
1297C
1298#include "vectorize.inc"
1299 DO n = 1, nindex1
1300 nn = index1(n)
1301 IF(xmag(1,nn)==xmaxv(nn)) THEN
1302 lmaxv(nn) = 1
1303 ELSEIF(xmag(2,nn)==xmaxv(nn)) THEN
1304 lmaxv(nn) = 2
1305 ELSE
1306 lmaxv(nn) = 3
1307 ENDIF
1308C
1309 vmag=sqrt(v(1,1,nn)**2+v(2,1,nn)**2)
1310 lmax = lmaxv(nn)
1311 IF(xmaxv(nn)>tol2(nn))THEN
1312 xmaxinv = one/xmaxv(nn)
1313 v(1,3,nn)=b(1,lmax,nn)*xmaxinv
1314 v(2,3,nn)=b(2,lmax,nn)*xmaxinv
1315 v(3,3,nn)=b(3,lmax,nn)*xmaxinv
1316 v(1,2,nn)=v(2,3,nn)*v(3,1,nn)-v(2,1,nn)*v(3,3,nn)
1317 v(2,2,nn)=v(3,3,nn)*v(1,1,nn)-v(3,1,nn)*v(1,3,nn)
1318 v(3,2,nn)=v(1,3,nn)*v(2,1,nn)-v(1,1,nn)*v(2,3,nn)
1319 vmag=one/sqrt(v(1,2,nn)**2+v(2,2,nn)**2+v(3,2,nn)**2)
1320 v(1,2,nn)=v(1,2,nn)*vmag
1321 v(2,2,nn)=v(2,2,nn)*vmag
1322 v(3,2,nn)=v(3,2,nn)*vmag
1323 ELSEIF(vmag>tol2(nn))THEN
1324 v(1,2,nn)=-v(2,1,nn)/vmag
1325 v(2,2,nn)=v(1,1,nn)/vmag
1326 v(3,2,nn)=zero
1327 ELSE
1328 v(1,2,nn)=one
1329 v(2,2,nn)=zero
1330 v(3,2,nn)=zero
1331 ENDIF
1332 ENDDO
1333C . . . . . . . . . . . . .
1334C SOLUTION DOUBLE
1335C . . . . . . . . . . . . .
1336 DO n = 1, nindex2
1337 nn = index2(n)
1338 xmag(1,nn)=sqrt(a(1,1,nn)**2+a(2,1,nn)**2)
1339 xmag(2,nn)=sqrt(a(1,2,nn)**2+a(2,2,nn)**2)
1340 xmag(3,nn)=sqrt(a(1,3,nn)**2+a(2,3,nn)**2)
1341C
1342 xmaxv(nn)=max(xmag(1,nn),xmag(2,nn),xmag(3,nn))
1343 ENDDO
1344C
1345#include "vectorize.inc"
1346 DO n = 1, nindex2
1347 nn = index2(n)
1348 IF(xmag(1,nn)==xmaxv(nn)) THEN
1349 lmaxv(nn) = 1
1350 ELSEIF(xmag(2,nn)==xmaxv(nn)) THEN
1351 lmaxv(nn) = 2
1352 ELSE
1353 lmaxv(nn) = 3
1354 ENDIF
1355C
1356 lmax = lmaxv(nn)
1357 IF(max(abs(a(3,1,nn)),abs(a(3,2,nn)),abs(a(3,3,nn)))
1358 & <tol2(nn))THEN
1359 xmaxinv = one/xmaxv(nn)
1360 v(1,1,nn)= zero
1361 v(2,1,nn)= zero
1362 v(3,1,nn)= one
1363 v(1,2,nn)=-a(2,lmax,nn)*xmaxinv
1364 v(2,2,nn)= a(1,lmax,nn)*xmaxinv
1365 v(3,2,nn)= zero
1366C
1367 ELSEIF(xmaxv(nn)>tol2(nn))THEN
1368 xmaxinv = one/xmaxv(nn)
1369 v(1,1,nn)=-a(2,lmax,nn)*xmaxinv
1370 v(2,1,nn)= a(1,lmax,nn)*xmaxinv
1371 v(3,1,nn)= zero
1372 v(1,2,nn)=-a(3,lmax,nn)*v(2,1,nn)
1373 v(2,2,nn)= a(3,lmax,nn)*v(1,1,nn)
1374 v(3,2,nn)= a(1,lmax,nn)*v(2,1,nn)-a(2,lmax,nn)*v(1,1,nn)
1375 vmag=one/sqrt(v(1,2,nn)**2+v(2,2,nn)**2+v(3,2,nn)**2)
1376 v(1,2,nn)=v(1,2,nn)*vmag
1377 v(2,2,nn)=v(2,2,nn)*vmag
1378 v(3,2,nn)=v(3,2,nn)*vmag
1379 ELSE
1380 v(1,1,nn) = one
1381 v(2,1,nn) = zero
1382 v(3,1,nn) = zero
1383 v(1,2,nn) = zero
1384 v(2,2,nn) = one
1385 v(3,2,nn) = zero
1386 ENDIF
1387 ENDDO
1388C
1389 DO nn = 1, nel
1390 vecdp(1,nn)=v(1,1,nn)
1391 vecdp(2,nn)=v(2,1,nn)
1392 vecdp(3,nn)=v(3,1,nn)
1393 vecdp(4,nn)=v(1,2,nn)
1394 vecdp(5,nn)=v(2,2,nn)
1395 vecdp(6,nn)=v(3,2,nn)
1396 vecdp(7,nn)=vecdp(2,nn)*vecdp(6,nn)-vecdp(3,nn)*vecdp(5,nn)
1397 vecdp(8,nn)=vecdp(3,nn)*vecdp(4,nn)-vecdp(1,nn)*vecdp(6,nn)
1398 vecdp(9,nn)=vecdp(1,nn)*vecdp(5,nn)-vecdp(2,nn)*vecdp(4,nn)
1399 ENDDO
1400C
1401 DO nn = 1, nel
1402 val(1,nn)=valdp(1,nn)
1403 val(2,nn)=valdp(2,nn)
1404 val(3,nn)=valdp(3,nn)
1405 vec(1,nn)=vecdp(1,nn)
1406 vec(2,nn)=vecdp(2,nn)
1407 vec(3,nn)=vecdp(3,nn)
1408 vec(4,nn)=vecdp(4,nn)
1409 vec(5,nn)=vecdp(5,nn)
1410 vec(6,nn)=vecdp(6,nn)
1411 vec(7,nn)=vecdp(7,nn)
1412 vec(8,nn)=vecdp(8,nn)
1413 vec(9,nn)=vecdp(9,nn)
1414 ENDDO
1415C reecriture pour contourner probleme sur itanium2 comp 9. + latency=16
1416 DO n = 1, nindex3
1417 nn = index3(n)
1418C str utilise com tableau temporaire au lieu de scalaires temporaires
1419 str(1,nn)=vec(7,nn)
1420 str(2,nn)=vec(8,nn)
1421 str(3,nn)=vec(9,nn)
1422 ENDDO
1423 DO n = 1, nindex3
1424 nn = index3(n)
1425 vec(7,nn)=vec(1,nn)
1426 vec(8,nn)=vec(2,nn)
1427 vec(9,nn)=vec(3,nn)
1428 vec(1,nn)=-str(1,nn)
1429 vec(2,nn)=-str(2,nn)
1430 vec(3,nn)=-str(3,nn)
1431 ENDDO
1432 RETURN

◆ valpvecdp_v()

subroutine valpvecdp_v ( sig,
val,
vec,
integer nel )

Definition at line 2230 of file sigeps33.F.

2231C-----------------------------------------------
2232C I m p l i c i t T y p e s
2233C-----------------------------------------------
2234#include "implicit_f.inc"
2235C-----------------------------------------------
2236C G l o b a l P a r a m e t e r s
2237C-----------------------------------------------
2238#include "mvsiz_p.inc"
2239C-----------------------------------------------
2240C D u m m y A r g u m e n t s
2241C-----------------------------------------------
2242 my_real
2243 . sig(mvsiz,6), val(mvsiz,3), vec(mvsiz,9)
2244 INTEGER NEL
2245C-----------------------------------------------
2246C L o c a l V a r i a b l e s
2247C-----------------------------------------------
2248 INTEGER I, L, N, NN, LMAX, LMAXV(MVSIZ),NINDEX3, INDEX3(MVSIZ),
2249 . NINDEX1, NINDEX2, INDEX1(MVSIZ), INDEX2(MVSIZ)
2250 double precision
2251 . cs(6), str(mvsiz,3), a(mvsiz,3,3), v(mvsiz,3,3), b(mvsiz,3,3),
2252 . xmag(mvsiz,3), pr, aa, bb, aaa(mvsiz),
2253 . cc, angp, dd, ftpi, ttpi, strmax,
2254 . tol1(mvsiz), tol2(mvsiz), xmaxv(mvsiz), norminf(mvsiz),
2255 . vmag, s11,
2256 . s21, s31, s12, s22, s32, s13, s23, s33, a11, a12, a13, a21,
2257 . a22, a23, a31, a32, a33,
2258 . mdemi, xmaxinv, flm,
2259 . valdp(mvsiz,3),vecdp(mvsiz,9)
2260 REAL FLMIN
2261C-----------------------------------------------
2262C DATA FTPI,TTPI / 4.188790205, 2.094395102 /
2263C FTPI=(4/3)*PI, TTPI=(2/3)*PI
2264C
2265C DEVIATEUR PRINCIPAL DE CONTRAINTE
2266C . . . . . . . . . . . . . . . . . . .
2267 mdemi = -half
2268 ttpi = acos(mdemi)
2269 ftpi = two*ttpi
2270C precision minimum dependant du type REAL ou DOUBLE
2271 CALL floatmin(cs(1),cs(2),flmin)
2272 flm = two*sqrt(flmin)
2273 nindex3=0
2274 DO nn = 1, nel
2275 cs(1) = sig(nn,1)
2276 cs(2) = sig(nn,2)
2277 cs(3) = sig(nn,3)
2278 cs(4) = sig(nn,4)
2279 cs(5) = sig(nn,5)
2280 cs(6) = sig(nn,6)
2281 pr = -(cs(1)+cs(2)+cs(3)) * third
2282 cs(1) = cs(1) + pr
2283 cs(2) = cs(2) + pr
2284 cs(3) = cs(3) + pr
2285 aaa(nn)=cs(4)**2 + cs(5)**2 + cs(6)**2 - cs(1)*cs(2)
2286 & - cs(2)*cs(3) - cs(1)*cs(3)
2287 norminf(nn) = max(abs(cs(1)),abs(cs(2)),abs(cs(3)),
2288 & abs(cs(4)),abs(cs(5)),abs(cs(6)))
2289 norminf(nn) = em10*norminf(nn)
2290C cas racine triple
2291c AA = MAX(AAA(NN),NORMINF(NN),EM20)
2292 aa = max(aaa(nn),norminf(nn))
2293C
2294 bb=cs(1)*cs(5)**2 + cs(2)*cs(6)**2
2295 & + cs(3)*cs(4)**2 - cs(1)*cs(2)*cs(3)
2296 & - two*cs(4)*cs(5)*cs(6)
2297C
2298 cc=-sqrt(twenty7/max(em20,aa))*bb*half/max(em20,aa)
2299 cc= min(cc,one)
2300 cc= max(cc,-one)
2301 angp=acos(cc) * third
2302 dd=two*sqrt(aa * third)
2303 str(nn,1)=dd*cos(angp)
2304 str(nn,2)=dd*cos(angp+ftpi)
2305 str(nn,3)=dd*cos(angp+ttpi)
2306C
2307 valdp(nn,1) = str(nn,1)-pr
2308 valdp(nn,2) = str(nn,2)-pr
2309 valdp(nn,3) = str(nn,3)-pr
2310C renforcement de precision en compression simple----
2311 IF(abs(str(nn,3))>abs(str(nn,1))
2312 & .AND.aaa(nn)>norminf(nn)) THEN
2313 aa=str(nn,1)
2314 str(nn,1)=str(nn,3)
2315 str(nn,3)=aa
2316 nindex3 = nindex3+1
2317 index3(nindex3) = nn
2318 ENDIF
2319C . . . . . . . . . . .
2320C VECTEURS PROPRES
2321C . . . . . . . . . . .
2322 strmax= max(abs(str(nn,1)),abs(str(nn,3)))
2323 tol1(nn)= max(em20,flm*strmax**2)
2324 tol2(nn)=flm*strmax * third
2325 a(nn,1,1)=cs(1)-str(nn,1)
2326 a(nn,2,2)=cs(2)-str(nn,1)
2327 a(nn,3,3)=cs(3)-str(nn,1)
2328 a(nn,1,2)=cs(4)
2329 a(nn,2,1)=cs(4)
2330 a(nn,2,3)=cs(5)
2331 a(nn,3,2)=cs(5)
2332 a(nn,1,3)=cs(6)
2333 a(nn,3,1)=cs(6)
2334C
2335 b(nn,1,1)=a(nn,2,1)*a(nn,3,2)-a(nn,3,1)
2336 . *a(nn,2,2)
2337 b(nn,1,2)=a(nn,2,2)*a(nn,3,3)-a(nn,3,2)
2338 . *a(nn,2,3)
2339 b(nn,1,3)=a(nn,2,3)*a(nn,3,1)-a(nn,3,3)
2340 . *a(nn,2,1)
2341 b(nn,2,1)=a(nn,3,1)*a(nn,1,2)-a(nn,1,1)
2342 . *a(nn,3,2)
2343 b(nn,2,2)=a(nn,3,2)*a(nn,1,3)-a(nn,1,2)
2344 . *a(nn,3,3)
2345 b(nn,2,3)=a(nn,3,3)*a(nn,1,1)-a(nn,1,3)
2346 . *a(nn,3,1)
2347 b(nn,3,1)=a(nn,1,1)*a(nn,2,2)-a(nn,2,1)
2348 . *a(nn,1,2)
2349 b(nn,3,2)=a(nn,1,2)*a(nn,2,3)-a(nn,2,2)
2350 . *a(nn,1,3)
2351 b(nn,3,3)=a(nn,1,3)*a(nn,2,1)-a(nn,2,3)
2352 . *a(nn,1,1)
2353 xmag(nn,1)=sqrt(b(nn,1,1)**2+b(nn,2,1)**2+b(nn,3,1)**2)
2354 xmag(nn,2)=sqrt(b(nn,1,2)**2+b(nn,2,2)**2+b(nn,3,2)**2)
2355 xmag(nn,3)=sqrt(b(nn,1,3)**2+b(nn,2,3)**2+b(nn,3,3)**2)
2356
2357 ENDDO
2358C
2359 nindex1 = 0
2360 nindex2 = 0
2361 DO nn = 1, nel
2362 xmaxv(nn)=max(xmag(nn,1),xmag(nn,2),xmag(nn,3))
2363 IF(xmag(nn,1)==xmaxv(nn)) THEN
2364 lmaxv(nn) = 1
2365 ELSEIF(xmag(nn,2)==xmaxv(nn)) THEN
2366 lmaxv(nn) = 2
2367 ELSE
2368 lmaxv(nn) = 3
2369 ENDIF
2370 IF(aaa(nn)<norminf(nn)) THEN
2371 valdp(nn,1) = sig(nn,1)
2372 valdp(nn,2) = sig(nn,2)
2373 valdp(nn,3) = sig(nn,3)
2374 v(nn,1,1) = one
2375 v(nn,2,1) = zero
2376 v(nn,3,1) = zero
2377 v(nn,1,2) = zero
2378 v(nn,2,2) = one
2379 v(nn,3,2) = zero
2380 ELSEIF(xmaxv(nn)>tol1(nn)) THEN
2381 nindex1 = nindex1 + 1
2382 index1(nindex1) = nn
2383 ELSE
2384 nindex2 = nindex2 + 1
2385 index2(nindex2) = nn
2386 ENDIF
2387 ENDDO
2388C
2389#include "vectorize.inc"
2390 DO n = 1, nindex1
2391 nn = index1(n)
2392 lmax = lmaxv(nn)
2393 xmaxinv = one/xmaxv(nn)
2394 v(nn,1,1)=b(nn,1,lmax)*xmaxinv
2395 v(nn,2,1)=b(nn,2,lmax)*xmaxinv
2396 v(nn,3,1)=b(nn,3,lmax)*xmaxinv
2397 a(nn,1,1)=a(nn,1,1)+str(nn,1)-str(nn,3)
2398 a(nn,2,2)=a(nn,2,2)+str(nn,1)-str(nn,3)
2399 a(nn,3,3)=a(nn,3,3)+str(nn,1)-str(nn,3)
2400C
2401 b(nn,1,1)=a(nn,2,1)*v(nn,3,1)-a(nn,3,1)*v(nn,2,1)
2402 b(nn,1,2)=a(nn,2,2)*v(nn,3,1)-a(nn,3,2)*v(nn,2,1)
2403 b(nn,1,3)=a(nn,2,3)*v(nn,3,1)-a(nn,3,3)*v(nn,2,1)
2404 b(nn,2,1)=a(nn,3,1)*v(nn,1,1)-a(nn,1,1)*v(nn,3,1)
2405 b(nn,2,2)=a(nn,3,2)*v(nn,1,1)-a(nn,1,2)*v(nn,3,1)
2406 b(nn,2,3)=a(nn,3,3)*v(nn,1,1)-a(nn,1,3)*v(nn,3,1)
2407 b(nn,3,1)=a(nn,1,1)*v(nn,2,1)-a(nn,2,1)*v(nn,1,1)
2408 b(nn,3,2)=a(nn,1,2)*v(nn,2,1)-a(nn,2,2)*v(nn,1,1)
2409 b(nn,3,3)=a(nn,1,3)*v(nn,2,1)-a(nn,2,3)*v(nn,1,1)
2410 xmag(nn,1)=sqrt(b(nn,1,1)**2+b(nn,2,1)**2+b(nn,3,1)**2)
2411 xmag(nn,2)=sqrt(b(nn,1,2)**2+b(nn,2,2)**2+b(nn,3,2)**2)
2412 xmag(nn,3)=sqrt(b(nn,1,3)**2+b(nn,2,3)**2+b(nn,3,3)**2)
2413C
2414 xmaxv(nn)=max(xmag(nn,1),xmag(nn,2),xmag(nn,3))
2415 ENDDO
2416C
2417#include "vectorize.inc"
2418 DO n = 1, nindex1
2419 nn = index1(n)
2420 IF(xmag(nn,1)==xmaxv(nn)) THEN
2421 lmaxv(nn) = 1
2422 ELSEIF(xmag(nn,2)==xmaxv(nn)) THEN
2423 lmaxv(nn) = 2
2424 ELSE
2425 lmaxv(nn) = 3
2426 ENDIF
2427C
2428 vmag=sqrt(v(nn,1,1)**2+v(nn,2,1)**2)
2429 lmax = lmaxv(nn)
2430 IF(xmaxv(nn)>tol2(nn))THEN
2431 xmaxinv = one/xmaxv(nn)
2432 v(nn,1,3)=b(nn,1,lmax)*xmaxinv
2433 v(nn,2,3)=b(nn,2,lmax)*xmaxinv
2434 v(nn,3,3)=b(nn,3,lmax)*xmaxinv
2435 v(nn,1,2)=v(nn,2,3)*v(nn,3,1)-v(nn,2,1)*v(nn,3,3)
2436 v(nn,2,2)=v(nn,3,3)*v(nn,1,1)-v(nn,3,1)*v(nn,1,3)
2437 v(nn,3,2)=v(nn,1,3)*v(nn,2,1)-v(nn,1,1)*v(nn,2,3)
2438 vmag=one/sqrt(v(nn,1,2)**2+v(nn,2,2)**2+v(nn,3,2)**2)
2439 v(nn,1,2)=v(nn,1,2)*vmag
2440 v(nn,2,2)=v(nn,2,2)*vmag
2441 v(nn,3,2)=v(nn,3,2)*vmag
2442 ELSEIF(vmag>tol2(nn))THEN
2443 v(nn,1,2)=-v(nn,2,1)/vmag
2444 v(nn,2,2)=v(nn,1,1)/vmag
2445 v(nn,3,2)=zero
2446 ELSE
2447 v(nn,1,2)=one
2448 v(nn,2,2)=zero
2449 v(nn,3,2)=zero
2450 ENDIF
2451 ENDDO
2452C . . . . . . . . . . . . .
2453C SOLUTION DOUBLE
2454C . . . . . . . . . . . . .
2455#include "vectorize.inc"
2456 DO n = 1, nindex2
2457 nn = index2(n)
2458 xmag(nn,1)=sqrt(a(nn,1,1)**2+a(nn,2,1)**2)
2459 xmag(nn,2)=sqrt(a(nn,1,2)**2+a(nn,2,2)**2)
2460 xmag(nn,3)=sqrt(a(nn,1,3)**2+a(nn,2,3)**2)
2461C
2462 xmaxv(nn)=max(xmag(nn,1),xmag(nn,2),xmag(nn,3))
2463 ENDDO
2464C
2465#include "vectorize.inc"
2466 DO n = 1, nindex2
2467 nn = index2(n)
2468 IF(xmag(nn,1)==xmaxv(nn)) THEN
2469 lmaxv(nn) = 1
2470 ELSEIF(xmag(nn,2)==xmaxv(nn)) THEN
2471 lmaxv(nn) = 2
2472 ELSE
2473 lmaxv(nn) = 3
2474 ENDIF
2475C
2476 lmax = lmaxv(nn)
2477 IF(max(abs(a(nn,3,1)),abs(a(nn,3,2)),abs(a(nn,3,3)))
2478 & <tol2(nn))THEN
2479 xmaxinv = one/xmaxv(nn)
2480 v(nn,1,1)= zero
2481 v(nn,2,1)= zero
2482 v(nn,3,1)= one
2483 v(nn,1,2)=-a(nn,2,lmax)*xmaxinv
2484 v(nn,2,2)= a(nn,1,lmax)*xmaxinv
2485 v(nn,3,2)= zero
2486C
2487 ELSEIF(xmaxv(nn)>tol2(nn))THEN
2488 xmaxinv = one/xmaxv(nn)
2489 v(nn,1,1)=-a(nn,2,lmax)*xmaxinv
2490 v(nn,2,1)= a(nn,1,lmax)*xmaxinv
2491 v(nn,3,1)= zero
2492 v(nn,1,2)=-a(nn,3,lmax)*v(nn,2,1)
2493 v(nn,2,2)= a(nn,3,lmax)*v(nn,1,1)
2494 v(nn,3,2)= a(nn,1,lmax)*v(nn,2,1)-a(nn,2,lmax)*v(nn,1,1)
2495 vmag=one/sqrt(v(nn,1,2)**2+v(nn,2,2)**2+v(nn,3,2)**2)
2496 v(nn,1,2)=v(nn,1,2)*vmag
2497 v(nn,2,2)=v(nn,2,2)*vmag
2498 v(nn,3,2)=v(nn,3,2)*vmag
2499 ELSE
2500 v(nn,1,1) = one
2501 v(nn,2,1) = zero
2502 v(nn,3,1) = zero
2503 v(nn,1,2) = zero
2504 v(nn,2,2) = one
2505 v(nn,3,2) = zero
2506 ENDIF
2507 ENDDO
2508C
2509 DO nn = 1, nel
2510 vecdp(nn,1)=v(nn,1,1)
2511 vecdp(nn,2)=v(nn,2,1)
2512 vecdp(nn,3)=v(nn,3,1)
2513 vecdp(nn,4)=v(nn,1,2)
2514 vecdp(nn,5)=v(nn,2,2)
2515 vecdp(nn,6)=v(nn,3,2)
2516 vecdp(nn,7)=vecdp(nn,2)*vecdp(nn,6)-vecdp(nn,3)*vecdp(nn,5)
2517 vecdp(nn,8)=vecdp(nn,3)*vecdp(nn,4)-vecdp(nn,1)*vecdp(nn,6)
2518 vecdp(nn,9)=vecdp(nn,1)*vecdp(nn,5)-vecdp(nn,2)*vecdp(nn,4)
2519 ENDDO
2520C
2521 DO nn = 1, nel
2522 val(nn,1)=valdp(nn,1)
2523 val(nn,2)=valdp(nn,2)
2524 val(nn,3)=valdp(nn,3)
2525 vec(nn,1)=vecdp(nn,1)
2526 vec(nn,2)=vecdp(nn,2)
2527 vec(nn,3)=vecdp(nn,3)
2528 vec(nn,4)=vecdp(nn,4)
2529 vec(nn,5)=vecdp(nn,5)
2530 vec(nn,6)=vecdp(nn,6)
2531 vec(nn,7)=vecdp(nn,7)
2532 vec(nn,8)=vecdp(nn,8)
2533 vec(nn,9)=vecdp(nn,9)
2534 ENDDO
2535C reecriture pour contourner probleme sur itanium2 comp 9. + latency=16
2536#include "vectorize.inc"
2537 DO n = 1, nindex3
2538 nn = index3(n)
2539C str utilise com tableau temporaire au lieu de scalaires temporaires
2540 str(nn,1)=vec(nn,7)
2541 str(nn,2)=vec(nn,8)
2542 str(nn,3)=vec(nn,9)
2543 ENDDO
2544#include "vectorize.inc"
2545 DO n = 1, nindex3
2546 nn = index3(n)
2547 vec(nn,7)=vec(nn,1)
2548 vec(nn,8)=vec(nn,2)
2549 vec(nn,9)=vec(nn,3)
2550 vec(nn,1)=-str(nn,1)
2551 vec(nn,2)=-str(nn,2)
2552 vec(nn,3)=-str(nn,3)
2553 ENDDO
2554 RETURN

◆ valpvecop()

subroutine valpvecop ( sig,
val,
vec,
integer nel )

Definition at line 1539 of file sigeps33.F.

1540C-----------------------------------------------
1541C I m p l i c i t T y p e s
1542C-----------------------------------------------
1543#include "implicit_f.inc"
1544C-----------------------------------------------
1545C G l o b a l P a r a m e t e r s
1546C-----------------------------------------------
1547#include "mvsiz_p.inc"
1548C-----------------------------------------------
1549C D u m m y A r g u m e n t s
1550C-----------------------------------------------
1551 my_real
1552 . sig(6,*), val(3,*), vec(9,*)
1553 INTEGER NEL
1554C-----------------------------------------------
1555C L o c a l V a r i a b l e s
1556C-----------------------------------------------
1557 INTEGER I, L, N, NN, LMAX, LMAXV(MVSIZ),NINDEX3, INDEX3(MVSIZ),
1558 . NINDEX1, NINDEX2, INDEX1(MVSIZ), INDEX2(MVSIZ)
1559 double precision
1560 . cs(6), str(3,mvsiz), a(3,3,mvsiz), v(3,3,mvsiz), b(3,3,mvsiz),
1561 . xmag(3,mvsiz), pr, aa, bb, aaa(mvsiz),
1562 . cc, angp, dd, ftpi, ttpi, strmax,
1563 . tol1(mvsiz), tol2(mvsiz), xmaxv(mvsiz), norminf(mvsiz),
1564 . vmag, s11,
1565 . s21, s31, s12, s22, s32, s13, s23, s33, a11, a12, a13, a21,
1566 . a22, a23, a31, a32, a33,
1567 . mdemi, xmaxinv, flm,
1568 . valdp(3,mvsiz),vecdp(9,mvsiz),s4_2,s5_2,s6_2,c1,c2,sq_aa
1569 REAL FLMIN
1570C-----------------------------------------------
1571C DATA FTPI,TTPI / 4.188790205, 2.094395102 /
1572C FTPI=(4/3)*PI, TTPI=(2/3)*PI
1573C
1574C DEVIATEUR PRINCIPAL DE CONTRAINTE
1575C . . . . . . . . . . . . . . . . . . .
1576 mdemi = -half
1577 ttpi = acos(mdemi)
1578 ftpi = two*ttpi
1579C precision minimum dependant du type REAL ou DOUBLE
1580 CALL floatmin(cs(1),cs(2),flmin)
1581 flm = two*sqrt(flmin)
1582 nindex3=0
1583 c1 = half*sqrt(twenty7)
1584 c2 = two*sqrt(third)
1585 DO nn = 1, nel
1586 cs(1) = sig(1,nn)
1587 cs(2) = sig(2,nn)
1588 cs(3) = sig(3,nn)
1589 cs(4) = sig(4,nn)
1590 cs(5) = sig(5,nn)
1591 cs(6) = sig(6,nn)
1592 pr = -(cs(1)+cs(2)+cs(3)) * third
1593 cs(1) = cs(1) + pr
1594 cs(2) = cs(2) + pr
1595 cs(3) = cs(3) + pr
1596 s4_2 = cs(4)*cs(4)
1597 s5_2 = cs(5)*cs(5)
1598 s6_2 = cs(6)*cs(6)
1599 aaa(nn)=s4_2 + s5_2 + s6_2 - cs(1)*cs(2)
1600 & - cs(2)*cs(3) - cs(1)*cs(3)
1601 norminf(nn) = max(abs(cs(1)),abs(cs(2)),abs(cs(3)),
1602 & abs(cs(4)),abs(cs(5)),abs(cs(6)))
1603 norminf(nn) = em10*norminf(nn)
1604C cas racine triple
1605 aa = max(aaa(nn),norminf(nn))
1606C
1607 bb=cs(1)*s5_2 + cs(2)*s6_2
1608 & + cs(3)*s4_2 - cs(1)*cs(2)*cs(3)
1609 & - two*cs(4)*cs(5)*cs(6)
1610C
1611 sq_aa = sqrt(aa)
1612 cc=-c1/max(em10,sq_aa)*bb/max(em20,aa)
1613c CC=-SQRT(TWENTY7/MAX(EM20,AA))*BB*HALF/MAX(EM20,AA)
1614 cc= min(cc,one)
1615 cc= max(cc,-one)
1616 angp=acos(cc) * third
1617 dd=c2*sq_aa
1618 str(1,nn)=dd*cos(angp)
1619 str(2,nn)=dd*cos(angp+ftpi)
1620 str(3,nn)=dd*cos(angp+ttpi)
1621C
1622 valdp(1,nn) = str(1,nn)-pr
1623 valdp(2,nn) = str(2,nn)-pr
1624 valdp(3,nn) = str(3,nn)-pr
1625C renforcement de precision en compression simple----
1626 IF(abs(str(3,nn))>abs(str(1,nn))
1627 & .AND.aaa(nn)>norminf(nn)) THEN
1628 aa=str(1,nn)
1629 str(1,nn)=str(3,nn)
1630 str(3,nn)=aa
1631 nindex3 = nindex3+1
1632 index3(nindex3) = nn
1633 ENDIF
1634C . . . . . . . . . . .
1635C VECTEURS PROPRES
1636C . . . . . . . . . . .
1637 strmax= max(abs(str(1,nn)),abs(str(3,nn)))
1638 tol1(nn)= max(em20,flm*strmax*strmax)
1639 tol2(nn)=flm*strmax * third
1640 a(1,1,nn)=cs(1)-str(1,nn)
1641 a(2,2,nn)=cs(2)-str(1,nn)
1642 a(3,3,nn)=cs(3)-str(1,nn)
1643 a(1,2,nn)=cs(4)
1644 a(2,1,nn)=cs(4)
1645 a(2,3,nn)=cs(5)
1646 a(3,2,nn)=cs(5)
1647 a(1,3,nn)=cs(6)
1648 a(3,1,nn)=cs(6)
1649C
1650 b(1,1,nn)=a(2,1,nn)*a(3,2,nn)-a(3,1,nn)
1651 . *a(2,2,nn)
1652 b(1,2,nn)=a(2,2,nn)*a(3,3,nn)-a(3,2,nn)
1653 . *a(2,3,nn)
1654 b(1,3,nn)=a(2,3,nn)*a(3,1,nn)-a(3,3,nn)
1655 . *a(2,1,nn)
1656 b(2,1,nn)=a(3,1,nn)*a(1,2,nn)-a(1,1,nn)
1657 . *a(3,2,nn)
1658 b(2,2,nn)=a(3,2,nn)*a(1,3,nn)-a(1,2,nn)
1659 . *a(3,3,nn)
1660 b(2,3,nn)=a(3,3,nn)*a(1,1,nn)-a(1,3,nn)
1661 . *a(3,1,nn)
1662 b(3,1,nn)=a(1,1,nn)*a(2,2,nn)-a(2,1,nn)
1663 . *a(1,2,nn)
1664 b(3,2,nn)=a(1,2,nn)*a(2,3,nn)-a(2,2,nn)
1665 . *a(1,3,nn)
1666 b(3,3,nn)=a(1,3,nn)*a(2,1,nn)-a(2,3,nn)
1667 . *a(1,1,nn)
1668 xmag(1,nn)=b(1,1,nn)*b(1,1,nn)+b(2,1,nn)*b(2,1,nn)+
1669 . b(3,1,nn)*b(3,1,nn)
1670 xmag(2,nn)=b(1,2,nn)*b(1,2,nn)+b(2,2,nn)*b(2,2,nn)+
1671 . b(3,2,nn)*b(3,2,nn)
1672 xmag(3,nn)=b(1,3,nn)*b(1,3,nn)+b(2,3,nn)*b(2,3,nn)+
1673 . b(3,3,nn)*b(3,3,nn)
1674
1675 ENDDO
1676C
1677 nindex1 = 0
1678 nindex2 = 0
1679 DO nn = 1, nel
1680 xmaxv(nn)=max(xmag(1,nn),xmag(2,nn),xmag(3,nn))
1681 IF(xmag(1,nn)==xmaxv(nn)) THEN
1682 lmaxv(nn) = 1
1683 ELSEIF(xmag(2,nn)==xmaxv(nn)) THEN
1684 lmaxv(nn) = 2
1685 ELSE
1686 lmaxv(nn) = 3
1687 ENDIF
1688 IF(aaa(nn)<norminf(nn)) THEN
1689 valdp(1,nn) = sig(1,nn)
1690 valdp(2,nn) = sig(2,nn)
1691 valdp(3,nn) = sig(3,nn)
1692 v(1,1,nn) = one
1693 v(2,1,nn) = zero
1694 v(3,1,nn) = zero
1695 v(1,2,nn) = zero
1696 v(2,2,nn) = one
1697 v(3,2,nn) = zero
1698 ELSEIF(xmaxv(nn)>tol1(nn)*tol1(nn)) THEN
1699 nindex1 = nindex1 + 1
1700 index1(nindex1) = nn
1701 ELSE
1702 nindex2 = nindex2 + 1
1703 index2(nindex2) = nn
1704 ENDIF
1705 ENDDO
1706C
1707#include "vectorize.inc"
1708 DO n = 1, nindex1
1709 nn = index1(n)
1710 lmax = lmaxv(nn)
1711 xmaxinv = one/sqrt(xmaxv(nn))
1712 v(1,1,nn)=b(1,lmax,nn)*xmaxinv
1713 v(2,1,nn)=b(2,lmax,nn)*xmaxinv
1714 v(3,1,nn)=b(3,lmax,nn)*xmaxinv
1715 a(1,1,nn)=a(1,1,nn)+str(1,nn)-str(3,nn)
1716 a(2,2,nn)=a(2,2,nn)+str(1,nn)-str(3,nn)
1717 a(3,3,nn)=a(3,3,nn)+str(1,nn)-str(3,nn)
1718C
1719 b(1,1,nn)=a(2,1,nn)*v(3,1,nn)-a(3,1,nn)*v(2,1,nn)
1720 b(1,2,nn)=a(2,2,nn)*v(3,1,nn)-a(3,2,nn)*v(2,1,nn)
1721 b(1,3,nn)=a(2,3,nn)*v(3,1,nn)-a(3,3,nn)*v(2,1,nn)
1722 b(2,1,nn)=a(3,1,nn)*v(1,1,nn)-a(1,1,nn)*v(3,1,nn)
1723 b(2,2,nn)=a(3,2,nn)*v(1,1,nn)-a(1,2,nn)*v(3,1,nn)
1724 b(2,3,nn)=a(3,3,nn)*v(1,1,nn)-a(1,3,nn)*v(3,1,nn)
1725 b(3,1,nn)=a(1,1,nn)*v(2,1,nn)-a(2,1,nn)*v(1,1,nn)
1726 b(3,2,nn)=a(1,2,nn)*v(2,1,nn)-a(2,2,nn)*v(1,1,nn)
1727 b(3,3,nn)=a(1,3,nn)*v(2,1,nn)-a(2,3,nn)*v(1,1,nn)
1728 xmag(1,nn)=b(1,1,nn)*b(1,1,nn)+b(2,1,nn)*b(2,1,nn)+
1729 . b(3,1,nn)*b(3,1,nn)
1730 xmag(2,nn)=b(1,2,nn)*b(1,2,nn)+b(2,2,nn)*b(2,2,nn)+
1731 . b(3,2,nn)*b(3,2,nn)
1732 xmag(3,nn)=b(1,3,nn)*b(1,3,nn)+b(2,3,nn)*b(2,3,nn)+
1733 . b(3,3,nn)*b(3,3,nn)
1734C
1735 xmaxv(nn)=max(xmag(1,nn),xmag(2,nn),xmag(3,nn))
1736 ENDDO
1737C
1738#include "vectorize.inc"
1739 DO n = 1, nindex1
1740 nn = index1(n)
1741 IF(xmag(1,nn)==xmaxv(nn)) THEN
1742 lmaxv(nn) = 1
1743 ELSEIF(xmag(2,nn)==xmaxv(nn)) THEN
1744 lmaxv(nn) = 2
1745 ELSE
1746 lmaxv(nn) = 3
1747 ENDIF
1748C
1749 vmag=sqrt(v(1,1,nn)**2+v(2,1,nn)**2)
1750 lmax = lmaxv(nn)
1751 xmaxv(nn)=sqrt(xmaxv(nn))
1752 IF(xmaxv(nn)>tol2(nn))THEN
1753 xmaxinv = one/xmaxv(nn)
1754 v(1,3,nn)=b(1,lmax,nn)*xmaxinv
1755 v(2,3,nn)=b(2,lmax,nn)*xmaxinv
1756 v(3,3,nn)=b(3,lmax,nn)*xmaxinv
1757 v(1,2,nn)=v(2,3,nn)*v(3,1,nn)-v(2,1,nn)*v(3,3,nn)
1758 v(2,2,nn)=v(3,3,nn)*v(1,1,nn)-v(3,1,nn)*v(1,3,nn)
1759 v(3,2,nn)=v(1,3,nn)*v(2,1,nn)-v(1,1,nn)*v(2,3,nn)
1760 c1 = v(1,2,nn)*v(1,2,nn)+v(2,2,nn)*v(2,2,nn)+
1761 & v(3,2,nn)*v(3,2,nn)
1762 vmag=one/sqrt(c1)
1763 v(1,2,nn)=v(1,2,nn)*vmag
1764 v(2,2,nn)=v(2,2,nn)*vmag
1765 v(3,2,nn)=v(3,2,nn)*vmag
1766 ELSEIF(vmag>tol2(nn))THEN
1767 v(1,2,nn)=-v(2,1,nn)/vmag
1768 v(2,2,nn)=v(1,1,nn)/vmag
1769 v(3,2,nn)=zero
1770 ELSE
1771 v(1,2,nn)=one
1772 v(2,2,nn)=zero
1773 v(3,2,nn)=zero
1774 ENDIF
1775 ENDDO
1776C . . . . . . . . . . . . .
1777C SOLUTION DOUBLE
1778C . . . . . . . . . . . . .
1779 DO n = 1, nindex2
1780 nn = index2(n)
1781 xmag(1,nn)=a(1,1,nn)*a(1,1,nn)+a(2,1,nn)*a(2,1,nn)
1782 xmag(2,nn)=a(1,2,nn)*a(1,2,nn)+a(2,2,nn)*a(2,2,nn)
1783 xmag(3,nn)=a(1,3,nn)*a(1,3,nn)+a(2,3,nn)*a(2,3,nn)
1784C
1785 xmaxv(nn)=max(xmag(1,nn),xmag(2,nn),xmag(3,nn))
1786 ENDDO
1787C
1788#include "vectorize.inc"
1789 DO n = 1, nindex2
1790 nn = index2(n)
1791 IF(xmag(1,nn)==xmaxv(nn)) THEN
1792 lmaxv(nn) = 1
1793 ELSEIF(xmag(2,nn)==xmaxv(nn)) THEN
1794 lmaxv(nn) = 2
1795 ELSE
1796 lmaxv(nn) = 3
1797 ENDIF
1798C
1799 lmax = lmaxv(nn)
1800 xmaxv(nn)=sqrt(xmaxv(nn))
1801 IF(max(abs(a(3,1,nn)),abs(a(3,2,nn)),abs(a(3,3,nn)))
1802 & <tol2(nn))THEN
1803 xmaxinv = one/xmaxv(nn)
1804 v(1,1,nn)= zero
1805 v(2,1,nn)= zero
1806 v(3,1,nn)= one
1807 v(1,2,nn)=-a(2,lmax,nn)*xmaxinv
1808 v(2,2,nn)= a(1,lmax,nn)*xmaxinv
1809 v(3,2,nn)= zero
1810C
1811 ELSEIF(xmaxv(nn)>tol2(nn))THEN
1812 xmaxinv = one/xmaxv(nn)
1813 v(1,1,nn)=-a(2,lmax,nn)*xmaxinv
1814 v(2,1,nn)= a(1,lmax,nn)*xmaxinv
1815 v(3,1,nn)= zero
1816 v(1,2,nn)=-a(3,lmax,nn)*v(2,1,nn)
1817 v(2,2,nn)= a(3,lmax,nn)*v(1,1,nn)
1818 v(3,2,nn)= a(1,lmax,nn)*v(2,1,nn)-a(2,lmax,nn)*v(1,1,nn)
1819 c1 = v(1,2,nn)*v(1,2,nn)+v(2,2,nn)*v(2,2,nn)+
1820 & v(3,2,nn)*v(3,2,nn)
1821 vmag=one/sqrt(c1)
1822 v(1,2,nn)=v(1,2,nn)*vmag
1823 v(2,2,nn)=v(2,2,nn)*vmag
1824 v(3,2,nn)=v(3,2,nn)*vmag
1825 ELSE
1826 v(1,1,nn) = one
1827 v(2,1,nn) = zero
1828 v(3,1,nn) = zero
1829 v(1,2,nn) = zero
1830 v(2,2,nn) = one
1831 v(3,2,nn) = zero
1832 ENDIF
1833 ENDDO
1834C
1835 DO nn = 1, nel
1836 vecdp(1,nn)=v(1,1,nn)
1837 vecdp(2,nn)=v(2,1,nn)
1838 vecdp(3,nn)=v(3,1,nn)
1839 vecdp(4,nn)=v(1,2,nn)
1840 vecdp(5,nn)=v(2,2,nn)
1841 vecdp(6,nn)=v(3,2,nn)
1842 vecdp(7,nn)=vecdp(2,nn)*vecdp(6,nn)-vecdp(3,nn)*vecdp(5,nn)
1843 vecdp(8,nn)=vecdp(3,nn)*vecdp(4,nn)-vecdp(1,nn)*vecdp(6,nn)
1844 vecdp(9,nn)=vecdp(1,nn)*vecdp(5,nn)-vecdp(2,nn)*vecdp(4,nn)
1845 ENDDO
1846C
1847 DO nn = 1, nel
1848 val(1,nn)=valdp(1,nn)
1849 val(2,nn)=valdp(2,nn)
1850 val(3,nn)=valdp(3,nn)
1851 vec(1,nn)=vecdp(1,nn)
1852 vec(2,nn)=vecdp(2,nn)
1853 vec(3,nn)=vecdp(3,nn)
1854 vec(4,nn)=vecdp(4,nn)
1855 vec(5,nn)=vecdp(5,nn)
1856 vec(6,nn)=vecdp(6,nn)
1857 vec(7,nn)=vecdp(7,nn)
1858 vec(8,nn)=vecdp(8,nn)
1859 vec(9,nn)=vecdp(9,nn)
1860 ENDDO
1861C reecriture pour contourner probleme sur itanium2 comp 9. + latency=16
1862 DO n = 1, nindex3
1863 nn = index3(n)
1864C str utilise com tableau temporaire au lieu de scalaires temporaires
1865 str(1,nn)=vec(7,nn)
1866 str(2,nn)=vec(8,nn)
1867 str(3,nn)=vec(9,nn)
1868 ENDDO
1869 DO n = 1, nindex3
1870 nn = index3(n)
1871 vec(7,nn)=vec(1,nn)
1872 vec(8,nn)=vec(2,nn)
1873 vec(9,nn)=vec(3,nn)
1874 vec(1,nn)=-str(1,nn)
1875 vec(2,nn)=-str(2,nn)
1876 vec(3,nn)=-str(3,nn)
1877 ENDDO
1878 RETURN

◆ valpvecop_v()

subroutine valpvecop_v ( sig,
val,
vec,
integer nel )

Definition at line 2564 of file sigeps33.F.

2565C-----------------------------------------------
2566C I m p l i c i t T y p e s
2567C-----------------------------------------------
2568#include "implicit_f.inc"
2569C-----------------------------------------------
2570C G l o b a l P a r a m e t e r s
2571C-----------------------------------------------
2572#include "mvsiz_p.inc"
2573C-----------------------------------------------
2574C D u m m y A r g u m e n t s
2575C-----------------------------------------------
2576 my_real
2577 . sig(mvsiz,6), val(mvsiz,3), vec(mvsiz,9)
2578 INTEGER NEL
2579C-----------------------------------------------
2580C L o c a l V a r i a b l e s
2581C-----------------------------------------------
2582 INTEGER I, L, N, NN, LMAX, LMAXV(MVSIZ),NINDEX3, INDEX3(MVSIZ),
2583 . NINDEX1, NINDEX2, INDEX1(MVSIZ), INDEX2(MVSIZ)
2584 double precision
2585 . cs1,cs2,cs3,cs4,cs5,cs6,
2586 . a(mvsiz,3,3), v(mvsiz,3,3), b(mvsiz,3,3),
2587 . pr, aa, bb, aaa(mvsiz),
2588 . cc, angp, dd, ftpi, ttpi, strmax,
2589 . tol1(mvsiz), tol2(mvsiz), xmaxv(mvsiz), norminf(mvsiz),
2590 . vmag, s11,
2591 . s21, s31, s12, s22, s32, s13, s23, s33, a11, a12, a13, a21,
2592 . s1,s2,s3,
2593 . a22, a23, a31, a32, a33,
2594 . mdemi, xmaxinv, flm,
2595 . s4_2,s5_2,s6_2,c1,c2,sq_aa,
2596 . valdp1(mvsiz),valdp2(mvsiz),valdp3(mvsiz),
2597 . str1(mvsiz),str2(mvsiz),str3(mvsiz),
2598 . xmag1(mvsiz),xmag2(mvsiz),xmag3(mvsiz),
2599 . vecdp1(mvsiz),vecdp2(mvsiz),vecdp3(mvsiz),
2600 . vecdp4(mvsiz),vecdp5(mvsiz),vecdp6(mvsiz),
2601 . vecdp7(mvsiz),vecdp8(mvsiz),vecdp9(mvsiz)
2602
2603 DOUBLE PRECISION :: SQRT3DMI, SIN_ANGP, COS_ANGP
2604 parameter(sqrt3dmi=sqrt(3.0d0) / 2.0d0 )
2605
2606 REAL FLMIN
2607C-----------------------------------------------
2608C DATA FTPI,TTPI / 4.188790205, 2.094395102 /
2609C FTPI=(4/3)*PI, TTPI=(2/3)*PI
2610C
2611C DEVIATEUR PRINCIPAL DE CONTRAINTE
2612C . . . . . . . . . . . . . . . . . . .
2613 mdemi = -half
2614 ttpi = acos(mdemi)
2615 ftpi = two*ttpi
2616C precision minimum dependant du type REAL ou DOUBLE
2617 CALL floatmin(cs1,cs2,flmin)
2618 flm = two*sqrt(flmin)
2619 nindex3=0
2620 index3(1:nel) =0
2621 c1 = half*sqrt(twenty7)
2622 c2 = two*sqrt(third)
2623 DO nn = 1, nel
2624 cs1 = sig(nn,1)
2625 cs2 = sig(nn,2)
2626 cs3 = sig(nn,3)
2627 cs4 = sig(nn,4)
2628 cs5 = sig(nn,5)
2629 cs6 = sig(nn,6)
2630 pr = -(cs1+cs2+cs3) * third
2631 cs1 = cs1 + pr
2632 cs2 = cs2 + pr
2633 cs3 = cs3 + pr
2634 s4_2 = cs4*cs4
2635 s5_2 = cs5*cs5
2636 s6_2 = cs6*cs6
2637 aaa(nn)=s4_2 + s5_2 + s6_2 - cs1*cs2
2638 & - cs2*cs3 - cs1*cs3
2639 norminf(nn) = max(abs(cs1),abs(cs2),abs(cs3),
2640 & abs(cs4),abs(cs5),abs(cs6))
2641 norminf(nn) = em10*norminf(nn)
2642C cas racine triple
2643 aa = max(aaa(nn),norminf(nn))
2644 sq_aa = sqrt(aa)
2645C
2646 bb=cs1*s5_2 + cs2*s6_2
2647 & + cs3*s4_2 - cs1*cs2*cs3
2648 & - two*cs4*cs5*cs6
2649C
2650 cc=-c1/max(em10,sq_aa)*bb/max(em20,aa)
2651c CC=-SQRT(TWENTY7/MAX(EM20,AA))*BB*HALF/MAX(EM20,AA)
2652 cc= min(cc,one)
2653 cc= max(cc,-one)
2654 angp=acos(cc) * third
2655 dd=c2*sq_aa
2656
2657 cos_angp = cos(angp)
2658 sin_angp = sin(angp)
2659c STR1(NN)=DD*COS_ANGP
2660c STR2(NN)=DD*((-HALF*COS_ANGP) - (SIN_ANGP*SQRT3DMI))
2661c STR3(NN)=DD*((-HALF*COS_ANGP) + (SIN_ANGP*SQRT3DMI))
2662 str3(nn)=dd*cos_angp
2663 str1(nn)=dd*((-half*cos_angp) - (sin_angp*sqrt3dmi))
2664 str2(nn)=-(str1(nn)+str3(nn))
2665C
2666 valdp1(nn) = str1(nn)-pr
2667 valdp2(nn) = str2(nn)-pr
2668 valdp3(nn) = str3(nn)-pr
2669C renforcement de precision en compression simple----
2670 IF(abs(str3(nn))>abs(str1(nn))
2671 & .AND.aaa(nn)>norminf(nn)) THEN
2672 aa=str1(nn)
2673 str1(nn)=str3(nn)
2674 str3(nn)=aa
2675 index3(nn) = 1
2676 nindex3 = 1
2677 ENDIF
2678
2679C . . . . . . . . . . .
2680C VECTEURS PROPRES
2681C . . . . . . . . . . .
2682 strmax= max(abs(str1(nn)),abs(str3(nn)))
2683 tol1(nn)= max(em20,flm*strmax*strmax)
2684 tol2(nn)=flm*strmax * third
2685 a(nn,1,1)=cs1-str1(nn)
2686 a(nn,2,2)=cs2-str1(nn)
2687 a(nn,3,3)=cs3-str1(nn)
2688 a(nn,1,2)=cs4
2689 a(nn,2,1)=cs4
2690 a(nn,2,3)=cs5
2691 a(nn,3,2)=cs5
2692 a(nn,1,3)=cs6
2693 a(nn,3,1)=cs6
2694C
2695 b(nn,1,1)=a(nn,2,1)*a(nn,3,2)-a(nn,3,1)
2696 . *a(nn,2,2)
2697 b(nn,1,2)=a(nn,2,2)*a(nn,3,3)-a(nn,3,2)
2698 . *a(nn,2,3)
2699 b(nn,1,3)=a(nn,2,3)*a(nn,3,1)-a(nn,3,3)
2700 . *a(nn,2,1)
2701 b(nn,2,1)=a(nn,3,1)*a(nn,1,2)-a(nn,1,1)
2702 . *a(nn,3,2)
2703 b(nn,2,2)=a(nn,3,2)*a(nn,1,3)-a(nn,1,2)
2704 . *a(nn,3,3)
2705 b(nn,2,3)=a(nn,3,3)*a(nn,1,1)-a(nn,1,3)
2706 . *a(nn,3,1)
2707 b(nn,3,1)=a(nn,1,1)*a(nn,2,2)-a(nn,2,1)
2708 . *a(nn,1,2)
2709 b(nn,3,2)=a(nn,1,2)*a(nn,2,3)-a(nn,2,2)
2710 . *a(nn,1,3)
2711 b(nn,3,3)=a(nn,1,3)*a(nn,2,1)-a(nn,2,3)
2712 . *a(nn,1,1)
2713 xmag1(nn)=b(nn,1,1)*b(nn,1,1)+b(nn,2,1)*b(nn,2,1)+
2714 . b(nn,3,1)*b(nn,3,1)
2715 xmag2(nn)=b(nn,1,2)*b(nn,1,2)+b(nn,2,2)*b(nn,2,2)+
2716 . b(nn,3,2)*b(nn,3,2)
2717 xmag3(nn)=b(nn,1,3)*b(nn,1,3)+b(nn,2,3)*b(nn,2,3)+
2718 . b(nn,3,3)*b(nn,3,3)
2719
2720 ENDDO
2721C
2722 nindex1 = 0
2723 nindex2 = 0
2724 DO nn = 1, nel
2725 xmaxv(nn)=max(xmag1(nn),xmag2(nn),xmag3(nn))
2726 IF(xmag1(nn)==xmaxv(nn)) THEN
2727 lmaxv(nn) = 1
2728 ELSEIF(xmag2(nn)==xmaxv(nn)) THEN
2729 lmaxv(nn) = 2
2730 ELSE
2731 lmaxv(nn) = 3
2732 ENDIF
2733 IF(aaa(nn)<norminf(nn)) THEN
2734 valdp1(nn) = sig(nn,1)
2735 valdp2(nn) = sig(nn,2)
2736 valdp3(nn) = sig(nn,3)
2737 v(nn,1,1) = one
2738 v(nn,2,1) = zero
2739 v(nn,3,1) = zero
2740 v(nn,1,2) = zero
2741 v(nn,2,2) = one
2742 v(nn,3,2) = zero
2743 ELSEIF(xmaxv(nn)>tol1(nn)*tol1(nn)) THEN
2744 nindex1 = nindex1 + 1
2745 index1(nindex1) = nn
2746 ELSE
2747 nindex2 = nindex2 + 1
2748 index2(nindex2) = nn
2749 ENDIF
2750 ENDDO
2751C
2752#include "vectorize.inc"
2753 DO n = 1, nindex1
2754 nn = index1(n)
2755 lmax = lmaxv(nn)
2756 xmaxinv = one/sqrt(xmaxv(nn))
2757 v(nn,1,1)=b(nn,1,lmax)*xmaxinv
2758 v(nn,2,1)=b(nn,2,lmax)*xmaxinv
2759 v(nn,3,1)=b(nn,3,lmax)*xmaxinv
2760 a(nn,1,1)=a(nn,1,1)+str1(nn)-str3(nn)
2761 a(nn,2,2)=a(nn,2,2)+str1(nn)-str3(nn)
2762 a(nn,3,3)=a(nn,3,3)+str1(nn)-str3(nn)
2763C
2764 b(nn,1,1)=a(nn,2,1)*v(nn,3,1)-a(nn,3,1)*v(nn,2,1)
2765 b(nn,1,2)=a(nn,2,2)*v(nn,3,1)-a(nn,3,2)*v(nn,2,1)
2766 b(nn,1,3)=a(nn,2,3)*v(nn,3,1)-a(nn,3,3)*v(nn,2,1)
2767 b(nn,2,1)=a(nn,3,1)*v(nn,1,1)-a(nn,1,1)*v(nn,3,1)
2768 b(nn,2,2)=a(nn,3,2)*v(nn,1,1)-a(nn,1,2)*v(nn,3,1)
2769 b(nn,2,3)=a(nn,3,3)*v(nn,1,1)-a(nn,1,3)*v(nn,3,1)
2770 b(nn,3,1)=a(nn,1,1)*v(nn,2,1)-a(nn,2,1)*v(nn,1,1)
2771 b(nn,3,2)=a(nn,1,2)*v(nn,2,1)-a(nn,2,2)*v(nn,1,1)
2772 b(nn,3,3)=a(nn,1,3)*v(nn,2,1)-a(nn,2,3)*v(nn,1,1)
2773 xmag1(nn)=b(nn,1,1)*b(nn,1,1)+b(nn,2,1)*b(nn,2,1)+
2774 . b(nn,3,1)*b(nn,3,1)
2775 xmag2(nn)=b(nn,1,2)*b(nn,1,2)+b(nn,2,2)*b(nn,2,2)+
2776 . b(nn,3,2)*b(nn,3,2)
2777 xmag3(nn)=b(nn,1,3)*b(nn,1,3)+b(nn,2,3)*b(nn,2,3)+
2778 . b(nn,3,3)*b(nn,3,3)
2779C
2780 xmaxv(nn)=max(xmag1(nn),xmag2(nn),xmag3(nn))
2781 ENDDO
2782C
2783#include "vectorize.inc"
2784 DO n = 1, nindex1
2785 nn = index1(n)
2786 IF(xmag1(nn)==xmaxv(nn)) THEN
2787 lmaxv(nn) = 1
2788 ELSEIF(xmag2(nn)==xmaxv(nn)) THEN
2789 lmaxv(nn) = 2
2790 ELSE
2791 lmaxv(nn) = 3
2792 ENDIF
2793C
2794 vmag=sqrt(v(nn,1,1)**2+v(nn,2,1)**2)
2795 lmax = lmaxv(nn)
2796 xmaxv(nn)=sqrt(xmaxv(nn))
2797 IF(xmaxv(nn)>tol2(nn))THEN
2798 xmaxinv = one/xmaxv(nn)
2799 v(nn,1,3)=b(nn,1,lmax)*xmaxinv
2800 v(nn,2,3)=b(nn,2,lmax)*xmaxinv
2801 v(nn,3,3)=b(nn,3,lmax)*xmaxinv
2802 v(nn,1,2)=v(nn,2,3)*v(nn,3,1)-v(nn,2,1)*v(nn,3,3)
2803 v(nn,2,2)=v(nn,3,3)*v(nn,1,1)-v(nn,3,1)*v(nn,1,3)
2804 v(nn,3,2)=v(nn,1,3)*v(nn,2,1)-v(nn,1,1)*v(nn,2,3)
2805 c1 = v(nn,1,2)*v(nn,1,2)+v(nn,2,2)*v(nn,2,2)+
2806 & v(nn,3,2)*v(nn,3,2)
2807 vmag=one/sqrt(c1)
2808 v(nn,1,2)=v(nn,1,2)*vmag
2809 v(nn,2,2)=v(nn,2,2)*vmag
2810 v(nn,3,2)=v(nn,3,2)*vmag
2811 ELSEIF(vmag>tol2(nn))THEN
2812 v(nn,1,2)=-v(nn,2,1)/vmag
2813 v(nn,2,2)=v(nn,1,1)/vmag
2814 v(nn,3,2)=zero
2815 ELSE
2816 v(nn,1,2)=one
2817 v(nn,2,2)=zero
2818 v(nn,3,2)=zero
2819 ENDIF
2820 ENDDO
2821C . . . . . . . . . . . . .
2822C SOLUTION DOUBLE
2823C . . . . . . . . . . . . .
2824#include "vectorize.inc"
2825 DO n = 1, nindex2
2826 nn = index2(n)
2827 xmag1(nn)=a(nn,1,1)*a(nn,1,1)+a(nn,2,1)*a(nn,2,1)
2828 xmag2(nn)=a(nn,1,2)*a(nn,1,2)+a(nn,2,2)*a(nn,2,2)
2829 xmag3(nn)=a(nn,1,3)*a(nn,1,3)+a(nn,2,3)*a(nn,2,3)
2830C
2831 xmaxv(nn)=max(xmag1(nn),xmag2(nn),xmag3(nn))
2832 ENDDO
2833C
2834#include "vectorize.inc"
2835 DO n = 1, nindex2
2836 nn = index2(n)
2837 IF(xmag1(nn)==xmaxv(nn)) THEN
2838 lmaxv(nn) = 1
2839 ELSEIF(xmag2(nn)==xmaxv(nn)) THEN
2840 lmaxv(nn) = 2
2841 ELSE
2842 lmaxv(nn) = 3
2843 ENDIF
2844C
2845 lmax = lmaxv(nn)
2846 xmaxv(nn)=sqrt(xmaxv(nn))
2847 IF(max(abs(a(nn,3,1)),abs(a(nn,3,2)),abs(a(nn,3,3)))
2848 & <tol2(nn))THEN
2849 xmaxinv = one/xmaxv(nn)
2850 v(nn,1,1)= zero
2851 v(nn,2,1)= zero
2852 v(nn,3,1)= one
2853 v(nn,1,2)=-a(nn,2,lmax)*xmaxinv
2854 v(nn,2,2)= a(nn,1,lmax)*xmaxinv
2855 v(nn,3,2)= zero
2856C
2857 ELSEIF(xmaxv(nn)>tol2(nn))THEN
2858 xmaxinv = one/xmaxv(nn)
2859 v(nn,1,1)=-a(nn,2,lmax)*xmaxinv
2860 v(nn,2,1)= a(nn,1,lmax)*xmaxinv
2861 v(nn,3,1)= zero
2862 v(nn,1,2)=-a(nn,3,lmax)*v(nn,2,1)
2863 v(nn,2,2)= a(nn,3,lmax)*v(nn,1,1)
2864 v(nn,3,2)= a(nn,1,lmax)*v(nn,2,1)-a(nn,2,lmax)*v(nn,1,1)
2865 c1 = v(nn,1,2)*v(nn,1,2)+v(nn,2,2)*v(nn,2,2)+
2866 & v(nn,3,2)*v(nn,3,2)
2867 vmag=one/sqrt(c1)
2868 v(nn,1,2)=v(nn,1,2)*vmag
2869 v(nn,2,2)=v(nn,2,2)*vmag
2870 v(nn,3,2)=v(nn,3,2)*vmag
2871 ELSE
2872 valdp1(nn) = sig(nn,1)
2873 valdp2(nn) = sig(nn,2)
2874 valdp3(nn) = sig(nn,3)
2875 v(nn,1,1) = one
2876 v(nn,2,1) = zero
2877 v(nn,3,1) = zero
2878 v(nn,1,2) = zero
2879 v(nn,2,2) = one
2880 v(nn,3,2) = zero
2881 index3(nn) = 0
2882 ENDIF
2883 ENDDO
2884C
2885 DO nn = 1, nel
2886 vecdp1(nn)=v(nn,1,1)
2887 vecdp2(nn)=v(nn,2,1)
2888 vecdp3(nn)=v(nn,3,1)
2889 vecdp4(nn)=v(nn,1,2)
2890 vecdp5(nn)=v(nn,2,2)
2891 vecdp6(nn)=v(nn,3,2)
2892 vecdp7(nn)=vecdp2(nn)*vecdp6(nn)-vecdp3(nn)*vecdp5(nn)
2893 vecdp8(nn)=vecdp3(nn)*vecdp4(nn)-vecdp1(nn)*vecdp6(nn)
2894 vecdp9(nn)=vecdp1(nn)*vecdp5(nn)-vecdp2(nn)*vecdp4(nn)
2895 ENDDO
2896C
2897 DO nn = 1, nel
2898 val(nn,1)=valdp1(nn)
2899 val(nn,2)=valdp2(nn)
2900 val(nn,3)=valdp3(nn)
2901 vec(nn,1)=vecdp1(nn)
2902 vec(nn,2)=vecdp2(nn)
2903 vec(nn,3)=vecdp3(nn)
2904 vec(nn,4)=vecdp4(nn)
2905 vec(nn,5)=vecdp5(nn)
2906 vec(nn,6)=vecdp6(nn)
2907 vec(nn,7)=vecdp7(nn)
2908 vec(nn,8)=vecdp8(nn)
2909 vec(nn,9)=vecdp9(nn)
2910 ENDDO
2911C reecriture pour contourner probleme sur itanium2 comp 9. + latency=16
2912 IF (nindex3 > 0) THEN
2913 DO nn = 1, nel
2914 IF(index3(nn) == 1) THEN
2915C strutilise com tableau temporaire au lieu de scalaires temporaires
2916 s1=vec(nn,7)
2917 s2=vec(nn,8)
2918 s3=vec(nn,9)
2919 vec(nn,7)=vec(nn,1)
2920 vec(nn,8)=vec(nn,2)
2921 vec(nn,9)=vec(nn,3)
2922 vec(nn,1)=-s1
2923 vec(nn,2)=-s2
2924 vec(nn,3)=-s3
2925 ENDIF
2926 ENDDO
2927 END IF !(NINDEX3 > 0) THEN
2928 RETURN