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