1 : C YDUDW SOURCE MAGN 10/05/18 21:16:57 6675 2 : SUBROUTINE YDUDW 3 : IMPLICIT INTEGER(I-N) 4 : IMPLICIT REAL*8 (A-H,O-Z) 5 : C*********************************************************************** 6 : C 7 : C CET OPERATEUR DISCRETISE L'OPERATEUR DE PENALISATION DIV(U)=EPS*P 8 : C EN 2D SUR LES ELEMENTS QUA4 TRI3 TRI7 et QUA9 PLAN OU AXI 9 : C EN 3D SUR LES ELEMENTS CUB8 ET PRI6 10 : C 11 : C SYNTAXE : 12 : C --------- 13 : C 14 : C DUDW(EPS) INCO UN : 15 : C 16 : C COEFFICIENTS : 17 : C -------------- 18 : C 19 : C 20 : C EPS (SCAL DOMA) PARAMETRE DE PENALISATION 21 : C (SCAL ELEM) 22 : C 23 : C INCONNUES : 24 : C ----------- 25 : C 26 : C UN CHAMP DE VITESSE 27 : C 28 : C 29 : C OPTIONS : POROSITE DIV(PU)=EPS*P 30 : C SOURCE DIV(U)-Q=EPS*P 31 : C 32 : C************************************************************************ 33 : 34 : -INC CCOPTIO 35 : -INC CCGEOME 36 : -INC SMCOORD 37 : -INC SMLENTI 38 : -INC SMELEME 39 : POINTEUR MELEM1.MELEME,MELEMS.MELEME,MELEMI.MELEME,MELEMP.MELEME 40 : POINTEUR MELEMC.MELEME 41 : -INC SMCHPOI 42 : POINTEUR IZGG1.MPOVAL,SRCE.MPOVAL,IZTGG1.MPOVAL 43 : 44 : -INC SIZFFB 45 : POINTEUR IZF1.IZFFM 46 : -INC SMMATRIK 47 : 48 : -INC SMLMOTS 49 : POINTEUR LINCO.MLMOTS 50 : SEGMENT TRAV 51 : REAL*8 AF(NP,NP,NINC,NINC),AS(NP,NINC),CT(MP1,NP,NINC),PQ(MP1) 52 : ENDSEGMENT 53 : CHARACTER*8 NOMZ,NOMI,TYPE,TYPC 54 : CHARACTER*8 NOM0,NOM1,NOM2,NOM3 55 : DIMENSION IXV(3) 56 : PARAMETER (NTB=1) 57 : CHARACTER*8 LTAB(NTB) 58 : DIMENSION KTAB(NTB) 59 : DATA LTAB/'KIZX '/ 60 : C***************************************************************************** 61 : CDUDW 62 : C write(6,*)' Operateur DUDW ' 63 : 63 : CALL LITABS(LTAB,KTAB,NTB,1,IRET) 64 : IF(IERR.NE.0)RETURN 65 : MTABX=KTAB(1) 66 : C 67 : C- Récupération de la table EQEX (pointeur MTAB1) 68 : C 68 : CALL LEKTAB(MTABX,'EQEX',MTAB1) 69 : IF(MTAB1.EQ.0)THEN 70 : C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24 71 : MOTERR( 1: 8) = ' EQEX ' 72 : MOTERR( 9:16) = ' EQEX ' 73 : MOTERR(17:24) = ' KIZX ' 73 : CALL ERREUR(786) 74 : RETURN 75 : ENDIF 75 : CALL ACME(MTAB1,'NAVISTOK',NASTOK) 76 : IF(NASTOK.EQ.0)THEN 76 : CALL ZDUDW(MTABX,MTAB1) 77 : RETURN 78 : ENDIF 79 : C 80 : C- Récupération de la table INCO (pointeur KINC) 81 : C 81 : CALL LEKTAB(MTAB1,'INCO',KINC) 82 : IF(KINC.EQ.0)THEN 83 : C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24 84 : MOTERR( 1: 8) = ' INCO ' 85 : MOTERR( 9:16) = ' INCO ' 86 : MOTERR(17:24) = ' EQEX ' 86 : CALL ERREUR(786) 87 : RETURN 88 : ENDIF 89 : 90 : C***************************************************************************** 91 : C OPTIONS 92 : C KIMPL = 0 -> EXPL 1 -> IMPL 2 -> SEMI 93 : C KFORM = 0 -> SI 1 -> EF 2 -> VF 3 -> EFMC 94 : C IDCEN = 0-> rien 1-> CENTREE 2-> SUPGDC 3-> SUPG 4-> TVISQUEU 5-> CNG 95 : 96 : IAXI=0 97 : IF(IFOMOD.EQ.0)IAXI=2 98 : C 99 : C- Récupération de la table des options KOPT (pointeur KOPTI) 100 : C 100 : CALL LEKTAB(MTABX,'KOPT',KOPTI) 101 : IF (KOPTI.EQ.0) THEN 102 : C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24 103 : MOTERR( 1: 8) = ' KOPT ' 104 : MOTERR( 9:16) = ' KOPT ' 105 : MOTERR(17:24) = ' KIZX ' 105 : CALL ERREUR(786) 106 : RETURN 107 : ENDIF 108 : 108 : CALL ACME(KOPTI,'KIMPL',KIMPL) 108 : CALL ACME(KOPTI,'KPOIN',KPRE) 108 : CALL ACME(KOPTI,'KFORM',KFORM) 109 : 110 : IF(KFORM.NE.0.AND.KFORM.NE.1)THEN 111 : C Option %m1:8 incompatible avec les données 112 : MOTERR( 1: 8) = 'EF/EFM1 ' 112 : CALL ERREUR(803) 113 : RETURN 114 : ENDIF 115 : 116 : C write(6,*)' Apres les options ' 117 : C***************************************************************************** 118 : C 119 : C- Récupération de la table DOMAINE associée au domaine local 120 : C 120 : CALL ACMM(MTABX,'NOMZONE',NOMZ) 120 : CALL LEKTAB(MTABX,'DOMZ',MTABZ) 121 : IF(MTABZ.EQ.0)THEN 122 : C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24 123 : MOTERR( 1: 8) = ' DOMZ ' 124 : MOTERR( 9:16) = ' DOMZ ' 125 : MOTERR(17:24) = ' KIZX ' 125 : CALL ERREUR(786) 126 : RETURN 127 : ENDIF 128 : 128 : CALL LEKTAB(MTABZ,'SOMMET',MELEMS) 128 : CALL LEKTAB(MTABZ,'MAILLAGE',MELEME) 128 : CALL LEKTAB(MTABZ,'MACRO1',MACRO1) 129 : MELEMI=MACRO1 129 : CALL LEKTAB(MTABZ,'QUADRATI',MQUAD) 130 : IF (IERR.NE.0) RETURN 131 : 132 : IF(MQUAD.EQ.0.AND.MACRO1.EQ.0.AND.KPRE.NE.2)THEN 133 : WRITE(6,*)'Operateur DUDW ' 134 : WRITE(6,*)'Type d''éléments non prévu' 135 : RETURN 136 : ENDIF 137 : 138 : IF(KPRE.EQ.2)THEN 138 : CALL LEKTAB(MTABZ,'CENTRE',MELEMC) 139 : MELEMP=MELEMC 140 : MELEMI=MELEME 141 : ELSEIF(KPRE.EQ.3)THEN 141 : CALL LEKTAB(MTABZ,'CENTREP0',MELEMC) 142 : MELEMP=MELEMC 143 : ELSEIF(KPRE.EQ.4)THEN 143 : CALL LEKTAB(MTABZ,'CENTREP1',MELEMC) 143 : CALL LEKTAB(MTABZ,'ELTP1NC ',MELEMP) 144 : ENDIF 145 : 146 : SEGACT MELEME 147 : SEGACT MELEMP 148 : C************************************************************************* 149 : C VERIFICATIONS SUR LES INCONNUES 150 : C write(6,*)' Verification des inconnues ' 151 : 152 : TYPE='LISTMOTS' 152 : CALL ACMO(MTABX,'LISTINCO',TYPE,LINCO) 153 : SEGACT LINCO 154 : 155 : NBINC=LINCO.MOTS(/2) 156 : IF(NBINC.NE.1)THEN 157 : C Indice %m1:8 : contient plus de %i1 %m9:16 158 : MOTERR( 1:8) = 'LISTINCO' 159 : INTERR(1) = 1 160 : MOTERR(9:16) = ' MOTS ' 160 : CALL ERREUR(799) 161 : RETURN 162 : ENDIF 163 : 164 : NOMI=LINCO.MOTS(1) 165 : 166 : TYPE=' ' 166 : CALL ACMO(KINC,NOMI,TYPE,MCHPOI) 167 : IF(TYPE.NE.'CHPOINT ')THEN 168 : C Indice %m1:8 : ne contient pas un objet de type %m9:16 169 : MOTERR( 1: 8) = 'INC '//NOMI 170 : MOTERR( 9:16) = 'CHPOINT ' 170 : CALL ERREUR(800) 171 : RETURN 172 : ELSE 172 : CALL LICHT(MCHPOI,IZTU1,TYPC,MELEM1) 173 : ENDIF 174 : 175 : C************************************************************************* 176 : C Le domaine de definition est donne par le SPG de la premiere inconnue 177 : C Les inconnues suivantes devront posseder ce meme pointeur 178 : C On verifie que les points de la zone sont tous inclus dans ce SPG 179 : 179 : CALL KRIPAD(MELEM1,MLENTI) 180 : 180 : CALL VERPAD(MLENTI,MELEME,IRET) 181 : IF(IRET.NE.0)THEN 182 : C Indice %m1:8 : L'objet %m9:16 n'a pas le bon support géométrique 183 : MOTERR(1: 8) = 'INC '//NOMI 184 : MOTERR(9:16) = 'CHPOINT ' 184 : CALL ERREUR(788) 185 : RETURN 186 : ENDIF 187 : 188 : NINC=IDIM 189 : 190 : MLENT1=MLENTI 191 : 192 : C***************************************************************************** 193 : C Lecture du ou des coefficients 194 : C write(6,*)' Lecture des coefficients' 195 : 195 : CALL ACME(MTABX,'IARG',IARG) 196 : 197 : C? IXV(1)=MELEMC 198 : IXV(1)=0 199 : IXV(2)=1 200 : IXV(3)=0 200 : CALL LEKCOF('Opérateur DUDW :', 201 : & MTABX,KINC,1,IXV,IZTG1,IZTGG1,NPT1,NC1,IK1,IRET) 202 : IF(IRET.EQ.0)RETURN 203 : KITT=2 204 : 205 : INDGS=0 206 : SRCE=IZTGG1 207 : MPOVA1 =IZTGG1 208 : MPOVA3 =IZTGG1 209 : 210 : IF(IARG.EQ.2)THEN 211 : IXV(1)=MELEMC 212 : IXV(2)=1 213 : IXV(3)=0 213 : CALL LEKCOF('Opérateur DUDW :', 214 : & MTABX,KINC,2,IXV,IZS,SRCE,NPTS,NCS,IKS,IRET) 215 : IF(IRET.EQ.0)RETURN 216 : INDGS=1 216 : CALL KRIPAD(MELEMC,MLENT1) 217 : ENDIF 218 : 219 : C************************************************************************* 220 : IF(KFORM.EQ.0)THEN 221 : C CAS FORMULATION EF SI (GRESHO) 222 : 223 : IF(KIMPL.NE.0)THEN 224 : GO TO 90 225 : ENDIF 226 : 227 : WRITE(6,*)' Operateur DUDW ' 228 : WRITE(6,*)' Cas Formulation EF SI ' 229 : WRITE(6,*)' Cas invalide ' 230 : GO TO 90 231 : 232 : C************************************************************************* 233 : ELSE IF(KFORM.EQ.1)THEN 234 : C CAS FORMULATION EF 235 : 236 : NUTOEL=0 237 : 238 : SEGACT MELEMI 239 : NBSOUS=MELEMI.LISOUS(/1) 240 : IF(NBSOUS.EQ.0)NBSOUS=1 241 : 242 : NRIGE=7 243 : NKID =9 244 : NKMT =7 245 : NMATRI=1 246 : SEGINI MATRIK 247 : IRIGEL(1,1)=MELEMI 248 : IRIGEL(2,1)=MELEMI 249 : NBOP=0 250 : NBME=NINC*NINC 251 : C NBME=2 252 : NBELC=0 253 : SEGINI IMATRI 254 : IRIGEL(4,1)=IMATRI 255 : KSPGP=MELEMS 256 : KSPGD=MELEMS 257 : 258 : IF(NBME.EQ.4)THEN 259 : NOM1='1'//NOMI(1:7) 260 : NOM2='2'//NOMI(1:7) 261 : LISPRI(1)=NOM1 262 : LISDUA(1)=NOM1 263 : LISPRI(2)=NOM2 264 : LISDUA(2)=NOM1 265 : LISPRI(3)=NOM1 266 : LISDUA(3)=NOM2 267 : LISPRI(4)=NOM2 268 : LISDUA(4)=NOM2 269 : ELSEIF(NBME.EQ.9)THEN 270 : NOM1='1'//NOMI(1:7) 271 : NOM2='2'//NOMI(1:7) 272 : NOM3='3'//NOMI(1:7) 273 : LISPRI(1)=NOM1 274 : LISDUA(1)=NOM1 275 : LISPRI(2)=NOM2 276 : LISDUA(2)=NOM1 277 : LISPRI(3)=NOM3 278 : LISDUA(3)=NOM1 279 : 280 : LISPRI(4)=NOM1 281 : LISDUA(4)=NOM2 282 : LISPRI(5)=NOM2 283 : LISDUA(5)=NOM2 284 : LISPRI(6)=NOM3 285 : LISDUA(6)=NOM2 286 : 287 : LISPRI(7)=NOM1 288 : LISDUA(7)=NOM3 289 : LISPRI(8)=NOM2 290 : LISDUA(8)=NOM3 291 : LISPRI(9)=NOM3 292 : LISDUA(9)=NOM3 293 : 294 : ELSEIF(NBME.EQ.2)THEN 295 : 296 : LISPRI(1)='1'//NOMI(1:7) 297 : LISDUA(1)='1'//NOMI(1:7) 298 : LISPRI(2)='2'//NOMI(1:7) 299 : LISDUA(2)='2'//NOMI(1:7) 300 : 301 : ELSE 302 : WRITE(6,*)' Operateur DUDW ' 303 : WRITE(6,*)' Cas invalide' 304 : GO TO 90 305 : ENDIF 306 : 307 : 308 : IF(INDGS.NE.0)THEN 309 : NAT=2 310 : NSOUPO=1 311 : SEGACT MELEMS 312 : N=MELEMS.NUM(/2) 313 : NC=NINC 314 : SEGINI MCHPO1,MSOUP1,MPOVA1 315 : MCHPO1.IFOPOI=IFOMOD 316 : MCHPO1.MOCHDE=TITREE 317 : MCHPO1.MTYPOI='SMBR' 318 : MCHPO1.JATTRI(1)=2 319 : MCHPO1.IPCHP(1)=MSOUP1 320 : MSOUP1.NOCOMP(1)=NOM1 321 : MSOUP1.NOCOMP(2)=NOM2 322 : IF(NINC.EQ.3)MSOUP1.NOCOMP(3)=NOM3 323 : MSOUP1.IGEOC=MELEMS 324 : MSOUP1.IPOVAL=MPOVA1 325 : 326 : SEGACT MELEMC 327 : N=MELEMC.NUM(/2) 328 : NC=1 329 : SEGINI MCHPO3,MSOUP3,MPOVA3 330 : MCHPO3.IFOPOI=IFOMOD 331 : MCHPO3.MOCHDE=TITREE 332 : MCHPO3.MTYPOI='SMBR' 333 : MCHPO3.JATTRI(1)=2 334 : MCHPO3.IPCHP(1)=MSOUP3 335 : MSOUP3.NOCOMP(1)='SCAL' 336 : MSOUP3.IGEOC=MELEMC 337 : MSOUP3.IPOVAL=MPOVA3 338 : 339 : ENDIF 340 : 341 : DO 101 L=1,NBSOUS 342 : IPT1=MELEMI 343 : IPT2=MELEMI 344 : IF(NBSOUS.NE.1)IPT1=MELEMI.LISOUS(L) 345 : SEGACT IPT1 346 : IF(INDGS.NE.0)THEN 347 : IPT2=MELEMP 348 : IF(NBSOUS.NE.1)IPT2=MELEMP.LISOUS(L) 349 : SEGACT IPT2 350 : ENDIF 351 : 352 : IF(MQUAD.NE.0)THEN 353 : IF(KPRE.EQ.2)NOM0=NOMS(IPT1.ITYPEL)//'PRP0' 354 : IF(KPRE.EQ.3)NOM0=NOMS(IPT1.ITYPEL)//'PRP0' 355 : IF(KPRE.EQ.4)NOM0=NOMS(IPT1.ITYPEL)//'PRP1' 356 : ELSEIF(MACRO1.NE.0)THEN 357 : IF(KPRE.EQ.2)NOM0=NOMS(IPT1.ITYPEL)//' ' 358 : IF(KPRE.EQ.3)NOM0=NOMS(IPT1.ITYPEL)//'MCP0' 359 : IF(KPRE.EQ.4)NOM0=NOMS(IPT1.ITYPEL)//'MCP1' 360 : ELSE 361 : IF(KPRE.EQ.2)NOM0=NOMS(IPT1.ITYPEL)//' ' 362 : ENDIF 363 : 363 : CALL KALPBG(NOM0,'FONFORM ',IZFFM) 364 : 365 : 366 : SEGACT IZFFM*MOD 367 : IZHR=KZHR(1) 368 : SEGACT IZHR*MOD 369 : NES=GR(/1) 370 : NPG=GR(/3) 371 : IZF1=KTP(1) 372 : SEGACT IZF1*MOD 373 : MP1=IZF1.FN(/1) 374 : NP = IPT1.NUM(/1) 375 : 376 : MP = NP 377 : NBEL=IPT1.NUM(/2) 378 : SEGINI IPM1,IPM2,IPM3,IPM4 379 : LIZAFM(L,1)=IPM1 380 : C LIZAFM(L,2)=IPM4 381 : LIZAFM(L,2)=IPM2 382 : LIZAFM(L,3)=IPM3 383 : LIZAFM(L,4)=IPM4 384 : IPM5=IPM1 385 : IPM6=IPM1 386 : IPM7=IPM1 387 : IPM8=IPM1 388 : IPM9=IPM1 389 : IF(NBME.EQ.9)THEN 390 : SEGINI IPM5,IPM6,IPM7,IPM8,IPM9 391 : LIZAFM(L,5)=IPM5 392 : LIZAFM(L,6)=IPM6 393 : LIZAFM(L,7)=IPM7 394 : LIZAFM(L,8)=IPM8 395 : LIZAFM(L,9)=IPM9 396 : ENDIF 397 : 398 : C Pour l'instant les sources et puits de masse ne sont pas actifs 399 : INDG2=0 400 : KI2=0 401 : KJ2=0 402 : 403 : SEGINI TRAV 404 : 405 : NPT=MPOVA1.VPOCHA(/1) 405 : CALL XDUDW(FN,IZF1.FN,GR,PG,XYZ,HR,PGSQ,RPG,AJ, 406 : & NES,IDIM,NP,MP1,NPG,IAXI,NINC, 407 : & IZTGG1.VPOCHA,IK1,SRCE.VPOCHA,INDGS,IKS, 408 : & IPT1.NUM,NBEL,NUTOEL,XCOOR,AF,AS,CT,PQ, 409 : & IPM1.AM,IPM2.AM,IPM3.AM,IPM4.AM,IPM5.AM,IPM6.AM,IPM7.AM, 410 : & IPM8.AM,IPM9.AM,MPOVA1.VPOCHA,NPT,LECT,IPT2.NUM,MLENT1.LECT, 411 : & MPOVA3.VPOCHA) 412 : 413 : SEGSUP TRAV 414 : SEGDES IPT1 415 : NUTOEL=NUTOEL+NBEL 416 : 101 CONTINUE 417 : SEGDES IMATRI 418 : SEGDES MELEME,MATRIK 419 : SEGDES IZTGG1 420 : 420 : CALL ECROBJ('MATRIK',MATRIK) 420 : C? CALL ECMO(MTABX,'MATELM','MATRIK',MATRIK) 421 : 422 : NAT=2 423 : NSOUPO=0 424 : SEGINI MCHPOI 425 : JATTRI(1)=2 425 : CALL ECROBJ('CHPOINT',MCHPOI) 426 : 427 : 428 : 429 : IF(INDGS.NE.0)THEN 430 : 431 : TYPE=' ' 431 : CALL ACMO(MTAB1,'SMBR',TYPE,MCHPO2) 432 : IF(TYPE.NE.'CHPOINT')THEN 432 : CALL ECMO(MTAB1,'SMBR','CHPOINT',MCHPO1) 433 : ELSE 433 : CALL ECROBJ('CHPOINT',MCHPO2) 433 : CALL ECROBJ('CHPOINT',MCHPO1) 433 : CALL PRFUSE 433 : CALL LIROBJ('CHPOINT',MCHPOI,1,IRET) 433 : CALL ECMO(MTAB1,'SMBR','CHPOINT',MCHPOI) 434 : ENDIF 435 : 436 : TYPE=' ' 436 : C CALL ACMO(MTAB1,'DUDWSRCE',TYPE,MCHPO4) 437 : C IF(TYPE.NE.'CHPOINT')THEN 437 : CALL ECMO(MTAB1,'DUDWSRCE','CHPOINT',MCHPO3) 438 : C ELSE 438 : C CALL ECROBJ('CHPOINT',MCHPO3) 438 : C CALL ECROBJ('CHPOINT',MCHPO4) 438 : C CALL PRFUSE 438 : C CALL LIROBJ('CHPOINT',MCHPOI,1,IRET) 438 : C CALL ECMO(MTAB1,'DUDWSRCE','CHPOINT',MCHPOI) 439 : C ENDIF 440 : 441 : ENDIF 442 : 443 : SEGSUP MLENTI 444 : IF(INDGS.NE.0)SEGSUP MLENT1 445 : 446 : RETURN 447 : C************************************************************************* 448 : ELSE IF(KFORM.EQ.2)THEN 449 : C CAS FORMULATION VF 450 : 451 : RETURN 452 : ENDIF 453 : C************************************************************************* 454 : 1002 FORMAT(10(1X,1PE11.4)) 455 : 90 CONTINUE 456 : WRITE(6,*)' Interuption anormale de DUDW' 457 : RETURN 458 : END 459 : 460 : 461 : 462 : 463 : 464 : 465 : 466 : 467 :
© Cast3M 2003 - All rights reserved.
Disclaimer