1 : C ZCLPLS SOURCE CHAT 05/01/13 04:21:27 5004 2 : SUBROUTINE ZCLPLS(HR,RPG,LE,NEL,K0,NPTD,NX,IES,NP,IAXI,IPADL, 3 : & COEFF,IK1, 4 : & TN,G, 5 : & VOLU,COTE,NELZ,DME,DIAM, 6 : & DT,DTT2,NUEL,DIAEL) 7 : C 8 : C VERSION VECTORISEE (CF XCVTIT POUR PLUS DE DETAILS) 9 : C 10 : IMPLICIT INTEGER(I-N) 11 : IMPLICIT REAL*8 (A-H,O-Z) 12 : C NOMBRE MAXI DE POINTS PAR ELEMENT : NPX 13 : PARAMETER(NPX=9) 14 : C LONGUEUR DES REGISTRES VECTORIELS DE LA MACHINE CIBLE 15 : PARAMETER(LRV=64) 16 : C*********************************************************************** 17 : C 18 : C CE SP DISCRETISE LE LAPLACIEN D UNE VARIABLE SCALAIRE 19 : C 20 : C EN 1D SUR L ELEMENT SEG2 21 : C EN 2D SUR L ELEMENT QUA4 PLAN OU AXI 22 : C EN 3D SUR L ELEMENT CUB8 PRI6 et TET4 23 : C LES OPERATEURS SONT SOUS INTEGRES 24 : C 25 : C COEFF(SCAL DOMA) LE COEFFICIENT DU LAPLACIEN ( NU ) 26 : C (SCAL ELEM) 27 : C 28 : C NUELD(NELZ): NUMERO DE L'ELEMENT NK DE LA ZONE DANS LE DOMAINE 29 : C (NELZ : NOMBRE D'ELEMENTS DE LA ZONE) 30 : C ROC(NELD) : ROC PAR ELEMENT DONNE SUR TOUT LE DOMAINE 31 : C 32 : C 33 : C*********************************************************************** 34 : -INC CCVQUA4 35 : -INC CCREEL 36 : C 37 : DIMENSION TN(NPTD,NX),GG(8) 38 : DIMENSION COEFF(*) 39 : DIMENSION COTE(NELZ,IES),DME(*),VOLU(*),DIAM(*) 40 : 41 : DIMENSION IPADL(1),LE(NP,1) 42 : DIMENSION HR(NEL,NP,IES),RPG(1) 43 : 44 : REAL*8 G(NPTD,NX) 45 : REAL*8 ROC(1) 46 : DIMENSION NUELD(1) 47 : DIMENSION QGGT(8,8),Q1(8,8),Q2(8,8),Q3(8,8) 48 : C 49 : C 50 : 51 : C* -TABLEAUX ADDITIONNELS POUR LA VECTORISATION ********************** 52 : DIMENSION AIRE(LRV),ALF (LRV),COEF(LRV),CLSR(LRV),GINC(LRV) 53 : DIMENSION AL (LRV),AH (LRV),AP (LRV),CFM (LRV),XMA (LRV) 54 : DIMENSION XMB(LRV),XMD(LRV),XMH(LRV),XMI(LRV),AHL(LRV) 55 : DIMENSION ALH(LRV),DPR(LRV),TETA(LRV,NPX,3),F(LRV,NPX),BF(LRV,NPX) 56 : C 57 : C --- TABLEAU POUR L'OPTION RAPIDE --- 58 : DIMENSION ALP(LRV) 59 : C*** 60 : SAVE IPAS,QGGT,Q1,Q2,Q3 61 : DATA IPAS/0/ 62 : 63 : C 64 : C INITIALISATIONS DIVERSES 65 : C 66 : C WRITE(6,*)' ROC=',(ROC(MM),MM=1,NEL) 67 : C WRITE(6,*)' NUELD=',(NUELD(MM),MM=1,NEL) 68 : 69 : NK=K0 69 : CALL INITD(BF,64,0.D0) 70 : 71 : IF(IES.EQ.2)GO TO 20 72 : IF(IES.EQ.3)GO TO 30 73 : C********** 74 : C * 1D* 75 : C********** 76 : C 77 : C 78 : C K EST LE NUMERO DE L'ELEMENT 79 : C KP PERMET DE SE SITUER A L'INTERIEUR DES TABLEAUX DE LRV ELEMENTS 80 : NPACK=INT(NEL/FLOAT(LRV))+1 81 : KPACKD=1 82 : KPACKF=NPACK 83 : C 84 : C*** BOUCLE SUR LES ELEMENTS *** 85 : C 86 : DO 70001 KPACK=KPACKD,KPACKF 87 : C DO 40 K=1,NEL 88 : C 89 : C POUR CHAQUE PAQUET DE LRV ELEMENTS 90 : KDEB=1+(KPACK-1)*LRV 91 : KFIN=MIN(NEL,KDEB+LRV-1) 92 : DO 70002 K=KDEB,KFIN 93 : KP=K-KDEB+1 94 : NK=K+K0 95 : K1=1+(1-IK1)*(NK-1) 96 : COEF(KP)=COEFF(K1) 97 : CLSR(KP)=COEFF(K1) 98 : AIRE(KP)=VOLU(NK) 99 : XMA(KP)=AIRE(KP) 100 : 70002 CONTINUE 101 : DO 4 I=1,NP 102 : DO 4 N=1,NX 103 : DO 4 K=KDEB,KFIN 104 : KP=K-KDEB+1 105 : NF=IPADL(LE(I,K)) 106 : TETA(KP,I,N)=TN(NF,N) 107 : 4 CONTINUE 108 : C 109 : C 110 : C 111 : DO 70003 K=KDEB,KFIN 112 : NK=K+K0 113 : KP=K-KDEB+1 114 : DT0=DT 115 : DT2=0.5*XMA(KP)*XMA(KP)/CLSR(KP) 116 : IF(DT2.LT.DT)DT=DT2 117 : IF(DT.EQ.DT0)GO TO 152 118 : DTT2=DT2 119 : DIAEL=XMA(KP) 120 : NUEL=NK 121 : 152 CONTINUE 122 : 70003 CONTINUE 123 : C 124 : C 125 : N=1 126 : C --- CALCUL DE L'INCREMENT GINC --------------------------- 127 : DO 70104 K=KDEB,KFIN 128 : KP=K-KDEB+1 129 : GINC(KP)=COEF(KP)*(TETA(KP,1,N)-TETA(KP,2,N))/AIRE(KP) 130 : 70104 CONTINUE 131 : C --- ACCUMULATION DANS LE TABLEAU G ----------------------- 132 : DO 70004 K=KDEB,KFIN 133 : KP=K-KDEB+1 134 : NF=IPADL(LE(1,K)) 135 : G(NF,N)=G(NF,N) + GINC(KP) 136 : NF=IPADL(LE(2,K)) 137 : G(NF,N)=G(NF,N) - GINC(KP) 138 : 70004 CONTINUE 139 : 70001 CONTINUE 140 : 40 CONTINUE 141 : IPAS=1 142 : RETURN 143 : C *********** 144 : C * 2D * 145 : C *********** 146 : 20 CONTINUE 147 : C 148 : 149 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 150 : 151 : C DIFFERENCES TRIANGLE / QUADRANGLE 152 : QUA4=0.D0 153 : IF(NP.EQ.4)QUA4=1.D0 154 : 155 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 156 : 157 : 158 : C 159 : C K EST LE NUMERO DE L'ELEMENT 160 : C KP PERMET DE SE SITUER A L'INTERIEUR DES TABLEAUX DE LRV ELEMENTS 161 : NPACK=INT(NEL/FLOAT(LRV))+1 162 : KPACKD=1 163 : KPACKF=NPACK 164 : C 165 : C*** BOUCLE SUR LES ELEMENTS *** 166 : C 167 : DO 80001 KPACK=KPACKD,KPACKF 168 : C 169 : C POUR CHAQUE PAQUET DE LRV ELEMENTS 170 : KDEB=1+(KPACK-1)*LRV 171 : KFIN=MIN(NEL,KDEB+LRV-1) 172 : DO 80002 K=KDEB,KFIN 173 : KP=K-KDEB+1 174 : NK=K+K0 175 : K1=1+(1-IK1)*(NK-1) 176 : XMI(KP)=DIAM(NK) 177 : XMA(KP)=DME(NK) 178 : COEF(KP)=COEFF(K1) 179 : CLSR(KP)=COEFF(K1) 180 : AIRE(KP)=VOLU(NK) 181 : CC 182 : AL(KP)=COTE(NK,1) 183 : AH(KP)=COTE(NK,2) 184 : XMB(KP)=(AL(KP)+AH(KP))/2. 185 : XMD(KP)=((AL(KP)*AH(KP))**2)/(AL(KP)**2+AH(KP)**2) 186 : AHL(KP)=AH(KP)/AL(KP) 187 : ALH(KP)=AL(KP)/AH(KP) 188 : DPR(KP)=1.D0 189 : IF(IAXI.NE.0)DPR(KP)=2.D0*XPI*RPG(K) 190 : 80002 CONTINUE 191 : CC 192 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 193 : C 194 : DO 5 I=1,NP 195 : DO 5 N=1,NX 196 : DO 5 K=KDEB,KFIN 197 : KP=K-KDEB+1 198 : NF=IPADL(LE(I,K)) 199 : TETA(KP,I,N)=TN(NF,N) 200 : 5 CONTINUE 201 : C 202 : C 203 : C 204 : DO 80004 K=KDEB,KFIN 205 : KP=K-KDEB+1 206 : NK=K+K0 207 : DT0=DT 208 : DT2=0.5*XMD(KP)/CLSR(KP) 209 : IF(DT2.LT.DT)DT=DT2 210 : IF(DT.EQ.DT0)GO TO 252 211 : DTT2=DT2 212 : DIAEL=XMB(KP) 213 : NUEL=NK 214 : 252 CONTINUE 215 : 80004 CONTINUE 216 : C 217 : C --- VERSION RAPIDE ----------------------- 218 : CF12=1.D0/12.D0 219 : DO 80011 K=KDEB,KFIN 220 : KP=K-KDEB+1 221 : ALP(KP)=CF12*(AHL(KP)+ALH(KP))*DPR(KP)*QUA4 222 : 80011 CONTINUE 223 : C ------------------------------------------ 224 : C 225 : DO 1 N=1,NX 226 : DO 80010 I= 1,NP 227 : DO 80010 K=KDEB,KFIN 228 : KP=K-KDEB+1 229 : BF(KP,I)=0.D0 230 : 80010 CONTINUE 231 : 232 : DO 1 I=1,NP 233 : DO 3 J= 1,NP 234 : DO 3 K=KDEB,KFIN 235 : KP=K-KDEB+1 236 : BF(KP,I)=BF(KP,I)+TETA(KP,J,N)* 237 : & ((HR(K,I,1)*HR(K,J,1)+HR(K,I,2)*HR(K,J,2))*AIRE(KP) 238 : & +ALP(KP)*VGGT(J,I))*COEF(KP) 239 : 3 CONTINUE 240 : 241 : C --- ACCUMULATION DANS LE TABLEAU G ------------------------------ 242 : DO 80006 K=KDEB,KFIN 243 : KP=K-KDEB+1 244 : NF=IPADL(LE(I,K)) 245 : G(NF,N)=G(NF,N)+BF(KP,I) 246 : 80006 CONTINUE 247 : 1 CONTINUE 248 : 80001 CONTINUE 249 : 50 CONTINUE 250 : 251 : C WRITE(6,*)' SUB XCLPLS G(1,=' 252 : C WRITE(6,1002) (G(MM,1),MM=1,NPTD) 253 : C WRITE(6,*)' SUB XCLPLS G(2,=' 254 : C WRITE(6,1002) (G(MM,2),MM=1,NPTD) 255 : C WRITE(6,*)' FIN **** ' 256 : 257 : IPAS=1 258 : RETURN 259 : 30 CONTINUE 260 : C *********** 261 : C * 3D * 262 : C *********** 263 : 263 : IF(IPAS.EQ.0)CALL CALHRH(QGGT,Q1,Q2,Q3,IES) 264 : CUB8=0.D0 265 : IF(NP.EQ.8)CUB8=1.D0 266 : 267 : 1003 FORMAT(' XCVTI ',10I10) 268 : C 269 : C K EST LE NUMERO DE L'ELEMENT 270 : C KP PERMET DE SE SITUER A L'INTERIEUR DES TABLEAUX DE LRV ELEMENTS 271 : NPACK=INT(NEL/FLOAT(LRV))+1 272 : KPACKD=1 273 : KPACKF=NPACK 274 : C 275 : C*** BOUCLE SUR LES ELEMENTS *** 276 : C 277 : DO 90001 KPACK=KPACKD,KPACKF 278 : C DO 60 K=1,NEL 279 : C 280 : C POUR CHAQUE PAQUET DE LRV ELEMENTS 281 : KDEB=1+(KPACK-1)*LRV 282 : KFIN=MIN(NEL,KDEB+LRV-1) 283 : DO 90002 K=KDEB,KFIN 284 : KP=K-KDEB+1 285 : NK=K+K0 286 : K1=1+(1-IK1)*(NK-1) 287 : AL(KP)=COTE(NK,1) 288 : AH(KP)=COTE(NK,2) 289 : AP(KP)=COTE(NK,3) 290 : CFM(KP)=AL(KP)*AH(KP)/AP(KP)+AL(KP)*AP(KP)/AH(KP)+ 291 : & AP(KP)*AH(KP)/AL(KP) 292 : XMI(KP)=DIAM(NK) 293 : XMA(KP)=DME(NK) 294 : COEF(KP)=COEFF(K1) 295 : C CLSR(KP)=COEFF(K1)/ROC(NUELD(NK)) 296 : CLSR(KP)=COEFF(K1) 297 : AIRE(KP)=VOLU(NK) 298 : XMB(KP)=(XMA(KP)+XMI(KP))*0.5 299 : XMD(KP)=((XMA(KP)*XMI(KP))**2)/(2.*XMA(KP)**2+XMI(KP)**2) 300 : XMH(KP)=AIRE(KP)**0.33 301 : 90002 CONTINUE 302 : C 303 : DO 15 I=1,NP 304 : DO 15 N=1,NX 305 : DO 15 K=KDEB,KFIN 306 : KP=K-KDEB+1 307 : NF=IPADL(LE(I,K)) 308 : TETA(KP,I,N)=TN(NF,N) 309 : 15 CONTINUE 310 : C WRITE(6,*)' K ** ',K,' TETA ',' VOLU=',AIRE,' COEF=',COEF 311 : C &,' XMH=',XMH 312 : C WRITE(6,1002)(TETA(MM,1),MM=1,8) 313 : C 314 : C 315 : DO 90003 K=KDEB,KFIN 316 : KP=K-KDEB+1 317 : NK=K+K0 318 : DT0=DT 319 : DT2=0.5*XMD(KP)/CLSR(KP) 320 : IF(DT2.LT.DT)DT=DT2 321 : IF(DT.EQ.DT0)GO TO 353 322 : DTT2=DT2 323 : DIAEL=XMB(KP) 324 : NUEL=NK 325 : 353 CONTINUE 326 : 90003 CONTINUE 327 : C 328 : DO 11 N=1,NX 329 : DO 90004 I=1,NP 330 : DO 90004 K=KDEB,KFIN 331 : KP=K-KDEB+1 332 : BF(KP,I)=0. 333 : 90004 CONTINUE 334 : DO 11 I=1,NP 335 : DO 13 J=1,NP 336 : DO 13 K=KDEB,KFIN 337 : KP=K-KDEB+1 338 : GEO1=CFM(KP)*QGGT(J,I)*CUB8 339 : BF(KP,I)=BF(KP,I)+TETA(KP,J,N)* 340 : &(HR(K,I,1)*HR(K,J,1)+HR(K,I,2)*HR(K,J,2)+HR(K,I,3)*HR(K,J,3)) 341 : & *AIRE(KP)+COEF(KP)*XMH(KP)*GEO1 342 : 13 CONTINUE 343 : DO 90006 K=KDEB,KFIN 344 : KP=K-KDEB+1 345 : NF=IPADL(LE(I,K)) 346 : G(NF,N)=G(NF,N)+BF(KP,I) 347 : 90006 CONTINUE 348 : 11 CONTINUE 349 : 90001 CONTINUE 350 : C WRITE(6,*)' GG(I)=' 351 : C WRITE(6,1002)GG 352 : 60 CONTINUE 352 : C CALL ARRET(0) 353 : IPAS=1 354 : RETURN 355 : 1001 FORMAT(' XCVTI',I10,6E12.5) 356 : 1002 FORMAT(10(1X,1PE11.4)) 357 : END 358 : 359 : 360 : 361 : 362 : 363 : 364 : 365 :
© Cast3M 2003 - All rights reserved.
Disclaimer