1 : C KFCE SOURCE CHAT 07/10/29 21:21:49 5940 2 : SUBROUTINE KFCE(IQUAD,MELEMQ) 3 : IMPLICIT INTEGER(I-N) 4 : IMPLICIT REAL*8 (A-H,O-Z) 5 : C************************************************************************* 6 : C 7 : C Objet : calcule les connectivites faces -> elements (SEG3) 8 : C (Ptc 1 Face Pc2) 9 : C les supports geometriques des faces (POI1) 10 : C 11 : C Syntaxe : A , B , C , D , E = 'KFCE' GEO GEO1 <'INCL' TABD> ; 12 : C 13 : C A : Connectivites elements -> faces (ELTFA) MELAF 14 : C B : Support geometrique des centres (CENTRE) MELEM1 15 : C C : Connectivites faces -> sommets (FACEP) MPFD 16 : C D : Connectivites faces -> elements (FACEL) MFD 17 : C E : Support geometrique des faces (FACE) MF1 18 : C 19 : C 20 : C GEO : objet maillage 21 : C GEO1 : support geometrique de l'objet maillage 22 : C 23 : C INCL : directive indiquant qu'on cherche a inclure les 24 : C points crees dans le domaine TABD 25 : C_________________________________________________________________________ 26 : C 27 : C NPF = IKAS(5,I,NTYEL)= le nombre de point constituant une face pour 28 : C le type d'element NTYEL 29 : C 30 : C 31 : C************************************************************************* 32 : PARAMETER (NBTYEL=10) 33 : CHARACTER*8 LTYPEL(NBTYEL),NOM8,TYPE 34 : DIMENSION IKAS(5,6,NBTYEL),NPS(4),LTPL(4),MPS(4) 35 : DIMENSION LNBFAC(NBTYEL),NBTYF(NBTYEL),LNBSO(25) 36 : C 37 : C CORRESPONDANCE 38 : C maillage SEG2 TRI3 QUA4 CUB8 PRI6 PYR5 TET4 TRI6 SEG3 QUA9 39 : C | | | | | | | | | | 40 : C V V V V V V V V V V 41 : C faces SEG2 TRI3 QUA4 PRI6 PYR5 QUA5 TET4 TRI3 SEG2 QUA4 42 : C 43 : C 44 : C ALIAS seg2 tri3 qua4 pri6 pyr5 qua5 tet4 cub8 tri6 seg3 qua8 45 : C ALIAS numero 2 4 8 16 25 9 23 14 6 3 10 46 : 47 : -INC CCOPTIO 48 : -INC SMCOORD 49 : -INC SMTABLE 50 : POINTEUR MTABD.MTABLE,MTBT.MTABLE 51 : -INC SMLENTI 52 : -INC SMELEME 53 : POINTEUR MFD.MELEME,MF1.MELEME,MELEM1.MELEME 54 : POINTEUR MPFD.MELEME,MELAF.MELEME 55 : POINTEUR MELEMC.MELEME,MELEF1.MELEME 56 : POINTEUR MELES1.MELEME,MELEMQ.MELEME 57 : POINTEUR MELEM0.MELEME 58 : -INC CCGEOME 59 : 60 : SEGMENT TRAV 61 : INTEGER ITRAV(NBFA,4) 62 : ENDSEGMENT 63 : 64 : CHARACTER*4 LISMO(1) 65 : PARAMETER (NTB=1) 66 : DIMENSION KTAB(NTB) 67 : CHARACTER*8 LTAB(NTB) 68 : DATA LTYPEL/'SEG2 ','TRI3 ','QUA4 ','CUB8 ', 69 : & 'PRI6 ','PYR5 ','TET4 ', 70 : & 'TRI6 ','SEG3 ','QUA8 '/ 71 : DATA LNBFAC/ 2 , 3 , 4 , 6 , 5 , 5 , 4 , 72 : & 3 , 2 , 4 / 73 : DATA NBTYF/ 2 , 4 , 8 , 16 , 25 , 9 , 23 , 74 : & 4 , 2 , 8 / 75 : C LNBSO pointe sur la ligne de IKAS en fct du type d'element (position 76 : C dans LNBSO) 77 : DATA LNBSO /0 ,1 ,9 ,2 ,0 ,8 ,0 ,3 , 78 : & 0 ,10,0 ,0 ,0 ,4 ,0 ,5 , 79 : & 0 ,0 ,0 ,0 ,0 ,0 ,7 ,0 ,6/ 80 : DATA IKAS/ 81 : & 1,0,0,0,1,2,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 82 : & 1,2,0,0,2,2,3,0,0,2,3,1,0,0,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 83 : & 1,2,0,0,2,2,3,0,0,2,3,4,0,0,2,4,1,0,0,2,0,0,0,0,0,0,0,0,0,0, 84 : & 1,2,3,4,4,5,6,7,8,4,1,2,6,5,4,2,3,7,6,4,3,4,8,7,4,4,1,5,8,4, 85 : & 1,2,3,0,3,4,5,6,0,3,1,2,5,4,4,2,3,6,5,4,3,1,4,6,4,0,0,0,0,0, 86 : & 1,2,3,4,4,1,2,5,0,3,2,3,5,0,3,3,4,5,0,3,4,1,5,0,3,0,0,0,0,0, 87 : & 1,2,4,0,3,1,2,3,0,3,4,2,3,0,3,1,4,3,0,3,0,0,0,0,0,0,0,0,0,0, 88 : & 1,3,0,0,2,3,5,0,0,2,5,1,0,0,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 89 : & 1,0,0,0,1,3,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 90 : & 1,3,0,0,2,3,5,0,0,2,5,7,0,0,2,7,1,0,0,2,0,0,0,0,0,0,0,0,0,0/ 91 : DATA LTPL /1,2,4,8/ 92 : DATA LISMO /'INCL'/ 93 : DATA LTAB /'DOMAINE '/ 94 : 95 : C*** 95 : CALL LIROBJ('MAILLAGE',MELEME,1,IRET) 96 : MELEM0=MELEME 97 : IF(IRET.EQ.0)GO TO 90 97 : CALL LIRREE(XVAL,0,IXV) 98 : XVAL=ABS(XVAL) 99 : 100 : MTABD=0 101 : KINC=0 101 : CALL LIRMOT(LISMO,1,IP,0) 102 : IF(IP.NE.0)THEN 103 : KINC=1 104 : NTO=1 104 : CALL LITABS(LTAB,KTAB,NTB,NTO,IRET) 105 : IF(IRET.EQ.0)RETURN 105 : CALL LIRREE(XVAL,1,IXV) 106 : IF(IRET.EQ.0)RETURN 107 : XVAL=ABS(XVAL) 108 : MTABD=KTAB(1) 109 : TYPE=' ' 109 : CALL ACMO(MTABD,'CENTRE',TYPE,MELEMC) 110 : IF(TYPE.NE.'MAILLAGE')RETURN 111 : TYPE=' ' 111 : CALL ACMO(MTABD,'FACE',TYPE,MELEF1) 112 : IF(TYPE.NE.'MAILLAGE')RETURN 113 : TYPE=' ' 113 : CALL ACMO(MTABD,'SOMMET',TYPE,MELES1) 114 : IF(TYPE.NE.'MAILLAGE')RETURN 115 : ENDIF 116 : 117 : 118 : SEGACT MELEME 119 : NBSOUS=LISOUS(/1) 120 : 121 : IF(NBSOUS.NE.0)THEN 122 : NBNN=0 123 : NBELEM=0 124 : NBREF=0 125 : SEGINI MELAF 126 : ELSE 127 : NBSOUS=1 128 : ENDIF 129 : 130 : NBFA=0 131 : IELIM=0 132 : NBS=NBSOUS 133 : 134 : DO 1 L=1,NBS 135 : IF(NBS.NE.1)THEN 136 : IPT1=LISOUS(L) 137 : SEGACT IPT1 138 : ELSE 139 : IPT1=MELEME 140 : ENDIF 141 : NOM8=NOMS(IPT1.ITYPEL)//' ' 141 : CALL OPTLI(IP,LTYPEL,NOM8,NBTYEL) 142 : C write(6,*)' NOM8,IP=',NOM8,IP 143 : IF(IP.EQ.0)THEN 144 : WRITE(6,*)' On ne sait rien faire pour cet element ',NOM8 145 : GO TO 90 146 : ENDIF 147 : IF((IDIM.EQ.2.AND.IP.EQ.1).OR. 148 : $ (IDIM.EQ.3.AND.IP.LE.3)) IELIM=1 149 : 150 : NBNN=LNBFAC(IP) 151 : NBELEM=IPT1.NUM(/2) 152 : NBFA=NBFA+NBELEM*NBNN 153 : NBSOUS=0 154 : NBREF=0 155 : SEGINI IPT2 156 : IPT2.ITYPEL=NBTYF(IP) 157 : IF(NBS.NE.1) THEN 158 : MELAF.LISOUS(L)=IPT2 159 : CG Pas besoin de désactiver MELAF, MELEME, IPT1, IPT2 160 : CG SEGDES IPT2 161 : CG SEGDES IPT1 162 : ELSE 163 : MELAF=IPT2 164 : ENDIF 165 : 1 CONTINUE 166 : CG SEGDES MELEME 167 : CG SEGDES MELAF 168 : C write(6,*)' Nombre de face total =',NBFA 169 : 170 : IDIM1=IDIM+1 171 : CG Pas besoin d'activer MELAF, MELEME, IPT1, IPT2 172 : CG SEGACT MELEME,MELAF*MOD 173 : DO 50 L=1,NBS 174 : IF(NBS.NE.1)THEN 175 : IPT1=LISOUS(L) 176 : IPT2=MELAF.LISOUS(L) 177 : CG SEGACT IPT1 178 : CG SEGACT IPT2*MOD 179 : ELSE 180 : IPT1=MELEME 181 : IPT2=MELAF 182 : ENDIF 183 : NOM8=NOMS(IPT1.ITYPEL)//' ' 183 : CALL OPTLI(IP,LTYPEL,NOM8,NBTYEL) 184 : IF(IP.GE.8.AND.IP.LE.10)THEN 185 : IQUAD=1 186 : ELSE 187 : IQUAD=0 188 : ENDIF 189 : 190 : NBFAC=LNBFAC(IP) 191 : NEL=IPT1.NUM(/2) 192 : NP=IPT1.NUM(/1) 193 : 194 : C ici on reserve la place pour les faces (ne concerne pas les TRI6...) 195 : NELN=0 196 : IF(IQUAD.EQ.0)THEN 197 : NELN=NEL*NBFAC 198 : C write(6,*)' NBFAC=',NBFAC,' IP=',IP,' NP=',NP,' NF=',NF 199 : NBPTI=XCOOR(/1)/IDIM1 200 : KF=NBPTI 201 : NBPTS=NBPTI+NELN 202 : SEGADJ MCOORD 203 : ENDIF 204 : 205 : DO 52 I=1,NBFAC 206 : NPF=IKAS(5,I,IP) 207 : I1=IKAS(1,I,IP) 208 : I2=IKAS(2,I,IP) 209 : I3=IKAS(3,I,IP) 210 : I4=IKAS(4,I,IP) 211 : 212 : GO TO (631,632,633,634),NPF 213 : 214 : 631 CONTINUE 215 : 216 : IF(IDIM.EQ.2)THEN 217 : C write(6,*)' Cas Segment la face est un poi1 ' 218 : 219 : IF(IQUAD.EQ.0)THEN 220 : DO 851 K=1,NEL 221 : KF1=KF 222 : KF=KF+1 223 : IPT2.NUM(I,K)=KF 224 : N1=IPT1.NUM(I1,K) 225 : N11=N1-1 226 : MD1=N11*IDIM1+1 227 : MD2=N11*IDIM1+2 228 : MG1=KF1*IDIM1+1 229 : MG2=KF1*IDIM1+2 230 : XCOOR(MG1)=XCOOR(MD1) 231 : XCOOR(MG2)=XCOOR(MD2) 232 : XCOOR(KF*IDIM1)=XCOOR(N1*IDIM1) 233 : 851 CONTINUE 234 : ELSE 235 : DO 751 K=1,NEL 236 : IPT2.NUM(I,K)=IPT1.NUM(I1,K) 237 : 751 CONTINUE 238 : ENDIF 239 : 240 : ELSE 241 : 242 : DO 951 K=1,NEL 243 : KF1=KF 244 : KF=KF+1 245 : IPT2.NUM(I,K)=KF 246 : N1=IPT1.NUM(I1,K) 247 : N11=N1-1 248 : MD1=N11*IDIM1+1 249 : MD2=N11*IDIM1+2 250 : MD3=N11*IDIM1+3 251 : MG1=KF1*IDIM1+1 252 : MG2=KF1*IDIM1+2 253 : MG3=KF1*IDIM1+3 254 : XCOOR(MG1)=XCOOR(MD1) 255 : XCOOR(MG2)=XCOOR(MD2) 256 : XCOOR(MG3)=XCOOR(MD3) 257 : XCOOR(KF*IDIM1)=XCOOR(N1*IDIM1) 258 : 951 CONTINUE 259 : ENDIF 260 : 261 : GO TO 640 262 : 263 : 632 CONTINUE 264 : 265 : IF(IDIM.EQ.2)THEN 266 : C write(6,*)' Cas Triangle ou quadrangle la face est un seg2 ' 267 : 268 : IF(IQUAD.EQ.0)THEN 269 : DO 852 K=1,NEL 270 : KF1=KF 271 : KF=KF+1 272 : IPT2.NUM(I,K)=KF 273 : N1=IPT1.NUM(I1,K) 274 : N2=IPT1.NUM(I2,K) 275 : MG1=KF1*IDIM1+1 276 : MG2=KF1*IDIM1+2 277 : XCOOR(MG1)= 278 : & (XCOOR((N1-1)*IDIM1+1)+XCOOR((N2-1)*IDIM1+1)) 279 : $ /2.D0 280 : XCOOR(MG2)= 281 : & (XCOOR((N1-1)*IDIM1+2)+XCOOR((N2-1)*IDIM1+2)) 282 : $ /2.D0 283 : XCOOR(KF*IDIM1)=(XCOOR(N1*IDIM1)+XCOOR(N2*IDIM1))/2 284 : $ .D0 285 : C write(6,*)' KF=',KF,XCOOR(MG1),XCOOR(MG2) 286 : 852 CONTINUE 287 : ELSE 288 : C write(6,*)' I1,I=',i1,i 289 : DO 752 K=1,NEL 290 : IPT2.NUM(I,K)=IPT1.NUM(I1+1,K) 291 : 752 CONTINUE 292 : ENDIF 293 : 294 : ELSE 295 : 296 : DO 952 K=1,NEL 297 : KF1=KF 298 : KF=KF+1 299 : IPT2.NUM(I,K)=KF 300 : N1=IPT1.NUM(I1,K) 301 : N2=IPT1.NUM(I2,K) 302 : MG1=KF1*IDIM1+1 303 : MG2=KF1*IDIM1+2 304 : MG3=KF1*IDIM1+3 305 : XCOOR(MG1)= 306 : & (XCOOR((N1-1)*IDIM1+1)+XCOOR((N2-1)*IDIM1+1))/2 307 : $ .D0 308 : XCOOR(MG2)= 309 : & (XCOOR((N1-1)*IDIM1+2)+XCOOR((N2-1)*IDIM1+2))/2 310 : $ .D0 311 : XCOOR(MG3)= 312 : & (XCOOR((N1-1)*IDIM1+3)+XCOOR((N2-1)*IDIM1+3))/2 313 : $ .D0 314 : XCOOR(KF*IDIM1)=(XCOOR(N1*IDIM1)+XCOOR(N2*IDIM1))/2.D0 315 : 952 CONTINUE 316 : ENDIF 317 : 318 : GO TO 640 319 : 320 : 633 CONTINUE 321 : DO 853 K=1,NEL 322 : KF1=KF 323 : KF=KF+1 324 : IPT2.NUM(I,K)=KF 325 : N1=IPT1.NUM(I1,K) 326 : N2=IPT1.NUM(I2,K) 327 : N3=IPT1.NUM(I3,K) 328 : MG1=KF1*IDIM1+1 329 : MG2=MG1+1 330 : MG3=MG1+2 331 : XCOOR(MG1)= 332 : & (XCOOR((N1-1)*IDIM1+1)+XCOOR((N2-1)*IDIM1+1) 333 : & +XCOOR((N3-1)*IDIM1+1))/3.D0 334 : XCOOR(MG2)= 335 : & (XCOOR((N1-1)*IDIM1+2)+XCOOR((N2-1)*IDIM1+2) 336 : & +XCOOR((N3-1)*IDIM1+2))/3.D0 337 : XCOOR(MG3)= 338 : & (XCOOR((N1-1)*IDIM1+3)+XCOOR((N2-1)*IDIM1+3) 339 : & +XCOOR((N3-1)*IDIM1+3))/3.D0 340 : XCOOR(KF*IDIM1)=(XCOOR(N1*IDIM1)+XCOOR(N2*IDIM1) 341 : & +XCOOR(N3*IDIM1))/3.D0 342 : 853 CONTINUE 343 : 344 : GO TO 640 345 : 346 : 634 CONTINUE 347 : DO 854 K=1,NEL 348 : KF1=KF 349 : KF=KF+1 350 : IPT2.NUM(I,K)=KF 351 : N1=IPT1.NUM(I1,K) 352 : N2=IPT1.NUM(I2,K) 353 : N3=IPT1.NUM(I3,K) 354 : N4=IPT1.NUM(I4,K) 355 : MG1=KF1*IDIM1+1 356 : MG2=MG1+1 357 : MG3=MG1+2 358 : XCOOR(MG1)= 359 : & (XCOOR((N1-1)*IDIM1+1)+XCOOR((N2-1)*IDIM1+1) 360 : & +XCOOR((N3-1)*IDIM1+1)+XCOOR((N4-1)*IDIM1+1))/4.D0 361 : XCOOR(MG2)= 362 : & (XCOOR((N1-1)*IDIM1+2)+XCOOR((N2-1)*IDIM1+2) 363 : & +XCOOR((N3-1)*IDIM1+2)+XCOOR((N4-1)*IDIM1+2))/4.D0 364 : XCOOR(MG3)= 365 : & (XCOOR((N1-1)*IDIM1+3)+XCOOR((N2-1)*IDIM1+3) 366 : & +XCOOR((N3-1)*IDIM1+3)+XCOOR((N4-1)*IDIM1+3))/4.D0 367 : XCOOR(KF*IDIM1)=(XCOOR(N1*IDIM1)+XCOOR(N2*IDIM1) 368 : & +XCOOR(N3*IDIM1)+XCOOR(N4*IDIM1))/4.D0 369 : 854 CONTINUE 370 : 371 : 640 CONTINUE 372 : 373 : 52 CONTINUE 374 : CG Normalement, il n'y a pas besoin de désactiver MELAF et MELEME 375 : CG vu qu'ils sont réutilisés juste après 376 : CG IF(NBS.NE.1)THEN 377 : CG SEGDES IPT2 378 : CG SEGDES IPT1 379 : CG ENDIF 380 : 50 CONTINUE 381 : CG SEGDES MELAF 382 : CG SEGDES MELEME 383 : 383 : CALL CRTABL(MTBT) 383 : CALL NOMOBJ('TABLE','tabletmp',MTBT) 384 : 384 : CALL ECMO(MTBT,'MELAF','MAILLAGE',MELAF) 385 : IF(IXV.EQ.0)XVAL=XCOOR(KF*IDIM1)*(1.D-4)+1.D-8 385 : CALL ECROBJ('MAILLAGE',MELAF) 386 : WRITE(6,1951) 387 : 1951 FORMAT(1H+,' Creation des points face ') 387 : CALL ECRREE(XVAL) 387 : CALL PRELIM(0) 387 : CALL LIROBJ('MAILLAGE',MELAF,1,IRET) 387 : CALL ECMO(MTBT,'MELAF','MAILLAGE',MELAF) 387 : CALL ECRCHA('POI1') 387 : CALL ECROBJ('MAILLAGE',MELAF) 387 : CALL PRCHAN 387 : CALL LIROBJ('MAILLAGE',MF1,1,IRET) 388 : IF(IRET.EQ.0)RETURN 388 : CALL ECMO(MTBT,'MF1','MAILLAGE',MF1) 388 : CALL KNBEL(MF1,NBFA) 389 : 390 : C Ici On cree les points centre 391 : C write(6,*)' Ici On cree les points centre' 392 : 392 : CALL ECROBJ('MAILLAGE',MELEM0) 392 : CALL CRECTR 392 : CALL LIROBJ('MAILLAGE',MELEM1,1,IRET) 393 : IF(IRET.EQ.0)RETURN 394 : 395 : IF(IQUAD.EQ.1)THEN 396 : 397 : C Cas des éléments quadratiques on cree des TRI7,QUA9... 398 : C write(6,*)' Cas des éléments quadratiques' 399 : 400 : SEGACT MELEME,MELEM1 401 : NBSOUS=LISOUS(/1) 402 : NBSOU2=NBSOUS 403 : IF(NBSOUS.EQ.0)NBSOU2=1 404 : KP=0 405 : DO 1880 L=1,NBSOU2 406 : IPT1=MELEME 407 : IF(NBSOU2.NE.1)IPT1=LISOUS(L) 408 : SEGACT IPT1 409 : NBNN=IPT1.NUM(/1)+1 410 : IF(IPT1.ITYPEL.EQ.3)NBNN=IPT1.NUM(/1) 411 : NBELEM=IPT1.NUM(/2) 412 : NBREF=0 413 : IF(L.EQ.1.AND.NBSOU2.EQ.1)THEN 414 : SEGINI MELEMQ 415 : IPT2=MELEMQ 416 : ELSEIF(L.EQ.1.AND.NBSOU2.NE.1)THEN 417 : NBNN=0 418 : NBELEM=0 419 : NBREF=0 420 : SEGINI MELEMQ 421 : NBNN=IPT1.NUM(/1)+1 422 : IF(IPT1.ITYPEL.EQ.3)NBNN=IPT1.NUM(/1) 423 : NBELEM=IPT1.NUM(/2) 424 : NBSOUS=0 425 : NBREF=0 426 : SEGINI IPT2 427 : MELEMQ.LISOUS(1)=IPT2 428 : NBSOUS=NBSOU2 429 : ELSE 430 : NBSOUS=0 431 : NBREF=0 432 : SEGINI IPT2 433 : MELEMQ.LISOUS(L)=IPT2 434 : NBSOUS=NBSOU2 435 : ENDIF 436 : 437 : IF(IPT1.ITYPEL.EQ.6)THEN 438 : IPT2.ITYPEL=7 439 : DO 1883 K=1,NBELEM 440 : KP=KP+1 441 : DO 1882 I=1,NBNN-1 442 : IPT2.NUM(I,K)=IPT1.NUM(I,K) 443 : 1882 CONTINUE 444 : IPT2.NUM(NBNN,K)=MELEM1.NUM(1,KP) 445 : 1883 CONTINUE 446 : ELSEIF(IPT1.ITYPEL.EQ.3)THEN 447 : IPT2.ITYPEL=3 448 : DO 1885 K=1,NBELEM 449 : DO 18851 I=1,NBNN 450 : IPT2.NUM(I,K)=IPT1.NUM(I,K) 451 : 18851 CONTINUE 452 : 1885 CONTINUE 453 : ELSEIF(IPT1.ITYPEL.EQ.10)THEN 454 : IPT2.ITYPEL=11 455 : DO 1886 K=1,NBELEM 456 : KP=KP+1 457 : DO 1887 I=1,NBNN-1 458 : IPT2.NUM(I,K)=IPT1.NUM(I,K) 459 : 1887 CONTINUE 460 : IPT2.NUM(NBNN,K)=MELEM1.NUM(1,KP) 461 : 1886 CONTINUE 462 : ELSE 463 : WRITE(6,*)' Type d''élément non encore prévu' 464 : ENDIF 465 : 466 : 1880 CONTINUE 466 : CALL ECMO(MTBT,'MELEMQ','MAILLAGE',MELEMQ) 467 : ENDIF 468 : 469 : 470 : C<<<<<<<<<<<<< 470 : CALL ECMO(MTBT,'MELEM1','MAILLAGE',MELEM1) 470 : C CALL ECMO(MTBT,'MELEM0','MAILLAGE',MELEM0) 470 : C CALL ECMO(MTBT,'MELAF','MAILLAGE',MELAF) 471 : 472 : IF(IELIM.NE.0)THEN 473 : 474 : WRITE(6,1952) 475 : 1952 FORMAT(1H+,' Elimination entre points centre et points face') 475 : CALL ECROBJ('MAILLAGE',MELEM1) 475 : CALL ECRREE(XVAL) 475 : CALL ECROBJ('MAILLAGE',MELAF) 475 : CALL PRELIM(0) 475 : CALL LIROBJ('MAILLAGE',MELEM1,1,IRET) 476 : IF(IRET.EQ.0)RETURN 477 : 478 : WRITE(6,1954) 479 : 1954 FORMAT(1H+,' Elimination entre points face et points sommet') 479 : CALL ECROBJ('MAILLAGE',MELEM0) 479 : CALL ECRREE(XVAL) 479 : CALL ECROBJ('MAILLAGE',MELAF) 479 : CALL PRELIM(0) 479 : CALL LIROBJ('MAILLAGE',MELEM0,1,IRET) 480 : IF(IRET.EQ.0)RETURN 481 : 482 : ENDIF 483 : 484 : C Ici On cree les connectivites FACEL (SEG3) 485 : C -> In KRIPAD : SEGINI MLENTI 486 : C -> In KRIPAD : SEGACT MF1 486 : CALL KRIPAD(MF1,MLENTI) 487 : SEGDES MF1 488 : SEGACT MELEM1 489 : 490 : NBNN=3 491 : NBELEM=NBFA 492 : NBSOUS=0 493 : NBREF=0 494 : SEGINI MFD 495 : MFD.ITYPEL=3 496 : SEGACT MELAF 497 : NBSOUS=MELAF.LISOUS(/1) 498 : IF(NBSOUS.EQ.0)NBSOUS=1 499 : NK=0 500 : DO 784 NS=1,NBSOUS 501 : IF(NBSOUS.NE.1)THEN 502 : MELEME=MELAF.LISOUS(NS) 503 : SEGACT MELEME 504 : ELSE 505 : MELEME=MELAF 506 : ENDIF 507 : NEL=NUM(/2) 508 : NP =NUM(/1) 509 : DO 783 K=1,NEL 510 : NK=NK+1 511 : NNK=MELEM1.NUM(1,NK) 512 : DO 782 I=1,NP 513 : NUF=NUM(I,K) 514 : NUF1=LECT(NUF) 515 : MFD.NUM(1,NUF1)=NNK 516 : MFD.NUM(2,NUF1)=NUF 517 : 782 CONTINUE 518 : 783 CONTINUE 519 : CG Pas besoin de désactiver 520 : CG IF(NBSOUS.NE.1)SEGDES MELEME 521 : CG SEGDES MELEME 522 : 784 CONTINUE 523 : 524 : NK=NK+1 525 : DO 785 NS=NBSOUS,1,-1 526 : IF(NBSOUS.NE.1)THEN 527 : MELEME=MELAF.LISOUS(NS) 528 : CG Pas besoin de réactiver 529 : CG SEGACT MELEME 530 : ELSE 531 : MELEME=MELAF 532 : ENDIF 533 : NEL=NUM(/2) 534 : NP =NUM(/1) 535 : DO 786 K=NEL,1,-1 536 : NK=NK-1 537 : NNK=MELEM1.NUM(1,NK) 538 : DO 787 I=1,NP 539 : NUF=NUM(I,K) 540 : NUF1=LECT(NUF) 541 : MFD.NUM(3,NUF1)=NNK 542 : 787 CONTINUE 543 : 786 CONTINUE 544 : CG Pas besoin de désactiver mais après, il s'appelle IPT1 ! 545 : CG IF(NBSOUS.NE.1)SEGDES MELEME 546 : 785 CONTINUE 547 : CG SEGDES MELAF 548 : SEGDES MFD 549 : SEGDES MELEM1 550 : 551 : 552 : C ici on cree les connectivites face --> sommet 553 : 554 : CG Sous-entendu : NBFA=nb_éléments(MF1) 555 : SEGINI TRAV 556 : SEGACT MELEM0 557 : CG SEGACT MELAF 558 : NBSOUS=MELAF.LISOUS(/1) 559 : IF(NBSOUS.EQ.0)NBSOUS=1 560 : DO 338 L=1,NBSOUS 561 : IF(NBSOUS.NE.1)THEN 562 : IPT1=MELAF.LISOUS(L) 563 : IPT2=MELEM0.LISOUS(L) 564 : CG Déja activé plus haut 565 : CG SEGACT IPT1 566 : SEGACT IPT2 567 : ELSE 568 : IPT1=MELAF 569 : IPT2=MELEM0 570 : ENDIF 571 : NP =IPT1.NUM(/1) 572 : NEL=IPT1.NUM(/2) 573 : IP=LNBSO(IPT2.ITYPEL) 574 : DO 339 K=1,NEL 575 : DO 3391 I=1,NP 576 : NUF=IPT1.NUM(I,K) 577 : NUF=LECT(NUF) 578 : ITRAV(NUF,1)=IPT2 579 : ITRAV(NUF,2)=IP 580 : ITRAV(NUF,3)=I 581 : ITRAV(NUF,4)=K 582 : 3391 CONTINUE 583 : 339 CONTINUE 584 : IF(NBSOUS.NE.1)SEGDES IPT1 585 : 338 CONTINUE 586 : SEGDES MELAF 587 : SEGSUP MLENTI 588 : 589 : C write(6,*)' Fin premiere partie ' 589 : CALL INITI(NPS,4,0) 590 : DO 340 L=1,NBFA 591 : IP=ITRAV(L,2) 592 : I =ITRAV(L,3) 593 : NBS=IKAS(5,I,IP) 594 : C write(6,*)' IP,I,nbs=',IP,I,nbs 595 : NPS(NBS)=NPS(NBS)+1 596 : 340 CONTINUE 597 : C write(6,*)' NPS=',NPS 598 : NT=0 599 : DO 341 I=1,4 600 : IF(NPS(I).NE.0)THEN 601 : NT=NT+1 602 : IT=I 603 : ENDIF 604 : 341 CONTINUE 605 : NBREF=0 606 : IF(NT.EQ.1)THEN 607 : NBSOUS=0 608 : NBELEM=NPS(IT) 609 : NBNN=IT 610 : SEGINI MELEME 611 : ITYPEL=LTPL(IT) 612 : NPS(IT)=MELEME 613 : ELSE 614 : NBSOUS=NT 615 : NBELEM=0 616 : NBNN=0 617 : SEGINI MELEME 618 : JJ=0 619 : DO 342 J=1,4 620 : IF(NPS(J).EQ.0)GO TO 342 621 : JJ=JJ+1 622 : NBSOUS=0 623 : NBREF=0 624 : NBNN=J 625 : NBELEM=NPS(J) 626 : SEGINI IPT3 627 : LISOUS(JJ)=IPT3 628 : IPT3.ITYPEL=LTPL(J) 629 : NPS(J)=IPT3 630 : C write(6,*)' ipt3=',ipt3,' j=',j,' jj=',jj 631 : 342 CONTINUE 632 : ENDIF 633 : MPFD=MELEME 634 : C write(6,*)' Fin deuxieme partie MPFD=',MPFD 635 : 635 : CALL INITI(MPS,4,0) 636 : 637 : DO 343 L=1,NBFA 638 : IPT2=ITRAV(L,1) 639 : IP=ITRAV(L,2) 640 : I =ITRAV(L,3) 641 : K =ITRAV(L,4) 642 : NBS=IKAS(5,I,IP) 643 : IPT3=NPS(NBS) 644 : MPS(NBS)=MPS(NBS)+1 645 : KK=MPS(NBS) 646 : DO 344 J=1,NBS 647 : JJ=IKAS(J,I,IP) 648 : C write(6,*)'j,kk=',j,kk,' jj,k',jj,k,' nbs=',nbs 649 : IPT3.NUM(J,KK)=IPT2.NUM(JJ,K) 650 : 344 CONTINUE 651 : 343 CONTINUE 652 : 653 : IF(MPFD.LISOUS(/1).NE.0)THEN 654 : NSO=MPFD.LISOUS(/1) 655 : DO 765 I=1,NSO 656 : IPT3=MPFD.LISOUS(I) 657 : SEGDES IPT3 658 : 765 CONTINUE 659 : ENDIF 660 : SEGDES MPFD 661 : 662 : IF(MELEM0.LISOUS(/1).NE.0)THEN 663 : NSO=MELEM0.LISOUS(/1) 664 : DO 766 I=1,NSO 665 : IPT3=MELEM0.LISOUS(I) 666 : SEGDES IPT3 667 : 766 CONTINUE 668 : SEGDES MELEM0 669 : ENDIF 670 : 671 : SEGSUP TRAV 672 : IF(KINC.EQ.1)THEN 673 : WRITE(6,1953) 674 : 1953 FORMAT(1H+ 675 : $ ,' Elimination avec les points d''un autre domaine') 675 : CALL ECMO(MTBT,'MFD','MAILLAGE',MFD) 675 : CALL ECMO(MTBT,'MPFD','MAILLAGE',MPFD) 675 : CALL ECROBJ('MAILLAGE',MELEF1) 675 : CALL ECROBJ('MAILLAGE',MELEMC) 675 : CALL PRFUSE 675 : CALL ECROBJ('MAILLAGE',MELES1) 675 : CALL PRFUSE 675 : CALL LIROBJ('MAILLAGE',IPT1,1,IRET) 676 : IF(IRET.EQ.0)RETURN 676 : CALL ECMO(MTBT,'IPT1','MAILLAGE',IPT1) 676 : CALL ECROBJ('MAILLAGE',MF1) 676 : CALL ECROBJ('MAILLAGE',MELEM1) 676 : CALL PRFUSE 676 : CALL LIROBJ('MAILLAGE',IPT2,1,IRET) 677 : IF(IRET.EQ.0)RETURN 677 : CALL ECMO(MTBT,'IPT2','MAILLAGE',IPT2) 677 : CALL ECROBJ('MAILLAGE',IPT1) 677 : CALL ECRREE(XVAL) 677 : CALL ECROBJ('MAILLAGE',IPT2) 678 : C write(6,*)' Avant prelim ipt1=',ipt1,' ipt2',ipt2 678 : CALL PRELIM(0) 678 : CALL LIROBJ('MAILLAGE',IPT3,1,IRET) 679 : C write(6,*)' retour PRELIM IPT3=',IPT3 680 : SEGSUP IPT1 681 : SEGSUP IPT2 682 : C SEGSUP IPT3 683 : ENDIF 684 : SEGSUP MTBT 684 : CALL ECROBJ('MAILLAGE',MF1) 684 : CALL ECROBJ('MAILLAGE',MFD) 684 : CALL ECROBJ('MAILLAGE',MPFD) 684 : CALL ECROBJ('MAILLAGE',MELEM1) 684 : CALL ECROBJ('MAILLAGE',MELAF) 685 : 686 : RETURN 687 : 688 : 90 CONTINUE 689 : WRITE(6,*)' Retour anormal de KFCE ' 690 : RETURN 691 : END 692 : 693 : 694 : 695 : 696 : 697 : 698 :
© Cast3M 2003 - All rights reserved.
Disclaimer