1 : C YDFDT SOURCE MAGN 10/05/18 21:16:53 6675 2 : SUBROUTINE YDFDT 3 : C--------------------------------------------------------------------- 4 : C Cet opérateur discrétise le terme de dérivée temporelle. 5 : C--------------------------------------------------------------------- 6 : C Les formulations numériques sont EF Galerkin, EF Petrov-Galerkin et 7 : C VF, la matrice masse pouvant etre condensée en EF. 8 : C Les schema en temps disponibles CRANK NICOLSON (DT**2) CNG (DT**4) 9 : C Les schema en temps disponibles BDF2 (DT**2) BDF4 (DT**4) 10 : C--------------------------------------------------------------------- 11 : C 12 : C------------------------- 13 : C Phrase d'appel GIBIANE : 14 : C------------------------- 15 : C 16 : C 'DFDT' TAB1 ; 17 : C 18 : C TAB1 : Table de sous-type 'KIZX' contenant les données pour DFDT 19 : C 20 : C-------------------------------- 21 : C Construction de TAB1 via EQEX : 22 : C-------------------------------- 23 : C 24 : C 'ZONE' TAB2 'OPER' 'DFDT' CHP1 CHP2 FLOT1 (CHP3) (CHP4) 'INCO' MOT1 25 : C 26 : C TAB2 : TABLE DOMAINE associée à la zone d'action de l'opérateur 27 : C CHP1 : Coefficient multiplicateur de la dérivée temporelle 28 : C (densité dans le cas de la qdm ou rho*cp en thermique) 29 : C (FLOTTANT ou MOT ou CHPO) - (spg : CENTRE ou SOMMET) 30 : C CHP2 : Inconnue au pas de temps précédant 31 : C (MOT ou CHPO SCAL ou CHPO VECT) - (spg : CENTRE ou SOMMET) 32 : C FLOT1 : Valeur du pas de temps (MOT ou REEL ou ENTIER) 33 : C CHP3 : Champ de vitesse - Optionnel (utilisé en Petrov-Galerkin) 34 : C (POINT ou MOT ou CHPO) - (spg : CENTRE ou SOMMET) 35 : C CHP4 : Diffusion - Optionnel (utilisé en Petrov-Galerkin) 36 : C (FLOTTANT ou MOT ou CHPO) - (spg : CENTRE ou SOMMET) 37 : C MOT1 : Nom de l'inconnue sur laquelle agit l'opérateur (MOT) 38 : C 39 : C------------ 40 : C Résultats : 41 : C------------ 42 : C 43 : C Les matrices élémentaires sont stockés dans un objet MATRIK qui 44 : C est rangé à l'indice de type MOT MATELM de la table KIZX. 45 : C Le second membre est stocké dans un CHPO et assemblé dans la table 46 : C EQEX à l'indice de type MOT SMBR. 47 : C--------------------------------------------------------------------- 48 : C 49 : IMPLICIT INTEGER(I-N) 50 : IMPLICIT REAL*8 (A-H,O-Z) 51 : C 52 : -INC CCOPTIO 53 : -INC CCREEL 54 : *- 55 : -INC CCGEOME 56 : -INC SIZFFB 57 : POINTEUR IZF1.IZFFM 58 : 59 : PARAMETER (LRV=64) 60 : SEGMENT PETROV 61 : REAL*8 WT(LRV,NP,NPG,KDIM),WS(LRV,NP,NPG,KDIM),HK(LRV,IDIM,NP,NPG) 62 : REAL*8 UIL(LRV,IDIM,NPG),DUIL(LRV,IDIM,NPG) 63 : REAL*8 PGSK(LRV,NPG),RPGK(LRV,NPG),AIRE(LRV),ANUK(LRV),COEFK(LRV) 64 : REAL*8 AJK(LRV,IDIM,IDIM,NPG) 65 : ENDSEGMENT 66 : 67 : -INC SMCOORD 68 : -INC SMLENTI 69 : POINTEUR IPADU.MLENTI,IPADS.MLENTI,IPADI.MLENTI 70 : -INC SMELEME 71 : POINTEUR MELEM1.MELEME,MELEMS.MELEME,MELEMZ.MELEME 72 : POINTEUR MELEMO.MELEME 73 : -INC SMCHAML 74 : -INC SMCHPOI 75 : POINTEUR IZRO.MPOVAL,MZDT.MPOVAL,IZTGG2.MPOVAL,IZTU1.MPOVAL 76 : POINTEUR IZTGN2.MPOVAL 77 : POINTEUR IZTCO.MPOVAL,MZUN.MPOVAL,MZMU.MPOVAL 78 : -INC SMMATRIK 79 : -INC SMLMOTS 80 : C 81 : CHARACTER*8 NOMI,TYPE,TYPC,NOM,NOM0,NOMDT,MTYP,CHAI 82 : PARAMETER (NTB=1) 83 : CHARACTER*8 LTAB(NTB) 84 : DIMENSION KTAB(NTB),IXV(4) 85 : DATA LTAB/'KIZX '/ 86 : C 87 : C DFDT 88 : C write(6,*)' Opérateur DFDT ' 89 : 90 : C Nouvelle directive EQUA de EQEX 91 : MTYP=' ' 91 : CALL QUETYP(MTYP,0,IRET) 92 : IF(IRET.EQ.0)THEN 93 : C% On attend un des objets : %m1:8 %m9:16 %m17:24 %m25:32 %m33:40 94 : MOTERR( 1: 8) = 'CHAI ' 95 : MOTERR( 9:16) = 'MMODEL ' 96 : MOTERR(17:24) = 'TABLE ' 96 : CALL ERREUR(472) 97 : RETURN 98 : ENDIF 99 : 100 : IF(MTYP.EQ.'MMODEL')THEN 100 : CALL YTCLSF('T ','DFDT ') 101 : RETURN 102 : ELSEIF(MTYP.EQ.'MOT ')THEN 102 : CALL LIRCHA(CHAI,1,IRET) 102 : CALL YTCLSF(CHAI,'DFDT ') 103 : RETURN 104 : ENDIF 105 : C Fin Nouvelle directive EQUA de EQEX 106 : 107 : DTM1 = 0.D0 108 : DTX = 0.D0 109 : MLENTI = 0 110 : IAXI = 0 111 : NPTU = 0 112 : IF (IFOMOD.EQ.0) IAXI=2 113 : DEUPI=1.D0 114 : IF(IAXI.NE.0)DEUPI=2.D0*XPI 115 : C 116 : C- Lecture de la table KIZX (pointeur MTABX) associée à DFDT 117 : C 117 : CALL LITABS(LTAB,KTAB,NTB,1,IRET) 118 : IF (IERR.NE.0) RETURN 119 : MTABX = KTAB(1) 120 : C 121 : C- Récupération de la table EQEX (pointeur MTAB1) 122 : C 122 : CALL LEKTAB(MTABX,'EQEX',MTAB1) 123 : IF (MTAB1.EQ.0)THEN 124 : C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24 125 : MOTERR( 1: 8) = ' EQEX ' 126 : MOTERR( 9:16) = ' EQEX ' 127 : MOTERR(17:24) = ' KIZX ' 127 : CALL ERREUR(786) 128 : RETURN 129 : ENDIF 129 : CALL ACME(MTAB1,'NAVISTOK',NASTOK) 130 : IF(NASTOK.EQ.0)THEN 130 : CALL ZDFDT(MTABX,MTAB1) 131 : RETURN 132 : ENDIF 133 : C 134 : C- Récupération de la table DOMAINE associée au domaine local 135 : C 136 : TYPE=' ' 136 : CALL ACMO(MTABX,'DOMZ',TYPE,MMDZ) 137 : IF(TYPE.NE.'MMODEL')THEN 138 : C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24 139 : MOTERR( 1: 8) = ' DOMZ ' 140 : MOTERR( 9:16) = ' DOMZ ' 141 : MOTERR(17:24) = ' KIZX ' 141 : CALL ERREUR(786) 142 : RETURN 143 : ENDIF 144 : 144 : CALL LEKMOD(MMDZ,MTABZ,INEFMD) 145 : C 146 : C- Récupération de la table INCO (pointeur KINC) 147 : C 147 : CALL LEKTAB(MTAB1,'INCO',KINC) 148 : IF (KINC.EQ.0) THEN 149 : C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24 150 : MOTERR( 1: 8) = ' INCO ' 151 : MOTERR( 9:16) = ' INCO ' 152 : MOTERR(17:24) = ' EQEX ' 152 : CALL ERREUR(786) 153 : RETURN 154 : ENDIF 155 : C 156 : C- Récupération de la table PASDETPS (pointeur KINT) 157 : C 157 : CALL LEKTAB(MTAB1,'PASDETPS',KINT) 158 : C?2 IF(KINT.EQ.0)KINT=KINC 159 : C 160 : C- Récupération de la table des options KOPT (pointeur KOPTI) 161 : C 161 : CALL LEKTAB(MTABX,'KOPT',KOPTI) 162 : IF (KOPTI.EQ.0) THEN 163 : C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24 164 : MOTERR( 1: 8) = ' KOPT ' 165 : MOTERR( 9:16) = ' KOPT ' 166 : MOTERR(17:24) = ' KIZX ' 166 : CALL ERREUR(786) 167 : RETURN 168 : ELSE 169 : IKOMP=0 169 : CALL ACME(KOPTI,'ISCHT',ISCHT) 170 : C- Initialisation des constantes en fonction du schema 171 : CDT=1.D0 172 : CN=1.D0 173 : CNM=1.D0 174 : CNMM=0.D0 175 : IF(ISCHT.EQ.1)THEN 176 : CDT=2.D0 177 : CN=3.D0 178 : CNM=4.D0 179 : CNMM=-1.D0 180 : ELSEIF(ISCHT.EQ.2)THEN 181 : CDT=2.D0 182 : CN=6.D0 183 : CNM=0.D0 184 : CNMM=6.D0 185 : ENDIF 186 : 187 : KDIM=1 187 : CALL ACME(KOPTI,'IDCEN',IDCEN) 188 : IF(IDCEN.EQ.2)KDIM=IDIM 188 : CALL ACME(KOPTI,'KFORM',KFORM) 188 : CALL ACME(KOPTI,'KPOIND',KPOIND) 188 : CALL ACME(KOPTI,'KMACO',KMACO) 188 : CALL ACMF(KOPTI,'CMD',CMD) 189 : ENDIF 190 : 191 : IDCENV = IDCEN 192 : IF (IDCEN.EQ.5.OR.IDCEN.EQ.4) IDCENV=1 193 : C 194 : C- Récupération des indices CENTRE, FACE, SOMMET et MAILLAGE 195 : C- de la table DOMAINE 196 : C 196 : CALL LEKTAB(MTABZ,'SOMMET',MELEMS) 196 : CALL LEKTAB(MTABZ,'FACE',MELEMF) 196 : CALL LEKTAB(MTABZ,'CENTRE',MELEMC) 196 : CALL LEKTAB(MTABZ,'MAILLAGE',MELEMZ) 196 : CALL LEKTAB(MTABZ,'MACRO1',MELEMI) 196 : CALL LEKTAB(MTABZ,'MACRO',MACRO) 196 : CALL LEKTAB(MTABZ,'QUADRATI',MQUAD) 197 : IF (IERR.NE.0) RETURN 198 : 199 : C 200 : C- Vérification des compatiblités Formulation/SPG 201 : C- Identification du spg de l'inconnue 202 : C- MELEM1=spg ; MELEME=connectivité associée ; 203 : C 204 : 205 : IF(KFORM.EQ.0)THEN 206 : IF(KPOIND.EQ.99)THEN 207 : KPOIND=0 208 : MELEM1=MELEMS 209 : MELEME=MELEMS 210 : ELSEIF(KPOIND.EQ.2)THEN 211 : MELEM1=MELEMC 212 : MELEME=MELEMC 213 : ELSEIF(KPOIND.EQ.0)THEN 214 : MELEM1=MELEMS 215 : MELEME=MELEMS 216 : ELSEIF(KPOIND.NE.2.AND.KPOIND.NE.0)THEN 217 : C Option %m1:8 incompatible avec les données 218 : MOTERR( 1: 8) = ' EF ' 218 : CALL ERREUR(803) 219 : RETURN 220 : ENDIF 221 : 222 : ELSEIF(KFORM.EQ.1)THEN 223 : IF(KPOIND.EQ.99)THEN 224 : KPOIND=0 225 : MELEM1=MELEMS 226 : MELEME=MELEMZ 227 : ELSEIF(KPOIND.EQ.0)THEN 228 : MELEM1=MELEMS 229 : MELEME=MELEMZ 230 : ELSEIF(KPOIND.EQ.2)THEN 231 : MELEM1=MELEMC 232 : MELEME=MELEMC 233 : ELSEIF(KPOIND.NE.2.AND.KPOIND.NE.0)THEN 234 : C Option %m1:8 incompatible avec les données 235 : MOTERR( 1: 8) = ' EF ' 235 : CALL ERREUR(803) 236 : RETURN 237 : ENDIF 238 : 239 : ELSEIF(KFORM.EQ.2)THEN 240 : IF(KPOIND.EQ.99)THEN 241 : KPOIND=2 242 : MELEM1=MELEMC 243 : MELEME=MELEMC 244 : ELSEIF(KPOIND.EQ.2)THEN 245 : MELEM1=MELEMC 246 : MELEME=MELEMC 247 : ELSEIF(KPOIND.NE.2)THEN 248 : C Option %m1:8 incompatible avec les données 249 : MOTERR( 1: 8) = ' VF ' 249 : CALL ERREUR(803) 250 : RETURN 251 : ENDIF 252 : 253 : ELSEIF(KFORM.EQ.3)THEN 254 : IF(KPOIND.EQ.99.OR.KPOIND.EQ.2)THEN 255 : KPOIND=2 256 : MELEM1=MELEMC 257 : MELEME=MELEMC 258 : ELSEIF(KPOIND.EQ.3.AND.(MACRO.NE.0.OR.MQUAD.NE.0))THEN 258 : CALL LEKTAB(MTABZ,'CENTREP0',MELEMC) 259 : MELEM1=MELEMC 260 : MELEME=MELEMC 261 : ELSEIF(KPOIND.EQ.4.AND.(MACRO.NE.0.OR.MQUAD.NE.0))THEN 261 : CALL LEKTAB(MTABZ,'CENTREP1',MELEMC) 261 : CALL LEKTAB(MTABZ,'ELTP1NC',MELEMQ) 262 : MELEM1=MELEMC 263 : MELEME=MELEMQ 264 : IF(MACRO.NE.0)MELEMO=MELEMI 265 : IF(MQUAD.NE.0)MELEMO=MELEMZ 266 : C ELSEIF(KPOIND.NE.2.AND.KPOIND.NE.3.AND.KPOIND.NE.4)THEN 267 : ELSE 268 : C Option %m1:8 incompatible avec les données 269 : MOTERR( 1: 8) = ' EFMC ' 269 : CALL ERREUR(803) 270 : RETURN 271 : ENDIF 272 : 273 : ENDIF 274 : C 275 : C- Récupération du nombre d'inconnues et du nom de l'inconnue NOMI 276 : C 277 : TYPE = 'LISTMOTS' 277 : CALL ACMO(MTABX,'LISTINCO',TYPE,MLMOTS) 278 : IF (IERR.NE.0) RETURN 279 : SEGACT MLMOTS 280 : NBINC = MOTS(/2) 281 : IF (NBINC.LE.0.OR.NBINC.GE.3) THEN 282 : C Indice %m1:8 : contient plus de %i1 %m9:16 283 : MOTERR( 1:8) = 'LISTINCO' 284 : INTERR(1) = 2 285 : MOTERR(9:16) = ' MOTS ' 285 : CALL ERREUR(799) 286 : RETURN 287 : ENDIF 288 : NOMI = MOTS(1) 289 : IF (NBINC.EQ.2) THEN 290 : IF (MOTS(1).NE.MOTS(2)) THEN 291 : C Indice %m1:8 : contient plus de %i1 %m9:16 292 : MOTERR( 1:8) = 'LISTINCO' 293 : INTERR(1) = 1 294 : MOTERR(9:16) = ' MOT ' 294 : CALL ERREUR(799) 295 : RETURN 296 : ENDIF 297 : ENDIF 298 : SEGDES MLMOTS 299 : C 300 : C- Récupération de l'inconnue 301 : C 302 : TYPE = ' ' 302 : CALL ACMO(KINC,NOMI,TYPE,MCHPOI) 303 : IF (TYPE.NE.'CHPOINT ') THEN 304 : C Indice %m1:8 : ne contient pas un objet de type %m9:16 305 : MOTERR( 1: 8) = 'INC '//NOMI 306 : MOTERR( 9:16) = 'CHPOINT ' 306 : CALL ERREUR(800) 307 : RETURN 308 : ELSE 308 : CALL LICHT(MCHPOI,IZTU1,TYPC,MELEMD) 309 : NINKO = IZTU1.VPOCHA(/2) 310 : IF (NINKO.NE.1.AND.NINKO.NE.IDIM) THEN 311 : C Indice %m1:8 : Le %m9:16 n'a pas le bon nombre de composantes 312 : MOTERR( 1: 8) = 'INC '//NOMI 313 : MOTERR( 9:16) = 'CHPOINT ' 313 : CALL ERREUR(784) 314 : RETURN 315 : ENDIF 316 : ENDIF 317 : C 318 : C************************************************************************* 319 : C Le domaine de definition est donne par le SPG de la premiere inconnue 320 : C Les inconnues suivantes devront posseder ce meme pointeur 321 : C On verifie que les points de la zone sont tous inclus dans ce SPG 322 : 322 : CALL KRIPAD(MELEMD,IPADI) 323 : IPADS=IPADI 324 : IPADU=IPADI 324 : IF(MELEMD.NE.MELEMS)CALL KRIPAD(MELEMS,IPADS) 325 : MLENTI=IPADI 325 : IF(MELEM1.NE.MELEMD)CALL KRIPAD(MELEM1,MLENTI) 326 : 326 : CALL VERPAD(IPADI,MELEM1,IRET) 327 : IF (IRET.NE.0) THEN 328 : C Indice %m1:8 : L'objet %m9:16 n'a pas le bon support géométrique 329 : MOTERR(1: 8) = 'INC '//NOMI 330 : MOTERR(9:16) = 'CHPOINT ' 330 : CALL ERREUR(788) 331 : RETURN 332 : ENDIF 333 : C 334 : C- Test du nombre d'arguments 335 : C 336 : ICD=0 337 : IF(ISCHT.NE.0)ICD=1 337 : CALL ACME(MTABX,'IARG',IARG) 338 : IF (IDCENV.EQ.1) THEN 339 : IF (IARG.NE.(3+ICD).AND.IARG.NE.(5+ICD)) THEN 340 : C Indice %m1:8 : nombre d'arguments incorrect 341 : MOTERR(1:8) = 'IARG ' 341 : CALL ERREUR(804) 342 : RETURN 343 : ENDIF 344 : ELSE 345 : IF (IARG.NE.(5+ICD)) THEN 346 : C Indice %m1:8 : nombre d'arguments incorrect 347 : MOTERR(1:8) = 'IARG ' 347 : CALL ERREUR(804) 348 : RETURN 349 : ENDIF 350 : ENDIF 351 : C 352 : C- Lecture du coefficient multiplicateur SCAL CENTRE ou SCAL SOMMET 353 : C 354 : IXV(1) = MELEMC 355 : IXV(2) = 1 356 : IXV(3) = 0 357 : IXV(4) = MELEMS 358 : IRET = 4 358 : CALL LEKCOF(' Operateur DFDT :', 359 : & MTABX,KINC,1,IXV,MRO,IZRO,NPT1,NC1,IK1,IRET) 360 : IF (IRET.EQ.0) RETURN 361 : C 362 : C- Lecture de l'inconnue au pas de temps précédant 363 : C 364 : IF (NINKO.EQ.1) THEN 365 : IXV(1) = MELEMD 366 : ELSE 367 : IXV(1) =-MELEMD 368 : ENDIF 369 : IXV(2) = 0 370 : IXV(3) = 0 370 : CALL LEKCOF(' Operateur DFDT :', 371 : & MTABX,KINC,2,IXV,IZTG2,IZTGG2,NPT2,NC2,IK2,IRET) 372 : IF (IRET.EQ.0) RETURN 373 : C 374 : C- Lecture de l'inconnue au pas de temps n-2 (ISCHT NE 0) 375 : C 376 : IF (NINKO.EQ.1) THEN 377 : IXV(1) = MELEMD 378 : ELSE 379 : IXV(1) =-MELEMD 380 : ENDIF 381 : IXV(2) = 0 382 : IXV(3) = 0 382 : CALL LEKCOF(' Operateur DFDT :', 383 : & MTABX,KINC,(2+ICD),IXV,IZTG22,IZTGN2,NPT21,NC21, 384 : & IK21,IRET) 385 : IF (IRET.EQ.0) RETURN 386 : C 387 : C- Lecture du pas de temps 388 : C 389 : 390 : TYPE=' ' 390 : CALL ACMO(MTABX,'ARG3',TYPE,KK) 391 : IF(TYPE.EQ.'FLOTTANT')THEN 391 : CALL ACMF(MTABX,'ARG3',DTR) 391 : IF(KINT.NE.0)CALL ECMF(KINT,'DELTAT',DTR) 392 : ELSEIF(TYPE.EQ.'MOT')THEN 392 : CALL ACMM(MTABX,'ARG3',NOMDT) 393 : IF(NOMDT.NE.'DELTAT')THEN 394 : 395 : IXV(1) = 0 396 : IXV(2) = 1 397 : IXV(3) = 0 397 : CALL LEKCOF(' Operateur DFDT :', 398 : & MTABX,KINC,(3+ICD),IXV,IZTG3,MZDT,NPT3,NC3,IK3,IRET) 399 : IF (IRET.EQ.0) RETURN 400 : DTR = MZDT.VPOCHA(1,1) 400 : IF(KINT.NE.0)CALL ECMF(KINT,'DELTAT',DTR) 401 : 402 : ELSE 403 : IF(KINT.EQ.0)THEN 404 : C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24 405 : MOTERR( 1: 8) = 'PASDETPS' 406 : MOTERR( 9:16) = 'PASDETPS' 407 : MOTERR(17:24) = ' EQEX ' 407 : CALL ERREUR(786) 408 : RETURN 409 : ENDIF 409 : CALL ACMF(KINT,'DELTAT',DTR) 409 : CALL ACMF(MTAB1,'ALFA',ALFA) 410 : C write(6,*)' DFDT nomdt=',nomdt,' DTR=',dtr,alfa,(DTR*ALFA) 411 : DTR=DTR*ALFA 412 : ENDIF 413 : ENDIF 414 : 415 : DTR = DTR * CDT 416 : 417 : C 418 : C- Lecture des données pour le décentrement 419 : C- Champ de vitesse et viscosité 420 : C 421 : IKU=1 422 : MZUN=IZTU1 423 : MZMU=IZRO 424 : IKMU=1 425 : 426 : IF (IARG.EQ.5) THEN 427 : IXV(1) =-MELEMS 428 : IXV(2) = 0 429 : IXV(3) = 1 429 : CALL LEKCOF(' Operateur DFDT :', 430 : & MTABX,KINC,(4+ICD),IXV,MUN,MZUN,NPTU,NCU,IKU,IRET) 431 : IF (IRET.EQ.0) RETURN 432 : IPADU=IPADS 433 : IF (IKU.EQ.2) IKU=1 434 : C 435 : IXV(1) = MELEMC 436 : IXV(2) = 1 437 : IXV(3) = 0 438 : IXV(4) = MELEMS 439 : IRET = 4 439 : CALL LEKCOF(' Operateur DFDT :', 440 : & MTABX,KINC,(5+ICD),IXV,MMU,MZMU,NPTMU,NCMU,IKMU,IRET) 441 : IF (IRET.EQ.0) RETURN 442 : ENDIF 443 : C 444 : C- Transformation éventuelle des arguments du SOMMET vers le CENTRE 445 : C 446 : IF (IK1.EQ.4) THEN 446 : CALL COSOCE(MELEMC,MELEMS,MELEMZ,IZRO,MPOVAL) 447 : IZRO = MPOVAL 448 : IK1 = 0 449 : IK1DETR = 1 450 : ELSE 451 : IK1DETR = 0 452 : ENDIF 453 : IF (IKMU.EQ.4) THEN 453 : CALL COSOCE(MELEMC,MELEMS,MELEMZ,MZMU,MPOVAL) 454 : MZMU = MPOVAL 455 : IKMU = 0 456 : IKMUDETR = 1 457 : ELSE 458 : IKMUDETR = 0 459 : ENDIF 460 : C 461 : C- Récupération de caractéristiques géométriques dans le cas EF 462 : C- ou du volume de chaque élément dans le cas VF 463 : C 464 : 465 : IF (KPOIND.EQ.0) THEN 466 : ELSEIF (KPOIND.EQ.2) THEN 466 : CALL LEKTAB(MTABZ,'XXVOLUM',MCHPOI) 467 : IF (IERR.NE.0) RETURN 467 : CALL LICHT(MCHPOI,MPOVAL,TYPC,IGEOM) 468 : ELSEIF (KPOIND.EQ.3.AND.MACRO.EQ.0) THEN 468 : CALL LEKTAB(MTABZ,'XXVOLUM',MCHPOI) 469 : IF (IERR.NE.0) RETURN 469 : CALL LICHT(MCHPOI,MPOVAL,TYPC,IGEOM) 470 : ELSEIF (KPOIND.EQ.3.AND.MACRO.NE.0) THEN 470 : CALL LEKTAB(MTABZ,'VOLUMAC',MCHPOI) 471 : IF (IERR.NE.0) RETURN 471 : CALL LICHT(MCHPOI,MPOVAL,TYPC,IGEOM) 472 : ENDIF 473 : C 474 : C- Récupération/Construction du chapeau de l'objet MATRIK 475 : C 476 : TYPE = ' ' 476 : CALL ACMO(MTABX,'MATELM',TYPE,MATRIK) 477 : IF (TYPE.EQ.'MATRIK'.AND.KMACO.NE.0) THEN 478 : SEGACT MATRIK 479 : NMATRI = IRIGEL(/2) 480 : MELEME = IRIGEL(1,1) 481 : SEGACT MELEME 482 : IMATRI = IRIGEL(4,1) 483 : SEGACT IMATRI 484 : NBME = LIZAFM(/2) 485 : NINKO = NBME 486 : NBSOUS = LIZAFM(/1) 487 : MELEM1 = KSPGP 488 : ELSE 489 : NRIGE = 7 490 : NKID = 9 491 : NKMT = 7 492 : NMATRI = 1 493 : SEGINI MATRIK 494 : IRIGEL(1,1) = MELEME 495 : IRIGEL(2,1) = MELEME 496 : 497 : IF (KFORM.EQ.0.OR.KFORM.EQ.2) THEN 498 : C EFM1 ou VF -> Diagonale 499 : IRIGEL(7,1) = 5 500 : ELSEIF (KFORM.EQ.1) THEN 501 : C EF matrice pleine sym si centree non sinon 502 : IF (IDCENV.EQ.1) THEN 503 : IRIGEL(7,1) = 0 504 : ELSE 505 : IRIGEL(7,1) = 2 506 : ENDIF 507 : ELSEIF (KFORM.EQ.3) THEN 508 : IF(KPOIND.EQ.2.OR.KPOIND.EQ.3)THEN 509 : IRIGEL(7,1) = 5 510 : ELSEIF(KPOIND.EQ.4)THEN 511 : IRIGEL(7,1) = 0 512 : ENDIF 513 : ENDIF 514 : 515 : NBME = NINKO 516 : SEGACT MELEME 517 : NBSOUS = LISOUS(/1) 518 : IF (NBSOUS.EQ.0) NBSOUS=1 519 : SEGINI IMATRI 520 : IRIGEL(4,1) = IMATRI 521 : KSPGP = MELEM1 522 : KSPGD = MELEM1 523 : IF (NBME.EQ.1) THEN 524 : LISPRI(1) = NOMI(1:4)//' ' 525 : LISDUA(1) = NOMI(1:4)//' ' 526 : ELSE 527 : DO 10 I=1,NBME 528 : WRITE(NOM,FMT='(I1,A7)')I,NOMI(1:7) 529 : LISPRI(I) = NOM(1:4)//' ' 530 : LISDUA(I) = NOM(1:4)//' ' 531 : 10 CONTINUE 532 : ENDIF 533 : ENDIF 534 : C 535 : C- Construction des segments IZAFM contenant les matrices élémentaires 536 : C 537 : IF (KMACO.NE.2) THEN 538 : DO 30 L=1,NBSOUS 539 : IPT1 = MELEME 540 : IF (NBSOUS.NE.1) IPT1=LISOUS(L) 541 : SEGACT IPT1 542 : NP = IPT1.NUM(/1) 543 : MP = NP 544 : NBEL = IPT1.NUM(/2) 545 : SEGDES IPT1 546 : DO 20 I=1,NBME 547 : SEGINI IZAFM 548 : LIZAFM(L,I) = IZAFM 549 : 20 CONTINUE 550 : 30 CONTINUE 551 : ENDIF 552 : C 553 : C--------------------- 554 : C Calcul des matrices 555 : C--------------------- 556 : C 557 : IF (KMACO.NE.2) THEN 558 : C 559 : C-------------------------------------------------- 560 : C-- Cas ou l'inconnue est au SOMMET (EF ou EFM1) -- 561 : C-------------------------------------------------- 562 : C 563 : IF (KPOIND.EQ.0) THEN 564 : K0 = 0 565 : IF (KFORM.EQ.0) THEN 566 : MELEME = MELEMZ 567 : SEGACT MELEME 568 : NBSOUS = LISOUS(/1) 569 : IF (NBSOUS.EQ.0) NBSOUS=1 570 : ENDIF 571 : DO 140 L=1,NBSOUS 572 : IPT1 = MELEME 573 : IF (NBSOUS.NE.1) IPT1=LISOUS(L) 574 : SEGACT IPT1 575 : NP = IPT1.NUM(/1) 576 : NBEL = IPT1.NUM(/2) 577 : NOM0 = NOMS(IPT1.ITYPEL)//' ' 577 : CALL KALPBG(NOM0,'FONFORM ',IZFFM) 578 : IF (IZFFM.EQ.0) THEN 579 : C Echec lors de la lecture des fonctions de forme d'un élément 579 : CALL ERREUR(662) 580 : RETURN 581 : ENDIF 582 : SEGACT IZFFM 583 : IZHR = KZHR(1) 584 : SEGACT IZHR*MOD 585 : NPG = GR(/3) 586 : NES = GR(/1) 587 : C 588 : C- Allocation des tableaux de travail pour PETROV-GALERKIN 589 : C 590 : SEGINI PETROV 591 : C 592 : C- Traitement des éléments par paquets de LRV éléments 593 : C 594 : NNN = MOD(NBEL,LRV) 595 : IF (NNN.EQ.0) THEN 596 : NPACK = NBEL / LRV 597 : ELSE 598 : NPACK = 1 + (NBEL-NNN)/LRV 599 : ENDIF 600 : KPACKD = 1 601 : KPACKF = NPACK 602 : DO 130 KPACK=KPACKD,KPACKF 603 : KDEB = 1 + (KPACK-1)*LRV 604 : KFIN = MIN(NBEL,KDEB+LRV-1) 605 : C 606 : C- Réduction des coefficients sur les éléments du paquet 607 : C 608 : DO 40 K=KDEB,KFIN 609 : KP = K - KDEB + 1 610 : NK = K + K0 611 : NKR = (1-IK1)*(NK-1) + 1 612 : COEFK(KP) = IZRO.VPOCHA(NKR,1) / DTR 613 : 40 CONTINUE 614 : IF(IDCENV.NE.1)THEN 615 : DO 41 K=KDEB,KFIN 616 : KP = K - KDEB + 1 617 : NK = K + K0 618 : NKR = (1-IK1)*(NK-1) + 1 619 : NKM = (1-IKMU)*(NK-1) + 1 620 : ANUK(KP) = MZMU.VPOCHA(NKM,1)/(IZRO.VPOCHA(NKR,1)+XPETIT) 621 : 41 CONTINUE 622 : ENDIF 623 : C 624 : C- Récupération de la fonction de projection suivant IDCEN 625 : C 626 : 627 : IF(IDCENV.EQ.2)THEN 628 : 629 : DO 119 N=1,NBME 629 : CALL SUPGEF(FN,GR,PG,XYZ,HR,PGSQ,RPG,AJ, 630 : & NES,NP,NPG,IAXI,XCOOR, 631 : & IPT1.NUM,KDEB,KFIN,LRV,IDCENV,CMD,IKOMP, 632 : & IZTU1.VPOCHA(1,N),IPADI.LECT,MZUN.VPOCHA,IPADU.LECT,NPTU,ANUK, 633 : & WT(1,1,1,N),WS(1,1,1,N),HK,PGSK,RPGK,AJK,AIRE,UIL,DUIL, 634 : & DTM1,DT,DTT1,DTT2,DIAEL,NUEL) 635 : 119 CONTINUE 636 : 637 : ELSE 638 : 638 : CALL SUPGEF(FN,GR,PG,XYZ,HR,PGSQ,RPG,AJ, 639 : & NES,NP,NPG,IAXI,XCOOR, 640 : & IPT1.NUM,KDEB,KFIN,LRV,IDCENV,CMD,IKOMP, 641 : & IZTU1.VPOCHA,IPADI.LECT,MZUN.VPOCHA,IPADU.LECT,NPTU,ANUK, 642 : & WT(1,1,1,1),WS(1,1,1,1),HK,PGSK,RPGK,AJK,AIRE,UIL,DUIL, 643 : & DTM1,DT,DTT1,DTT2,DIAEL,NUEL) 644 : ENDIF 645 : C 646 : DO 120 N=1,NBME 647 : N1=1 648 : IF(IDCEN.EQ.2)N1=N 649 : C 650 : C- Matrice masse consistante 651 : C 652 : IF (KFORM.EQ.1) THEN 653 : IZAFM = LIZAFM(L,N) 654 : DO 80 K=KDEB,KFIN 655 : KP = K - KDEB + 1 656 : DO 70 I=1,NP 657 : DO 60 J=1,NP 658 : UU = 0.D0 659 : DO 50 LG=1,NPG 660 : UU=UU+FN(I,LG)*WT(KP,J,LG,N1)*PGSK(KP,LG) 661 : & *DEUPI*RPGK(KP,LG) 662 : 50 CONTINUE 663 : AM(K,I,J)=UU*COEFK(KP)*CN 664 : 60 CONTINUE 665 : 70 CONTINUE 666 : 80 CONTINUE 667 : C 668 : C- Matrice masse condensée 669 : C 670 : ELSE 671 : IZAFM = LIZAFM(1,N) 672 : DO 110 K=KDEB,KFIN 673 : KP = K - KDEB + 1 674 : 675 : DO 100 J=1,NP 676 : UU = 0.D0 677 : DO 90 LG=1,NPG 678 : UU = UU + WT(KP,J,LG,N1)*PGSK(KP,LG) 679 : & *DEUPI*RPGK(KP,LG) 680 : 90 CONTINUE 681 : 682 : KJ = MLENTI.LECT(IPT1.NUM(J,K)) 683 : AM(KJ,1,1)=AM(KJ,1,1) + (UU*COEFK(KP)*CN) 684 : 100 CONTINUE 685 : 110 CONTINUE 686 : ENDIF 687 : 120 CONTINUE 688 : 130 CONTINUE 689 : K0 = K0 + NBEL 690 : SEGSUP PETROV 691 : 140 CONTINUE 692 : C 693 : C-------------------------------------------------- 694 : C-- Cas ou l'inconnue est au CENTRE (VF ou EFM1) -- 695 : C-------------------------------------------------- 696 : C 697 : ELSEIF (KPOIND.EQ.2) THEN 698 : DO 160 N=1,NBME 699 : IZAFM = LIZAFM(1,N) 700 : DO 150 I=1,NBEL 701 : NKR = (1-IK1)*(I-1) + 1 702 : AM(I,1,1) = (VPOCHA(I,1) * IZRO.VPOCHA(NKR,1)*CN)/DTR 703 : 150 CONTINUE 704 : 160 CONTINUE 705 : C 706 : C------------------------------------------------------------ 707 : C-- Cas ou l'inconnue est au CENTREP0 (EFMC MACRO ou QUAD) -- 708 : C------------------------------------------------------------ 709 : C 710 : ELSEIF (KPOIND.EQ.3) THEN 711 : DO 163 N=1,NBME 712 : IZAFM = LIZAFM(1,N) 713 : DO 153 I=1,NBEL 714 : NKR = (1-IK1)*(I-1) + 1 715 : AM(I,1,1) = (VPOCHA(I,1) * IZRO.VPOCHA(NKR,1)*CN)/DTR 716 : 153 CONTINUE 717 : 163 CONTINUE 718 : C 719 : C------------------------------------------------------------ 720 : C-- Cas ou l'inconnue est au CENTREP1 (EFMC MACRO ou QUAD) -- 721 : C------------------------------------------------------------ 722 : C 723 : ELSEIF (KPOIND.EQ.4) THEN 724 : 725 : DO 164 N=1,NBME 726 : IZAFM = LIZAFM(1,N) 727 : 728 : SEGACT MELEMO 729 : NBSOUS=MELEMO.LISOUS(/1) 730 : IF(NBSOUS.EQ.0)NBSOUS=1 731 : DO 164 LS=1,NBSOUS 732 : IPT1=MELEMO 733 : IF(NBSOUS.NE.1)IPT1=MELEMO.LISOUS(LS) 734 : SEGACT IPT1 735 : NP1=IPT1.NUM(/1) 736 : NBEL1=IPT1.NUM(/2) 737 : IF(MQUAD.NE.0)THEN 738 : WRITE(NOM0,FMT='(A8)')NOMS(IPT1.ITYPEL)//'PRP1' 739 : ELSEIF(MACRO.NE.0)THEN 740 : WRITE(NOM0,FMT='(A8)')NOMS(IPT1.ITYPEL)//'MCP1' 741 : ENDIF 741 : CALL KALPBG(NOM0,'FONFORM ',IZFFM) 742 : 743 : SEGACT IZFFM*MOD 744 : IZHR=KZHR(1) 745 : SEGACT IZHR*MOD 746 : NES=GR(/1) 747 : NPG=GR(/3) 748 : IZF1=KTP(1) 749 : SEGACT IZF1*MOD 750 : MP=IZF1.FN(/1) 751 : 752 : NKR=1 753 : DO 164 K=1,NBEL 754 : 755 : DO 159 I1=1,NP1 756 : J1=IPT1.NUM(I1,K) 757 : DO 159 ND=1,IDIM 758 : XYZ(ND,I1)=XCOOR((J1-1)*(IDIM+1)+ND) 759 : 159 CONTINUE 759 : CALL CALJBC(FN,GR,PG,XYZ,HR,PGSQ,RPG,NES, 760 : & IDIM,NP1,NPG,IAXI,AIR1) 761 : 762 : DO 154 I=1,MP 763 : DO 155 J=1,MP 764 : U=0.D0 765 : DO 157 L=1,NPG 766 : U=U+IZF1.FN(I,L)*IZF1.FN(J,L)*PGSQ(L) 767 : 157 CONTINUE 768 : AM(K,I,J)=(U*IZRO.VPOCHA(NKR,1) * CN) / DTR 769 : 155 CONTINUE 770 : NKR = NKR + (1-IK1) 771 : 154 CONTINUE 772 : 164 CONTINUE 773 : ENDIF 774 : ENDIF 775 : C 776 : C-------------------------- 777 : C Création du second-membre 778 : C-------------------------- 779 : C Saturation de la matrice calculée ci-dessus par l'inconnue 780 : C au pas de temps précédent 781 : C 782 : C 783 : C- Construction du chapeau de l'objet CHAMPOIN qui contiendra le RHS 784 : C 785 : NAT = 2 786 : NSOUPO = 1 787 : SEGINI MCHPO1 788 : MCHPO1.IFOPOI = IFOMOD 789 : MCHPO1.MOCHDE = TITREE 790 : MCHPO1.MTYPOI = 'SMBR' 791 : MCHPO1.JATTRI(1) = 2 792 : NC = NINKO 793 : SEGINI MSOUP1 794 : MCHPO1.IPCHP(1) = MSOUP1 795 : DO 170 N=1,NINKO 796 : MSOUP1.NOCOMP(N) = LISDUA(N)(1:4) 797 : 170 CONTINUE 798 : MSOUP1.IGEOC = MELEM1 799 : SEGACT MELEM1 800 : N = MELEM1.NUM(/2) 801 : SEGINI MPOVA1 802 : MSOUP1.IPOVAL = MPOVA1 803 : SEGACT MELEM1 804 : NBSOUS = LIZAFM(/1) 805 : C 806 : MELEME = IRIGEL(1,1) 807 : IF(ISCHT.NE.2)THEN 808 : DO 220 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 210 N=1,NINKO 815 : IZAFM = LIZAFM(L,N) 816 : SEGACT IZAFM 817 : DO 200 K=1,NBEL 818 : DO 190 J=1,NP 819 : UU = 0.D0 820 : IU = LECT(IPT1.NUM(J,K)) 821 : DO 180 I=1,NP 822 : IK = IPADI.LECT(IPT1.NUM(I,K)) 823 : UU = UU + AM(K,I,J)*IZTGG2.VPOCHA(IK,N)/CN*CNM 824 : 180 CONTINUE 825 : MPOVA1.VPOCHA(IU,N) = MPOVA1.VPOCHA(IU,N) + UU 826 : 190 CONTINUE 827 : 200 CONTINUE 828 : 210 CONTINUE 829 : 220 CONTINUE 830 : ENDIF 831 : 832 : C Cas BDF2 M Tn-2 833 : IF(ISCHT.EQ.1)THEN 834 : DO 222 L=1,NBSOUS 835 : IPT1 = MELEME 836 : IF (NBSOUS.NE.1) IPT1=LISOUS(L) 837 : SEGACT IPT1 838 : NP = IPT1.NUM(/1) 839 : NBEL = IPT1.NUM(/2) 840 : DO 212 N=1,NINKO 841 : IZAFM = LIZAFM(L,N) 842 : SEGACT IZAFM 843 : DO 202 K=1,NBEL 844 : DO 192 J=1,NP 845 : UU = 0.D0 846 : IU = LECT(IPT1.NUM(J,K)) 847 : DO 182 I=1,NP 848 : IK = IPADI.LECT(IPT1.NUM(I,K)) 849 : UU = UU + AM(K,I,J)*IZTGN2.VPOCHA(IK,N)/CN*CNMM 850 : 182 CONTINUE 851 : MPOVA1.VPOCHA(IU,N) = MPOVA1.VPOCHA(IU,N) + UU 852 : 192 CONTINUE 853 : 202 CONTINUE 854 : 212 CONTINUE 855 : 222 CONTINUE 856 : ENDIF 857 : C 858 : C- Sauvegarde de l'objet MATRIK à l'indice MATELM de la table KIZX et du 859 : C- second-membre à l'indice SMBR de la table EQEX (assemblage éventuel) 860 : C 860 : CALL ECROBJ('MATRIK',MATRIK) 860 : CALL ECROBJ('CHPOINT ',MCHPO1) 861 : 862 : C 863 : C- Désactivation et ménage 864 : C 865 : SEGDES MCHPO1,MSOUP1,MPOVA1 866 : SEGDES IZTGG2 867 : SEGDES IMATRI,MATRIK*NOMOD 868 : IF (IK1DETR.EQ.1) SEGSUP IZRO 869 : IF (IKMUDETR.EQ.1) SEGSUP MZMU 870 : 871 : IF(IPADS.NE.IPADI)SEGSUP IPADS 872 : IF(MLENTI.NE.IPADI)SEGSUP IPADI 873 : SEGSUP IPADI 874 : 875 : RETURN 876 : 1002 FORMAT(10(1X,1PE11.4)) 877 : END 878 : 879 : 880 : 881 : 882 : 883 : 884 : 885 : 886 : 887 : 888 : 889 : 890 : 891 : 892 : 893 : 894 : 895 : 896 : 897 : 898 : 899 : 900 :
© Cast3M 2003 - All rights reserved.
Disclaimer