1 : C PENTE3 SOURCE CHAT 05/01/13 02:12:24 5004 2 : SUBROUTINE PENTE3(NFAC,MELEFL,MLECEN,NCOMP,MPOCHP, 3 : & MPOGRA,MPOMIN,MPOMAX,MPOALP) 4 : C 5 : C************************************************************************ 6 : C 7 : C PROJET : CASTEM 2000 8 : C 9 : C NOM : PENTE3 10 : C 11 : C DESCRIPTION : Cette subroutine est appellée par la subroutine 12 : C PENTE1 (calcul du gradient d'un CHPOINT 2D de type 13 : C CENTRE) 14 : C Elle contient la partie du calcul du limiteur. 15 : C 16 : C LANGUAGE : ESOPE 2000 avec extensions CISI 17 : C 18 : C AUTEUR : A. BECCANTINI, DRN/DMT/SEMT/TTMF 19 : C AUTEUR (Modif.) : R. MOREL, DRN/DMT/SEMT/TTMF 20 : C 21 : C************************************************************************ 22 : C 23 : C 24 : C APPELES (E/S) : none 25 : C 26 : C APPELES (Calcul) : none 27 : C 28 : C 29 : C************************************************************************ 30 : C 31 : C ENTREES : NFAC : nombre de faces 32 : C 33 : C MELEFL : pointeur du MELEME 'FACEL' 34 : C 35 : C MLECEN : pointeur de MLENTI qui contient la table 36 : C numerotation global/local de CENTREs 37 : C 38 : C NCOMP : nombre de composantes de CHPOINT dont on veut 39 : C calculer les pentes 40 : C 41 : C MPOCHP : pointeur de MPOVAL de CHPOINT dont on veut le 42 : C gradient 43 : C 44 : C MPOGRA : pointeur de MPOVAL du gradient 45 : C 46 : C MPOMIN : pointeur de MPOVAL du minimum sur le stencil 47 : C 48 : C MPOMAX : pointeur de MPOVAL du maximum sur le stencil 49 : C 50 : C SORTIES : MPOALP : pointeur de MELVAL du limiteur 51 : C 52 : C************************************************************************ 53 : C 54 : C HISTORIQUE (Anomalies et modifications éventuelles) 55 : C 56 : C HISTORIQUE : Cree le 6-4-1998 57 : C 58 : C HISTORIQUE : Modifie le 20-10-1998 pour extension 3D 59 : C 60 : C************************************************************************ 61 : C 62 : IMPLICIT INTEGER(I-N) 63 : -INC CCOPTIO 64 : C 65 : C**** Variables de CCOPTIO 66 : C 67 : C INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR 68 : C & ,IOTER, IOLEC, IOIMP, IOCAR, IOACQ 69 : C & ,IOPER, IOSGB, IOGRA, IOSAU, IORES 70 : C & ,IECHO, IIMPI, IOSPI 71 : C & ,IDIM 72 : CC & ,MCOORD 73 : C & ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE 74 : C & ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU 75 : C & ,NORINC,NORVAL,NORIND,NORVAD 76 : C & ,NUCROU, IPSAUV 77 : C 78 : INTEGER NFAC,NCOMP 79 : & ,NLCF,NGCF,NGCEG,NLCEG,NGCED,NLCED 80 : & ,INDCEL,I1,IGR,IDIMP1 81 : C 82 : C**** Pour l'instant 2D 83 : C 84 : REAL*8 XG,YG,XD,YD,XC,YC,DXG,DYG,DXD,DYD,DG,DD,DT 85 : & ,PHI,PHIMAX,PHIMIN,DPHIST,GRADX,GRADY,DPHI0,DPHI1,ALPHA 86 : & ,VALEUR 87 : C 88 : C**** Extension 3D 89 : C 90 : REAL*8 ZG,ZD,ZC,DZG,DZD,GRADZ 91 : C 92 : C 93 : -INC SMCOORD 94 : -INC SMCHPOI 95 : -INC SMELEME 96 : -INC SMLENTI 97 : POINTEUR MPOMIN.MPOVAL, MPOMAX.MPOVAL, MPOALP.MPOVAL, 98 : & MPOCHP.MPOVAL, MPOGRA.MPOVAL 99 : POINTEUR MELEFL.MELEME, MLECEN.MLENTI 100 : C 101 : C**** N.B. Tous les pointeurs ici sont déjà activés! 102 : C 103 : C**** Maillage FACEL 104 : C Boucle sur les faces 105 : C 106 : C 107 : C**** NGCEG = Numero global centre d'elt "gauche" 108 : C NLCEG = Numero local centre d'elt "gauche" 109 : C NGCF = " global centre de face 110 : C NLCF = numero local centre de face 111 : C NGCEG = Numero global centre d'elt "gauche" 112 : C NLCEG = Numero local centre d'elt "gauche" 113 : C 114 : IDIMP1 = IDIM+1 115 : C 116 : C Cas de la dimension 2 117 : C 118 : IF (IDIM .EQ. 2) THEN 119 : DO NLCF = 1, NFAC 120 : NGCEG = MELEFL.NUM(1,NLCF) 121 : NLCEG = MLECEN.LECT(NGCEG) 122 : NGCF = MELEFL.NUM(2,NLCF) 123 : NGCED = MELEFL.NUM(3,NLCF) 124 : NLCED = MLECEN.LECT(NGCED) 125 : IF(NGCEG .EQ. NGCED)THEN 126 : C 127 : C********** Limiteur dans le cas mur 128 : C 129 : C 130 : C********** Coordonees et parametres geometriques 131 : C 132 : C 133 : C XCOOR : table de coordonnes de points; 134 : C pour le point de numero global NG 135 : C X_NG = XCOOR((NG-1)*(IDIM+1)+1) 136 : C Y_NG = XCOOR((NG-1)*(IDIM+1)+2) 137 : C Z_NG = XCOOR((NG-1)*(IDIM+1)+3) 138 : C (dans notre cas IDIM=2, i.e. 2D) 139 : C 140 : INDCEL = (NGCEG-1)*IDIMP1 141 : XG = XCOOR(INDCEL+1) 142 : YG = XCOOR(INDCEL+2) 143 : INDCEL = (NGCF-1)*IDIMP1 144 : XC = XCOOR(INDCEL+1) 145 : YC = XCOOR(INDCEL+2) 146 : DXG = XC - XG 147 : DYG = YC - YG 148 : C 149 : C********** Boucle sur le composantes 150 : C 151 : DO I1 = 1, NCOMP 152 : IGR = 2*I1 - 1 153 : C 154 : C************* Limiteur 155 : C 156 : PHI = MPOCHP.VPOCHA(NLCEG,I1) 157 : PHIMAX = MPOMAX.VPOCHA(NLCEG,I1) 158 : PHIMIN = MPOMIN.VPOCHA(NLCEG,I1) 159 : DPHIST = PHIMAX - PHIMIN 160 : GRADX = MPOGRA.VPOCHA(NLCEG,IGR) 161 : GRADY = MPOGRA.VPOCHA(NLCEG,IGR+1) 162 : DPHI0 = GRADX * DXG + GRADY * DYG 163 : ALPHA = MPOALP.VPOCHA(NLCEG,I1) 164 : IF(ABS(DPHI0) .LE. (DPHIST*1.0D-8))THEN 165 : VALEUR = 1.0D0 166 : ELSEIF(DPHI0 .GT. 0.0D0)THEN 167 : DPHI1 = PHIMAX - PHI 168 : VALEUR = DPHI1/DPHI0*0.5D0 169 : ELSE 170 : DPHI1 = PHIMIN - PHI 171 : VALEUR = DPHI1/DPHI0*0.5D0 172 : ENDIF 173 : ALPHA = MIN(ALPHA, VALEUR) 174 : MPOALP.VPOCHA(NLCEG,I1) = ALPHA 175 : ENDDO 176 : ELSE 177 : C 178 : C******* NGCEG != NGCED 179 : C 180 : INDCEL = (NGCEG-1)*IDIMP1 181 : XG = XCOOR(INDCEL+1) 182 : YG = XCOOR(INDCEL+2) 183 : INDCEL = (NGCED-1)*IDIMP1 184 : XD = XCOOR(INDCEL+1) 185 : YD = XCOOR(INDCEL+2) 186 : INDCEL = (NGCF-1)*IDIMP1 187 : XC = XCOOR(INDCEL+1) 188 : YC = XCOOR(INDCEL+2) 189 : C 190 : DXG = XC - XG 191 : DYG = YC - YG 192 : DXD = XC - XD 193 : DYD = YC - YD 194 : DG = DXG * DXG + DYG * DYG 195 : DG = SQRT(DG) 196 : DD = DXD * DXD + DYD * DYD 197 : DD = SQRT(DD) 198 : DT = DG + DD 199 : C 200 : C********** Boucle sur le composantes 201 : C 202 : DO I1 = 1, NCOMP 203 : IGR = 2*I1 - 1 204 : C 205 : C************* Limiteur a gauche 206 : C 207 : PHI = MPOCHP.VPOCHA(NLCEG,I1) 208 : PHIMAX = MPOMAX.VPOCHA(NLCEG,I1) 209 : PHIMIN = MPOMIN.VPOCHA(NLCEG,I1) 210 : DPHIST = PHIMAX - PHIMIN 211 : GRADX = MPOGRA.VPOCHA(NLCEG,IGR) 212 : GRADY = MPOGRA.VPOCHA(NLCEG,IGR+1) 213 : DPHI0 = GRADX * DXG + GRADY * DYG 214 : ALPHA = MPOALP.VPOCHA(NLCEG,I1) 215 : IF(ABS(DPHI0) .LE. (DPHIST*1.0D-8))THEN 216 : VALEUR = 1.0D0 217 : ELSEIF(DPHI0 .GT. 0.0D0)THEN 218 : DPHI1 = PHIMAX - PHI 219 : VALEUR = DPHI1/DPHI0*DG/DT 220 : ELSE 221 : DPHI1 = PHIMIN - PHI 222 : VALEUR = DPHI1/DPHI0*DG/DT 223 : ENDIF 224 : 225 : ALPHA = MIN(ALPHA, VALEUR) 226 : MPOALP.VPOCHA(NLCEG,I1) = ALPHA 227 : C 228 : C************* Limiteur a droite 229 : C 230 : PHI = MPOCHP.VPOCHA(NLCED,I1) 231 : PHIMAX = MPOMAX.VPOCHA(NLCED,I1) 232 : PHIMIN = MPOMIN.VPOCHA(NLCED,I1) 233 : DPHIST = PHIMAX - PHIMIN 234 : GRADX = MPOGRA.VPOCHA(NLCED,IGR) 235 : GRADY = MPOGRA.VPOCHA(NLCED,IGR+1) 236 : DPHI0 = GRADX * DXD + GRADY * DYD 237 : ALPHA = MPOALP.VPOCHA(NLCED,I1) 238 : IF(ABS(DPHI0) .LE. (DPHIST*1.0D-8))THEN 239 : VALEUR = 1.0D0 240 : ELSEIF(DPHI0 .GT. 0.0D0)THEN 241 : DPHI1 = PHIMAX - PHI 242 : VALEUR = DPHI1/DPHI0*DD/DT 243 : ELSE 244 : DPHI1 = PHIMIN - PHI 245 : VALEUR = DPHI1/DPHI0*DD/DT 246 : ENDIF 247 : ALPHA = MIN(ALPHA, VALEUR) 248 : MPOALP.VPOCHA(NLCED,I1) = ALPHA 249 : ENDDO 250 : ENDIF 251 : ENDDO 252 : C 253 : C Cas de la dimension 3 254 : C 255 : ELSE 256 : DO NLCF = 1, NFAC 257 : NGCEG = MELEFL.NUM(1,NLCF) 258 : NLCEG = MLECEN.LECT(NGCEG) 259 : NGCF = MELEFL.NUM(2,NLCF) 260 : NGCED = MELEFL.NUM(3,NLCF) 261 : NLCED = MLECEN.LECT(NGCED) 262 : IF(NGCEG .EQ. NGCED)THEN 263 : C 264 : C********** Limiteur dans le cas mur 265 : C 266 : C 267 : C********** Coordonees et parametres geometriques 268 : C 269 : C 270 : C XCOOR : table de coordonnes de points; 271 : C pour le point de numero global NG 272 : C X_NG = XCOOR((NG-1)*(IDIM+1)+1) 273 : C Y_NG = XCOOR((NG-1)*(IDIM+1)+2) 274 : C Z_NG = XCOOR((NG-1)*(IDIM+1)+3) 275 : C 276 : INDCEL = (NGCEG-1)*IDIMP1 277 : XG = XCOOR(INDCEL+1) 278 : YG = XCOOR(INDCEL+2) 279 : ZG = XCOOR(INDCEL+3) 280 : INDCEL = (NGCF-1)*IDIMP1 281 : XC = XCOOR(INDCEL+1) 282 : YC = XCOOR(INDCEL+2) 283 : ZC = XCOOR(INDCEL+3) 284 : DXG = XC - XG 285 : DYG = YC - YG 286 : DZG = ZC - ZG 287 : C 288 : C********** Boucle sur le composantes 289 : C 290 : DO I1 = 1, NCOMP 291 : IGR = IDIM*(I1-1)+1 292 : C 293 : C************* Limiteur 294 : C 295 : PHI = MPOCHP.VPOCHA(NLCEG,I1) 296 : PHIMAX = MPOMAX.VPOCHA(NLCEG,I1) 297 : PHIMIN = MPOMIN.VPOCHA(NLCEG,I1) 298 : DPHIST = PHIMAX - PHIMIN 299 : GRADX = MPOGRA.VPOCHA(NLCEG,IGR) 300 : GRADY = MPOGRA.VPOCHA(NLCEG,IGR+1) 301 : GRADZ = MPOGRA.VPOCHA(NLCEG,IGR+2) 302 : DPHI0 = GRADX * DXG + GRADY * DYG + GRADZ * DZG 303 : ALPHA = MPOALP.VPOCHA(NLCEG,I1) 304 : IF(ABS(DPHI0) .LE. (DPHIST*1.0D-8))THEN 305 : VALEUR = 1.0D0 306 : ELSEIF(DPHI0 .GT. 0.0D0)THEN 307 : DPHI1 = PHIMAX - PHI 308 : VALEUR = DPHI1/DPHI0*0.5D0 309 : ELSE 310 : DPHI1 = PHIMIN - PHI 311 : VALEUR = DPHI1/DPHI0*0.5D0 312 : ENDIF 313 : ALPHA = MIN(ALPHA, VALEUR) 314 : MPOALP.VPOCHA(NLCEG,I1) = ALPHA 315 : ENDDO 316 : ELSE 317 : C 318 : C******* NGCEG != NGCED 319 : C 320 : INDCEL = (NGCEG-1)*IDIMP1 321 : XG = XCOOR(INDCEL+1) 322 : YG = XCOOR(INDCEL+2) 323 : ZG = XCOOR(INDCEL+3) 324 : INDCEL = (NGCED-1)*IDIMP1 325 : XD = XCOOR(INDCEL+1) 326 : YD = XCOOR(INDCEL+2) 327 : ZD = XCOOR(INDCEL+3) 328 : INDCEL = (NGCF-1)*IDIMP1 329 : XC = XCOOR(INDCEL+1) 330 : YC = XCOOR(INDCEL+2) 331 : ZC = XCOOR(INDCEL+3) 332 : C 333 : DXG = XC - XG 334 : DYG = YC - YG 335 : DZG = ZC - ZG 336 : DXD = XC - XD 337 : DYD = YC - YD 338 : DZD = ZC - ZD 339 : DG = DXG * DXG + DYG * DYG + DZG * DZG 340 : DG = SQRT(DG) 341 : DD = DXD * DXD + DYD * DYD + DZD * DZD 342 : DD = SQRT(DD) 343 : DT = DG + DD 344 : C 345 : C********** Boucle sur le composantes 346 : C 347 : DO I1 = 1, NCOMP 348 : IGR = IDIM*(I1-1)+1 349 : C 350 : C************* Limiteur a gauche 351 : C 352 : PHI = MPOCHP.VPOCHA(NLCEG,I1) 353 : PHIMAX = MPOMAX.VPOCHA(NLCEG,I1) 354 : PHIMIN = MPOMIN.VPOCHA(NLCEG,I1) 355 : DPHIST = PHIMAX - PHIMIN 356 : GRADX = MPOGRA.VPOCHA(NLCEG,IGR) 357 : GRADY = MPOGRA.VPOCHA(NLCEG,IGR+1) 358 : GRADZ = MPOGRA.VPOCHA(NLCEG,IGR+2) 359 : DPHI0 = GRADX * DXG + GRADY * DYG + GRADZ * DZG 360 : ALPHA = MPOALP.VPOCHA(NLCEG,I1) 361 : IF(ABS(DPHI0) .LE. (DPHIST*1.0D-8))THEN 362 : VALEUR = 1.0D0 363 : ELSEIF(DPHI0 .GT. 0.0D0)THEN 364 : DPHI1 = PHIMAX - PHI 365 : VALEUR = DPHI1/DPHI0*DG/DT 366 : ELSE 367 : DPHI1 = PHIMIN - PHI 368 : VALEUR = DPHI1/DPHI0*DG/DT 369 : ENDIF 370 : 371 : ALPHA = MIN(ALPHA, VALEUR) 372 : MPOALP.VPOCHA(NLCEG,I1) = ALPHA 373 : C 374 : C************* Limiteur a droite 375 : C 376 : PHI = MPOCHP.VPOCHA(NLCED,I1) 377 : PHIMAX = MPOMAX.VPOCHA(NLCED,I1) 378 : PHIMIN = MPOMIN.VPOCHA(NLCED,I1) 379 : DPHIST = PHIMAX - PHIMIN 380 : GRADX = MPOGRA.VPOCHA(NLCED,IGR) 381 : GRADY = MPOGRA.VPOCHA(NLCED,IGR+1) 382 : GRADZ = MPOGRA.VPOCHA(NLCED,IGR+2) 383 : DPHI0 = GRADX * DXD + GRADY * DYD + GRADZ * DZD 384 : ALPHA = MPOALP.VPOCHA(NLCED,I1) 385 : IF(ABS(DPHI0) .LE. (DPHIST*1.0D-8))THEN 386 : VALEUR = 1.0D0 387 : ELSEIF(DPHI0 .GT. 0.0D0)THEN 388 : DPHI1 = PHIMAX - PHI 389 : VALEUR = DPHI1/DPHI0*DD/DT 390 : ELSE 391 : DPHI1 = PHIMIN - PHI 392 : VALEUR = DPHI1/DPHI0*DD/DT 393 : ENDIF 394 : ALPHA = MIN(ALPHA, VALEUR) 395 : MPOALP.VPOCHA(NLCED,I1) = ALPHA 396 : ENDDO 397 : ENDIF 398 : ENDDO 399 : ENDIF 400 : C 401 : RETURN 402 : END 403 : 404 : 405 : 406 : 407 : 408 : 409 : 410 : 411 : 412 : 413 : 414 : 415 : 416 :
© Cast3M 2003 - All rights reserved.
Disclaimer