1 : C RAFF3D SOURCE CHAT 11/03/16 21:30:04 6902 2 : SUBROUTINE RAFF3D(IPT2,ICPR,KARPOS,KARETE,KMILIE,MELVA2,NACREE, 3 : .KARAF,IPT4,JPLANS,JPLAN3,JPLCOM,JNOEFA,IPT7,JFARAF,KARET2) 4 : IMPLICIT REAL*8(A-H,O-Z) 5 : IMPLICIT INTEGER(I-N) 6 : -INC CCOPTIO 7 : -INC SMCOORD 8 : -INC CCGEOME 9 : -INC SMELEME 10 : -INC SMCHPOI 11 : -INC SMMODEL 12 : -INC SMCHAML 13 : C 14 : C====================================================================== 15 : C Declarations 16 : C====================================================================== 17 : SEGMENT ICPR((XCOOR(/1)/(IDIM+1)),2) 18 : SEGMENT KARETE(NBNDS,NCOL) 19 : SEGMENT KARET2(NBNDS,NCOL) 20 : SEGMENT KMILIE(NBNDS,NCOL) 21 : SEGMENT KARPOS(NBNDS) 22 : SEGMENT JPLANS(JPLA1,JPLA2) 23 : SEGMENT JPLAN3(JPLA1,JPLA2) 24 : SEGMENT JPLCOM(JPLA1) 25 : SEGMENT JNOEFA(JNBFA,5) 26 : SEGMENT JFARAF(JNBFA,LLLL) 27 : SEGMENT NUMNOE(INUMNO) 28 : SEGMENT IWORK2(JNBSOM) 29 : SEGMENT IWORK1(IJKLMN) 30 : 31 : REAL*8 LON68,LON510,LON47,LON29 32 : C 33 : C====================================================================== 34 : C Initialisations 35 : C====================================================================== 36 : SEGACT JPLANS,JPLCOM,JNOEFA,JPLAN3*MOD 37 : SEGACT IPT2,ICPR,KARPOS,KARETE,KMILIE*MOD,MELVA2,JFARAF*MOD 38 : IJKLMN=4 39 : ipt7=0 40 : SEGINI IWORK1 41 : IWORK1(1)=4 42 : IWORK1(2)=8 43 : IWORK1(3)=6 44 : IWORK1(4)=10 45 : C JNBFA=JNOEFA(/1) 46 : C LLLL=2 47 : C IF (IPT2.ITYPEL.EQ.15) LLLL=10 48 : C IF ((IPT2.ITYPEL.EQ.17).OR.(IPT2.ITYPEL.EQ.24)) LLLL=6 49 : C LLLL=10 50 : C SEGINI JFARAF 51 : NCOMPL=LNELM(1,(IPT2.ITYPEL-1)*2+2) 52 : IF (NCOMPL.NE.0) GOTO 999 53 : NBNN=NBNNE(LNELM(2,(IPT2.ITYPEL-1)*2+1)) 54 : INELM=LNELM(1,(IPT2.ITYPEL-1)*2+1) 55 : C 56 : C Creation du squelette du maillage resultat 57 : NBELEM=IPT2.NUM(/2)-KARAF+INELM*KARAF 58 : NBSOUS=0 59 : NBREF=0 60 : SEGINI IPT4 61 : IPT4.ITYPEL=LNELM(2,(IPT2.ITYPEL-1)*2+1) 62 : LPO2=LPOS2(IPT2.ITYPEL) 63 : NBPT0=XCOOR(/1)/(IDIM+1) 64 : NBPT1=XCOOR(/1)/(IDIM+1) 65 : NBPTS=NBPT1+NACREE+KARAF*NBINTE(IPT2.ITYPEL) 66 : SEGADJ MCOORD 67 : INUMNO=NBRAF(IPT2.ITYPEL) 68 : SEGINI NUMNOE 69 : LCOMP=1 70 : C 71 : C====================================================================== 72 : C Phase de raffinement 3D 73 : C====================================================================== 74 : C======================================= 75 : C A) Boucle sur les elements a raffiner ! 76 : C======================================= 77 : DO 6 IARAF=1,MELVA2.VELCHE(/2) 78 : IF (MELVA2.VELCHE(1,IARAF).NE.1) THEN 79 : DO 1 IJKL=1,IPT2.NUM(/1) 80 : IPT4.NUM(IJKL,NBELEM)=IPT2.NUM(IJKL,IARAF) 81 : 1 CONTINUE 82 : NBELEM=NBELEM-1 83 : GOTO 6 84 : ENDIF 85 : JCOMPT=0 86 : JPOS5=LPOS5(IPT2.ITYPEL) 87 : C 88 : C========================================================== 89 : C B) Boucle sur les noeuds a creer pour raffiner l'element ! 90 : C========================================================== 91 : DO 4 I=1,NBRAF(IPT2.ITYPEL) 92 : JPOS1=LPOS1(1,I+LPOS2(IPT2.ITYPEL)-1) 93 : JLONG=LPOS1(2,I+LPOS2(IPT2.ITYPEL)-1) 94 : JLISN=LPOS3(IPT2.ITYPEL)+JCOMPT 95 : LTYPNO=JTYPNO(JPOS5-1+I) 96 : C 97 : C------------------------------- 98 : C ** B.1 / On est sur une arete ! 99 : C------------------------------- 100 : IF (LTYPNO.EQ.0) THEN 101 : NPTA=IPT2.NUM(LISNOE(JLISN),IARAF) 102 : NPTB=IPT2.NUM(LISNOE(JLISN+1),IARAF) 103 : NMIN=MIN(ICPR(NPTA,1),ICPR(NPTB,1)) 104 : NMAX=MAX(ICPR(NPTA,1),ICPR(NPTB,1)) 105 : DO 2 K=1,MAX(1,KARPOS(NMIN)) 106 : IF (KARETE(NMIN,K).EQ.NMAX) NEXIST=K 107 : 2 CONTINUE 108 : IF (KMILIE(NMIN,NEXIST).GT.0) THEN 109 : NUMNOE(I)=KMILIE(NMIN,NEXIST) 110 : JCOMPT=JCOMPT+JLONG 111 : GOTO 4 112 : ELSE 113 : NBPT1=NBPT1+1 114 : NUMNOE(I)=NBPT1 115 : KMILIE(NMIN,NEXIST)=NBPT1 116 : ENDIF 117 : C 118 : C------------------------------ 119 : C ** B.2 / On est sur une face ! 120 : C------------------------------ 121 : C Il faut identifier cette face (numero, sommets...) 122 : ELSEIF ((LTYPNO.GT.0).AND.(LTYPNO.LT.7)) THEN 123 : C => a) On initialise toutes les donnees relatives a cette face 124 : JTYPEL=IPT2.ITYPEL 125 : JLTEL2=LTEL(1,JTYPEL)-1+LTYPNO 126 : JLTEL2=LTEL(2,JTYPEL)-1+LTYPNO 127 : JLDEL1=LDEL(1,JLTEL2) 128 : JTYFAC=IWORK1(JLDEL1) 129 : JLDEL2=LDEL(2,JLTEL2) 130 : JNBSOM=NBSOM(JTYFAC) 131 : JSPOS=NSPOS(JTYFAC) 132 : C => b) On identifie les sommets de la face (n° global) 133 : SEGINI IWORK2 134 : DO 10 IAA=1,JNBSOM 135 : NGLOBA=IPT2.NUM(LFAC(JLDEL2-1+IBSOM(JSPOS-1+IAA)),IARAF) 136 : IWORK2(IAA)=NGLOBA 137 : 10 CONTINUE 138 : C => c) On classe ces sommets par ordre croissant (NPTA < NPTB < NPTC) 139 : NPTA=(XCOOR(/1)/(IDIM+1))+1 140 : NPTB=NPTA+1 141 : NPTC=NPTB+1 142 : DO 11 ICC=1,JNBSOM 143 : IF (IWORK2(ICC).LT.NPTA) THEN 144 : NPTC=NPTB 145 : NPTB=NPTA 146 : NPTA=IWORK2(ICC) 147 : ELSEIF (IWORK2(ICC).LT.NPTB) THEN 148 : NPTC=NPTB 149 : NPTB=IWORK2(ICC) 150 : ELSEIF (IWORK2(ICC).LT.NPTC) THEN 151 : NPTC=IWORK2(ICC) 152 : ENDIF 153 : 11 CONTINUE 154 : C => d) On passe ces sommets en n° locale 155 : NPTA2=ICPR(NPTA,1) 156 : NPTB2=ICPR(NPTB,1) 157 : NPTC2=ICPR(NPTC,1) 158 : IF ((NPTA2.LT.NPTB2).AND.(NPTA2.LT.NPTC2)) THEN 159 : NPTA=NPTA2 160 : NPTB=MIN(NPTB2,NPTC2) 161 : NPTC=MAX(NPTB2,NPTC2) 162 : ENDIF 163 : IF ((NPTB2.LT.NPTA2).AND.(NPTB2.LT.NPTC2)) THEN 164 : NPTA=NPTB2 165 : NPTB=MIN(NPTA2,NPTC2) 166 : NPTC=MAX(NPTA2,NPTC2) 167 : ENDIF 168 : IF ((NPTC2.LT.NPTA2).AND.(NPTC2.LT.NPTB2)) THEN 169 : NPTA=NPTC2 170 : NPTB=MIN(NPTA2,NPTB2) 171 : NPTC=MAX(NPTA2,NPTB2) 172 : ENDIF 173 : C => e) On cherche le numero de la face 174 : NEXIS2=0 175 : DO 12 IEE=1,JPLCOM(NPTA) 176 : MTMP=JPLANS(NPTA,IEE) 177 : JJ1=JNOEFA(MTMP,1) 178 : JJ2=JNOEFA(MTMP,2) 179 : JJ3=JNOEFA(MTMP,3) 180 : IF(JJ1.EQ.NPTA.AND.JJ2.EQ.NPTB.AND.JJ3.EQ.NPTC) THEN 181 : NEXIS2=IEE 182 : ENDIF 183 : 12 CONTINUE 184 : JNUMFA=JPLANS(NPTA,NEXIS2) 185 : C 186 : C--------------------------------- 187 : C ** B.3 / Raffinement de la face ! 188 : C--------------------------------- 189 : KTEST1=JPLAN3(NPTA,NEXIS2) 190 : KTEST2=NBINTE(JTYFAC) 191 : C => a) Si la face est de type QUA4 (un seul noeud a creer, au milieu) 192 : IF (JTYFAC.EQ.8) THEN 193 : IF (KTEST1.LT.KTEST2) THEN 194 : NBPT1=NBPT1+1 195 : NUMNOE(I)=NBPT1 196 : JPLAN3(NPTA,NEXIS2)=JPLAN3(NPTA,NEXIS2)+1 197 : JFARAF(JNUMFA,1)=NBPT1 198 : ELSEIF (KTEST1.EQ.KTEST2) THEN 199 : NUMNOE(I)=JFARAF(JNUMFA,1) 200 : JFARAF(JNUMFA,1)=0 201 : JCOMPT=JCOMPT+JLONG 202 : GOTO 4 203 : ELSE 204 : WRITE(*,*) 'ERREUR' 205 : ENDIF 206 : ENDIF 207 : C => b) Si la face n'est pas de type QUA4 (donc de type TRI6 ou QUA8) 208 : IF (JTYFAC.NE.8) THEN 209 : IF (KTEST1.LT.KTEST2) THEN 210 : NBPT1=NBPT1+1 211 : NUMNOE(I)=NBPT1 212 : JPLAN3(NPTA,NEXIS2)=JPLAN3(NPTA,NEXIS2)+1 213 : NEXIS3=0 214 : XCO2=0.25 215 : DO 13 KBB=1,JLONG 216 : XCO1=XCOEFF(JPOS1-1+KBB) 217 : IF (XCO1.EQ.XCO2) NEXIS3=KBB 218 : 13 CONTINUE 219 : JFARAF(JNUMFA,2*(KTEST1)+1)=NBPT1 220 : IF (NEXIS3.NE.0) THEN 221 : NGLOB=IPT2.NUM(LISNOE(JLISN-1+NEXIS3),IARAF) 222 : JFARAF(JNUMFA,2*(KTEST1)+2)=NGLOB 223 : ELSE 224 : JFARAF(JNUMFA,2*(KTEST1)+2)=0 225 : ENDIF 226 : ELSEIF (KTEST1.EQ.KTEST2) THEN 227 : NEXIS3=0 228 : NEXIS4=0 229 : XCO2=0.25 230 : DO 14 KBB=1,JLONG 231 : XCO1=XCOEFF(JPOS1-1+KBB) 232 : IF (XCO1.EQ.XCO2) NEXIS3=KBB 233 : 14 CONTINUE 234 : IF (NEXIS3.NE.0) THEN 235 : NGLOB=IPT2.NUM(LISNOE(JLISN-1+NEXIS3),IARAF) 236 : DO 15 KAA=2,JFARAF(/2),2 237 : IF (JFARAF(JNUMFA,KAA).EQ.NGLOB) NEXIS4=KAA 238 : 15 CONTINUE 239 : NUMNOE(I)=JFARAF(JNUMFA,NEXIS4-1) 240 : JFARAF(JNUMFA,NEXIS4-1)=0 241 : ELSE 242 : DO 16 KAA=2,JFARAF(/2),2 243 : IF (JFARAF(JNUMFA,KAA).EQ.0) NEXIS4=KAA 244 : 16 CONTINUE 245 : NUMNOE(I)=JFARAF(JNUMFA,NEXIS4-1) 246 : JFARAF(JNUMFA,NEXIS4-1)=0 247 : ENDIF 248 : JCOMPT=JCOMPT+JLONG 249 : GOTO 4 250 : ENDIF 251 : ENDIF 252 : C 253 : C------------------------------------------------------ 254 : C ** B.4 / On est a l'interieur du volume de l'element ! 255 : C------------------------------------------------------ 256 : ELSEIF (LTYPNO.EQ.7) THEN 257 : NBPT1=NBPT1+1 258 : NUMNOE(I)=NBPT1 259 : ENDIF 260 : C 261 : IF (NBPT1.EQ.NBPTS) THEN 262 : NBPTS=NBPTS+200 263 : SEGADJ MCOORD 264 : ENDIF 265 : C 266 : C============================== 267 : C C) Creation du nouveau point ! 268 : C============================== 269 : C On continue ici que lorsque l'on doit creer un nouveau point 270 : XPT=0. 271 : YPT=0. 272 : ZPT=0. 273 : DO 3 J=1,JLONG 274 : NGLOB=IPT2.NUM(LISNOE(JLISN-1+J),IARAF) 275 : XINI=XCOOR((NGLOB-1)*(IDIM+1)+1) 276 : YINI=XCOOR((NGLOB-1)*(IDIM+1)+2) 277 : ZINI=XCOOR((NGLOB-1)*(IDIM+1)+3) 278 : XPT=XPT+XINI*XCOEFF(JPOS1-1+J) 279 : YPT=YPT+YINI*XCOEFF(JPOS1-1+J) 280 : ZPT=ZPT+ZINI*XCOEFF(JPOS1-1+J) 281 : 3 CONTINUE 282 : XCOOR((NBPT1-1)*(IDIM+1)+1)=XPT 283 : XCOOR((NBPT1-1)*(IDIM+1)+2)=YPT 284 : XCOOR((NBPT1-1)*(IDIM+1)+3)=ZPT 285 : JCOMPT=JCOMPT+JLONG 286 : 4 CONTINUE 287 : JPOS4=LPOS4(IPT2.ITYPEL) 288 : C 289 : C Cas des tetraedres 290 : LONDIA=0 291 : JTYPEL=IPT2.ITYPEL 292 : IF (JTYPEL.EQ.23) THEN 293 : LPT6=NUMNOE(6) 294 : LPT5=NUMNOE(5) 295 : LPT8=NUMNOE(8) 296 : LPT10=NUMNOE(10) 297 : XPT6=XCOOR((LPT6-1)*(IDIM+1)+1) 298 : XPT5=XCOOR((LPT5-1)*(IDIM+1)+1) 299 : XPT8=XCOOR((LPT8-1)*(IDIM+1)+1) 300 : XPT10=XCOOR((LPT10-1)*(IDIM+1)+1) 301 : YPT6=XCOOR((LPT6-1)*(IDIM+1)+2) 302 : YPT5=XCOOR((LPT5-1)*(IDIM+1)+2) 303 : YPT8=XCOOR((LPT8-1)*(IDIM+1)+2) 304 : YPT10=XCOOR((LPT10-1)*(IDIM+1)+2) 305 : ZPT6=XCOOR((LPT6-1)*(IDIM+1)+3) 306 : ZPT5=XCOOR((LPT5-1)*(IDIM+1)+3) 307 : ZPT8=XCOOR((LPT8-1)*(IDIM+1)+3) 308 : ZPT10=XCOOR((LPT10-1)*(IDIM+1)+3) 309 : LON68=SQRT((XPT8-XPT6)**2+(YPT8-YPT6)**2+(ZPT8-ZPT6)**2) 310 : LON510=SQRT((XPT10-XPT5)**2+(YPT10-YPT5)**2+(ZPT10-ZPT5)**2) 311 : IF (LON510.LT.LON68) LONDIA=4*4 312 : ENDIF 313 : IF (JTYPEL.EQ.24) THEN 314 : LPT2=NUMNOE(2) 315 : LPT4=NUMNOE(4) 316 : LPT7=NUMNOE(7) 317 : LPT9=NUMNOE(9) 318 : XPT2=XCOOR((LPT2-1)*(IDIM+1)+1) 319 : XPT4=XCOOR((LPT4-1)*(IDIM+1)+1) 320 : XPT7=XCOOR((LPT7-1)*(IDIM+1)+1) 321 : XPT9=XCOOR((LPT9-1)*(IDIM+1)+1) 322 : YPT2=XCOOR((LPT2-1)*(IDIM+1)+2) 323 : YPT4=XCOOR((LPT4-1)*(IDIM+1)+2) 324 : YPT7=XCOOR((LPT7-1)*(IDIM+1)+2) 325 : YPT9=XCOOR((LPT9-1)*(IDIM+1)+2) 326 : ZPT2=XCOOR((LPT2-1)*(IDIM+1)+3) 327 : ZPT4=XCOOR((LPT4-1)*(IDIM+1)+3) 328 : ZPT7=XCOOR((LPT7-1)*(IDIM+1)+3) 329 : ZPT9=XCOOR((LPT9-1)*(IDIM+1)+3) 330 : LON47=SQRT((XPT7-XPT4)**2+(YPT7-YPT4)**2+(ZPT7-ZPT4)**2) 331 : LON29=SQRT((XPT9-XPT2)**2+(YPT9-YPT2)**2+(ZPT9-ZPT2)**2) 332 : IF (LON29.LT.LON47) LONDIA=4*10 333 : ENDIF 334 : C 335 : C=================================== 336 : C D) Creation des nouveaux elements ! 337 : C=================================== 338 : C On remplit la portion de IPT4 relative aux elements crees a partir 339 : C de la division de l'element IARAF (indice de boucle 1). 340 : C Cette portion de IPT4 contient les colonnes dont la valeur s'etend 341 : C de INELM*(LCOMP-1)+1 a INELM*LCOMP. 342 : DO 5 J=1,INELM 343 : DO 5 I=1,IPT4.NUM(/1) 344 : NTEMP=LIELM(JPOS4-1+NBNN*(J-1)+I) 345 : IF (((JTYPEL.EQ.23).OR.(JTYPEL.EQ.24)).AND.(J.GT.4)) THEN 346 : NTEMP=LIELM(JPOS4-1+NBNN*(J-1)+I+LONDIA) 347 : ENDIF 348 : IF (NTEMP.GT.NBNN) THEN 349 : IPT4.NUM(I,INELM*(LCOMP-1)+J)=NUMNOE(NTEMP-NBNN) 350 : ELSE 351 : IPT4.NUM(I,INELM*(LCOMP-1)+J)=IPT2.NUM(NTEMP,IARAF) 352 : ENDIF 353 : 5 CONTINUE 354 : LCOMP=LCOMP+1 355 : 6 CONTINUE 356 : C 357 : NBPTS=NBPT1 358 : SEGADJ MCOORD 359 : C 360 : C======================================================================= 361 : C Preparation du maillage de relations 362 : C======================================================================= 363 : C================================== 364 : C A) Relations dues aux faces (3D) ! 365 : C================================== 366 : C Tous les noeuds qui restent dans JFARAF sont a creer en tant que 367 : C relations de conformite 368 : C--------------------------------- 369 : C ** A.1 / Initialisation de IPT5 ! 370 : C--------------------------------- 371 : C 1/ Comptage du nombre de noeuds soumis a des relations 372 : NBRELA=0 373 : DO 114 IHF=1,JFARAF(/1) 374 : IF (JNOEFA(IHF,5).EQ.0) GOTO 114 375 : IF (JFARAF(IHF,1).EQ.0) GOTO 114 376 : DO 113 JF=1,JFARAF(/2),2 377 : IF (JFARAF(IHF,JF).NE.0) NBRELA=NBRELA+1 378 : 113 CONTINUE 379 : 114 CONTINUE 380 : C 381 : C 2/ Creation de IPT5 382 : IF (NBRELA.EQ.0) THEN 383 : IPT5=0 384 : GOTO 210 385 : ENDIF 386 : C IF (IPT2.ITYPEL.EQ.14) NBNN=4 387 : C IF (IPT2.ITYPEL.EQ.15) NBNN=8 388 : C IF (IPT2.ITYPEL.EQ.17) NBNN=5 389 : C IF (IPT2.ITYPEL.EQ.24) NBNN=5 390 : NBNN=10 391 : NBELEM=NBRELA 392 : NBSOUS=0 393 : NBREF=0 394 : SEGINI IPT5 395 : IPT5.ITYPEL=48 396 : C 397 : C 3/ Renseignement des noeuds support des relations 398 : NBRELA=0 399 : DO 116 IPF=1,JFARAF(/1) 400 : IF (JNOEFA(IPF,5).EQ.0) GOTO 116 401 : IF (JFARAF(IPF,1).EQ.0) GOTO 116 402 : DO 115 JF=1,JFARAF(/2),2 403 : IF (JFARAF(IPF,JF).EQ.0) GOTO 115 404 : NBRELA=NBRELA+1 405 : IPT5.NUM(1,NBRELA)=JFARAF(IPF,JF) 406 : 115 CONTINUE 407 : 116 CONTINUE 408 : C 409 : C----------------------------------------------------- 410 : C ** A.2 / Recherche des noeuds formant les relations ! 411 : C----------------------------------------------------- 412 : C 1/ Boucle sur l'ensemble des noeuds a creer 413 : DO 200 IARAF=1,MELVA2.VELCHE(/2) 414 : IF (MELVA2.VELCHE(1,IARAF).NE.1) GOTO 200 415 : JCOMPT=0 416 : JPOS5=LPOS5(IPT2.ITYPEL) 417 : DO 190 I=1,NBRAF(IPT2.ITYPEL) 418 : JPOS1=LPOS1(1,I+LPOS2(IPT2.ITYPEL)-1) 419 : JLONG=LPOS1(2,I+LPOS2(IPT2.ITYPEL)-1) 420 : JLISN=LPOS3(IPT2.ITYPEL)+JCOMPT 421 : LTYPNO=JTYPNO(JPOS5-1+I) 422 : IF ((LTYPNO.EQ.0).OR.(LTYPNO.EQ.7)) GOTO 189 423 : C 424 : C 2/ Preparation pour trouver le noeud et la face en question 425 : JTYPEL=IPT2.ITYPEL 426 : JLTEL2=LTEL(2,JTYPEL)-1+LTYPNO 427 : JLDEL1=LDEL(1,JLTEL2) 428 : JTYFAC=IWORK1(JLDEL1) 429 : JLDEL2=LDEL(2,JLTEL2) 430 : JNBSOM=NBSOM(JTYFAC) 431 : JSPOS=NSPOS(JTYFAC) 432 : SEGINI IWORK2 433 : C 434 : C 3/ Classement des 3 sommets par ordre croissant de n° globale 435 : DO 100 IAA=1,JNBSOM 436 : NGLOBA=IPT2.NUM(LFAC(JLDEL2-1+IBSOM(JSPOS-1+IAA)),IARAF) 437 : IWORK2(IAA)=NGLOBA 438 : 100 CONTINUE 439 : NPTA=(XCOOR(/1)/(IDIM+1))+1 440 : NPTB=NPTA+1 441 : NPTC=NPTB+1 442 : DO 110 ICC=1,JNBSOM 443 : IF (IWORK2(ICC).LT.NPTA) THEN 444 : NPTC=NPTB 445 : NPTB=NPTA 446 : NPTA=IWORK2(ICC) 447 : ELSEIF (IWORK2(ICC).LT.NPTB) THEN 448 : NPTC=NPTB 449 : NPTB=IWORK2(ICC) 450 : ELSEIF (IWORK2(ICC).LT.NPTC) THEN 451 : NPTC=IWORK2(ICC) 452 : ENDIF 453 : 110 CONTINUE 454 : C 455 : C 4/ Classement des 3 sommets par ordre croissant de n° locale 456 : NPTA2=ICPR(NPTA,1) 457 : NPTB2=ICPR(NPTB,1) 458 : NPTC2=ICPR(NPTC,1) 459 : IF ((NPTA2.LT.NPTB2).AND.(NPTA2.LT.NPTC2)) THEN 460 : NPTA=NPTA2 461 : NPTB=MIN(NPTB2,NPTC2) 462 : NPTC=MAX(NPTB2,NPTC2) 463 : ENDIF 464 : IF ((NPTB2.LT.NPTA2).AND.(NPTB2.LT.NPTC2)) THEN 465 : NPTA=NPTB2 466 : NPTB=MIN(NPTA2,NPTC2) 467 : NPTC=MAX(NPTA2,NPTC2) 468 : ENDIF 469 : IF ((NPTC2.LT.NPTA2).AND.(NPTC2.LT.NPTB2)) THEN 470 : NPTA=NPTC2 471 : NPTB=MIN(NPTA2,NPTB2) 472 : NPTC=MAX(NPTA2,NPTB2) 473 : ENDIF 474 : C 475 : C 5/ Recherche du numero de la face 476 : NEXIS2=0 477 : DO 120 IEE=1,JPLCOM(NPTA) 478 : MTMP=JPLANS(NPTA,IEE) 479 : JJ1=JNOEFA(MTMP,1) 480 : JJ2=JNOEFA(MTMP,2) 481 : JJ3=JNOEFA(MTMP,3) 482 : IF(JJ1.EQ.NPTA.AND.JJ2.EQ.NPTB.AND.JJ3.EQ.NPTC) THEN 483 : NEXIS2=IEE 484 : ENDIF 485 : 120 CONTINUE 486 : JNUMFA=JPLANS(NPTA,NEXIS2) 487 : C 488 : C 6/ Recherche du numero global du point 489 : IF (JNOEFA(JNUMFA,5).EQ.0) GOTO 189 490 : IF (JTYFAC.EQ.8) INOEGL=JFARAF(JNUMFA,1) 491 : IF (JTYFAC.NE.8) THEN 492 : NEXIS3=0 493 : NEXIS4=0 494 : XCO2=0.25 495 : DO 140 KBB=1,JLONG 496 : XCO1=XCOEFF(JPOS1-1+KBB) 497 : IF (XCO1.EQ.XCO2) NEXIS3=KBB 498 : 140 CONTINUE 499 : IF (NEXIS3.NE.0) THEN 500 : NGLOB=IPT2.NUM(LISNOE(JLISN-1+NEXIS3),IARAF) 501 : DO 150 KAA=2,JFARAF(/2),2 502 : IF (JFARAF(JNUMFA,KAA).EQ.NGLOB) NEXIS4=KAA 503 : 150 CONTINUE 504 : INOEGL=JFARAF(JNUMFA,NEXIS4-1) 505 : ELSE 506 : DO 160 KAA=2,JFARAF(/2),2 507 : IF (JFARAF(JNUMFA,KAA).EQ.0) NEXIS4=KAA 508 : 160 CONTINUE 509 : INOEGL=JFARAF(JNUMFA,NEXIS4-1) 510 : ENDIF 511 : ENDIF 512 : C 513 : C------------------------------ 514 : C ** A.3 / Remplissage de IPT5 ! 515 : C------------------------------ 516 : C 1/ Recherche de la position du point dans IPT5 517 : NEXIS5=0 518 : DO 170 IGG=1,IPT5.NUM(/2) 519 : IF (INOEGL.EQ.IPT5.NUM(1,IGG)) NEXIS5=IGG 520 : 170 CONTINUE 521 : IF (NEXIS5.EQ.0) GOTO 189 522 : C 523 : C 2/ Renseignement des points formant les relations 524 : DO 180 IHH=1,JLONG 525 : IPT5.NUM(1+IHH,NEXIS5)=IPT2.NUM(LISNOE(JLISN-1+IHH),IARAF) 526 : 180 CONTINUE 527 : IF (JLONG.EQ.4) IPT5.NUM(10,NEXIS5)=3 528 : IF (JLONG.EQ.5) IPT5.NUM(10,NEXIS5)=4 529 : IF (JLONG.EQ.8) THEN 530 : IF (JPOS1.EQ.16) IPT5.NUM(10,NEXIS5)=5 531 : IF (JPOS1.EQ.24) IPT5.NUM(10,NEXIS5)=6 532 : ENDIF 533 : 189 CONTINUE 534 : JCOMPT=JCOMPT+JLONG 535 : 190 CONTINUE 536 : 200 CONTINUE 537 : 210 CONTINUE 538 : C 539 : C=================================== 540 : C B) Relations dues aux aretes (2D) ! 541 : C=================================== 542 : C On cree un maillage IPT6 contenant tous les noeuds soumis a des 543 : C relations 544 : C 545 : SEGSUP KARET2 546 : NBNDS=KARETE(/1) 547 : NCOL=KARETE(/2) 548 : SEGINI KARET2 549 : ILPL=LPL(IPT2.ITYPEL) 550 : ILPT=LPT(IPT2.ITYPEL) 551 : DO 317 J=1,IPT4.NUM(/2) 552 : DO 317 K=1,ILPL*2-1,2 553 : NPTA=IPT4.NUM(KSEGM(ILPT+K-1),J) 554 : NPTB=IPT4.NUM(KSEGM(ILPT+K),J) 555 : IF((NPTA.GT.NBPT0).OR.(NPTB.GT.NBPT0)) THEN 556 : GOTO 317 557 : ENDIF 558 : NMIN=MIN(ICPR(NPTA,1),ICPR(NPTB,1)) 559 : NMAX=MAX(ICPR(NPTA,1),ICPR(NPTB,1)) 560 : NEXIST=0 561 : DO 316 I=1,MAX(1,KARPOS(NMIN)) 562 : IF (KARETE(NMIN,I).EQ.NMAX) THEN 563 : KARET2(NMIN,I)=KARET2(NMIN,I)+1 564 : ENDIF 565 : 316 CONTINUE 566 : 317 CONTINUE 567 : 568 : C 1/ Comptage du nombre de noeuds soumis a des relations 569 : NBELEM=0 570 : DO 27 J=1,KMILIE(/2) 571 : DO 27 I=1,KMILIE(/1) 572 : IF (KARET2(I,J).EQ.0) GOTO 27 573 : IF (KMILIE(I,J).GT.0) NBELEM=NBELEM+1 574 : 27 CONTINUE 575 : C 576 : C 2/ Creation de IPT6 577 : IPT6=0 578 : IF (NBELEM.EQ.0) GOTO 999 579 : NBNN=5 580 : NBREF=0 581 : NBSOUS=0 582 : SEGINI IPT6 583 : IPT6.ITYPEL=48 584 : C 585 : C 3/ Renseignement des noeuds support des relations 586 : DO 28 J=1,KMILIE(/2) 587 : DO 28 I=1,KMILIE(/1) 588 : IF (KARET2(I,J).EQ.0) GOTO 28 589 : IF (KMILIE(I,J).GT.0) THEN 590 : NBREF=NBREF+1 591 : IPT6.NUM(1,NBREF)=KMILIE(I,J) 592 : ENDIF 593 : 28 CONTINUE 594 : C 595 : C 4/ Recherche des noeuds formant les relations 596 : DO 24 IARAF=1,MELVA2.VELCHE(/2) 597 : IF (MELVA2.VELCHE(1,IARAF).NE.1) GOTO 24 598 : JCOMPT=0 599 : JPOS5=LPOS5(IPT2.ITYPEL) 600 : DO 23 I=1,NBRAF(IPT2.ITYPEL) 601 : JPOS1=LPOS1(1,I+LPOS2(IPT2.ITYPEL)-1) 602 : JLONG=LPOS1(2,I+LPOS2(IPT2.ITYPEL)-1) 603 : JLISN=LPOS3(IPT2.ITYPEL)+JCOMPT 604 : LTYPNO=JTYPNO(JPOS5-1+I) 605 : IF (LTYPNO.NE.0) GOTO 22 606 : NPTA=IPT2.NUM(LISNOE(JLISN),IARAF) 607 : NPTB=IPT2.NUM(LISNOE(JLISN+1),IARAF) 608 : NMIN=MIN(ICPR(NPTA,1),ICPR(NPTB,1)) 609 : NMAX=MAX(ICPR(NPTA,1),ICPR(NPTB,1)) 610 : DO 29 K=1,MAX(1,KARPOS(NMIN)) 611 : IF (KARETE(NMIN,K).EQ.NMAX) NEXIST=K 612 : 29 CONTINUE 613 : IF (KARET2(NMIN,NEXIST).EQ.0) GOTO 22 614 : IF (KMILIE(NMIN,NEXIST).EQ.0) GOTO 22 615 : NEXIS5=0 616 : DO 20 MM=1,IPT6.NUM(/2) 617 : INOEGL=KMILIE(NMIN,NEXIST) 618 : INRELA=IPT6.NUM(1,MM) 619 : IF (INOEGL.EQ.INRELA) NEXIS5=MM 620 : 20 CONTINUE 621 : C 622 : C 5/ Renseignement des noeuds formant les relations 623 : DO 21 IHH=1,JLONG 624 : IPT6.NUM(1+IHH,NEXIS5)=IPT2.NUM(LISNOE(JLISN-1+IHH),IARAF) 625 : 21 CONTINUE 626 : IF (JLONG.EQ.2) IPT6.NUM(5,NEXIS5)=1 627 : IF (JLONG.EQ.3) IPT6.NUM(5,NEXIS5)=2 628 : 22 CONTINUE 629 : JCOMPT=JCOMPT+JLONG 630 : 23 CONTINUE 631 : 24 CONTINUE 632 : 633 : 444 CONTINUE 634 : 635 : C 636 : C============================================ 637 : C C) Creation du maillage de relations final ! 638 : C============================================ 639 : IF (IPT5.EQ.0) THEN 640 : IPT7=IPT6 641 : GOTO 999 642 : ENDIF 643 : NBELEM=IPT5.NUM(/2)+IPT6.NUM(/2) 644 : C NBNN=MAX(IPT5.NUM(/1),IPT6.NUM(/1)) 645 : NBNN=10 646 : NBREF=0 647 : NBSOUS=0 648 : SEGINI IPT7 649 : IPT7.ITYPEL=48 650 : DO 42 NEO=1,IPT5.NUM(/2) 651 : DO 42 MOR=1,10 652 : IPT7.NUM(MOR,NEO)=IPT5.NUM(MOR,NEO) 653 : IPT7.ICOLOR(NEO)=IPT5.ICOLOR(NEO) 654 : 42 CONTINUE 655 : NN5 = IPT5.NUM(/2) 656 : DO 43 NEO=IPT5.NUM(/2)+1,IPT6.NUM(/2)+IPT5.NUM(/2) 657 : DO 43 MOR=1,IPT6.NUM(/1) 658 : IF (MOR.LT.IPT6.NUM(/1)) THEN 659 : IPT7.NUM(MOR,NEO)=IPT6.NUM(MOR,NEO-NN5) 660 : ELSE 661 : IPT7.NUM(10,NEO)=IPT6.NUM(MOR,NEO-NN5) 662 : ENDIF 663 : IPT7.ICOLOR(NEO)=IPT6.ICOLOR(NEO-NN5) 664 : 43 CONTINUE 665 : C======================================================================= 666 : C Fin du programme 667 : C======================================================================= 668 : 999 CONTINUE 669 : RETURN 670 : END 671 : 672 : 673 : 674 : 675 : 676 : 677 :
© Cast3M 2003 - All rights reserved.
Disclaimer