1 : C YTOIMP SOURCE MAGN 10/05/18 21:17:33 6675 2 : SUBROUTINE YTOIMP 3 : IMPLICIT INTEGER(I-N) 4 : IMPLICIT REAL*8 (A-H,O-Z) 5 : C*********************************************************************** 6 : C ----------------------------------------------------------- 7 : C --------- TOIMP ---------------------------------------- 8 : C ----------------------------------------------------------- 9 : C --------- PARAMETRI DELLO OPERATORE (NELLO ORDINE) : ----- 10 : C ----------------------------------------------------------- 11 : C --------- TENSION ( tau et pression ) --------- 12 : C ----------------------------------------------------------- 13 : C 14 : C SYNTAXE : 15 : C 16 : C TOIMP (tau pression) 17 : C 18 : C 1------2 19 : C (R1,AL1) LEF FLUIDE NOEUDS 1 2 20 : C 21 : C 22 : C 23 : C CAS TRIDIMENSIONNEL 24 : C 4 ________ 3 25 : C / FLUIDE / 26 : C 1 /________/2 27 : C 28 : C 29 : C*********************************************************************** 30 : 31 : -INC CCOPTIO 32 : -INC CCGEOME 33 : -INC CCREEL 34 : -INC SMMATRIK 35 : -INC SIZFFB 36 : PARAMETER (LRV=64) 37 : SEGMENT PETROV 38 : REAL*8 WT(LRV,NP,NPG,KDIM),WS(LRV,NP,NPG,KDIM),HK(LRV,IDIM,NP,NPG) 39 : REAL*8 UIL(LRV,IDIM,NPG),DUIL(LRV,IDIM,NPG) 40 : REAL*8 PGSK(LRV,NPG),RPGK(LRV,NPG),AIRE(LRV),ANUK(LRV) 41 : REAL*8 AJK(LRV,IDIM,IDIM,NPG) 42 : ENDSEGMENT 43 : 44 : -INC SMCOORD 45 : -INC SMLENTI 46 : POINTEUR IPADU.MLENTI 47 : -INC SMELEME 48 : POINTEUR MELEM1.MELEME 49 : -INC SMCHPOI 50 : POINTEUR IZG1.MCHPOI, IZGG1.MPOVAL 51 : POINTEUR IZTU1.MPOVAL,MZMU.MPOVAL,MZUN.MPOVAL 52 : POINTEUR MZTO.MPOVAL,MZDT.MPOVAL 53 : 54 : -INC SMLMOTS 55 : POINTEUR LINCO.MLMOTS 56 : CHARACTER*8 NOMI,NOMA,TYPE,TYPC,NOM0 57 : CHARACTER*4 NOM4(3) 58 : PARAMETER (NTB=1) 59 : CHARACTER*8 LTAB(NTB) 60 : DIMENSION KTAB(NTB),IXV(3) 61 : SAVE IPAS 62 : DATA LTAB/'KIZX '/,IPAS/0/ 63 : C***************************************************************************** 64 : CTOIMP 65 : C write(6,*)' Debut TOIMP ' 65 : CALL LITABS(LTAB,KTAB,NTB,1,IRET) 66 : IF (IERR.NE.0) RETURN 67 : MTABX=KTAB(1) 68 : C 69 : C- Récupération de la table EQEX (pointeur MTAB1) 70 : C 70 : CALL LEKTAB(MTABX,'EQEX',MTAB1) 71 : IF(MTAB1.EQ.0)THEN 72 : C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24 73 : MOTERR( 1: 8) = ' EQEX ' 74 : MOTERR( 9:16) = ' EQEX ' 75 : MOTERR(17:24) = ' KIZX ' 75 : CALL ERREUR(786) 76 : RETURN 77 : ENDIF 77 : CALL ACME(MTAB1,'NAVISTOK',NASTOK) 78 : IF(NASTOK.EQ.0)THEN 78 : CALL ZZOIMP(MTABX,MTAB1) 79 : RETURN 80 : ENDIF 81 : C 82 : C- Récupération de la table INCO (pointeur KINC) 83 : C 83 : CALL LEKTAB(MTAB1,'INCO',KINC) 84 : IF(KINC.EQ.0)THEN 85 : C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24 86 : MOTERR( 1: 8) = ' INCO ' 87 : MOTERR( 9:16) = ' INCO ' 88 : MOTERR(17:24) = ' EQEX ' 88 : CALL ERREUR(786) 89 : RETURN 90 : ENDIF 91 : 92 : C***************************************************************************** 93 : C OPTIONS 94 : C KIMPL = 0 -> EXPL 1 -> IMPL 2 -> SEMI 95 : C KFORM = 0 -> SI 1 -> EF 2 -> VF 3 -> EFMC 96 : C IDCEN = 0-> rien 1-> CENTREE 2-> SUPGDC 3-> SUPG 4-> TVISQUEU 5-> CNG 97 : C Pour IDCEN on autorise 1,2 et 3 98 : 99 : IAXI=0 100 : IF(IFOMOD.EQ.0)IAXI=2 101 : C 102 : C- Récupération de la table des options KOPT (pointeur KOPTI) 103 : C 103 : CALL LEKTAB(MTABX,'KOPT',KOPTI) 104 : IF (KOPTI.EQ.0) THEN 105 : C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24 106 : MOTERR( 1: 8) = ' KOPT ' 107 : MOTERR( 9:16) = ' KOPT ' 108 : MOTERR(17:24) = ' KIZX ' 108 : CALL ERREUR(786) 109 : RETURN 110 : ENDIF 111 : 112 : IKOMP=0 112 : CALL ACME(KOPTI,'KIMPL',KIMPL) 112 : CALL ACME(KOPTI,'KFORM',KFORM) 112 : CALL ACMF(KOPTI,'CMD',CMD) 112 : CALL ACME(KOPTI,'IDCEN',IDCEN) 113 : 114 : KDIM=1 115 : IF(IDCEN.EQ.2)KDIM=IDIM 116 : KDCEN=(IDCEN-1)*(IDCEN-2)*(IDCEN-3) 117 : IF(KDCEN.NE.0)THEN 118 : MOTERR( 1: 8) = 'IDCEN' 118 : CALL ERREUR(803) 119 : RETURN 120 : ENDIF 121 : 122 : IF(KFORM.NE.0.AND.KFORM.NE.1)THEN 123 : C Option %m1:8 incompatible avec les données 124 : MOTERR( 1: 8) = 'EF/EFM1 ' 124 : CALL ERREUR(803) 125 : RETURN 126 : ENDIF 127 : IF (IERR.NE.0) RETURN 128 : 129 : C write(6,*)' Apres les options ' 130 : C***************************************************************************** 131 : C 132 : C- Récupération de la table DOMAINE associée au domaine local 133 : C 133 : CALL LEKTAB(MTABX,'DOMZ',MTABZ) 134 : IF(MTABZ.EQ.0)THEN 135 : C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24 136 : MOTERR( 1: 8) = ' DOMZ ' 137 : MOTERR( 9:16) = ' DOMZ ' 138 : MOTERR(17:24) = ' KIZX ' 138 : CALL ERREUR(786) 139 : RETURN 140 : ENDIF 141 : 141 : CALL LEKTAB(MTABZ,'MAILLAGE',MELEME) 141 : CALL LEKTAB(MTABZ,'CENTRE',MELEMC) 142 : IF (IERR.NE.0) RETURN 143 : 144 : SEGACT MELEME 145 : 146 : C************************************************************************* 147 : C VERIFICATIONS SUR LES INCONNUES 148 : C 149 : C- Récupération du nombre d'inconnues et du nom de l'inconnue NOMI 150 : C 151 : TYPE='LISTMOTS' 151 : CALL ACMO(MTABX,'LISTINCO',TYPE,LINCO) 152 : IF (IERR.NE.0) RETURN 153 : SEGACT LINCO 154 : NBINC=LINCO.MOTS(/2) 155 : IF(NBINC.NE.1)THEN 156 : C Indice %m1:8 : contient plus de %i1 %m9:16 157 : MOTERR( 1:8) = 'LISTINCO' 158 : INTERR(1) = 1 159 : MOTERR(9:16) = ' MOTS ' 159 : CALL ERREUR(799) 160 : RETURN 161 : ENDIF 162 : 163 : C --> 1 ere Inconnue 164 : 165 : NOMI=LINCO.MOTS(1) 166 : DO 15 I=1,IDIM 167 : WRITE(NOM4(I),FMT='(I1,A3)')I,NOMI(1:3) 168 : 15 CONTINUE 169 : 170 : TYPE=' ' 170 : CALL ACMO(KINC,NOMI,TYPE,MCHPOI) 171 : IF(TYPE.NE.'CHPOINT ')THEN 172 : C Indice %m1:8 : ne contient pas un objet de type %m9:16 173 : MOTERR( 1: 8) = 'INC '//NOMI 174 : MOTERR( 9:16) = 'CHPOINT ' 174 : CALL ERREUR(800) 175 : RETURN 176 : ELSE 176 : CALL LICHT(MCHPOI,IZTU1,TYPC,MELEM1) 177 : ENDIF 178 : NPTU=IZTU1.VPOCHA(/2) 179 : 180 : C************************************************************************* 181 : C Le domaine de definition est donne par le SPG de la premiere inconnue 182 : C Les inconnues suivantes devront posseder ce meme pointeur 183 : C On verifie que les points de la zone sont tous inclus dans ce SPG 184 : 184 : CALL KRIPAD(MELEM1,MLENTI) 185 : IPADU=MLENTI 186 : 187 : IF(IPAS.EQ.0)THEN 187 : CALL VERPAD(MLENTI,MELEME,IRET) 188 : IF(IRET.NE.0)THEN 189 : C Indice %m1:8 : L'objet %m9:16 n'a pas le bon support géométrique 190 : MOTERR(1: 8) = 'INC '//NOMI 191 : MOTERR(9:16) = 'CHPOINT ' 191 : CALL ERREUR(788) 192 : IPAS=0 193 : RETURN 194 : ENDIF 195 : ENDIF 196 : 197 : C************************************************************************* 198 : C Lecture du coefficient 199 : C Type du coefficient : 200 : C IK1=0 CHPOINT IK1=1 scalaire IK1=2 vecteur 201 : 201 : CALL ACME(MTABX,'IARG',IARG) 202 : IF(IARG.NE.1.AND.IDCEN.EQ.1.AND. 203 : & IARG.NE.3.AND.(IDCEN.EQ.2.OR.IDCEN.EQ.2) 204 : &)THEN 205 : WRITE(6,*)'Opérateur TOIMP : nombre d''arguments incorrect' 206 : C Indice %m1:8 : nombre d'arguments incorrect 207 : MOTERR(1:8) = 'IARG ' 207 : CALL ERREUR(804) 208 : RETURN 209 : ENDIF 210 : 211 : MZMU=IZTU1 212 : MZUN=IZTU1 213 : MZDT=IZTU1 214 : 215 : IXV(1)=-MELEMC 216 : IXV(2)=0 217 : IXV(3)=1 217 : CALL LEKCOF('Opérateur TOIMP :', 218 : & MTABX,KINC,1,IXV,MTO,MZTO,NTAU,NC1,IKS,IRET) 219 : IF(IRET.EQ.0)RETURN 220 : 221 : IF(IKS.EQ.2)IKS=1 222 : 223 : 224 : IF(IARG.GE.2)THEN 225 : 226 : C 2ème coefficient UN , champ de vitesse transportant 227 : IXV(1)=-MELEM1 228 : IXV(2)=0 229 : IXV(3)=1 229 : CALL LEKCOF('Opérateur TOIM :', 230 : & MTABX,KINC,2,IXV,MUN,MZUN,NPTU,NC2,IKU,IRET) 231 : IF(IRET.EQ.0)GO TO 90 232 : IF(IKU.EQ.2)IKU=1 233 : 234 : C 3ème coefficient viscosité 235 : IXV(1)=0 236 : IXV(2)=1 237 : IXV(3)=0 237 : CALL LEKCOF('Opérateur TOIM :', 238 : & MTABX,KINC,3,IXV,MMU,MZMU,NPT3,NC3,IKM,IRET) 239 : IF(IRET.EQ.0)GO TO 90 240 : ELSE 240 : CALL LEKTAB(MTABX,'ARGS3',MMU) 240 : CALL LICHT(MMU,MZMU,TYPC,IGEOM0) 241 : MZMU.VPOCHA(1,1)=1.D0 242 : IKM=1 243 : 244 : ENDIF 245 : 246 : IF(IARG.EQ.4)THEN 247 : C 4ème coefficient Dt 248 : IXV(1)=0 249 : IXV(2)=1 250 : IXV(3)=0 250 : CALL LEKCOF('Opérateur TOIM :', 251 : & MTABX,KINC,4,IXV,MDT,MZDT,NPT4,NC4,IKT,IRET) 252 : IF(IRET.EQ.0)RETURN 253 : DT=MZDT.VPOCHA(1,1) 254 : ELSE 255 : DT=0.D0 256 : 257 : ENDIF 258 : 259 : C write(6,*)' Fin lecture Arguments ' 260 : C Fin lecture Arguments ************************************************ 261 : 262 : 263 : IF(KIMPL.EQ.0)THEN 264 : NC=IZTU1.VPOCHA(/2) 265 : TYPE='SOMMET' 265 : CALL KRCHPT(TYPE,MELEM1,NC,IZG1,NOM4) 266 : 267 : ELSE 268 : 269 : NAT=2 270 : NSOUPO=1 271 : SEGACT MELEM1 272 : N=MELEM1.NUM(/2) 273 : NC=IZTU1.VPOCHA(/2) 274 : NINKO=NC 275 : SEGINI MCHPO1,MSOUP1,MPOVA1 276 : MCHPO1.IFOPOI=IFOMOD 277 : MCHPO1.MOCHDE=TITREE 278 : MCHPO1.MTYPOI='SMBR' 279 : MCHPO1.JATTRI(1)=2 280 : MCHPO1.IPCHP(1)=MSOUP1 281 : DO 177 N=1,NINKO 282 : WRITE(NOM0,FMT='(I1,A7)')N,NOMI(1:7) 283 : MSOUP1.NOCOMP(N)=NOM0(1:4) 284 : 177 CONTINUE 285 : MSOUP1.IGEOC=MELEM1 286 : MSOUP1.IPOVAL=MPOVA1 287 : IZG1=MCHPO1 288 : 289 : ENDIF 290 : 290 : CALL LICHT(IZG1,IZGG1,TYPC,IGEOM) 291 : 292 : IF(IGEOM.NE.MELEM1)THEN 293 : WRITE(6,*)' Opérateur TOIM' 294 : WRITE(6,*)' Incompatibilité de SPG entre 1ères inconnues' 295 : RETURN 296 : ENDIF 297 : 298 : SEGACT MELEME 299 : NBSOUS=LISOUS(/1) 300 : IF(NBSOUS.EQ.0)NBSOUS=1 301 : 302 : NPTD=IZGG1.VPOCHA(/1) 303 : K0=0 304 : DO 1 LS=1,NBSOUS 305 : IPT1=MELEME 306 : IF(NBSOUS.NE.1)IPT1=LISOUS(LS) 307 : SEGACT IPT1 308 : 309 : NOM0 = NOMS(IPT1.ITYPEL)//' ' 309 : CALL KALPBG(NOM0,'FONFORM ',IZFFM) 310 : SEGACT IZFFM 311 : IZHR=KZHR(1) 312 : SEGACT IZHR*MOD 313 : NES=GR(/1) 314 : NPG=GR(/3) 315 : 316 : NP =IPT1.NUM(/1) 317 : NBEL=IPT1.NUM(/2) 318 : 319 : SEGINI PETROV 320 : 321 : DEUPI=1.D0 322 : IF(IAXI.NE.0)DEUPI=2.D0*XPI 323 : 324 : C Calcul du nombre de paquets de LRV éléments 325 : C 326 : NNN=MOD(NBEL,LRV) 327 : IF(NNN.EQ.0) THEN 328 : NPACK=NBEL/LRV 329 : ELSE 330 : NPACK=1+(NBEL-NNN)/LRV 331 : ENDIF 332 : KPACKD=1 333 : KPACKF=NPACK 334 : C 335 : C ******* BOUCLE SUR LES PAQUETS DE LRV ELEMENTS ********** 336 : C 337 : DO 7001 KPACK=KPACKD,KPACKF 338 : C 339 : C ======= A L'INTERIEUR DE CHAQUE PAQUET DE LRV ELEMENTS ======= 340 : C 341 : C 1. Calcul des limites du paquet courant. 342 : KDEB=1+(KPACK-1)*LRV 343 : KFIN=MIN(NBEL,KDEB+LRV-1) 344 : C 345 : DO 7002 K=KDEB,KFIN 346 : KP=K-KDEB+1 347 : NK=K0+K 348 : NKM=(1-IKM)*(NK-1)+1 349 : ANUK(KP)=MZMU.VPOCHA(NKM,1) 350 : 7002 CONTINUE 351 : 352 : IF(IDCEN.EQ.2)THEN 353 : 354 : DO 108 NC=1,IDIM 354 : CALL SUPGEF(FN,GR,PG,XYZ,HR,PGSQ,RPG,AJ, 355 : & NES,NP,NPG,IAXI,XCOOR, 356 : & IPT1.NUM,KDEB,KFIN,LRV,IDCEN,CMD,IKOMP, 357 : & MZUN.VPOCHA(1,NC),IPADU.LECT,MZUN.VPOCHA,IPADU.LECT,NPTU,ANUK, 358 : & WT(1,1,1,NC),WS(1,1,1,NC),HK,PGSK,RPGK,AJK,AIRE,UIL,DUIL, 359 : & DTM1,DT,DTT1,DTT2,DIAEL,NUEL) 360 : 108 CONTINUE 361 : 362 : ELSE 363 : 363 : CALL SUPGEF(FN,GR,PG,XYZ,HR,PGSQ,RPG,AJ, 364 : & NES,NP,NPG,IAXI,XCOOR, 365 : & IPT1.NUM,KDEB,KFIN,LRV,IDCEN,CMD,IKOMP, 366 : & MZUN.VPOCHA,IPADU.LECT,MZUN.VPOCHA,IPADU.LECT,NPTU,ANUK, 367 : & WT,WS,HK,PGSK,RPGK,AJK,AIRE,UIL,DUIL, 368 : & DTM1,DT,DTT1,DTT2,DIAEL,NUEL) 369 : 370 : ENDIF 371 : 372 : 373 : IF(NES.NE.IDIM)THEN 374 : DO 7003 K=KDEB,KFIN 375 : 376 : KP=K-KDEB+1 377 : NK=K0+K 378 : KA=1+(1-IKS)*(NK-1) 379 : 380 : DO 50 I=1,NP 381 : DO 50 N=1,IDIM 382 : N1=1 383 : IF(IDCEN.EQ.2)N1=N 384 : FF1=0.D0 385 : DO 51 L=1,NPG 386 : DO 52 M=1,IDIM 387 : FF1=FF1 + MZTO.VPOCHA(KA,M)*AJK(KP,N,M,L)*WT(KP,I,L,N1)*PGSK(KP,L) 388 : &*DEUPI*RPGK(KP,L) 389 : 52 CONTINUE 390 : 51 CONTINUE 391 : 392 : NF=LECT(IPT1.NUM(I,K)) 393 : IZGG1.VPOCHA(NF,N)=IZGG1.VPOCHA(NF,N)+FF1 394 : 395 : 50 CONTINUE 396 : 7003 CONTINUE 397 : 398 : ELSEIF(NES.EQ.IDIM)THEN 399 : DO 7004 K=KDEB,KFIN 400 : 401 : KP=K-KDEB+1 402 : NK=K0+K 403 : KA=1+(1-IKS)*(NK-1) 404 : 405 : DO 60 I=1,NP 406 : DO 60 N=1,IDIM 407 : N1=1 408 : IF(IDCEN.EQ.2)N1=N 409 : FF1=0.D0 410 : DO 61 L=1,NPG 411 : FF1=FF1 + MZTO.VPOCHA(KA,N)*WT(KP,I,L,N1)*PGSK(KP,L) 412 : &*DEUPI*RPGK(KP,L) 413 : 61 CONTINUE 414 : 415 : NF=LECT(IPT1.NUM(I,K)) 416 : IZGG1.VPOCHA(NF,N)=IZGG1.VPOCHA(NF,N)+FF1 417 : 418 : 60 CONTINUE 419 : 7004 CONTINUE 420 : ENDIF 421 : 422 : 7001 CONTINUE 423 : 424 : SEGDES MZTO,IPT1 425 : K0=K0+NBEL 426 : 427 : SEGSUP PETROV 428 : 1 CONTINUE 429 : SEGDES MELEME 430 : 431 : NRIGE=7 432 : NKID =9 433 : NKMT =7 434 : NMATRI=0 435 : SEGINI MATRIK 435 : CALL ECROBJ('MATRIK',MATRIK) 435 : CALL ECROBJ('CHPOINT',IZG1) 436 : 437 : 438 : SEGDES MELEME 439 : 440 : SEGDES IZTU1 441 : SEGDES IZG1,IZGG1 442 : SEGDES LINCO 443 : SEGSUP MLENTI 444 : IPAS=1 445 : C write(6,*)' FIN TOIMP ' 446 : RETURN 447 : 90 CONTINUE 448 : WRITE(6,*)' Interuption anormale de TOIMP ' 449 : C Option %m1:8 incompatible avec les données 450 : MOTERR( 1: 8) = ' EF ' 450 : CALL ERREUR(803) 451 : RETURN 452 : 1002 FORMAT(10(1X,1PE11.4)) 453 : END 454 : 455 : 456 : 457 : 458 : 459 : 460 : 461 : 462 :
© Cast3M 2003 - All rights reserved.
Disclaimer