1 : C YNSKE SOURCE PV 13/04/16 21:15:28 7765 2 : SUBROUTINE YNSKE 3 : IMPLICIT INTEGER(I-N) 4 : IMPLICIT REAL*8 (A-H,O-Z) 5 : C*********************************************************************** 6 : C 7 : C CET OPERATEUR DISCRETISE LES EQUATIONS DE NAVIER STOKES COUPLEES 8 : C AUX DEUX EQUATIONS DU MODELE K - EPSILON 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 LES OPERATEURS SONT "SOUS-INTEGRES" 12 : C 13 : C SYNTAXE : 14 : C --------- 15 : C 1/ Cas incompréssible 16 : C 17 : C du/dt + u Grad u = (nu+nut) Lapl u - 1/ro Grad p < + S > 18 : C < + g beta (T-Tref) > 19 : C 20 : C 'OPER' 'NSKE' nu nut 'INCO' UN : 21 : C 'OPER' 'NSKE' nu nut s 'INCO' UN : 22 : C 'OPER' 'NSKE' nu nut gb tn tref 'INCO' UN : 23 : C 24 : C 25 : C 2/ Cas compréssible 26 : C 27 : C dG/dt + u Grad G + G Div u = (mu+mut)(Lapl u + 1/3 Grad Div u) 28 : C - Grad p < + S > 29 : C 30 : C 'OPER' 'NSKE' ro mu mut un 'INCO' GN : 31 : C 'OPER' 'NSKE' ro mu mut un S 'INCO' GN : 32 : C 33 : C ro densité 34 : C FLOTTANT ou CHPOINT SCAL CENTRE 35 : C nu,mu viscosité cinématique (resp. dynamique) moléculaire 36 : C FLOTTANT ou CHPOINT SCAL CENTRE 37 : C nut,mut viscosité cinématique (resp. dynamique) turbulente 38 : C CHPOINT SCAL CENTRE 39 : C s source volumique de qdm 40 : C POINT ou CHPOINT VECT CENTRE 41 : C gb coéfficient de flottabilité (g*beta où g est l'accéllération 42 : C de la pesanteur et beta le coéfficient de dilatabilité) 43 : C POINT ou CHPOINT VECT CENTRE 44 : C tn Champ de température CHPOINT SCAL SOMMET 45 : C tref température de référence 46 : C FLOTTANT ou CHPOINT SCAL SOMMET 47 : C 48 : C Champ de vitesse -> VITESS 49 : C un Champ de vitesse transportant -> UTRANS 50 : C CHPOINT VECT SOMMET 51 : C gn Champ de vitesse massique (transporté) -> IZTU1 (Inconnue) 52 : C CHPOINT VECT SOMMET 53 : C kn Energie cinétique turbulente 54 : C CHPOINT SCAL SOMMET 55 : C en Taux de disparition de k 56 : C CHPOINT SCAL SOMMET 57 : C Constantes du modèle K - Epsilon 58 : C cnu = 0.09 c1 = 1.41 c2 = 1.92 59 : C 60 : C************************************************************************ 61 : 62 : -INC CCOPTIO 63 : -INC CCGEOME 64 : -INC SIZFFB 65 : POINTEUR IZF1.IZFFM 66 : SEGMENT PETROV 67 : REAL*8 WT(LRV,NP,NPG,KDIM),WS(LRV,NP,NPG,KDIM),HK(LRV,IDIM,NP,NPG) 68 : REAL*8 UIL(LRV,IDIM,NPG),DUIL(LRV,IDIM,NPG) 69 : REAL*8 PGSK(LRV,NPG),RPGK(LRV,NPG),AIRE(LRV),ANUK(LRV),COEFK(LRV) 70 : REAL*8 AJK(LRV,IDIM,IDIM,NPG) 71 : ENDSEGMENT 72 : 73 : -INC SMCHAML 74 : -INC SMCOORD 75 : -INC SMLENTI 76 : POINTEUR IPADI.MLENTI,IPADU.MLENTI,IPADF.MLENTI,IPADS.MLENTI 77 : POINTEUR IPADQ.MLENTI 78 : -INC SMELEME 79 : POINTEUR MELEMC.MELEME,MELEMS.MELEME,MELEMI.MELEME 80 : POINTEUR MELEMK.MELEME,MELEP1.MELEME,IGEOM0.MELEME 81 : -INC SMCHPOI 82 : POINTEUR IZG1.MCHPOI,IZGG1.MPOVAL 83 : POINTEUR IZTU1.MPOVAL,IZTU2.MPOVAL,IZTU3.MPOVAL 84 : POINTEUR VITESS.MPOVAL,UTRANS.MPOVAL 85 : POINTEUR MZNU.MPOVAL,MZGB.MPOVAL,MZTN.MPOVAL,MZTR.MPOVAL 86 : POINTEUR MZRO.MPOVAL 87 : POINTEUR MZNT.MPOVAL 88 : POINTEUR IZVOL.MPOVAL,IZTCO.MPOVAL 89 : POINTEUR MCHVOL.MCHPOI 90 : -INC SMMATRIK 91 : POINTEUR KMATRI.IMATRI,IPM.IZAFM 92 : 93 : SEGMENT IMATRS 94 : INTEGER LIZAFS(NBSOUS,NBME) 95 : ENDSEGMENT 96 : POINTEUR IPMS.IZAFM,IPS1.IZAFM,IPS2.IZAFM,IPS3.IZAFM 97 : 98 : -INC SMLMOTS 99 : POINTEUR LINCO.MLMOTS 100 : CHARACTER*8 NOMZ,NOMI,NOMA,TYPE,TYPC,NOM,NOM0 101 : CHARACTER*4 NOM4(5) 102 : PARAMETER (NTB=1) 103 : CHARACTER*8 LTAB(NTB) 104 : DIMENSION KTAB(NTB),IXV(3),RO(1) 105 : SAVE IPAS 106 : DATA LTAB/'KIZX '/,IPAS/0/ 107 : C***************************************************************************** 108 : CNSKE 109 : C write(6,*)' DEBUT NSKE ' 110 : RO(1)=1.d0 110 : CALL LITABS(LTAB,KTAB,NTB,1,IRET) 111 : IF (IERR.NE.0) RETURN 112 : MTABX=KTAB(1) 113 : C 114 : C- Récupération de la table EQEX (pointeur MTAB1) 115 : C 115 : CALL LEKTAB(MTABX,'EQEX',MTAB1) 116 : IF(MTAB1.EQ.0)THEN 117 : C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24 118 : MOTERR( 1: 8) = ' EQEX ' 119 : MOTERR( 9:16) = ' EQEX ' 120 : MOTERR(17:24) = ' KIZX ' 120 : CALL ERREUR(786) 121 : RETURN 122 : ENDIF 122 : CALL ACME(MTAB1,'NAVISTOK',NASTOK) 123 : IF(NASTOK.EQ.0)THEN 123 : CALL ZNSKE(MTABX,MTAB1) 124 : RETURN 125 : ENDIF 126 : C 127 : C- Récupération de la table INCO (pointeur KINC) 128 : C 128 : CALL LEKTAB(MTAB1,'INCO',KINC) 129 : IF(KINC.EQ.0)THEN 130 : C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24 131 : MOTERR( 1: 8) = ' INCO ' 132 : MOTERR( 9:16) = ' INCO ' 133 : MOTERR(17:24) = ' EQEX ' 133 : CALL ERREUR(786) 134 : RETURN 135 : ENDIF 136 : 137 : C***************************************************************************** 138 : C OPTIONS 139 : C KIMPL = 0 -> EXPL 1 -> IMPL 2 -> SEMI 140 : C KFORM = 0 -> SI 1 -> EF 2 -> VF 3 -> EFMC 141 : C IDCEN = 0-> rien 1-> CENTREE 2-> SUPGDC 3-> SUPG 4-> TVISQUEU 5-> CNG 142 : 143 : IAXI=0 144 : IF(IFOMOD.EQ.0)IAXI=2 145 : C 146 : C- Récupération de la table des options KOPT (pointeur KOPTI) 147 : C 147 : CALL LEKTAB(MTABX,'KOPT',KOPTI) 148 : IF (KOPTI.EQ.0) THEN 149 : C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24 150 : MOTERR( 1: 8) = ' KOPT ' 151 : MOTERR( 9:16) = ' KOPT ' 152 : MOTERR(17:24) = ' KIZX ' 152 : CALL ERREUR(786) 153 : RETURN 154 : ENDIF 155 : 155 : CALL ACME(KOPTI,'MTRMASS ',MMPG) 156 : IPG=0 157 : IF(MMPG.EQ.3)IPG=1 157 : CALL ACME(KOPTI,'IDCEN',IDCEN) 158 : KDIM=1 159 : IF(IDCEN.EQ.2)KDIM=IDIM 159 : CALL ACME(KOPTI,'RNG ',IMODEL) 159 : CALL ACME(KOPTI,'IKOMP',IKOMP) 159 : CALL ACME(KOPTI,'KIMPL',KIMPL) 159 : CALL ACME(KOPTI,'KPOIN',KPRE) 159 : CALL ACME(KOPTI,'KFORM',KFORM) 159 : CALL ACME(KOPTI,'KMACO',KMACO) 159 : CALL ACME(KOPTI,'IDIV',IDIV) 159 : CALL ACMF(KOPTI,'CMD',CMD) 159 : CALL ACMF(KOPTI,'AIMPL',AIMPL) 160 : AG=AIMPL 161 : AD=AIMPL-1.D0 162 : 163 : IF (IERR.NE.0) RETURN 164 : 165 : IF(KFORM.NE.0.AND.KFORM.NE.1)THEN 166 : C Option %m1:8 incompatible avec les données 167 : MOTERR( 1: 8) = 'EF/EFM1 ' 167 : CALL ERREUR(803) 168 : RETURN 169 : ENDIF 170 : 171 : IF(KFORM.NE.0.AND.KPRE.NE.2)THEN 172 : C Option %m1:8 incompatible avec les données 173 : MOTERR( 1: 8) = 'EFM1 ' 173 : CALL ERREUR(803) 174 : RETURN 175 : ENDIF 176 : 177 : C write(6,*)' Apres les options ' 178 : C***************************************************************************** 179 : C 180 : C- Récupération de la table DOMAINE associée au domaine local 181 : C 181 : CALL ACMM(MTABX,'NOMZONE',NOMZ) 181 : CALL LEKTAB(MTABX,'DOMZ',MTABZ) 182 : IF(MTABZ.EQ.0)THEN 183 : C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24 184 : MOTERR( 1: 8) = ' DOMZ ' 185 : MOTERR( 9:16) = ' DOMZ ' 186 : MOTERR(17:24) = ' KIZX ' 186 : CALL ERREUR(786) 187 : RETURN 188 : ENDIF 189 : 189 : CALL LEKTAB(MTABZ,'MAILLAGE',MELEME) 189 : CALL LEKTAB(MTABZ,'CENTRE',MELEMC) 189 : CALL LEKTAB(MTABZ,'SOMMET',MELEMS) 190 : IF (IERR.NE.0) RETURN 191 : 192 : MELEMQ=MELEMC 193 : MELEP1=MELEMC 194 : 195 : IF(KPRE.NE.2)THEN 195 : CALL LEKTAB(MTABZ,'MACRO',MACRO) 195 : CALL LEKTAB(MTABZ,'MACRO1',MELEMK) 195 : CALL LEKTAB(MTABZ,'QUADRATI',MQUAD) 196 : IF(KPRE.EQ.3)THEN 196 : CALL LEKTAB(MTABZ,'CENTREP0',MELEMQ) 197 : MELEP1=MELEMQ 198 : ELSEIF(KPRE.EQ.4)THEN 198 : CALL LEKTAB(MTABZ,'CENTREP1',MELEMQ) 198 : CALL LEKTAB(MTABZ,'ELTP1NC ',MELEP1) 199 : ENDIF 200 : ENDIF 201 : 202 : C************************************************************************* 203 : C VERIFICATIONS SUR LES INCONNUES 204 : C 205 : C- Récupération du nombre d'inconnues et du nom de l'inconnue NOMI 206 : C 207 : TYPE='LISTMOTS' 207 : CALL ACMO(MTABX,'LISTINCO',TYPE,LINCO) 208 : IF (IERR.NE.0) RETURN 209 : SEGACT LINCO 210 : NBINC=LINCO.MOTS(/2) 211 : IF(NBINC.NE.3)THEN 212 : C Indice %m1:8 : contient plus de %i1 %m9:16 213 : MOTERR( 1:8) = 'LISTINCO' 214 : INTERR(1) = 3 215 : MOTERR(9:16) = ' MOTS ' 215 : CALL ERREUR(799) 216 : RETURN 217 : ENDIF 218 : 219 : C 220 : C- Récupération des inconnues 221 : C 222 : C 223 : C 1ère inconnue UN vitesse 224 : C 225 : NOMI=LINCO.MOTS(1) 226 : NOMA=NOMI 227 : DO 15 I=1,IDIM 228 : WRITE(NOM4(I),FMT='(I1,A3)')I,NOMI(1:3) 229 : 15 CONTINUE 230 : 231 : TYPE=' ' 231 : CALL ACMO(KINC,NOMI,TYPE,MCHPOI) 232 : IF(TYPE.NE.'CHPOINT ')THEN 233 : C Indice %m1:8 : ne contient pas un objet de type %m9:16 234 : MOTERR( 1: 8) = 'INC '//NOMI 235 : MOTERR( 9:16) = 'CHPOINT ' 235 : CALL ERREUR(800) 236 : RETURN 237 : ELSE 237 : CALL LICHT(MCHPOI,IZTU1,TYPC,MELEMI) 238 : NINKO = IZTU1.VPOCHA(/2) 239 : NPTI = IZTU1.VPOCHA(/1) 240 : IF (NINKO.NE.IDIM) THEN 241 : C Indice %m1:8 : Le %m9:16 n'a pas le bon nombre de composantes 242 : MOTERR( 1: 8) = 'INC '//NOMI 243 : MOTERR( 9:16) = 'CHPOINT ' 243 : CALL ERREUR(784) 244 : RETURN 245 : ENDIF 246 : C write(6,*)' MCHPOI,MELEMI=',MCHPOI,MELEMI 247 : C On fait pointer ces deux tableaux sur le champ U inconu (tjs présent) pour 248 : C eviter de les enlever lors de l'appel FORTRAN si les options sont absentes 249 : UTRANS=IZTU1 250 : IKW=0 251 : VITESS=IZTU1 252 : ENDIF 253 : C 254 : C 2ème inconnue KN énergie turbulente 255 : C 256 : NOMI=LINCO.MOTS(2) 257 : NKN=IDIM+1 258 : NOM4(NKN)=NOMI(1:4) 259 : 260 : TYPE=' ' 260 : CALL ACMO(KINC,NOMI,TYPE,MCHPOI) 261 : IF(TYPE.NE.'CHPOINT ')THEN 262 : WRITE(6,*)' L objet CHPOINT ',NOMI,' n existe pas dans la table' 263 : C Indice %m1:8 : ne contient pas un objet de type %m9:16 264 : MOTERR( 1: 8) = 'INC '//NOMI 265 : MOTERR( 9:16) = 'CHPOINT ' 265 : CALL ERREUR(800) 266 : RETURN 267 : ELSE 267 : CALL LICHT(MCHPOI,IZTU2,TYPC,IGEOM0) 268 : IF(IGEOM0.NE.MELEMI)THEN 269 : segact igeom0,melemi 269 : call crech1(igeom0,0) 269 : call crech1(melemi,0) 270 : endif 271 : IF(IGEOM0.NE.MELEMI)THEN 272 : WRITE(6,*)' Opérateur NSKE ' 273 : WRITE(6,*)' Le SPG de l''inconnue 2 n''est pas convenable' 274 : C Indice %m1:8 : L'objet %m9:16 n'a pas le bon support géométrique 275 : MOTERR(1: 8) = 'INC '//NOMI 276 : MOTERR(9:16) = 'CHPOINT ' 276 : CALL ERREUR(788) 277 : RETURN 278 : ENDIF 279 : ENDIF 280 : C 281 : C 3ème inconnue EN dissipation de KN 282 : C 283 : NOMI=LINCO.MOTS(3) 284 : NEN=IDIM+2 285 : NOM4(NEN)=NOMI(1:4) 286 : 287 : TYPE=' ' 287 : CALL ACMO(KINC,NOMI,TYPE,MCHPOI) 288 : IF(TYPE.NE.'CHPOINT ')THEN 289 : WRITE(6,*)' L objet CHPOINT ',NOMI,' n existe pas dans la table' 290 : C Indice %m1:8 : ne contient pas un objet de type %m9:16 291 : MOTERR( 1: 8) = 'INC '//NOMI 292 : MOTERR( 9:16) = 'CHPOINT ' 292 : CALL ERREUR(800) 293 : RETURN 294 : ELSE 294 : CALL LICHT(MCHPOI,IZTU3,TYPC,IGEOM0) 295 : IF(IGEOM0.NE.MELEMI)THEN 296 : segact igeom0,melemi 296 : call crech1(igeom0,0) 296 : call crech1(melemi,0) 297 : endif 298 : IF(IGEOM0.NE.MELEMI)THEN 299 : WRITE(6,*)' Opérateur NSKE ' 300 : WRITE(6,*)' Le SPG de l''inconnue 3 n''est pas convenable' 301 : C Indice %m1:8 : L'objet %m9:16 n'a pas le bon support géométrique 302 : MOTERR(1: 8) = 'INC '//NOMI 303 : MOTERR(9:16) = 'CHPOINT ' 303 : CALL ERREUR(788) 304 : RETURN 305 : ENDIF 306 : ENDIF 307 : 308 : C***************************************************************************** 309 : C Le domaine de definition est donne par le SPG de la premiere inconnue 310 : C Les inconnues suivantes devront posseder ce meme pointeur 311 : C On verifie que les points de la zone sont tous inclus dans ce SPG 312 : 312 : CALL KRIPAD(MELEMI,IPADI) 313 : IPADS=IPADI 314 : NPTS = NPTI 315 : IF(MELEMI.NE.MELEMS)THEN 315 : CALL KRIPAD(MELEMS,IPADS) 316 : NPTS=MELEMS.NUM(/2) 317 : ENDIF 318 : IPADU=IPADI 319 : NPTU=NPTI 320 : 321 : IF(IPAS.EQ.0)THEN 322 : IPAS=1 322 : CALL VERPAD(IPADI,MELEME,IRET) 323 : IF(IRET.NE.0)THEN 324 : C Indice %m1:8 : L'objet %m9:16 n'a pas le bon support géométrique 325 : MOTERR(1: 8) = 'INC '//NOMI 326 : MOTERR(9:16) = 'CHPOINT ' 326 : CALL ERREUR(788) 327 : RETURN 328 : ENDIF 329 : ENDIF 330 : 331 : C************************************************************************* 332 : C Lecture du coefficient 333 : C Type du coefficient : 334 : C IK1=0 CHPOINT IK1=1 scalaire IK1=2 vecteur 335 : 336 : C write(6,*)' Opérateur NSKE lecture des coefficients' 336 : CALL ACME(MTABX,'IARG',IARG) 337 : IF(IKOMP.EQ.0)THEN 338 : IF(IARG.NE.2.AND.IARG.NE.3.AND.IARG.NE.5)THEN 339 : WRITE(6,*)' Opérateur NSKE : option incompréssible ' 340 : WRITE(6,*)' Nombre d''arguments incorrect : ',IARG 341 : WRITE(6,*)' On attend 2 , 3 ou 5 ' 342 : C Indice %m1:8 : nombre d'arguments incorrect 343 : MOTERR(1:8) = 'IARG ' 343 : CALL ERREUR(804) 344 : RETURN 345 : ENDIF 346 : ELSEIF(IKOMP.EQ.1)THEN 347 : IF(IARG.NE.4.AND.IARG.NE.5)THEN 348 : WRITE(6,*)' Opérateur NSKE : option compréssible ' 349 : WRITE(6,*)' Nombre d''arguments incorrect : ',IARG 350 : WRITE(6,*)' On attend 4 ou 5 ' 351 : C Indice %m1:8 : nombre d'arguments incorrect 352 : MOTERR(1:8) = 'IARG ' 352 : CALL ERREUR(804) 353 : RETURN 354 : ENDIF 355 : ENDIF 356 : 357 : C--Cas incompréssible 358 : IF(IKOMP.EQ.0)THEN 359 : 360 : MZRO=IZTU1 361 : 362 : C 1er coefficient : nu 363 : IXV(1)=MELEMC 364 : IXV(2)=1 365 : IXV(3)=0 365 : CALL LEKCOF('Opérateur NSKE :', 366 : & MTABX,KINC,1,IXV,MNU,MZNU,NPT1,NC1,IKN,IRET) 367 : IF(IRET.EQ.0)RETURN 368 : 369 : MZGB=MZNU 370 : MZTN=MZNU 371 : MZTR=MZNU 372 : 373 : C 2ème coefficient : nut 374 : IXV(1)=MELEMC 375 : IXV(2)=0 376 : IXV(3)=0 376 : CALL LEKCOF('Opérateur NSKE :', 377 : & MTABX,KINC,2,IXV,MNT,MZNT,NPT2,NC2,IKNT,IRET) 378 : IF(IRET.EQ.0)RETURN 379 : 380 : IF(IARG.GE.3)THEN 381 : C 3ème coefficient : gbeta 382 : IXV(1)=-MELEMC 383 : IXV(2)=0 384 : IXV(3)=1 384 : CALL LEKCOF('Opérateur NSKE :', 385 : & MTABX,KINC,3,IXV,MGB,MZGB,NELG,NC3,IKG,IRET) 386 : IF(IRET.EQ.0)RETURN 387 : 388 : IF(IKG.EQ.2)IKG=1 389 : 390 : IF(IARG.EQ.5)THEN 391 : 392 : C 4ème coefficient : tn 393 : IXV(1)=MELEMS 394 : IXV(2)=1 395 : IXV(3)=0 395 : CALL LEKCOF('Opérateur NSKE :', 396 : & MTABX,KINC,4,IXV,MTN,MZTN,NPT4,NC4,IKTN,IRET) 397 : IF(IRET.EQ.0)RETURN 398 : 399 : C 5ème coefficient : tref 400 : IXV(1)=MELEMS 401 : IXV(2)=1 402 : IXV(3)=0 402 : CALL LEKCOF('Opérateur NSKE :', 403 : & MTABX,KINC,5,IXV,MTR,MZTR,NPT5,NC5,IKTR,IRET) 404 : IF(IRET.EQ.0)RETURN 405 : 406 : ENDIF 407 : ENDIF 408 : 409 : C--Cas compréssible 410 : ELSEIF(IKOMP.EQ.1)THEN 411 : C 1er coefficient : ro 412 : IXV(1)=MELEMC 413 : IXV(2)=1 414 : IXV(3)=0 414 : CALL LEKCOF('Opérateur NSKE :', 415 : & MTABX,KINC,1,IXV,MRO,MZRO,NPT1,NCR,IKR,IRET) 416 : IF(IRET.EQ.0)RETURN 417 : 418 : MZGB=MZRO 419 : MZTN=MZRO 420 : MZTR=MZRO 421 : 422 : C 2ème coefficient : mu 423 : IXV(1)=MELEMC 424 : IXV(2)=1 425 : IXV(3)=0 425 : CALL LEKCOF('Opérateur NSKE :', 426 : & MTABX,KINC,2,IXV,MNU,MZNU,NPT2,NC2,IKN,IRET) 427 : IF(IRET.EQ.0)RETURN 428 : 429 : C 3ème coefficient : mut 430 : IXV(1)=MELEMC 431 : IXV(2)=0 432 : IXV(3)=0 432 : CALL LEKCOF('Opérateur NSKE :', 433 : & MTABX,KINC,3,IXV,MNT,MZNT,NPT3,NC3,IKNT,IRET) 434 : IF(IRET.EQ.0)RETURN 435 : 436 : C 4ème coefficient : un (en compréssible) 437 : IXV(1)=-MELEMS 438 : IXV(2)=0 439 : IXV(3)=0 439 : CALL LEKCOF('Opérateur NSKE :', 440 : & MTABX,KINC,4,IXV,MUN,UTRANS,NPT4,NC4,IK4,IRET) 441 : IF(IRET.EQ.0)RETURN 442 : IPADU=IPADS 443 : NPTU=NPTS 444 : 445 : IF(IARG.EQ.5)THEN 446 : C 5ème coefficient : s (en compréssible) 447 : IXV(1)=-MELEMC 448 : IXV(2)=0 449 : IXV(3)=1 449 : CALL LEKCOF('Opérateur NSKE :', 450 : & MTABX,KINC,5,IXV,MGB,MZGB,NELG,NC5,IKG,IRET) 451 : IF(IRET.EQ.0)RETURN 452 : 453 : IF(IKG.EQ.2)IKG=1 454 : ENDIF 455 : 456 : ENDIF 457 : 458 : C write(6,*)' Opérateur NSKE : Fin lecture Arguments ' 459 : C Fin lecture Arguments ************************************************ 460 : 461 : 461 : CALL LEKTAB(MTABZ,'XXDXDY',MCHPDX) 462 : IF(IRET.EQ.0)RETURN 462 : CALL LICHT(MCHPDX,IZTCO,TYPC,IGEOM) 463 : NELZ=IZTCO.VPOCHA(/1) 464 : 465 : 466 : C write(6,*)' FORMULATION ',kform 467 : C*********** FORMULATIONS ********** 468 : 469 : IF(KFORM.EQ.0)THEN 470 : C Formulation EFM1 471 : IF(KIMPL.NE.0)THEN 472 : WRITE(6,*)' Operateur NSKE' 473 : C Option %m1:8 incompatible avec les données 474 : MOTERR( 1: 8) = ' EFM1 ' 474 : CALL ERREUR(803) 475 : RETURN 476 : ENDIF 477 : 477 : CALL LEKTAB(MTABZ,'MATESI',MATRIK) 478 : IF(MATRIK.EQ.0)GO TO 90 479 : SEGACT MATRIK 480 : KMATRI=IRIGEL(4,1) 481 : SEGACT KMATRI 482 : SEGDES MATRIK 483 : 483 : CALL LEKTAB(MTABZ,'XXPSOML',MCHELM) 483 : CALL LEKTAB(MTABZ,'XXVOLUM',MCHVOL) 484 : IF(IRET.EQ.0)RETURN 485 : 486 : SEGACT MCHELM 486 : CALL LICHT(MCHVOL,IZVOL,TYPC,IGEOM) 487 : 488 : C***************************************************************************** 489 : C 490 : C Création du Chpoint second membre recevant l'incrément en explicite 491 : 492 : NC=IDIM+2 493 : TYPE='SOMMET' 493 : CALL KRCHPT(TYPE,MELEMS,NC,IZG1,NOM4) 494 : 494 : CALL LICHT(IZG1,IZGG1,TYPC,IGEOM) 495 : 496 : SEGDES LINCO 497 : 498 : 499 : C###################################################################### 500 : 501 : NRIGE=7 502 : NKID =9 503 : NKMT =7 504 : NMATRI=1 505 : SEGINI MATRIK 506 : IRIGEL(1,1)=MELEMS 507 : IRIGEL(2,1)=MELEMS 508 : IRIGEL(7,1)=5 509 : NBME=1 510 : NBSOUS=1 511 : SEGINI IMATRI 512 : IRIGEL(4,1)=IMATRI 513 : SEGACT MELEMS 514 : KSPGP=MELEMS 515 : KSPGD=MELEMS 516 : LISPRI(1)=NOM4(NEN)//' ' 517 : LISDUA(1)=NOM4(NEN)//' ' 518 : NP=1 519 : MP=1 520 : NBEL=MELEMS.NUM(/2) 521 : SEGINI IZAFM 522 : LIZAFM(1,1)=IZAFM 523 : SEGDES MATRIK,IMATRI 524 : 525 : SEGACT MELEME 526 : NBSOUS=LISOUS(/1) 527 : IF(NBSOUS.EQ.0)NBSOUS=1 528 : NUTOEL=0 529 : DT=1.D30 530 : IES=IDIM 531 : 532 : DO 1 L=1,NBSOUS 533 : IPT1=MELEME 534 : IF(NBSOUS.NE.1)IPT1=LISOUS(L) 535 : SEGACT IPT1 536 : IPM=KMATRI.LIZAFM(L,1) 537 : SEGACT IPM 538 : IPM1=IPM 539 : IF(IAXI.NE.0)THEN 540 : IPM1=KMATRI.LIZAFM(L,2) 541 : SEGACT IPM1 542 : ENDIF 543 : 544 : MCHAML=ICHAML(L) 545 : SEGACT MCHAML 546 : MELVAL=IELVAL(1) 547 : SEGACT MELVAL 548 : 549 : IF(IMACHE(L).NE.IPT1)THEN 550 : write(*,*)'IPT1,IMACHE ',IPT1,IMACHE(L) 551 : goto 90 552 : ENDIF 553 : 554 : NP =IPT1.NUM(/1) 555 : NBEL=IPT1.NUM(/2) 556 : 556 : CALL YCVIT(IPM.AM,IPM1.AM,VELCHE, 557 : & IPT1.NUM,NBEL,NUTOEL,NPTI,NPTS,NPTU,IES,NP,IAXI, 558 : & IPADI.LECT,IPADS.LECT,IPADU.LECT,IKOMP,IARG, 559 : & MZRO.VPOCHA,IKR, 560 : & MZNU.VPOCHA,IKN,MZGB.VPOCHA,IKG,NELG,MZTN.VPOCHA,IKTN, 561 : & MZTR.VPOCHA,IKTR,MZNT.VPOCHA, 562 : & IZTU1.VPOCHA,UTRANS.VPOCHA,IZTU2.VPOCHA,IZTU3.VPOCHA, 563 : & IZGG1.VPOCHA,IZGG1.VPOCHA(1,NKN),IZGG1.VPOCHA(1,NEN),AM, 564 : & IZVOL.VPOCHA,IZTCO.VPOCHA,NELZ,IDCEN,IMODEL, 565 : & DT,DTT1,DTT2,NUEL,DIAEL,ECHLM) 566 : 567 : SEGDES IPM,IPT1,MCHAML,MELVAL 568 : NUTOEL=NUTOEL+NBEL 569 : 570 : 1 CONTINUE 571 : SEGDES MELEME 572 : SEGDES KMATRI,MCHELM,IZAFM 573 : 574 : SEGDES MZTN 575 : SEGDES IZTU1,IZTU2,IZTU3 576 : SEGDES IZGG1 577 : SEGDES IZVOL,IZTCO 578 : 578 : CALL ECROBJ('MATRIK',MATRIK) 579 : 579 : CALL ECROBJ('CHPOINT',IZG1) 580 : 580 : CALL LEKTAB(MTAB1,'PASDETPS',MTABT) 581 : IF(MTABT.EQ.0)THEN 581 : CALL CRTABL(MTABT) 581 : CALL ECMM(MTABT,'SOUSTYPE','PASDETPS') 581 : CALL ECMO(MTAB1,'PASDETPS','TABLE ',MTABT) 582 : DTP=1.D30+DT 583 : IPAT=1 583 : CALL ECME(MTABT,'NUPASDT',IPAT) 584 : ELSE 584 : CALL ACMF(MTABT,'DELTAT',DTP) 585 : ENDIF 586 : 587 : IF(DT.LT.DTP)THEN 587 : CALL ECMF(MTABT,'DELTAT',DT) 587 : CALL ECMM(MTABT,'OPER','NSKE') 587 : CALL ACMM(MTABX,'NOMZONE',NOMZ) 587 : CALL ECMM(MTABT,'ZONE',NOMZ) 587 : CALL ECMF(MTABT,'DTCONV',DTT1) 587 : CALL ECMF(MTABT,'DTDIFU',DTT2) 587 : CALL ECMF(MTABT,'DIAEL',DIAEL) 587 : CALL ECME(MTABT,'NUEL',NUEL) 588 : ENDIF 589 : C************************************************************************* 590 : ELSE IF(KFORM.EQ.1)THEN 591 : C CAS FORMULATION EF 592 : 593 : IF(KIMPL.EQ.0)THEN 594 : C Option %m1:8 incompatible avec les données 595 : MOTERR( 1: 8) = ' EF ' 595 : CALL ERREUR(803) 596 : RETURN 597 : ENDIF 598 : 599 : DT=0.D0 600 : IF(IDCEN.EQ.4)THEN 601 : TYPE=' ' 601 : CALL ACMO(MTABT,'DELTAT',TYPE,IENT) 602 : IF(TYPE.NE.'ENTIER')THEN 603 : WRITE(6,*)' Opérateur NS ' 604 : WRITE(6,*)' On reclame un pas de temps' 605 : C Option %m1:8 incompatible avec les données 606 : MOTERR( 1: 8) = ' EF ' 606 : CALL ERREUR(803) 607 : RETURN 608 : ELSE 608 : CALL ACMF(MTABT,'DELTAT',DT) 609 : ENDIF 610 : ENDIF 611 : 612 : IHV=1 613 : NUTOEL=0 614 : NINKO=VITESS.VPOCHA(/2) 615 : SEGACT MELEME 616 : NBSOUS=LISOUS(/1) 617 : IF(NBSOUS.EQ.0)NBSOUS=1 618 : 619 : TYPE=' ' 619 : CALL ACMO(MTABX,'MATELM',TYPE,MATRIK) 620 : IF(TYPE.EQ.'MATRIK'.AND.KMACO.NE.0)THEN 621 : SEGACT MATRIK 622 : NMATRI=IRIGEL(/2) 623 : MELEME=IRIGEL(1,1) 624 : SEGACT MELEME 625 : IMATRI=IRIGEL(4,1) 626 : SEGACT IMATRI 627 : NBME=LIZAFM(/2) 628 : NINKO=NBME 629 : MELEMS=KSPGP 630 : 631 : 632 : ELSE 633 : 634 : NRIGE=7 635 : NKID =9 636 : NKMT =7 637 : NMATRI=1 638 : SEGINI MATRIK 639 : IRIGEL(1,1)=MELEME 640 : IRIGEL(2,1)=MELEME 641 : 642 : IRIGEL(7,1)=2 643 : NBOP=0 644 : NBME=NINKO 645 : NBELC=0 646 : SEGINI IMATRI,IMATRS 647 : IRIGEL(4,1)=IMATRI 648 : KSPGP=MELEMS 649 : KSPGD=MELEMS 650 : 651 : IF(NBME.EQ.1)THEN 652 : LISPRI(1)=NOMI(1:4)//' ' 653 : LISDUA(1)=NOMA(1:4)//' ' 654 : ELSE 655 : DO 102 I=1,NBME 656 : WRITE(NOM,FMT='(I1,A7)')I,NOMI(1:7) 657 : LISPRI(I)=NOM(1:4)//' ' 658 : WRITE(NOM,FMT='(I1,A7)')I,NOMA(1:7) 659 : LISDUA(I)=NOM(1:4)//' ' 660 : 102 CONTINUE 661 : ENDIF 662 : 663 : NUTOEL=0 664 : DO 101 L=1,NBSOUS 665 : IPT1=MELEME 666 : IF(NBSOUS.NE.1)IPT1=LISOUS(L) 667 : SEGACT IPT1 668 : WRITE(NOM0,FMT='(A8)')NOMS(IPT1.ITYPEL)//' ' 668 : CALL KALPBG(NOM0,'FONFORM ',IZFFM) 669 : SEGACT IZFFM*MOD 670 : IZHR=KZHR(1) 671 : SEGACT IZHR*MOD 672 : NES=GR(/1) 673 : NPG=GR(/3) 674 : 675 : NP = IPT1.NUM(/1) 676 : MP = NP 677 : NBEL=IPT1.NUM(/2) 678 : SEGINI IPM1,IPS1 679 : LIZAFM(L,1)=IPM1 680 : LIZAFS(L,1)=IPS1 681 : IPM2=IPM1 682 : IPM3=IPM1 683 : IPS2=IPS1 684 : IPS3=IPS1 685 : IF(NBME.GE.2)THEN 686 : SEGINI IPM2,IPS2 687 : LIZAFM(L,2)=IPM2 688 : LIZAFS(L,2)=IPS2 689 : ENDIF 690 : IF(NBME.GE.3)THEN 691 : SEGINI IPM3,IPS3 692 : LIZAFM(L,3)=IPM3 693 : LIZAFS(L,3)=IPS3 694 : ENDIF 695 : 696 : 697 : KITT=2 698 : KJTT=IKN 699 : NPT=UTRANS.VPOCHA(/1) 700 : 701 : SEGINI PETROV 701 : CALL XCONV(FN,GR,PG,XYZ,HR,PGSQ,RPG,AJ, 702 : & NES,IDIM,NP,NPG,IAXI,AG,AD,IDIV,CMD,IKOMP,LRV, 703 : & WT,WS,HK,PGSK,RPGK,AIRE,AJK,UIL,DUIL,COEFK,ANUK, 704 : & RO,1,UTRANS.VPOCHA,NPTU,IPADU.LECT,MZNU.VPOCHA,IK1, 705 : & IPT1.NUM,NBEL,NUTOEL,XCOOR, 706 : & IPM1.AM,IPM2.AM,IPM3.AM, 707 : & IPS1.AM,IPS2.AM,IPS3.AM, 708 : & NINKO,IDCEN,DT, 709 : & IZTU1.VPOCHA,NPTI,IPADI.LECT) 710 : SEGSUP PETROV 711 : 711 : CALL XLAPL(FN,GR,PG,XYZ,HR,PGSQ,RPG,NES,IDIM,NP,NPG,IAXI, 712 : & MZNU.VPOCHA,MZNU.VPOCHA,MZNU.VPOCHA,KITT,KJTT,IKN, 713 : & IPT1.NUM,NBEL,NUTOEL,XCOOR,AIMPL,IKOMP, 714 : & IPM1.AM,IPM2.AM,IPM3.AM, 715 : & IPS1.AM,IPS2.AM,IPS3.AM, 716 : & NINKO,IHV,IARG,MZNU.VPOCHA) 717 : 718 : NUTOEL=NUTOEL+NBEL 719 : 101 CONTINUE 720 : 721 : ENDIF 722 : 723 : 724 : IF(KIMPL.EQ.2.OR.KIMPL.EQ.0.OR.IARG.GT.1)THEN 725 : NAT=2 726 : NSOUPO=1 727 : SEGACT MELEMS 728 : N=MELEMS.NUM(/2) 729 : NC=NINKO 730 : SEGINI MCHPO1,MSOUP1,MPOVA1 731 : MCHPO1.IFOPOI=IFOMOD 732 : MCHPO1.MOCHDE=TITREE 733 : MCHPO1.MTYPOI='SMBR' 734 : MCHPO1.JATTRI(1)=2 735 : MCHPO1.IPCHP(1)=MSOUP1 736 : DO 177 N=1,NINKO 737 : MSOUP1.NOCOMP(N)=LISDUA(N)(1:4) 738 : 177 CONTINUE 739 : MSOUP1.IGEOC=MELEMS 740 : MSOUP1.IPOVAL=MPOVA1 741 : ENDIF 742 : 743 : IF(IARG.EQ.3.OR.IARG.EQ.5)THEN 743 : CALL KRIPAD(MELEMQ,IPADQ) 744 : IF(IARG.EQ.3)THEN 745 : IKAS=2 746 : ELSEIF(IARG.EQ.5)THEN 747 : IKAS=3 748 : ENDIF 749 : 750 : 751 : 752 : IF(MACRO.NE.0.AND.KPRE.EQ.2)MELEMK=MELEME 753 : SEGACT MELEMK 754 : NBSOUS=MELEMK.LISOUS(/1) 755 : IF(NBSOUS.EQ.0)NBSOUS=1 756 : 757 : NUTOEL=0 758 : DO 1102 L=1,NBSOUS 759 : IPT1=MELEMK 760 : IF(NBSOUS.NE.1)IPT1=MELEMK.LISOUS(L) 761 : SEGACT IPT1 762 : 763 : IF(MQUAD.NE.0)THEN 764 : IF(KPRE.EQ.2)WRITE(NOM0,FMT='(A8)')NOMS(IPT1.ITYPEL)//'PRP0' 765 : IF(KPRE.EQ.3)WRITE(NOM0,FMT='(A8)')NOMS(IPT1.ITYPEL)//'PRP0' 766 : IF(KPRE.EQ.4)WRITE(NOM0,FMT='(A8)')NOMS(IPT1.ITYPEL)//'PRP1' 767 : ELSEIF(MACRO.NE.0)THEN 768 : IF(KPRE.EQ.2)WRITE(NOM0,FMT='(A8)')NOMS(IPT1.ITYPEL)//' ' 769 : IF(KPRE.EQ.3)WRITE(NOM0,FMT='(A8)')NOMS(IPT1.ITYPEL)//'MCP0' 770 : IF(KPRE.EQ.4)WRITE(NOM0,FMT='(A8)')NOMS(IPT1.ITYPEL)//'MCP1' 771 : ELSE 772 : IF(KPRE.EQ.2)WRITE(NOM0,FMT='(A8)')NOMS(IPT1.ITYPEL)//' ' 773 : ENDIF 773 : CALL KALPBG(NOM0,'FONFORM ',IZFFM) 774 : 775 : 776 : SEGACT IZFFM*MOD 777 : IZHR=KZHR(1) 778 : SEGACT IZHR*MOD 779 : NES=GR(/1) 780 : NPG=GR(/3) 781 : IZF1=KTP(1) 782 : SEGACT IZF1*MOD 783 : MP1=IZF1.FN(/1) 784 : NP = IPT1.NUM(/1) 785 : MP = NP 786 : NBEL=IPT1.NUM(/2) 787 : NELG=MZGB.VPOCHA(/1) 788 : NPT =MPOVA1.VPOCHA(/1) 789 : 790 : SEGACT MELEP1 791 : 792 : 792 : CALL XSOUR(FN,IZF1.FN,GR,PG,XYZ,HR,PGSQ,RPG, 793 : & NES,IDIM,NP,MP1,NPG,IAXI,IPT1.NUM,IKAS,KPRE, 794 : & MZGB.VPOCHA,IKG,NELG,IPADQ.LECT,MELEP1.NUM, 795 : & MZTN.VPOCHA,IKTN,MZTR.VPOCHA,IKTR,IPADS.LECT, 796 : & NBEL,NUTOEL,XCOOR,MPOVA1.VPOCHA,NPT) 797 : 798 : NUTOEL=NUTOEL+NBEL 799 : 1102 CONTINUE 800 : 801 : ENDIF 802 : 803 : IF(KIMPL.EQ.2.OR.KIMPL.EQ.0)THEN 804 : 805 : NBSOUS=LISOUS(/1) 806 : IF(NBSOUS.EQ.0)NBSOUS=1 807 : 808 : DO 1533 L=1,NBSOUS 809 : IPT1=MELEME 810 : IF(NBSOUS.NE.1)IPT1=LISOUS(L) 811 : SEGACT IPT1 812 : NP=IPT1.NUM(/1) 813 : NBEL=IPT1.NUM(/2) 814 : DO 2 N=1,NINKO 815 : IPMS=LIZAFS(L,N) 816 : SEGACT IPMS 817 : DO 12 K=1,NBEL 818 : DO 13 J=1,NP 819 : UU=0.D0 820 : IU=IPADS.LECT(IPT1.NUM(J,K)) 821 : DO 14 I=1,NP 822 : IK=IPADI.LECT(IPT1.NUM(I,K)) 823 : UU=UU+IPMS.AM(K,I,J)*IZTU1.VPOCHA(IK,N) 824 : 14 CONTINUE 825 : MPOVA1.VPOCHA(IU,N)=MPOVA1.VPOCHA(IU,N)+UU 826 : 13 CONTINUE 827 : 12 CONTINUE 828 : 829 : 2 CONTINUE 830 : 831 : 1533 CONTINUE 832 : ENDIF 833 : 834 : SEGDES IPM1,IPM2,IPM3 835 : IPS=IPS1 836 : SEGSUP IPS1 837 : IF(IPS2.NE.IPS)SEGSUP IPS2 838 : IF(IPS3.NE.IPS)SEGSUP IPS3 839 : SEGDES IZTCO 840 : 840 : CALL ECROBJ('MATRIK',MATRIK) 841 : 842 : IF(KIMPL.EQ.2.OR.KIMPL.EQ.0.OR.IARG.GT.1)THEN 843 : 844 : C? TYPE=' ' 844 : C? CALL ACMO(MTAB1,'SMBR',TYPE,MCHPO2) 845 : C? IF(TYPE.NE.'CHPOINT')THEN 845 : C? CALL ECMO(MTAB1,'SMBR','CHPOINT',MCHPO1) 846 : C? ELSE 846 : C? CALL ECROBJ('CHPOINT',MCHPO2) 846 : C? CALL ECROBJ('CHPOINT',MCHPO1) 846 : C?? CALL OPERAD 846 : C? CALL PRFUSE 846 : C? CALL LIROBJ('CHPOINT',MCHPOI,1,IRET) 846 : C? CALL ECMO(MTAB1,'SMBR','CHPOINT',MCHPOI) 847 : C? ENDIF 848 : 848 : CALL ECROBJ('CHPOINT',MCHPO1) 849 : ELSE 850 : NAT=2 851 : NSOUPO=0 852 : SEGINI MCHPOI 853 : JATTRI(1)=2 853 : CALL ECROBJ('CHPOINT',MCHPOI) 854 : ENDIF 855 : 856 : SEGDES IMATRI 857 : SEGDES MELEME,MATRIK 858 : IF(IKN.EQ.0)THEN 859 : SEGDES MZNU 860 : ENDIF 860 : C? CALL ECMO(MTABX,'MATELM','MATRIK',MATRIK) 861 : 862 : C************************************************************************* 863 : ELSE IF(KFORM.EQ.2)THEN 864 : C CAS FORMULATION VF 865 : WRITE(6,*)' FORMULATION VF NON DISPONIBLE ' 866 : C Option %m1:8 incompatible avec les données 867 : MOTERR( 1: 8) = ' EF ' 867 : CALL ERREUR(803) 868 : RETURN 869 : ENDIF 870 : C************************************************************************* 871 : 872 : 873 : 874 : SEGSUP IPADI 875 : IPAS=1 876 : C write(6,*)' FIN NSKE ' 877 : RETURN 878 : 879 : 90 CONTINUE 880 : WRITE(6,*)' Interuption anormale de NSKE' 881 : C Option %m1:8 incompatible avec les données 882 : MOTERR( 1: 8) = ' EF ' 882 : CALL ERREUR(803) 883 : RETURN 884 : 1002 FORMAT(10(1X,1PE11.4)) 885 : END 886 : 887 : 888 : 889 : 890 : 891 : 892 : 893 : 894 : 895 : 896 : 897 : 898 : 899 : 900 : 901 : 902 : 903 : 904 : 905 :
© Cast3M 2003 - All rights reserved.
Disclaimer