1 : C DIVS SOURCE MAGN 10/05/18 21:15:21 6675 2 : SUBROUTINE DIVS(IKAS,MTABX) 3 : IMPLICIT INTEGER(I-N) 4 : IMPLICIT REAL*8 (A-H,O-Z) 5 : C************************************************************************ 6 : C Operateur KMAC 7 : C 8 : C OBJET : Cree un objet de type MATRIK 9 : C 10 : C 11 : C IKAS=1 KMAB calcul de B (DIV U) 12 : C IKAS=2 KMBT calcul de Bt uniquement (GRAD P) 13 : C IKAS=3 KBBT calcul de B assemblage pour B et Bt 14 : C 15 : C*********************************************************************** 16 : -INC CCOPTIO 17 : -INC CCGEOME 18 : -INC CCREEL 19 : -INC SMCOORD 20 : -INC SIZFFB 21 : -INC SMCHPOI 22 : POINTEUR IZCH2.MCHPOI,IZCCH2.MPOVAL 23 : POINTEUR IZDV.MCHPOI,IZDDV.MPOVAL,IZTU1.MPOVAL,TETAN.MPOVAL 24 : POINTEUR IZTG1.MCHPOI,IZTGG1.MPOVAL,IZBETA.MPOVAL 25 : POINTEUR IZTGG.MCHPOI,IZANU.MPOVAL 26 : -INC SMMATRIK 27 : 28 : -INC SMLENTI 29 : POINTEUR IZIPAD.MLENTI,MLENTI1.MLENTI,MLENTI2.MLENTI,IPADU.MLENTI 30 : -INC SMLMOTS 31 : POINTEUR LINCO.MLMOTS 32 : -INC SMELEME 33 : POINTEUR MELEMZ.MELEME,MELEMB.MELEME,MELSTB.MELEME 34 : POINTEUR MELEM1.MELEME,MELES1.MELEME,MCTREI.MELEME 35 : POINTEUR IGEOM.MELEME,MELEMA.MELEME 36 : POINTEUR IZLEMC.MELEME,MELEMS.MELEME,MELEMC.MELEME 37 : POINTEUR MELEMI.MELEME,MELEMP.MELEME 38 : POINTEUR MACRO1.MELEME 39 : 40 : CHARACTER*8 TYPE,TYPC,NOMZ,NOMP,NOMD,NOM0 41 : CHARACTER*8 NOMPP,NOMDD 42 : CHARACTER*4 NOM 43 : INTEGER IPAD,IPAD2,IK 44 : REAL*8 KAUX,TETA1 45 : DIMENSION IXV(5) 46 : C 47 : C************************************************************************* 48 : CKMIC 49 : C write(6,*)' Operateur KMIC MTABX=',MTABX 50 : C/MELEMS 51 : 52 : C- Récupération de la table EQEX (pointeur MTAB1) 53 : C 53 : CALL LEKTAB(MTABX,'EQEX',MTAB1) 54 : IF(MTAB1.EQ.0)THEN 55 : C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24 56 : MOTERR( 1: 8) = ' EQEX ' 57 : MOTERR( 9:16) = ' EQEX ' 58 : MOTERR(17:24) = ' KIZX ' 58 : CALL ERREUR(786) 59 : RETURN 60 : ENDIF 60 : CALL ACME(MTAB1,'NAVISTOK',NASTOK) 61 : IF(NASTOK.EQ.0)THEN 61 : CALL ZKMIC(IKAS,MTABX,MTAB1) 62 : RETURN 63 : ENDIF 64 : C 65 : C- Récupération de la table INCO (pointeur KINC) 66 : C 66 : CALL LEKTAB(MTAB1,'INCO',KINC) 67 : IF(KINC.EQ.0)THEN 68 : C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24 69 : MOTERR( 1: 8) = ' INCO ' 70 : MOTERR( 9:16) = ' INCO ' 71 : MOTERR(17:24) = ' EQEX ' 71 : CALL ERREUR(786) 72 : RETURN 73 : ENDIF 74 : 75 : C************************************************************************* 76 : C OPTIONS 77 : C KIMPL = 0 -> EXPL 1 -> IMPL 2 -> CN 78 : C KFORM = 0 -> EFSI 1 -> EF 2 -> VF 3 -> EFMC 79 : C KPRE=3 pression P0 KPRE=4 pression P1 KPRE=2 cas macro 1ère génération 80 : C KPRE=5 pression MSOMMET 81 : 82 : IAXI=0 83 : IK=0 84 : IF(IFOMOD.EQ.0)IAXI=2 85 : C 86 : C- Récupération de la table des options KOPT (pointeur KOPTI) 87 : C 87 : CALL LEKTAB(MTABX,'KOPT',KOPTI) 88 : IF (KOPTI.EQ.0) THEN 89 : C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24 90 : MOTERR( 1: 8) = ' KOPT ' 91 : MOTERR( 9:16) = ' KOPT ' 92 : MOTERR(17:24) = ' KIZX ' 92 : CALL ERREUR(786) 93 : RETURN 94 : ENDIF 95 : 95 : CALL ACME(KOPTI,'KIMPL',KIMPL) 95 : CALL ACME(KOPTI,'KPOIN',KPRE) 95 : CALL ACMF(KOPTI,'AIMPL',TETA1) 95 : CALL ACME(KOPTI,'KFORM',KFORM) 96 : C gbm rajoute option décentrement 18/10/99 96 : CALL ACME(KOPTI,'IDCEN',IDCEN) 97 : C gbm rajoute IKOMP = 2 pour l'opérateur (div) et 98 : C on ne modifie pas IKOMP pour grad 99 : C pour appel ultérieur à SUPGEF. 100 : C IKOMP = 2 veut dire conservatif mais on a div(flux) 101 : C au lieu de div(ro*U). 102 : if (IKAS .EQ. 1 .OR. IKAS .EQ. 3) THEN 103 : IKOMP = 2 104 : endif 105 : * GBM modif 22/11/99 105 : CALL ACMF(KOPTI,'CMD',CMD) 106 : 107 : 108 : IF (IERR.NE.0) RETURN 109 : C write(6,*)' Apres les options ' 110 : C************************************************************************* 111 : C 112 : C- Récupération de la table DOMAINE associée au domaine local 113 : C 113 : CALL ACMM(MTABX,'NOMZONE',NOMZ) 114 : TYPE=' ' 114 : CALL ACMO(MTABX,'DOMZ',TYPE,MMODEL) 115 : IF(TYPE.NE.'MMODEL')THEN 116 : C On attend un des objets : %m1:8 %m9:16 %m17:24 %m25:32 %m33:40 117 : MOTERR( 1: 8) = ' MMODEL ' 118 : MOTERR( 8:16) = ' MMODEL ' 119 : MOTERR(17:24) = ' MMODEL ' 120 : MOTERR(25:32) = ' MMODEL ' 121 : MOTERR(33:40) = ' MMODEL ' 121 : CALL ERREUR(471) 122 : RETURN 123 : ENDIF 124 : 124 : CALL LEKMOD(MMODEL,MTABZ,INEFMD) 125 : 125 : CALL LEKTAB(MTABZ,'MAILLAGE',MELEME) 125 : CALL LEKTAB(MTABZ,'SOMMET',MELEMS) 125 : CALL LEKTAB(MTABZ,'CENTRE',MELEMC) 125 : C?? CALL LEKTAB(MTABZ,'MACRO',MACRO) 126 : MACRO1=0 126 : IF(INEFMD.EQ.2)CALL LEKTAB(MTABZ,'MACRO1',MACRO1) 126 : C?? IF(INEFMD.EQ.3)CALL LEKTAB(MTABZ,'QUADRATI',MQUAD) 127 : IF (IERR.NE.0) RETURN 128 : IF(INEFMD.EQ.4.AND.KPRE.NE.5)THEN 129 : C% Données incompatibles 129 : CALL ERREUR(21) 130 : RETURN 131 : ENDIF 132 : 133 : MELEMI=MELEME 134 : IF(MACRO1.NE.0.AND.KPRE.NE.2)THEN 135 : MELEMI=MELEME 136 : ENDIF 137 : 138 : IF(KPRE.EQ.2.AND.INEFMD.EQ.3)KPRE=3 139 : IF(INEFMD.EQ.1)KPRE=2 140 : 141 : * ON MODIFIE POUR SOMMET -> SOMMET GBM 18/10/99 142 : MELEMP=MELEME 142 : CALL LEKTAB(MTABZ,'CENTRE',MELEMC) 143 : 144 : C************************************************************************* 145 : C VERIFICATIONS SUR LES INCONNUES 146 : 147 : C write(6,*)' Verification sur les inconnues ' 148 : TYPE='LISTMOTS' 148 : CALL ACMO(MTABX,'LISTINCO',TYPE,LINCO) 149 : IF (IERR.NE.0) RETURN 150 : SEGACT LINCO 151 : NBINC=LINCO.MOTS(/2) 152 : IF(NBINC.NE.2)THEN 153 : WRITE(6,*)'Operateur KMAC ' 154 : WRITE(6,*)'Nombre d''inconnues incorrecte : ',NBINC,' On attend 2' 155 : C Indice %m1:8 : contient plus de %i1 %m9:16 156 : MOTERR( 1:8) = 'LISTINCO' 157 : INTERR(1) = 2 158 : MOTERR(9:16) = ' MOTS ' 158 : CALL ERREUR(799) 159 : RETURN 160 : ENDIF 161 : 162 : C On recupere PHI n et TETA n pour Cranck-Nicholson 163 : NOMP=LINCO.MOTS(1) 164 : TYPE=' ' 164 : CALL ACMO(KINC,NOMP,TYPE,MCHPOI) 165 : IF(TYPE.NE.'CHPOINT ')THEN 166 : WRITE(6,*)' Opérateur KMAC :' 167 : WRITE(6,*)' L objet CHPOINT ',NOMP, 168 : & ' n existe pas dans la table' 169 : C Indice %m1:8 : ne contient pas un objet de type %m9:16 170 : MOTERR( 1: 8) = 'INC '//NOMP 171 : MOTERR( 9:16) = 'CHPOINT ' 171 : CALL ERREUR(800) 172 : RETURN 173 : ELSE 173 : CALL LICHT(MCHPOI,IZTU1,TYPC,IGEOM0) 174 : ENDIF 175 : C************************************************************************* 176 : C Le domaine de definition est donne par le SPG de la premiere inconnue 177 : C Les inconnues suivantes devront posseder ce meme pointeur 178 : C On verifie que les points de la zone sont tous inclus dans ce SPG 179 : C Inconnue Primale 180 : 181 : C write(6,*)' Verification inconnue primale ' 181 : CALL KRIPAD(IGEOM0,IPADU) 182 : IF(IKAS.EQ.1.OR.IKAS.EQ.3)THEN 183 : MELEMK=MELEMS 184 : ELSE 185 : MELEMK=MELEMC 186 : ENDIF 187 : 187 : CALL VERPAD(IPADU,MELEMK,IRET) 188 : IF(IRET.NE.0)THEN 189 : WRITE(6,*)' Opérateur KMAC ' 190 : WRITE(6,*)' La zone ',NOMZ,' n''est pas incluse dans le domaine' 191 : & , ' de définition de l''inconnue ',NOMP 192 : WRITE(6,*)' MELEMK=',melemk,' IGEOM0=',IGEOM0 193 : C Indice %m1:8 : L'objet %m9:16 n'a pas le bon support géométrique 194 : MOTERR(1: 8) = 'INC '//NOMP 195 : MOTERR(9:16) = 'CHPOINT ' 195 : CALL ERREUR(788) 196 : IPAS=0 197 : RETURN 198 : ENDIF 199 : 200 : 201 : C************************************************************************* 202 : 203 : NOMD=LINCO.MOTS(2) 204 : TYPE=' ' 204 : CALL ACMO(KINC,NOMD,TYPE,MCHPOI) 205 : IF(TYPE.NE.'CHPOINT ')THEN 206 : WRITE(6,*)' Opérateur KMAC :' 207 : WRITE(6,*)' L objet CHPOINT ',NOMD, 208 : & ' n existe pas dans la table' 209 : C Indice %m1:8 : ne contient pas un objet de type %m9:16 210 : MOTERR( 1: 8) = 'INC '//NOMD 211 : MOTERR( 9:16) = 'CHPOINT ' 211 : CALL ERREUR(800) 212 : RETURN 213 : ELSE 213 : CALL LICHT(MCHPOI,TETAN,TYPC,IGEOM0) 214 : ENDIF 215 : 216 : NC=TETAN.VPOCHA(/2) 217 : C************************************************************************* 218 : C Le domaine de definition est donne par le SPG de la premiere inconnue 219 : C Les inconnues suivantes devront posseder ce meme pointeur 220 : C On verifie que les points de la zone sont tous inclus dans ce SPG 221 : C Inconnue Duale 222 : 223 : C write(6,*)' IGEOM0=',igeom0 223 : CALL KRIPAD(IGEOM0,MLENTI) 224 : IF(IKAS.EQ.1.OR.IKAS.EQ.3)THEN 225 : C ON MODIFIE POUR CAS MEME ESPACE DIV U SOMMET -> SOMMET GBM 18/10 226 : MELEMK=MELEMS 227 : ELSE 228 : MELEMK=MELEMK 229 : ENDIF 230 : 231 : C write(6,*)' Verification inconnue duale ',MELEMK 231 : CALL VERPAD(MLENTI,MELEMK,IRET) 232 : IF(IRET.NE.0)THEN 233 : WRITE(6,*)' Opérateur KMAC ' 234 : WRITE(6,*)' La zone ',NOMZ,' n''est pas incluse dans le domaine' 235 : & ,' de définition de l''inconnue ',NOMD 236 : WRITE(6,*)' MELEMK=',melemk,' IGEOM0=',IGEOM0 237 : C Indice %m1:8 : L'objet %m9:16 n'a pas le bon support géométrique 238 : MOTERR(1: 8) = 'INC '//NOMD 239 : MOTERR(9:16) = 'CHPOINT ' 239 : CALL ERREUR(788) 240 : IPAS=0 241 : RETURN 242 : ENDIF 243 : 244 : C GBM rajoute NPTU 18/10 245 : NPTU = IZTU1.VPOCHA(/1) 246 : 247 : C************************************************************************* 248 : C Lecture du ou des coefficients 249 : C Type du coefficient : 250 : C IK1=0 CHPOINT IK1=1 scalaire IK1=2 vecteur 251 : 252 : C write(6,*)' Verification sur les coefficients ' 252 : CALL ACME(MTABX,'IARG',IARG) 253 : C write(6,*)' IARG=',iarg 254 : 255 : IF(MACRO1.NE.0.AND.IKAS.NE.2.AND.KPRE.EQ.2)THEN 256 : TYPE=' ' 256 : CALL ACMO(MTABZ,'MELSTB',TYPE,MELSTB) 257 : SEGACT MELSTB 258 : IF(IDIM.EQ.2)NBELEM=MELSTB.NUM(/2)/4 259 : IF(IDIM.EQ.3)NBELEM=MELSTB.NUM(/2)/8 260 : NBNN=MELSTB.NUM(/1) 261 : NBSOUS=0 262 : NBREF=0 263 : SEGINI MELEMA 264 : MELEMA.ITYPEL=MELSTB.ITYPEL 265 : 266 : NKPE=4 267 : IF(IDIM.EQ.3)NKPE=8 268 : do 4878 k=1,nbelem 269 : mi=(k-1)*NKPE+1 270 : do 4879 i=1,nbnn 271 : MELEMA.num(i,k)=melstb.num(i,mi) 272 : 4879 continue 273 : C write(6,*)k,(MELEMA.num(i,k),i=1,nbnn) 274 : 4878 continue 275 : 276 : TYPE=' ' 276 : CALL ACMO(MTABZ,'MCHPOC',TYPE,MCHPOC) 277 : TYPE=' ' 277 : CALL ACMO(MTABZ,'CENTRE',TYPE,MCTREI) 278 : ENDIF 279 : 280 : 281 : C 1er COEF 282 : 283 : IXV(1)=MELEMC 284 : IXV(2)=1 285 : IXV(3)=0 286 : IXV(4)=MELEMS 287 : IXV(5)=-MELEMS 288 : IRET=5 288 : CALL LEKCOF('Opérateur KMAC :', 289 : & MTABX,KINC,1,IXV,IZTG1,IZTGG1,NPT1,NC1,IK1,IRET) 290 : IF(IRET.EQ.0)RETURN 290 : IF(IK1.GE.4)CALL KRIPAD(MELEMS,MLENTI) 291 : 292 : 293 : BETA0=1.D-1 294 : IF(IARG.EQ.2)THEN 295 : C 2ème COEF 296 : IXV(1)=0 297 : IXV(2)=1 298 : IXV(3)=0 298 : CALL LEKCOF('Opérateur KMAC :', 299 : & MTABX,KINC,2,IXV,IZTG2,IZBETA,NPT2,NC2,IK2,IRET) 300 : IF(IRET.EQ.0)RETURN 301 : BETA0=IZBETA.VPOCHA(1,1) 302 : ENDIF 303 : 304 : C GBM rajoute 3eme coef mot clef 'SOMM' 305 : 306 : C GBM rajoute 4eme coef viscosité pour SUPGEF 307 : 308 : ANUK=1.D-15 309 : IF(IARG.EQ.4)THEN 310 : C 4ème COEF 311 : IXV(1)=0 312 : IXV(2)=1 313 : IXV(3)=0 313 : CALL LEKCOF('Opérateur KMAC :', 314 : & MTABX,KINC,4,IXV,IZTGG,IZANU,NPT3,NC3,IK3,IRET) 315 : IF(IRET.EQ.0)RETURN 316 : ANUK=IZANU.VPOCHA(1,1) 317 : ENDIF 318 : 319 : NOMP=LINCO.MOTS(1) 320 : NOMD=LINCO.MOTS(2)(1:4) 321 : 322 : NRIGE=7 323 : NKID =9 324 : NKMT =7 325 : NMATRI=1 326 : IF(MACRO1.NE.0.AND.IKAS.NE.2.AND.KPRE.EQ.2)NMATRI=2 327 : IF(MACRO1.NE.0.AND.IKAS.NE.2.AND.KPRE.EQ.4.AND.IDIM.EQ.2)NMATRI=2 328 : SEGINI MATRIK 329 : 330 : C Fin Stabilisation de toutes sortes 331 : C write(6,*)'C Fin Stabilisation de toutes sortes' 332 : 333 : NBME=IDIM 333 : CALL KRIPAD(MELEMI,MLENTI1) 334 : SEGACT MELEMI 335 : NBSOUS=MELEMI.LISOUS(/1) 336 : IF(NBSOUS.EQ.0)NBSOUS=1 337 : SEGINI IMATRI 338 : 339 : IF(IKAS.EQ.2)THEN 340 : KSPGD=MELEMS 341 : KSPGP=MELEMC 342 : IRIGEL(2,1)=MELEME 343 : IRIGEL(1,1)=MELEME 344 : ELSE 345 : KSPGP=MELEMS 346 : KSPGD=MELEMS 347 : IRIGEL(1,1)=MELEME 348 : IRIGEL(2,1)=MELEME 349 : ENDIF 350 : SEGACT MELEMP 351 : 352 : C write(6,*)' ds kmac melemp=',IRIGEL(1,1) 353 : C write(6,*)' ds kmac melemd=',IRIGEL(2,1) 354 : 355 : IRIGEL(4,1)=IMATRI 356 : * je modifie celui là pour retirer les multiplicateurs 357 : IF(IKAS.EQ.1)IRIGEL(7,1)= 5 358 : IF(IKAS.EQ.2)IRIGEL(7,1)=5 359 : IF(IKAS.EQ.3)IRIGEL(7,1)=5 360 : 361 : K0=0 362 : DO 11 L=1,NBSOUS 363 : IPT1=MELEMI 364 : IF(NBSOUS.NE.1)IPT1=MELEMI.LISOUS(L) 365 : SEGACT IPT1 366 : NBEL=IPT1.NUM(/2) 367 : 368 : IF(IKAS.EQ.2)THEN 369 : MP=IPT1.NUM(/1) 370 : NP=MELEMP.NUM(/1) 371 : ELSE 372 : NP=IPT1.NUM(/1) 373 : MP=MELEMP.NUM(/1) 374 : ENDIF 375 : 376 : DO 12 I=1,NBME 377 : SEGINI IZAFM 378 : C write(6,*)' ni izafm np=',np,' mp=',mp,' nbel=',nbel,izafm,l,i 379 : LIZAFM(L,I)=IZAFM 380 : IF(IKAS.EQ.2)THEN 381 : WRITE(NOM,FMT='(I1,A3)')I,NOMD(1:3) 382 : LISDUA(I)=NOM//' ' 383 : LISPRI(I)=NOMP 384 : ELSE 385 : WRITE(NOM,FMT='(I1,A3)')I,NOMP(1:3) 386 : LISPRI(I)=NOM//' ' 387 : LISDUA(I)=NOMD 388 : ENDIF 389 : 12 CONTINUE 390 : IPM1=LIZAFM(L,1) 391 : IPM2=LIZAFM(L,2) 392 : IPM3=LIZAFM(L,2) 393 : IF(IDIM.EQ.3)IPM3=LIZAFM(L,3) 394 : 395 : IF (IKAS .EQ. 1 .OR. IKAS .EQ. 3) THEN 396 : c on code div roU -> T, 397 : IF(IK1.LT.4)THEN 398 : c cas où le coef multiplicateur est scalaire, point 399 : c ou champ au centre. 399 : CALL KSPRUS 400 : & (IPT1,IPM1,IPM2,IPM3,IAXI,IKAS,INEFMD,KPRE,IZTGG1,IK1,KO, 401 : & ANUK,IDCEN,IKOMP,NPTU,MLENTI,IZTU1,TETAN,NP,NBEL,CMD) 402 : ELSE 403 : C MODIFé par GBM, cas coef au sommet. 403 : CALL KSPRJS 404 : & (IPT1,IPM1,IPM2,IPM3,IAXI,IKAS,INEFMD,KPRE,IZTGG1 405 : & ,IPADU,IK1) 406 : ENDIF 407 : ELSE 408 : c on code gradT -> U 409 : c on inverse la place de IZTU1 et de TETAN dans les arguments 410 : c ca l'utilisateur rentre la pression en premier, IZTU1 est donc 411 : c une pression (nom mal choisi). GBM 18/10/99 412 : IF(IK1.LT.4)THEN 413 : c cas où le coef multiplicateur est scalaire, point 414 : c ou champ au centre. 414 : CALL KSPRUS 415 : & (IPT1,IPM1,IPM2,IPM3,IAXI,IKAS,INEFMD,KPRE,IZTGG1,IK1,KO, 416 : & ANUK,IDCEN,IKOMP,NPTU,MLENTI,TETAN,IZTU1,NP,NBEL,CMD) 417 : ELSE 418 : C MODIFé par GBM, cas coef au sommet. 418 : CALL KSPRJS 419 : & (IPT1,IPM1,IPM2,IPM3,IAXI,IKAS,INEFMD,KPRE,IZTGG1, 420 : & IPADU,IK1) 421 : ENDIF 422 : ENDIF 423 : C 424 : K0=K0+NBEL 425 : SEGDES IPM1,IPM2,IPM3 426 : SEGDES IPT1 427 : 11 CONTINUE 428 : SEGDES MELEMI 429 : IF(IK1.EQ.4)SEGSUP MLENTI 430 : 430 : CALL ECROBJ('MATRIK',MATRIK) 431 : 432 : NAT=2 433 : NSOUPO=0 434 : SEGINI MCHPOI 435 : JATTRI(1)=2 435 : CALL ECROBJ('CHPOINT',MCHPOI) 436 : 437 : C write(6,*)' Fin operateur KMIC' 438 : SEGDES IMATRI,MATRIK 439 : RETURN 440 : 1001 FORMAT(20(1X,I5)) 441 : END 442 : 443 : 444 : 445 : 446 : 447 : 448 : 449 : 450 :
© Cast3M 2003 - All rights reserved.
Disclaimer