1 : C KMIC SOURCE MAGN 10/05/18 21:15:44 6675 2 : SUBROUTINE KMIC(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,VISCO.MPOVAL 25 : -INC SMMATRIK 26 : POINTEUR IPM.IZAFM 27 : SEGMENT IMATRS 28 : ENDSEGMENT 29 : POINTEUR IPMS.IZAFM,IPS1.IZAFM,IPS2.IZAFM,IPS3.IZAFM 30 : 31 : PARAMETER (LRV=64) 32 : DIMENSION VISC(3) 33 : 34 : -INC SMLENTI 35 : POINTEUR IZIPAD.MLENTI,MLENTI1.MLENTI,MLENTI2.MLENTI 36 : -INC SMLMOTS 37 : POINTEUR LINCO.MLMOTS 38 : -INC SMELEME 39 : POINTEUR MELEMZ.MELEME,MELEMB.MELEME,MELSTB.MELEME 40 : POINTEUR MELEM1.MELEME,MELES1.MELEME,MCTREI.MELEME 41 : POINTEUR IGEOM.MELEME,MELEMA.MELEME 42 : POINTEUR IZLEMC.MELEME,MELEMS.MELEME,MELEMC.MELEME 43 : POINTEUR MELEMI.MELEME,MELEMP.MELEME 44 : POINTEUR MACRO1.MELEME 45 : 46 : CHARACTER*8 TYPE,TYPC,NOMZ,NOMP,NOMD,NOM0,NOMK 47 : CHARACTER*8 NOMPP,NOMDD 48 : CHARACTER*4 NOM 49 : INTEGER IPAD,IPAD2,IK 50 : REAL*8 KAUX,TETA1 51 : DIMENSION IXV(5) 52 : C 53 : C************************************************************************* 54 : CKMIC 55 : C write(6,*)' Operateur KMIC MTABX=',MTABX 56 : C 57 : C- Récupération de la table EQEX (pointeur MTAB1) 58 : C 58 : CALL LEKTAB(MTABX,'EQEX',MTAB1) 59 : IF(MTAB1.EQ.0)THEN 60 : C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24 61 : MOTERR( 1: 8) = ' EQEX ' 62 : MOTERR( 9:16) = ' EQEX ' 63 : MOTERR(17:24) = ' KIZX ' 63 : CALL ERREUR(786) 64 : RETURN 65 : ENDIF 66 : 66 : CALL ACME(MTABX,'IARG',IARG) 67 : IF(IARG.GE.3)THEN 68 : C 3ème COEF 68 : CALL ACMM(MTABX,'ARG3',NOMK) 69 : 70 : IF(NOMK.EQ.'SOMM')THEN 71 : * Cet opérateur calcul un gradient ou div de sommets vers sommets 72 : * sans contrainte de type multiplicateur de Lagrange. 72 : CALL DIVS(IKAS,MTABX) 73 : RETURN 74 : * 75 : ELSEIF(NOMK.EQ.'NOMU')THEN 76 : * Cet opérateur calcule la divergence (ou grad) de sommets vers ctp1 77 : * sans contrainte de cmultiplicateurs de Lagrange. 77 : CALL DIVC(IKAS,MTABX) 78 : C write(6,*)' RETOUR KKMIC' 79 : RETURN 80 : 81 : ENDIF 82 : ENDIF 83 : 83 : CALL ACME(MTAB1,'DISCPRES',IDISCP) 84 : 84 : CALL ACME(MTAB1,'NAVISTOK',NASTOK) 85 : IF(NASTOK.EQ.0)THEN 85 : CALL ZKMIC(IKAS,MTABX,MTAB1) 86 : RETURN 87 : ENDIF 88 : C 89 : C- Récupération de la table INCO (pointeur KINC) 90 : C 90 : CALL LEKTAB(MTAB1,'INCO',KINC) 91 : IF(KINC.EQ.0)THEN 92 : C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24 93 : MOTERR( 1: 8) = ' INCO ' 94 : MOTERR( 9:16) = ' INCO ' 95 : MOTERR(17:24) = ' EQEX ' 95 : CALL ERREUR(786) 96 : RETURN 97 : ENDIF 98 : 99 : C************************************************************************* 100 : C OPTIONS 101 : C KIMPL = 0 -> EXPL 1 -> IMPL 2 -> CN 102 : C KFORM = 0 -> EFSI 1 -> EF 2 -> VF 3 -> EFMC 103 : C KPRE=3 pression P0 KPRE=4 pression P1 KPRE=2 cas macro 1ère génération 104 : C KPRE=5 pression MSOMMET 105 : 106 : NOMUL=1 107 : IAXI=0 108 : IK=0 109 : IF(IFOMOD.EQ.0)IAXI=2 110 : C 111 : C- Récupération de la table des options KOPT (pointeur KOPTI) 112 : C 112 : CALL LEKTAB(MTABX,'KOPT',KOPTI) 113 : IF (KOPTI.EQ.0) THEN 114 : C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24 115 : MOTERR( 1: 8) = ' KOPT ' 116 : MOTERR( 9:16) = ' KOPT ' 117 : MOTERR(17:24) = ' KIZX ' 117 : CALL ERREUR(786) 118 : RETURN 119 : ENDIF 120 : 120 : CALL ACME(KOPTI,'KIMPL',KIMPL) 120 : CALL ACME(KOPTI,'KPOIN',KPRE) 121 : C write(6,*)'KBBT KPRE=',KPRE 121 : CALL ACMF(KOPTI,'AIMPL',TETA1) 121 : CALL ACME(KOPTI,'KFORM',KFORM) 121 : C CALL ACME(KOPTI,'IKOMP',IKOMP) 122 : 123 : IF (IERR.NE.0) RETURN 124 : C write(6,*)' Apres les options ' 125 : C************************************************************************* 126 : C 127 : C- Récupération de la table DOMAINE associée au domaine local 128 : C 128 : CALL ACMM(MTABX,'NOMZONE',NOMZ) 129 : TYPE=' ' 129 : CALL ACMO(MTABX,'DOMZ',TYPE,MMODEL) 130 : IF(TYPE.NE.'MMODEL')THEN 131 : C On attend un des objets : %m1:8 %m9:16 %m17:24 %m25:32 %m33:40 132 : MOTERR( 1: 8) = ' MMODEL ' 133 : MOTERR( 8:16) = ' MMODEL ' 134 : MOTERR(17:24) = ' MMODEL ' 135 : MOTERR(25:32) = ' MMODEL ' 136 : MOTERR(33:40) = ' MMODEL ' 136 : CALL ERREUR(471) 137 : RETURN 138 : ENDIF 139 : 139 : CALL LEKMOD(MMODEL,MTABZ,INEFMD) 140 : 140 : CALL LEKTAB(MTABZ,'MAILLAGE',MELEME) 140 : CALL LEKTAB(MTABZ,'SOMMET',MELEMS) 140 : CALL LEKTAB(MTABZ,'CENTRE',MELEMC) 140 : C?? CALL LEKTAB(MTABZ,'MACRO',MACRO) 141 : MACRO1=0 141 : IF(INEFMD.EQ.2)CALL LEKTAB(MTABZ,'MACRO1',MACRO1) 141 : C?? IF(INEFMD.EQ.3)CALL LEKTAB(MTABZ,'QUADRATI',MQUAD) 142 : IF (IERR.NE.0) RETURN 143 : IF(INEFMD.EQ.4.AND.KPRE.NE.5)THEN 144 : C% Données incompatibles 144 : CALL ERREUR(21) 145 : RETURN 146 : ENDIF 147 : 148 : MELEMI=MELEME 149 : IF(MACRO1.NE.0.AND.KPRE.NE.2)THEN 150 : MELEMI=MACRO1 151 : ENDIF 152 : 153 : IF(KPRE.EQ.2.AND.INEFMD.EQ.3)KPRE=3 154 : IF(INEFMD.EQ.1.AND.KPRE.NE.5)KPRE=2 155 : 156 : IF(KPRE.EQ.3)THEN 156 : CALL LEKTAB(MTABZ,'CENTREP0',MELEMC) 157 : MELEMP=MELEMC 158 : ELSEIF(KPRE.EQ.4)THEN 158 : CALL LEKTAB(MTABZ,'CENTREP1',MELEMC) 158 : CALL LEKTAB(MTABZ,'ELTP1NC ',MELEMP) 159 : ELSEIF(KPRE.EQ.2)THEN 159 : CALL LEKTAB(MTABZ,'CENTRE',MELEMC) 160 : MELEMP=MELEMC 161 : ELSEIF(KPRE.EQ.5)THEN 161 : CALL LEKTAB(MTABZ,'MSOMMET',MELEMC) 161 : CALL LEKTAB(MTABZ,'MMAIL ',MELEMP) 162 : ENDIF 163 : 164 : IF(IDISCP.EQ.0)THEN 164 : CALL ECME(MTAB1,'DISCPRES',KPRE) 165 : ELSEIF(IDISCP.NE.KPRE)THEN 166 : C% Données incompatibles 166 : CALL ERREUR(21) 167 : RETURN 168 : ENDIF 169 : 170 : IKOMP=1 171 : C IF(KPRE.EQ.5)IKOMP=0 172 : C IF(KPRE.EQ.5)IKOMP=1 173 : 174 : C************************************************************************* 175 : C VERIFICATIONS SUR LES INCONNUES 176 : 177 : C write(6,*)' Verification sur les inconnues ' 178 : TYPE='LISTMOTS' 178 : CALL ACMO(MTABX,'LISTINCO',TYPE,LINCO) 179 : IF (IERR.NE.0) RETURN 180 : SEGACT LINCO 181 : NBINC=LINCO.MOTS(/2) 182 : IF(NBINC.NE.2)THEN 183 : WRITE(6,*)'Operateur KMAB ' 184 : WRITE(6,*)'Nombre d''inconnues incorrecte : ',NBINC,' On attend 2' 185 : C Indice %m1:8 : contient plus de %i1 %m9:16 186 : MOTERR( 1:8) = 'LISTINCO' 187 : INTERR(1) = 2 188 : MOTERR(9:16) = ' MOTS ' 188 : CALL ERREUR(799) 189 : RETURN 190 : ENDIF 191 : 192 : C On recupere PHI n et TETA n pour Cranck-Nicholson 193 : NOMP=LINCO.MOTS(1) 194 : TYPE=' ' 194 : CALL ACMO(KINC,NOMP,TYPE,MCHPOI) 195 : IF(TYPE.NE.'CHPOINT ')THEN 196 : WRITE(6,*)' Opérateur KMAB :' 197 : WRITE(6,*)' L objet CHPOINT ',NOMP, 198 : & ' n existe pas dans la table' 199 : C Indice %m1:8 : ne contient pas un objet de type %m9:16 200 : MOTERR( 1: 8) = 'INC '//NOMP 201 : MOTERR( 9:16) = 'CHPOINT ' 201 : CALL ERREUR(800) 202 : RETURN 203 : ELSE 203 : CALL LICHT(MCHPOI,IZTU1,TYPC,IGEOM0) 204 : ENDIF 205 : C************************************************************************* 206 : C Le domaine de definition est donne par le SPG de la premiere inconnue 207 : C Les inconnues suivantes devront posseder ce meme pointeur 208 : C On verifie que les points de la zone sont tous inclus dans ce SPG 209 : C Inconnue Primale 210 : 211 : C write(6,*)' Verification inconnue primale ' 211 : CALL KRIPAD(IGEOM0,MLENTI) 212 : IF(IKAS.EQ.1.OR.IKAS.EQ.3)THEN 213 : MELEMK=MELEMS 214 : ELSE 215 : MELEMK=MELEMC 216 : ENDIF 217 : 218 : C write(6,*)' MELEMK=',melemk 218 : CALL VERPAD(MLENTI,MELEMK,IRET) 219 : IF(IRET.NE.0)THEN 220 : WRITE(6,*)' Opérateur KMAB ' 221 : WRITE(6,*)' La zone ',NOMZ,' n''est pas incluse dans le domaine' 222 : & , ' de définition de l''inconnue ',NOMP 223 : WRITE(6,*)' MELEMK=',melemk,' IGEOM0=',IGEOM0 224 : C Indice %m1:8 : L'objet %m9:16 n'a pas le bon support géométrique 225 : MOTERR(1: 8) = 'INC '//NOMP 226 : MOTERR(9:16) = 'CHPOINT ' 226 : CALL ERREUR(788) 227 : IPAS=0 228 : RETURN 229 : ENDIF 230 : 231 : SEGSUP MLENTI 232 : 233 : C************************************************************************* 234 : 235 : NOMD=LINCO.MOTS(2) 236 : TYPE=' ' 236 : CALL ACMO(KINC,NOMD,TYPE,MCHPOI) 237 : IF(TYPE.NE.'CHPOINT ')THEN 238 : WRITE(6,*)' Opérateur KMAB :' 239 : WRITE(6,*)' L objet CHPOINT ',NOMD, 240 : & ' n existe pas dans la table' 241 : C Indice %m1:8 : ne contient pas un objet de type %m9:16 242 : MOTERR( 1: 8) = 'INC '//NOMD 243 : MOTERR( 9:16) = 'CHPOINT ' 243 : CALL ERREUR(800) 244 : RETURN 245 : ELSE 245 : CALL LICHT(MCHPOI,TETAN,TYPC,IGEOM0) 246 : ENDIF 247 : 248 : NC=TETAN.VPOCHA(/2) 249 : C************************************************************************* 250 : C Le domaine de definition est donne par le SPG de la premiere inconnue 251 : C Les inconnues suivantes devront posseder ce meme pointeur 252 : C On verifie que les points de la zone sont tous inclus dans ce SPG 253 : C Inconnue Duale 254 : 255 : C write(6,*)' IGEOM0=',igeom0 255 : CALL KRIPAD(IGEOM0,MLENTI) 256 : IF(IKAS.EQ.1.OR.IKAS.EQ.3)THEN 257 : MELEMK=MELEMC 258 : ELSE 259 : MELEMK=MELEMS 260 : ENDIF 261 : 262 : C write(6,*)' Verification inconnue duale ',MELEMK 262 : CALL VERPAD(MLENTI,MELEMK,IRET) 263 : IF(IRET.NE.0)THEN 264 : WRITE(6,*)' Opérateur KMAB ' 265 : WRITE(6,*)' La zone ',NOMZ,' n''est pas incluse dans le domaine' 266 : & ,' de définition de l''inconnue ',NOMD 267 : WRITE(6,*)' MELEMK=',melemk,' IGEOM0=',IGEOM0 268 : C Indice %m1:8 : L'objet %m9:16 n'a pas le bon support géométrique 269 : MOTERR(1: 8) = 'INC '//NOMD 270 : MOTERR(9:16) = 'CHPOINT ' 270 : CALL ERREUR(788) 271 : IPAS=0 272 : RETURN 273 : ENDIF 274 : 275 : SEGSUP MLENTI 276 : 277 : C************************************************************************* 278 : C Lecture du ou des coefficients 279 : C Type du coefficient : 280 : C IK1=0 CHPOINT IK1=1 scalaire IK1=2 vecteur 281 : 282 : C write(6,*)' Verification sur les coefficients ' 282 : CALL ACME(MTABX,'IARG',IARG) 283 : C write(6,*)' IARG=',iarg 284 : 285 : IF(MACRO1.NE.0.AND.IKAS.NE.2.AND.KPRE.EQ.2)THEN 286 : TYPE=' ' 286 : CALL ACMO(MTABZ,'MELSTB',TYPE,MELSTB) 287 : SEGACT MELSTB 288 : IF(IDIM.EQ.2)NBELEM=MELSTB.NUM(/2)/4 289 : IF(IDIM.EQ.3)NBELEM=MELSTB.NUM(/2)/8 290 : NBNN=MELSTB.NUM(/1) 291 : NBSOUS=0 292 : NBREF=0 293 : SEGINI MELEMA 294 : MELEMA.ITYPEL=MELSTB.ITYPEL 295 : 296 : NKPE=4 297 : IF(IDIM.EQ.3)NKPE=8 298 : do 4878 k=1,nbelem 299 : mi=(k-1)*NKPE+1 300 : do 4879 i=1,nbnn 301 : MELEMA.num(i,k)=melstb.num(i,mi) 302 : 4879 continue 303 : C write(6,*)k,(MELEMA.num(i,k),i=1,nbnn) 304 : 4878 continue 305 : 306 : TYPE=' ' 306 : CALL ACMO(MTABZ,'MCHPOC',TYPE,MCHPOC) 307 : TYPE=' ' 307 : CALL ACMO(MTABZ,'CENTRE',TYPE,MCTREI) 308 : ENDIF 309 : 310 : 311 : C 1er COEF 312 : 313 : IXV(1)=MELEMC 314 : IXV(2)=1 315 : IXV(3)=0 316 : IXV(4)=MELEMS 317 : IXV(5)=-MELEMS 318 : IRET=5 318 : CALL LEKCOF('Opérateur KMAB :', 319 : & MTABX,KINC,1,IXV,IZTG1,IZTGG1,NPT1,NC1,IK1,IRET) 320 : IF(IRET.EQ.0)RETURN 320 : IF(IK1.GE.4)CALL KRIPAD(MELEMS,MLENTI) 321 : 322 : 323 : BETA0=1.D-1 324 : IF(IARG.EQ.2)THEN 325 : C 2ème COEF 326 : IXV(1)=0 327 : IXV(2)=1 328 : IXV(3)=0 328 : CALL LEKCOF('Opérateur KMAB :', 329 : & MTABX,KINC,2,IXV,IZTG2,IZBETA,NPT2,NC2,IK2,IRET) 330 : IF(IRET.EQ.0)RETURN 331 : BETA0=IZBETA.VPOCHA(1,1) 332 : C? IF(IZBETA.VPOCHA(1,1).LT.0.D0)THEN 333 : C% Données incompatibles 333 : C? CALL ERREUR(21) 334 : C? RETURN 335 : C? ENDIF 336 : ENDIF 337 : 338 : IF(IARG.EQ.3)THEN 339 : C 3ème COEF 339 : CALL ACMM(MTABX,'ARG3',NOMK) 340 : C write(6,*)' NOMK=',nomk 341 : IF(NOMK.EQ.'CAVIFERM')IKOMP=0 342 : IF(NOMK.EQ.'NOMUL ')NOMUL=0 343 : ENDIF 344 : 345 : NOMP=LINCO.MOTS(1) 346 : NOMD=LINCO.MOTS(2)(1:4) 347 : 348 : NRIGE=7 349 : NKID =9 350 : NKMT =7 351 : NMATRI=1 352 : IF(MACRO1.NE.0.AND.IKAS.NE.2.AND.KPRE.EQ.2)NMATRI=2 353 : IF(MACRO1.NE.0.AND.IKAS.NE.2.AND.KPRE.EQ.4.AND.IDIM.EQ.2)NMATRI=2 354 : IF(KPRE.EQ.5)NMATRI=NMATRI+1 355 : C write(6,*)' NMATRI=',nmatri,' KPRE=',kpre 356 : SEGINI MATRIK 357 : 358 : C CAS Stabilisation via MACRO CENTRE 359 : IF(MACRO1.NE.0.AND.IKAS.NE.2.AND.KPRE.EQ.2)THEN 360 : IF(IKAS.EQ.3)BETA0=-BETA0 361 : C write(6,*)' CAS Stabilisation via MACRO CENTRE ',BETA0,IKAS 362 : NK=0 363 : I2=2 364 : NBME=1 365 : NBSOUS=1 366 : SEGINI IMATRI 367 : IRIGEL(4,i2)=IMATRI 368 : KSPGP=MCTREI 369 : KSPGD=MCTREI 370 : IRIGEL(1,i2)=MELEMA 371 : IRIGEL(2,i2)=MELEMA 372 : IRIGEL(7,i2)=0 372 : CALL LICHT(MCHPOC,MPOVAL,TYPC,IGEOM) 373 : 374 : C write(6,*)' MELSTB=',melstb,' MCTREI=',MCTREI 375 : 376 : SEGACT MELSTB 377 : 378 : NBSOUS=MELSTB.LISOUS(/1) 379 : IF(NBSOUS.NE.0)THEN 379 : CALL ERREUR(5) 380 : ENDIF 381 : 382 : NBEL=MELEMA.NUM(/2) 383 : NBCI=MELSTB.NUM(/2) 384 : NP =MELSTB.NUM(/1) 385 : MP =NP 386 : 387 : C write(6,*)' nbel,nbci,np,mp=',nbel,nbci,np,mp 388 : SEGINI IZAFM 389 : LIZAFM(1,1)=IZAFM 390 : LISPRI(1)=NOMD 391 : LISDUA(1)=NOMD 392 : 392 : CALL KRIPAD(MCTREI,MLENT1) 393 : 394 : C write(6,*)' IK1=',ik1 395 : IF(IK1.LT.4)THEN 396 : DO 33 K=1,NBEL 397 : 398 : NK=NK+1 399 : KC=1+(1-IK1)*(NK-1) 400 : C write(6,*)' K=',K,' NK=',nk,np,kc 401 : 402 : DO 32 J=1,NP 403 : K1=MLENT1.LECT(MELEMA.NUM(J,K)) 404 : II=J 405 : DO 321 I=1,NP 406 : U=VPOCHA(K1,I)*BETA0*IZTGG1.VPOCHA(KC,1) 407 : IF(I.EQ.1)U=ABS(VPOCHA(K1,I))*BETA0*IZTGG1.VPOCHA(KC,1) 408 : IF(II.LE.NP)THEN 409 : AM(K,II,J)=U 410 : ELSE 411 : AM(K,II-NP,J)=U 412 : ENDIF 413 : II=II+1 414 : 321 CONTINUE 415 : 32 CONTINUE 416 : 33 CONTINUE 417 : C write(6,*)' FIN BCL 33 ' 418 : ELSE 419 : 420 : DO 34 K=1,NBEL 421 : 422 : NK=NK+1 423 : UA=1.D0 424 : 425 : DO 35 J=1,NP 426 : K1=MLENT1.LECT(MELEMA.NUM(J,K)) 427 : II=J 428 : DO 351 I=1,NP 429 : U=VPOCHA(K1,I)*BETA0*UA 430 : IF(I.EQ.1)U=ABS(VPOCHA(K1,I))*BETA0*UA 431 : IF(II.LE.NP)THEN 432 : AM(K,II,J)=U 433 : ELSE 434 : AM(K,II-NP,J)=U 435 : ENDIF 436 : II=II+1 437 : 351 CONTINUE 438 : 35 CONTINUE 439 : 34 CONTINUE 440 : 441 : 442 : ENDIF 443 : SEGSUP MLENT1 444 : ENDIF 445 : C write(6,*)' Apres 34 ' 446 : 447 : C CAS Stabilisation via MACRO CENTREP1 448 : DEUPI=1.D0 449 : IF(IAXI.NE.0)DEUPI=2.D0*XPI 450 : IF(IKAS.EQ.3)DEUPI=-DEUPI 451 : 452 : IF(MACRO1.NE.0.AND.IKAS.NE.2.AND.KPRE.EQ.4.AND.IDIM.EQ.2)THEN 453 : I2=2 454 : NBME=1 455 : NBSOUS=1 456 : SEGINI IMATRI 457 : IRIGEL(4,i2)=IMATRI 458 : KSPGP=MELEMC 459 : KSPGD=MELEMC 460 : IRIGEL(1,i2)=MELEMP 461 : IRIGEL(2,i2)=MELEMP 462 : IRIGEL(7,i2)=0 463 : 464 : SEGACT MELEMP 465 : NBSOUS=MELEMP.LISOUS(/1) 466 : C write(6,*)' MELEMP2 NBSOUS=',NBSOUS,' MELEMP=',MELEMP 467 : IF(NBSOUS.NE.0)THEN 467 : CALL ERREUR(5) 468 : ENDIF 469 : 470 : SEGACT MELEMI 471 : 472 : NBEL=MELEMP.NUM(/2) 473 : NP =MELEMP.NUM(/1) 474 : MP =NP 475 : NOM0=NOMS(MELEMP.ITYPEL)//' ' 476 : 477 : SEGINI IZAFM 478 : LIZAFM(1,1)=IZAFM 479 : LISPRI(1)=NOMD 480 : LISDUA(1)=NOMD 481 : 481 : CALL KALPBG(NOM0,'FONFORM ',IZFFM) 482 : 483 : SEGACT IZFFM*MOD 484 : IZHR=KZHR(1) 485 : SEGACT IZHR*MOD 486 : NES=GR(/1) 487 : NPG=GR(/3) 488 : 489 : NK=0 490 : DO 430 LM=1,MAX(1,MELEMI.LISOUS(/1)) 491 : IPT1=MELEMI 492 : IF(MELEMI.LISOUS(/1).NE.0)IPT1=MELEMI.LISOUS(LM) 493 : SEGACT IPT1 494 : ITYP=IPT1.ITYPEL 495 : BETA=0.D0 496 : IF(NOMS(ITYP).EQ.'TRI6')BETA=BETA0 497 : C write(6,*)' NOMS(ITYP)=',NOMS(ITYP),BETA 498 : NBEL=IPT1.NUM(/2) 499 : 500 : C Cas coefficient FLOTTANT ou CENTRE 501 : IF(IK1.LT.4)THEN 502 : DO 43 KK=1,NBEL 503 : NK=NK+1 504 : DO 42 I=1,NP 505 : J=MELEMP.NUM(I,NK) 506 : DO 42 N=1,IDIM 507 : XYZ(N,I)=XCOOR((J-1)*(IDIM+1) +N) 508 : 42 CONTINUE 509 : 509 : CALL CALJBR 510 : &(FN,GR,PG,XYZ,HR,PGSQ,RPG,NES,IDIM,NP,NPG,IAXI,AIRE,AJ,SGN) 511 : KC=1+(1-IK1)*(NK-1) 512 : 513 : DO 44 I=1,NP 514 : DO 44 J=1,NP 515 : U=0.D0 516 : DO 433 L=1,NPG 517 : UH=0.D0 518 : DO 4333 KH=1,IDIM 519 : UH=UH+HR(KH,J,L)*HR(KH,I,L) 520 : 4333 CONTINUE 521 : U=U-UH*PGSQ(L)*DEUPI*BETA*IZTGG1.VPOCHA(KC,1) 522 : 433 CONTINUE 523 : AM(NK,I,J)=U 524 : 44 CONTINUE 525 : 526 : 43 CONTINUE 527 : 528 : ELSEIF(IK1.EQ.4.OR.IK1.EQ.5)THEN 529 : C Cas coefficient SOMMET 530 : 531 : DO 45 KK=1,NBEL 532 : NK=NK+1 533 : DO 46 I=1,NP 534 : J=MELEMP.NUM(I,NK) 535 : DO 46 N=1,IDIM 536 : XYZ(N,I)=XCOOR((J-1)*(IDIM+1) +N) 537 : 46 CONTINUE 538 : 538 : CALL CALJBR 539 : &(FN,GR,PG,XYZ,HR,PGSQ,RPG,NES,IDIM,NP,NPG,IAXI,AIRE,AJ,SGN) 540 : 541 : DO 47 I=1,NP 542 : DO 47 J=1,NP 543 : U=0.D0 544 : DO 453 L=1,NPG 545 : UH=0.D0 546 : DO 4533 KH=1,IDIM 547 : UH=UH+HR(KH,J,L)*HR(KH,I,L) 548 : 4533 CONTINUE 549 : 550 : UA=0.D0 551 : DO 454 JJ=1,NP 552 : J1=LECT(IPT1.NUM(JJ,NK)) 553 : DO 454 KH=1,IDIM 554 : KG=(IK1-4)*(KH-1)+1 555 : UA=UA+FN(JJ,L)*IZTGG1.VPOCHA(J1,KG) 556 : 454 CONTINUE 557 : UA=UA/KG 558 : 559 : U=U-UH*PGSQ(L)*DEUPI*BETA*UA 560 : 453 CONTINUE 561 : AM(NK,I,J)=U 562 : 47 CONTINUE 563 : 564 : 45 CONTINUE 565 : ENDIF 566 : 567 : SEGDES IPT1 568 : 430 CONTINUE 569 : 570 : C write(6,*)'C Fin Stabilisation de toutes sortes' 571 : ENDIF 572 : C Fin Stabilisation de toutes sortes 573 : 574 : NBME=IDIM 574 : CALL KRIPAD(MELEMI,MLENTI1) 575 : SEGACT MELEMI 576 : NBSOUS=MELEMI.LISOUS(/1) 577 : IF(NBSOUS.EQ.0)NBSOUS=1 578 : SEGINI IMATRI 579 : 580 : IF(IKAS.EQ.2)THEN 581 : KSPGD=MELEMS 582 : KSPGP=MELEMC 583 : IRIGEL(2,1)=MELEMI 584 : IRIGEL(1,1)=MELEMP 585 : ELSE 586 : KSPGP=MELEMS 587 : KSPGD=MELEMC 588 : IRIGEL(1,1)=MELEMI 589 : IRIGEL(2,1)=MELEMP 590 : ENDIF 591 : SEGACT MELEMP 592 : 593 : C write(6,*)' ds kmac melemp=',IRIGEL(1,1) 594 : C write(6,*)' ds kmac melemd=',IRIGEL(2,1) 595 : 596 : IRIGEL(4,1)=IMATRI 597 : IF(NOMUL.EQ.0)THEN 598 : IRIGEL(7,1)=2 599 : ELSE 600 : IF(IKAS.EQ.1)IRIGEL(7,1)=-3 601 : IF(IKAS.EQ.2)IRIGEL(7,1)=3 602 : IF(IKAS.EQ.3)IRIGEL(7,1)=4 603 : ENDIF 604 : 605 : SEGACT MELEMP 606 : 607 : K0=0 608 : DO 11 L=1,MAX(1,MELEMI.LISOUS(/1)) 609 : IPT1=MELEMI 610 : IPT2=MELEMP 611 : IF(MELEMI.LISOUS(/1).NE.0)IPT1=MELEMI.LISOUS(L) 612 : IF(MELEMP.LISOUS(/1).NE.0)IPT2=MELEMP.LISOUS(L) 613 : SEGACT IPT1,IPT2 614 : NBEL=IPT1.NUM(/2) 615 : 616 : IF(IKAS.EQ.2)THEN 617 : MP=IPT1.NUM(/1) 618 : C avt NP=MELEMP.NUM(/1) 619 : NP=IPT2.NUM(/1) 620 : ELSE 621 : NP=IPT1.NUM(/1) 622 : C avt MP=MELEMP.NUM(/1) 623 : MP=IPT2.NUM(/1) 624 : ENDIF 625 : 626 : DO 12 I=1,NBME 627 : SEGINI IZAFM 628 : LIZAFM(L,I)=IZAFM 629 : IF(IKAS.EQ.2)THEN 630 : WRITE(NOM,FMT='(I1,A3)')I,NOMD(1:3) 631 : LISDUA(I)=NOM//' ' 632 : LISPRI(I)=NOMP 633 : ELSE 634 : WRITE(NOM,FMT='(I1,A3)')I,NOMP(1:3) 635 : LISPRI(I)=NOM//' ' 636 : LISDUA(I)=NOMD 637 : ENDIF 638 : 12 CONTINUE 639 : IPM1=LIZAFM(L,1) 640 : IPM2=LIZAFM(L,2) 641 : IPM3=LIZAFM(L,2) 642 : IF(IDIM.EQ.3)IPM3=LIZAFM(L,3) 643 : 644 : IF(IK1.LT.4)THEN 645 : C write(6,*)' Appel KPRUSS' 645 : CALL KPRUSS(IPT1,IPT2, 646 : &IPM1,IPM2,IPM3,IAXI,IKAS,INEFMD,KPRE,IZTGG1,IK1,K0,IKOMP) 647 : C write(6,*)'Retour KPRUSS' 648 : ELSE 649 : C write(6,*)' Appel KPRJSS' 649 : CALL KPRJSS(IPT1,IPT2,IPM1,IPM2,IPM3, 650 : & IAXI,IKAS,INEFMD,KPRE,IZTGG1,MLENTI,IK1,IKOMP) 651 : ENDIF 652 : 653 : K0=K0+NBEL 654 : SEGDES IPM1,IPM2,IPM3 655 : SEGDES IPT1 656 : 11 CONTINUE 657 : SEGDES MELEMI 658 : IF(IK1.EQ.4)SEGSUP MLENTI 659 : 660 : 661 : C%%%%%%%%%%%%%%%%% CAS DES PRESSIONS CONTINUES %%%%%%%%%%%%%%%%%%%%%%% 662 : 663 : IF(KPRE.EQ.5)THEN 663 : CALL LEKTAB(MTABZ,'MMAIL',MELEME) 663 : CALL LEKTAB(MTABZ,'MSOMMET',MELEMS) 664 : 665 : SEGACT MELEME 666 : NBSOU1=LISOUS(/1) 667 : IF(NBSOU1.EQ.0)NBSOU1=1 668 : 669 : IRIGEL(1,2)=MELEME 670 : IRIGEL(2,2)=MELEME 671 : IRIGEL(7,2)=0 672 : 673 : c%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 674 : 675 : IHV=0 676 : AIMPL=1.D0 677 : C --- CALCUL DU NOMBRE DE PAQUETS DE LRV ELEMENTS --- 677 : CALL KPLRVM(MELEME,LRV,NBSOUS) 678 : NBME=1 679 : 680 : SEGINI IMATRI 681 : IRIGEL(4,2)=IMATRI 682 : KSPGP=MELEMS 683 : KSPGD=MELEMS 684 : LISPRI(1)=NOMD//' ' 685 : LISDUA(1)=NOMD//' ' 686 : 687 : LL=0 688 : NUTOEL=0 689 : SEGACT MELEME 690 : 691 : DO 101 L=1,MAX(1,LISOUS(/1)) 692 : IPT1=MELEME 693 : IF(LISOUS(/1).NE.0)IPT1=LISOUS(L) 694 : SEGACT IPT1 695 : NOM0=NOMS(IPT1.ITYPEL)//' ' 695 : CALL KALPBG(NOM0,'FONFORM ',IZFFM) 696 : SEGACT IZFFM*MOD 697 : IZHR=KZHR(1) 698 : SEGACT IZHR*MOD 699 : NES=GR(/1) 700 : NPG=GR(/3) 701 : 702 : NP = IPT1.NUM(/1) 703 : MP = NP 704 : NBELEM=IPT1.NUM(/2) 705 : 706 : NNN=MOD(NBELEM,LRV) 707 : IF(NNN.EQ.0) NPACK=NBELEM/LRV 708 : IF(NNN.NE.0) NPACK=1+(NBELEM-NNN)/LRV 709 : 710 : DO 701 KPACK=1,NPACK 711 : 712 : C --- POUR CHAQUE PAQUET DE LRV ELEMENTS ou moins 713 : KDEB=1+(KPACK-1)*LRV 714 : KFIN=MIN(NBELEM,KDEB+LRV-1) 715 : 716 : NBEL=KFIN-KDEB+1 717 : LL=LL+1 718 : 719 : SEGINI IPM1,IPS1 720 : C write(6,*)' IPM1,LL=',IPM1,LL 721 : LIZAFM(LL,1)=IPM1 722 : 723 : IPM2=IPM1 724 : IPM3=IPM1 725 : IPS2=IPS1 726 : IPS3=IPS1 727 : 728 : KITT=2 729 : KJTT=IK1 730 : VISCO=IZTGG1 730 : CALL INITD(VISC,3, 1.D0) 731 : 731 : CALL XLAPL(FN,GR,PG,XYZ,HR,PGSQ,RPG,NES,IDIM,NP,NPG,IAXI, 732 : & VISC,VISC,VISC,KITT,KJTT,IK1, 733 : & IPT1.NUM(1,KDEB),NBEL,NUTOEL,XCOOR,AIMPL,IKOMP, 734 : & IPM1.AM,IPM2.AM,IPM3.AM, 735 : & IPS1.AM,IPS2.AM,IPS3.AM, 736 : & NINKO,IHV,IARG,VISC) 737 : 738 : NUTOEL=NUTOEL+NBEL 739 : SEGDES IPM1 740 : 741 : 701 CONTINUE 742 : SEGDES IPT1 743 : 101 CONTINUE 744 : 745 : 746 : 747 : 748 : 749 : 750 : c%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 751 : 752 : ENDIF 753 : 753 : CALL ECROBJ('MATRIK',MATRIK) 754 : 755 : NAT=2 756 : NSOUPO=0 757 : SEGINI MCHPOI 758 : JATTRI(1)=2 758 : CALL ECROBJ('CHPOINT',MCHPOI) 759 : 760 : C write(6,*)' Fin operateur KMIC' 761 : SEGDES IMATRI,MATRIK 762 : RETURN 763 : 1001 FORMAT(20(1X,I5)) 764 : 1002 FORMAT(10(1X,1PE11.4)) 765 : END 766 : 767 : 768 : 769 : 770 : 771 : 772 : 773 : 774 : 775 : 776 : 777 : 778 :
© Cast3M 2003 - All rights reserved.
Disclaimer