C PRIS15 SOURCE CHAT 05/01/13 02:27:41 5004 SUBROUTINE PRIS15(ICM,XE,XEC,XX,SHPC,SGA,QSI,IG) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(3,15),AJ(3,3),X(3),XE(3,*),XEC(3,2),XX(3,2), & BJ(3,3),F(3),SHPC(6,20,*),SGA(3,5,5),QSI(3,2),DT(3),SIB(3) DATA EPSI/1.D-7/, ZE/0.D0/,UN/1.D0/,TWO/2.D0/,DE/0.5D0/ DO 10 I=1,3 S11=XE(I,1)+XE(I,10) S22=XE(I,1)-XE(I,10) T11=XE(I,2)+XE(I,11) T22=XE(I,2)-XE(I,11) U11=XE(I,6)+XE(I,15) U22=XE(I,6)-XE(I,15) W11=XE(I,3)+XE(I,12) Z11=XE(I,5)+XE(I,14) A(I,1)= XE(I,7)-XX(I,IG) A(I,2)=-S11+TWO*T11-W11-XE(I,7)+XE(I,8) A(I,3)=-S11-Z11+TWO*U11-XE(I,7)+XE(I,9) A(I,4)=-DE*S22 A(I,5)= S11-TWO*T11+W11 A(I,6)=TWO*(S11-T11+XE(I,4)+XE(I,13)-U11) A(I,7)=S11+Z11-TWO*U11 A(I,8)=1.5D0*S22+DE*(XE(I,5)-XE(I,14))-TWO*U22 A(I,9)=DE*S11-XE(I,7) A(I,10)=1.5D0*S22-TWO*T22+DE*(XE(I,3)-XE(I,12)) A(I,11)=-S22-XE(I,5)+XE(I,14)+TWO*U22 A(I,12)=DE*(-S11+Z11)+XE(I,7)-XE(I,9) A(I,13)=DE*(-S11+W11)+XE(I,7)-XE(I,8) A(I,14)=-S22+TWO*T22-XE(I,3)+XE(I,12) 10 A(I,15)=TWO*(-S22+T22-XE(I,4)+XE(I,13)+U22) C C --- Initialisation C X(1)=0.25 X(2)=0.25 X(3)=0.00 C ITER=0 ICRI=0 100 CONTINUE X11=X(1)*X(1) X22=X(2)*X(2) X33=X(3)*X(3) X12=X(1)*X(2) X23=X(2)*X(3) X31=X(3)*X(1) X123=X12*X(3) DO 20 I=1,3 AJ(I,1)=A(I,2)+TWO*A(I,5)*X(1)+A(I,6)*X(2)+A(I,10)*X(3) & +A(I,13)*X33+TWO*A(I,14)*X31+A(I,15)*X23 C AJ(I,2)=A(I,3)+A(I,6)*X(1)+TWO*A(I,7)*X(2)+A(I,8)*X(3) & +TWO*A(I,11)*X23+A(I,12)*X33+A(I,15)*X31 C 20 AJ(I,3)=A(I,4)+A(I,8)*X(2)+TWO*A(I,9)*X(3)+A(I,10)*X(1) & +A(I,11)*X22+TWO*A(I,12)*X23+TWO*A(I,13)*X31 & +A(I,14)*X11+A(I,15)*X12 C --- Calcul de l'inverse de la matrice AJ et son determinant :BJ=inverse(AJ) CALL INVERT(AJ,BJ,DET,ICRI) DO 40 I=1,3 40 F(I)=A(I,1)+A(I,2)*X(1)+A(I,3)*X(2)+A(I,4)*X(3)+A(I,5)*X11 & +A(I,6)*X12+A(I,7)*X22+A(I,8)*X23+A(I,9)*X33 & +A(I,10)*X31+A(I,11)*X22*X(3)+A(I,12)*X(2)*X33 & +A(I,13)*X(1)*X33+A(I,14)*X11*X(3)+A(I,15)*X123 DD=SQRT(F(1)*F(1)+F(2)*F(2)+F(3)*F(3)) IF(DD.GT.EPSI) THEN ITER=ITER+1 DO 50 I=1,3 50 X(I)=X(I)-BJ(I,1)*F(1)-BJ(I,2)*F(2)-BJ(I,3)*F(3) GOTO 100 END IF C DO 60 I=1,3 60 QSI(I,IG)=X(I) C --- Calcul des fonctions d'interpolation leur derivees AUX=UN-X(1)-X(2) DZEM=(UN-X(3)) DZEP=(UN+X(3)) DZMP=DZEM*DZEP PAUX =X(1)-DE*X(3)-UN PAUX1=X(2)-DE*X(3)-UN PAUX2=X(1)+X(2)+DE*X(3) PAUX3=X(1)+DE*X(3)-UN PAUX4=X(2)+DE*X(3)-UN PAUX5=X(1)+X(2)-DE*X(3) TX1=TWO*X(1) TX2=TWO*X(2) TXM2=TX2*DZEM TXP2=TX2*DZEP C --- FONCTIONS D'INTERPOLATION SHPC(1,1,IG)=-AUX*DZEM*PAUX2 SHPC(1,2,IG)=TX1*AUX*DZEM SHPC(1,3,IG)=X(1)*DZEM*PAUX SHPC(1,4,IG)=TX1*X(2)*DZEM SHPC(1,5,IG)=X(2)*DZEM*PAUX1 SHPC(1,6,IG)=AUX*TXM2 SHPC(1,7,IG)=AUX*DZMP SHPC(1,8,IG)=X(1)*DZMP SHPC(1,9,IG)=X(2)*DZMP SHPC(1,10,IG)=-AUX*DZEP*PAUX5 SHPC(1,11,IG)=TX1*AUX*DZEP SHPC(1,12,IG)=X(1)*DZEP*PAUX3 SHPC(1,13,IG)=TX1*X(2)*DZEP SHPC(1,14,IG)=X(2)*DZEP*PAUX4 SHPC(1,15,IG)=TXP2*AUX IF(ICM.EQ.0) GOTO 3000 SHPC(2,1,IG)=(PAUX2-AUX)*DZEM SHPC(2,2,IG)=TWO*(AUX-X(1))*DZEM SHPC(2,3,IG)=(PAUX+X(1))*DZEM SHPC(3,1,IG)=(PAUX2-AUX)*DZEM SHPC(3,2,IG)=-TX1*DZEM SHPC(3,3,IG)=ZE SHPC(4,1,IG)=AUX*(PAUX2-DE*DZEM) SHPC(4,2,IG)=-TX1*AUX SHPC(4,3,IG)=-X(1)*(DE*DZEM+PAUX) SHPC(2,4,IG)=TXM2 SHPC(2,5,IG)=ZE SHPC(2,6,IG)=-TXM2 SHPC(3,4,IG)=TX1*DZEM SHPC(3,5,IG)=(X(2)+PAUX1)*DZEM SHPC(3,6,IG)=TWO*(AUX-X(2))*DZEM SHPC(4,4,IG)=-TX1*X(2) SHPC(4,5,IG)=-X(2)*(PAUX1+DE*DZEM) SHPC(4,6,IG)=-TX2*AUX SHPC(2,7,IG)=-DZMP SHPC(2,8,IG)=DZMP SHPC(2,9,IG)=ZE SHPC(3,7,IG)=-DZMP SHPC(3,8,IG)=ZE SHPC(3,9,IG)=DZMP SHPC(4,7,IG)=-TWO*X(3)*AUX SHPC(4,8,IG)=-TX1*X(3) SHPC(4,9,IG)=-TX2*X(3) SHPC(2,10,IG)=(PAUX5-AUX)*DZEP SHPC(2,11,IG)=TWO*(AUX-X(1))*DZEP SHPC(2,12,IG)=(PAUX3+X(1))*DZEP SHPC(3,10,IG)=(PAUX5-AUX)*DZEP SHPC(3,11,IG)=-TX1*DZEP SHPC(3,12,IG)=ZE SHPC(4,10,IG)=AUX*(DE*DZEP-PAUX5) SHPC(4,11,IG)=TX1*AUX SHPC(4,12,IG)=X(1)*(PAUX3+DE*DZEP) SHPC(2,13,IG)=TXP2 SHPC(2,14,IG)=ZE SHPC(2,15,IG)=-TXP2 SHPC(3,13,IG)=TX1*DZEP SHPC(3,14,IG)=(X(2)+PAUX4)*DZEP SHPC(3,15,IG)=TWO*(AUX-X(2))*DZEP SHPC(4,13,IG)=TX1*X(2) SHPC(4,14,IG)=X(2)*(PAUX4+DE*DZEP) SHPC(4,15,IG)=TX2*AUX 3000 CONTINUE CALL PROSGA(ICM,XE,XEC,XX,X,SHPC,SGA,QSI,IG,15) END