1 : C QUADRI SOURCE CHAT 05/01/13 02:40:46 5004 2 : SUBROUTINE QUADRI(ICM,IDIM,XE,XEC,XX,SHPC,SGA,QSI,IG) 3 : IMPLICIT INTEGER(I-N) 4 : IMPLICIT REAL*8 (A-H,O-Z) 5 : DIMENSION A(2,8),AJ(2,2),X(2),XE(3,*),XEC(3,2),XX(3,2), 6 : & BJ(2,2),F(2),SHPC(6,20,*),SGA(3,5,5),QSI(3,2),DT(2),SIB(2) 7 : DATA EPSI/1.D-7/ 8 : C --- IG: NUMERO DU POINT DE GAUSS 9 : DO 10 I=1,2 10 : A(I,1)=-XE(I,1)-XE(I,3)-XE(I,5)-XE(I,7) 11 : & +2.D0*(XE(I,2)+XE(I,4)+XE(I,6)+XE(I,8))-4.D0*XX(I,IG) 12 : A(I,2)= 2.D0*( XE(I,4)-XE(I,8)) 13 : A(I,3)= 2.D0*(-XE(I,2)+XE(I,6)) 14 : A(I,4)= XE(I,1)-2.D0*XE(I,2)+XE(I,3)+XE(I,5)-2.D0*XE(I,6)+XE(I,7) 15 : A(I,5)= XE(I,1)-XE(I,3)+XE(I,5)-XE(I,7) 16 : A(I,6)= XE(I,1)+XE(I,3)-2.D0*XE(I,4)+XE(I,5)+XE(I,7)-2.D0*XE(I,8) 17 : A(I,7)=-XE(I,1)+2.D0*XE(I,2)-XE(I,3)+XE(I,5)-2.D0*XE(I,6)+XE(I,7) 18 : 10 A(I,8)=-XE(I,1)+XE(I,3)-2.D0*XE(I,4)+XE(I,5)-XE(I,7)+2.D0*XE(I,8) 19 : C --- Initialisation 20 : DO 15 I=1,2 21 : 15 X(I)=0.D0 22 : C 23 : ITER=0 24 : ICRI=0 25 : 100 CONTINUE 26 : C --- LES DERIVEES DES FONCTIONS DE TRANSFORMATION 27 : QP=X(1)*X(2) 28 : Q1=X(1)*X(1) 29 : Q2=X(2)*X(2) 30 : DO 20 I=1,IDIM 31 : AJ(I,1)=A(I,2)+2.D0*A(I,4)*X(1)+A(I,5)*X(2) 32 : & +2.D0*A(I,7)*QP+A(I,8)*Q2 33 : 20 AJ(I,2)=A(I,3)+A(I,5)*X(1)+2.D0*A(I,6)*X(2) 34 : & +A(I,7)*Q1+2.D0*A(I,8)*QP 35 : C --- Calcul de l'inverse de la matrice AJ et son determinant :BJ=inverse(AJ) 35 : CALL INVER2(AJ,BJ,DET,ICRI) 36 : DO 40 I=1,IDIM 37 : 40 F(I)=A(I,1)+A(I,2)*X(1)+A(I,3)*X(2)+A(I,4)*Q1+A(I,5)*QP 38 : & +A(I,6)*Q2+A(I,7)*Q1*X(2)+A(I,8)*Q2*X(1) 39 : DD=SQRT(F(1)*F(1)+F(2)*F(2)) 40 : IF(DD.GT.EPSI) THEN 41 : ITER=ITER+1 42 : DO 50 I=1,IDIM 43 : 50 X(I)=X(I)-BJ(I,1)*F(1)-BJ(I,2)*F(2) 44 : GOTO 100 45 : END IF 46 : DO 60 I=1,IDIM 47 : 60 QSI(I,IG)=X(I) 48 : C=0.25D0 49 : D=0.5D0 50 : YK1=1.D0+X(1) 51 : YK2=1.D0-X(1) 52 : YN1=1.D0+X(2) 53 : YN2=1.D0-X(2) 54 : YKN=X(1)+X(2) 55 : YNK=X(1)-X(2) 56 : C --- Calcul des fonctions d'interpolation leur derivees 57 : C --- Fonctions d'interpolation Ni 58 : SHPC(1,1,IG)=-C*YK2*YN2*(YK1+X(2)) 59 : SHPC(1,2,IG)= D*YK1*YK2*YN2 60 : SHPC(1,3,IG)=-C*YK1*YN2*(YK2+X(2)) 61 : SHPC(1,4,IG)= D*YK1*YN1*YN2 62 : SHPC(1,5,IG)=-C*YK1*YN1*(YK2-X(2)) 63 : SHPC(1,6,IG)= D*YK1*YK2*YN1 64 : SHPC(1,7,IG)=-C*YK2*YN1*(YK1-X(2)) 65 : SHPC(1,8,IG)= D*YK2*YN1*YN2 66 : IF(ICM.EQ.0) GOTO 1000 67 : C --- Derivee de Ni par rapport a QSI 68 : SHPC(2,1,IG)= C*YN2*(X(1)+YKN) 69 : SHPC(2,2,IG)=-YN2*X(1) 70 : SHPC(2,3,IG)= C*YN2*(X(1)+YNK) 71 : SHPC(2,4,IG)= D*YN1*YN2 72 : SHPC(2,5,IG)= C*YN1*(X(1)+YKN) 73 : SHPC(2,6,IG)=-YN1*X(1) 74 : SHPC(2,7,IG)= C*YN1*(X(1)+YNK) 75 : SHPC(2,8,IG)=-SHPC(2,4,IG) 76 : C --- Derivee de Ni par rapport a ETA 77 : SHPC(3,1,IG)= C*YK2*(YKN+X(2)) 78 : SHPC(3,2,IG)=-D*YK1*YK2 79 : SHPC(3,3,IG)=-C*YK1*(YNK-X(2)) 80 : SHPC(3,4,IG)=-YK1*X(2) 81 : SHPC(3,5,IG)= C*YK1*(YKN+X(2)) 82 : SHPC(3,6,IG)=-SHPC(3,2,IG) 83 : SHPC(3,7,IG)=-C*YK2*(YNK-X(2)) 84 : SHPC(3,8,IG)=-YK2*X(2) 85 : C 86 : DJ=0.D0 87 : DO 70 I=1,IDIM 88 : DT(I)=0.50D0*(XEC(I,2)-XEC(I,1)) 89 : SIB(I)=DT(I) 90 : DJ=DJ+DT(I)*DT(I) 91 : DO 70 J=1,IDIM 92 : 70 BJ(I,J)=0.D0 93 : DO 90 I=1,IDIM 94 : DO 90 J=1,IDIM 95 : DO 80 K=1,8 96 : 80 BJ(I,J)=BJ(I,J)+XE(I,K)*SHPC(J+1,K,IG) 97 : 90 CONTINUE 98 : 99 : C --- LES DRIVEES 99 : CALL LINEQ1(BJ,DT,X) 100 : DJ=SQRT(DJ) 101 : DO 120 I=1,IDIM 102 : 120 SIB(I)=SIB(I)/DJ 103 : 1000 CONTINUE 104 : DO 2000 IS=1,IDIM 105 : C --- Coordonnees du point d'integration dans le systeme reel (x,y,z) 106 : SGA(IS,1,IG)=XX(IS,IG) 107 : C --- Coordonnees du point d'integration dans le systeme de reference 108 : C (qsi,eta,dze) 109 : SGA(IS,2,IG)=QSI(IS,IG) 110 : C --- Derivee de (qsi,eta) par rapport a miu 111 : SGA(IS,3,IG)=X(IS) 112 : C --- Cosinus du cable avec les axes 113 : SGA(IS,4,IG)=SIB(IS) 114 : C --- Jacobien 115 : SGA(IS,5,IG)=DJ 116 : 2000 CONTINUE 117 : END 118 : 119 :
© Cast3M 2003 - All rights reserved.
Disclaimer