1 : C PENTE2 SOURCE PV 09/03/12 21:29:58 6325 2 : SUBROUTINE PENTE2(IOP2,INORM,MLECEN,MLEFAC, 3 : & MLENCL,IMCHAM, 4 : & NCOMP,MPOCHP,MPOVCL,MPOGRA, 5 : & MPOMIN,MPOMAX) 6 : C 7 : C************************************************************************ 8 : C 9 : C PROJET : CASTEM 2000 10 : C 11 : C NOM : PENTE2 12 : C 13 : C DESCRIPTION : Cette subroutine est appellée par la subroutine 14 : C PENTE1 (calcul du gradient d'un CHPOINT de type 15 : C CENTRE) 16 : C Elle contient la partie du calcul du gradient. 17 : C 18 : C LANGUAGE : ESOPE 2000 avec extensions CISI 19 : C 20 : C AUTEUR : A. BECCANTINI, DRN/DMT/SEMT/TTMF 21 : C AUTEUR (modif.) : R. MOREL, DRN/DMT/SEMT/TTMF 22 : C 23 : C************************************************************************ 24 : C 25 : C 26 : C APPELES (Outils) : LICHT, ERREUR 27 : C 28 : C 29 : C************************************************************************ 30 : C 31 : C ENTREES : IOP2 : Entier indiquant le type de reconstruction (cf PENT) 32 : C 33 : C INORM : Pointeur sur le CHPOINT des normales aux faces 34 : C 35 : C MLECEN : table numerotation global/local pour les 36 : C points CENTRE (segment MLENTI) 37 : C 38 : C MLEFAC : table numerotation global/local pour les 39 : C points FACE 40 : C 41 : C MLENCL : table numerotation global/local pour les 42 : C SPG du CHPOINT de C.L. ou 43 : C table de 0 si ce CHPOINT n'existe pas 44 : C 45 : C IMCHAM : MCHAML de type 'GRADGEO' 46 : C 47 : C NCOMP : Nombre de composantes du CHPOINT dont 48 : C on veut le gradient 49 : C 50 : C MPOCHP : MPOVAL du CHPOINT dont on veut le gradient 51 : C 52 : C MPOVCL : MPOVAL du CHPOINT de C.L. 53 : C 54 : C SORTIES : MPOGRA : MPOVAL du gradient 55 : C 56 : C MPOMAX : MPOVAL du max sur le stencil 57 : C 58 : C MPOMIN : MPOVAL du min sur le stencil 59 : C 60 : C 61 : C************************************************************************ 62 : C 63 : C HISTORIQUE (Anomalies et modifications éventuelles) 64 : C 65 : C HISTORIQUE : Cree le 4-6-1998 66 : C 67 : C HISTORIQUE : Modifie le 20-10-1998 pour extension 3D 68 : C Modifie le 28-04-2000 pour reconstruction quadatique 69 : C exacte (A. Beccantini) 70 : C HISTORIQUE : Modifie le 05.07.2003: elimination reconstruction 71 : C quadratique exacte. 72 : C Algorithme de BArth-Jespersen remplacé par une 73 : C reconstruction lineaire exacte qui se base sur 74 : C LSM 75 : C************************************************************************ 76 : C 77 : IMPLICIT INTEGER(I-N) 78 : -INC CCOPTIO 79 : -INC SMCHAML 80 : -INC SMELEME 81 : -INC SMCHPOI 82 : -INC SMLENTI 83 : C 84 : POINTEUR INORM.MCHPOI 85 : POINTEUR MLENCL.MLENTI, MLECEN.MLENTI, MLEFAC.MLENTI 86 : POINTEUR MPOMAX.MPOVAL, MPOMIN.MPOVAL, MPONOR.MPOVAL, 87 : & MPOVCL.MPOVAL, MPOCHP.MPOVAL, MPOGRA.MPOVAL 88 : POINTEUR MELVX.MELVAL, MELVY.MELVAL, MELVZ.MELVAL, 89 : & MELVXX.MELVAL, MELVYY.MELVAL, MELVZZ.MELVAL, 90 : & MELVXY.MELVAL, MELVXZ.MELVAL, MELVYZ.MELVAL 91 : C 92 : C**** Variables de COOPTIO 93 : C 94 : INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR 95 : & ,IOTER, IOLEC, IOIMP, IOCAR, IOACQ 96 : & ,IOPER, IOSGB, IOGRA, IOSAU, IORES 97 : & ,IECHO, IIMPI, IOSPI 98 : & ,IDIM 99 : & ,MCOORD 100 : & ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE 101 : & ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU 102 : & ,NORINC,NORVAL,NORIND,NORVAD 103 : & ,NUCROU, IPSAUV, IFICLE, IPREFI, ISAFOR, IREFOR 104 : C 105 : INTEGER IOP2,NCOMP,NLELEM,NVOI,I3,IMCHAM 106 : & ,NGCE,NLCE,NTOTV,NGCV,NLCV,ICOM,IGEOM,K,NBELEM,NBNN 107 : & ,NLCF,NLCL,NMAI 108 : REAL*8 VAL,VCOMP(3),VN 109 : REAL*8 NX,NY,NZ 110 : CHARACTER*8 TYPE 111 : C 112 : C**** Segments déjà activés 113 : C 114 : C MPOVCL 115 : C MLENCL 116 : C MLECEN 117 : C MLEFAC 118 : C MPOCHP 119 : C MPOGRA*MOD 120 : C MPOMIN*MOD 121 : C MPOMAX*MOD 122 : C 123 : C**** N.B. MPOMAX.VPOCHA, MPOMIN.VPOCHA déjà initialisé: 124 : C 125 : C MPOMAX.VPOCHA = MPOCHP.VPOCHA 126 : C MPOMIN.VPOCHA = MPOCHP.VPOCHA 127 : C 127 : CALL LICHT(INORM,MPONOR,TYPE,IGEOM) 128 : C 129 : C**** En LICHT 130 : C SEGACT*MOD MPONOR 131 : C MPONOR: pointeur sur les normales aux faces 132 : C 133 : MCHELM = IMCHAM 134 : SEGACT MCHELM 135 : NMAI=MCHELM.ICHAML(/1) 136 : C 137 : C********************************************* 138 : C Premier cas : on traite un champ scalaire en 139 : C reconstruction non-quadratique 140 : C********************************************* 141 : C 142 : C 143 : IF(IOP2.EQ.1) THEN 144 : C 145 : C**** N.B: IOP2 = 1 -> EULESCAL 146 : C 147 : DO 10 K=1,NMAI 148 : C 149 : C********** Boucle sur les differents maillages elementaires 150 : C IPT1 = maillage duale 151 : C 152 : IPT1=MCHELM.IMACHE(K) 153 : SEGACT IPT1 154 : NBNN=IPT1.NUM(/1) 155 : NTOTV=NBNN 156 : C 157 : C********** NTOTV = Nombre de voisins 158 : C 159 : NBELEM = IPT1.NUM(/2) 160 : MCHAM1 = MCHELM.ICHAML(K) 161 : SEGACT MCHAM1 162 : MELVX = MCHAM1.IELVAL(1) 163 : MELVY = MCHAM1.IELVAL(2) 164 : SEGACT MELVX 165 : SEGACT MELVY 166 : IF (IDIM .EQ. 3) THEN 167 : MELVZ = MCHAM1.IELVAL(3) 168 : SEGACT MELVZ 169 : ENDIF 170 : SEGDES MCHAM1 171 : DO NLELEM = 1, NBELEM, 1 172 : C 173 : C************* Boucle sur les ELTs 174 : C NLELEM = Numero local d'ELT dans le maillage 1 175 : C 176 : C NGCE = numero global centre d'elt 177 : C NLCE = numero local centre d'elt 178 : C 179 : C MLECEN = table numerotation global/local pour les 180 : C points CENTRE (segment MLENTI) 181 : C 182 : C MLEFAC = table numerotation global/local pour les 183 : C points FACE 184 : C 185 : C MLENCL = table numerotation global/local pour les 186 : C SPG du CHPOINT de C.L. ou 187 : C table de 0 si ce CHPOINT n'existe pas 188 : C 189 : NGCE = IPT1.NUM(NBNN,NLELEM) 190 : NLCE = MLECEN.LECT(NGCE) 191 : DO NVOI = 1, NTOTV 192 : C 193 : C**************** Boucle sur les voisins de NGCE 194 : C 195 : C NGCV = Numero global voisin 196 : C NLCV = " local de NGCV dans le maillage CENTRE 197 : C NLCL = " local de NGCV dans le SPG du CHPOINT 198 : C des C.L. 199 : C 200 : NGCV = IPT1.NUM(NVOI,NLELEM) 201 : NLCV = MLECEN.LECT(NGCV) 202 : IF(NLCV.EQ.0)THEN 203 : NLCL = MLENCL.LECT(NGCV) 204 : C 205 : C**************** NGCV n'est pas un point centre, mais un point face, 206 : C et on est sur les bords 207 : C 208 : DO I3 = 1, NCOMP 209 : C 210 : C********************** I3 : numero de composante du CHPOI 211 : C ICOM : numero de composante du gradient 212 : C 213 : ICOM = IDIM*(I3 -1)+1 214 : IF (NLCL .EQ. 0) THEN 215 : C 216 : C************************* NGCV n'est pas un point du SPG de C.L. 217 : C 218 : VAL = MPOCHP.VPOCHA(NLCE,I3) 219 : ELSE 220 : C 221 : C************************* NGCV est un point du SPG de C.L. 222 : C 223 : VAL = MPOVCL.VPOCHA(NLCL,I3) 224 : ENDIF 225 : MPOGRA.VPOCHA(NLCE,ICOM)= MPOGRA.VPOCHA(NLCE 226 : $ ,ICOM)+VAL * MELVX.VELCHE(NVOI,NLELEM) 227 : MPOGRA.VPOCHA(NLCE,ICOM+1)=MPOGRA.VPOCHA(NLCE 228 : $ ,ICOM+1) +VAL * MELVY.VELCHE(NVOI,NLELEM) 229 : IF (IDIM .EQ. 3) THEN 230 : MPOGRA.VPOCHA(NLCE,ICOM+2)=MPOGRA.VPOCHA(NLCE 231 : $ ,ICOM+2) +VAL * MELVZ.VELCHE(NVOI 232 : $ ,NLELEM) 233 : ENDIF 234 : MPOMAX.VPOCHA(NLCE,I3)=MAX(MPOMAX.VPOCHA(NLCE,I3 235 : $ ),VAL) 236 : MPOMIN.VPOCHA(NLCE,I3)=MIN(MPOMIN.VPOCHA(NLCE,I3 237 : $ ),VAL) 238 : ENDDO 239 : ELSE 240 : C 241 : C********** NGCV est un point centre 242 : C 243 : DO I3 = 1, NCOMP 244 : ICOM = IDIM*(I3 -1)+1 245 : VAL = MPOCHP.VPOCHA(NLCV,I3) 246 : MPOGRA.VPOCHA(NLCE,ICOM)=MPOGRA.VPOCHA(NLCE,ICOM 247 : $ )+VAL * MELVX.VELCHE(NVOI,NLELEM) 248 : MPOGRA.VPOCHA(NLCE,ICOM+1)=MPOGRA.VPOCHA(NLCE 249 : $ ,ICOM+1)+VAL * MELVY.VELCHE(NVOI,NLELEM) 250 : IF (IDIM .EQ. 3) THEN 251 : MPOGRA.VPOCHA(NLCE,ICOM+2)=MPOGRA.VPOCHA(NLCE 252 : $ ,ICOM+2)+VAL * MELVZ.VELCHE(NVOI,NLELEM 253 : $ ) 254 : ENDIF 255 : MPOMAX.VPOCHA(NLCE,I3)=MAX(MPOMAX.VPOCHA(NLCE,I3 256 : $ ),VAL) 257 : MPOMIN.VPOCHA(NLCE,I3)=MIN(MPOMIN.VPOCHA(NLCE,I3 258 : $ ),VAL) 259 : ENDDO 260 : ENDIF 261 : ENDDO 262 : C 263 : C******* Fin boucle sur les voisins de NGCE 264 : C 265 : ENDDO 266 : SEGDES MELVX 267 : SEGDES MELVY 268 : IF (IDIM .EQ. 3) SEGDES MELVZ 269 : SEGDES IPT1 270 : 10 CONTINUE 271 : C 272 : C********************************************************** 273 : C Deuxieme cas : on traite un champ vectoriel (IOP2=2) 274 : C********************************************************** 275 : C 276 : C 277 : ELSEIF(IOP2 .EQ. 2)THEN 278 : DO 20 K=1,NMAI 279 : IPT1=MCHELM.IMACHE(K) 280 : SEGACT IPT1 281 : NBNN=IPT1.NUM(/1) 282 : NTOTV=NBNN 283 : NBELEM=IPT1.NUM(/2) 284 : MCHAM1=MCHELM.ICHAML(K) 285 : SEGACT MCHAM1 286 : MELVX=MCHAM1.IELVAL(1) 287 : MELVY=MCHAM1.IELVAL(2) 288 : SEGACT MELVX 289 : SEGACT MELVY 290 : IF (IDIM .EQ. 3) THEN 291 : MELVZ=MCHAM1.IELVAL(3) 292 : SEGACT MELVZ 293 : ENDIF 294 : SEGDES MCHAM1 295 : DO NLELEM = 1, NBELEM 296 : NGCE = IPT1.NUM(NBNN,NLELEM) 297 : NLCE = MLECEN.LECT(NGCE) 298 : DO NVOI = 1, NTOTV 299 : NGCV = IPT1.NUM(NVOI,NLELEM) 300 : NLCV = MLECEN.LECT(NGCV) 301 : IF(NLCV .EQ. 0)THEN 302 : NLCL = MLENCL.LECT(NGCV) 303 : C 304 : C************* NGCV n'est pas un point centre, mais un point face. 305 : C 306 : IF(NLCL .NE. 0)THEN 307 : C 308 : C************************* NGCV est un point du SPG de C.L. 309 : C 310 : DO I3 = 1, NCOMP 311 : VCOMP(I3) = MPOVCL.VPOCHA(NLCL,I3) 312 : ENDDO 313 : ELSE 314 : NLCF = MLEFAC.LECT(NGCV) 315 : IF (IDIM .EQ.2) THEN 316 : VCOMP(1) = MPOCHP.VPOCHA(NLCE,1) 317 : VCOMP(2) = MPOCHP.VPOCHA(NLCE,2) 318 : NX = MPONOR.VPOCHA(NLCF,1) 319 : NY = MPONOR.VPOCHA(NLCF,2) 320 : VN = VCOMP(1)*NX+VCOMP(2)*NY 321 : VCOMP(1) = VCOMP(1)-2*NX*VN 322 : VCOMP(2) = VCOMP(2)-2*NY*VN 323 : ELSE 324 : VCOMP(1) = MPOCHP.VPOCHA(NLCE,1) 325 : VCOMP(2) = MPOCHP.VPOCHA(NLCE,2) 326 : VCOMP(3) = MPOCHP.VPOCHA(NLCE,3) 327 : NX = MPONOR.VPOCHA(NLCF,1) 328 : NY = MPONOR.VPOCHA(NLCF,2) 329 : NZ = MPONOR.VPOCHA(NLCF,3) 330 : VN = VCOMP(1)*NX+VCOMP(2)*NY+VCOMP(3)*NZ 331 : VCOMP(1) = VCOMP(1)-2*NX*VN 332 : VCOMP(2) = VCOMP(2)-2*NY*VN 333 : VCOMP(3) = VCOMP(3)-2*NZ*VN 334 : ENDIF 335 : ENDIF 336 : DO I3 = 1, NCOMP 337 : ICOM = IDIM*(I3 -1)+1 338 : VAL = VCOMP(I3) 339 : MPOGRA.VPOCHA(NLCE,ICOM)= MPOGRA.VPOCHA(NLCE 340 : $ ,ICOM)+VAL * MELVX.VELCHE(NVOI 341 : $ ,NLELEM) 342 : MPOGRA.VPOCHA(NLCE,ICOM+1)=MPOGRA.VPOCHA(NLCE 343 : $ ,ICOM+1) +VAL * MELVY.VELCHE(NVOI 344 : $ ,NLELEM) 345 : IF (IDIM .EQ. 3) THEN 346 : MPOGRA.VPOCHA(NLCE,ICOM+2)=MPOGRA.VPOCHA(NLCE 347 : $ ,ICOM+2) +VAL * MELVZ.VELCHE(NVOI 348 : $ ,NLELEM) 349 : ENDIF 350 : MPOMAX.VPOCHA(NLCE,I3)=MAX(MPOMAX.VPOCHA(NLCE,I3 351 : $ ),VAL) 352 : MPOMIN.VPOCHA(NLCE,I3)=MIN(MPOMIN.VPOCHA(NLCE,I3 353 : $ ),VAL) 354 : ENDDO 355 : ELSE 356 : C 357 : C********** NGCV est un point centre 358 : C 359 : DO I3 = 1, NCOMP 360 : ICOM = IDIM*(I3 -1)+1 361 : VAL = MPOCHP.VPOCHA(NLCV,I3) 362 : MPOGRA.VPOCHA(NLCE,ICOM)=MPOGRA.VPOCHA(NLCE,ICOM 363 : $ )+VAL * MELVX.VELCHE(NVOI,NLELEM) 364 : MPOGRA.VPOCHA(NLCE,ICOM+1)=MPOGRA.VPOCHA(NLCE 365 : $ ,ICOM+1)+VAL * MELVY.VELCHE(NVOI,NLELEM) 366 : IF (IDIM .EQ. 3) THEN 367 : MPOGRA.VPOCHA(NLCE,ICOM+2)=MPOGRA.VPOCHA(NLCE 368 : $ ,ICOM+2)+VAL * MELVZ.VELCHE(NVOI,NLELEM 369 : $ ) 370 : ENDIF 371 : MPOMAX.VPOCHA(NLCE,I3)=MAX(MPOMAX.VPOCHA(NLCE,I3 372 : $ ),VAL) 373 : MPOMIN.VPOCHA(NLCE,I3)=MIN(MPOMIN.VPOCHA(NLCE,I3 374 : $ ),VAL) 375 : ENDDO 376 : ENDIF 377 : ENDDO 378 : C 379 : C******* Fin boucle sur les voisins de NGCE 380 : C 381 : ENDDO 382 : SEGDES MELVX 383 : SEGDES MELVY 384 : IF (IDIM .EQ. 3) SEGDES MELVZ 385 : SEGDES IPT1 386 : 20 CONTINUE 387 : ENDIF 388 : SEGDES MPONOR 389 : SEGDES MCHELM 390 : RETURN 391 : END 392 : 393 : 394 :
© Cast3M 2003 - All rights reserved.
Disclaimer