1 : C VOLUME SOURCE FANDEUR 10/08/30 21:15:51 6732 2 : C MODIF : O.STAB / 25.03.97 / APPEL A VOLOS QUI 3 : C AUTORISE UN RACCORD DE 2 GRILLES NON IDENTIQUE 4 : C FABRICATION DE CUBES ET PRISMES PAR TRANSLATION ET ROTATION ET 5 : C ENTRE SURFACES OPPOSEES 6 : C MODIFICATION AOUT 1984 : MAILLAGE AUTOMATIQUE A L'INTERIEUR D'UNE 7 : C SURFACE ENVELOPPE 8 : C DECEMBRE 1984 VERIFICATION QUE LES TOPOLOGIES DU HAUT ET DU BAS 9 : C SONT SIMILAIRES (MEMES ICPR) 10 : C JANVIER 1985 NOMBRE DE COUCHES IMPOSE SI INBR NEGATIF 11 : C NOVEMBRE 1985 OPTION GENERATRICE (APPEL VOLUMG) 12 : c 12/97 KICH modif en liaison avec evolution de PROPER 13 : 14 : SUBROUTINE VOLUME 15 : IMPLICIT INTEGER(I-N) 16 : IMPLICIT REAL*8 (A-H,O-Z) 17 : -INC SMELEME 18 : -INC CCOPTIO 19 : -INC SMCOORD 20 : -INC CCGEOME 21 : -INC CCREEL 22 : *- 23 : -INC SMLREEL 24 : logical ltelq 25 : SEGMENT TABPAR(NCOUCH) 26 : SEGMENT ICPR(NBNNEL,NBELEC) 27 : CHARACTER*4 MCLE(6) 28 : DATA MCLE /'TRAN','ROTA','DINI','DFIN','GENE','PROG'/ 29 : 30 : IDIMP1 = IDIM + 1 31 : c 32 : MLREEL=0 33 : IMPOI=0 34 : IMPOF=0 35 : c 36 : IF (ILCOUR.LT.14.OR.ILCOUR.GT.17) then 36 : IF (ILCOUR.LT.23.OR.ILCOUR.GT.26) CALL ERREUR(16) 37 : ENDIF 38 : IF (IERR.NE.0) RETURN 39 : 40 : INBR=0 41 : ICLE=3 41 : CALL LIRENT(INBR,0,IREINB) 41 : 80 CALL LIRMOT(MCLE,6,JCLE,0) 42 : IF (JCLE.EQ.0) GOTO 87 43 : GOTO (81,82,83,84,79,78),JCLE 44 : C ... mot-clé 'TRAN' ... 45 : 81 CONTINUE 46 : C ... mot-clé 'ROTA' ... 47 : 82 CONTINUE 48 : C POUR LE CAS OU LE DECOUPAGE EST DERRIERE TRAN 48 : IF (INBR.EQ.0) CALL LIRENT(INBR,0,IREINB) 49 : IF (ICLE.NE.3) GOTO 86 50 : ICLE=JCLE 51 : GOTO 80 52 : C ... mot-clé 'DINI' ... 53 : 83 IF(IMPOI.EQ.1) GOTO 86 54 : IMPOI=1 54 : CALL LIRREE(XXX,1,IRETOU) 55 : DEN1D=XXX 56 : IF (DEN1D.LE.0.) THEN 56 : CALL ERREUR(17) 57 : RETURN 58 : ENDIF 59 : IF (IERR.NE.0) RETURN 60 : GOTO 80 61 : C ... mot-clé 'DFIN' ... 62 : 84 IF (IMPOF.EQ.1) GOTO 86 63 : IMPOF=1 63 : CALL LIRREE(XXX,1,IRETOU) 64 : DEN2D=XXX 65 : IF (DEN2D.LE.0.) THEN 65 : CALL ERREUR(17) 66 : RETURN 67 : ENDIF 68 : IF (IERR.NE.0) RETURN 69 : GOTO 80 70 : C ... mot-clé 'GENE' ... 71 : 79 CONTINUE 71 : CALL VOLUMG 72 : RETURN 73 : C ... mot-clé 'PROG' ... 74 : 78 CONTINUE 74 : CALL LIROBJ('LISTREEL',MLREEL,1,IRETOU) 75 : IF (IERR.NE.0) RETURN 76 : GOTO 80 77 : 86 CONTINUE 77 : CALL REFUS 78 : 79 : C ... Fin de lecture des mots clés ... 80 : 87 CONTINUE 80 : CALL LIROBJ('MAILLAGE',IPT1,1,IRETOU) 81 : C ... Si pas d'option TRAN ni ROTA on veut lire un deuxième maillage ... 81 : IF (ICLE.EQ.3) CALL LIROBJ('MAILLAGE',IPT2,0,IRETOU) 82 : C ... S'il n'y en a pas, on va remplir l'enveloppe .(ou épaisseur).. 83 : IF(IRETOU.EQ.0) GOTO 4400 84 : C ================================== 85 : C ------- DEBUT DE MODIF - O.STAB 05.12.96 -------- 86 : C ================================== 87 : C 88 : c WRITE(IOIMP,*) ' ICLE = ',ICLE 89 : IF((ICLE.NE.1).AND.(ICLE.NE.2).AND.(ICLE.NE.5))THEN 89 : CALL LIROBJ('POINT ',IPO1,0,IRETOU) 90 : IF( IRETOU.EQ.0)GOTO 871 90 : CALL LIROBJ('POINT ',IPO2,0,IRETOU) 91 : IF( IRETOU.EQ.0)GOTO 871 92 : C 93 : C calcul des densités moyennes si pas données 94 : C 95 : DEN1D=1. 96 : IF(IMPOI.EQ.0.AND.INBR.EQ.0) THEN 97 : MELEME = IPT1 98 : SEGACT MELEME 99 : NBNN=NUM(/1) 100 : NBELEM=NUM(/2) 101 : NPR=0 102 : DEN1=0.D0 103 : DO 710 I=1,NBNN 104 : DO 710 J=1,NBELEM 105 : IR=NUM(I,J) 106 : IF (IR.EQ.0) GOTO 710 107 : IREF=(IR-1)*IDIMP1 108 : NPR=NPR+1 109 : DEN1=DEN1+XCOOR(IREF+4) 110 : 710 CONTINUE 111 : DEN1D=DEN1/NPR 112 : SEGDES MELEME 113 : ENDIF 114 : DEN2D=1. 115 : IF(IMPOF.EQ.0.AND.INBR.EQ.0) THEN 116 : MELEME = IPT2 117 : SEGACT MELEME 118 : NBNN=NUM(/1) 119 : NBELEM=NUM(/2) 120 : NPR=0 121 : DEN2=0.D0 122 : DO 711 I=1,NBNN 123 : DO 711 J=1,NBELEM 124 : IR=NUM(I,J) 125 : IF (IR.EQ.0) GOTO 711 126 : IREF=(IR-1)*IDIMP1 127 : NPR=NPR+1 128 : DEN2=DEN2+XCOOR(IREF+4) 129 : 711 CONTINUE 130 : DEN2D=DEN2/NPR 131 : SEGDES MELEME 132 : ENDIF 132 : CALL VOLOS(IPT1,IPT2,IPO1,IPO2,DEN1D,DEN2D,INBR) 133 : RETURN 134 : ENDIF 135 : C ================================== 136 : C ------- FIN DE MODIF - O.STAB 05.12.96 -------- 137 : C ================================== 138 : C 139 : 871 CONTINUE 140 : 141 : C ... Début du traitement ... 142 : ISVOL1=0 143 : ISVOL2=0 144 : SEGACT IPT1 145 : C SI IPT1 VOLUME IL FAUT EN EXTRAIRE LA FACE 1 146 : C ... dans la pratique le maillage initial soit n'a pas de 147 : C sous-maillages, soit il en a 2 (triangles et quadrilatères), 148 : C sinon la programmation ci-dessous poserait des problèmes : 149 : C en cas de IPT1.LISOUS(/1) > 2 un saut immédiat vers 3101 150 : C provoquerait SEGDES IPT3,IPT4 qui ne sont pas encore initialisés ... 151 : 3100 IF (IPT1.LISOUS(/1).EQ.0) GOTO 1000 152 : IF (IPT1.LISOUS(/1).NE.2) GOTO 3101 153 : IDEUX=2 154 : IPT3=IPT1.LISOUS(1) 155 : IPT4=IPT1.LISOUS(2) 156 : SEGACT IPT3,IPT4 157 : IP=IPT3.ITYPEL*IPT4.ITYPEL 158 : c ... TRI3*QUA4 TRI6*QUA8 ... 159 : IF (IP.NE.32.AND.IP.NE.60) GOTO 3101 160 : IS=IPT3.ITYPEL+IPT4.ITYPEL 161 : c ... TRI3+QUA4 TRI6+QUA8 ... 162 : IF (IS.NE.12.AND.IS.NE.16) GOTO 3101 163 : c ... ici on a deux maillages : un composé de triangles, l'autre de quadrilatères ... 164 : INCR=1 165 : IF (IS.EQ.16) INCR=2 166 : c ... NBNNEL = nombre de noeuds / élément du maillage total ... 167 : NBNNEL=4*INCR 168 : C EN FAIT ON CREE UN SEGMENT QUI CONTIENT LES CUBES ET LES TRIANGLES 169 : C 0 DANS LA DERNIERE POSITION DU TRIANGLE 170 : NBSOUS=0 171 : NBREF=0 172 : NBNN=NBNNEL 173 : NBELE3=IPT3.NUM(/2) 174 : IF (IPT3.ITYPEL.LE.6) NBTRI=NBELE3 175 : IF (IPT3.ITYPEL.GE.8) NBQUA=NBELE3 176 : NBELE4=IPT4.NUM(/2) 177 : IF (IPT4.ITYPEL.LE.6) NBTRI=NBELE4 178 : IF (IPT4.ITYPEL.GE.8) NBQUA=NBELE4 179 : NBELEM=NBELE3+NBELE4 180 : SEGINI MELEME 181 : C*C ... Initialisation du nouveau maillage à 0 (est-ce nécessaire ?) ... 182 : C* DO 1100 I=1,NBNN 183 : C* DO 1100 J=1,NBELEM 184 : C* NUM(I,J)=0 185 : C* 1100 CONTINUE 186 : C ... On transvase IPT3 dans le nouveau maillage ... 187 : DO 1101 J=1,NBELE3 188 : ICOLOR(J)=IPT3.ICOLOR(J) 189 : DO 1101 I=1,IPT3.NUM(/1) 190 : NUM(I,J)=IPT3.NUM(I,J) 191 : 1101 CONTINUE 192 : C ... On transvase IPT4 dans le nouveau maillage ... 193 : DO 1102 J=1,NBELE4 194 : K=J+NBELE3 195 : ICOLOR(K)=IPT4.ICOLOR(J) 196 : DO 1102 I=1,IPT4.NUM(/1) 197 : NUM(I,K)=IPT4.NUM(I,J) 198 : 1102 CONTINUE 199 : GOTO 1001 200 : 201 : C RECHERCHE DE LA PREMIERE FACE DE IPT1 202 : 3101 SEGDES IPT3,IPT4 202 : 3102 IF (IPT1.LISREF(/1).LT.2) CALL ERREUR(16) 203 : IF (IERR.NE.0) RETURN 204 : ISVOL1=IPT1 205 : IAUX=IPT1.LISREF(2) 206 : SEGDES IPT1 207 : IPT1=IAUX 208 : SEGACT IPT1 209 : GOTO 3100 210 : 211 : c ... On vient ici si le maillage est simple (homogène) ... 212 : 1000 CONTINUE 213 : IDEUX=1 214 : NBNNEL=IPT1.NUM(/1) 215 : IF (IPT1.ITYPEL.NE.8.AND.IPT1.ITYPEL.NE.10.AND.IPT1.ITYPEL.NE.4 216 : #.AND.IPT1.ITYPEL.NE.6) GOTO 3102 217 : INCR=1 218 : IF (KDEGRE(IPT1.ITYPEL).EQ.3) INCR=2 219 : MELEME=IPT1 220 : 221 : c ... ici MELEME est le maillage contanant tous les éléments de la surface initiale ... 222 : 1001 SEGACT MCOORD 223 : c ... SI c'est 'ROTA' ... 224 : IF (ICLE.EQ.2) GOTO 1 225 : c ... si ce n'est ni ROTA ni TRAN ... 226 : IF (ICLE.EQ.3) GOTO 2 227 : c ... alors c'est TRAN ... 227 : CALL LIROBJ('POINT ',IVEC,1,IRETOU) 228 : IF (IERR.NE.0) RETURN 229 : IREF=(IVEC-1)*IDIMP1 230 : XTRAN=XCOOR(IREF+1) 231 : YTRAN=XCOOR(IREF+2) 232 : ZTRAN=XCOOR(IREF+3) 233 : DEN2=XCOOR(IREF+4) 234 : DLONG=SQRT(XTRAN*XTRAN+YTRAN*YTRAN+ZTRAN*ZTRAN) 235 : XDIS=XTRAN 236 : YDIS=YTRAN 237 : ZDIS=ZTRAN 238 : GOTO 2 239 : 240 : C ... cas 'ROTA' ... 241 : c ... lecture de l'angle ... 241 : 1 CALL LIRREE(XXX,0,IRETOU) 242 : ANGLI=XXX 243 : IF (IRETOU.EQ.1) GOTO 85 244 : IF (IREINB.NE.0) THEN 245 : ANGLI=INBR 246 : INBR=0 247 : ELSE 247 : CALL LIRREE(XXX,1,IRETOU) 248 : ANGLI=XXX 249 : IF (IERR.NE.0) RETURN 250 : ENDIF 251 : 85 CONTINUE 252 : ANGLE=ANGLI*XPI/180.D0 253 : c ... et des deux points définissant l'axe de rotation ... 253 : CALL LIROBJ('POINT ',IP1,1,IRETOU) 253 : CALL LIROBJ('POINT ',IP2,1,IRETOU) 254 : IF (IERR.NE.0) RETURN 255 : IREF1=(IP1-1)*IDIMP1 256 : IREF2=(IP2-1)*IDIMP1 257 : XPT1=XCOOR(IREF1+1) 258 : YPT1=XCOOR(IREF1+2) 259 : ZPT1=XCOOR(IREF1+3) 260 : XVEC=XCOOR(IREF2+1)-XCOOR(IREF1+1) 261 : YVEC=XCOOR(IREF2+2)-XCOOR(IREF1+2) 262 : ZVEC=XCOOR(IREF2+3)-XCOOR(IREF1+3) 263 : DEN2=(XCOOR(IREF2+4)+XCOOR(IREF1+4))*0.5 264 : RAY=SQRT(XVEC*XVEC+YVEC*YVEC+ZVEC*ZVEC) 265 : XVEC=XVEC/RAY 266 : YVEC=YVEC/RAY 267 : ZVEC=ZVEC/RAY 268 : IF (ANGLE.GE.0.) GOTO 2 269 : ANGLE=-ANGLE 270 : XVEC=-XVEC 271 : YVEC=-YVEC 272 : ZVEC=-ZVEC 273 : 274 : c ... calcul des moyennes des densités et des coordonnées ... 275 : 2 CONTINUE 276 : NBNN=NUM(/1) 277 : NBELEM=NUM(/2) 278 : NPR=0 279 : DEN1=0. 280 : XG=0. 281 : YG=0. 282 : ZG=0. 283 : DO 5 I=1,NBNN 284 : DO 5 J=1,NBELEM 285 : IR=NUM(I,J) 286 : IF (IR.EQ.0) GOTO 5 287 : IREF=(IR-1)*IDIMP1 288 : NPR=NPR+1 289 : DEN1=DEN1+XCOOR(IREF+4) 290 : XG=XG+XCOOR(IREF+1) 291 : YG=YG+XCOOR(IREF+2) 292 : ZG=ZG+XCOOR(IREF+3) 293 : 5 CONTINUE 294 : DEN1=DEN1/NPR 295 : XG=XG/NPR 296 : YG=YG/NPR 297 : ZG=ZG/NPR 298 : 299 : c ... cas 'TRAN' => GOTO 6 ... 300 : IF (ICLE.EQ.1) GOTO 6 301 : c ... cas 'ROTA' => GOTO 3 ... 302 : IF (ICLE.EQ.2) GOTO 3 303 : c ... cas du volume entre deux surfaces ... 304 : C COMPATIBILITE DU 2EME OBJET ET RECHERCHE DU CENTRE DE GRAVITE 305 : SEGACT IPT2 306 : 3150 IF (IPT2.LISOUS(/1).EQ.0) GOTO 1020 306 : IF (IDEUX.NE.2) CALL ERREUR(21) 307 : IF (IERR.NE.0) RETURN 308 : IF (IPT2.LISOUS(/1).NE.2) GOTO 3151 309 : IPT5=IPT2.LISOUS(1) 310 : IPT6=IPT2.LISOUS(2) 311 : SEGACT IPT5,IPT6 312 : IF (IPT5.ITYPEL.NE.IPT3.ITYPEL) GOTO 3151 313 : IF (IPT6.ITYPEL.NE.IPT4.ITYPEL) GOTO 3151 313 : IF (IPT5.NUM(/2).NE.IPT3.NUM(/2)) CALL ERREUR(21) 313 : IF (IPT6.NUM(/2).NE.IPT4.NUM(/2)) CALL ERREUR(21) 313 : IF (IPT5.NUM(/1).NE.IPT3.NUM(/1)) CALL ERREUR(21) 313 : IF (IPT6.NUM(/1).NE.IPT4.NUM(/1)) CALL ERREUR(21) 314 : IF (IERR.NE.0) RETURN 315 : GOTO 1021 316 : 3151 SEGDES IPT5,IPT6 316 : 3152 IF (IPT2.LISREF(/1).LT.2) CALL ERREUR(21) 317 : IF (IERR.NE.0) RETURN 318 : IAUX=IPT2.LISREF(1) 319 : ISVOL2=IPT2 320 : SEGDES IPT2 321 : IPT2=IAUX 322 : SEGACT IPT2 323 : GOTO 3150 324 : c ... les deux maillages doivent être simples ... 325 : 1020 IF (IPT1.ITYPEL.NE.IPT2.ITYPEL) GOTO 3152 326 : IF (IDEUX.NE.1) GOTO 3152 326 : IF (IPT2.NUM(/1).NE.IPT1.NUM(/1)) CALL ERREUR(21) 326 : IF (IPT2.NUM(/2).NE.IPT1.NUM(/2)) CALL ERREUR(21) 327 : 1021 CONTINUE 328 : c ... calcul des densités et coordonnées moyennes de la deuxième surface ... 329 : NPR=0 330 : XG2=0. 331 : YG2=0. 332 : ZG2=0. 333 : DEN2=0. 334 : IPT7=IPT2 335 : M1031=2 336 : IF (IDEUX.EQ.1) M1031=1 337 : DO 1031 M=1,M1031 338 : IF (M1031.NE.1) IPT7=IPT2.LISOUS(M) 339 : DO 4 I=1,IPT7.NUM(/1) 340 : DO 4 J=1,IPT7.NUM(/2) 341 : IREF=(IPT7.NUM(I,J)-1)*IDIMP1 342 : DEN2=DEN2+XCOOR(IREF+4) 343 : XG2=XG2+XCOOR(IREF+1) 344 : YG2=YG2+XCOOR(IREF+2) 345 : ZG2=ZG2+XCOOR(IREF+3) 346 : 4 CONTINUE 347 : NPR=NPR+IPT7.NUM(/1)*IPT7.NUM(/2) 348 : 1031 CONTINUE 349 : DEN2=DEN2/NPR 350 : XG2=XG2/NPR 351 : YG2=YG2/NPR 352 : ZG2=ZG2/NPR 353 : DLONG=SQRT((XG2-XG)**2+(YG2-YG)**2+(ZG2-ZG)**2) 354 : GOTO 6 355 : 356 : c ... cas 'ROTA' ... 357 : 3 CONTINUE 358 : XV1=XG-XPT1 359 : YV1=YG-YPT1 360 : ZV1=ZG-ZPT1 361 : PV1=XV1*XVEC+YV1*YVEC+ZV1*ZVEC 362 : XV1=XV1-PV1*XVEC 363 : YV1=YV1-PV1*YVEC 364 : ZV1=ZV1-PV1*ZVEC 365 : RAY=SQRT(XV1*XV1+YV1*YV1+ZV1*ZV1) 366 : XV1=XV1/RAY 367 : YV1=YV1/RAY 368 : ZV1=ZV1/RAY 369 : XV2=YVEC*ZV1-ZVEC*YV1 370 : YV2=ZVEC*XV1-XVEC*ZV1 371 : ZV2=XVEC*YV1-YVEC*XV1 372 : C RAYON MOYEN 373 : C ANGLE EN RADIANS D'OU LONGUEUR MOYENNE 374 : DLONG=RAY*ABS(ANGLE) 375 : 376 : c ... partie commune, recherche du nombre de couches à créér et des densités ... 377 : 6 CONTINUE 378 : IF (IMPOI.EQ.1) DEN1=DEN1D 379 : IF (IMPOF.EQ.1) DEN2=DEN2D 380 : C JE NE VOIS PAS DANS QUELS CAS CA INTERVIENT 380 : C CALL LIRENT(INBR,0,IRETOU) 381 : IPT3=MELEME 382 : DENI=DEN1 383 : DECA=DEN2-DEN1 384 : DEN1=DEN1/DLONG 385 : DEN2=DEN2/DLONG 386 : IF (MLREEL.NE.0)THEN 387 : SEGACT,MLREEL 388 : INBR=PROG(/1)-1 389 : SEGDES,MLREEL 390 : ENDIF 390 : CALL DECOUP(INBR,DEN1,DEN2,APROG,NCOUCH,DENI,DECA,DLONG) 391 : NX=NCOUCH-1 392 : 393 : IF (IIMPI.EQ.1) WRITE(IOIMP,9000) NCOUCH,APROG 394 : 9000 FORMAT(/' COUCHES ',I6,' RAISON ',G12.5) 395 : 396 : C ... Initialisation du nouveau maillage ... 397 : C ON FAIT TOUJOURS COMME SI IL N'Y AVAIT QU'UN TYPE D'ELEMENT 398 : NBSOUS=0 399 : C MODIF POUR CONSTRUIRE TOUJOURS LE POURTOUR 400 : NBREF=3 401 : IF (IPT1.LISREF(/1).NE.0) NBREF=3 402 : NBNN=2*NBNNEL+(INCR-1)*(NBNNEL/2) 403 : NBNNV=NBNN 404 : NBASE=NBELEM 405 : NBELEM=NBELEM*NCOUCH 406 : SEGINI IPT7 407 : IF (NBNNV.EQ.6 ) IPT7.ITYPEL=16 408 : IF (NBNNV.EQ.15) IPT7.ITYPEL=17 409 : IF (NBNNV.EQ.8 ) IPT7.ITYPEL=14 410 : IF (NBNNV.EQ.20) IPT7.ITYPEL=15 411 : IPT7.LISREF(1)=IPT1 412 : C*c ... Mise à 0 des connectivités ... 413 : C* DO 1040 I=1,NBNN 414 : C* DO 1040 J=1,NBELEM 415 : C* IPT7.NUM(I,J)=0 416 : C* 1040 CONTINUE 417 : SEGINI TABPAR 418 : c ... si ce n'est ni TRAN ni ROTA on saute ... 419 : IF (ICLE.EQ.3) GOTO 16 420 : IF (ICLE.EQ.2) GOTO 10 421 : c ... cas TRAN, on fait appel à l'opérateur PLUS ... 422 : IOPTG=1 422 : CALL ECROBJ('POINT ',IVEC) 422 : CALL ECROBJ('MAILLAGE',IPT1) 422 : CALL PROPER(IOPTG) 423 : GOTO 11 424 : c ... cas ROTA, on fait appel à l'opérateur TOURner ... 425 : 10 XXX=ANGLI 425 : CALL ECRREE(XXX) 425 : CALL ECROBJ('POINT ',IP2) 425 : CALL ECROBJ('POINT ',IP1) 425 : CALL ECROBJ('MAILLAGE',IPT1) 425 : CALL TOURNE 426 : 11 CONTINUE 427 : c ... puis on lit le second MAILLAGE ... 427 : CALL LIROBJ('MAILLAGE',IPT2,1,IRETOU) 428 : IF (IERR.NE.0) RETURN 429 : C IPT3 ET IPT4 ONT ETE DESCENDU DANS L'OPERATION AINSI QUE MCOORD/REFPO 430 : 16 SEGACT IPT1,IPT2,MCOORD 431 : IPT4=IPT2 432 : c ... si le 1er maillage est simple on suppose que le deuxième le sera aussi ... 433 : IF (IDEUX.EQ.1) GOTO 15 434 : IPT5=IPT2.LISOUS(1) 435 : IPT6=IPT2.LISOUS(2) 436 : SEGACT IPT5,IPT6 437 : C ON FAIT COMME POUR LE BAS 438 : NBSOUS=0 439 : NBREF=0 440 : NBNN=4*INCR 441 : c ... qui ont les mêmes nombres d'éléments que les sous-maillages du premier ... 442 : NBELEM=NBELE3+NBELE4 443 : SEGINI MELEME 444 : c ... et que l'on transvase succesivement dans une nouvelle entité (MELEME) ... 445 : C* DO 1110 J=1,NBELEM 446 : C* DO 1110 I=1,NBNN 447 : C* NUM(I,J)=0 448 : C* 1110 CONTINUE 449 : c ... d'abord IPT5 ... 450 : DO 1111 J=1,NBELE3 451 : ICOLOR(J)=IPT5.ICOLOR(J) 452 : DO 1111 I=1,IPT5.NUM(/1) 453 : NUM(I,J)=IPT5.NUM(I,J) 454 : 1111 CONTINUE 455 : c ... puis IPT6 ... 456 : DO 1112 J=1,NBELE4 457 : K=J+NBELE3 458 : ICOLOR(K)=IPT6.ICOLOR(J) 459 : DO 1112 I=1,IPT6.NUM(/1) 460 : NUM(I,K)=IPT6.NUM(I,J) 461 : 1112 CONTINUE 462 : SEGDES IPT5,IPT6,IPT2 463 : c ... et qui remplace le maillage lu ... 464 : IPT4=MELEME 465 : 15 IPT7.LISREF(2)=IPT2 466 : 467 : C CONSTRUCTION DE LA TABLE DES POINTS EFFECTIFS 468 : c ... IPT3 = maillage (parfois bâtard) contenant toutes les facettes 469 : c de la surface initiale (?) ... 470 : NBELEC=IPT3.NUM(/2) 471 : c ... ICPR(ligne = nombre maxi de noeuds / facette, colonne = nb facettes) ... 472 : SEGINI ICPR 473 : C*c ... mise à 0 ... 474 : C* DO 12 I=1,NBNNEL 475 : C* DO 12 J=1,NBELEC 476 : C* 12 ICPR(I,J)=0 477 : 478 : c ... on parcourt les 2 maillages ... 479 : DO 13 J=1,NBELEC 480 : DO 13 I=1,NBNNEL 481 : c ... IR = N° du noeud (ou 0) du 1er ... 482 : IR=IPT3.NUM(I,J) 483 : c ... IR2 = N° du noeud (ou 0) equivalent du 2nd ... 484 : IR2=IPT4.NUM(I,J) 485 : c ... si le 1er est absent ... 486 : IF (IR.EQ.0) GOTO 1120 487 : c ... sinon, si son equivalent est nul => kk !!!!!!! ... 488 : IF (IR2.EQ.0) GOTO 8833 489 : I1=IR 490 : I1R2=IR2 491 : c ... si ce n'est pas le 1er élément ... 492 : IF (J.EQ.1) GOTO 13 493 : c ... on va vérifier que l'equivalence est la même pour 494 : c tous les éléments précédents ... 495 : JM1=J-1 496 : DO 14 JJ=1,JM1 497 : DO 14 II=1,NBNNEL 498 : IR=IPT3.NUM(II,JJ) 499 : IR2=IPT4.NUM(II,JJ) 500 : IF (IR.EQ.0) GOTO 14 501 : IF (IR.NE.I1) GOTO 8834 502 : IF (IR2.NE.I1R2) GOTO 8833 503 : c ... on met dans ICPR(n° noeud,n° élt) la valeur de II+(JJ-1)*8 504 : c tq. le noeud II de l'élt JJ de IPT3 est le même que 505 : c le noeud I de l'élt J (toutefois JJ < J donc aucun noeud ne 506 : c pointe sur lui même) ... 507 : ICPR(I,J)=II+(JJ-1)*8 508 : GOTO 13 509 : 8834 IF (IR2.EQ.I1R2) GOTO 8833 510 : 14 CONTINUE 511 : GOTO 13 512 : c ... si le 1er est absent (suite), on met ICPR correspondant à -1 ... 513 : 1120 ICPR(I,J)=-1 514 : c ... si son equivalent est non nul => kk !!!!!!! ... 515 : IF (IR2.NE.0) GOTO 8833 516 : 13 CONTINUE 517 : GOTO 8835 518 : 8833 CONTINUE 519 : 520 : C LES TOPOLOGIES SONT DIFFERENTES 521 : SEGSUP ICPR 521 : CALL ERREUR(21) 522 : RETURN 523 : 8835 CONTINUE 524 : 525 : C ON FABRIQUE POUR LE MOMENT DES CUBES A 8 OU 20 NOEUDS ET DES PRISMES 526 : C A 6 OU 15 NOEUDS 527 : C D'ABORD LES POINTS DU BAS 528 : 529 : DIN=DEN1 530 : DO 20 I=1,NBELEC 531 : c ... On commence par donner la couleur ... 532 : c ... Si c'est TRAN ou ROTA ce sera celle du maillage d'origine ... 533 : IF (ICLE.NE.3) THEN 534 : IPT7.ICOLOR(I)=IPT3.ICOLOR(I) 535 : c ... sinon, une <> au sens de ITABM ... 536 : ELSE 537 : ICOLI=IPT3.ICOLOR(I) 538 : C ... CORRECTION PROBLEME SAUSSAIS 29 NOVEMBRE 1985 539 : ICOLJ=IPT4.ICOLOR(I) 540 : IPT7.ICOLOR(I)=ITABM(ICOLI,ICOLJ) 541 : ENDIF 542 : c ... Puis on commence par transvaser les connectivités de la surface initiale ... 543 : DO 20 J=1,NBNNEL 544 : IR=IPT3.NUM(J,I) 545 : IF (IR.EQ.0) GOTO 20 546 : IPT7.NUM(J,I)=IR 547 : 20 CONTINUE 548 : 549 : IBASE=XCOOR(/1)/IDIMP1 550 : 551 : C ON FABRIQUE ENSUITE LES COUCHES 552 : C ON AFFECTE SEULEMENT LES NUMEROS DE NOEUDS 553 : 554 : c ... IDIF = nombre de noeuds placés < > les couches ... 555 : IDIF=(INCR-1)*(NBNNEL/2) 556 : NX=NCOUCH-1 557 : DO 21 ICOUCH=1,NCOUCH 558 : DIN=DIN*APROG 559 : TABPAR(ICOUCH)=DIN 560 : IF (ICOUCH.EQ.NCOUCH) GOTO 21 561 : JBASE=(ICOUCH-1)*NBELEC 562 : IF (INCR.EQ.1) GOTO 2000 563 : 564 : C ON FABRIQUE D'ABORD LA COUCHE INTERMEDIAIRE 565 : 566 : DO 2001 J=1,NBELEC 567 : IPT7.ICOLOR(J+JBASE)=IPT7.ICOLOR(J) 568 : DO 2001 IA=1,(NBNNEL/2) 569 : I=2*IA-1 570 : IF (ICPR(I,J).EQ.-1) GOTO 2001 571 : IF (ICPR(I,J).NE.0) GOTO 2002 572 : IBASE=IBASE+1 573 : IPT7.NUM(IA+NBNNEL,J+JBASE)=IBASE 574 : GOTO 2001 575 : 2002 IAUX=ICPR(I,J) 576 : JJ=(IAUX-1)/8+1 577 : II=IAUX-8*JJ+8 578 : IIA=(II+1)/2 579 : IPT7.NUM(IA+NBNNEL,J+JBASE)=IPT7.NUM(IIA+NBNNEL,JJ+JBASE) 580 : 2001 CONTINUE 581 : 2000 CONTINUE 582 : DO 22 J=1,NBELEC 583 : IPT7.ICOLOR(J+JBASE)=IPT7.ICOLOR(J) 584 : DO 22 I=1,NBNNEL 585 : IF (ICPR(I,J).EQ.-1) GOTO 22 586 : IF (ICPR(I,J).NE.0) GOTO 23 587 : IBASE=IBASE+1 588 : IPT7.NUM(I,J+JBASE+NBELEC)=IBASE 589 : IPT7.NUM(I+NBNNEL+IDIF,J+JBASE)=IBASE 590 : GOTO 22 591 : 23 IAUX=ICPR(I,J) 592 : JJ=(IAUX-1)/8+1 593 : II=IAUX-8*JJ+8 594 : IPT7.NUM(I,J+JBASE+NBELEC)=IPT7.NUM(II,JJ+JBASE+NBELEC) 595 : IPT7.NUM(I+NBNNEL+IDIF,J+JBASE)=IPT7.NUM(II+NBNNEL+IDIF, 596 : & JJ+JBASE) 597 : 22 CONTINUE 598 : 21 CONTINUE 599 : IF (MLREEL.NE.0)THEN 600 : SEGACT,MLREEL 601 : DPROG=PROG(NCOUCH+1)-PROG(1) 602 : DO 12345 ICOUCH=1,NCOUCH 603 : TABPAR(ICOUCH)=(PROG(ICOUCH+1)-PROG(ICOUCH))/DPROG 604 : 12345 CONTINUE 605 : SEGDES,MLREEL 606 : ENDIF 607 : 25 CONTINUE 608 : C ON FAIT LES POINTS DU HAUT ET EVENTUELLEMENT LA COUCHE INTERMEDIAIRE 609 : C PRECEDENTE 610 : JBASE=NBELEC*NX 611 : IF (INCR.EQ.1) GOTO 2003 612 : DO 2004 J=1,NBELEC 613 : IPT7.ICOLOR(J+JBASE)=IPT7.ICOLOR(J) 614 : DO 2004 IA=1,(NBNNEL/2) 615 : I=2*IA-1 616 : IF (ICPR(I,J).EQ.-1) GOTO 2004 617 : IF (ICPR(I,J).NE.0) GOTO 2005 618 : IBASE=IBASE+1 619 : IPT7.NUM(IA+NBNNEL,J+JBASE)=IBASE 620 : GOTO 2004 621 : 2005 IAUX=ICPR(I,J) 622 : JJ=(IAUX-1)/8+1 623 : II=IAUX-8*JJ+8 624 : IIA=(II+1)/2 625 : IPT7.NUM(IA+NBNNEL,J+JBASE)=IPT7.NUM(IIA+NBNNEL,JJ+JBASE) 626 : 2004 CONTINUE 627 : 2003 CONTINUE 628 : DO 30 J=1,NBELEC 629 : IPT7.ICOLOR(J+JBASE)=IPT7.ICOLOR(J) 630 : DO 30 I=1,NBNNEL 631 : IF (ICPR(I,J).EQ.-1) GOTO 30 632 : IPT7.NUM(I+NBNNEL+IDIF,J+JBASE)=IPT4.NUM(I,J) 633 : 30 CONTINUE 634 : DPAR=0. 635 : C CREATION DES POINTS 636 : IADR=XCOOR(/1)/4 637 : NBPTS=IADR+NCOUCH*INCR*NBELEC*NBNNEL 638 : SEGADJ MCOORD 639 : DO 61 ICOUCH=1,NCOUCH 640 : DIN=TABPAR(ICOUCH) 641 : DO 610 IC=1,INCR 642 : IC1=INCR+1-IC 643 : DPAR=DPAR+DIN/INCR 644 : UMDPAR=1.-DPAR 645 : DEN610=DENI+DECA*DPAR 646 : IF (ICOUCH.EQ.NCOUCH.AND.IC.EQ.INCR) GOTO 610 647 : IF (ICLE.NE.2) GOTO 63 648 : ANG=DPAR*DLONG/RAY 649 : SI=SIN(ANG) 650 : CO=COS(ANG) 651 : 63 CONTINUE 652 : DO 620 J=1,NBELEC 653 : DO 62 I=1,NBNNEL,IC1 654 : IF (ICPR(I,J).NE.0) GOTO 62 655 : IREF=4*IPT3.NUM(I,J)-4 656 : GOTO (67,64,66),ICLE 657 : 67 XCOOR(IADR*IDIMP1+1)=XCOOR(IREF+1)+DPAR*XDIS 658 : XCOOR(IADR*IDIMP1+2)=XCOOR(IREF+2)+DPAR*YDIS 659 : XCOOR(IADR*IDIMP1+3)=XCOOR(IREF+3)+DPAR*ZDIS 660 : GOTO 65 661 : 66 IREF2=4*IPT4.NUM(I,J)-4 662 : XCOOR(IADR*IDIMP1+1)=UMDPAR*XCOOR(IREF+1)+DPAR*XCOOR(IREF2+1) 663 : XCOOR(IADR*IDIMP1+2)=UMDPAR*XCOOR(IREF+2)+DPAR*XCOOR(IREF2+2) 664 : XCOOR(IADR*IDIMP1+3)=UMDPAR*XCOOR(IREF+3)+DPAR*XCOOR(IREF2+3) 665 : GOTO 65 666 : 64 X1=XCOOR(IREF+1)-XPT1 667 : Y1=XCOOR(IREF+2)-YPT1 668 : Z1=XCOOR(IREF+3)-ZPT1 669 : XV=X1*XV1+Y1*YV1+Z1*ZV1 670 : YV=X1*XV2+Y1*YV2+Z1*ZV2 671 : ZV=X1*XVEC+Y1*YVEC+Z1*ZVEC 672 : XD=XV*CO-YV*SI 673 : YD=XV*SI+YV*CO 674 : ZD=ZV 675 : XCOOR(IADR*IDIMP1+1)=XD*XV1+YD*XV2+ZD*XVEC+XPT1 676 : XCOOR(IADR*IDIMP1+2)=XD*YV1+YD*YV2+ZD*YVEC+YPT1 677 : XCOOR(IADR*IDIMP1+3)=XD*ZV1+YD*ZV2+ZD*ZVEC+ZPT1 678 : GOTO 65 679 : 65 CONTINUE 680 : XCOOR((IADR+1)*IDIMP1)=DEN610 681 : IADR=IADR+1 682 : 62 CONTINUE 683 : 620 CONTINUE 684 : 610 CONTINUE 685 : 61 CONTINUE 686 : NBPTS=IADR 687 : SEGADJ MCOORD 688 : 60 CONTINUE 689 : C C'EST FINI 690 : C IL RESTE DANS LE CAS OU ON A DES CUBES ET DES PRISMES A LES SEPARER 691 : C ET A SUPPRIMER LES SEGMENTS SUPPLEMENTAIRES DE TRAVAIL 692 : C D'ABORD FAIRE LE POURTOUR A PARTIR DU CONTOUR 693 : IF (IPT7.LISREF(/1).EQ.2) GOTO 3000 693 : CALL ECROBJ('MAILLAGE',IPT1) 693 : CALL PRCONT 693 : CALL LIROBJ('MAILLAGE',IPT5,1,IRETOU) 694 : IF (IERR.NE.0) GOTO 3000 695 : C IPT5 LE CONTOUR IPT6 SERA LE POURTOUR 696 : SEGACT IPT5 697 : NBASE=IPT5.NUM(/2) 698 : NBNN=INCR*4 699 : NBELEM=NBASE*NCOUCH 700 : NBSOUS=0 701 : NBREF=0 702 : SEGINI IPT6 703 : IPT6.ITYPEL=6+2*INCR 704 : SEGACT IPT3 705 : DO 3001 IEL=1,NBASE 706 : DO 3001 IP=1,INCR+1 707 : INP=IPT5.NUM(IP,IEL) 708 : DO 3003 IELS=1,NBELEC 709 : DO 3003 IPS=1,NBNNEL 710 : IPSP=IPT3.NUM(IPS,IELS) 711 : IF (IPSP.EQ.0) GOTO 3003 712 : IF (IPSP.EQ.INP) GOTO 3002 713 : 3003 CONTINUE 714 : GOTO 3000 715 : 3002 CONTINUE 716 : DO 3004 IC=1,NCOUCH 717 : IBASE=(IC-1)*NBASE 718 : JBASE=(IC-1)*NBELEC 719 : C PTS DU BAS 720 : IPT6.NUM(IP,IEL+IBASE)=IPT7.NUM(IPS,IELS+JBASE) 721 : C PTS DU HAUT 722 : IPT6.NUM(NBNN+2-INCR-IP,IEL+IBASE)= 723 : # IPT7.NUM(IPS+NBNNEL+IDIF,IELS+JBASE) 724 : C EVENTUELLEMENT PTS MILIEUX 725 : IF (INCR.EQ.1.OR.IP.EQ.2) GOTO 3004 726 : IPT6.NUM(10-2*IP,IEL+IBASE)=IPT7.NUM((IPS+1)/2+NBNNEL,IELS+JBASE) 727 : 3004 CONTINUE 728 : 3001 CONTINUE 729 : DO 3005 I=1,NCOUCH 730 : DO 3005 J=1,NBASE 731 : IPT6.ICOLOR(J+(I-1)*NBASE)=IPT5.ICOLOR(J) 732 : 3005 CONTINUE 733 : SEGDES IPT5,IPT6 734 : IPT7.LISREF(3)=IPT6 735 : 3000 CONTINUE 736 : IF (IDEUX.EQ.1) GOTO 1500 737 : SEGSUP IPT3,IPT4 738 : MELEME=IPT7 739 : NBSOUS=2 740 : NBREF=LISREF(/1) 741 : NBNN=0 742 : NBELEM=0 743 : SEGINI IPT7 744 : IPT7.LISREF(1)=LISREF(1) 745 : IPT7.LISREF(2)=LISREF(2) 746 : IF (NBREF.EQ.3) IPT7.LISREF(3)=LISREF(3) 747 : NBSOUS=0 748 : NBREF=0 749 : NBNN=6 750 : IF (INCR.EQ.2) NBNN=15 751 : NBELEM=NBTRI*NCOUCH 752 : SEGINI IPT3 753 : IPT3.ITYPEL=16 754 : IF (INCR.EQ.2) IPT3.ITYPEL=17 755 : IPT7.LISOUS(1)=IPT3 756 : NBNN=8 757 : IF (INCR.EQ.2) NBNN=20 758 : NBELEM=NBQUA*NCOUCH 759 : SEGINI IPT4 760 : IPT4.ITYPEL=14 761 : IF (INCR.EQ.2) IPT4.ITYPEL=15 762 : IPT7.LISOUS(2)=IPT4 763 : IT=0 764 : IQ=0 765 : DO 1501 J=1,NUM(/2) 766 : IF (NUM(NBNNV,J).EQ.0) GOTO 1502 767 : C C'EST UN CUBE 768 : IQ=IQ+1 769 : IPT4.ICOLOR(IQ)=ICOLOR(J) 770 : DO 1503 K=1,IPT4.NUM(/1) 771 : IPT4.NUM(K,IQ)=NUM(K,J) 772 : 1503 CONTINUE 773 : GOTO 1501 774 : 1502 IT=IT+1 775 : IPT3.ICOLOR(IT)=ICOLOR(J) 776 : C C'EST UN PRISME 777 : IF (INCR.EQ.2) GOTO 2020 778 : IPT3.NUM(1,IT)=NUM(1,J) 779 : IPT3.NUM(2,IT)=NUM(2,J) 780 : IPT3.NUM(3,IT)=NUM(3,J) 781 : IPT3.NUM(4,IT)=NUM(NBNNEL+1,J) 782 : IPT3.NUM(5,IT)=NUM(NBNNEL+2,J) 783 : IPT3.NUM(6,IT)=NUM(NBNNEL+3,J) 784 : GOTO 1501 785 : 2020 CONTINUE 786 : DO 2021 L=1,6 787 : IPT3.NUM(L,IT)=NUM(L,J) 788 : 2021 CONTINUE 789 : IPT3.NUM(7,IT)=NUM(NBNNEL+1,J) 790 : IPT3.NUM(8,IT)=NUM(NBNNEL+2,J) 791 : IPT3.NUM(9,IT)=NUM(NBNNEL+3,J) 792 : DO 2022 L=1,6 793 : IPT3.NUM(L+9,IT)=NUM(NBNNEL+IDIF+L,J) 794 : 2022 CONTINUE 795 : 1501 CONTINUE 796 : SEGDES IPT3,IPT4 797 : SEGSUP MELEME 798 : 1500 SEGDES IPT1,IPT2 799 : SEGSUP ICPR,TABPAR 800 : IF (ISVOL1.EQ.0) GOTO 3200 801 : IPT8=ISVOL1 802 : SEGACT IPT8 803 : ltelq=.false. 803 : CALL FUSE(IPT8,IPT7,IRET,ltelq) 804 : SEGDES IPT7,IPT8 805 : IPT7=IRET 806 : 3200 CONTINUE 807 : IF (ISVOL2.EQ.0) GOTO 3201 808 : IPT8=ISVOL2 809 : SEGACT IPT8 810 : ltelq=.false. 810 : CALL FUSE(IPT7,IPT8,IRET,ltelq) 811 : SEGDES IPT7,IPT8 812 : IPT7=IRET 813 : 3201 CONTINUE 814 : SEGDES IPT7 814 : CALL ECROBJ('MAILLAGE',IPT7) 815 : RETURN 816 : 817 : 4400 CONTINUE 818 : mchpoi=0 819 : epai=0.d0 819 : Call LIRREE(EPAI,0,iretou) 819 : if(iretou.eq.0) call lirobj('CHPOINT ' , MCHPOI,0,iretch) 820 : if(iretou+iretch.eq.0) then 821 : C MAILLAGE AUTOMATIQUE DE VOLUME 822 : write(IOIMP,*) ' appel a demete' 822 : CALL DEMETE(IPT1) 823 : IF (IERR.NE.0) RETURN 824 : IPT7=IPT1 825 : GOTO 3201 826 : else 826 : call lirobj('POINT ',ip1,1,iretpt) 827 : if(ierr.ne.0) return 827 : call volshb(ipt1,epai,mchpoi,ip1,ipt7) 828 : if(ierr.ne.0) return 829 : go to 3201 830 : endif 831 : END 832 : 833 :
© Cast3M 2003 - All rights reserved.
Disclaimer