1 : C YCLS SOURCE PV 12/07/09 21:15:05 7430 2 : SUBROUTINE YCLS(NOMPR,MTABX,IHV,MZIN,NOMI,MATRIK,MCHPO1) 3 : IMPLICIT INTEGER(I-N) 4 : IMPLICIT REAL*8 (A-H,O-Z) 5 : C*********************************************************************** 6 : C 7 : C CET OPERATEUR DISCRETISE Convection Diffusion + source 8 : C 9 : C SYNTAXE : 10 : C --------- 11 : C 12 : C 13 : C COEFFICIENTS : 14 : C -------------- 15 : C 16 : C 17 : C ALF (SCAL DOMA) DIFFUSIVITE THERMIQUE 18 : C (SCAL ELEM) 19 : C (SCAL NOEU) 20 : C 21 : C INCONNUES : 22 : C ----------- 23 : C 24 : C TN CHAMP DE TEMPERATURE 25 : C UN CHAMP DE VITESSE TRANSPORTANT 26 : C 27 : C NOMPR Nom de l'opérateur appelant (NS,TSCA,KONV ou LAPN) 28 : C suivant l'opérateur les coefs sont différents 29 : C Schéma en temps ou l'opérateur intervient : 30 : C SEMI CNG TVISQ 31 : C teta - schéma Crank Nicolson généralisé Tenseur visqueux 32 : C Ces schéma nécessitent que DFDT soit en CENTREE 33 : C Le Schéma BDF2 ne fait intervenir que DFDT 34 : C 35 : C E/ MTABX Table de l'opérateur -> coefficients 36 : C E/ IHV=0 SCALAIRE 37 : C IHV=1 VITESSE 38 : C E/ MZIN pointeur du CHPOINT de l'inconnue 39 : C E/ NOMI nom de l'inconnue 40 : C************************************************************************ 41 : 42 : -INC CCOPTIO 43 : -INC CCREEL 44 : -INC CCGEOME 45 : -INC SMCOORD 46 : -INC SMLENTI 47 : POINTEUR MLENT4.MLENTI 48 : -INC SMLMOTS 49 : -INC SMMODEL 50 : -INC SMELEME 51 : POINTEUR IGEOM.MELEME,MELEMS.MELEME,MELEMC.MELEME,MELEMP.MELEME 52 : POINTEUR MELEM2.MELEME 53 : -INC SMCHPOI 54 : POINTEUR MZIN.MCHPOI 55 : POINTEUR MTETA1.MPOVAL,MTETA2.MPOVAL,MTETA3.MPOVAL,MTETA4.MPOVAL 56 : -INC SMTABLE 57 : POINTEUR MTABZ.MTABLE 58 : -INC SIZFFB 59 : POINTEUR IZF1.IZFFM,IZH2.IZHR,IZW.IZFFM,IZWH.IZHR 60 : SEGMENT SAJT 61 : REAL*8 AJT(IDIM,IDIM,NPG),RF1(NP,MP,IDIM),SM1(NP,IDIM) 62 : REAL*8 TN1(NP,IDIM),TN2(NP,IDIM) 63 : REAL*8 WT(MP,NPG), WS(MP,NPG) 64 : REAL*8 WTC(IDIM*NPG,MP),WSC(IDIM*NPG,MP),HRD(IDIM*NPG,NP) 65 : REAL*8 FND(NPG,NP) ,FNE(NPG,NP,IDIM),HRE(NPG,NP,IDIM) 66 : REAL*8 GMU(NPG,IDIM),USI(IDIM) 67 : ENDSEGMENT 68 : -INC SMRIGID 69 : C-INC SMMATRIK 70 : C******************************************************************* 71 : C 72 : SEGMENT MATRIK 73 : REAL*8 COEMTK(NMATRI) 74 : INTEGER JRIGEL(NRIGE,NMATRI) 75 : INTEGER KSYM,KMINC,KMINCP,KMINCD,KIZM 76 : INTEGER KISPGT,KISPGP,KISPGD 77 : INTEGER KNTTT,KNTTP,KNTTD 78 : INTEGER KIDMAT(NKID) 79 : INTEGER KKMMT(NKMT) 80 : ENDSEGMENT 81 : 82 : SEGMENT JMATRI 83 : CHARACTER*8 LISPRI(NBMF),LJSDUA(NBMF) 84 : INTEGER LIZAFM(NBSOUS,NBMF) 85 : INTEGER KSPGP,KSPGD 86 : ENDSEGMENT 87 : POINTEUR JMATRS.JMATRI,JMATR1.JMATRI,JMATR2.JMATRI,JMATR3.JMATRI 88 : 89 : C Stokage matrices elementaires non assemblees (valeurs) 90 : SEGMENT IZAFM 91 : REAL*8 AM(NBEL,NP,MP) 92 : ENDSEGMENT 93 : POINTEUR IPM1.IZAFM,IPM2.IZAFM,IPM3.IZAFM,IPM4.IZAFM 94 : 95 : C******************************************************************* 96 : POINTEUR IPM.IZAFM 97 : -INC SMCHAML 98 : POINTEUR MCHELG.MCHELM,MCHAMG.MCHAML,MELVAG.MELVAL 99 : 100 : LOGICAL TDFDT,TCONV,TSOUR,TSOURB,TMATRI 101 : LOGICAL ICAL,XPG,XRIG 102 : LOGICAL XDIAG,XTV,XTG,XBDF,XCONS 103 : CHARACTER*8 CHHH,TYPE,NOM,NOM0,CHAI,SCHT,NOMPER,NOMPR,MPRE 104 : CHARACTER*4 INCOD 105 : CHARACTER*4 NOMI(*),NOMA(3),NOMP(3),NOMD(3),NMACO 106 : DIMENSION ICOEF(4),SU(3),XPOI(3) 107 : C***************************************************************************** 108 : CTCLSF 109 : 110 : NOMPER=NOMPR 111 : c write(6,*)' DEBUT YCLS appele par ',NOMPER 112 : 113 : C- Table de l'opérateur (pointeur MTABX) 114 : 115 : C- Récupération de la table RV (pointeur MTAB1) 115 : CALL LEKTAB(MTABX,'EQEX',MTAB1) 116 : C- Récupération de la table des options de l'opérateur (pointeur KOPTI) 116 : CALL LEKTAB(MTABX,'KOPT',KOPTI) 117 : C- Récupération de la table DOMAINE de la zone (pointeur MTABZ) 117 : CALL LEKTAB(MTABX,'DOMZ',MTABZ) 118 : C- Récupération de la table INCO (pointeur KINC) 118 : CALL LEKTAB(MTAB1,'INCO',KINC) 119 : IF(IERR.NE.0) RETURN 120 : 121 : C***************************************************************************** 122 : C Traitement des options 123 : C KIMPL = 0 -> EXPL 1 -> IMPL 2 -> CN 124 : C AIMPL=0 explicite 125 : C IKOMP = 0 Formulation non conservative des termes de convection 126 : C IKOMP = 1 Formulation faible conservative des termes de convection 127 : C i.e. integre par partie on fait apparaitre l'integrale de surface du 128 : C flux convectif ce qui permet d'imposer le cas echeant le flux total 129 : C (convectif + diffusif) 130 : C IKOMP = 2 Correction de Temam pour la stabilite L2 131 : C KMU = 0 QDM : mu = cste sur l'element 132 : C KMU = 1 QDM : Formulation en Tau (contrainte visqueuse) 133 : C KMU = 2 QDM : Formulation en Grad de Mu 134 : C KKALE=1 Indice indiquant que l'on est en ALE 135 : C sinon KKALE=0 136 : C KFORM = 0 -> SI 1 -> EF 2 -> VF 3 -> EFMC 137 : C IDCEN : entier indiquant le type de décentrement souhaité 138 : C IDCEN 1-> CENTREE 2-> SUPGDC 3-> SUPG 4-> TVISQUEU 5-> CNG 139 : C E/ CMD : Coefficient multiplicateur du décentrement 140 : C Si IDCEN=4 ou =5 CMD=DT 141 : 141 : CALL ACME(KOPTI,'KFORM',KFORM) 141 : CALL ACME(KOPTI,'KIMPL',KIMPL) 141 : CALL ACMF(KOPTI,'AIMPL',AIMPL) 141 : CALL ACME(KOPTI,'IKOMP',IKOMP) 141 : CALL ACME(KOPTI,'KMACO',KMACO) 142 : KKMACO=KMACO 142 : CALL ACME(KOPTI,'IDCEN',IDCEN) 142 : CALL ACMF(KOPTI,'CMD',CMD) 142 : CALL ACME(KOPTI,'ALE',KKALE) 142 : CALL ACME(KOPTI,'KMU',KMU) 143 : IF(KMU.EQ.2.AND.IHV.EQ.1)THEN 144 : IHM=1 145 : ELSE 146 : IHM=0 147 : ENDIF 148 : 149 : C Restrictions des options 150 : IF(NOMPER.NE.'NS')KKALE=0 151 : IF(NOMPER.EQ.'LAPN')IDCEN=1 152 : 153 : IF(KIMPL.EQ.0)AIMPL=0.D0 154 : IF(KIMPL.EQ.1)AIMPL=1.D0 155 : 156 : c write(6,*)' YCLS KMACO=',KMACO,'IDCEN=',IDCEN 157 : IF(KMACO.EQ.1.AND.IDCEN.NE.5)THEN 157 : CALL ACMM(KOPTI,'NMACO',NMACO) 158 : TYPE=' ' 158 : CALL ACMO(MTAB1,NMACO,TYPE,MATRIK) 159 : IF(TYPE.EQ.'MATRIK')THEN 160 : KKMACO=2 161 : ELSE 162 : KKMACO=1 163 : ENDIF 164 : ENDIF 165 : IF(IDCEN.EQ.5)KKMACO=0 166 : c write(6,*)' YCLS KMACO=',kmaco,' KKMACO=',KKMACO,'IKOMP=',ikomp 167 : 168 : C***************************************************************************** 169 : 170 : 171 : 172 : KPOIND=0 173 : NBMF=1 174 : XRIG=.FALSE. 175 : TMATRI=.TRUE. 176 : C C2LAPN=1.D0 Laplacien C2LAPN=0.D0 pas de laplacien 177 : C2LAPN=1.D0 178 : TCONV=.TRUE. 179 : TSOUR=.FALSE. 180 : TSOURB=.FALSE. 181 : XPG=.FALSE. 182 : IF(IDCEN.GT.1)XPG=.TRUE. 183 : IF(IDCEN.EQ.4)AIMPL=0.D0 184 : C C1CNG cas particulier de Crank Nicolson généralisé 185 : C1CNG=1.D0 186 : IF(IDCEN.EQ.5)THEN 187 : AIMPL=1.D0 188 : C1CNG=0.5D0 189 : ENDIF 190 : 191 : AIEX=AIMPL 192 : IF(AIMPL.EQ.0.D0.AND.KKMACO.NE.0)THEN 193 : TMATRI=.TRUE. 194 : AIEX=1.D0 195 : ENDIF 196 : 197 : IF(AIMPL.EQ.0.D0.AND.KKMACO.EQ.0)THEN 198 : TMATRI=.FALSE. 199 : ENDIF 200 : 201 : IAXI=0 202 : IF(IFOMOD.EQ.0)IAXI=2 203 : DEUPI=1.D0 204 : IF(IAXI.NE.0)DEUPI=2.D0*XPI 205 : 206 : 207 : IF(XRIG)THEN 208 : IF(IHV.EQ.0)THEN 209 : NOMP(1)='T ' 210 : NOMD(1)='Q ' 211 : ELSEIF(IHV.EQ.1)THEN 212 : NOMP(1)='UX ' 213 : NOMP(2)='UY ' 214 : NOMP(3)='UZ ' 215 : NOMD(1)='FX ' 216 : NOMD(2)='FY ' 217 : NOMD(3)='FZ ' 218 : ENDIF 219 : ELSE 220 : IF(IHV.EQ.0)THEN 221 : NOMP(1)=NOMI(1) 222 : NOMD(1)=NOMI(1) 223 : ELSEIF(IHV.EQ.1)THEN 224 : WRITE(NOMP(1),FMT='(I1,A3)')1,NOMI(1)(1:3) 225 : WRITE(NOMD(1),FMT='(I1,A3)')1,NOMI(1)(1:3) 226 : WRITE(NOMP(2),FMT='(I1,A3)')2,NOMI(1)(1:3) 227 : WRITE(NOMD(2),FMT='(I1,A3)')2,NOMI(1)(1:3) 228 : WRITE(NOMP(3),FMT='(I1,A3)')3,NOMI(1)(1:3) 229 : WRITE(NOMD(3),FMT='(I1,A3)')3,NOMI(1)(1:3) 230 : ENDIF 231 : ENDIF 232 : c write(6,*)' NOMP=',nomp 233 : c write(6,*)' NOMD=',nomd 234 : 234 : CALL LEKTAB(MTABZ,'SOMMET',MELEMS) 234 : CALL LEKTAB(MTABZ,'MAILLAGE',MELEME) 235 : MELEM2=MELEME 236 : 237 : C On cree un second membre non vide 238 : SEGACT MELEMS 239 : N=MELEMS.NUM(/2) 240 : IF(IHV.EQ.0)NC=1 241 : IF(IHV.EQ.1)NC=IDIM 242 : SEGDES MELEMS 243 : NSOUPO=1 244 : NAT=2 245 : SEGINI MCHPO1,MSOUP1,MPOVA1 246 : MCHPO1.JATTRI(1)=2 247 : MCHPO1.IFOPOI=IFOMOD 248 : MCHPO1.MTYPOI=' ' 249 : MCHPO1.MOCHDE=' ' 250 : MCHPO1.IPCHP(1)=MSOUP1 251 : MSOUP1.IGEOC=MELEMS 252 : MSOUP1.IPOVAL=MPOVA1 253 : IF(IHV.EQ.0)THEN 254 : MSOUP1.NOCOMP(1)=NOMD(1) 255 : ELSE 256 : DO 91 N=1,NC 257 : MSOUP1.NOCOMP(N)=NOMD(N) 258 : 91 CONTINUE 259 : ENDIF 260 : SEGDES MSOUP1,MCHPO1 261 : 262 : C************************************************************************* 263 : C Lecture des Inco(s) aux temps precedents si Transitoire 264 : c write(6,*)'Lecture des INCO ' 265 : 265 : CALL LRCHT(MZIN,MTETA1,TYPE,IGEOM) 265 : CALL KRIPAD(IGEOM,MLENT1) 265 : CALL VERPAD(MLENT1,MELEMS,IRET) 266 : IF(IRET.NE.0)THEN 267 : C% Indice %m1:8 : L'objet %m9:16 n'a pas le bon support géométrique 268 : MOTERR(1: 8) = 'TETA1' 269 : MOTERR(9:16) = 'CHPOINT ' 270 : WRITE(IOIMP,*)'Opérateur : ',NOMPER 270 : CALL ERREUR(788) 271 : RETURN 272 : ENDIF 273 : 274 : NC=MTETA1.VPOCHA(/2) 275 : 276 : IF((IHV.EQ.0.AND.NC.NE.1).OR. 277 : & (IHV.EQ.1.AND.NC.NE.IDIM))THEN 278 : C% Indice %m1:8 : L'objet %m9:16 n'a pas le bon nombre de composantes 279 : MOTERR( 1: 8) = 'Inconnue' 280 : MOTERR( 9:16) = 'CHPOINT' 281 : WRITE(IOIMP,*)'Operateur : ',NOMPER 281 : CALL ERREUR(784) 282 : RETURN 283 : ENDIF 284 : 285 : SEGDES IGEOM,MELEMS 286 : 287 : C***************************************************************************** 288 : C Lecture du coefficient 289 : c write(6,*)'Lecture des COEF ' 290 : NMUVARI=0 291 : MCHELG =0 292 : 292 : CALL ACME(MTABX,'IARG',IARG) 293 : c write(6,*)' YCLS IARG=',IARG,' KKALE=',KKALE 294 : 295 : C--------- CAS NS ou TSCA Formulation non conservative ou conservative ------- 296 : IF((NOMPER.EQ.'NS'.OR.NOMPER.EQ.'TSCA').AND. 297 : & (IKOMP.EQ.0.OR.IKOMP.EQ.1))THEN 298 : 299 : c_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._ 300 : IF( (IHV.EQ.0.AND.(IARG.EQ.3.OR.IARG.EQ.4)).OR. 301 : c (IHV.EQ.1.AND.(IARG.EQ.3.OR.IARG.EQ.4.OR.IARG.EQ.6)) 302 : c )THEN 303 : C2LAPN=1.D0 304 : TCONV=.TRUE. 305 : TSOUR=.FALSE. 306 : TSOURB=.FALSE. 307 : IF(IARG.EQ.4)TSOUR=.TRUE. 308 : IF(IHV.EQ.1.AND.IARG.EQ.6)THEN 309 : TSOUR=.TRUE. 310 : TSOURB=.TRUE. 311 : ENDIF 312 : C Ro 313 : NUCOEF=1 313 : CALL LEKMOF(MTABZ,NUCOEF,MTABX,KINC,0,MCHEL1,KPOIND,0,MUVARI) 314 : IF (IERR.NE.0) RETURN 315 : C Un 316 : NUCOEF=NUCOEF+1 316 : CALL LEKMOF(MTABZ,NUCOEF,MTABX,KINC,1,MCHEL3,KPOIND,0,MUVARI) 317 : IF (IERR.NE.0) RETURN 318 : C Mu 319 : NUCOEF=NUCOEF+1 319 : CALL LEKMOF(MTABZ,NUCOEF,MTABX,KINC,0,MCHEL2,KPOIND,IHM,MUVARI) 320 : IF(MUVARI.NE.0)THEN 321 : MCHELG=MUVARI 322 : NMUVARI=NMUVARI+1 323 : ENDIF 324 : IF (IERR.NE.0) RETURN 325 : 326 : ELSE 327 : WRITE(6,*)' Nombre d''arguments incorrect : ',IARG 328 : IF(IHV.EQ.0)THEN 329 : WRITE(6,*)' On attend 3 ou 4' 330 : ENDIF 331 : IF(IHV.EQ.1)THEN 332 : WRITE(6,*)' On attend 3 ou 6' 333 : ENDIF 334 : C Indice %m1:8 : nombre d'arguments incorrect 335 : MOTERR(1:8) = NOMPER 335 : CALL ERREUR(804) 336 : RETURN 337 : ENDIF 338 : c_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._ 339 : 340 : 341 : IF(TSOUR)THEN 342 : C Source -> MCHEL4 343 : C Gb 344 : NUCOEF=NUCOEF+1 345 : c write(6,*)' Coef 3 ',NUCOEF,' MTABZ=',MTABZ 346 : IVC=IHV 346 : CALL LEKMOF(MTABZ,NUCOEF,MTABX,KINC,IVC,MCHEL4,KPOIND,0,MUVARI) 347 : IF (IERR.NE.0) RETURN 348 : 349 : IF(TSOURB)THEN 350 : C Source Boussinesq pour la QDM : Gb(T0 - T) -> MCHEL4 351 : C Tn 352 : c write(6,*)' Coef 4 BOU',NUCOEF 353 : NUCOEF=NUCOEF+1 353 : CALL LEKMOF(MTABZ,NUCOEF,MTABX,KINC,0,MCHEL5,KPOIND,0,MUVARI) 354 : IF (IERR.NE.0) RETURN 355 : C To 356 : NUCOEF=NUCOEF+1 356 : CALL LEKMOF(MTABZ,NUCOEF,MTABX,KINC,0,MCHEL6,KPOIND,0,MUVARI) 357 : IF (IERR.NE.0) RETURN 358 : C -> MCHEL4 = Gb*(Tn - T0) 358 : CALL MELBOU(MTABZ,MCHEL4,MCHEL4,MCHEL5,MCHEL6) 359 : ENDIF 360 : ENDIF 361 : 362 : C write(6,*)'FIN CAS NS ou TSCA Formulation non conservative' 363 : C----- FIN CAS NS ou TSCA Formulation non conservative ou non conservative --- 364 : 365 : C--------- CAS NS ou TSCA Formulation Temam ---------------------------------- 366 : ELSEIF(( NOMPER.EQ.'TSCA').AND.IKOMP.EQ.2)THEN 367 : 368 : TCONV=.TRUE. 369 : TSOUR=.TRUE. 370 : TSOURB=.FALSE. 371 : C On convient d'avoir 3 coefs tout le temps (mu , u , s) 372 : 373 : C Nu 374 : NUCOEF=1 375 : c write(6,*)' Coef 1 ',NUCOEF,' MTABZ=',MTABZ 375 : CALL LEKMOF(MTABZ,NUCOEF,MTABX,KINC,0,MCHEL2,KPOIND,IHM,MUVARI) 376 : IF(MUVARI.NE.0)THEN 377 : MCHELG=MUVARI 378 : NMUVARI=NMUVARI+1 379 : ENDIF 380 : IF (IERR.NE.0) RETURN 381 : 382 : SEGACT MELEME 383 : L1=72 384 : N1=MAX(1,LISOUS(/1)) 385 : N2=1 386 : N3=6 387 : SEGINI MCHEL1 388 : SEGACT MCHEL2 389 : DO 75 L=1,N1 390 : MCHAM2=MCHEL2.ICHAML(L) 391 : SEGACT MCHAM2 392 : MELVA2=MCHAM2.IELVAL(1) 393 : SEGACT MELVA2 394 : N1PTEL=MELVA2.VELCHE(/1) 395 : N1EL=1 396 : N2PTEL=0 397 : N2EL=0 398 : SEGINI MCHAM1,MELVA1 399 : MCHEL1.ICHAML(L)=MCHAM1 400 : MCHAM1.IELVAL(1)=MELVA1 401 : DO 76 LG=1,N1PTEL 402 : MELVA1.VELCHE(LG,1)=1.D0 403 : 76 CONTINUE 404 : 75 CONTINUE 405 : 406 : C Un 407 : NUCOEF=NUCOEF+1 408 : c write(6,*)' Coef 2 ',NUCOEF 408 : CALL LEKMOF(MTABZ,NUCOEF,MTABX,KINC,1,MCHEL3,KPOIND,0,MUVARI) 409 : IF (IERR.NE.0) RETURN 409 : CALL LEKDIV(MTABZ,NUCOEF,MTABX,KINC,MCHEL5,KPOIND) 410 : IF (IERR.NE.0) RETURN 411 : 412 : IF(TSOUR)THEN 413 : C Source -> MCHEL4 414 : C Gb 415 : NUCOEF=NUCOEF+1 416 : c write(6,*)' Coef 3 ',NUCOEF,' MTABZ=',MTABZ 417 : IVC=IHV 417 : CALL LEKMOF(MTABZ,NUCOEF,MTABX,KINC,IVC,MCHEL4,KPOIND,0,MUVARI) 418 : IF (IERR.NE.0) RETURN 419 : ENDIF 420 : 421 : C----- FIN CAS NS ou TSCA Formulation Temam ---------------------------------- 422 : 423 : C--------- CAS LAPN ---------------------------------------------------------- 424 : ELSEIF(NOMPER.EQ.'LAPN')THEN 425 : 426 : TSOUR=.FALSE. 427 : TSOURB=.FALSE. 428 : TCONV =.FALSE. 429 : C Alfa 430 : NUCOEF=1 431 : c write(6,*)' Coef 1 ',NUCOEF 431 : CALL LEKMOF(MTABZ,NUCOEF,MTABX,KINC,0,MCHEL2,KPOIND,IHM,MUVARI) 432 : IF(MUVARI.NE.0)THEN 433 : MCHELG=MUVARI 434 : NMUVARI=NMUVARI+1 435 : ENDIF 436 : IF (IERR.NE.0) RETURN 437 : 438 : SEGACT MELEME 439 : L1=72 440 : N1=MAX(1,LISOUS(/1)) 441 : N2=1 442 : N3=6 443 : SEGINI MCHEL1,MCHEL3 444 : SEGACT MCHEL2 445 : DO 74 L=1,N1 446 : MCHAM2=MCHEL2.ICHAML(L) 447 : SEGACT MCHAM2 448 : MELVA2=MCHAM2.IELVAL(1) 449 : SEGACT MELVA2 450 : N1PTEL=MELVA2.VELCHE(/1) 451 : N1EL=1 452 : N2PTEL=0 453 : N2EL=0 454 : SEGINI MCHAM1,MELVA1 455 : MCHEL1.ICHAML(L)=MCHAM1 456 : MCHAM1.IELVAL(1)=MELVA1 457 : DO 73 LG=1,N1PTEL 458 : MELVA1.VELCHE(LG,1)=0.D0 459 : 73 CONTINUE 460 : N1PTEL=N1PTEL*IDIM 461 : SEGINI MCHAM3,MELVA3 462 : MCHEL3.ICHAML(L)=MCHAM3 463 : MCHAM3.IELVAL(1)=MELVA3 464 : DO 731 LG=1,N1PTEL 465 : MELVA3.VELCHE(LG,1)=0.D0 466 : 731 CONTINUE 467 : 74 CONTINUE 468 : C----- FIN CAS LAPN ---------------------------------------------------------- 469 : 470 : C--------- CAS KONV ---------------------------------------------------------- 471 : ELSEIF(NOMPER.EQ.'KONV')THEN 472 : 473 : IF(IARG.NE.3.AND.IARG.NE.4)THEN 474 : WRITE(6,*)' Nombre d''arguments incorrect : ',IARG 475 : WRITE(6,*)' On attend 3 ou 4' 476 : C Indice %m1:8 : nombre d'arguments incorrect 477 : MOTERR(1:8) = NOMPER 477 : CALL ERREUR(804) 478 : RETURN 479 : ENDIF 480 : 481 : C2LAPN=0.D0 482 : TCONV=.TRUE. 483 : TSOUR=.FALSE. 484 : TSOURB=.FALSE. 485 : C Ro 486 : NUCOEF=1 487 : c write(6,*)' Coef 1 ',NUCOEF 487 : CALL LEKMOF(MTABZ,NUCOEF,MTABX,KINC,0,MCHEL1,KPOIND,0,MUVARI) 488 : IF (IERR.NE.0) RETURN 489 : 490 : C Un 491 : NUCOEF=NUCOEF+1 492 : c write(6,*)' Coef 2 ',NUCOEF 492 : CALL LEKMOF(MTABZ,NUCOEF,MTABX,KINC,1,MCHEL3,KPOIND,0,MUVARI) 493 : IF (IERR.NE.0) RETURN 494 : C Mu 495 : NUCOEF=NUCOEF+1 496 : c write(6,*)' Coef 3 ',NUCOEF 496 : CALL LEKMOF(MTABZ,NUCOEF,MTABX,KINC,0,MCHEL2,KPOIND,0,MUVARI) 497 : IF (IERR.NE.0) RETURN 498 : 499 : C Dt 500 : IF(IDCEN.EQ.4.OR.IDCEN.EQ.5)THEN 501 : NUCOEF=NUCOEF+1 502 : c write(6,*)' Coef 4 ',NUCOEF 502 : CALL LEKMOF(MTABZ,NUCOEF,MTABX,KINC,0,MCHEL4,KPOIND,0,MUVARI) 503 : IF (IERR.NE.0) RETURN 504 : SEGACT MCHEL4 505 : MCHAM4=MCHEL4.ICHAML(1) 506 : SEGACT MCHAM4 507 : MELVA4=MCHAM4.IELVAL(1) 508 : SEGACT MELVA4 509 : N1PTEL=MELVA4.VELCHE(/1) 510 : N1EL=MELVA4.VELCHE(/2) 511 : IF(N1EL.EQ.1)THEN 512 : DT= MELVA4.VELCHE(1,1) 513 : SEGSUP MELVA4,MCHAM4,MCHEL4 514 : ELSE 514 : call erreur(991) 515 : RETURN 516 : ENDIF 517 : 518 : CMD=DT 519 : ENDIF 520 : 521 : C----- FIN CAS KONV ---------------------------------------------------------- 522 : 523 : ENDIF 524 : c write(6,*)' Fin lect COEF ' 525 : IF (IERR.NE.0) RETURN 526 : C************************************************************************* 527 : C******* CALCUL ********************************************************** 528 : C************************************************************************* 529 : 530 : SEGACT MTETA1 531 : 531 : CALL KRIPAD(MELEMS,MLENTI) 532 : SEGDES MELEMS 533 : 534 : SEGACT MELEME,MELEM2 535 : NBSOUS=LISOUS(/1) 536 : IF(NBSOUS.EQ.0)NBSOUS=1 537 : 538 : C MATRIK . MATRIK . MATRIK . MATRIK . MATRIK . MATRIK MATRIK . MATRIK . 539 : IF(KKMACO.NE.2)THEN 540 : C ... Création / Ajustement des matrices élémentaires 541 : IF(TMATRI)THEN 542 : C Cas RIGIDITE 543 : IF(XRIG)THEN 544 : 545 : NRIGE=8 546 : NRIGEL=0 547 : SEGINI MRIGID 548 : 549 : NRIGE=8 550 : NRI=IRIGEL(/2) 551 : NRIGEL=NBSOUS+NRI 552 : SEGADJ MRIGID 553 : c write(6,*)' NRIGE,NRIGEL,MRIGID=',NRIGE,NRIGEL,MRIGID 554 : 555 : ELSE 556 : C Cas MATRIK 557 : 558 : NRIGE=7 559 : NKID =9 560 : NKMT =7 561 : NMATRI=0 562 : SEGINI MATRIK 563 : c write(6,*)' YCLS creation MATRIK =',matrik 564 : 565 : NRIGE=7 566 : NKID =9 567 : NKMT =7 568 : NMATR0=JRIGEL(/2) 569 : 570 : NBME=1 571 : IF(IHV.EQ.1)NBME=IDIM 572 : NMATRI=NBME 573 : SEGADJ MATRIK 574 : 575 : DO 41 M=1,NBME 576 : JRIGEL(1,M)=MELEME 577 : JRIGEL(2,M)=MELEM2 578 : JRIGEL(7,M)=0 579 : IF(TCONV)JRIGEL(7,M)=2 580 : SEGINI JMATRI 581 : JRIGEL(4,M)=JMATRI 582 : KSPGP=MELEMS 583 : KSPGD=MELEMS 584 : LISPRI(1)=NOMP(M)//' ' 585 : LJSDUA(1)=NOMD(M)//' ' 586 : 41 CONTINUE 587 : ENDIF 588 : ENDIF 589 : C ... Fin Création / Ajustement des matrices élémentaires 590 : 591 : 592 : ELSE 593 : C Activation RIGIDITE 594 : IF(XRIG)THEN 595 : SEGACT MRIGID 596 : ELSE 597 : C Activation MATRIK 598 : SEGACT MATRIK 599 : ENDIF 600 : 601 : C ... Activation des matrices élémentaires 602 : IF(TMATRI)THEN 603 : C Cas RIGIDITE 604 : IF(XRIG)THEN 605 : NRIGE=8 606 : NRI=IRIGEL(/2) 607 : NRIGEL=NBSOUS+NRI 608 : SEGAct MRIGID 609 : c write(6,*)' NRIGE,NRIGEL,MRIGID=',NRIGE,NRIGEL,MRIGID 610 : 611 : ELSE 612 : C Cas MATRIK Pleine 613 : NMATRI=JRIGEL(/2) 614 : NBME=NMATRI 615 : 616 : DO 42 M=1,NBME 617 : JMATRI=JRIGEL(4,M) 618 : SEGACT JMATRI*MOD 619 : LISPRI(1)=NOMP(M)//' ' 620 : LJSDUA(1)=NOMD(M)//' ' 621 : 42 CONTINUE 622 : 623 : ENDIF 624 : ENDIF 625 : C ... Fin Activation des matrices élémentaires 626 : 627 : ENDIF 628 : C FIN MATRIK . FIN MATRIK . FIN MATRIK . FIN MATRIK . FIN MATRIK . FIN M 629 : C ...................................................................... 630 : 631 : SEGACT MCHEL1 632 : SEGACT MCHEL2 633 : IF(MCHELG.NE.0)THEN 634 : IF(NMUVARI.GT.1)THEN 635 : WRITE(6,*)' Trop de coefficients pour la diffusion ?!?' 636 : C% Données incompatibles 636 : CALL ERREUR(22) 637 : RETURN 638 : ENDIF 639 : SEGACT MCHELG 640 : ENDIF 641 : SEGACT MCHEL3 642 : IF(TSOUR)THEN 643 : SEGACT MCHEL4 644 : ENDIF 645 : IF(IKOMP.EQ.2)THEN 646 : SEGACT MCHEL5 647 : ENDIF 648 : 649 : IF(MAX(1,MELEM2.LISOUS(/1)).NE.MAX(1,LISOUS(/1)))THEN 650 : WRITE(6,*)' Geometries incompatibles dans ',nomper 651 : C% Données incompatibles 651 : CALL ERREUR(22) 652 : RETURN 653 : ENDIF 654 : 655 : NKD=0 656 : DO 101 L=1,MAX(1,LISOUS(/1)) 657 : IPT1=MELEME 658 : IPT2=MELEM2 659 : IF(LISOUS(/1).NE.0)IPT1=LISOUS(L) 660 : IF(MELEM2.LISOUS(/1).NE.0)THEN 661 : IPT2=MELEM2.LISOUS(L) 662 : NKD=0 663 : ENDIF 664 : SEGACT IPT1,IPT2 665 : 666 : NOM0 = NOMS(IPT1.ITYPEL)//' ' 666 : CALL KALPBG(NOM0,'FONFORM ',IZFFM) 667 : IF(IZFFM.EQ.0)RETURN 668 : SEGACT IZFFM*MOD 669 : IZHR=KZHR(1) 670 : SEGACT IZHR*MOD 671 : IZF1 = KTP(1) 672 : IZH2 = KZHR(2) 673 : IZW = IZFFM 674 : IZWH = IZHR 675 : 676 : NES=GR(/1) 677 : NPG=GR(/3) 678 : 679 : NP = IPT1.NUM(/1) 680 : MP = IPT2.NUM(/1) 681 : C? MP = IZW.FN(/1) ceci doit etre identique 682 : 683 : NBEL=IPT1.NUM(/2) 684 : 685 : SEGINI SAJT 686 : 687 : c....................................................................... 688 : C MATRIK . MATRIK . MATRIK . MATRIK . MATRIK . MATRIK MATRIK . MATRIK . 689 : c write(6,*)' KKMACO=',KKMACO,'TMATRI=',TMATRI 690 : IF(KKMACO.NE.2)THEN 691 : IF(TMATRI)THEN 692 : C Cas RIGIDITE 693 : IF(XRIG)THEN 694 : IRIGEL(1,NRI+L)=MELEME 695 : COERIG(L)=1.D0 696 : 697 : IRIGEL(7,NRI+L)=0 698 : IF(TCONV)IRIGEL(7,NRI+L)=2 699 : 700 : NBME=1 701 : IF(IHV.EQ.1)NBME=IDIM 702 : NLIGRP=NP 703 : NLIGRD=MP 704 : SEGINI DESCR 705 : IRIGEL(3,NRI+L)=DESCR 706 : IF(NBME.EQ.1)THEN 707 : DO 102 I=1,NLIGRP 708 : LISINC(I)=NOMP(1)//' ' 709 : NOELEP(I)=I 710 : 102 CONTINUE 711 : DO 103 I=1,NLIGRD 712 : LISDUA(I)=NOMD(1)//' ' 713 : NOELED(I)=I 714 : 103 CONTINUE 715 : ELSE 716 : ENDIF 717 : SEGDES DESCR 718 : 719 : NELRIG=NBEL 720 : SEGINI xMATRI 721 : IRIGEL(4,NRI+L)=xMATRI 722 : c write(6,*)'NELRIG,IMATRI=',NELRIG,IMATRI 723 : 724 : * DO 104 K=1,NELRIG 725 : * SEGINI XMATRI 726 : c write(6,*)'NLIGRD,NLIGRP,XMATRI=',NLIGRD,NLIGRP,XMATRI 727 : * IMATTT(K)=XMATRI 728 : * 104 CONTINUE 729 : 730 : ELSE 731 : C Cas MATRIK 732 : SEGINI IPM1 733 : mtail=(IPM1.AM(/1))*(ipm1.am(/2))*(ipm1.am(/3)) 734 : JMATR1=JRIGEL(4,NMATR0+1) 735 : JMATR1.LIZAFM(L,1)=IPM1 736 : IPM2=IPM1 737 : IPM3=IPM1 738 : IF(NBME.GE.2)THEN 739 : JMATR2=JRIGEL(4,NMATR0+2) 740 : IF(IAXI.NE.0)THEN 741 : SEGINI IPM2 742 : mtail=(IPM2.AM(/1))*(ipm2.am(/2))*(ipm2.am(/3)) 743 : JMATR2.LIZAFM(L,1)=IPM2 744 : ICAL=.TRUE. 745 : ELSE 746 : IPM2=IPM1 747 : JMATR2.LIZAFM(L,1)=IPM2 748 : ICAL=.FALSE. 749 : ENDIF 750 : ENDIF 751 : IF(NBME.GE.3)THEN 752 : JMATR3=JRIGEL(4,NMATR0+3) 753 : IPM3=IPM1 754 : JMATR3.LIZAFM(L,1)=IPM3 755 : ICAL=.FALSE. 756 : ENDIF 757 : ENDIF 758 : ENDIF 759 : c write(6,*)' YCLS appelé par ',NOMPER,' Taille matrice=',mtail 760 : 761 : ELSE 762 : 763 : C Cas RIGIDITE 764 : IF(XRIG)THEN 765 : 766 : ELSE 767 : C Cas MATRIK 768 : JMATR1=JRIGEL(4,1) 769 : SEGACT JMATR1 770 : IPM1=JMATR1.LIZAFM(L,1) 771 : IPM2=IPM1 772 : IPM3=IPM1 773 : SEGACT IPM1 774 : IF(NBME.GE.2)THEN 775 : JMATR2=JRIGEL(4,2) 776 : SEGACT JMATR2 777 : IPM2=JMATR2.LIZAFM(L,1) 778 : SEGACT IPM2 779 : ENDIF 780 : IF(NBME.GE.3)THEN 781 : JMATR3=JRIGEL(4,3) 782 : IPM3=JMATR3.LIZAFM(L,1) 783 : SEGACT IPM3 784 : ENDIF 785 : ENDIF 786 : ENDIF 787 : C FIN MATRIK . FIN MATRIK . FIN MATRIK . FIN MATRIK . FIN MATRIK . FIN M 788 : C ...................................................................... 789 : 790 : 791 : C----Ro = 1. 792 : IK1=1 793 : MCHAM1=MCHEL1.ICHAML(L) 794 : SEGACT MCHAM1 795 : 796 : MELVA1=MCHAM1.IELVAL(1) 797 : SEGACT MELVA1 798 : N1PTEL=MELVA1.VELCHE(/1) 799 : N1EL=MELVA1.VELCHE(/2) 800 : IF(N1EL.EQ.1)THEN 801 : IK1=1 802 : ELSEIF(N1EL.EQ.NBEL)THEN 803 : IK1=0 804 : ENDIF 805 : 806 : C----Nu ou Lambda 807 : IK2=1 808 : MCHAM2=MCHEL2.ICHAML(L) 809 : SEGACT MCHAM2 810 : 811 : MELVA2=MCHAM2.IELVAL(1) 812 : SEGACT MELVA2 813 : N1PTEL=MELVA2.VELCHE(/1) 814 : N1EL=MELVA2.VELCHE(/2) 815 : IF(N1EL.EQ.1)THEN 816 : IK2=1 817 : ELSEIF(N1EL.EQ.NBEL)THEN 818 : IK2=0 819 : ENDIF 820 : 821 : IF(MCHELG.NE.0)THEN 822 : MCHAMG=MCHELG.ICHAML(L) 823 : SEGACT MCHAMG 824 : 825 : MELVAG=MCHAMG.IELVAL(1) 826 : SEGACT MELVAG 827 : ENDIF 828 : 829 : C----U 830 : IK3=1 831 : MCHAM3=MCHEL3.ICHAML(L) 832 : SEGACT MCHAM3 833 : 834 : MELVA3=MCHAM3.IELVAL(1) 835 : SEGACT MELVA3 836 : N1PTEL=MELVA3.VELCHE(/1) 837 : N1EL=MELVA3.VELCHE(/2) 838 : IF(N1EL.EQ.1)THEN 839 : IK3=1 840 : ELSEIF(N1EL.EQ.NBEL)THEN 841 : IK3=0 842 : ENDIF 843 : 844 : C----sources 845 : IK4=1 846 : IF(TSOUR)THEN 847 : MCHAM4=MCHEL4.ICHAML(L) 848 : SEGACT MCHAM4 849 : 850 : MELVA4=MCHAM4.IELVAL(1) 851 : SEGACT MELVA4 852 : N1PTEL=MELVA4.VELCHE(/1) 853 : N1EL=MELVA4.VELCHE(/2) 854 : IF(N1EL.EQ.1)THEN 855 : IK4=1 856 : ELSEIF(N1EL.EQ.NBEL)THEN 857 : IK4=0 858 : ENDIF 859 : ENDIF 860 : 861 : C----Div (rocp u ) 862 : IK5=1 863 : IF(IKOMP.EQ.2)THEN 864 : MCHAM5=MCHEL5.ICHAML(L) 865 : SEGACT MCHAM5 866 : 867 : MELVA5=MCHAM5.IELVAL(1) 868 : SEGACT MELVA5 869 : N1PTEL=MELVA5.VELCHE(/1) 870 : N1EL=MELVA5.VELCHE(/2) 871 : IF(N1EL.EQ.1)THEN 872 : IK5=1 873 : ELSEIF(N1EL.EQ.NBEL)THEN 874 : IK5=0 875 : ENDIF 876 : ENDIF 877 : 878 : c write(6,*)' AVANT 108 NC=',NC,' NBEL=',NBEL,MP,NP,NC 879 : C=============================================== 880 : AI1=AIMPL-1.D0 881 : IF(AIMPL.EQ.0.D0)THEN 882 : AI2=-1.D0 883 : ELSE 884 : AI2=AI1/AIMPL 885 : ENDIF 886 : 887 : DO 108 KE=1,NBEL 888 : 889 : NK1=KE + IK1*(1 - KE) 890 : NK2=KE + IK2*(1 - KE) 891 : NK3=KE + IK3*(1 - KE) 892 : NK4=KE + IK4*(1 - KE) 893 : NK5=KE + IK5*(1 - KE) 894 : 895 : DO 109 I=1,NP 896 : J=IPT1.NUM(I,KE) 897 : DO 109 N=1,IDIM 898 : XYZ(N,I)=XCOOR((J-1)*(IDIM+1)+N) 899 : 109 CONTINUE 900 : 901 : DO 107 I=1,NP 902 : I1=MLENT1.LECT(IPT1.NUM(I,KE)) 903 : DO 107 N=1,NC 904 : TN1(I,N)=MTETA1.VPOCHA(I1,N) 905 : 107 CONTINUE 906 : 906 : CALL CALJBR(FN,GR,PG,XYZ,HR,PGSQ,RPG,NES, 907 : * IDIM,NP,NPG,IAXI,AIRE,AJ,SGN) 908 : 909 : IF (IDCEN.EQ.0.OR.IDCEN.EQ.1) THEN 909 : CALL RSETD(WT,FN,(NP*NPG)) 909 : CALL RSETD(WS,FN,(NP*NPG)) 910 : ELSE 910 : CALL CALDEC(WT,WS,XYZ,GR,HR,FN,NES,IDIM,NP,NPG,AJT, 911 : & IDCEN,CMD,MELVA1.VELCHE(1,NK1),MELVA2.VELCHE(1,NK2), 912 : & MELVA3.VELCHE(1,NK3),TN1,NC,IKOMP,XREF,AIRE,KE) 913 : ENDIF 914 : 914 : CALL INITD(SM1,(MP*IDIM),0.D0) 915 : 916 : C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: 917 : IF(KKMACO.EQ.2)THEN 918 : 919 : C ...... Chargement Rigidite ou Matrik 920 : C Cas RIGIDITE 921 : IF(XRIG)THEN 922 : * XMATRI=IMATTT(KE) 923 : DO 911 I=1,MP 924 : DO 911 J=1,NP 925 : RE(I,J,ke)=RF1(J,I,1) 926 : 911 CONTINUE 927 : * SEGDES XMATRI 928 : ELSE 929 : C Cas MATRIK 930 : DO 912 N=1,NBME 931 : JMATR1=JRIGEL(4,N) 932 : IPM4=JMATR1.LIZAFM(L,1) 933 : DO 912 I=1,MP 934 : DO 912 J=1,NP 935 : RF1(J,I,N)=IPM4.AM(KE,J,I) 936 : 912 CONTINUE 937 : ENDIF 938 : C...... Multiplication Matrice Vecteur 939 : 940 : DO 710 I=1,MP 941 : DO 711 J=1,NP 942 : 943 : DO 716 N=1,NC 944 : SM1(I,N) = SM1(I,N) + AI2*RF1(J,I,N)*TN1(J,N) 945 : 716 CONTINUE 946 : 947 : 711 CONTINUE 948 : 710 CONTINUE 949 : 950 : ELSEIF(KKMACO.NE.2)THEN 951 : C................................................................. 952 : C...... Convection Laplacien .... préparation forme conservative 953 : C...... Convection Laplacien .... et non conservative 954 : C préparation 955 : AKOMP=1.D0 956 : IF(IKOMP.EQ.1)AKOMP=-1.D0 957 : NIPG=NPG*IDIM 958 : I1=0 959 : DO 420 LG=1,NPG 960 : PDR=PGSQ(LG)*DEUPI*RPG(LG) 961 : C1=C1CNG*MELVA1.VELCHE(LG,NK1)*PDR*AKOMP 962 : C2=C2LAPN*MELVA2.VELCHE(LG,NK2)*PDR 963 : DO 420 N=1,IDIM 964 : I1=I1+1 965 : C3=MELVA3.VELCHE((N-1)*NPG+LG,NK3)*C1 966 : DO 420 I=1,MP 967 : WTC(I1,I)=WT(I,LG)*C3 + IZWH.HR(N,I,LG)*C2 968 : 420 CONTINUE 969 : 970 : I1=0 971 : DO 421 LG=1,NPG 972 : DO 421 N=1,IDIM 973 : I1=I1+1 974 : DO 421 J=1,NP 975 : HRD(I1,J)=HR(N,J,LG) 976 : 421 CONTINUE 977 : 978 : IF(IKOMP.EQ.0.OR.IKOMP.EQ.2)THEN 979 : C...... Convection .... forme non conservative 980 : DO 401 I=1,MP 981 : DO 401 J=1,NP 982 : U3=0.D0 983 : DO 402 I1=1,NIPG 984 : U3=U3+(WTC(I1,I) * HRD(I1,J)) 985 : 402 CONTINUE 986 : 987 : DO 403 N=1,NC 988 : SM1(I,N) = SM1(I,N) + AI1*U3*TN1(J,N) 989 : 403 CONTINUE 990 : RF1(J,I,1) = AIEX*U3 991 : 992 : 401 CONTINUE 993 : ENDIF 994 : 995 : IF(IKOMP.EQ.1)THEN 996 : C...... Convection .... forme conservative 997 : DO 411 I=1,MP 998 : DO 411 J=1,NP 999 : U3=0.D0 1000 : DO 412 I1=1,NIPG 1001 : U3=U3+(WTC(I1,I) * HRD(I1,J)) 1002 : 412 CONTINUE 1003 : 1004 : DO 413 N=1,NC 1005 : SM1(J,N) = SM1(J,N) + AI1*U3*TN1(I,N) 1006 : 413 CONTINUE 1007 : RF1(I,J,1) = AIEX*U3 1008 : 1009 : 411 CONTINUE 1010 : ENDIF 1011 : 1012 : C...... Fin Convection .... 1013 : 1014 : C.......... Cas KMU=1 et FTAU pour QDM (Formulation en To) 1015 : IF(KMU.EQ.1)THEN 1016 : c write(6,*)' Cas FTAU' 1017 : 1018 : DO 4162 LG=1,NPG 1019 : PDR=PGSQ(LG)*DEUPI*RPG(LG) 1020 : C2=C2LAPN*MELVA2.VELCHE(LG,NK2)*PDR 1021 : DO 4162 N=1,IDIM 1022 : DO 4162 J=1,NP 1023 : FNE(LG,J,N)=HR(N,J,LG)*C2 1024 : HRE(LG,J,N)=HR(N,J,LG) 1025 : 4162 CONTINUE 1026 : 1027 : IF(IDIM.EQ.2)THEN 1028 : DO 4192 I=1,NP 1029 : US1=0.D0 1030 : US2=0.D0 1031 : DO 4182 J=1,NP 1032 : U11=0.D0 1033 : U12=0.D0 1034 : U21=0.D0 1035 : U22=0.D0 1036 : DO 4172 LG=1,NPG 1037 : U11=U11+HRE(LG,J,1)*FNE(LG,I,1) 1038 : U12=U12+HRE(LG,J,1)*FNE(LG,I,2) 1039 : U21=U21+HRE(LG,J,2)*FNE(LG,I,1) 1040 : U22=U22+HRE(LG,J,2)*FNE(LG,I,2) 1041 : 4172 CONTINUE 1042 : US1 = US1 - U11*TN1(J,1) - U12*TN1(J,2) 1043 : US2 = US2 - U21*TN1(J,1) - U22*TN1(J,2) 1044 : 4182 CONTINUE 1045 : SM1(I,1) = SM1(I,1) + US1 1046 : SM1(I,2) = SM1(I,2) + US2 1047 : 4192 CONTINUE 1048 : ELSE 1049 : DO 5193 I=1,NP 1050 : US1=0.D0 1051 : US2=0.D0 1052 : US3=0.D0 1053 : DO 5183 J=1,NP 1054 : U11=0.D0 1055 : U12=0.D0 1056 : U13=0.D0 1057 : U21=0.D0 1058 : U22=0.D0 1059 : U23=0.D0 1060 : U31=0.D0 1061 : U32=0.D0 1062 : U33=0.D0 1063 : DO 5173 LG=1,NPG 1064 : U11=U11+HRE(LG,J,1)*FNE(LG,I,1) 1065 : U12=U12+HRE(LG,J,1)*FNE(LG,I,2) 1066 : U13=U13+HRE(LG,J,1)*FNE(LG,I,3) 1067 : 1068 : U21=U21+HRE(LG,J,2)*FNE(LG,I,1) 1069 : U22=U22+HRE(LG,J,2)*FNE(LG,I,2) 1070 : U23=U23+HRE(LG,J,2)*FNE(LG,I,3) 1071 : 1072 : U31=U31+HRE(LG,J,3)*FNE(LG,I,1) 1073 : U32=U32+HRE(LG,J,3)*FNE(LG,I,2) 1074 : U33=U33+HRE(LG,J,3)*FNE(LG,I,3) 1075 : 5173 CONTINUE 1076 : US1 = US1 - U11*TN1(J,1) - U12*TN1(J,2) - U13*TN1(J,3) 1077 : US2 = US2 - U21*TN1(J,1) - U22*TN1(J,2) - U23*TN1(J,3) 1078 : US3 = US3 - U31*TN1(J,1) - U32*TN1(J,2) - U33*TN1(J,3) 1079 : 5183 CONTINUE 1080 : SM1(I,1) = SM1(I,1) + US1 1081 : SM1(I,2) = SM1(I,2) + US2 1082 : SM1(I,3) = SM1(I,3) + US3 1083 : 5193 CONTINUE 1084 : ENDIF 1085 : 1086 : 1087 : ENDIF 1088 : C...... FIN Cas KMU=1 et FTAU pour QDM (Formulation en To) 1089 : 1090 : C.......... Cas KMU=2 et MUVARI pour QDM (Formulation en Grad de Mu) 1091 : IF(KMU.EQ.2.AND.MCHELG.NE.0)THEN 1092 : c write(6,*)' Cas MUVARI' 1093 : 1093 : CALL INITD(GMU,(NPG*IDIM),0.D0) 1094 : DO 415 LG=1,NPG 1095 : DO 415 I=1,NP 1096 : DO 415 N=1,NC 1097 : GMU(LG,N)=GMU(LG,N)+HR(N,I,LG)*MELVAG.VELCHE(I,NK2) 1098 : 415 CONTINUE 1099 : 1100 : 1101 : DO 416 LG=1,NPG 1102 : PDR=PGSQ(LG)*DEUPI*RPG(LG) 1103 : DO 416 N=1,IDIM 1104 : DO 416 J=1,NP 1105 : FNE(LG,J,N)=FN(J,LG)*GMU(LG,N) 1106 : HRE(LG,J,N)=HR(N,J,LG)*PDR 1107 : 416 CONTINUE 1108 : 1109 : IF(IDIM.EQ.2)THEN 1110 : DO 419 I=1,NP 1111 : US1=0.D0 1112 : US2=0.D0 1113 : DO 418 J=1,NP 1114 : U11=0.D0 1115 : U12=0.D0 1116 : U21=0.D0 1117 : U22=0.D0 1118 : DO 417 LG=1,NPG 1119 : U11=U11+HRE(LG,J,1)*FNE(LG,I,1) 1120 : U12=U12+HRE(LG,J,1)*FNE(LG,I,2) 1121 : U21=U21+HRE(LG,J,2)*FNE(LG,I,1) 1122 : U22=U22+HRE(LG,J,2)*FNE(LG,I,2) 1123 : 417 CONTINUE 1124 : US1 = US1 - U11*TN1(J,1) - U12*TN1(J,2) 1125 : US2 = US2 - U21*TN1(J,1) - U22*TN1(J,2) 1126 : 418 CONTINUE 1127 : SM1(I,1) = SM1(I,1) + US1 1128 : SM1(I,2) = SM1(I,2) + US2 1129 : 419 CONTINUE 1130 : ELSE 1131 : DO 4193 I=1,NP 1132 : US1=0.D0 1133 : US2=0.D0 1134 : US3=0.D0 1135 : DO 4183 J=1,NP 1136 : U11=0.D0 1137 : U12=0.D0 1138 : U13=0.D0 1139 : U21=0.D0 1140 : U22=0.D0 1141 : U23=0.D0 1142 : U31=0.D0 1143 : U32=0.D0 1144 : U33=0.D0 1145 : DO 4173 LG=1,NPG 1146 : c U11=U11+FN(I,LG)*HRE(LG,J,1)*GMU(LG,1) 1147 : c U12=U12+FN(I,LG)*HRE(LG,J,1)*GMU(LG,2) 1148 : c U21=U21+FN(I,LG)*HRE(LG,J,2)*GMU(LG,1) 1149 : c U22=U22+FN(I,LG)*HRE(LG,J,2)*GMU(LG,2) 1150 : 1151 : U11=U11+HRE(LG,J,1)*FNE(LG,I,1) 1152 : U12=U12+HRE(LG,J,1)*FNE(LG,I,2) 1153 : U13=U13+HRE(LG,J,1)*FNE(LG,I,3) 1154 : 1155 : U21=U21+HRE(LG,J,2)*FNE(LG,I,1) 1156 : U22=U22+HRE(LG,J,2)*FNE(LG,I,2) 1157 : U23=U23+HRE(LG,J,2)*FNE(LG,I,3) 1158 : 1159 : U31=U31+HRE(LG,J,3)*FNE(LG,I,1) 1160 : U32=U32+HRE(LG,J,3)*FNE(LG,I,2) 1161 : U33=U33+HRE(LG,J,3)*FNE(LG,I,3) 1162 : 4173 CONTINUE 1163 : US1 = US1 - U11*TN1(J,1) - U12*TN1(J,2) - U13*TN1(J,3) 1164 : US2 = US2 - U21*TN1(J,1) - U22*TN1(J,2) - U23*TN1(J,3) 1165 : US3 = US3 - U31*TN1(J,1) - U32*TN1(J,2) - U33*TN1(J,3) 1166 : 4183 CONTINUE 1167 : SM1(I,1) = SM1(I,1) + US1 1168 : SM1(I,2) = SM1(I,2) + US2 1169 : SM1(I,3) = SM1(I,3) + US3 1170 : 4193 CONTINUE 1171 : ENDIF 1172 : 1173 : 1174 : ENDIF 1175 : C...... FIN Cas KMU=2 et MUVARI pour QDM (Formulation en Grad de Mu) 1176 : 1177 : C.......... Cas CNG uniquement (complément) 1178 : IF(IDCEN.EQ.5)THEN 1179 : C préparation 1180 : I1=0 1181 : DO 520 LG=1,NPG 1182 : PDR=PGSQ(LG)*DEUPI*RPG(LG) 1183 : C1=C1CNG*MELVA1.VELCHE(LG,NK1)*PDR 1184 : DO 520 N=1,IDIM 1185 : I1=I1+1 1186 : C3=MELVA3.VELCHE((N-1)*NPG+LG,NK3)*C1 1187 : DO 520 I=1,MP 1188 : WSC(I1,I)=WS(I,LG)*C3 1189 : 520 CONTINUE 1190 : 1191 : DO 510 I=1,MP 1192 : DO 511 J=1,NP 1193 : U3=0.D0 1194 : DO 512 I1=1,NIPG 1195 : U3=U3+WSC(I1,I) * HRD(I1,J) 1196 : 512 CONTINUE 1197 : 1198 : DO 516 N=1,NC 1199 : SM1(I,N) = SM1(I,N) - U3*TN1(J,N) 1200 : 516 CONTINUE 1201 : 1202 : 511 CONTINUE 1203 : 510 CONTINUE 1204 : ENDIF 1205 : C..........Fin Cas CNG uniquement (complément) 1206 : 1207 : C.............................................................. 1208 : C.............................................................. 1209 : IF(IKOMP.EQ.2)THEN 1210 : C...... Convection Laplacien .... Formulation conservative 1211 : C préparation 1212 : DO 430 LG=1,NPG 1213 : PDR=PGSQ(LG)*DEUPI*RPG(LG) 1214 : C5=C1CNG*MELVA5.VELCHE(LG,NK5)*PDR 1215 : DO 430 I=1,MP 1216 : WTC(LG,I)=WT(I,LG)*C5 1217 : 430 CONTINUE 1218 : 1219 : DO 431 LG=1,NPG 1220 : DO 431 J=1,NP 1221 : FND(LG,J)=FN(J,LG) 1222 : 431 CONTINUE 1223 : 1224 : DO 440 I=1,MP 1225 : DO 441 J=1,NP 1226 : U3=0.D0 1227 : DO 442 LG=1,NPG 1228 : U3=U3+WTC(LG,I) * FND(LG,J) 1229 : 442 CONTINUE 1230 : 1231 : DO 446 N=1,NC 1232 : SM1(I,N) = SM1(I,N) + AI1*U3*TN1(J,N) 1233 : 446 CONTINUE 1234 : RF1(J,I,1)=RF1(J,I,1)+(AIEX*U3) 1235 : 1236 : 441 CONTINUE 1237 : 440 CONTINUE 1238 : 1239 : C.......... Cas CNG uniquement (complément) 1240 : IF(IDCEN.EQ.5)THEN 1241 : C préparation 1242 : I1=0 1243 : DO 530 LG=1,NPG 1244 : PDR=PGSQ(LG)*DEUPI*RPG(LG) 1245 : C1=C1CNG*MELVA1.VELCHE(LG,NK1)*PDR 1246 : DO 530 N=1,IDIM 1247 : I1=I1+1 1248 : C3=MELVA3.VELCHE((N-1)*NPG+LG,NK3)*C1 1249 : DO 530 I=1,MP 1250 : WSC(I1,I)=WS(I,LG)*C3 1251 : 530 CONTINUE 1252 : 1253 : DO 540 I=1,MP 1254 : DO 541 J=1,NP 1255 : U3=0.D0 1256 : DO 542 I1=1,NIPG 1257 : U3=U3+WSC(I1,I) * HRD(I1,J) 1258 : 542 CONTINUE 1259 : 1260 : DO 546 N=1,NC 1261 : SM1(I,N) = SM1(I,N) - U3*TN1(J,N) 1262 : 546 CONTINUE 1263 : 1264 : 541 CONTINUE 1265 : 540 CONTINUE 1266 : ENDIF 1267 : C..........Fin Cas CNG uniquement (complément) 1268 : ENDIF 1269 : C..........Fin Complément pour IKOMP=2 ...................... 1270 : 1271 : C======================================================================= 1272 : c Cas 2D axi pour la QDM 1273 : IF(IHV.EQ.1.AND.IAXI.EQ.2)THEN 1274 : DO 320 LG=1,NPG 1275 : PDR=PGSQ(LG)*DEUPI/RPG(LG) 1276 : C2=C2LAPN*MELVA2.VELCHE(LG,NK2)*PDR 1277 : * pv semble inutile et faux C3=MELVA3.VELCHE((N-1)*NPG+LG,NK3) 1278 : DO 320 I=1,MP 1279 : WTC(LG,I)=WT(I,LG)*C2 1280 : 320 CONTINUE 1281 : 1282 : DO 310 I=1,MP 1283 : DO 311 J=1,NP 1284 : UR=0.D0 1285 : DO 312 LG=1,NPG 1286 : UR=UR+WTC(LG,I)*FN(J,LG) 1287 : 312 CONTINUE 1288 : SM1(I,1) = SM1(I,1) + AI1*UR*TN1(J,1) 1289 : RF1(J,I,2) = RF1(J,I,1) 1290 : RF1(J,I,1)=RF1(J,I,1)+(AIEX*UR) 1291 : 311 CONTINUE 1292 : 310 CONTINUE 1293 : ENDIF 1294 : c Cas 2D axi pour la QDM Fin 1295 : C======================================================================= 1296 : 1297 : ENDIF 1298 : 1299 : C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: 1300 : IF(TSOUR) THEN 1301 : C...... Source 1302 : DO 610 I=1,MP 1303 : DO 617 N=1,NC 1304 : U4=0.D0 1305 : DO 615 LG=1,NPG 1306 : C4=MELVA4.VELCHE((N-1)*NPG+LG,NK4) 1307 : U4=U4+WT(I,LG)*PGSQ(LG)*C4*DEUPI*RPG(LG) 1308 : 615 CONTINUE 1309 : SM1(I,N)=SM1(I,N)+ U4 1310 : 617 CONTINUE 1311 : 610 CONTINUE 1312 : C...... Source Fin 1313 : ENDIF 1314 : C======================================================================= 1315 : 1316 : C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 1317 : IF(TMATRI.AND.KKMACO.NE.2)THEN 1318 : C ...... Chargement Rigidite ou Matrik 1319 : C Cas RIGIDITE 1320 : IF(XRIG)THEN 1321 : * XMATRI=IMATTT(KE) 1322 : DO 901 I=1,MP 1323 : DO 901 J=1,NP 1324 : RE(I,J,ke)=RF1(J,I,1) 1325 : 901 CONTINUE 1326 : * SEGDES XMATRI 1327 : ELSE 1328 : C Cas MATRIK 1329 : NBMM=1 1330 : IF(ICAL)NBMM=NBME 1331 : DO 902 N=1,NBMM 1332 : JMATR1=JRIGEL(4,NMATR0+N) 1333 : IPM4=JMATR1.LIZAFM(L,1) 1334 : DO 902 I=1,NP 1335 : DO 902 J=1,NP 1336 : IPM4.AM(KE,J,I)=RF1(J,I,N) 1337 : 902 CONTINUE 1338 : ENDIF 1339 : ENDIF 1340 : 1341 : C ...... Chargement Second membre 1342 : DO 910 I=1,NP 1343 : I1=LECT(IPT1.NUM(I,KE)) 1344 : DO 910 N=1,NC 1345 : MPOVA1.VPOCHA(I1,N)=MPOVA1.VPOCHA(I1,N)+SM1(I,N) 1346 : 910 CONTINUE 1347 : C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 1348 : 1349 : 108 CONTINUE 1350 : 1351 : SEGSUP MCHAM1,MELVA1 1352 : SEGSUP MCHAM2,MELVA2 1353 : SEGSUP MCHAM3,MELVA3 1354 : 1355 : IF(TSOUR)THEN 1356 : SEGSUP MCHAM4,MELVA4 1357 : ENDIF 1358 : 1359 : IF(IKOMP.EQ.2)THEN 1360 : SEGSUP MCHAM5,MELVA5 1361 : ENDIF 1362 : 1363 : 1364 : 1365 : SEGDES IPT1,IPT2 1366 : 1367 : IF(TMATRI)THEN 1368 : C Cas RIGIDITE 1369 : IF(XRIG)THEN 1370 : SEGDES xMATRI 1371 : ELSE 1372 : C Cas MATRIK 1373 : SEGDES IPM1 1374 : IF(NBME.GE.2)SEGDES IPM2 1375 : IF(NBME.GE.3)SEGDES IPM3 1376 : ENDIF 1377 : ENDIF 1378 : 1379 : SEGSUP IZFFM,IZHR,IZF1,IZH2 1380 : SEGSUP SAJT 1381 : 1382 : 101 CONTINUE 1383 : 1384 : IF(TMATRI)THEN 1385 : IF(.NOT.XRIG)THEN 1386 : NBMM=JRIGEL(/2) 1387 : IF(NBMM.NE.0)THEN 1388 : DO 141 M=1,NBMM 1389 : JMATRI=JRIGEL(4,M) 1390 : SEGDES JMATRI 1391 : 141 CONTINUE 1392 : ENDIF 1393 : ENDIF 1394 : ENDIF 1395 : 1396 : C""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" 1397 : C""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" 1398 : C""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" 1399 : 1400 : 1401 : SEGSUP MCHEL1 1402 : SEGSUP MCHEL2 1403 : SEGSUP MCHEL3 1404 : 1405 : IF(TSOUR)THEN 1406 : SEGSUP MCHEL4 1407 : ENDIF 1408 : 1409 : IF(IKOMP.EQ.2)THEN 1410 : SEGSUP MCHEL5 1411 : ENDIF 1412 : 1413 : SEGDES MCHPO1,MPOVA1 1414 : SEGDES MELEME,MELEM2 1415 : 1416 : SEGSUP MLENTI 1417 : c IF(TLINCO)THEN 1418 : SEGSUP MLENT1 1419 : SEGDES MTETA1 1420 : c ENDIF 1421 : 1422 : IF(TMATRI)THEN 1423 : C Cas RIGIDITE 1424 : IF(XRIG)THEN 1425 : SEGDES MRIGID 1426 : ELSE 1427 : C Cas MATRIK 1428 : SEGDES MATRIK 1429 : ENDIF 1430 : ENDIF 1431 : 1432 : 1433 : c segact MCHPO1 1434 : c MSOUP1=MCHPO1.IPCHP(1) 1435 : c segact MSOUP1 1436 : c write(6,*)MCHPO1.IPCHP(/1),MSOUP1.NOCOMP(/2) 1437 : c write(6,*)(MSOUP1.NOCOMP(i),i=1,MSOUP1.NOCOMP(/2)) 1437 : c if(ihv.eq.1)call ecmo(mtabz,'TOTO','CHPOINT',MCHPO1) 1438 : 1439 : C************************************************************************* 1440 : 1441 : IF(KKMACO.EQ.1)THEN 1442 : TYPE='MATRIK' 1442 : CALL ECMO(MTAB1,NMACO,TYPE,MATRIK) 1443 : ENDIF 1444 : 1445 : IF(AIMPL.EQ.0.D0)THEN 1446 : NRIGE=7 1447 : NKID =9 1448 : NKMT =7 1449 : NMATRI=0 1450 : SEGINI MATRIK 1451 : ENDIF 1452 : 1453 : c write(6,*)' FIN YCLS' 1454 : RETURN 1455 : 1001 FORMAT(20(1X,I5)) 1456 : 1002 FORMAT(10(1X,1PE11.4)) 1457 : END 1458 : 1459 : 1460 : 1461 : 1462 : 1463 : 1464 : 1465 : 1466 : 1467 : 1468 : 1469 : 1470 : 1471 : 1472 : 1473 : 1474 : 1475 : 1476 : 1477 : 1478 : 1479 : 1480 : 1481 : 1482 : 1483 : 1484 : 1485 : 1486 : 1487 : 1488 : 1489 :
© Cast3M 2003 - All rights reserved.
Disclaimer