1 : C YLAPL SOURCE MAGN 10/05/18 21:17:07 6675 2 : SUBROUTINE YLAPL 3 : IMPLICIT INTEGER(I-N) 4 : IMPLICIT REAL*8 (A-H,O-Z) 5 : C*********************************************************************** 6 : C 7 : C CET OPERATEUR DISCRETISE LE LAPLACIEN 8 : C EN 2D SUR LES ELEMENTS QUA4 ET TRI3 PLAN OU AXI 9 : C EN 3D SUR LES ELEMENTS CUB8 ET PRI6 10 : C EN 0D SUR LES ELEMENTS POI1 (discrétisation en 1D de l'équation 11 : C de conduction de la chaleur : 12 : C dérivée temporelle + laplacien) 13 : C 14 : C CET OPERATEUR EST "SOUS-INTEGRES" 15 : C 16 : C SYNTAXE : 17 : C --------- 18 : C 19 : C LAPL(ALF) INCO TN : 20 : C 21 : C COEFFICIENTS : 22 : C -------------- 23 : C 24 : C 25 : C ALF (SCAL DOMA) DIFFUSIVITE THERMIQUE 26 : C (SCAL ELEM) 27 : C 28 : C INCONNUES : 29 : C ----------- 30 : C 31 : C TN CHAMP DE TEMPERATURE 32 : C 33 : C************************************************************************ 34 : 35 : -INC CCVQUA4 36 : -INC CCOPTIO 37 : -INC CCGEOME 38 : -INC SMCOORD 39 : -INC SMLENTI 40 : POINTEUR IPADI.MLENTI,IPADS.MLENTI 41 : -INC SMELEME 42 : POINTEUR MELEMI.MELEME,IGEOM0.MELEME,MELEMS.MELEME 43 : -INC SMCHPOI 44 : POINTEUR MCHPIN.MCHPOI 45 : POINTEUR IZG1.MCHPOI, IZG2.MCHPOI 46 : POINTEUR MZLAM.MPOVAL 47 : POINTEUR IZGG1.MPOVAL,IZGG2.MPOVAL 48 : POINTEUR IZTU1.MPOVAL,IZTU2.MPOVAL,IZTU3.MPOVAL,IZTU4.MPOVAL 49 : POINTEUR IZTG5.MPOVAL 50 : POINTEUR IZVOL.MPOVAL,IZTCO.MPOVAL,IZDIAE.MPOVAL,IZTDI.MPOVAL 51 : 52 : -INC SIZFFB 53 : -INC SMMATRIK 54 : POINTEUR IPM.IZAFM 55 : 56 : C SEGMENT IMATRS 57 : C INTEGER LIZAFS(NBSOUS,NBME) 58 : C ENDSEGMENT 59 : POINTEUR IPMS.IZAFM,IPS1.IZAFM,IPS2.IZAFM,IPS3.IZAFM 60 : 61 : -INC SMLMOTS 62 : POINTEUR LINCO.MLMOTS 63 : CHARACTER*8 NOMZ,NOMI,NOMA,TYPE,TYPC,NOM,NOM0,CHAI,MTYP 64 : CHARACTER*4 NOM4 65 : REAL*8 XVAL(9) 66 : LOGICAL LOGI 67 : PARAMETER (NTB=2) 68 : CHARACTER*8 LTAB(NTB) 69 : DIMENSION KTAB(NTB),IXV(4) 70 : * initialisons kform a 0 PV 71 : kform=0 72 : C***************************************************************************** 73 : CLAPL 74 : C write(6,*)' DEBUT YLAPL' 75 : C 76 : C Trois traitements différents suivant la discrétisation 2D/3D EF, VF, ou 0D 77 : C (respectivement, table d'entrée de soustype KIZX 78 : C ou de soustype OPER_0D) 79 : C 80 : C 81 : C**** EN VF, LAPN est un operatéur normal, 82 : C 83 : C JACO RESI DELTAT = 'LAPN' 'VF' ... 84 : C 84 : CALL LIRCHA(NOM4,0,IRET) 85 : IF(IRET .NE. 0)THEN 86 : IF(NOM4 .EQ. 'VF ')THEN 86 : CALL YLAPL1 87 : GOTO 9999 88 : ELSE 89 : C 90 : C********** Je m'excuse et je le remets dans la pile 91 : C 91 : CALL REFUS 92 : ENDIF 93 : ENDIF 94 : C 95 : C Nouvelle directive EQUA de EQEX 96 : C 97 : MTYP=' ' 97 : CALL QUETYP(MTYP,0,IRET) 98 : IF(IRET.EQ.0)THEN 99 : C% On attend un des objets : %m1:8 %m9:16 %m17:24 %m25:32 %m33:40 100 : MOTERR( 1: 8) = 'CHAI ' 101 : MOTERR( 9:16) = 'MMODEL ' 102 : MOTERR(17:24) = 'TABLE ' 102 : CALL ERREUR(472) 103 : RETURN 104 : ENDIF 105 : 106 : IF(MTYP.EQ.'MMODEL')THEN 106 : CALL YTCLSF(' L ','LAPN ') 107 : RETURN 108 : ELSEIF(MTYP.EQ.'MOT ')THEN 108 : CALL LIRCHA(CHAI,1,IRET) 108 : CALL YTCLSF(CHAI,'LAPN ') 109 : RETURN 110 : ENDIF 111 : C Fin Nouvelle directive EQUA de EQEX 112 : 113 : LTAB(1) = 'KIZX ' 114 : LTAB(2) = 'OPER_0D ' 115 : KTAB(1) = 0 116 : KTAB(2) = 0 116 : CALL LITABS(LTAB,KTAB,NTB,0,IRET) 117 : IF(IRET.EQ.0)THEN 118 : WRITE(6,*)' Opérateur LAPN :' 119 : WRITE(6,*)' On attend un ensemble de table soustypes' 120 : RETURN 121 : ENDIF 122 : C 123 : C Bifurcation en cas de discrétisation 0D 124 : C 125 : IF (KTAB(1).NE.0) THEN 126 : MTABX = KTAB(1) 127 : ELSEIF (KTAB(2).NE.0) THEN 128 : IPTAB1 = KTAB(2) 128 : CALL TOWA (IPTAB1) 129 : RETURN 130 : ELSE 131 : WRITE(6,*)' Opérateur LAPN :' 132 : WRITE(6,*)' On attend une table de soustype KIZX ou OPER_0D' 133 : RETURN 134 : ENDIF 135 : 136 : C....................................................................... 137 : C 138 : C- Récupération de la table EQEX (pointeur MTAB1) 139 : C 139 : CALL LEKTAB(MTABX,'EQEX',MTAB1) 140 : IF(MTAB1.EQ.0)THEN 141 : C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24 142 : MOTERR( 1: 8) = ' EQEX ' 143 : MOTERR( 9:16) = ' EQEX ' 144 : MOTERR(17:24) = ' KIZX ' 144 : CALL ERREUR(786) 145 : RETURN 146 : ENDIF 146 : CALL ACME(MTAB1,'NAVISTOK',NASTOK) 147 : IF(NASTOK.EQ.0)THEN 147 : CALL ZLAPL(MTABX,MTAB1) 148 : RETURN 149 : ENDIF 150 : C 151 : C- Récupération de la table INCO (pointeur KINC) 152 : C 152 : CALL LEKTAB(MTAB1,'INCO',KINC) 153 : IF(KINC.EQ.0)THEN 154 : C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24 155 : MOTERR( 1: 8) = ' INCO ' 156 : MOTERR( 9:16) = ' INCO ' 157 : MOTERR(17:24) = ' EQEX ' 157 : CALL ERREUR(786) 158 : RETURN 159 : ENDIF 160 : C....................................................................... 161 : 161 : CALL ACMM(MTABX,'NOMZONE',NOMZ) 162 : 162 : CALL LIROBJ('MMODEL',MMDZ,0,IRET) 163 : IF(IRET.EQ.0)THEN 164 : TYPE=' ' 164 : CALL ACMO(MTABX,'DOMZ',TYPE,MMDZ) 165 : IF(TYPE.NE.'MMODEL')THEN 166 : C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24 167 : MOTERR( 1: 8) = ' DOMZ ' 168 : MOTERR( 9:16) = ' DOMZ ' 169 : MOTERR(17:24) = ' KIZX ' 169 : CALL ERREUR(786) 170 : RETURN 171 : ENDIF 172 : ENDIF 173 : 174 : C***************************************************************************** 175 : C OPTIONS 176 : C KIMPL = 0 -> EXPL 1 -> IMPL 2 -> SEMI 177 : C KFORM = 0 -> SI 1 -> EF 2 -> VF 3 -> EFMC 178 : C IDCEN = 0-> rien 1-> CENTREE 2-> SUPGDC 3-> SUPG 4-> TVISQUEU 5-> CNG 179 : 180 : IAXI=0 181 : IF(IFOMOD.EQ.0)IAXI=2 182 : C 183 : C- Récupération de la table des options KOPT (pointeur KOPTI) 184 : C 184 : CALL LEKTAB(MTABX,'KOPT',KOPTI) 185 : IF (KOPTI.EQ.0) THEN 186 : C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24 187 : MOTERR( 1: 8) = ' KOPT ' 188 : MOTERR( 9:16) = ' KOPT ' 189 : MOTERR(17:24) = ' KIZX ' 189 : CALL ERREUR(786) 190 : RETURN 191 : ENDIF 192 : 193 : C***************************************************************************** 194 : C 195 : C- Récupération de la table DOMAINE associée au domaine local 196 : C 197 : 198 : C E/ MMODEL : Pointeur de la table contenant l'information cherchée 199 : C /S IPOINT : Pointeur sur la table DOMAINE 200 : C /S INEFMD : Type formulation INEFMD=1 LINE,=2 MACRO,=3 QUADRATIQUE 201 : C INEFMD=4 LINB 202 : 202 : CALL LEKMOD(MMDZ,MTABZ,INEFMD) 203 : 203 : CALL LEKTAB(MTABZ,'MAILLAGE',MELEME) 203 : CALL LEKTAB(MTABZ,'SOMMET',MELEMS) 203 : CALL LEKTAB(MTABZ,'CENTRE',MELEMC) 204 : IF (IERR.NE.0) RETURN 205 : 206 : SEGACT MELEME 207 : 208 : C***************************************************************************** 209 : C VERIFICATIONS SUR LES INCONNUES 210 : C 211 : C- Récupération du nombre d'inconnues et du nom de l'inconnue NOMI 212 : C 213 : TYPE = 'LISTMOTS' 213 : CALL ACMO(MTABX,'LISTINCO',TYPE,LINCO) 214 : IF (IERR.NE.0) RETURN 215 : SEGACT LINCO 216 : NBINC = LINCO.MOTS(/2) 217 : IF (NBINC.LE.0.OR.NBINC.GE.3) THEN 218 : C Indice %m1:8 : contient plus de %i1 %m9:16 219 : MOTERR( 1:8) = 'LISTINCO' 220 : INTERR(1) = 2 221 : MOTERR(9:16) = ' MOTS ' 221 : CALL ERREUR(799) 222 : RETURN 223 : ENDIF 224 : NOMI = LINCO.MOTS(1) 225 : NOMA = NOMI 226 : NOM4 = NOMI(1:4) 227 : IF (NBINC.EQ.2) THEN 228 : IF (LINCO.MOTS(1).EQ.LINCO.MOTS(2)) THEN 229 : NINCO = 1 230 : ELSE 231 : IF (KFORM.EQ.0) 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) = ' MOTS ' 235 : CALL ERREUR(799) 236 : RETURN 237 : ELSE 238 : NOMA = LINCO.MOTS(2) 239 : ENDIF 240 : ENDIF 241 : ENDIF 242 : C 243 : C- Récupération de l'inconnue 244 : C 245 : TYPE = ' ' 245 : CALL ACMO(KINC,NOMI,TYPE,MCHPOI) 246 : IF (TYPE.NE.'CHPOINT ') THEN 247 : C Indice %m1:8 : ne contient pas un objet de type %m9:16 248 : MOTERR( 1: 8) = 'INC '//NOMI 249 : MOTERR( 9:16) = 'CHPOINT ' 249 : CALL ERREUR(800) 250 : RETURN 251 : ELSE 251 : CALL LICHT(MCHPOI,IZTU1,TYPC,MELEMI) 252 : MCHPIN=MCHPOI 253 : IF (NBINC.EQ.2) THEN 254 : TYPE = ' ' 254 : CALL ACMO(KINC,NOMA,TYPE,MCHPO2) 255 : IF (TYPE.NE.'CHPOINT ') THEN 256 : C Indice %m1:8 : ne contient pas un objet de type %m9:16 257 : MOTERR( 1: 8) = 'INC '//NOMA 258 : MOTERR( 9:16) = 'CHPOINT ' 258 : CALL ERREUR(800) 259 : RETURN 260 : ELSE 260 : CALL LICHT(MCHPO2,IZTU2,TYPC,IGEOM2) 261 : IF (IGEOM2.NE.MELEMI) THEN 262 : C Indice %m1:8 : L'objet %m9:16 n'a pas le bon support géométrique 263 : MOTERR(1: 8) = 'INC '//NOMA 264 : MOTERR(9:16) = 'CHPOINT ' 264 : CALL ERREUR(788) 265 : RETURN 266 : ENDIF 267 : ENDIF 268 : ENDIF 269 : C***************************************************************************** 270 : C Le domaine de definition est donne par le SPG de la premiere inconnue 271 : C Les inconnues suivantes devront posseder ce meme pointeur 272 : C On verifie que les points de la zone sont tous inclus dans ce SPG 273 : IF (KFORM.NE.2) THEN 273 : CALL KRIPAD(MELEMI,IPADI) 274 : IPADS=IPADI 274 : IF(MELEMI.NE.MELEMS)CALL KRIPAD(MELEMS,IPADS) 275 : 275 : CALL VERPAD(IPADI,MELEME,IRET) 276 : IF (IRET.NE.0) THEN 277 : C Indice %m1:8 : L'objet %m9:16 n'a pas le bon support géométrique 278 : MOTERR(1: 8) = 'INC '//NOMI 279 : MOTERR(9:16) = 'CHPOINT ' 279 : CALL ERREUR(788) 280 : RETURN 281 : ENDIF 282 : ENDIF 283 : C***************************************************************************** 284 : ENDIF 285 : 286 : C***************************************************************************** 287 : C KIMPL = 0 -> EXPL 1 -> IMPL 2 -> SEMI 288 : C KFORM = 0 -> SI 1 -> EF 2 -> VF 289 : 289 : CALL ACME(KOPTI,'KFORM',KFORM) 290 : C************************************************************************* 291 : C Je vous arrete tout de suite ! si KFORM=1 on ne fai rien d'autre ici ! 292 : IF(KFORM.EQ.1)THEN 293 : C CAS FORMULATION EF 294 : 295 : NINKO=IZTU1.VPOCHA(/2) 296 : IHV=0 297 : IF(NINKO.EQ.IDIM)IHV=1 298 : 298 : CALL YCLS('LAPN ',MTABX,IHV,MCHPIN,NOMI,MATRIK,MCHPO1) 299 : 299 : CALL ECROBJ('MATRIK',MATRIK) 299 : CALL ECROBJ('CHPOINT',MCHPO1) 300 : C write(6,*)' FIN YLAPL' 301 : RETURN 302 : ENDIF 303 : C************************************************************************* 304 : 304 : CALL ACME(KOPTI,'IKOMP',IKOMP) 304 : CALL ACME(KOPTI,'KIMPL',KIMPL) 304 : CALL ACME(KOPTI,'KMACO',KMACO) 304 : CALL ACMF(KOPTI,'AIMPL',AIMPL) 305 : 306 : IF(KIMPL.EQ.0)THEN 306 : CALL LEKTAB(MTABZ,'MATESI',MATRIK) 307 : SEGACT MATRIK 308 : IMATRI=IRIGEL(4,1) 309 : SEGACT IMATRI 310 : 310 : CALL LEKTAB(MTABZ,'XXDXDY',MCHPDX) 310 : CALL LEKTAB(MTABZ,'XXVOLUM',MCHPVO) 310 : CALL LEKTAB(MTABZ,'XXDIAME',MCHPDE) 310 : CALL LEKTAB(MTABZ,'XXDIEMIN',MCHPDI) 311 : IF (IERR.NE.0) RETURN 312 : 312 : CALL LICHT(MCHPDX,IZTCO,TYPC,IGEOM) 312 : CALL LICHT(MCHPDI,IZTDI,TYPC,IGEOM) 312 : CALL LICHT(MCHPDE,IZDIAE,TYPC,IGEOM) 312 : CALL LICHT(MCHPVO,IZVOL,TYPC,IGEOM) 313 : 314 : ENDIF 315 : 316 : C***************************************************************************** 317 : C Lecture du coefficient 318 : C Type du coefficient : 319 : C IK1=0 CHPOINT IK1=1 scalaire IK1=2 vecteur 320 : 320 : CALL ACME(MTABX,'IARG',IARG) 321 : IF(IARG.GT.1)THEN 322 : WRITE(6,*)' Operateur LAPN ' 323 : WRITE(6,*)' Nombre d''arguments incorrect : ',IARG 324 : WRITE(6,*)' On attend 1 ' 325 : RETURN 326 : ENDIF 327 : 328 : IXV(1)=MELEMC 329 : IXV(2)=1 330 : IXV(3)=0 331 : IXV(4)=MELEMS 332 : 333 : IRET=0 334 : IF(KFORM.EQ.1)IRET=4 335 : 335 : CALL LEKCOF('Opérateur LAPN :', 336 : & MTABX,KINC,1,IXV,MLAM,MZLAM,NPT1,NC1,IK1,IRET) 337 : IF(IRET.EQ.0)RETURN 338 : 339 : C write(6,*)' KFORM,KIMPL=',KFORM,KIMPL 340 : C************************************************************************* 341 : IF(KFORM.EQ.0)THEN 342 : C CAS FORMULATION EF SI (GRESHO) 343 : 344 : IF(KIMPL.NE.0)THEN 345 : WRITE(6,*)' Opérateur LAPN :' 346 : C Option %m1:8 incompatible avec les données 347 : MOTERR( 1: 8) = ' EFM1 ' 348 : WRITE(6,*)' Options incompatibles : EFM1 et (IMPL ou SEMI) ' 348 : CALL ERREUR(803) 349 : RETURN 350 : ENDIF 351 : 352 : NC=IZTU1.VPOCHA(/2) 353 : NPTI=IZTU1.VPOCHA(/1) 354 : TYPE='SOMMET' 354 : CALL KRCHPT(TYPE,MELEMS,NC,IZG1,NOM4) 355 : 355 : CALL LICHT(IZG1,IZGG1,TYPC,IGEOM) 356 : 357 : NPTS=IZGG1.VPOCHA(/1) 358 : 359 : NCOT=IZTCO.VPOCHA(/1) 360 : 361 : SEGACT MELEME 362 : NBSOUS=LISOUS(/1) 363 : IF(NBSOUS.EQ.0)NBSOUS=1 364 : NUTOEL=0 365 : DT=1.E30 366 : 367 : 368 : DO 1 L=1,NBSOUS 369 : IPT1=MELEME 370 : IF(NBSOUS.NE.1)IPT1=LISOUS(L) 371 : SEGACT IPT1 372 : 373 : IZAFM=LIZAFM(L,1) 374 : IPM1=IZAFM 375 : SEGACT IZAFM 376 : IF(IAXI.NE.0)THEN 377 : IPM1=LIZAFM(L,2) 378 : SEGACT IPM1 379 : ENDIF 380 : 381 : NP =IPT1.NUM(/1) 382 : NBEL=IPT1.NUM(/2) 383 : IES=IDIM 384 : NINKO=IZTU1.VPOCHA(/2) 385 : 385 : CALL YCLPLS(AM,IPM1.AM,IPT1.NUM, 386 : & NBEL,NUTOEL,NPTI,NPTS,NINKO,IES,NP,IAXI,IPADI.LECT,IPADS.LECT, 387 : & MZLAM.VPOCHA,IK1, 388 : & IZTU1.VPOCHA,IZGG1.VPOCHA, 389 : & IZVOL.VPOCHA,IZTCO.VPOCHA,NCOT,IZDIAE.VPOCHA,IZTDI.VPOCHA, 390 : & DT,DTT2,NUEL,DIAEL) 391 : 392 : SEGDES IZAFM,IPT1 393 : IF(IAXI.NE.0)SEGDES IPM1 394 : NUTOEL=NUTOEL+NBEL 395 : 1 CONTINUE 396 : SEGDES MATRIK,IMATRI,MELEME 397 : DTT1=0. 398 : 398 : CALL LEKTAB(MTAB1,'PASDETPS',MTABT) 399 : 400 : IF(MTABT.EQ.0)THEN 400 : CALL CRTABL(MTABT) 400 : CALL ECMM(MTABT,'SOUSTYPE','PASDETPS') 400 : CALL ECMO(MTAB1,'PASDETPS','TABLE ',MTABT) 401 : DTP=1.E30+DT 402 : IPAT=1 403 : DTT1=0. 403 : CALL ECME(MTABT,'NUPASDT',IPAT) 404 : ELSE 404 : CALL ACMF(MTABT,'DELTAT',DTP) 405 : ENDIF 406 : 407 : IF(DT.LT.DTP)THEN 407 : CALL ECMF(MTABT,'DELTAT',DT) 407 : CALL ECMM(MTABT,'OPER','LAPL') 407 : CALL ECMM(MTABT,'ZONE',NOMZ) 407 : CALL ECMF(MTABT,'DTCONV',DTT1) 407 : CALL ECMF(MTABT,'DTDIFU',DTT2) 407 : CALL ECMF(MTABT,'DIAEL',DIAEL) 407 : CALL ECME(MTABT,'NUEL',NUEL) 408 : ENDIF 409 : 410 : SEGDES IZTU1 411 : SEGDES IZGG1 412 : IF(IK1.EQ.0)THEN 413 : SEGDES MZLAM 414 : ENDIF 415 : SEGDES IZVOL,IZTCO,IZDIAE,IZTDI 416 : SEGDES LINCO 417 : 418 : NRIGE=7 419 : NKID =9 420 : NKMT =7 421 : NMATRI=0 422 : SEGINI MATRIK 422 : CALL ECROBJ('MATRIK',MATRIK) 422 : CALL ECROBJ('CHPOINT',IZG1) 423 : 424 : C************************************************************************* 425 : ELSE IF(KFORM.EQ.2)THEN 426 : C CAS FORMULATION VF 426 : CALL LAPLVF(MTABX) 427 : ENDIF 428 : C************************************************************************* 429 : 430 : C write(6,*)' FIN YLAPL' 431 : 9999 CONTINUE 432 : RETURN 433 : 1001 FORMAT(20(1X,I5)) 434 : 1002 FORMAT(10(1X,1PE11.4)) 435 : END 436 : 437 : 438 : 439 : 440 : 441 : 442 : 443 : 444 : 445 : 446 : 447 : 448 : 449 : 450 : 451 : 452 : 453 : 454 : 455 :
© Cast3M 2003 - All rights reserved.
Disclaimer