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