1 : C SHAPX SOURCE BP208322 13/11/13 21:15:23 7865 2 : C SHAPX 3 : c 4 : SUBROUTINE SHAPX(XG,YG,ZG,MELE,LV7,NBRMAX,NBENRJ, 5 : &TLREEL,JX,SHP,IRET) 6 : c SUBROUTINE SHAPX(XG,YG,MELE,NFTIP,NSING,ITABX,JX,SHP,IRET) 7 : C 8 : C======================================================================= 9 : C 10 : C Entree: SHP = FONCTIONS DE FORME STANDARD ET LEUR DERIVEES 11 : C LV7 = FONCTIONS DE NIVEAU ET LEUR DERIVEES 12 : C LV7(ienr, PSI/PHI, valeur/deriveeparqsi/pareta, iNOEUD) 13 : C Sortie: SHP = TOUTES LES FONCTIONS DE FORME DE L'EF (std + enrichissement) 14 : c ET LEUR DERIVEES DANS LE REPERE GLOBAL 15 : C 16 : C 17 : C=====PARTIE DECLARATIVE================================================ 18 : C 19 : implicit real*8(a-h,o-z) 20 : implicit integer(i-n) 21 : -INC CCOPTIO 22 : -INC CCREEL 23 : -INC SMLREEL 24 : C 25 : DIMENSION SHP(6,*) 26 : c PARAMETER (NBRMAX=5) 27 : REAL*8 LV7(NBRMAX,2,6,*) 28 : INTEGER TLREEL(NBRMAX,*) 29 : c DIMENSION LV7(NBRMAX,2,6,*) 30 : c DIMENSION TLREEL(NBRMAX,*) 31 : REAL*8 PSI,PHI,H,R,THETA,XG,YG,ZG 32 : C 33 : C 34 : c 35 : c WRITE(*,*) '#### ENTREE DANS SHAPX ####' 36 : C 37 : c 38 : C=====AIGUILLAGE selon element========================================== 39 : c 40 : if(MELE.eq.263) GOTO 263 41 : if(MELE.eq.264) GOTO 264 42 : c Erreur pour les elements non programmés 42 : call erreur (21) 43 : RETURN 44 : C 45 : CTY if(MELE.eq.263) GOTO 263 46 : CTY Erreur pour les elements non programmés 46 : CTY call erreur (21) 47 : CTY RETURN 48 : c 49 : c 50 : C=====XQ4R============================================================== 51 : c 52 : 263 continue 53 : NBNN = 4 54 : c KMAX =LV7(/1) 55 : c 56 : c 57 : C++++ ENRICHISSEMENT H ++++++++++++++++++++++++++++++ 58 : KENR = 1 59 : 60 : C>>>> boucle sur les noeuds 61 : DO 200 I=1,NBNN 62 : 63 : MLREEL=TLREEL(KENR,I) 64 : if(MLREEL.eq.0) goto 200 65 : 66 : C on calcule les fonctions utiles 67 : PHI = LV7(KENR,2,1,I) 68 : H = SIGN(1.D0,PHI) 69 : C 70 : c on modifie les fonctions de forme et leurs derivees 71 : II = (KENR*NBNN) + I 72 : DO IND=1,3 73 : SHP(IND,II) = H*SHP(IND,I) 74 : ENDDO 75 : C 76 : 200 CONTINUE 77 : C<<<<<<<<<<<<< 78 : C 79 : C 80 : C++++ ENRICHISSEMENT F1 et F2 ++++++++++++++++++++++++ 81 : C 82 : c DO 300 KENR=2,NBRMAX 83 : DO 300 KENR=2,NBENRJ 84 : IENR = 4*(KENR-2) + 2 85 : 86 : C>>>>>>> boucle sur les noeuds 87 : DO 305 I=1,NBNN 88 : C 89 : MLREEL=TLREEL(KENR,I) 90 : c write(*,*) 'JX,KENR,I=',JX,KENR,I,' MLREEL=',MLREEL 91 : if(MLREEL.eq.0) goto 305 92 : 93 : C on recupere les level set et leur derivees 94 : PSI = LV7(KENR,1,1,I) 95 : PSIX = LV7(KENR,1,2,I) 96 : PSIY = LV7(KENR,1,3,I) 97 : PHI = LV7(KENR,2,1,I) 98 : PHIX = LV7(KENR,2,2,I) 99 : PHIY = LV7(KENR,2,3,I) 100 : 101 : C on calcule les fonctions utiles 102 : H = SIGN(1.D0,PHI) 103 : R = ( (PSI**2) + (PHI**2) )**0.5 104 : THETA = H * ((XPI/2.) - (ATAN(PSI/ABS(PHI)))) 105 : SIN1T = SIN(THETA) 106 : COS1T = COS(THETA) 107 : SIN05T = SIN(THETA/2.) 108 : COS05T = COS(THETA/2.) 109 : C 110 : C et leurs derivées 111 : C RX = ( (PSI*PSIX) + (PHI*PHIX) ) / (R**0.5) 112 : C RY = ( (PSI*PSIY) + (PHI*PHIY) ) / (R**0.5) 113 : C THETAX= ( ((-1.*ABS(PHI))*PSIX) + ((H*PSI)*PHIX) ) / (R**2) 114 : C THETAY= ( ((-1.*ABS(PHI))*PSIY) + ((H*PSI)*PHIY) ) / (R**2) 115 : 116 : c on calcule les fonctions de forme d ENRICHISSEMENT + derivees 117 : 118 : C+++ Fb = (r**0.5)*sin(theta/2) 119 : II = (IENR*NBNN) + I 120 : 121 : C Fonction de forme Ni*Fb 122 : FB = (R**0.5) * SIN05T 123 : SHP(1,II) = SHP(1,I) * FB 124 : C Derivees (Ni*Fb),X = (Ni,X*Fb) + (Ni*Fb,X) 125 : SHP(2,II) = ( SHP(2,I) * FB ) + 126 : $ ( SHP(1,I)*0.5/(R**0.5) * 127 : $ ( ( SIN05T * ((COS1T*PSIX)+(SIN1T*PHIX)) ) 128 : $ + ( COS05T * ((-1.*SIN1T*PSIX)+(COS1T*PHIX)) ) ) ) 129 : C Derivees (Ni*Fb),y = (Ni,y*Fb) + (Ni*Fb,y) 130 : SHP(3,II) = ( SHP(3,I) * FB ) + 131 : $ ( SHP(1,I)*0.5/(R**0.5) * 132 : $ ( ( SIN05T * ((COS1T*PSIY)+(SIN1T*PHIY)) ) 133 : $ + ( COS05T * ((-1.*SIN1T*PSIY)+(COS1T*PHIY)) ) ) ) 134 : 135 : 136 : C+++ Fc = (r**0.5)*cos(theta/2) 137 : II = ((IENR+1)*NBNN) + I 138 : 139 : C Fonction de forme Ni*Fc 140 : FC = (R**0.5) * COS05T 141 : SHP(1,II) = SHP(1,I) * FC 142 : C Derivees (Ni*Fc),x = (Ni,x*Fc) + (Ni*Fc,x) 143 : SHP(2,II) = ( SHP(2,I) * FC ) + 144 : $ ( SHP(1,I)*0.5/(R**0.5) * 145 : $ ( ( COS05T * ((COS1T*PSIX)+(SIN1T*PHIX)) ) 146 : $ - ( SIN05T * ((-1.*SIN1T*PSIX)+(COS1T*PHIX)) ) ) ) 147 : C Derivees (Ni*Fc),y = (Ni,y*Fc) + (Ni*Fc,y) 148 : SHP(3,II) = ( SHP(3,I) * FC) + 149 : $ ( SHP(1,I)*0.5/(R**0.5) * 150 : $ ( ( COS05T * ((COS1T*PSIY)+(SIN1T*PHIY)) ) 151 : $ - ( SIN05T * ((-1.*SIN1T*PSIY)+(COS1T*PHIY)) ) ) ) 152 : 153 : 154 : C+++ Fd = (r**0.5)*sin(theta/2)*sin(theta) 155 : II = ((IENR+2)*NBNN) + I 156 : 157 : C Fonction de forme Ni*Fd 158 : FD = (R**0.5) * SIN05T * SIN1T 159 : SHP(1,II)= SHP(1,I) * FD 160 : C Derivees (Ni*Fc),X = (Ni,X*Fc) + (Ni*Fc,X) 161 : SHP(2,II) = ( SHP(2,I) * FD ) + 162 : $ ( SHP(1,I)*0.5/(R**0.5) * 163 : $ ( ( SIN05T*SIN1T * ((COS1T*PSIX)+(SIN1T*PHIX)) ) 164 : $ + ( (COS05T*SIN1T + (2.*SIN05T*COS1T)) * 165 : $ ((-1.*SIN1T*PSIX)+(COS1T*PHIX)) ) ) ) 166 : C Derivees (Ni*Fc),Y = (Ni,Y*Fc) + (Ni*Fc,Y) 167 : SHP(3,II) = ( SHP(3,I) * FD ) + 168 : $ ( SHP(1,I)*0.5/(R**0.5) * 169 : $ ( ( SIN05T*SIN1T * ((COS1T*PSIY)+(SIN1T*PHIY)) ) 170 : $ + ( (COS05T*SIN1T + (2.*SIN05T*COS1T)) * 171 : $ ((-1.*SIN1T*PSIY)+(COS1T*PHIY)) ) ) ) 172 : 173 : 174 : C+++ Fe = (r**0.5)*cos(theta/2)*sin(theta) 175 : II = ((IENR+3)*NBNN) + I 176 : 177 : C Fonction de forme Ni*Fe 178 : FE = (R**0.5) * COS05T * SIN1T 179 : SHP(1,II)= SHP(1,I) * FE 180 : C Derivees (Ni*Fc),x = (Ni,x*Fe) + (Ni*Fe,x) 181 : SHP(2,II) = ( SHP(2,I) * FE ) + 182 : $ ( SHP(1,I)*0.5/(R**0.5) * 183 : $ ( ( COS05T*SIN1T * ((COS1T*PSIX)+(SIN1T*PHIX)) ) 184 : $ - ( (SIN05T*SIN1T - (2.*COS05T*COS1T)) * 185 : $ ((-1.*SIN1T*PSIX)+(COS1T*PHIX)) ) ) ) 186 : C Derivees (Ni*Fe),y = (Ni,y*Fe) + (Ni*Fe,y) 187 : SHP(3,II) = ( SHP(3,I) * FE ) + 188 : $ ( SHP(1,I)*0.5/(R**0.5) * 189 : $ ( ( COS05T*SIN1T * ((COS1T*PSIY)+(SIN1T*PHIY)) ) 190 : $ - ( (SIN05T*SIN1T - (2.*COS05T*COS1T)) * 191 : $ ((-1.*SIN1T*PSIY)+(COS1T*PHIY)) ) ) ) 192 : 193 : C 194 : C 195 : 305 CONTINUE 196 : C<<<<<<<<<<<<< 197 : 198 : 300 CONTINUE 199 : 200 : GOTO 999 201 : 202 : C_____________________________________________ 203 : C FONCTIONS DE FORME TRIDIMENSIONNELLES 204 : C 205 : C=====XC8R============================================================== 206 : c 207 : 264 continue 208 : C 209 : NBNN = 8 210 : c KMAX =LV7(/1) 211 : c 212 : c 213 : C++++ ENRICHISSEMENT H ++++++++++++++++++++++++++++++ 214 : KENR = 1 215 : 216 : C>>>> boucle sur les noeuds 217 : DO 400 I=1,NBNN 218 : 219 : MLREEL=TLREEL(KENR,I) 220 : if(MLREEL.eq.0) goto 400 221 : 222 : C on calcule les fonctions utiles 223 : PHI = LV7(KENR,2,1,I) 224 : H = SIGN(1.D0,PHI) 225 : C 226 : c on modifie les fonctions de forme et leurs derivees 227 : II = (KENR*NBNN) + I 228 : DO IND=1,4 229 : SHP(IND,II) = H*SHP(IND,I) 230 : ENDDO 231 : C 232 : 400 CONTINUE 233 : C<<<<<<<<<<<<< 234 : C 235 : C 236 : C++++ ENRICHISSEMENT F1 et F2 ++++++++++++++++++++++++ 237 : C 238 : c DO 500 KENR=2,NBRMAX 239 : DO 500 KENR=2,NBENRJ 240 : IENR = 4*(KENR-2) + 2 241 : 242 : C>>>>>>> boucle sur les noeuds 243 : DO 505 I=1,NBNN 244 : C 245 : MLREEL=TLREEL(KENR,I) 246 : c write(*,*) 'JX,KENR,I=',JX,KENR,I,' MLREEL=',MLREEL 247 : if(MLREEL.eq.0) goto 505 248 : 249 : C on recupere les level set et leur derivees 250 : PSI = LV7(KENR,1,1,I) 251 : PSIX = LV7(KENR,1,2,I) 252 : PSIY = LV7(KENR,1,3,I) 253 : PSIZ = LV7(KENR,1,4,I) 254 : PHI = LV7(KENR,2,1,I) 255 : PHIX = LV7(KENR,2,2,I) 256 : PHIY = LV7(KENR,2,3,I) 257 : PHIZ = LV7(KENR,2,4,I) 258 : 259 : C on calcule les fonctions utiles 260 : H = SIGN(1.D0,PHI) 261 : R = ( (PSI**2) + (PHI**2) )**0.5 262 : THETA = H * ((XPI/2.) - (ATAN(PSI/ABS(PHI)))) 263 : SIN1T = SIN(THETA) 264 : COS1T = COS(THETA) 265 : SIN05T = SIN(THETA/2.) 266 : COS05T = COS(THETA/2.) 267 : C 268 : C et leurs derivées 269 : C RX = ( (PSI*PSIX) + (PHI*PHIX) ) / (R**0.5) 270 : C RY = ( (PSI*PSIY) + (PHI*PHIY) ) / (R**0.5) 271 : C THETAX= ( ((-1.*ABS(PHI))*PSIX) + ((H*PSI)*PHIX) ) / (R**2) 272 : C THETAY= ( ((-1.*ABS(PHI))*PSIY) + ((H*PSI)*PHIY) ) / (R**2) 273 : 274 : c on calcule les fonctions de forme d ENRICHISSEMENT + derivees 275 : 276 : C+++ Fb = (r**0.5)*sin(theta/2) 277 : II = (IENR*NBNN) + I 278 : 279 : C Fonction de forme Ni*Fb 280 : FB = (R**0.5) * SIN05T 281 : SHP(1,II) = SHP(1,I) * FB 282 : C Derivees (Ni*Fb),X = (Ni,X*Fb) + (Ni*Fb,X) 283 : SHP(2,II) = ( SHP(2,I) * FB ) + 284 : $ ( SHP(1,I)*0.5/(R**0.5) * 285 : $ ( ( SIN05T * ((COS1T*PSIX)+(SIN1T*PHIX)) ) 286 : $ + ( COS05T * ((-1.*SIN1T*PSIX)+(COS1T*PHIX)) ) ) ) 287 : C Derivees (Ni*Fb),y = (Ni,y*Fb) + (Ni*Fb,y) 288 : SHP(3,II) = ( SHP(3,I) * FB ) + 289 : $ ( SHP(1,I)*0.5/(R**0.5) * 290 : $ ( ( SIN05T * ((COS1T*PSIY)+(SIN1T*PHIY)) ) 291 : $ + ( COS05T * ((-1.*SIN1T*PSIY)+(COS1T*PHIY)) ) ) ) 292 : C Derivees (Ni*Fb),z = (Ni,z*Fb) + (Ni*Fb,z) 293 : SHP(4,II) = ( SHP(4,I) * FB ) + 294 : $ ( SHP(1,I)*0.5/(R**0.5) * 295 : $ ( ( SIN05T * ((COS1T*PSIZ)+(SIN1T*PHIZ)) ) 296 : $ + ( COS05T * ((-1.*SIN1T*PSIZ)+(COS1T*PHIZ)) ) ) ) 297 : C 298 : C+++ Fc = (r**0.5)*cos(theta/2) 299 : II = ((IENR+1)*NBNN) + I 300 : 301 : C Fonction de forme Ni*Fc 302 : FC = (R**0.5) * COS05T 303 : SHP(1,II) = SHP(1,I) * FC 304 : C Derivees (Ni*Fc),x = (Ni,x*Fc) + (Ni*Fc,x) 305 : SHP(2,II) = ( SHP(2,I) * FC ) + 306 : $ ( SHP(1,I)*0.5/(R**0.5) * 307 : $ ( ( COS05T * ((COS1T*PSIX)+(SIN1T*PHIX)) ) 308 : $ - ( SIN05T * ((-1.*SIN1T*PSIX)+(COS1T*PHIX)) ) ) ) 309 : C Derivees (Ni*Fc),y = (Ni,y*Fc) + (Ni*Fc,y) 310 : SHP(3,II) = ( SHP(3,I) * FC) + 311 : $ ( SHP(1,I)*0.5/(R**0.5) * 312 : $ ( ( COS05T * ((COS1T*PSIY)+(SIN1T*PHIY)) ) 313 : $ - ( SIN05T * ((-1.*SIN1T*PSIY)+(COS1T*PHIY)) ) ) ) 314 : C Derivees (Ni*Fc),z = (Ni,z*Fc) + (Ni*Fc,z) 315 : SHP(4,II) = ( SHP(4,I) * FC) + 316 : $ ( SHP(1,I)*0.5/(R**0.5) * 317 : $ ( ( COS05T * ((COS1T*PSIZ)+(SIN1T*PHIZ)) ) 318 : $ - ( SIN05T * ((-1.*SIN1T*PSIZ)+(COS1T*PHIZ)) ) ) ) 319 : 320 : 321 : C+++ Fd = (r**0.5)*sin(theta/2)*sin(theta) 322 : II = ((IENR+2)*NBNN) + I 323 : 324 : C Fonction de forme Ni*Fd 325 : FD = (R**0.5) * SIN05T * SIN1T 326 : SHP(1,II)= SHP(1,I) * FD 327 : C Derivees (Ni*Fd),X = (Ni,X*Fd) + (Ni*Fd,X) 328 : SHP(2,II) = ( SHP(2,I) * FD ) + 329 : $ ( SHP(1,I)*0.5/(R**0.5) * 330 : $ ( ( SIN05T*SIN1T * ((COS1T*PSIX)+(SIN1T*PHIX)) ) 331 : $ + ( (COS05T*SIN1T + (2.*SIN05T*COS1T)) * 332 : $ ((-1.*SIN1T*PSIX)+(COS1T*PHIX)) ) ) ) 333 : C Derivees (Ni*Fd),Y = (Ni,Y*Fd) + (Ni*Fd,Y) 334 : SHP(3,II) = ( SHP(3,I) * FD ) + 335 : $ ( SHP(1,I)*0.5/(R**0.5) * 336 : $ ( ( SIN05T*SIN1T * ((COS1T*PSIY)+(SIN1T*PHIY)) ) 337 : $ + ( (COS05T*SIN1T + (2.*SIN05T*COS1T)) * 338 : $ ((-1.*SIN1T*PSIY)+(COS1T*PHIY)) ) ) ) 339 : C Derivees (Ni*Fd),z = (Ni,z*Fd) + (Ni*Fd,z) 340 : SHP(4,II) = ( SHP(4,I) * FD ) + 341 : $ ( SHP(1,I)*0.5/(R**0.5) * 342 : $ ( ( SIN05T*SIN1T * ((COS1T*PSIZ)+(SIN1T*PHIZ)) ) 343 : $ + ( (COS05T*SIN1T + (2.*SIN05T*COS1T)) * 344 : $ ((-1.*SIN1T*PSIZ)+(COS1T*PHIZ)) ) ) ) 345 : 346 : 347 : C+++ Fe = (r**0.5)*cos(theta/2)*sin(theta) 348 : II = ((IENR+3)*NBNN) + I 349 : 350 : C Fonction de forme Ni*Fe 351 : FE = (R**0.5) * COS05T * SIN1T 352 : SHP(1,II)= SHP(1,I) * FE 353 : C Derivees (Ni*Fe),x = (Ni,x*Fe) + (Ni*Fe,x) 354 : SHP(2,II) = ( SHP(2,I) * FE ) + 355 : $ ( SHP(1,I)*0.5/(R**0.5) * 356 : $ ( ( COS05T*SIN1T * ((COS1T*PSIX)+(SIN1T*PHIX)) ) 357 : $ - ( (SIN05T*SIN1T - (2.*COS05T*COS1T)) * 358 : $ ((-1.*SIN1T*PSIX)+(COS1T*PHIX)) ) ) ) 359 : C Derivees (Ni*Fe),y = (Ni,y*Fe) + (Ni*Fe,y) 360 : SHP(3,II) = ( SHP(3,I) * FE ) + 361 : $ ( SHP(1,I)*0.5/(R**0.5) * 362 : $ ( ( COS05T*SIN1T * ((COS1T*PSIY)+(SIN1T*PHIY)) ) 363 : $ - ( (SIN05T*SIN1T - (2.*COS05T*COS1T)) * 364 : $ ((-1.*SIN1T*PSIY)+(COS1T*PHIY)) ) ) ) 365 : C Derivees (Ni*Fe),z = (Ni,z*Fe) + (Ni*Fe,z) 366 : SHP(4,II) = ( SHP(4,I) * FE ) + 367 : $ ( SHP(1,I)*0.5/(R**0.5) * 368 : $ ( ( COS05T*SIN1T * ((COS1T*PSIZ)+(SIN1T*PHIZ)) ) 369 : $ - ( (SIN05T*SIN1T - (2.*COS05T*COS1T)) * 370 : $ ((-1.*SIN1T*PSIZ)+(COS1T*PHIZ)) ) ) ) 371 : 372 : C 373 : C 374 : 505 CONTINUE 375 : C<<<<<<<<<<<<< 376 : 377 : 500 CONTINUE 378 : 379 : 380 : C_____________________________________________ 381 : C FIN DU PROGRAMME 382 : 999 CONTINUE 383 : C 384 : 385 : IRET = 1 386 : c WRITE(*,*) 'C EST FINI POUR SHAPX' 387 : 388 : RETURN 389 : END 390 : 391 : 392 : 393 :
© Cast3M 2003 - All rights reserved.
Disclaimer