C CUBTWE SOURCE CHAT 05/01/12 22:32:47 5004 SUBROUTINE CUBTWE(ICM,XE,XEC,XX,SHPC,SGA,QSI,IG) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(3,20),A1(3,8),A2(2,8),Q(3,20),P(3,8), & AJ(3,3),X(3),XE(3,*),XEC(3,2),XX(3,2),AL(3),DS(6),XT(3), & BJ(3,3),F(3),SHPC(6,20,*),SGA(3,5,5),QSI(3,2),DT(3),SIB(3) DATA EPSI/1.D-7/, UN/1.D0/,TWO/2.D0/ DATA DS/-1.D0, 1.D0, -1.D0, 1.D0, -1.D0,-1.D0/ C --- ORDRE: 1 3 5 7 13 15 17 19; 2 6 14 18; 4 8 16 20; 9 10 11 12 DATA Q/-1.D0,-1.D0,-1.D0, 1.D0,-1.D0,-1.D0, 1.D0, 1.D0,-1.D0, & -1.D0, 1.D0,-1.D0, -1.D0,-1.D0, 1.D0, 1.D0,-1.D0, 1.D0, & 1.D0, 1.D0, 1.D0, -1.D0, 1.D0, 1.D0, 0.D0,-1.D0,-1.D0, & 0.D0, 1.D0,-1.D0, 0.D0,-1.D0, 1.D0, 0.D0, 1.D0, 1.D0, & 1.D0, 0.D0,-1.D0, -1.D0, 0.D0,-1.D0, 1.D0, 0.D0, 1.D0, & -1.D0, 0.D0, 1.D0, -1.D0,-1.D0, 0.D0, 1.D0,-1.D0, 0.D0, & 1.D0, 1.D0, 0.D0, -1.D0, 1.D0, 0.D0/ D12=0.D0 DO 5 I=1,3 D12=D12+(XEC(I,2)-XEC(I,1))**2 DO 5 J=1,20 5 A(I,J)=0.D0 D12=DSQRT(D12) DO 6 I=1,4 K=2*I-1 DO 6 J=1,3 P(J,I)=XE(J,K)*0.125D0 6 P(J,I+4)=XE(J,K+12)*0.125D0 C IG: NUMERO DU POINT DE GAUSS DO 10 J=1,8 DO 10 I=1,3 CO=P(I,J) DO 14 JH=1,3 14 X(JH)=CO*Q(JH,J) A(I,1)=A(I,1)-2.D0*CO A(I,2)=A(I,2)-X(1) A(I,3)=A(I,3)-X(2) A(I,4)=A(I,4)-X(3) A(I,5)=A(I,5)+CO A(I,7)=A(I,7)+CO A(I,9)=A(I,9)+CO A(I,11)=A(I,11)+X(2) A(I,12)=A(I,12)+X(1) A(I,13)=A(I,13)+X(3) A(I,14)=A(I,14)+X(2) A(I,15)=A(I,15)+X(1) A(I,16)=A(I,16)+X(3) A(I,17)=A(I,17)+X(1)*Q(2,J)*Q(3,J) A(I,18)=A(I,18)+X(2)*Q(3,J) A(I,19)=A(I,19)+X(1)*Q(3,J) 10 A(I,20)=A(I,20)+X(1)*Q(2,J) DO 11 I=1,3 P(I,1)=XE(I,2)*0.25D0 P(I,2)=XE(I,6)*0.25D0 P(I,3)=XE(I,14)*0.25D0 11 P(I,4)=XE(I,18)*0.25D0 DO 12 J=1,4 K=J+8 DO 12 I=1,3 CO=P(I,J) CC=CO*Q(2,K) A(I,1)=A(I,1)+CO A(I,3)=A(I,3)+CC A(I,4)=A(I,4)+CO*Q(3,K) A(I,5)=A(I,5)-CO A(I, 8)=A(I, 8)+CC*Q(3,K) A(I,11)=A(I,11)-CC A(I,16)=A(I,16)-CO*Q(3,K) 12 A(I,18)=A(I,18)-CC*Q(3,K) DO 21 I=1,3 P(I,1)=XE(I,4)*0.25D0 P(I,2)=XE(I,8)*0.25D0 P(I,3)=XE(I,16)*0.25D0 21 P(I,4)=XE(I,20)*0.25D0 DO 22 J=1,4 K=J+12 DO 22 I=1,3 CO=P(I,J) CC=CO*Q(1,K) A(I,1)=A(I,1)+CO A(I,2)=A(I,2)+CC A(I,4)=A(I,4)+CO*Q(3,K) A(I,7)=A(I,7)-CO A(I,10)=A(I,10)+CC*Q(3,K) A(I,12)=A(I,12)-CC A(I,13)=A(I,13)-CO*Q(3,K) 22 A(I,19)=A(I,19)-CC*Q(3,K) DO 31 I=1,3 DO 31 J=1,4 31 P(I,J)=XE(I,J+8)*0.25D0 DO 32 J=1,4 K=J+16 DO 32 I=1,3 CO=P(I,J) CC=CO*Q(1,K) A(I,1)=A(I,1)+CO A(I,2)=A(I,2)+CC A(I,3)=A(I,3)+CO*Q(2,K) A(I,6)=A(I,6)+CC*Q(2,K) A(I,9)=A(I,9)-CO A(I,14)=A(I,14)-CO*Q(2,K) A(I,15)=A(I,15)-CC 32 A(I,20)=A(I,20)-CC*Q(2,K) C DO 15 I=1,3 A(I,1)=A(I,1)-XX(I,IG) 15 X(I)=0.D0 C ITER=0 ICRI=0 100 CONTINUE C --- LES DERIVEES DES FONCTIONS DE TRANSFORMATION 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=X(1)*X(2)*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)+ & TWO*A(I,11)*X12+A(I,12)*X22+A(I,15)*X33+TWO*A(I,16)*X31+ & A(I,17)*X23+TWO*A(I,18)*X123+A(I,19)*X22*X(3)+A(I,20)*X(2)*X33 C AJ(I,2)=A(I,3)+A(I,6)*X(1)+TWO*A(I,7)*X(2)+A(I,8)*X(3)+ & A(I,11)*X11+TWO*A(I,12)*X12+TWO*A(I,13)*X23+A(I,14)*X33+ & A(I,17)*X31+A(I,18)*X11*X(3)+TWO*A(I,19)*X123+A(I,20)*X(1)*X33 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,13)*X22+TWO*A(I,14)*X23+TWO*A(I,15)*X31+A(I,16)*X11+ & A(I,17)*X12+A(I,18)*X11*X(2)+A(I,19)*X(1)*X22+TWO*A(I,20)*X123 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)*X11*X(2)+A(I,12)*X(1)*X22+A(I,13)*X22*X(3) & +A(I,14)*X(2)*X33+A(I,15)*X(1)*X33+A(I,16)*X11*X(3) & +(A(I,17)+A(I,18)*X(1)+A(I,19)*X(2)+A(I,20)*X(3))*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 DO 60 I=1,3 60 QSI(I,IG)=X(I) C --- Calcul des fonctions d'interpolation leur derivees C --- Fonctions d'interpolation Ni DO 310 I=1,4 I1=2*I-1 I2=I1+12 S1=1.D0+X(1)*Q(1,I) S2=1.D0+X(2)*Q(2,I) S3=1.D0+X(3)*Q(3,I) S4=-1.D0+X(1)*Q(1,I)+X(2)*Q(2,I)+X(3)*Q(3,I) J=I+4 T1=1.D0+X(1)*Q(1,J) T2=1.D0+X(2)*Q(2,J) T3=1.D0+X(3)*Q(3,J) T4=-1.D0+X(1)*Q(1,J)+X(2)*Q(2,J)+X(3)*Q(3,J) C --- NOEUDS 1 3 5 7 13 15 17 19 SHPC(1,I1,IG)=0.125D0*S1*S2*S3*(-UN+S4) SHPC(1,I2,IG)=0.125D0*T1*T2*T3*(-UN+T4) IF(ICM.NE.0) THEN SHPC(2,I1,IG)=0.125D0*Q(1,I)*S2*S3*(S4+X(1)*Q(1,I)) SHPC(2,I2,IG)=0.125D0*Q(1,J)*T2*T3*(T4+X(1)*Q(1,J)) SHPC(3,I1,IG)=0.125D0*Q(2,I)*S1*S3*(S4+X(2)*Q(2,I)) SHPC(3,I2,IG)=0.125D0*Q(2,J)*T1*T3*(T4+X(2)*Q(2,J)) SHPC(4,I1,IG)=0.125D0*Q(3,I)*S1*S2*(S4+X(3)*Q(3,I)) SHPC(4,I2,IG)=0.125D0*Q(3,J)*T1*T2*(T4+X(3)*Q(3,J)) END IF C C --- NOEUDS 2 6 14 18 S0=1.D0-X(1)*X(1) I3=(I-1)*4+2 IF(I.GE.3) I3=I*4+2 J=I+8 S1=1.D0+X(2)*Q(2,J) S2=1.D0+X(3)*Q(3,J) SHPC(1,I3,IG)=0.25D0*S0*S1*S2 IF(ICM.NE.0) THEN SHPC(2,I3,IG)=-0.5D0*X(1)*S1*S2 SHPC(3,I3,IG)=0.25D0*Q(2,J)*S0*S2 SHPC(4,I3,IG)=0.25D0*Q(3,J)*S0*S1 END IF C C --- NOEUDS 4 8 16 20 I4=I*4 IF(I.GE.3) I4=I*4+4 J=I+12 S0=1.D0-X(2)*X(2) S1=1.D0+X(1)*Q(1,J) S2=1.D0+X(3)*Q(3,J) SHPC(1,I4,IG)=0.25D0*S0*S1*S2 IF(ICM.NE.0) THEN SHPC(2,I4,IG)=0.25D0*Q(1,J)*S0*S2 SHPC(3,I4,IG)=-0.5D0*X(2)*S1*S2 SHPC(4,I4,IG)=0.25D0*Q(3,J)*S0*S1 END IF C C --- NOEUDS 9 10 11 12 I5=I+8 J=I+16 S0=1.D0-X(3)*X(3) S1=1.D0+X(1)*Q(1,J) S2=1.D0+X(2)*Q(2,J) SHPC(1,I5,IG)=0.25D0*S0*S1*S2 IF(ICM.NE.0) THEN SHPC(2,I5,IG)=0.25D0*Q(1,J)*S0*S2 SHPC(3,I5,IG)=0.25D0*Q(2,J)*S0*S1 SHPC(4,I5,IG)=-0.5D0*X(3)*S1*S2 END IF 310 CONTINUE CALL PROSGA(ICM,XE,XEC,XX,X,SHPC,SGA,QSI,IG,20) END