1 : C YFRT SOURCE CHAT 11/03/16 21:33:55 6902 2 : SUBROUTINE YFRT(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) 66 : ENDSEGMENT 67 : -INC SMRIGID 68 : C-INC SMMATRIK 69 : C******************************************************************* 70 : C 71 : SEGMENT MATRIK 72 : REAL*8 COEMTK(NMATRI) 73 : INTEGER JRIGEL(NRIGE,NMATRI) 74 : INTEGER KSYM,KMINC,KMINCP,KMINCD,KIZM 75 : INTEGER KISPGT,KISPGP,KISPGD 76 : INTEGER KNTTT,KNTTP,KNTTD 77 : INTEGER KIDMAT(NKID) 78 : INTEGER KKMMT(NKMT) 79 : ENDSEGMENT 80 : 81 : SEGMENT JMATRI 82 : CHARACTER*8 LISPRI(NBMF),LJSDUA(NBMF) 83 : INTEGER LIZAFM(NBSOUS,NBMF) 84 : INTEGER KSPGP,KSPGD 85 : ENDSEGMENT 86 : POINTEUR JMATRS.JMATRI,JMATR1.JMATRI,JMATR2.JMATRI,JMATR3.JMATRI 87 : 88 : C Stokage matrices elementaires non assemblees (valeurs) 89 : SEGMENT IZAFM 90 : REAL*8 AM(NBEL,NP,MP) 91 : ENDSEGMENT 92 : POINTEUR IPM1.IZAFM,IPM2.IZAFM,IPM3.IZAFM,IPM4.IZAFM 93 : 94 : C******************************************************************* 95 : POINTEUR IPM.IZAFM 96 : -INC SMCHAML 97 : 98 : LOGICAL TDFDT,TCONV,TSOUR,TSOURB,TMATRI 99 : LOGICAL ICAL,XPG,XRIG 100 : LOGICAL XDIAG,XTV,XTG,XBDF,XCONS 101 : CHARACTER*8 CHHH,TYPE,NOM,NOM0,CHAI,SCHT,NOMPER,NOMPR,MPRE 102 : CHARACTER*4 INCOD 103 : CHARACTER*4 NOMI(3),NOMA(3),NOMP(3),NOMD(3),NMACO 104 : DIMENSION ICOEF(4),SU(3),XPOI(3) 105 : C***************************************************************************** 106 : CTCLSF 107 : NOMPER=NOMPR 108 : write(6,*)' DEBUT YFRT appele par ',NOMPER 109 : 110 : C- Table de l'opérateur (pointeur MTABX) 111 : 112 : C- Récupération de la table RV (pointeur MTAB1) 112 : CALL LEKTAB(MTABX,'EQEX',MTAB1) 113 : C- Récupération de la table des options de l'opérateur (pointeur KOPTI) 113 : CALL LEKTAB(MTABX,'KOPT',KOPTI) 114 : C- Récupération de la table DOMAINE de la zone (pointeur MTABZ) 114 : CALL LEKTAB(MTABX,'DOMZ',MTABZ) 115 : C- Récupération de la table INCO (pointeur KINC) 115 : CALL LEKTAB(MTAB1,'INCO',KINC) 116 : IF(IERR.NE.0) RETURN 117 : 118 : C***************************************************************************** 119 : C Traitement des options 120 : C KIMPL = 0 -> EXPL 1 -> IMPL 2 -> CN 121 : C AIMPL=0 explicite 122 : C KKALE=1 Indice indiquant que l'on est en ALE 123 : C sinon KKALE=0 124 : C KFORM = 0 -> SI 1 -> EF 2 -> VF 3 -> EFMC 125 : C IDCEN : entier indiquant le type de décentrement souhaité 126 : C IDCEN 1-> CENTREE 2-> SUPGDC 3-> SUPG 4-> TVISQUEU 5-> CNG 127 : C E/ CMD : Coefficient multiplicateur du décentrement 128 : C Si IDCEN=4 ou =5 CMD=DT 129 : 129 : CALL ACME(KOPTI,'KFORM',KFORM) 129 : CALL ACME(KOPTI,'KIMPL',KIMPL) 129 : CALL ACMF(KOPTI,'AIMPL',AIMPL) 130 : IKOMP=0 130 : *? CALL ACME(KOPTI,'IKOMP',IKOMP) 130 : CALL ACME(KOPTI,'KMACO',KMACO) 131 : KKMACO=KMACO 131 : CALL ACME(KOPTI,'IDCEN',IDCEN) 131 : CALL ACMF(KOPTI,'CMD',CMD) 131 : *? CALL ACME(KOPTI,'ALE',KKALE) 132 : 133 : C Restrictions des options 134 : IF(NOMPER.NE.'NS')KKALE=0 135 : IF(NOMPER.EQ.'LAPN')IDCEN=1 136 : 137 : IF(KIMPL.EQ.0)AIMPL=0.D0 138 : IF(KIMPL.EQ.1)AIMPL=1.D0 139 : 140 : c write(6,*)' YFRT KMACO=',KMACO,'IDCEN=',IDCEN 141 : IF(KMACO.EQ.1.AND.IDCEN.NE.5)THEN 141 : CALL ACMM(KOPTI,'NMACO',NMACO) 142 : TYPE=' ' 142 : CALL ACMO(MTAB1,NMACO,TYPE,MATRIK) 143 : IF(TYPE.EQ.'MATRIK')THEN 144 : KKMACO=2 145 : ELSE 146 : KKMACO=1 147 : ENDIF 148 : ENDIF 149 : IF(IDCEN.EQ.5)KKMACO=0 150 : c write(6,*)' YFRT KMACO=',kmaco,' KKMACO=',KKMACO,'IKOMP=',ikomp 151 : 152 : C***************************************************************************** 153 : 154 : KPOIND=0 155 : NBMF=1 156 : XRIG=.FALSE. 157 : TMATRI=.TRUE. 158 : XPG=.FALSE. 159 : IF(IDCEN.GT.1)XPG=.TRUE. 160 : IF(IDCEN.EQ.4)AIMPL=0.D0 161 : IF(IDCEN.EQ.5)THEN 162 : AIMPL=1.D0 163 : ENDIF 164 : 165 : AIEX=AIMPL 166 : IF(AIMPL.EQ.0.D0.AND.KKMACO.NE.0)THEN 167 : TMATRI=.TRUE. 168 : AIEX=1.D0 169 : ENDIF 170 : 171 : IF(AIMPL.EQ.0.D0.AND.KKMACO.EQ.0)THEN 172 : TMATRI=.FALSE. 173 : ENDIF 174 : 175 : IAXI=0 176 : IF(IFOMOD.EQ.0)IAXI=2 177 : DEUPI=1.D0 178 : IF(IAXI.NE.0)DEUPI=2.D0*XPI 179 : 180 : IF(XRIG)THEN 181 : IF(IHV.EQ.0)THEN 182 : NOMP(1)='T ' 183 : NOMD(1)='Q ' 184 : ELSEIF(IHV.EQ.1)THEN 185 : NOMP(1)='UX ' 186 : NOMP(2)='UY ' 187 : NOMP(3)='UZ ' 188 : NOMD(1)='FX ' 189 : NOMD(2)='FY ' 190 : NOMD(3)='FZ ' 191 : ENDIF 192 : ELSE 193 : IF(IHV.EQ.0)THEN 194 : NOMP(1)=NOMI(1) 195 : NOMD(1)=NOMI(1) 196 : ELSEIF(IHV.EQ.1)THEN 197 : WRITE(NOMP(1),FMT='(I1,A3)')1,NOMI(1)(1:3) 198 : WRITE(NOMD(1),FMT='(I1,A3)')1,NOMI(1)(1:3) 199 : WRITE(NOMP(2),FMT='(I1,A3)')2,NOMI(1)(1:3) 200 : WRITE(NOMD(2),FMT='(I1,A3)')2,NOMI(1)(1:3) 201 : WRITE(NOMP(3),FMT='(I1,A3)')3,NOMI(1)(1:3) 202 : WRITE(NOMD(3),FMT='(I1,A3)')3,NOMI(1)(1:3) 203 : ENDIF 204 : ENDIF 205 : c write(6,*)' NOMP=',nomp 206 : c write(6,*)' NOMD=',nomd 207 : 207 : CALL LEKTAB(MTABZ,'SOMMET',MELEMS) 207 : CALL LEKTAB(MTABZ,'MAILLAGE',MELEME) 208 : MELEM2=MELEME 209 : 210 : C On cree un second membre non vide 211 : SEGACT MELEMS 212 : N=MELEMS.NUM(/2) 213 : IF(IHV.EQ.0)NC=1 214 : IF(IHV.EQ.1)NC=IDIM 215 : SEGDES MELEMS 216 : NSOUPO=1 217 : NAT=2 218 : SEGINI MCHPO1,MSOUP1,MPOVA1 219 : MCHPO1.JATTRI(1)=2 220 : MCHPO1.IFOPOI=IFOMOD 221 : MCHPO1.MTYPOI=' ' 222 : MCHPO1.MOCHDE=' ' 223 : MCHPO1.IPCHP(1)=MSOUP1 224 : MSOUP1.IGEOC=MELEMS 225 : MSOUP1.IPOVAL=MPOVA1 226 : IF(IHV.EQ.0)THEN 227 : MSOUP1.NOCOMP(1)=NOMD(1) 228 : ELSE 229 : DO 91 N=1,NC 230 : MSOUP1.NOCOMP(N)=NOMD(N) 231 : 91 CONTINUE 232 : ENDIF 233 : SEGDES MSOUP1,MCHPO1 234 : 235 : C************************************************************************* 236 : C Lecture des Inco(s) aux temps precedents si Transitoire 237 : c write(6,*)'Lecture des INCO MZIN=',MZIN 238 : 238 : CALL LRCHT(MZIN,MTETA1,TYPE,IGEOM) 238 : CALL KRIPAD(IGEOM,MLENT1) 238 : CALL VERPAD(MLENT1,MELEMS,IRET) 239 : IF(IRET.NE.0)THEN 240 : C% Indice %m1:8 : L'objet %m9:16 n'a pas le bon support géométrique 241 : MOTERR(1: 8) = 'TETA1' 242 : MOTERR(9:16) = 'CHPOINT ' 243 : WRITE(IOIMP,*)'Opérateur : ',NOMPER 243 : CALL ERREUR(788) 244 : RETURN 245 : ENDIF 246 : 247 : NC=MTETA1.VPOCHA(/2) 248 : 249 : IF((IHV.EQ.0.AND.NC.NE.1).OR. 250 : & (IHV.EQ.1.AND.NC.NE.IDIM))THEN 251 : C% Indice %m1:8 : L'objet %m9:16 n'a pas le bon nombre de composantes 252 : MOTERR( 1: 8) = 'Inconnue' 253 : MOTERR( 9:16) = 'CHPOINT' 254 : WRITE(IOIMP,*)'Operateur : ',NOMPER 254 : CALL ERREUR(784) 255 : RETURN 256 : ENDIF 257 : 258 : SEGDES IGEOM,MELEMS 259 : 260 : C***************************************************************************** 261 : C Lecture du coefficient 262 : c write(6,*)'Lecture des COEF ' 263 : 263 : CALL ACME(MTABX,'IARG',IARG) 264 : c write(6,*)' YFRT IARG=',IARG,' KKALE=',KKALE 265 : 266 : C--------- CAS FROT Formulation non conservative ---------------------------- 267 : IF(NOMPER.EQ.'FROT')THEN 268 : 269 : C Coef K={Kx Ky Kz} 270 : NUCOEF=1 271 : c write(6,*)' Coef 1 ',NUCOEF,' MTABZ=',MTABZ 271 : CALL LEKMOF(MTABZ,NUCOEF,MTABX,KINC,1,MCHEL1,KPOIND,0,MCHELG) 272 : IF (IERR.NE.0) RETURN 273 : 274 : C Coef Beta={Betax Betay Betaz} 275 : NUCOEF=NUCOEF+1 276 : c write(6,*)' Coef 2 ',NUCOEF,' MTABZ=',MTABZ 276 : CALL LEKMOF(MTABZ,NUCOEF,MTABX,KINC,1,MCHEL2,KPOIND,0,MCHELG) 277 : IF (IERR.NE.0) RETURN 278 : 279 : C Un 280 : NUCOEF=NUCOEF+1 281 : c write(6,*)' Coef 3 ',NUCOEF 281 : CALL LEKMOF(MTABZ,NUCOEF,MTABX,KINC,1,MCHEL3,KPOIND,0,MCHELG) 282 : IF (IERR.NE.0) RETURN 283 : 284 : 285 : c write(6,*)'FIN CAS FROT Formulation non conservative' 286 : C----- FIN CAS FROT Formulation non conservative ---------------------------- 287 : 288 : 289 : ENDIF 290 : c write(6,*)' Fin lect COEF ' 291 : C************************************************************************* 292 : C******* CALCUL ********************************************************** 293 : C************************************************************************* 294 : 295 : SEGACT MTETA1 296 : 296 : CALL KRIPAD(MELEMS,MLENTI) 297 : SEGDES MELEMS 298 : 299 : SEGACT MELEME,MELEM2 300 : NBSOUS=LISOUS(/1) 301 : IF(NBSOUS.EQ.0)NBSOUS=1 302 : 303 : C MATRIK . MATRIK . MATRIK . MATRIK . MATRIK . MATRIK MATRIK . MATRIK . 304 : IF(KKMACO.NE.2)THEN 305 : C ... Création / Ajustement des matrices élémentaires 306 : IF(TMATRI)THEN 307 : C Cas RIGIDITE 308 : IF(XRIG)THEN 309 : 310 : NRIGE=8 311 : NRIGEL=0 312 : SEGINI MRIGID 313 : 314 : NRIGE=8 315 : NRI=IRIGEL(/2) 316 : NRIGEL=NBSOUS+NRI 317 : SEGADJ MRIGID 318 : c write(6,*)' NRIGE,NRIGEL,MRIGID=',NRIGE,NRIGEL,MRIGID 319 : 320 : ELSE 321 : C Cas MATRIK 322 : 323 : NRIGE=7 324 : NKID =9 325 : NKMT =7 326 : NMATRI=0 327 : SEGINI MATRIK 328 : c write(6,*)' YFRT creation MATRIK =',matrik 329 : 330 : NRIGE=7 331 : NKID =9 332 : NKMT =7 333 : NMATR0=JRIGEL(/2) 334 : 335 : NBME=1 336 : IF(IHV.EQ.1)NBME=IDIM 337 : NMATRI=NBME 338 : SEGADJ MATRIK 339 : 340 : DO 41 M=1,NBME 341 : JRIGEL(1,M)=MELEME 342 : JRIGEL(2,M)=MELEM2 343 : JRIGEL(7,M)=0 344 : c? IF(TCONV)JRIGEL(7,M)=2 345 : SEGINI JMATRI 346 : JRIGEL(4,M)=JMATRI 347 : KSPGP=MELEMS 348 : KSPGD=MELEMS 349 : LISPRI(1)=NOMP(M)//' ' 350 : LJSDUA(1)=NOMD(M)//' ' 351 : 41 CONTINUE 352 : ENDIF 353 : ENDIF 354 : C ... Fin Création / Ajustement des matrices élémentaires 355 : 356 : 357 : ELSE 358 : C Activation RIGIDITE 359 : IF(XRIG)THEN 360 : SEGACT MRIGID 361 : ELSE 362 : C Activation MATRIK 363 : SEGACT MATRIK 364 : ENDIF 365 : 366 : C ... Activation des matrices élémentaires 367 : IF(TMATRI)THEN 368 : C Cas RIGIDITE 369 : IF(XRIG)THEN 370 : NRIGE=8 371 : NRI=IRIGEL(/2) 372 : NRIGEL=NBSOUS+NRI 373 : SEGAct MRIGID 374 : c write(6,*)' NRIGE,NRIGEL,MRIGID=',NRIGE,NRIGEL,MRIGID 375 : 376 : ELSE 377 : C Cas MATRIK Pleine 378 : NMATRI=JRIGEL(/2) 379 : NBME=NMATRI 380 : 381 : DO 42 M=1,NBME 382 : JMATRI=JRIGEL(4,M) 383 : SEGACT JMATRI*MOD 384 : LISPRI(1)=NOMP(M)//' ' 385 : LJSDUA(1)=NOMD(M)//' ' 386 : 42 CONTINUE 387 : 388 : ENDIF 389 : ENDIF 390 : C ... Fin Activation des matrices élémentaires 391 : 392 : ENDIF 393 : C FIN MATRIK . FIN MATRIK . FIN MATRIK . FIN MATRIK . FIN MATRIK . FIN M 394 : C ...................................................................... 395 : 396 : SEGACT MCHEL1 397 : SEGACT MCHEL2 398 : SEGACT MCHEL3 399 : 400 : IF(MAX(1,MELEM2.LISOUS(/1)).NE.MAX(1,LISOUS(/1)))THEN 401 : WRITE(6,*)' Geometries incompatibles dans ',nomper 402 : C% Données incompatibles 402 : CALL ERREUR(22) 403 : RETURN 404 : ENDIF 405 : 406 : NKD=0 407 : DO 101 L=1,MAX(1,LISOUS(/1)) 408 : IPT1=MELEME 409 : IPT2=MELEM2 410 : IF(LISOUS(/1).NE.0)IPT1=LISOUS(L) 411 : IF(MELEM2.LISOUS(/1).NE.0)THEN 412 : IPT2=MELEM2.LISOUS(L) 413 : NKD=0 414 : ENDIF 415 : SEGACT IPT1,IPT2 416 : 417 : NOM0 = NOMS(IPT1.ITYPEL)//' ' 417 : CALL KALPBG(NOM0,'FONFORM ',IZFFM) 418 : IF(IZFFM.EQ.0)RETURN 419 : SEGACT IZFFM*MOD 420 : IZHR=KZHR(1) 421 : SEGACT IZHR*MOD 422 : IZF1 = KTP(1) 423 : IZH2 = KZHR(2) 424 : IZW = IZFFM 425 : IZWH = IZHR 426 : 427 : NES=GR(/1) 428 : NPG=GR(/3) 429 : 430 : NP = IPT1.NUM(/1) 431 : MP = IPT2.NUM(/1) 432 : C? MP = IZW.FN(/1) ceci doit etre identique 433 : 434 : NBEL=IPT1.NUM(/2) 435 : 436 : SEGINI SAJT 437 : 438 : c....................................................................... 439 : C MATRIK . MATRIK . MATRIK . MATRIK . MATRIK . MATRIK MATRIK . MATRIK . 440 : c write(6,*)' KKMACO=',KKMACO,'TMATRI=',TMATRI 441 : IF(KKMACO.NE.2)THEN 442 : IF(TMATRI)THEN 443 : C Cas RIGIDITE 444 : IF(XRIG)THEN 445 : IRIGEL(1,NRI+L)=MELEME 446 : COERIG(L)=1.D0 447 : 448 : IRIGEL(7,NRI+L)=0 449 : c? IF(TCONV)IRIGEL(7,NRI+L)=2 450 : 451 : NBME=1 452 : IF(IHV.EQ.1)NBME=IDIM 453 : NLIGRP=NP 454 : NLIGRD=MP 455 : SEGINI DESCR 456 : IRIGEL(3,NRI+L)=DESCR 457 : IF(NBME.EQ.1)THEN 458 : DO 102 I=1,NLIGRP 459 : LISINC(I)=NOMP(1)//' ' 460 : NOELEP(I)=I 461 : 102 CONTINUE 462 : DO 103 I=1,NLIGRD 463 : LISDUA(I)=NOMD(1)//' ' 464 : NOELED(I)=I 465 : 103 CONTINUE 466 : ELSE 467 : ENDIF 468 : SEGDES DESCR 469 : 470 : NELRIG=NBEL 471 : SEGINI xMATRI 472 : IRIGEL(4,NRI+L)=xMATRI 473 : c write(6,*)'NELRIG,IMATRI=',NELRIG,IMATRI 474 : 475 : * DO 104 K=1,NELRIG 476 : * SEGINI XMATRI 477 : c write(6,*)'NLIGRD,NLIGRP,XMATRI=',NLIGRD,NLIGRP,XMATRI 478 : * IMATTT(K)=XMATRI 479 : * 104 CONTINUE 480 : 481 : ELSE 482 : C Cas MATRIK 483 : SEGINI IPM1 484 : c mtail=(IPM1.AM(/1))*(ipm1.am(/2))*(ipm1.am(/3)) 485 : JMATR1=JRIGEL(4,NMATR0+1) 486 : JMATR1.LIZAFM(L,1)=IPM1 487 : IPM2=IPM1 488 : IPM3=IPM1 489 : IF(NBME.GE.2)THEN 490 : JMATR2=JRIGEL(4,NMATR0+2) 491 : SEGINI IPM2 492 : c mtail=(IPM2.AM(/1))*(ipm2.am(/2))*(ipm2.am(/3)) 493 : JMATR2.LIZAFM(L,1)=IPM2 494 : ICAL=.TRUE. 495 : ENDIF 496 : IF(NBME.GE.3)THEN 497 : JMATR3=JRIGEL(4,NMATR0+3) 498 : SEGINI IPM3 499 : JMATR3.LIZAFM(L,1)=IPM3 500 : ICAL=.TRUE. 501 : ENDIF 502 : ENDIF 503 : ENDIF 504 : c write(6,*)' YFRT appelé par ',NOMPER,' Taille matrice=',mtail 505 : 506 : ELSE 507 : 508 : C Cas RIGIDITE 509 : IF(XRIG)THEN 510 : 511 : ELSE 512 : C Cas MATRIK 513 : JMATR1=JRIGEL(4,1) 514 : SEGACT JMATR1 515 : IPM1=JMATR1.LIZAFM(L,1) 516 : IPM2=IPM1 517 : IPM3=IPM1 518 : SEGACT IPM1 519 : IF(NBME.GE.2)THEN 520 : JMATR2=JRIGEL(4,2) 521 : SEGACT JMATR2 522 : IPM2=JMATR2.LIZAFM(L,1) 523 : SEGACT IPM2 524 : ENDIF 525 : IF(NBME.GE.3)THEN 526 : JMATR3=JRIGEL(4,3) 527 : IPM3=JMATR3.LIZAFM(L,1) 528 : SEGACT IPM3 529 : ENDIF 530 : ENDIF 531 : ENDIF 532 : C FIN MATRIK . FIN MATRIK . FIN MATRIK . FIN MATRIK . FIN MATRIK . FIN M 533 : C ...................................................................... 534 : 535 : 536 : C----K 537 : IK1=1 538 : MCHAM1=MCHEL1.ICHAML(L) 539 : SEGACT MCHAM1 540 : 541 : MELVA1=MCHAM1.IELVAL(1) 542 : SEGACT MELVA1 543 : N1PTEL=MELVA1.VELCHE(/1) 544 : N1EL=MELVA1.VELCHE(/2) 545 : IF(N1EL.EQ.1)THEN 546 : IK1=1 547 : ELSEIF(N1EL.EQ.NBEL)THEN 548 : IK1=0 549 : ENDIF 550 : 551 : C----Beta 552 : IK2=1 553 : MCHAM2=MCHEL2.ICHAML(L) 554 : SEGACT MCHAM2 555 : 556 : MELVA2=MCHAM2.IELVAL(1) 557 : SEGACT MELVA2 558 : N1PTEL=MELVA2.VELCHE(/1) 559 : N1EL=MELVA2.VELCHE(/2) 560 : IF(N1EL.EQ.1)THEN 561 : IK2=1 562 : ELSEIF(N1EL.EQ.NBEL)THEN 563 : IK2=0 564 : ENDIF 565 : 566 : C----U 567 : IK3=1 568 : MCHAM3=MCHEL3.ICHAML(L) 569 : SEGACT MCHAM3 570 : 571 : MELVA3=MCHAM3.IELVAL(1) 572 : SEGACT MELVA3 573 : N1PTEL=MELVA3.VELCHE(/1) 574 : N1EL=MELVA3.VELCHE(/2) 575 : IF(N1EL.EQ.1)THEN 576 : IK3=1 577 : ELSEIF(N1EL.EQ.NBEL)THEN 578 : IK3=0 579 : ENDIF 580 : 581 : 582 : c write(6,*)' AVANT 108 NC=',NC,' NBEL=',NBEL,MP,NP,NC 583 : C=============================================== 584 : AI1=AIMPL-1.D0 585 : IF(AIMPL.EQ.0.D0)THEN 586 : AI2=-1.D0 587 : ELSE 588 : AI2=AI1/AIMPL 589 : ENDIF 590 : 591 : DO 108 KE=1,NBEL 592 : 593 : NK1=KE + IK1*(1 - KE) 594 : NK2=KE + IK2*(1 - KE) 595 : NK3=KE + IK3*(1 - KE) 596 : 597 : DO 109 I=1,NP 598 : J=IPT1.NUM(I,KE) 599 : DO 109 N=1,IDIM 600 : XYZ(N,I)=XCOOR((J-1)*(IDIM+1)+N) 601 : 109 CONTINUE 602 : 603 : DO 107 I=1,NP 604 : I1=MLENT1.LECT(IPT1.NUM(I,KE)) 605 : DO 107 N=1,NC 606 : TN1(I,N)=MTETA1.VPOCHA(I1,N) 607 : 107 CONTINUE 608 : 608 : CALL CALJBR(FN,GR,PG,XYZ,HR,PGSQ,RPG,NES, 609 : * IDIM,NP,NPG,IAXI,AIRE,AJ,SGN) 610 : 611 : IF (IDCEN.EQ.0.OR.IDCEN.EQ.1) THEN 611 : CALL RSETD(WT,FN,(NP*NPG)) 611 : CALL RSETD(WS,FN,(NP*NPG)) 612 : ELSE 612 : CALL CALDEC(WT,WS,XYZ,GR,HR,FN,NES,IDIM,NP,NPG,AJT, 613 : & IDCEN,CMD,MELVA1.VELCHE(1,NK1),MELVA2.VELCHE(1,NK2), 614 : & MELVA3.VELCHE(1,NK3),TN1,NC,IKOMP,XREF,Aire,KE) 615 : ENDIF 616 : 616 : CALL INITD(SM1,(MP*IDIM),0.D0) 617 : 618 : C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: 619 : IF(KKMACO.EQ.2)THEN 620 : 621 : C ...... Chargement Rigidite ou Matrik 622 : C Cas RIGIDITE 623 : IF(XRIG)THEN 624 : * XMATRI=IMATTT(KE) 625 : DO 911 I=1,MP 626 : DO 911 J=1,NP 627 : RE(I,J,ke)=RF1(J,I,1) 628 : 911 CONTINUE 629 : * SEGDES XMATRI 630 : ELSE 631 : C Cas MATRIK 632 : DO 912 N=1,NBME 633 : JMATR1=JRIGEL(4,N) 634 : IPM4=JMATR1.LIZAFM(L,1) 635 : DO 912 I=1,MP 636 : DO 912 J=1,NP 637 : RF1(J,I,N)=IPM4.AM(KE,J,I) 638 : 912 CONTINUE 639 : ENDIF 640 : C...... Multiplication Matrice Vecteur 641 : 642 : DO 710 I=1,MP 643 : DO 711 J=1,NP 644 : 645 : DO 716 N=1,NC 646 : SM1(I,N) = SM1(I,N) + AI2*RF1(J,I,N)*TN1(J,N) 647 : 716 CONTINUE 648 : 649 : 711 CONTINUE 650 : 710 CONTINUE 651 : 652 : ELSEIF(KKMACO.NE.2)THEN 653 : C.............................................................. 654 : C...... Frot .... Formulation non conservative 655 : C préparation 656 : 657 : DO 416 N=1,NC 658 : DO 410 I=1,MP 659 : c DO 411 J=1,NP 660 : U3=0.D0 661 : DO 412 LG=1,NPG 662 : PDR=PGSQ(LG)*DEUPI*RPG(LG) 663 : C1=MELVA1.VELCHE((N-1)*NPG+LG,NK1)*PDR 664 : C2=MELVA2.VELCHE((N-1)*NPG+LG,NK2) 665 : C3=ABS ( MELVA3.VELCHE((N-1)*NPG+LG,NK3) ) + 1.e-30 666 : CK=C1*(C3**(C2-1.D0)) 667 : U3=U3+CK * FN(I,LG) 668 : 412 CONTINUE 669 : 670 : SM1(I,N) = SM1(I,N) + AI1*U3*TN1(I,N) 671 : RF1(I,I,N) = AIEX*U3 672 : c? SM1(I,N) = SM1(I,N) + AI1*U3*TN1(J,N) 673 : c? RF1(J,I,N) = AIEX*U3 674 : 675 : 411 CONTINUE 676 : 410 CONTINUE 677 : 416 CONTINUE 678 : 679 : 680 : C.............................................................. 681 : 682 : 683 : ENDIF 684 : 685 : C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: 686 : C======================================================================= 687 : 688 : C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 689 : IF(TMATRI.AND.KKMACO.NE.2)THEN 690 : C ...... Chargement Rigidite ou Matrik 691 : C Cas RIGIDITE 692 : IF(XRIG)THEN 693 : * XMATRI=IMATTT(KE) 694 : DO 901 I=1,NP 695 : DO 901 J=1,NP 696 : RE(I,J,ke)=RF1(J,I,1) 697 : 901 CONTINUE 698 : * SEGDES XMATRI 699 : ELSE 700 : C Cas MATRIK 701 : NBMM=1 702 : IF(ICAL)NBMM=NBME 703 : DO 902 N=1,NBMM 704 : JMATR1=JRIGEL(4,NMATR0+N) 705 : IPM4=JMATR1.LIZAFM(L,1) 706 : DO 902 I=1,NP 707 : DO 902 J=1,NP 708 : IPM4.AM(KE,J,I)=RF1(J,I,N) 709 : 902 CONTINUE 710 : c write(6,*)' RF1 ' 711 : c write(6,1002)(((RF1(J,I,N),j=1,np),i=1,np),n=1,nbmm) 712 : ENDIF 713 : ENDIF 714 : 715 : C ...... Chargement Second membre 716 : DO 910 I=1,NP 717 : I1=LECT(IPT1.NUM(I,KE)) 718 : DO 910 N=1,NC 719 : MPOVA1.VPOCHA(I1,N)=MPOVA1.VPOCHA(I1,N)+SM1(I,N) 720 : 910 CONTINUE 721 : C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 722 : 723 : 108 CONTINUE 724 : 725 : SEGSUP MCHAM1,MELVA1 726 : SEGSUP MCHAM2,MELVA2 727 : SEGSUP MCHAM3,MELVA3 728 : 729 : 730 : 731 : SEGDES IPT1,IPT2 732 : 733 : IF(TMATRI)THEN 734 : C Cas RIGIDITE 735 : IF(XRIG)THEN 736 : SEGDES xMATRI 737 : ELSE 738 : C Cas MATRIK 739 : SEGDES IPM1 740 : IF(NBME.GE.2)SEGDES IPM2 741 : IF(NBME.GE.3)SEGDES IPM3 742 : ENDIF 743 : ENDIF 744 : 745 : SEGSUP IZFFM,IZHR,IZF1,IZH2 746 : SEGSUP SAJT 747 : 748 : 101 CONTINUE 749 : 750 : IF(TMATRI)THEN 751 : IF(.NOT.XRIG)THEN 752 : NBMM=JRIGEL(/2) 753 : IF(NBMM.NE.0)THEN 754 : DO 141 M=1,NBMM 755 : JMATRI=JRIGEL(4,M) 756 : SEGDES JMATRI 757 : 141 CONTINUE 758 : ENDIF 759 : ENDIF 760 : ENDIF 761 : 762 : C""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" 763 : C""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" 764 : C""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" 765 : 766 : 767 : SEGSUP MCHEL1 768 : SEGSUP MCHEL2 769 : SEGSUP MCHEL3 770 : 771 : SEGDES MCHPO1,MPOVA1 772 : SEGDES MELEME,MELEM2 773 : 774 : SEGSUP MLENTI 775 : c IF(TLINCO)THEN 776 : SEGSUP MLENT1 777 : SEGDES MTETA1 778 : c ENDIF 779 : 780 : IF(TMATRI)THEN 781 : C Cas RIGIDITE 782 : IF(XRIG)THEN 783 : SEGDES MRIGID 784 : ELSE 785 : C Cas MATRIK 786 : SEGDES MATRIK 787 : ENDIF 788 : ENDIF 789 : 790 : 791 : c segact MCHPO1 792 : c MSOUP1=MCHPO1.IPCHP(1) 793 : c segact MSOUP1 794 : c write(6,*)MCHPO1.IPCHP(/1),MSOUP1.NOCOMP(/2) 795 : c write(6,*)(MSOUP1.NOCOMP(i),i=1,MSOUP1.NOCOMP(/2)) 795 : c if(ihv.eq.1)call ecmo(mtabz,'TOTO','CHPOINT',MCHPO1) 796 : 797 : C************************************************************************* 798 : 799 : IF(KKMACO.EQ.1)THEN 800 : TYPE='MATRIK' 800 : CALL ECMO(MTAB1,NMACO,TYPE,MATRIK) 801 : ENDIF 802 : 803 : IF(AIMPL.EQ.0.D0)THEN 804 : NRIGE=7 805 : NKID =9 806 : NKMT =7 807 : NMATRI=0 808 : SEGINI MATRIK 809 : ENDIF 810 : 811 : c write(6,*)' FIN YFRT' 812 : RETURN 813 : 1001 FORMAT(20(1X,I5)) 814 : 1002 FORMAT(10(1X,1PE11.4)) 815 : END 816 : 817 : 818 : 819 : 820 : 821 : 822 : 823 : 824 : 825 : 826 : 827 : 828 : 829 : 830 : 831 : 832 : 833 : 834 : 835 : 836 : 837 : 838 : 839 :
© Cast3M 2003 - All rights reserved.
Disclaimer