1 : C SHPCAB SOURCE CHAT 05/01/13 03:15:42 5004 2 : SUBROUTINE SHPCAB(IGAU,XEL,XEC,MELE,NBNO,IDIM, 3 : & QSIG,SHPC,XX,SGA,ICM) 4 : C======================================================================C 5 : C C 6 : C INPUT: C 7 : C C 8 : C MELE: numero de l'element beton C 9 : C QSIG: COORDONNÉES DES POINTS DE GAUSS UNIDIMENSIONNELS C 10 : C IGAU: NOMBRE DE POINTS DE GAUSS C 11 : C XEC : LES COORDONNEES DU CABLE C 12 : C XEL : LES COORDONNEES DE L'ELEMENT BETON C 13 : C IGAU: NOMBRE DE POINTS DE GAUSS C 14 : C NBNO: NOMBRE DE NOEUDS C 15 : C LRE : NOMBRE DE COLONNES DE LA MATRICE B C 16 : C QSIG: COORDONNEES DES POINTS DE GAUSS A UNE DIMENSION C 17 : C ICM = 0: CALCUL DES FONCTIONS DE FORME SEULEMENT C 18 : C ICM =NON 0: CALCUL DES FONCTIONS DE FORME ET LEURS DERIVEES C 19 : C C 20 : C OUTPUT: C 21 : C C 22 : C DJ : JACOBIEN C 23 : C XX(3,2) : COORDONNEES DES POINTS DE GAUSS DANS LE SYSTEME GLOBAL C 24 : C QSI(3,2): COORDONNEES DES POINTS DE GAUSS DANS L'ELEMENT DE C 25 : C REFERENCE C 26 : C SHPC(6,NBNO,*)=FONCTIONS DE FORMES ET DERIVEES (CABLE) C 27 : C SGAU: SIN, COS C 28 : C C 29 : C G. L. JIANG C 30 : C======================================================================C 31 : IMPLICIT INTEGER(I-N) 32 : IMPLICIT REAL*8 (A-H,O-Z) 33 : PARAMETER(XZER=0.D0) 34 : DIMENSION XX(3,2),QSI(3,2),XEC(3,2),XEL(3,*),SIB(3), 35 : & SHPC(6,20,5),SGA(3,5,5),QSIG(6),COE(3,3),BOE(3), 36 : & XA(6),YA(6),ZA(6) 37 : DATA EPSI/1.D-6/ 38 : DO 10 I1=1,6 39 : DO 10 I2=1,NBNO 40 : DO 10 I3=1,IGAU 41 : 10 SHPC(I1,I2,I3)=0.D0 42 : IF(ICM.NE.0) THEN 43 : DO 11 I1=1,3 44 : DO 11 I2=1,5 45 : DO 11 I3=1,5 46 : 11 SGA(I1,I2,I3)=0.0 47 : DO 15 IG=1,IGAU 48 : X1=0.5D0*(1.D0-QSIG(IG)) 49 : X2=0.5D0*(1.D0+QSIG(IG)) 50 : DO 12 IS=1,IDIM 51 : 12 XX(IS,IG)=X1*XEC(IS,1)+X2*XEC(IS,2) 52 : 15 CONTINUE 53 : END IF 54 : C 55 : C --- TRI3 56 : C 57 : IF(MELE.EQ.4) THEN 58 : A11=XEL(1,2)-XEL(1,1) 59 : A12=XEL(1,3)-XEL(1,1) 60 : A21=XEL(2,2)-XEL(2,1) 61 : A22=XEL(2,3)-XEL(2,1) 62 : DET=A11*A22-A12*A21 63 : DO 20 IG=1,IGAU 64 : B1=XX(1,IG)-XEL(1,1) 65 : B2=XX(2,IG)-XEL(2,1) 66 : IF(ABS(DET).LT.EPSI) THEN 67 : WRITE(6,*) 'SUBROUTINE SHPCAB DET=0' 68 : ELSE 69 : QSI(1,IG)= (B1*A22-B2*A12)/DET 70 : QSI(2,IG)=-(B1*A21-B2*A11)/DET 71 : END IF 72 : 20 CONTINUE 73 : IF(ICM.NE.0) THEN 74 : T1=0.5D0*(XEC(1,2)-XEC(1,1)) 75 : T2=0.5D0*(XEC(2,2)-XEC(2,1)) 76 : DJ=SQRT(T1*T1+T2*T2) 77 : C --- COS BETA 78 : SIB(1)=T1/DJ 79 : C --- SIN BETA 80 : SIB(2)=T2/DJ 81 : do 21 ig=1,igau 82 : sga(1,4,ig)=sib(1) 83 : sga(2,4,ig)=sib(2) 84 : 21 sga(1,5,ig)=dj 85 : END IF 86 : C --- Calcul de fonctions de forme 87 : DO 30 IG=1,IGAU 88 : SHPC(1,1,IG)=1.D0-QSI(1,IG)-QSI(2,IG) 89 : SHPC(1,2,IG)=QSI(1,IG) 90 : 30 SHPC(1,3,IG)=QSI(2,IG) 91 : IF(ICM.EQ.0) GOTO 4000 92 : C --- Derivees par rapport a QSI 93 : DO 40 IG=1,IGAU 94 : SHPC(2,1,IG)=-1.D0 95 : SHPC(2,2,IG)= 1.D0 96 : 40 SHPC(2,3,IG)= 0.D0 97 : C --- Derivees par rapport a ETA 98 : DO 50 IG=1,IGAU 99 : SHPC(3,1,IG)=-1.D0 100 : SHPC(3,2,IG)= 0.D0 101 : 50 SHPC(3,3,IG)= 1.D0 102 : S1= (T1*A22-T2*A12)/DET 103 : S2=-(T1*A21-T2*A11)/DET 104 : C --- Derivees de QSI par rapport a MIU 105 : DO 60 IG=1,IGAU 106 : 60 SGA(1,3,IG)=S1 107 : C --- Derivees de ETA par rapport a MIU 108 : DO 70 IG=1,IGAU 109 : 70 SGA(2,3,IG)=S2 110 : GOTO 4000 111 : END IF 112 : C 113 : C --- TRI 6 114 : C 115 : IF(MELE.EQ.6) THEN 116 : DO 72 IG=1,IGAU 116 : CALL TRINO6(ICM,2,XEL,XEC,XX,SHPC,SGA,QSI,IG) 117 : 72 CONTINUE 118 : GOTO 5000 119 : END IF 120 : C 121 : C --- QUA4 122 : C 123 : IF(MELE.EQ.8) THEN 124 : S1=XEL(1,1)-XEL(1,2) 125 : S2=XEL(1,3)-XEL(1,4) 126 : A11= S1+S2 127 : A12=-S1+S2 128 : A13=-XEL(1,1)-XEL(1,2)+XEL(1,3)+XEL(1,4) 129 : A21= XEL(2,1)-XEL(2,2)+XEL(2,3)-XEL(2,4) 130 : A22=-XEL(2,1)+XEL(2,2)+XEL(2,3)-XEL(2,4) 131 : A23=-XEL(2,1)-XEL(2,2)+XEL(2,3)+XEL(2,4) 132 : DO 320 IG=1,IGAU 133 : A14= XEL(1,1)+XEL(1,2)+XEL(1,3)+XEL(1,4)-4.D0*XX(1,IG) 134 : A24= XEL(2,1)+XEL(2,2)+XEL(2,3)+XEL(2,4)-4.D0*XX(2,IG) 135 : IF(ABS(A11).LT.EPSI.AND.ABS(A21).LT.EPSI) THEN 135 : CALL LINEQU(A12,A13,A22,A23,-A14,-A24,QSI(1,IG),QSI(2,IG)) 136 : GOTO 320 137 : END IF 138 : D1=A12*A21-A22*A11 139 : D2=A13*A21-A23*A11 140 : D3=A14*A21-A24*A11 141 : IF(ABS(D2).LT.EPSI) THEN 142 : QSI(1,IG)=-D3/D1 143 : S1=A11*QSI(1,IG)+A13 144 : IF(ABS(S1).GT.EPSI) THEN 145 : QSI(2,IG)=-A14/S1 146 : GOTO 320 147 : ELSE 148 : QSI(2,IG)=-A24/(A21*QSI(1,IG)+A23) 149 : GOTO 320 150 : END IF 151 : END IF 152 : A1=-(A12*A21-A22*A11)/D2 153 : A2=-(A14*A21-A24*A11)/D2 154 : A=A11*A1 155 : B=A11*A2+A12+A13*A1 156 : C=A13*A2+A14 157 : IF(ABS(A).GT.EPSI) THEN 158 : cs=B*B-4.0*A*C 159 : BC=SQRT(B*B-4.0*A*C) 160 : QS1=0.5*(-B+BC)/A 161 : QS2=0.5*(-B-BC)/A 162 : ET1=A1*QS1+A2 163 : ET2=A1*QS2+A2 164 : IF(ABS(QS1).LT.1.0.AND.ABS(ET1).LT.1.0) THEN 165 : QSI(1,IG)=QS1 166 : QSI(2,IG)=ET1 167 : ELSE 168 : QSI(1,IG)=QS2 169 : QSI(2,IG)=ET2 170 : END IF 171 : GOTO 320 172 : ELSE 173 : QSI(1,IG)=-C/B 174 : QSI(2,IG)=A1*QSI(1,IG)+A2 175 : GOTO 320 176 : END IF 177 : 320 CONTINUE 178 : IF(ICM.NE.0) THEN 179 : T1=0.5D0*(XEC(1,2)-XEC(1,1)) 180 : T2=0.5D0*(XEC(2,2)-XEC(2,1)) 181 : DJ=SQRT(T1*T1+T2*T2) 182 : C --- COS BETA 183 : SIB(1)=T1/DJ 184 : C --- SIN BETA 185 : SIB(2)=T2/DJ 186 : do 321 ig=1,igau 187 : sga(1,4,ig)=sib(1) 188 : sga(2,4,ig)=sib(2) 189 : 321 sga(1,5,ig)=dj 190 : END IF 191 : C --- Calcul de fonctions de forme 192 : DO 330 IG=1,IGAU 193 : SHPC(1,1,IG)=0.25D0*(1.D0-QSI(1,IG))*(1.D0-QSI(2,IG)) 194 : SHPC(1,2,IG)=0.25D0*(1.D0+QSI(1,IG))*(1.D0-QSI(2,IG)) 195 : SHPC(1,3,IG)=0.25D0*(1.D0+QSI(1,IG))*(1.D0+QSI(2,IG)) 196 : 330 SHPC(1,4,IG)=0.25D0*(1.D0-QSI(1,IG))*(1.D0+QSI(2,IG)) 197 : IF(ICM.EQ.0) GOTO 4000 198 : DO 350 IG=1,IGAU 199 : C --- Derivees par rapport a QSI 200 : SHPC(2,1,IG)=0.25D0*(-1.D0+QSI(2,IG)) 201 : SHPC(2,2,IG)=0.25D0*( 1.D0-QSI(2,IG)) 202 : SHPC(2,3,IG)=0.25D0*( 1.D0+QSI(2,IG)) 203 : SHPC(2,4,IG)=0.25D0*(-1.D0-QSI(2,IG)) 204 : C --- Derivees par rapport a ETA 205 : SHPC(3,1,IG)=0.25D0*(-1.D0+QSI(1,IG)) 206 : SHPC(3,2,IG)=0.25D0*(-1.D0-QSI(1,IG)) 207 : SHPC(3,3,IG)=0.25D0*( 1.D0+QSI(1,IG)) 208 : 350 SHPC(3,4,IG)=0.25D0*( 1.D0-QSI(1,IG)) 209 : DO 360 IG=1,IGAU 210 : A11=0.0 211 : A12=0.0 212 : A21=0.0 213 : A22=0.0 214 : DO 355 IS=1,4 215 : A11=A11+XEL(1,IS)*SHPC(2,IS,IG) 216 : A12=A12+XEL(1,IS)*SHPC(3,IS,IG) 217 : A21=A21+XEL(2,IS)*SHPC(2,IS,IG) 218 : 355 A22=A22+XEL(2,IS)*SHPC(3,IS,IG) 219 : DET=A11*A22-A12*A21 220 : S1= (T1*A22-T2*A12)/DET 221 : S2=-(T1*A21-T2*A11)/DET 222 : C --- Derivees de QSI par rapport a MIU 223 : SGA(1,3,IG)=S1 224 : C --- Derivees de ETA par rapport a MIU 225 : SGA(2,3,IG)=S2 226 : 360 CONTINUE 227 : GOTO 4000 228 : END IF 229 : C 230 : C --- QUA8 231 : C 232 : IF(MELE.EQ.10) THEN 233 : DO 374 IG=1,IGAU 233 : CALL QUADRI(ICM,2,XEL,XEC,XX,SHPC,SGA,QSI,IG) 234 : 374 CONTINUE 235 : GOTO 5000 236 : END IF 237 : C 238 : C --- CUB8 239 : C 240 : IF(MELE.EQ.14) THEN 241 : DO 381 IG=1,IGAU 241 : CALL CUBEIT(ICM,XEL,XEC,XX,SHPC,SGA,QSI,IG) 242 : 381 CONTINUE 243 : GOTO 5000 244 : END IF 245 : C 246 : C --- CU20 247 : C 248 : IF(MELE.EQ.15) THEN 249 : DO 391 IG=1,IGAU 249 : CALL CUBTWE(ICM,XEL,XEC,XX,SHPC,SGA,QSI,IG) 250 : 391 CONTINUE 251 : GOTO 5000 252 : END IF 253 : C 254 : C --- PRI6 255 : C 256 : IF(MELE.EQ.16) THEN 257 : DO 401 IG=1,IGAU 257 : CALL PRISIX(ICM,XEL,XEC,XX,SHPC,SGA,QSI,IG) 258 : 401 CONTINUE 259 : GOTO 5000 260 : END IF 261 : C 262 : C --- PR15 263 : C 264 : IF(MELE.EQ.17) THEN 265 : DO 411 IG=1,IGAU 265 : CALL PRIS15(ICM,XEL,XEC,XX,SHPC,SGA,QSI,IG) 266 : 411 CONTINUE 267 : GOTO 5000 268 : END IF 269 : C 270 : 4000 CONTINUE 271 : DO 2000 IS=1,IDIM 272 : C --- Coordonnees du point d'integration dans le systeme reel (x,y,z) 273 : SGA(IS,1,IG)=XX(IS,IG) 274 : C --- Coordonnees du point d'integration dans le systeme de reference 275 : C (qsi,eta, ) 276 : SGA(IS,2,IG)=QSI(IS,IG) 277 : 2000 CONTINUE 278 : 5000 CONTINUE 279 : END 280 : 281 :
© Cast3M 2003 - All rights reserved.
Disclaimer