1 SUBROUTINE ELIMIN 2C 3 INCLUDE 'com_coor.f' 4 INCLUDE 'com_faces.f' 5 INCLUDE 'com_options.f' 6C 7 INTEGER NODE(27),IPP(4,4,9),NUMF(6),IP8(4,6),NOD(8),IP827(8,8) 8 & ,IP9(4,4),IP4(2,4),IP4T(3,4),IP6(3,4) 9 & ,NFPR(6),IBID(28),NONO(8),IFAFA(6) 10 REAL*4 PRES(6),XBID(28) 11 EQUIVALENCE(IBID(2),NODE) 12 EQUIVALENCE(IBID,XBID) 13 EQUIVALENCE(NOD,NODE) 14 LOGICAL*4 ZTRI,EGAL 15 CHARACTER*128 CBID 16 DATA IPP / 25,18, 6,13, 25,17, 1, 9, 17 & 25, 9, 2,18, 25,13, 5,17, 18 & 23,14, 7,15, 23,16, 5,13, 19 & 23,13, 6,14, 23,15, 8,16, 20 & 22,19, 7,14, 22,18, 2,10, 21 & 22,10, 3,19, 22,14, 6,18, 22 & 21,11, 3,10, 21, 9, 1,12, 23 & 21,12, 4,11, 21,10, 2, 9, 24 & 26,15, 7,19, 26,11, 4,20, 25 & 26,20, 8,15, 26,19, 3,11, 26 & 24,16, 8,20, 24,12, 1,17, 27 & 24,17, 5,16, 24,20, 4,12, 28 & 27,23,14,22, 27,24,16,23, 29 & 27,21,12,24, 27,22,10,21, 30 & 27,26,15,23, 27,23,13,25, 31 & 27,25, 9,21, 27,21,11,26, 32 & 27,22,19,26, 27,26,20,24, 33 & 27,24,17,25, 27,25,18,22/ 34 DATA NUMF / 3,6,2,5,4,1 / 35 DATA IP4 / 4,1 , 2,3 , 1,2 , 3,4 / 36 DATA IP4T / 1,3,2 37 & ,2,3,4 38 & ,1,2,4 39 & ,1,4,3 / 40 DATA IP6 / 6,1,4 41 & ,4,2,5 42 & ,5,3,6 43 & ,4,5,6 / 44 DATA IP8 / 5,8,4,1 45 & ,2,3,7,6 46 & ,5,1,2,6 47 & ,4,8,7,3 48 & ,1,4,3,2 49 & ,8,5,6,7 / 50 DATA IP9 / 9,5,2,6 51 & ,9,6,3,7 52 & ,9,7,4,8 53 & ,9,8,1,5 / 54 DATA IP827 / 1, 9,21,12,17,25,27,24 55 & , 9, 2,10,21,25,18,22,27 56 & ,21,10, 3,11,27,22,19,26 57 & ,12,21,11, 4,24,27,26,20 58 & ,17,25,27,24, 5,13,23,16 59 & ,25,18,22,27,13, 6,14,23 60 & ,27,22,19,26,23,14, 7,15 61 & ,24,27,26,20,16,23,15, 8 / 62 DATA US3 / .33333333333333333333 / 63C 64 IF (ISTDOUT.EQ.0) THEN 65 IF (ILANG.EQ.0) THEN 66 PRINT*,'Fin de lecture des noeuds. Lecture des �l�ments...' 67 ELSE 68 PRINT*,'End of nodes reading. Beginning elements reading..' 69 ENDIF 70 ENDIF 71 DO I=1,NUMNP 72 ISD(I) = 0 73 INDEX(I) = 0 74 ENDDO 75 IF (NBCORN.NE.0) THEN 76 DO I=1,NUMNP 77 RAYON2(I) = 0. 78 ENDDO 79 ENDIF 80C 81 NDSMAX = 0 82 IFACEXT = 0 83 IVOL = 0 84 NEL = 0 85 ICPT = 0 86 NFC = 0 87 IGROUP = 0 88 NUMSD = 0 89 NUMEL = 0 90 NUMNP2 = NUMNP 91 MODTET = 0 92 ICOIN = 0 93 IVIEUX = 0 94 1000 CONTINUE 95 IF (ICOURB.EQ.-5) THEN 96 IGROUP = 1 97 ICPT0 = 1 98 NELEM = NELDEL 99 NEL = NELEM 100 IF (ICOURXYZ.EQ.0) THEN 101 NDS = 3 102 DO N=1,NEL 103 NN = (N-1)*NDS+1 104 CALL QUALITY(NODEL(NN),N,NDS) 105 ENDDO 106 ELSEIF(ICOURXYZ.EQ.1) THEN 107 NDS = 2 108 ELSEIF(ICOURXYZ.EQ.2) THEN 109C 110C trait� dans elimin3 111C 112 ICPT = NF 113 GOTO 9999 114 ENDIF 115 NDSEL = NDS 116 ELSE 117 IF (IARCH.EQ.0) THEN 118 READ(IFAC,END=9999,ERR=9999) N,NTYP,NELEM,NDS 119 ELSE 120 CALL litrecbin(IFAC,NODE,LLL,0) 121 IF (LLL.LE.0) GOTO 9999 122 N = NODE(1) 123 NTYP = NODE(2) 124 NELEM = NODE(3) 125 NDS = NODE(4) 126 ENDIF 127 IF (NELEM.LE.0) GOTO 1000 128 NDSEL = NDS 129 IF (NDS*(NEL+NELEM).GT.NEMAX) CALL TROPDEPOINTS(NEL+NELEM,NDS,0) 130CC maintenant on lit le no de sd directement 131CCC IGROUP = IGROUP+1 132 IGROUP = N 133 ICPT0 = ICPT+1 134 NEL = NEL + NELEM 135 ENDIF 136 NMODP = MAX(1,NELEM/50) 137 IF (NMODP.EQ.1) NMODP = NFMAX 138C 139C element a 27 noeuds 140C 141 IF (NDS.GE.27) THEN 142 IF (ILANG.EQ.0) THEN 143 ELEMENTS = 'Hexa�dres 27 noeuds' 144 ELSE 145 ELEMENTS = 'Hexaedrons 27 nodes' 146 ENDIF 147 MODTET = 40 148 LELEM = 19 149 IOPREF = 1 150 NBFAC = 24 151 NFEXT = 6 152 ICOIN = 1 153 IOPT = 0 154 DO N=1,NELEM 155 IF (MOD(N,NMODP).EQ.0.AND.ISTDOUT.EQ.0) 156 & CALL MYFLUSH(100.*REAL(N)/REAL(NELEM)) 157 IF (IARCH.EQ.0) THEN 158 READ(IFAC) NE,(NODE(J),J=1,NDS) 159 ELSE 160 CALL litrecbin(IFAC,IBID,LLL,0) 161 ENDIF 162 IMAT = IREF3(NODE(27)) 163 NUMEL = NUMEL+1 164 CALL ECRNOD(NODE,NODEL,NUMEL,NDS) 165 DO I=1,NDS 166 IF (ISD(NODE(I)).EQ.0) THEN 167 ISD(NODE(I)) = IGROUP 168 ELSEIF(ISD(NODE(I)).NE.IGROUP) THEN 169 ISD(NODE(I)) = -1 170 ENDIF 171 ENDDO 172 DO I=9,20 173 INDEX(NODE(I)) = INDEX(NODE(I))+1 174 ENDDO 175 DO I=21,26 176 INDEX(NODE(I)) = -1 177 ENDDO 178C 179C Decomposition de l'element a 27 noeuds en 24 faces exterieures 180C 181 DO IFCC = 1,NFEXT 182 IF (IELIMI.EQ.0) IOPT = 1 183 NUMFAC = NUMF(IFCC) 184 DO K=1,4 185 CALL TRIAGE(NODE(IPP(1,K,IFCC)) 186 & ,NODE(IPP(2,K,IFCC)) 187 & ,NODE(IPP(3,K,IFCC)) 188 & ,NODE(IPP(4,K,IFCC)) 189 & ,ICPT,ICPT0,NFC,IOPT,NUMEL,NUMFAC,NODE(27)) 190 ISD2(ICPT) = IGROUP 191 ISD3(ICPT) = IMAT 192 ENDDO 193 IF (NFC.GT.NFCM) CALL COMPAC(ICPT,ICPT0,NFC) 194 ENDDO 195C 196C Calcul des rayons de courbure et puissances dioptriques pour les yeux 197C 198 DO K=1,6 199 IF (ICOR(NODE(20+K)).NE.0) CALL RAYFAC(NODE,K) 200 ENDDO 201C 202C Decomposition en tetraedres pour isosurfaces 203C 204 DO K=1,8 205 DO KK=1,8 206 NONO(KK) = NODE(IP827(KK,K)) 207 ENDDO 208 CALL TETRA(NONO) 209 ENDDO 210 ENDDO 211C 212C 4 noeuds.... 213C 214 ELSEIF(NDS.EQ.4) THEN 215C 216C ... 2d --> quadrangles 217C 218 IF (I2D.NE.0.OR.NTYP.EQ.3) THEN 219 NBFAC = 1 220 IPREFC = -1 221 IOPREF = 0 222 IF (I2D.EQ.0) THEN 223 ELEMENTS = 'Quadrangles 3D Q1' 224 ELSE 225 ELEMENTS = 'Quadrangles 2D Q1' 226 ENDIF 227 LELEM = 17 228 DO N=1,NELEM 229 IF (MOD(N,NMODP).EQ.0.AND.ISTDOUT.EQ.0) 230 & CALL MYFLUSH(100.*REAL(N)/REAL(NELEM)) 231 ICPT = ICPT+1 232 IF (IARCH.EQ.0) THEN 233 IF (IVIEUX.EQ.0) THEN 234 READ(IFAC,END=3131,ERR=3131) NE,(NOD(J),J=1,NDS) 235 & ,NPR,(NFPR(K),PRES(K),K=1,NPR) 236 GOTO 3132 237 3131 NPR = 0 238 IVIEUX = 1 239 3132 NUMEL = NUMEL+1 240 ELSE 241 READ(IFAC) NE,(NOD(J),J=1,NDS) 242 NUMEL = NUMEL+1 243 ENDIF 244 ELSE 245 CALL litrecbin(IFAC,IBID,LLL,0) 246cc NE = IBID(1) 247 IF (LLL/4.GT.NDS+1) THEN 248 NPR = IBID(NDS+2) 249 IF (NPR.GT.0) THEN 250 DO J=1,NPR 251 NFPR(J) = IBID(NDS+1+2*J) 252 PRES(J) = XBID(NDS+2*(J+1)) 253 ENDDO 254 ENDIF 255 ELSE 256 NPR = 0 257 ENDIF 258 NUMEL = NUMEL+1 259 ENDIF 260 IF (NPR.GT.0) THEN 261 DO J=1,NPR 262 IF (PRES(J).EQ.9999.) THEN 263 IR = 5 264 ELSEIF(PRES(J).EQ.-9999.) THEN 265 IR = 1 266 ELSE 267 IR = MOD(NINT(PRES(J)),16) 268 IF (IR.LT.0) IR = IR+16 269 ENDIF 270 IREF(NOD(IP4(1,NFPR(J)))) = IR 271 IREF(NOD(IP4(2,NFPR(J)))) = IR 272 ENDDO 273 ENDIF 274 DO J=1,NDS 275 ICLAS(J,ICPT) = NOD(J) 276 ENDDO 277 CALL ECRNOD(ICLAS(1,ICPT),NODEL,NUMEL,NDS) 278 DO J=1,NDS 279 NN = ICLAS(J,ICPT) 280 INDEX(NN) = INDEX(NN)+1 281 IF (INDEX(NN).LE.NVMX) NVOI(INDEX(NN),NN) = ICPT 282 IF (ISD(NN).EQ.0) THEN 283 ISD(NN) = IGROUP 284 ELSEIF(ISD(NN).NE.IGROUP) THEN 285 ISD(NN) = -1 286 ENDIF 287 ENDDO 288 ICLAS(5,ICPT) = NUMEL 289 ICLAS(6,ICPT) = 0 290 NUMNP2 = NUMNP2+1 291 X(NUMNP2) = .25*( X(ICLAS(1,ICPT))+X(ICLAS(2,ICPT)) 292 & +X(ICLAS(3,ICPT))+X(ICLAS(4,ICPT)) ) 293 Y(NUMNP2) = .25*( Y(ICLAS(1,ICPT))+Y(ICLAS(2,ICPT)) 294 & +Y(ICLAS(3,ICPT))+Y(ICLAS(4,ICPT)) ) 295 Z(NUMNP2) = .25*( Z(ICLAS(1,ICPT))+Z(ICLAS(2,ICPT)) 296 & +Z(ICLAS(3,ICPT))+Z(ICLAS(4,ICPT)) ) 297 DEPX(NUMNP2) = .25*(DEPX(ICLAS(1,ICPT))+DEPX(ICLAS(2,ICPT)) 298 & +DEPX(ICLAS(3,ICPT))+DEPX(ICLAS(4,ICPT))) 299 DEPY(NUMNP2) = .25*(DEPY(ICLAS(1,ICPT))+DEPY(ICLAS(2,ICPT)) 300 & +DEPY(ICLAS(3,ICPT))+DEPY(ICLAS(4,ICPT))) 301 DEPZ(NUMNP2) = .25*(DEPZ(ICLAS(1,ICPT))+DEPZ(ICLAS(2,ICPT)) 302 & +DEPZ(ICLAS(3,ICPT))+DEPZ(ICLAS(4,ICPT))) 303 ICLAS(7,ICPT) = NUMNP2 304 ISD2(ICPT) = IGROUP 305 ENDDO 306C 307C ... 3d --> tetraedres 308C 309 ELSE 310 IOPREF = 1 311 NBFAC = 4 312 MODTET = 1 313 IF (ILANG.EQ.0) THEN 314 ELEMENTS = 'Tetra�dres' 315 LELEM = 10 316 ELSE 317 ELEMENTS = 'Tetraedrons' 318 LELEM = 11 319 ENDIF 320 IF (IELIMI.EQ.0) THEN 321 IOPT = 4 322 ELSE 323 IOPT = 0 324 ENDIF 325 IF (NEL.EQ.NELEM) CALL CHECKFACEXT(CBID,LBID,IFACEXT) 326 IBORLU = 0 327 4436 CONTINUE 328 DO N=1,NELEM 329 IF (MOD(N,NMODP).EQ.0.AND.ISTDOUT.EQ.0) 330 & CALL MYFLUSH(100.*REAL(N)/REAL(NELEM)) 331 IF (IFACEXT.GT.0.AND.IBORLU.EQ.0) THEN 332 READ(IFACEXT,*,END=4434,ERR=4434) 333 & NB,NFA,(IFAFA(I),I=1,NFA) 334 IBORLU = 1 335 GOTO 4435 336 4434 IF (N.EQ.1) THEN 337 IFACEXT = -1 338 IF (ILANG.EQ.0) THEN 339 PRINT*,'*** Erreur � la lecture de '//CBID(1:LBID)// 340 & ' - �l�ment',N 341 ELSE 342 PRINT*,'*** Error reading '//CBID(1:LBID)// 343 & ' - element',N 344 ENDIF 345 NUMEL = 0 346 NUMNP2 = NUMNP 347 GOTO 4436 348 ENDIF 349 4435 CONTINUE 350 ENDIF 351 IF (IARCH.EQ.0) THEN 352 READ(IFAC) NE,N1,N2,N3,N4 353 ELSE 354 CALL litrecbin(IFAC,IBID,LLL,0) 355cc NE = IBID(1) 356 N1 = IBID(2) 357 N2 = IBID(3) 358 N3 = IBID(4) 359 N4 = IBID(5) 360 ENDIF 361 IVOL = IVOL+1 362 NTET(1,IVOL) = N1 363 NTET(2,IVOL) = N2 364 NTET(3,IVOL) = N3 365 NTET(4,IVOL) = N4 366 NUMEL = NUMEL+1 367 CALL ECRNOD(NTET(1,IVOL),NODEL,NUMEL,4) 368 NUMNP2 = NUMNP2+1 369 X(NUMNP2) = .25*( X(N1)+X(N2)+X(N3)+X(N4) ) 370 Y(NUMNP2) = .25*( Y(N1)+Y(N2)+Y(N3)+Y(N4) ) 371 Z(NUMNP2) = .25*( Z(N1)+Z(N2)+Z(N3)+Z(N4) ) 372 DEPX(NUMNP2) = .25*( DEPX(N1)+DEPX(N2)+DEPX(N3)+DEPX(N4) ) 373 DEPY(NUMNP2) = .25*( DEPY(N1)+DEPY(N2)+DEPY(N3)+DEPY(N4) ) 374 DEPZ(NUMNP2) = .25*( DEPZ(N1)+DEPZ(N2)+DEPZ(N3)+DEPZ(N4) ) 375 IF (IFACEXT.GT.0) THEN 376 IF (NB.EQ.NUMEL) THEN 377 IBORLU = 0 378 DO K=1,4 379 IFLAG = 0 380 DO KK=1,NFA 381 IF (IFAFA(KK).EQ.K) IFLAG=KK 382 ENDDO 383 IF (IFLAG.NE.0) THEN 384 IOPT = 0 385 CALL ORDON2(NTET(IP4T(1,K),IVOL) 386 & ,NTET(IP4T(2,K),IVOL) 387 & ,NTET(IP4T(3,K),IVOL),M1,M2,M3) 388 CALL TRIAGE(M1,M2,M3,0,ICPT,ICPT0,NFC,IOPT,NUMEL,K 389 & ,NUMNP2) 390 ISD2(ICPT) = IGROUP 391 ENDIF 392 ENDDO 393 ENDIF 394 ELSE 395 ICPT00 = ICPT 396 CALL ORDON2(N1,N4,N3,M1,M2,M3) 397 CALL TRIAGE(M1,M2,M3,0,ICPT,ICPT0,NFC,IOPT,NUMEL,4,NUMNP2) 398 ISD2(ICPT) = IGROUP 399 CALL ORDON2(N1,N2,N4,M1,M2,M3) 400 CALL TRIAGE(M1,M2,M3,0,ICPT,ICPT0,NFC,IOPT,NUMEL,3,NUMNP2) 401 ISD2(ICPT) = IGROUP 402 CALL ORDON2(N2,N3,N4,M1,M2,M3) 403 CALL TRIAGE(M1,M2,M3,0,ICPT,ICPT0,NFC,IOPT,NUMEL,2,NUMNP2) 404 ISD2(ICPT) = IGROUP 405 CALL ORDON2(N1,N3,N2,M1,M2,M3) 406 CALL TRIAGE(M1,M2,M3,0,ICPT,ICPT0,NFC,IOPT,NUMEL,1,NUMNP2) 407 ISD2(ICPT) = IGROUP 408 DO K=ICPT00,ICPT 409 DO J=1,NDS 410 NN = ICLAS(J,K) 411 INDEX(NN) = INDEX(NN)+1 412 IF (INDEX(NN).LE.NVMX) NVOI(INDEX(NN),NN) = K 413 JJ = ICLAS(J,K) 414 IF (ISD(JJ).EQ.0) THEN 415 ISD(JJ) = IGROUP 416 ELSEIF(ISD(JJ).NE.IGROUP) THEN 417 ISD(JJ) = -1 418 ENDIF 419 ENDDO 420 ENDDO 421 IF (NFC.GT.NFCM) CALL COMPAC(ICPT,ICPT0,NFC) 422 ENDIF 423 ENDDO 424 ENDIF 425C 426C Triangles P2 427C 428 ELSEIF(NDS.EQ.6) THEN 429C 430C c'est comme ca qu'on dit qu'on fait du Q2 (ou P2). idiot 431C 432 MODTET = 40 433 IOPREF = 1 434 NBFAC = 4 435 ICOIN = 1 436 IF (I2D.EQ.0) THEN 437 ELEMENTS = 'Triangles 3D Q2' 438 ELSE 439 ELEMENTS = 'Triangles 2D Q2' 440 ENDIF 441 LELEM = 15 442 DO N=1,NELEM 443 IF (MOD(N,NMODP).EQ.0.AND.ISTDOUT.EQ.0) 444 & CALL MYFLUSH(100.*REAL(N)/REAL(NELEM)) 445 IF (IARCH.EQ.0) THEN 446 READ(IFAC) NE,(NODE(J),J=1,NDS) 447 ELSE 448 CALL litrecbin(IFAC,IBID,LLL,0) 449cc NE = IBID(1) 450 ENDIF 451 NUMEL = NUMEL+1 452 CALL ECRNOD(NODE,NODEL,NUMEL,NDS) 453 DO I=1,4 454 ICPT = ICPT+1 455 DO J=1,3 456 NN = NODE(IP6(J,I)) 457 INDEX(NN) = INDEX(NN)+1 458 IF (INDEX(NN).LE.NVMX) NVOI(INDEX(NN),NN) = ICPT 459 IF (ISD(NN).EQ.0) THEN 460 ISD(NN) = IGROUP 461 ELSEIF(ISD(NN).NE.IGROUP) THEN 462 ISD(NN) = -1 463 ENDIF 464 ICLAS(J,ICPT) = NN 465 ENDDO 466 ICLAS(4,ICPT) = ICLAS(3,ICPT) 467 ICLAS(5,ICPT) = NUMEL 468 ICLAS(6,ICPT) = 0 469 NUMNP2 = NUMNP2+1 470 X(NUMNP2) = X(ICLAS(1,ICPT)) 471 Y(NUMNP2) = Y(ICLAS(1,ICPT)) 472 Z(NUMNP2) = Z(ICLAS(1,ICPT)) 473 DEPX(NUMNP2) = DEPX(ICLAS(1,ICPT)) 474 DEPY(NUMNP2) = DEPY(ICLAS(1,ICPT)) 475 DEPZ(NUMNP2) = DEPZ(ICLAS(1,ICPT)) 476 ICLAS(7,ICPT) = NUMNP2 477 ISD2(ICPT) = IGROUP 478 ENDDO 479 ENDDO 480C 481C Triangles P1 et segments 482C 483 ELSEIF(NDS.EQ.3.OR.NDS.EQ.2) THEN 484ccc? IOPREF = 0 485 IOPREF = 1 486 NBFAC = 1 487 IF (NDS.EQ.3) THEN 488 IF (I2D.EQ.0) THEN 489 ELEMENTS = 'Triangles 3D P1' 490 ELSE 491 ELEMENTS = 'Triangles 2D P1' 492 ENDIF 493 LELEM = 15 494 ELSE 495 IF (I2D.EQ.0) THEN 496 ELEMENTS = 'Segments 3D' 497 ELSE 498 ELEMENTS = 'Segments 2D' 499 ENDIF 500 LELEM = 11 501 ENDIF 502 DO N=1,NELEM 503 IF (MOD(N,NMODP).EQ.0.AND.ISTDOUT.EQ.0) 504 & CALL MYFLUSH(100.*REAL(N)/REAL(NELEM)) 505 ICPT = ICPT+1 506 IF (ICOURB.NE.-5) THEN 507 IF (IARCH.EQ.0) THEN 508 READ(IFAC) NE,(ICLAS(J,ICPT),J=1,NDS) 509 ELSE 510 CALL litrecbin(IFAC,IBID,LLL,0) 511cc NE = IBID(1) 512 DO J=1,NDS 513 ICLAS(J,ICPT) = IBID(J+1) 514 ENDDO 515 ENDIF 516 ENDIF 517 NUMEL = NUMEL+1 518 CALL ECRNOD(ICLAS(1,ICPT),NODEL,NUMEL,NDS) 519 ICLAS(4,ICPT) = ICLAS(3,ICPT) 520 ICLAS(5,ICPT) = NUMEL 521 ICLAS(6,ICPT) = 0 522 NUMNP2 = NUMNP2+1 523 IF (NDS.EQ.3) THEN 524 X(NUMNP2) = ( X(ICLAS(1,ICPT))+X(ICLAS(2,ICPT)) 525 & + X(ICLAS(3,ICPT)) )*US3 526 Y(NUMNP2) = ( Y(ICLAS(1,ICPT))+Y(ICLAS(2,ICPT)) 527 & + Y(ICLAS(3,ICPT)) )*US3 528 Z(NUMNP2) = ( Z(ICLAS(1,ICPT))+Z(ICLAS(2,ICPT)) 529 & + Z(ICLAS(3,ICPT)) )*US3 530 DEPX(NUMNP2) = ( DEPX(ICLAS(1,ICPT))+DEPX(ICLAS(2,ICPT)) 531 & + DEPX(ICLAS(3,ICPT)) )*US3 532 DEPY(NUMNP2) = ( DEPY(ICLAS(1,ICPT))+DEPY(ICLAS(2,ICPT)) 533 & + DEPY(ICLAS(3,ICPT)) )*US3 534 DEPZ(NUMNP2) = ( DEPZ(ICLAS(1,ICPT))+DEPZ(ICLAS(2,ICPT)) 535 & + DEPZ(ICLAS(3,ICPT)) )*US3 536 ELSE 537 X(NUMNP2) = ( X(ICLAS(1,ICPT))+X(ICLAS(2,ICPT)) )*.5 538 Y(NUMNP2) = ( Y(ICLAS(1,ICPT))+Y(ICLAS(2,ICPT)) )*.5 539 Z(NUMNP2) = ( Z(ICLAS(1,ICPT))+Z(ICLAS(2,ICPT)) )*.5 540 DEPX(NUMNP2) = (DEPX(ICLAS(1,ICPT))+DEPX(ICLAS(2,ICPT)))*.5 541 DEPY(NUMNP2) = (DEPY(ICLAS(1,ICPT))+DEPY(ICLAS(2,ICPT)))*.5 542 DEPZ(NUMNP2) = (DEPZ(ICLAS(1,ICPT))+DEPZ(ICLAS(2,ICPT)))*.5 543 ICLAS(3,ICPT) = ICLAS(2,ICPT) 544 ENDIF 545 ICLAS(7,ICPT) = NUMNP2 546 ISD2(ICPT) = IGROUP 547 DO J=1,NDS 548 NN = ICLAS(J,ICPT) 549 INDEX(NN) = INDEX(NN)+1 550 IF (INDEX(NN).LE.NVMX) NVOI(INDEX(NN),NN) = ICPT 551 JJ = ICLAS(J,ICPT) 552 IF (ISD(JJ).EQ.0) THEN 553 ISD(JJ) = IGROUP 554 ELSEIF(ISD(JJ).NE.IGROUP) THEN 555 ISD(JJ) = -1 556 ENDIF 557 ENDDO 558 ENDDO 559C 560C Briques a 8 noeuds 561C 562 ELSEIF(NDS.EQ.8) THEN 563 IOPREF = 0 564 NBFAC = 6 565 MODTET = 5 566 IF (ILANG.EQ.0) THEN 567 ELEMENTS = 'Hexa�dres 8 noeuds' 568 ELSE 569 ELEMENTS = 'Hexaedrons 8 nodes' 570 ENDIF 571 LELEM = 18 572 IF (IELIMI.EQ.0) THEN 573 IOPT = 8 574 ELSE 575 IOPT = 0 576 ENDIF 577 IF (NEL.EQ.NELEM) CALL CHECKFACEXT(CBID,LBID,IFACEXT) 578 IBORLU = 0 579 3436 CONTINUE 580 DO N=1,NELEM 581 IF (MOD(N,NMODP).EQ.0.AND.ISTDOUT.EQ.0) 582 & CALL MYFLUSH(100.*REAL(N)/REAL(NELEM)) 583 IF (IFACEXT.GT.0.AND.IBORLU.EQ.0) THEN 584 READ(IFACEXT,*,END=3434,ERR=3434) NB,NFA,(IFAFA(I),I=1,NFA) 585 IBORLU = 1 586 GOTO 3435 587 3434 IF (N.EQ.1) THEN 588 IFACEXT = -1 589 IF (ILANG.EQ.0) THEN 590 PRINT*,'*** Erreur � la lecture de '//CBID(1:LBID)// 591 & ' - �l�ment',N 592 ELSE 593 PRINT*,'*** Error reading '//CBID(1:LBID)// 594 & ' - element',N 595 ENDIF 596 NUMEL = 0 597 NUMNP2 = NUMNP 598 GOTO 3436 599 ENDIF 600 3435 CONTINUE 601 ENDIF 602 IF (IARCH.EQ.0) THEN 603 IF (IVIEUX.EQ.0) THEN 604 READ(IFAC,END=3232,ERR=3232) NE,(NODE(J),J=1,NDS) 605 & ,NPR,(NFPR(J),PRES(J),J=1,NPR) 606 GOTO 3233 607 3232 NPR = 0 608 IVIEUX = 1 609 3233 NUMEL = NUMEL+1 610 ELSE 611 READ(IFAC) NE,(NODE(J),J=1,NDS) 612 NUMEL = NUMEL+1 613 ENDIF 614 ELSE 615 CALL litrecbin(IFAC,IBID,LLL,0) 616cc NE = IBID(1) 617 IF (LLL/4.GT.NDS+1) THEN 618 NPR = IBID(NDS+2) 619 IF (NPR.GT.0) THEN 620 DO J=1,NPR 621 NFPR(J) = IBID(NDS+1+2*J) 622 PRES(J) = XBID(NDS+2*(J+1)) 623 ENDDO 624 ENDIF 625 ELSE 626 NPR = 0 627 ENDIF 628 NUMEL = NUMEL+1 629 ENDIF 630 IF (NPR.GT.0) THEN 631 DO J=1,NPR 632 IF (PRES(J).EQ.9999.) THEN 633 IR = 5 634 ELSEIF(PRES(J).EQ.-9999.) THEN 635 IR = 1 636 ELSE 637 IR = MOD(NINT(PRES(J)),16) 638 IF (IR.LT.0) IR = IR+16 639 ENDIF 640 IREF(NODE(IP8(1,NFPR(J)))) = IR 641 IREF(NODE(IP8(2,NFPR(J)))) = IR 642 ENDDO 643 ENDIF 644 CALL ECRNOD(NODE,NODEL,NUMEL,NDS) 645 DO I=1,NDS 646 IF (ISD(NODE(I)).EQ.0) THEN 647 ISD(NODE(I)) = IGROUP 648 ELSEIF(ISD(NODE(I)).NE.IGROUP) THEN 649 ISD(NODE(I)) = -1 650 ENDIF 651 ENDDO 652 NUMNP2 = NUMNP2+1 653 X(NUMNP2) = 654 & .125*( X(NODE(1))+X(NODE(2))+X(NODE(3))+X(NODE(4)) 655 & + X(NODE(5))+X(NODE(6))+X(NODE(7))+X(NODE(8)) ) 656 Y(NUMNP2) = 657 & .125*( Y(NODE(1))+Y(NODE(2))+Y(NODE(3))+Y(NODE(4)) 658 & + Y(NODE(5))+Y(NODE(6))+Y(NODE(7))+Y(NODE(8)) ) 659 Z(NUMNP2) = 660 & .125*( Z(NODE(1))+Z(NODE(2))+Z(NODE(3))+Z(NODE(4)) 661 & + Z(NODE(5))+Z(NODE(6))+Z(NODE(7))+Z(NODE(8)) ) 662 DEPX(NUMNP2) = .125*( DEPX(NODE(1))+DEPX(NODE(2)) 663 & + DEPX(NODE(3))+DEPX(NODE(4)) 664 & + DEPX(NODE(5))+DEPX(NODE(6)) 665 & + DEPX(NODE(7))+DEPX(NODE(8)) ) 666 DEPY(NUMNP2) = .125*( DEPY(NODE(1))+DEPY(NODE(2)) 667 & + DEPY(NODE(3))+DEPY(NODE(4)) 668 & + DEPY(NODE(5))+DEPY(NODE(6)) 669 & + DEPY(NODE(7))+DEPY(NODE(8)) ) 670 DEPZ(NUMNP2) = .125*( DEPZ(NODE(1))+DEPZ(NODE(2)) 671 & + DEPZ(NODE(3))+DEPZ(NODE(4)) 672 & + DEPZ(NODE(5))+DEPZ(NODE(6)) 673 & + DEPZ(NODE(7))+DEPZ(NODE(8)) ) 674 IF (IFACEXT.GT.0) THEN 675 IF (NB.EQ.NUMEL) THEN 676 IBORLU = 0 677 DO K=1,6 678 IFLAG = 0 679 DO KK=1,NFA 680 IF (IFAFA(KK).EQ.K) IFLAG=KK 681 ENDDO 682 IF (IFLAG.NE.0) THEN 683 IOPT = 0 684 CALL ORDONN(NODE(IP8(1,K)),NODE(IP8(2,K)) 685 & ,NODE(IP8(3,K)),NODE(IP8(4,K)),M1,M2,M3,M4) 686 CALL TRIAGE(M1,M2,M3,M4,ICPT,ICPT0,NFC,IOPT,NUMEL,K 687 & ,NUMNP2) 688 ISD2(ICPT) = IGROUP 689 ENDIF 690 ENDDO 691 ENDIF 692 ELSE 693 DO K=1,6 694 CALL ORDONN(NODE(IP8(1,K)),NODE(IP8(2,K)) 695 & ,NODE(IP8(3,K)),NODE(IP8(4,K)),M1,M2,M3,M4) 696 CALL TRIAGE(M1,M2,M3,M4,ICPT,ICPT0,NFC,IOPT,NUMEL,K 697 & ,NUMNP2) 698 ISD2(ICPT) = IGROUP 699 ENDDO 700 IF (NFC.GT.NFCM) CALL COMPAC(ICPT,ICPT0,NFC) 701 ENDIF 702 DO K=1,8 703 INDEX(NODE(K)) = INDEX(NODE(K))+1 704 ENDDO 705 CALL TETRA(NODE) 706 ENDDO 707C 708C 9 noeuds.... 709C 710 ELSEIF(NDS.EQ.9) THEN 711C 712C c'est comme ca qu'on dit qu'on fait du Q2 (ou P2). idiot 713C 714 MODTET = 40 715 IOPREF = 1 716 NBFAC = 4 717 ICOIN = 1 718 IF (I2D.EQ.0) THEN 719 ELEMENTS = 'Quadrangles 3D Q2' 720 ELSE 721 ELEMENTS = 'Quadrangles 2D Q2' 722 ENDIF 723 LELEM = 17 724 DO N=1,NELEM 725 IF (MOD(N,NMODP).EQ.0.AND.ISTDOUT.EQ.0) 726 & CALL MYFLUSH(100.*REAL(N)/REAL(NELEM)) 727 IF (IARCH.EQ.0) THEN 728 READ(IFAC) NE,(NODE(J),J=1,NDS) 729 ELSE 730 CALL litrecbin(IFAC,IBID,LLL,0) 731cc NE = IBID(1) 732 ENDIF 733 IMAT = IREF3(NODE(9)) 734 NUMEL = NUMEL+1 735 CALL ECRNOD(NODE,NODEL,NUMEL,NDS) 736 DO I=1,4 737 ICPT = ICPT+1 738 DO J=1,4 739 NN = NODE(IP9(J,I)) 740 INDEX(NN) = INDEX(NN)+1 741 IF (INDEX(NN).LE.NVMX) NVOI(INDEX(NN),NN) = ICPT 742 IF (ISD(NN).EQ.0) THEN 743 ISD(NN) = IGROUP 744 ELSEIF(ISD(NN).NE.IGROUP) THEN 745 ISD(NN) = -1 746 ENDIF 747 ICLAS(J,ICPT) = NN 748 ENDDO 749 ICLAS(5,ICPT) = NUMEL 750 ICLAS(6,ICPT) = 0 751 NUMNP2 = NUMNP2+1 752 X(NUMNP2) = X(ICLAS(1,ICPT)) 753 Y(NUMNP2) = Y(ICLAS(1,ICPT)) 754 Z(NUMNP2) = Z(ICLAS(1,ICPT)) 755 DEPX(NUMNP2) = DEPX(ICLAS(1,ICPT)) 756 DEPY(NUMNP2) = DEPY(ICLAS(1,ICPT)) 757 DEPZ(NUMNP2) = DEPZ(ICLAS(1,ICPT)) 758 ICLAS(7,ICPT) = NUMNP2 759 ISD2(ICPT) = IGROUP 760 ISD3(ICPT) = IMAT 761 ENDDO 762 ENDDO 763 ELSE 764 IF (ILANG.EQ.0) THEN 765 PRINT*,'***** El�ments �',NDS,' noeuds de type inconnu *****' 766 ELSE 767 PRINT*,'*****',NDS,'-nodes elements are unsupported *****' 768 ENDIF 769 STOP 770 ENDIF 771 IF (ISTDOUT.EQ.0) CALL ENDFLUSH 772 NDSMAX = MAX(NDSMAX,NDS) 773 IF (ICOURB.NE.-5) GOTO 1000 774C 775C Recopie les coordonnees des points des faces selectionnees 776C 777 9999 IF (ICOURB.NE.-5) CLOSE(IFAC) 778 IF (IFACEXT.GT.0) CLOSE(IFACEXT) 779 ZTRI = NDSMAX.EQ.2.OR.NDSMAX.EQ.3.OR.NDSMAX.EQ.6 780 & .OR.(NDSMAX.EQ.4.AND.I2D.EQ.0.AND.NTYP.NE.3) 781 IF (ZTRI) THEN 782 NDSB = 3 783 ELSE 784 NDSB = 4 785 ENDIF 786 IF (INOPO.NE.0) THEN 787 IF (ISAUVEMESH.EQ.1) THEN 788 CALL EXEC('/bin/mv -vi '//NOM_FICH(1:LONG)//'* .') 789 ENDIF 790 CALL EXEC('/bin/rm -f '//NOM_FICH(1:LONG)//'*') 791 ENDIF 792 NFACE = 0 793 DO 30 N=1,ICPT 794 IF (ICLAS(1,N).NE.0) THEN 795 NFACE = NFACE + 1 796 ISD2(NFACE) = ISD2(N) 797 ISD3(NFACE) = ISD3(N) 798 NUMSD = MAX(NUMSD,ISD2(NFACE)) 799 NNUMFA(NFACE) = ICLAS(5,N) 800 NRFAC(NFACE) = ICLAS(6,N) 801 DO I=1,NDSB 802 NFAC(I,NFACE) = ICLAS(I,N) 803 XF(I,NFACE) = X(ICLAS(I,N)) 804 YF(I,NFACE) = Y(ICLAS(I,N)) 805 ZF(I,NFACE) = Z(ICLAS(I,N)) 806 ENDDO 807 NFAC(NDSB+1,NFACE) = ICLAS(7,N) 808 XF(NDSB+1,NFACE) = X(ICLAS(7,N)) 809 YF(NDSB+1,NFACE) = Y(ICLAS(7,N)) 810 ZF(NDSB+1,NFACE) = Z(ICLAS(7,N)) 811 IF (ZTRI) THEN 812 XMA = MAX(X(ICLAS(1,N)),X(ICLAS(2,N)) 813 & ,X(ICLAS(3,N))) 814 XMI = MIN(X(ICLAS(1,N)),X(ICLAS(2,N)) 815 & ,X(ICLAS(3,N))) 816 ISMIN = MIN( ISD(ICLAS(1,N)),ISD(ICLAS(2,N)) 817 & ,ISD(ICLAS(3,N)) ) 818 ISMAX = MAX( ISD(ICLAS(1,N)),ISD(ICLAS(2,N)) 819 & ,ISD(ICLAS(3,N)) ) 820 ELSE 821 XMA = MAX(X(ICLAS(1,N)),X(ICLAS(2,N)) 822 & ,X(ICLAS(3,N)),X(ICLAS(4,N))) 823 XMI = MIN(X(ICLAS(1,N)),X(ICLAS(2,N)) 824 & ,X(ICLAS(3,N)),X(ICLAS(4,N))) 825 ISMIN = MIN( ISD(ICLAS(1,N)),ISD(ICLAS(2,N)) 826 & ,ISD(ICLAS(3,N)),ISD(ICLAS(4,N)) ) 827 ISMAX = MAX( ISD(ICLAS(1,N)),ISD(ICLAS(2,N)) 828 & ,ISD(ICLAS(3,N)),ISD(ICLAS(4,N)) ) 829 ENDIF 830 IF (EGAL(XMI,XMA)) THEN 831 IF (EGAL(XMI,-XMED).AND.ISYM.EQ.4) THEN 832 IFPLAN(NFACE) = 3 833 ELSE 834 IFPLAN(NFACE) = 0 835 ENDIF 836 ELSE 837 IF (ZTRI) THEN 838 YMA = MAX(Y(ICLAS(1,N)),Y(ICLAS(2,N)) 839 & ,Y(ICLAS(3,N))) 840 YMI = MIN(Y(ICLAS(1,N)),Y(ICLAS(2,N)) 841 & ,Y(ICLAS(3,N))) 842 ELSE 843 YMA = MAX(Y(ICLAS(1,N)),Y(ICLAS(2,N)) 844 & ,Y(ICLAS(3,N)),Y(ICLAS(4,N))) 845 YMI = MIN(Y(ICLAS(1,N)),Y(ICLAS(2,N)) 846 & ,Y(ICLAS(3,N)),Y(ICLAS(4,N))) 847 ENDIF 848 IF (EGAL(YMI,YMA)) THEN 849 IF (EGAL(YMI,-YMED)) THEN 850 IFPLAN(NFACE) = 2 851 ELSE 852 IFPLAN(NFACE) = 0 853 ENDIF 854 ENDIF 855 ENDIF 856 IF (ISMIN.NE.ISMAX) THEN 857 ITOTO = ISD2(NFACE) 858 IP = 1000 859 IF (ISD(ICLAS(1,N)).LT.0.AND.ISD(ICLAS(2,N)).LT.0) THEN 860 ITOTO = ITOTO + IP 861 IP = IP*10 862 ENDIF 863 IF (ISD(ICLAS(2,N)).LT.0.AND.ISD(ICLAS(3,N)).LT.0) THEN 864 ITOTO = ITOTO + IP*2 865 IP = IP*10 866 ENDIF 867 IF (ZTRI) THEN 868 IF (ISD(ICLAS(3,N)).LT.0.AND.ISD(ICLAS(1,N)).LT.0) THEN 869 ITOTO = ITOTO + IP*3 870 ENDIF 871 ELSE 872 IF (ISD(ICLAS(3,N)).LT.0.AND.ISD(ICLAS(4,N)).LT.0) THEN 873 ITOTO = ITOTO + IP*3 874 IP = IP*10 875 ENDIF 876 IF (ISD(ICLAS(4,N)).LT.0.AND.ISD(ICLAS(1,N)).LT.0) THEN 877 ITOTO = ITOTO + IP*4 878 ENDIF 879 ENDIF 880 ISD2(NFACE) = ITOTO 881 ENDIF 882 ENDIF 883 30 CONTINUE 884 IGROUP = NUMSD 885C 886C Pourquoi c'etait commente ? 6/98 887C 888 IF (I2D.NE.0) NDS = NDSB 889C 890 IF (NDS.EQ.3.OR.NDS.EQ.2) IPREFC = -1 891 IF (IPREFC.NE.-1) IORIENT = 1 892 IF (NDS.EQ.2) IFC = 1 893 IF (ISTDOUT.EQ.0) THEN 894 IF (ILANG.EQ.0) THEN 895 IF (IVIEUX.NE.0) 896 & PRINT*,'*** Attention !! Vieux format de maillage.........' 897 PRINT*,'Fin de lecture. Recherche de la fronti�re.........' 898 ELSE 899 IF (IVIEUX.NE.0) 900 & PRINT*,'*** Warning !! Old-style mesh file................' 901 PRINT*,'End of reading. Computing the boundary............' 902 ENDIF 903 ENDIF 904 CALL MINDIAM(NFACE) 905 CALL FROETC(NTYP) 906 NUMBIS = NUMNP2 907 IF (ZTRI) NDS=3 908C 909C Determination des elements du bord 910C 911 DO I=1,NEL 912 IELBOR(I) = 0 913 ENDDO 914 DO I=1,NFACE 915 IF (IELBOR(NNUMFA(I)).EQ.0) THEN 916 IELBOR(NNUMFA(I)) = I 917 ELSE 918 IELBOR(NNUMFA(I)) = -I 919 ENDIF 920 ENDDO 921 NNN = 0 922 DO I=1,NEL 923 IF (IELBOR(I).NE.0) NNN = NNN+1 924 ENDDO 925 IF (I2D.EQ.0.AND.IPS2D.EQ.0) THEN 926 IF (ISTDOUT.EQ.0) THEN 927 IF (ILANG.EQ.0) THEN 928 PRINT*,NNN,' �l�ments de bord sur' 929 & ,NEL,' (',REAL(NINT(1000.*REAL(NNN)/REAL(NEL)))*.1,' %)' 930 ELSE 931 PRINT*,NNN,' boundary elements among' 932 & ,NEL,' (',REAL(NINT(1000.*REAL(NNN)/REAL(NEL)))*.1,' %)' 933 ENDIF 934 ENDIF 935 ENDIF 936C 937C Ecriture du fichier .facext 938C (8 noeuds et 4 noeuds tetra seulement pour l'instant) 939C 940 IF (IFACEXT.EQ.-1.AND.IELIMI.EQ.0) THEN 941 IF (CBID(1:5).NE.'/tmp/') THEN 942 CALL PREMIER_LIBRE(IFACEXT) 943 OPEN(IFACEXT,FILE=CBID(1:LBID),IOSTAT=IERR) 944 IF (IERR.EQ.0) THEN 945 DO I=1,NEL 946 ITAB(I) = 0 947 ENDDO 948 DO I=1,NFACE 949 ITAB(NNUMFA(I)) = ITAB(NNUMFA(I))+1 950 ENDDO 951 NN = 0 952 DO I=1,NEL 953 IF (ITAB(I).GT.0) THEN 954 DO J=1,ITAB(I) 955 NN = NN+1 956 IFAFA(J) = NRFAC(NN) 957 ENDDO 958 WRITE(IFACEXT,*) I,ITAB(I),(IFAFA(K),K=1,ITAB(I)) 959 ENDIF 960 ENDDO 961 CLOSE(IFACEXT) 962 ENDIF 963 ENDIF 964 ENDIF 965C 966 END 967C======================================================================= 968 SUBROUTINE ECRNOD(NODE,NODEL,NUMEL,NDS) 969C 970 INTEGER NODE(*),NODEL(*) 971C 972 II = (NUMEL-1)*NDS 973 DO I=1,NDS 974 NODEL(II+I) = NODE(I) 975 ENDDO 976C 977C calcul des criteres de qualite des EF 978C 979 CALL QUALITY(NODE,NUMEL,NDS) 980 END 981C======================================================================= 982 SUBROUTINE QUALITY(NODE,N,ND) 983 INCLUDE 'com_coor.f' 984 INCLUDE 'com_faces.f' 985C 986 DIMENSION ANG(4) 987 INTEGER NODE(*),ISUC3(3),ISUC4(4),IPRE4(4),IP8(4,6) 988 DATA ISUC3 / 2,3,1 / 989 DATA ISUC4 / 2,3,4,1 / 990 DATA IPRE4 / 4,1,2,3 / 991 DATA IP8 / 5,8,4,1, 992 & 2,3,7,6, 993 & 5,1,2,6, 994 & 4,8,7,3, 995 & 1,4,3,2, 996 & 8,5,6,7 / 997C 998 RMAX = 0. 999 IF (ND.EQ.3.OR.ND.EQ.6) THEN 1000C 1001C Triangles 1002C 1003 XC = (X(NODE(1))+X(NODE(2))+X(NODE(3)))/3. 1004 YC = (Y(NODE(1))+Y(NODE(2))+Y(NODE(3)))/3. 1005 ZC = (Z(NODE(1))+Z(NODE(2))+Z(NODE(3)))/3. 1006 DO J=1,3 1007 K = ISUC3(J) 1008 L = ISUC3(K) 1009 ANG(J) = CALCANG(X,Y,Z,NODE,J,K,L) 1010 RR = SQRT((X(NODE(J))-XC)**2 1011 & +(Y(NODE(J))-YC)**2 1012 & +(Z(NODE(J))-ZC)**2) 1013 RMAX = AMAX1(RMAX,RR) 1014 ENDDO 1015 QUALITE(N,1) = AMIN1(ANG(1),ANG(2),ANG(3)) 1016 & /AMAX1(ANG(1),ANG(2),ANG(3)) 1017 ELSEIF(ND.EQ.4.OR.ND.EQ.9) THEN 1018C 1019C quadrangles (4,9 noeuds) (a revoir pour les tetra) 1020C 1021 XC = (X(NODE(1))+X(NODE(2))+X(NODE(3))+X(NODE(4)))*0.25 1022 YC = (Y(NODE(1))+Y(NODE(2))+Y(NODE(3))+Y(NODE(4)))*0.25 1023 ZC = (Z(NODE(1))+Z(NODE(2))+Z(NODE(3))+Z(NODE(4)))*0.25 1024 DO J=1,4 1025 K = ISUC4(J) 1026 L = IPRE4(J) 1027 ANG(J) = CALCANG(X,Y,Z,NODE,J,K,L) 1028 RR = SQRT((X(NODE(J))-XC)**2 1029 & +(Y(NODE(J))-YC)**2 1030 & +(Z(NODE(J))-ZC)**2) 1031 RMAX = AMAX1(RMAX,RR) 1032 ENDDO 1033 QUALITE(N,1) = AMIN1(ANG(1),ANG(2),ANG(3),ANG(4)) 1034 & /AMAX1(ANG(1),ANG(2),ANG(3),ANG(4)) 1035 ELSE 1036C 1037C Hexa 8, 27 noeuds 1038C 1039 ANGMIN = 5. 1040 ANGMAX = 0. 1041 XC = (X(NODE(1))+X(NODE(2))+X(NODE(3))+X(NODE(4)) 1042 & +X(NODE(5))+X(NODE(6))+X(NODE(7))+X(NODE(8)))*0.125 1043 YC = (Y(NODE(1))+Y(NODE(2))+Y(NODE(3))+Y(NODE(4)) 1044 & +Y(NODE(5))+Y(NODE(6))+Y(NODE(7))+Y(NODE(8)))*0.125 1045 ZC = (Z(NODE(1))+Z(NODE(2))+Z(NODE(3))+Z(NODE(4)) 1046 & +Z(NODE(5))+Z(NODE(6))+Z(NODE(7))+Z(NODE(8)))*0.125 1047 DO I=1,6 1048 DO II=1,4 1049 J = IP8(II,I) 1050 K = IP8(ISUC4(II),I) 1051 L = IP8(IPRE4(II),I) 1052 AA = CALCANG(X,Y,Z,NODE,J,K,L) 1053 ANGMAX = AMAX1(ANGMAX,AA) 1054 ANGMIN = AMIN1(ANGMIN,AA) 1055 ENDDO 1056 ENDDO 1057 QUALITE(N,1) = ANGMIN/ANGMAX 1058 DO J=1,8 1059 RR = SQRT((X(NODE(J))-XC)**2 1060 & +(Y(NODE(J))-YC)**2 1061 & +(Z(NODE(J))-ZC)**2) 1062 RMAX = AMAX1(RMAX,RR) 1063 ENDDO 1064 ENDIF 1065C 1066C Pas fini... 1067C 1068 QUALITE(N,2) = RMAX 1069C 1070 END 1071C======================================================================= 1072 REAL*4 FUNCTION CALCANG(X,Y,Z,NODE,J,K,L) 1073 INTEGER NODE(*) 1074 DIMENSION X(*),Y(*),Z(*) 1075CC DATA CQVSPI / 57.29577951 / 1076 DATA PI / 3.14159265358979 / 1077C 1078 JJ = NODE(J) 1079 KK = NODE(K) 1080 LL = NODE(L) 1081 U1 = X(KK)-X(JJ) 1082 V1 = Y(KK)-Y(JJ) 1083 W1 = Z(KK)-Z(JJ) 1084 U2 = X(LL)-X(JJ) 1085 V2 = Y(LL)-Y(JJ) 1086 W2 = Z(LL)-Z(JJ) 1087 XN1 = SQRT(U1**2+V1**2+W1**2) 1088 XN2 = SQRT(U2**2+V2**2+W2**2) 1089 IF (XN1*XN2.NE.0.) THEN 1090 ANG = ASIN(MIN(1.,SQRT((U1*V2-U2*V1)**2 1091 & +(V1*W2-V2*W1)**2 1092 & +(W1*U2-W2*U1)**2)/(XN1*XN2))) 1093 SCAL = U1*U2 + V1*V2 + W1*W2 1094 IF (SCAL.GE.0.) THEN 1095 CALCANG = ANG 1096 ELSE 1097 CALCANG = PI-ANG 1098 ENDIF 1099 ELSE 1100 CALCANG = -1. 1101 ENDIF 1102 END 1103