1 : C YGMV SOURCE CHAT 05/01/13 04:17:43 5004 2 : SUBROUTINE YGMV 3 : IMPLICIT INTEGER(I-N) 4 : IMPLICIT REAL*8 (A-H,O-Z) 5 : C*********************************************************************** 6 : C 7 : C CE SP DISCRETISE LE TERME DE SOURCE DE QDM DANS LES EQUATIONS DE 8 : C NAVIER STOKES 9 : C EN 2D SUR LES ELEMENTS QUA4 ET TRI3 PLAN OU AXI 10 : C EN 3D SUR LES ELEMENTS CUB8 ET PRI6 11 : C L OPERATEUR EST "SOUS-INTEGRE" 12 : C 13 : C SYNTAXE : 14 : C --------- 15 : C 16 : C GMV TABGMV INCO UN : 17 : C 18 : C COEFFICIENT : 19 : C ------------- 20 : C 21 : C TABGMV Table contenant les entrees suivantes 22 : C 23 : C 24 : C INCONNUES : 25 : C ----------- 26 : C 27 : C UN CHAMPS DE VITESSE 28 : C 29 : C 30 : C*********************************************************************** 31 : 32 : -INC CCOPTIO 33 : -INC SMLREEL 34 : POINTEUR MQR.MLREEL,MPR.MLREEL 35 : -INC SMEVOLL 36 : 37 : -INC SMCOORD 38 : -INC SMLENTI 39 : POINTEUR IZIPAD.MLENTI 40 : -INC SMELEME 41 : POINTEUR MELEM1.MELEME,MELEMC.MELEME,MDEBI.MELEME 42 : POINTEUR MENTR.MELEME,MSORT.MELEME 43 : -INC SMCHPOI 44 : POINTEUR IZG1.MCHPOI,IZGG1.MPOVAL 45 : POINTEUR IZTU1.MPOVAL,IZPP.MPOVAL 46 : POINTEUR IZVOL.MPOVAL 47 : 48 : 49 : -INC SMTABLE 50 : POINTEUR KIZG.TABLE,MTABX.MTABLE,OPTI.MTABLE 51 : POINTEUR INCO.TABLE,KOPT.MTABLE 52 : POINTEUR MTABZ.MTABLE,MTABD.MTABLE,MTABP.MTABLE 53 : POINTEUR MTABA.MTABLE 54 : -INC SMLMOTS 55 : POINTEUR LINCO.MLMOTS 56 : CHARACTER*8 NOMZ,NOMI,NOMA,TYPE,TYPC 57 : REAL*8 XVEC1(3),XVEC2(3) 58 : PARAMETER (NTB=1) 59 : CHARACTER*8 LTAB(NTB) 60 : DIMENSION KTAB(NTB) 61 : SAVE IPAS 62 : DATA LTAB/'KIZX '/ 63 : DATA IPAS/0/ 64 : C***************************************************************************** 65 : 66 : CGMV 66 : CALL LITABS(LTAB,KTAB,NTB,1,IRET) 67 : IF(IRET.EQ.0)THEN 68 : WRITE(6,*)' Operateur GMV ' 69 : WRITE(6,*)' On attend un ensemble de table soustypes' 70 : RETURN 71 : ENDIF 72 : MTABX=KTAB(1) 73 : SEGACT MTABX 74 : 75 : C***************************************************************************** 76 : C OPTIONS 77 : C CES PARAMETRES SONT INITIALISES POUR ETRE EN DECENTRE ET NE PAS AVOIR DE 78 : C POROSITE : 79 : 80 : IOP4=0 81 : IAXI=0 82 : IF(IFOMOD.EQ.0)IAXI=2 82 : CALL LEKTAB(MTABX,'KOPT',OPTI) 83 : IF(OPTI.NE.0)THEN 84 : SEGACT OPTI 85 : TYPE=' ' 85 : CALL ACMO(OPTI,'POROSITE',TYPE,MCHPOI) 86 : IF(TYPE.EQ.'CHPOINT ')THEN 87 : IOP7=1 87 : CALL LICHT(MCHPOI,IZPORO,TYPC,IGEOM) 88 : ENDIF 89 : ENDIF 90 : C***************************************************************************** 91 : 91 : CALL ACMM(MTABX,'NOMZONE',NOMZ) 91 : CALL LEKTAB(MTABX,'DOMZ',MTABZ) 92 : IF(MTABZ.EQ.0)THEN 93 : WRITE(6,*)' Operateur GMV ' 94 : WRITE(6,*)' On ne trouve pas l''indice DOMZ ? ' 95 : GO TO 90 96 : ENDIF 97 : SEGACT MTABZ 98 : 98 : CALL LEKTAB(MTABZ,'MAILLAGE',MELEME) 99 : IF(MELEME.EQ.0)GO TO 90 100 : SEGACT MELEME 101 : 102 : 102 : CALL LEKTAB(MTABZ,'XXVOLUM',MCHPOI) 103 : IF(MCHPOI.EQ.0)GO TO 90 103 : CALL LICHT(MCHPOI,IZVOL,TYPC,IGEOM) 104 : 105 : C*** 106 : 107 : TYPE='LISTMOTS' 107 : CALL ACMO(MTABX,'LISTINCO',TYPE,LINCO) 108 : SEGACT LINCO 109 : 109 : CALL ACME(MTABX,'IARG',IARG) 110 : 111 : IF(IARG.NE.1)THEN 112 : WRITE(6,*)' Operateur GMV ' 113 : WRITE(6,*)' Nombre d''arguments ( ',IARG,' ) incorrect ' 114 : RETURN 115 : ENDIF 116 : 117 : TYPE=' ' 117 : CALL ACMO(MTABX,'ARG1',TYPE,MTABA) 118 : IF(TYPE.NE.'TABLE')THEN 119 : WRITE(6,*)' Operateur GMV ' 120 : WRITE(6,*) 121 : &' LE TYPE DE L''ARGUMENT (',TYPE,') N EST PAS CONVENABLE' 122 : WRITE(6,*)' On attend une table ' 123 : RETURN 124 : ELSE 125 : SEGACT MTABA 126 : 127 : TYPE=' ' 127 : CALL ACMO(MTABA,'DIR',TYPE,IPQ) 128 : 129 : IF(TYPE.NE.'POINT')THEN 130 : WRITE(6,*)' Operateur GMV ' 131 : WRITE(6,*)' Entree DIR erronee ' 132 : RETURN 133 : ELSE 134 : XVEC1(1)=XCOOR((IPQ-1)*(IDIM+1) +1) 135 : XVEC1(2)=XCOOR((IPQ-1)*(IDIM+1) +2) 136 : XNN=XVEC1(1)*XVEC1(1)+XVEC1(2)*XVEC1(2) 137 : IF(IDIM.EQ.3)THEN 138 : XVEC1(3)=XCOOR((IPQ-1)*(IDIM+1) +3) 139 : XNN=XNN+XVEC1(3)*XVEC1(3) 140 : ENDIF 141 : XVEC2(1)=XVEC1(1)/XNN 142 : XVEC2(2)=XVEC1(2)/XNN 143 : IF(IDIM.EQ.3)XVEC2(3)=XVEC1(3)/XNN 144 : 145 : ENDIF 146 : 147 : TYPE=' ' 147 : CALL ACMO(MTABA,'IMPR',TYPE,IRET) 148 : IF(TYPE.NE.'ENTIER')THEN 149 : IMPR=0 150 : ELSE 150 : CALL ACME(MTABA,'IMPR',IMPR) 151 : ENDIF 152 : 153 : TYPE=' ' 153 : CALL ACMO(MTABA,'KIMP',TYPE,IRET) 154 : IF(TYPE.NE.'FLOTTANT')THEN 155 : IKIMP=0 156 : ELSE 157 : IKIMP=1 157 : CALL ACMF(MTABA,'KIMP',AKIMP) 158 : ENDIF 159 : 160 : IF(IKIMP.EQ.0)THEN 161 : 162 : TYPE=' ' 162 : CALL ACMO(MTABA,'K0',TYPE,IRET) 163 : IF(TYPE.NE.'FLOTTANT')THEN 164 : AK0=1. 164 : CALL ECMF(MTABA,'K0',AK0) 165 : ELSE 165 : CALL ACMF(MTABA,'K0',AK0) 166 : ENDIF 167 : 168 : TYPE=' ' 168 : CALL ACMO(MTABA,'OMEGA',TYPE,IRET) 169 : IF(TYPE.NE.'FLOTTANT')THEN 170 : W1=0.1 170 : CALL ECMF(MTABA,'OMEGA',W1) 171 : ELSE 171 : CALL ACMF(MTABA,'OMEGA',W1) 172 : ENDIF 173 : 174 : TYPE=' ' 174 : CALL ACMO(MTABA,'GMV',TYPE,MEVOLL) 175 : IF(TYPE.EQ.'EVOLUTIO')THEN 176 : SEGACT MEVOLL 177 : KEVOLL=IEVOLL(1) 178 : SEGACT KEVOLL 179 : MQR=IPROGX 180 : MPR=IPROGY 181 : SEGDES KEVOLL,MEVOLL 182 : ELSE 183 : WRITE(6,*)' Operateur GMV ' 184 : WRITE(6,*)' Entree GMV erronee ' 185 : WRITE(6,*)' On attend un type EVOLUTION' 186 : RETURN 187 : ENDIF 188 : 189 : ENDIF 190 : 191 : 192 : TYPE=' ' 192 : CALL ACMO(MTABA,'LDEBIT',TYPE,MDEBI) 193 : IF(TYPE.NE.'MAILLAGE')THEN 194 : WRITE(6,*)' Operateur GMV ' 195 : WRITE(6,*)' Entree LDEBIT erronee ' 196 : WRITE(6,*)' On attend un type MAILLAGE' 197 : RETURN 198 : ENDIF 199 : 200 : TYPE=' ' 200 : CALL ACMO(MTABA,'PENTREE',TYPE,MENTR) 201 : IF(TYPE.NE.'MAILLAGE')THEN 202 : WRITE(6,*)' Operateur GMV ' 203 : WRITE(6,*)' Entree PENTREE erronee ' 204 : WRITE(6,*)' On attend un type MAILLAGE' 205 : RETURN 206 : ENDIF 207 : 208 : TYPE=' ' 208 : CALL ACMO(MTABA,'PSORTIE',TYPE,MSORT) 209 : IF(TYPE.NE.'MAILLAGE')THEN 210 : WRITE(6,*)' Operateur GMV ' 211 : WRITE(6,*)' Entree PSORTIE erronee ' 212 : WRITE(6,*)' On attend un type MAILLAGE' 213 : RETURN 214 : ENDIF 215 : 216 : ENDIF 217 : C****************** Fin lecture Table Arguments ***************** 218 : 218 : CALL LEKTAB(MTABX,'EQEX',MTAB1) 219 : IF(MTAB1.EQ.0)THEN 220 : WRITE(6,*)' Operateur GMV ' 221 : WRITE(6,*)' On ne trouve pas l''indice EQEX ? ' 222 : GO TO 90 223 : ENDIF 224 : SEGACT MTAB1 224 : CALL LEKTAB(MTAB1,'DOMAINE',MTABD) 225 : IF(MTABD.EQ.0)THEN 226 : WRITE(6,*)' Operateur GMV ' 227 : WRITE(6,*)' On ne trouve pas l''indice DOMAINE ?' 228 : GO TO 90 229 : ENDIF 229 : CALL LEKTAB(MTABD,'SOMMET',MELEM1) 230 : IF(MELEM1.EQ.0)THEN 231 : WRITE(6,*)' Operateur GMV ' 232 : WRITE(6,*)' On ne trouve pas l''indice SOMMET ?' 233 : GO TO 90 234 : ENDIF 235 : 235 : CALL LEKTAB(MTABD,'CENTRE',MELEMC) 236 : IF(MELEMC.EQ.0)THEN 237 : WRITE(6,*)' Operateur GMV ' 238 : WRITE(6,*)' On ne trouve pas l''indice CENTRE ?' 239 : GO TO 90 240 : ENDIF 241 : 241 : CALL LEKTAB(MTAB1,'INCO',INCO) 242 : IF(INCO.EQ.0)THEN 243 : WRITE(6,*)' Operateur GMV ' 244 : WRITE(6,*)'Il n''y a pas de table INCO ' 245 : RETURN 246 : ENDIF 247 : SEGACT INCO 248 : 248 : CALL KRIPAD(MELEM1,IZIPAD) 249 : SEGACT MELEMC 250 : NBC=MELEMC.NUM(/2) 250 : CALL RSETXI(IZIPAD.LECT,MELEMC.NUM,NBC) 251 : SEGDES MELEMC 252 : C***************************************************************************** 253 : 253 : CALL LEKTAB(MTAB1,'KIZG',KIZG) 254 : IF(KIZG.EQ.0)THEN 254 : CALL CRTABL(KIZG) 254 : CALL ECMM(KIZG,'SOUSTYPE','KIZG') 254 : CALL ECMO(MTAB1,'KIZG','TABLE ',KIZG) 255 : ELSE 256 : SEGACT KIZG 257 : ENDIF 258 : 259 : C VERIFICATIONS SUR LES INCONNUES 260 : NBINC=LINCO.MOTS(/2) 261 : IF(NBINC.NE.1)THEN 262 : WRITE(6,*)' Operateur GMV ' 263 : WRITE(6,*)'Nombre d''inconnues incorrecte : ',NBINC,' On attend 1' 264 : RETURN 265 : ENDIF 266 : 267 : NOMI=LINCO.MOTS(1) 268 : SEGACT INCO 269 : 270 : TYPE=' ' 270 : CALL ACMO(INCO,NOMI,TYPE,MCHPOI) 271 : IF(TYPE.NE.'CHPOINT ')THEN 272 : WRITE(6,*)' Operateur GMV ' 273 : WRITE(6,*)' L objet CHPOINT ',NOMI,' n existe pas dans la table' 274 : RETURN 275 : ELSE 275 : CALL LICHT(MCHPOI,IZTU1,TYPC,IGEOM0) 276 : ENDIF 277 : 278 : TYPE=' ' 278 : CALL ACMO(KIZG,NOMI,TYPE,IZG1) 279 : IF(TYPE.NE.'CHPOINT ')THEN 280 : NC=IZTU1.VPOCHA(/2) 281 : TYPE='SOMMET' 281 : CALL CRCHPT(TYPE,IGEOM0,NC,IZG1) 281 : CALL ECMO(KIZG,NOMI,'CHPOINT ',IZG1) 282 : ENDIF 283 : 283 : CALL LICHT(IZG1,IZGG1,TYPC,IGEOM) 284 : 285 : TYPE=' ' 285 : CALL ACMO(MTAB1,'PRESSION',TYPE,MTABP) 286 : IF(TYPE.NE.'TABLE')THEN 287 : WRITE(6,*)' Operateur GMV ' 288 : WRITE(6,*)' Il n''y a pas de table pression' 289 : RETURN 290 : ELSE 291 : SEGACT MTABP 292 : TYPE=' ' 292 : CALL ACMO(MTABP,'PRESSION',TYPE,MCHP) 293 : IF(TYPE.NE.'CHPOINT')THEN 294 : WRITE(6,*)' Operateur GMV ' 295 : WRITE(6,*)' Il n''y a pas de Champ de pression' 296 : RETURN 297 : ELSE 297 : CALL LICHT(MCHP,IZPP,TYPE,IGEOMP) 298 : IF(IZPP.EQ.0)GO TO 90 299 : ENDIF 300 : 301 : ENDIF 302 : 303 : SEGACT MELEME,MENTR 304 : NE=MENTR.NUM(/2) 305 : PE=0. 306 : DO 31 I=1,NE 307 : I1=MENTR.NUM(1,I) 308 : NPP1=IZIPAD.LECT(I1) 309 : PE=PE+IZPP.VPOCHA(NPP1,1) 310 : 31 CONTINUE 311 : PE=PE/FLOAT(NE) 312 : 313 : SEGACT MSORT 314 : NS=MSORT.NUM(/2) 315 : PS=0. 316 : DO 32 I=1,NS 317 : I1=MSORT.NUM(1,I) 318 : NPP1=IZIPAD.LECT(I1) 319 : PS=PS+IZPP.VPOCHA(NPP1,1) 320 : 32 CONTINUE 321 : PS=PS/FLOAT(NS) 322 : DELTAP=PS-PE 323 : 324 : SEGACT MDEBI 325 : NNP=MDEBI.NUM(/1) 326 : NNE=MDEBI.NUM(/2) 327 : DO 33 K=1,NNE 328 : DO 33 I=1,NNP 329 : I1=MDEBI.NUM(I,K) 330 : IF(IZIPAD.LECT(I1).EQ.0)THEN 331 : WRITE(6,*)' Operateur GMV ' 332 : WRITE(6,*)' LDEBIT SPG incompatible ' 333 : RETURN 334 : ENDIF 335 : 33 CONTINUE 336 : 337 : IMPRD=0 337 : CALL FFDBIT(IZTU1,MDEBI,IZIPAD,IAXI,Q,IMPRD) 338 : 339 : 340 : IF(IKIMP.EQ.0)THEN 341 : SEGACT MQR,MPR 342 : NPC=MPR.PROG(/1) 343 : 343 : CALL INTEPL(DELTAP,MQR.PROG,MPR.PROG,NPC,QTH,IRET) 344 : 344 : CALL ACMF(MTABA,'K0',AK0) 345 : IF(QTH.LT.Q)THEN 346 : AK1=AK0/(1.D0+(Q-QTH)/(Q+1.E-20)) 347 : ELSE 348 : IF(DELTAP.GE.0)THEN 349 : AK1=AK0*QTH/(Q+1.E-20) 350 : ELSE 351 : AK1=AK0*Q/(QTH+1.E-20) 352 : ENDIF 353 : ENDIF 354 : AK0=W1*AK1+(1.-W1)*AK0 355 : 356 : ELSE 357 : AK0=AKIMP 358 : ENDIF 359 : 360 : IF(IMPR.NE.0)THEN 361 : IMP=MOD(IPAS,IMPR) 362 : IF(IMP.EQ.0)THEN 363 : 364 : IF(IKIMP.NE.0)THEN 365 : WRITE(6,*)' GMV : CAS K IMPOSE ',AKIMP 366 : QTH=0. 367 : ENDIF 368 : WRITE(6,1888) DELTAP,Q,QTH,AK0 369 : 1888 FORMAT(1X,' GMV : DELTAP=',1PE12.4,' QR=',1PE11.4,' QTh=',1PE11.4, 370 : &' AK0=',1PE11.4) 371 : ENDIF 372 : ENDIF 373 : 374 : XVEC2(1)=XVEC2(1)*AK0 375 : XVEC2(2)=XVEC2(2)*AK0 376 : XVEC2(3)=XVEC2(3)*AK0 376 : CALL ECMF(MTABA,'K0',AK0) 377 : 378 : NPT=IZGG1.VPOCHA(/1) 379 : 380 : SEGACT MELEME 381 : NBSOUS=LISOUS(/1) 382 : IF(NBSOUS.NE.0)THEN 383 : WRITE(6,*)' Operateur GMV ' 384 : WRITE(6,*) 385 : & ' La zone doit etre compose d''un seul type d''element' 386 : RETURN 387 : ENDIF 388 : 389 : SEGACT IZTU1 390 : NBSOUS=1 391 : NUTOEL=0 392 : 393 : DO 1 L=1,NBSOUS 394 : IPT1=MELEME 395 : IF(NBSOUS.NE.1)IPT1=LISOUS(L) 396 : SEGACT IPT1 397 : 398 : NP =IPT1.NUM(/1) 399 : NBEL=IPT1.NUM(/2) 400 : 401 : IES=IDIM 402 : C 402 : CALL ZGMV(IPT1.NUM,NBEL,NUTOEL,NPT,IES,NP,IAXI, 403 : & IZIPAD.LECT,XVEC2, 404 : & IZTU1.VPOCHA, 405 : & IZGG1.VPOCHA, 406 : & IZVOL.VPOCHA) 407 : 408 : SEGDES IPT1 409 : NUTOEL=NUTOEL+NBEL 410 : 411 : 1 CONTINUE 412 : SEGDES MELEME 413 : 414 : IF(IKIMP.EQ.0)SEGDES MPR,MQR 415 : SEGDES IZTU1 416 : SEGDES IZGG1 417 : SEGDES IZVOL 418 : SEGDES LINCO 419 : SEGDES MTABX,MTAB1,INCO,KIZG,MTABA,MTABP,MTABD,MTABZ 420 : SEGDES IZPP,MDEBI,MSORT,MENTR 421 : SEGSUP IZIPAD 422 : 423 : 89 CONTINUE 424 : IPAS=IPAS+1 425 : RETURN 426 : 90 CONTINUE 427 : WRITE(6,*)' Interuption anormale de GMV' 428 : RETURN 429 : 1002 FORMAT(10(1X,1PE11.4)) 430 : END 431 : 432 : 433 :
© Cast3M 2003 - All rights reserved.
Disclaimer