C SHPCAB SOURCE CHAT 05/01/13 03:15:42 5004 SUBROUTINE SHPCAB(IGAU,XEL,XEC,MELE,NBNO,IDIM, & QSIG,SHPC,XX,SGA,ICM) C======================================================================C C C C INPUT: C C C C MELE: numero de l'element beton C C QSIG: COORDONNÉES DES POINTS DE GAUSS UNIDIMENSIONNELS C C IGAU: NOMBRE DE POINTS DE GAUSS C C XEC : LES COORDONNEES DU CABLE C C XEL : LES COORDONNEES DE L'ELEMENT BETON C C IGAU: NOMBRE DE POINTS DE GAUSS C C NBNO: NOMBRE DE NOEUDS C C LRE : NOMBRE DE COLONNES DE LA MATRICE B C C QSIG: COORDONNEES DES POINTS DE GAUSS A UNE DIMENSION C C ICM = 0: CALCUL DES FONCTIONS DE FORME SEULEMENT C C ICM =NON 0: CALCUL DES FONCTIONS DE FORME ET LEURS DERIVEES C C C C OUTPUT: C C C C DJ : JACOBIEN C C XX(3,2) : COORDONNEES DES POINTS DE GAUSS DANS LE SYSTEME GLOBAL C C QSI(3,2): COORDONNEES DES POINTS DE GAUSS DANS L'ELEMENT DE C C REFERENCE C C SHPC(6,NBNO,*)=FONCTIONS DE FORMES ET DERIVEES (CABLE) C C SGAU: SIN, COS C C C C G. L. JIANG C C======================================================================C IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) PARAMETER(XZER=0.D0) DIMENSION XX(3,2),QSI(3,2),XEC(3,2),XEL(3,*),SIB(3), & SHPC(6,20,5),SGA(3,5,5),QSIG(6),COE(3,3),BOE(3), & XA(6),YA(6),ZA(6) DATA EPSI/1.D-6/ DO 10 I1=1,6 DO 10 I2=1,NBNO DO 10 I3=1,IGAU 10 SHPC(I1,I2,I3)=0.D0 IF(ICM.NE.0) THEN DO 11 I1=1,3 DO 11 I2=1,5 DO 11 I3=1,5 11 SGA(I1,I2,I3)=0.0 DO 15 IG=1,IGAU X1=0.5D0*(1.D0-QSIG(IG)) X2=0.5D0*(1.D0+QSIG(IG)) DO 12 IS=1,IDIM 12 XX(IS,IG)=X1*XEC(IS,1)+X2*XEC(IS,2) 15 CONTINUE END IF C C --- TRI3 C IF(MELE.EQ.4) THEN A11=XEL(1,2)-XEL(1,1) A12=XEL(1,3)-XEL(1,1) A21=XEL(2,2)-XEL(2,1) A22=XEL(2,3)-XEL(2,1) DET=A11*A22-A12*A21 DO 20 IG=1,IGAU B1=XX(1,IG)-XEL(1,1) B2=XX(2,IG)-XEL(2,1) IF(ABS(DET).LT.EPSI) THEN WRITE(6,*) 'SUBROUTINE SHPCAB DET=0' ELSE QSI(1,IG)= (B1*A22-B2*A12)/DET QSI(2,IG)=-(B1*A21-B2*A11)/DET END IF 20 CONTINUE IF(ICM.NE.0) THEN T1=0.5D0*(XEC(1,2)-XEC(1,1)) T2=0.5D0*(XEC(2,2)-XEC(2,1)) DJ=SQRT(T1*T1+T2*T2) C --- COS BETA SIB(1)=T1/DJ C --- SIN BETA SIB(2)=T2/DJ do 21 ig=1,igau sga(1,4,ig)=sib(1) sga(2,4,ig)=sib(2) 21 sga(1,5,ig)=dj END IF C --- Calcul de fonctions de forme DO 30 IG=1,IGAU SHPC(1,1,IG)=1.D0-QSI(1,IG)-QSI(2,IG) SHPC(1,2,IG)=QSI(1,IG) 30 SHPC(1,3,IG)=QSI(2,IG) IF(ICM.EQ.0) GOTO 4000 C --- Derivees par rapport a QSI DO 40 IG=1,IGAU SHPC(2,1,IG)=-1.D0 SHPC(2,2,IG)= 1.D0 40 SHPC(2,3,IG)= 0.D0 C --- Derivees par rapport a ETA DO 50 IG=1,IGAU SHPC(3,1,IG)=-1.D0 SHPC(3,2,IG)= 0.D0 50 SHPC(3,3,IG)= 1.D0 S1= (T1*A22-T2*A12)/DET S2=-(T1*A21-T2*A11)/DET C --- Derivees de QSI par rapport a MIU DO 60 IG=1,IGAU 60 SGA(1,3,IG)=S1 C --- Derivees de ETA par rapport a MIU DO 70 IG=1,IGAU 70 SGA(2,3,IG)=S2 GOTO 4000 END IF C C --- TRI 6 C IF(MELE.EQ.6) THEN DO 72 IG=1,IGAU CALL TRINO6(ICM,2,XEL,XEC,XX,SHPC,SGA,QSI,IG) 72 CONTINUE GOTO 5000 END IF C C --- QUA4 C IF(MELE.EQ.8) THEN S1=XEL(1,1)-XEL(1,2) S2=XEL(1,3)-XEL(1,4) A11= S1+S2 A12=-S1+S2 A13=-XEL(1,1)-XEL(1,2)+XEL(1,3)+XEL(1,4) A21= XEL(2,1)-XEL(2,2)+XEL(2,3)-XEL(2,4) A22=-XEL(2,1)+XEL(2,2)+XEL(2,3)-XEL(2,4) A23=-XEL(2,1)-XEL(2,2)+XEL(2,3)+XEL(2,4) DO 320 IG=1,IGAU A14= XEL(1,1)+XEL(1,2)+XEL(1,3)+XEL(1,4)-4.D0*XX(1,IG) A24= XEL(2,1)+XEL(2,2)+XEL(2,3)+XEL(2,4)-4.D0*XX(2,IG) IF(ABS(A11).LT.EPSI.AND.ABS(A21).LT.EPSI) THEN CALL LINEQU(A12,A13,A22,A23,-A14,-A24,QSI(1,IG),QSI(2,IG)) GOTO 320 END IF D1=A12*A21-A22*A11 D2=A13*A21-A23*A11 D3=A14*A21-A24*A11 IF(ABS(D2).LT.EPSI) THEN QSI(1,IG)=-D3/D1 S1=A11*QSI(1,IG)+A13 IF(ABS(S1).GT.EPSI) THEN QSI(2,IG)=-A14/S1 GOTO 320 ELSE QSI(2,IG)=-A24/(A21*QSI(1,IG)+A23) GOTO 320 END IF END IF A1=-(A12*A21-A22*A11)/D2 A2=-(A14*A21-A24*A11)/D2 A=A11*A1 B=A11*A2+A12+A13*A1 C=A13*A2+A14 IF(ABS(A).GT.EPSI) THEN cs=B*B-4.0*A*C BC=SQRT(B*B-4.0*A*C) QS1=0.5*(-B+BC)/A QS2=0.5*(-B-BC)/A ET1=A1*QS1+A2 ET2=A1*QS2+A2 IF(ABS(QS1).LT.1.0.AND.ABS(ET1).LT.1.0) THEN QSI(1,IG)=QS1 QSI(2,IG)=ET1 ELSE QSI(1,IG)=QS2 QSI(2,IG)=ET2 END IF GOTO 320 ELSE QSI(1,IG)=-C/B QSI(2,IG)=A1*QSI(1,IG)+A2 GOTO 320 END IF 320 CONTINUE IF(ICM.NE.0) THEN T1=0.5D0*(XEC(1,2)-XEC(1,1)) T2=0.5D0*(XEC(2,2)-XEC(2,1)) DJ=SQRT(T1*T1+T2*T2) C --- COS BETA SIB(1)=T1/DJ C --- SIN BETA SIB(2)=T2/DJ do 321 ig=1,igau sga(1,4,ig)=sib(1) sga(2,4,ig)=sib(2) 321 sga(1,5,ig)=dj END IF C --- Calcul de fonctions de forme DO 330 IG=1,IGAU SHPC(1,1,IG)=0.25D0*(1.D0-QSI(1,IG))*(1.D0-QSI(2,IG)) SHPC(1,2,IG)=0.25D0*(1.D0+QSI(1,IG))*(1.D0-QSI(2,IG)) SHPC(1,3,IG)=0.25D0*(1.D0+QSI(1,IG))*(1.D0+QSI(2,IG)) 330 SHPC(1,4,IG)=0.25D0*(1.D0-QSI(1,IG))*(1.D0+QSI(2,IG)) IF(ICM.EQ.0) GOTO 4000 DO 350 IG=1,IGAU C --- Derivees par rapport a QSI SHPC(2,1,IG)=0.25D0*(-1.D0+QSI(2,IG)) SHPC(2,2,IG)=0.25D0*( 1.D0-QSI(2,IG)) SHPC(2,3,IG)=0.25D0*( 1.D0+QSI(2,IG)) SHPC(2,4,IG)=0.25D0*(-1.D0-QSI(2,IG)) C --- Derivees par rapport a ETA SHPC(3,1,IG)=0.25D0*(-1.D0+QSI(1,IG)) SHPC(3,2,IG)=0.25D0*(-1.D0-QSI(1,IG)) SHPC(3,3,IG)=0.25D0*( 1.D0+QSI(1,IG)) 350 SHPC(3,4,IG)=0.25D0*( 1.D0-QSI(1,IG)) DO 360 IG=1,IGAU A11=0.0 A12=0.0 A21=0.0 A22=0.0 DO 355 IS=1,4 A11=A11+XEL(1,IS)*SHPC(2,IS,IG) A12=A12+XEL(1,IS)*SHPC(3,IS,IG) A21=A21+XEL(2,IS)*SHPC(2,IS,IG) 355 A22=A22+XEL(2,IS)*SHPC(3,IS,IG) DET=A11*A22-A12*A21 S1= (T1*A22-T2*A12)/DET S2=-(T1*A21-T2*A11)/DET C --- Derivees de QSI par rapport a MIU SGA(1,3,IG)=S1 C --- Derivees de ETA par rapport a MIU SGA(2,3,IG)=S2 360 CONTINUE GOTO 4000 END IF C C --- QUA8 C IF(MELE.EQ.10) THEN DO 374 IG=1,IGAU CALL QUADRI(ICM,2,XEL,XEC,XX,SHPC,SGA,QSI,IG) 374 CONTINUE GOTO 5000 END IF C C --- CUB8 C IF(MELE.EQ.14) THEN DO 381 IG=1,IGAU CALL CUBEIT(ICM,XEL,XEC,XX,SHPC,SGA,QSI,IG) 381 CONTINUE GOTO 5000 END IF C C --- CU20 C IF(MELE.EQ.15) THEN DO 391 IG=1,IGAU CALL CUBTWE(ICM,XEL,XEC,XX,SHPC,SGA,QSI,IG) 391 CONTINUE GOTO 5000 END IF C C --- PRI6 C IF(MELE.EQ.16) THEN DO 401 IG=1,IGAU CALL PRISIX(ICM,XEL,XEC,XX,SHPC,SGA,QSI,IG) 401 CONTINUE GOTO 5000 END IF C C --- PR15 C IF(MELE.EQ.17) THEN DO 411 IG=1,IGAU CALL PRIS15(ICM,XEL,XEC,XX,SHPC,SGA,QSI,IG) 411 CONTINUE GOTO 5000 END IF C 4000 CONTINUE DO 2000 IS=1,IDIM C --- Coordonnees du point d'integration dans le systeme reel (x,y,z) SGA(IS,1,IG)=XX(IS,IG) C --- Coordonnees du point d'integration dans le systeme de reference C (qsi,eta, ) SGA(IS,2,IG)=QSI(IS,IG) 2000 CONTINUE 5000 CONTINUE END