1 : C INCLUS SOURCE PASCAL 09/01/22 21:15:36 6263 2 : C 3 : C CE SOUS PROGRAMME EXTRAIT DU PREMIER OBJET L'ENSEMBLE DES ELEMENTS 4 : C STRICTEMENT INCLUS DANS LE DEUXIEME OBJET 5 : C 6 : SUBROUTINE INCLUS 7 : IMPLICIT INTEGER(I-N) 8 : IMPLICIT REAL*8 (A-H,O-Z) 9 : -INC CCREEL 10 : *- 11 : -INC SMELEME 12 : -INC CCOPTIO 13 : -INC SMCOORD 14 : -INC CCGEOME 15 : CHARACTER*4 LEMOT(5) 16 : DATA LEMOT/'BARY','STRI','LARG','VOLU','NOID'/ 17 : SEGMENT /FER/(NFI(ITT),MAI(IPP),ITOUR) 18 : SEGMENT ICPR(XCOOR(/1)/(IDIM+1)) 19 : 19 : CALL LIROBJ('MAILLAGE',IPT1,1,IRETOU) 19 : CALL LIROBJ('MAILLAGE',IPT2,1,IRETOU) 20 : IF (IERR.NE.0) RETURN 21 : 22 : * On fixe les options 23 : * Par défaut : non barycentrique, appartenance large, et erreur si 24 : * maillage vide 25 : IBARY=0 26 : IFRON=0 27 : IVERI=0 28 : DO I=1,3 28 : CALL LIRMOT(LEMOT,5,IRET,0) 29 : IF (IRET.EQ.1) IBARY=1 30 : IF (IRET.EQ.2) IFRON=1 31 : IF (IRET.EQ.3) IFRON=0 32 : IF (IRET.EQ.4) THEN 33 : * Cas de maillage volumique 33 : CALL INCLU3(IPT1,IPT2) 34 : RETURN 35 : ENDIF 36 : IF (IRET.EQ.5) IVERI=1 37 : IF (IRET.EQ.0) GOTO 20 38 : ENDDO 39 : 40 : 20 CONTINUE 41 : 42 : C CRITERE D'INCLUSION PAS UTILISE EN 2D... 42 : CALL LIRREE(XCRITT,0,IRET) 43 : XCRITT=0.D0 44 : 45 : * Cas d'une surface plane. 46 : SEGACT IPT2 47 : 48 : C Appel éventuel à CONTOUR 49 : IF ((KSURF(IPT2.ITYPEL).EQ.0).AND.(IPT2.LISOUS(/1).EQ.0)) GOTO 25 49 : CALL ECROBJ('MAILLAGE',IPT2) 49 : CALL PRCONT 49 : CALL LIROBJ('MAILLAGE',IPT2,1,IRETOU) 50 : IF (IERR.NE.0) RETURN 51 : SEGACT IPT2 52 : 53 : 25 CONTINUE 54 : * Transformation éventuelle en SEG2 (SEG3 ou maillage complexe) 55 : IF(IPT2.ITYPEL.EQ.3.OR.IPT2.LISOUS(/1).NE.0) 55 : $ CALL CHANGE (IPT2,2) 56 : 57 : * Remplissage segment FER 58 : ICLE=0 58 : CALL AVTRSF(IPT2,FER,ICLE) 59 : IF (IERR.NE.0) RETURN 60 : 61 : * Détermination du plan par son centre de gravité et son vecteur 62 : * normal unitaire (pour changement de variable) 63 : SEGACT MCOORD 64 : SEGINI ICPR 65 : DO 1 I=1,ICPR(/1) 66 : ICPR(I)=0 67 : 1 CONTINUE 68 : XGRAV=0. 69 : YGRAV=0. 70 : ZGRAV=0. 71 : DO 30 I=MAI(1)+1,MAI(ITOUR+1) 72 : IREF=NFI(I)*(IDIM+1)-IDIM 73 : XGRAV=XGRAV+XCOOR(IREF) 74 : YGRAV=YGRAV+XCOOR(IREF+1) 75 : ZGRAV=ZGRAV+XCOOR(IREF+2) 76 : 30 CONTINUE 77 : XGRAV=XGRAV/FLOAT(MAI(ITOUR+1)-MAI(1)) 78 : YGRAV=YGRAV/FLOAT(MAI(ITOUR+1)-MAI(1)) 79 : ZGRAV=ZGRAV/FLOAT(MAI(ITOUR+1)-MAI(1)) 80 : IF (IDIM.EQ.2) ZGRAV=0. 81 : XNORM=0.D0 82 : YNORM=0.D0 83 : ZNORM=0.D0 84 : DO 42 IT=1,ITOUR 85 : IREF=NFI(MAI(IT+1))*(IDIM+1)-IDIM 86 : XV1=XCOOR(IREF)-XGRAV 87 : YV1=XCOOR(IREF+1)-YGRAV 88 : ZV1=XCOOR(IREF+2)-ZGRAV 89 : IF (IDIM.EQ.2) ZV1=0. 90 : DO 421 I=MAI(IT-1+1)+1,MAI(IT+1) 91 : IREF=(NFI(I))*(IDIM+1)-IDIM 92 : XV2=XCOOR(IREF)-XGRAV 93 : YV2=XCOOR(IREF+1)-YGRAV 94 : ZV2=0. 95 : IF (IDIM.GE.3) ZV2=XCOOR(IREF+2)-ZGRAV 96 : XNORM=XNORM+YV1*ZV2-ZV1*YV2 97 : YNORM=YNORM+ZV1*XV2-XV1*ZV2 98 : ZNORM=ZNORM+XV1*YV2-YV1*XV2 99 : XV1=XV2 100 : YV1=YV2 101 : ZV1=ZV2 102 : 421 CONTINUE 103 : 42 CONTINUE 104 : DNORM=SQRT(XNORM**2+YNORM**2+ZNORM**2) 105 : XNORM=XNORM/DNORM 106 : YNORM=YNORM/DNORM 107 : ZNORM=ZNORM/DNORM 108 : 109 : MELEME=IPT1 110 : SEGACT MELEME 111 : NBS =LISOUS(/1) 112 : NBSOUS=NBS 113 : IPT8=0 114 : IF (NBSOUS.NE.0) THEN 115 : NBREF=0 116 : NBNN=0 117 : NBELEM=0 118 : SEGINI IPT8 119 : IOB1=0 120 : ENDIF 121 : 122 : IF(IBARY.EQ.0) THEN 123 : * Eligibilité d'après chaque point de l'élément 124 : DO 7001 IOB=1,MAX(1,NBS) 125 : IF (NBS.NE.0) THEN 126 : IPT1=LISOUS(IOB) 127 : SEGACT IPT1 128 : ENDIF 129 : NBELEM=IPT1.NUM(/2) 130 : NBNN=IPT1.NUM(/1) 131 : NBREF=0 132 : NBSOUS=0 133 : SEGINI IPT4 134 : IPT4.ITYPEL=IPT1.ITYPEL 135 : DO 2 I=1,NBNN 136 : DO 2222 J=1,NBELEM 137 : IPT4.NUM(I,J)=0 138 : 2222 CONTINUE 139 : 2 CONTINUE 140 : 141 : ICO=0 142 : DO 10 J=1,NBELEM 143 : DO 11 I=1,NBNN 144 : IP=IPT1.NUM(I,J) 145 : IF (ICPR(IP).EQ.-1) GOTO 10 146 : IF (ICPR(IP).EQ.1) GOTO 11 147 : IREF=IP*(IDIM+1)-IDIM 148 : XP=XCOOR(IREF) 149 : YP=XCOOR(IREF+1) 150 : ZP=XCOOR(IREF+2) 151 : IF (IDIM.EQ.2) ZP=0. 152 : ANG=0.D0 153 : DO 12 IT=1,ITOUR 154 : IF (IP.EQ.NFI(MAI(IT+1)))THEN 155 : IF(IFRON.EQ.0) GO TO 13 156 : GO TO 18 157 : ENDIF 158 : IREF=NFI(MAI(IT+1))*(IDIM+1)-IDIM 159 : XV1=XCOOR(IREF)-XP 160 : YV1=XCOOR(IREF+1)-YP 161 : ZV1=XCOOR(IREF+2)-ZP 162 : IF (IDIM.EQ.2) ZV1=0. 163 : DO 120 M=MAI(IT-1+1)+1,MAI(IT+1) 164 : IF (IP.EQ.NFI(M))THEN 165 : IF(IFRON.EQ.0) GO TO 13 166 : GO TO 18 167 : ENDIF 168 : IREF=NFI(M)*(IDIM+1)-IDIM 169 : XV2=XCOOR(IREF)-XP 170 : YV2=XCOOR(IREF+1)-YP 171 : ZV2=XCOOR(IREF+2)-ZP 172 : IF (IDIM.EQ.2) ZV2=0. 173 : V1=XNORM*(YV1*ZV2-ZV1*YV2)+YNORM*(ZV1*XV2-XV1 174 : $ *ZV2)+ZNORM*(XV1*YV2-YV1*XV2) 175 : V2=XV1*XV2+YV1*YV2+ZV1*ZV2 176 : IF( IFRON.EQ.1) THEN 177 : V6 = ( XV1 -XV2)*( XV1 -XV2) + ( YV1 -YV2)* ( 178 : $ YV1 -YV2)+ ( ZV1-ZV2)* ( ZV1-ZV2) 179 : IF( V2.LE.0. . OR. ( V2 . LT . V6* 1.D-4)) 180 : $ THEN 181 : V4 =SQRT (XV1* XV1 + YV1 * YV1 + ZV1 * ZV1 182 : $ ) 183 : V5 =SQRT (XV2* XV2 + YV2 * YV2 + ZV2 184 : $ * ZV2 ) 185 : V3 = ((YV1*ZV2-ZV1*YV2)*(YV1*ZV2-ZV1*YV2) 186 : $ +(ZV1*XV2-XV1*ZV2)*(ZV1*XV2-XV1*ZV2) 187 : $ +(XV1*YV2-YV1*XV2) *(XV1*YV2-YV1*XV2) 188 : $ ) **0.5D0 189 : IF( V4* 1.D4. LT. V5 . OR . V5*1.D4 .LT.V4 190 : $ ) GO TO 18 191 : IF( V3*1.D4 . LT . ( ABS ( V2) ) ) 192 : # go to 18 193 : ENDIF 194 : ENDIF 195 : IF (V1.EQ.0.D0.AND.V2.EQ.0.D0) GOTO 121 196 : ANG=ANG+ATAN2(V1,V2) 197 : 121 CONTINUE 198 : XV1=XV2 199 : YV1=YV2 200 : ZV1=ZV2 201 : 120 CONTINUE 202 : 12 CONTINUE 203 : IF (ABS(ANG).GT.XPI) GOTO 13 204 : 18 ICPR(IP)=-1 205 : GOTO 10 206 : 13 ICPR(IP)=1 207 : 11 CONTINUE 208 : * Cet élément est ok, on le met de côté. 209 : DO 14 L=1,NBNN 210 : IPT4.NUM(L,J)=IPT1.NUM(L,J) 211 : 14 CONTINUE 212 : IPT4.ICOLOR(J)=IPT1.ICOLOR(J) 213 : ICO=ICO+1 214 : 10 CONTINUE 215 : 216 : NBELEM=ICO 217 : IF (NBELEM.EQ.0) GOTO 7001 218 : 219 : * sauvegarde partition ok 220 : SEGINI IPT3 221 : IPT3.ITYPEL=IPT1.ITYPEL 222 : JJ=0 223 : DO 21 J=1,IPT4.NUM(/2) 224 : IF (IPT4.NUM(1,J).EQ.0) GOTO 21 225 : JJ=JJ+1 226 : IPT3.ICOLOR(JJ)=IPT4.ICOLOR(J) 227 : DO 22 I=1,NBNN 228 : IPT3.NUM(I,JJ)=IPT4.NUM(I,J) 229 : 22 CONTINUE 230 : 21 CONTINUE 231 : SEGSUP IPT4 232 : 233 : IF (NBS.NE.0) THEN 234 : SEGDES IPT1,IPT3 235 : IOB1=IOB1+1 236 : IPT8.LISOUS(IOB1)=IPT3 237 : ENDIF 238 : 7001 CONTINUE 239 : ELSE 240 : * Eligibilité d'après la position du barycentre 241 : DO 7051 IOB=1,MAX(1,NBS) 242 : IF (NBS.NE.0) THEN 243 : IPT1=LISOUS(IOB) 244 : SEGACT IPT1 245 : ENDIF 246 : NBELEM=IPT1.NUM(/2) 247 : NBNN=IPT1.NUM(/1) 248 : NBREF=0 249 : NBSOUS=0 250 : SEGINI IPT4 251 : IPT4.ITYPEL=IPT1.ITYPEL 252 : DO 52 I=1,NBNN 253 : DO 521 J=1,NBELEM 254 : IPT4.NUM(I,J)=0 255 : 521 CONTINUE 256 : 52 CONTINUE 257 : ICO=0 258 : DO 60 J=1,NBELEM 259 : XP=0.D0 260 : YP=0.D0 261 : ZP=0.D0 262 : DO 61 I=1,NBNN 263 : IP=IPT1.NUM(I,J) 264 : IREF=IP*(IDIM+1)-IDIM 265 : XP=XCOOR(IREF) +XP 266 : YP=XCOOR(IREF+1)+YP 267 : ZP=XCOOR(IREF+2) +ZP 268 : 61 CONTINUE 269 : IF (IDIM.EQ.2) ZP=0.D0 270 : XP = XP / FLOAT(NBNN) 271 : YP = YP / FLOAT(NBNN) 272 : ZP = ZP / FLOAT(NBNN) 273 : ANG=0.D0 274 : DO 62 IT=1,ITOUR 275 : IREF=NFI(MAI(IT+1))*(IDIM+1)-IDIM 276 : XV1=XCOOR(IREF)-XP 277 : YV1=XCOOR(IREF+1)-YP 278 : ZV1=XCOOR(IREF+2)-ZP 279 : IF (IDIM.EQ.2) ZV1=0. 280 : DO 621 M=MAI(IT-1+1)+1,MAI(IT+1) 281 : IREF=NFI(M)*(IDIM+1)-IDIM 282 : XV2=XCOOR(IREF)-XP 283 : YV2=XCOOR(IREF+1)-YP 284 : ZV2=XCOOR(IREF+2)-ZP 285 : IF (IDIM.EQ.2) ZV2=0. 286 : V1=XNORM*(YV1*ZV2-ZV1*YV2)+YNORM*(ZV1*XV2-XV1*ZV2) 287 : $ +ZNORM*(XV1*YV2-YV1*XV2) 288 : V2=XV1*XV2+YV1*YV2+ZV1*ZV2 289 : IF (V1.EQ.0.D0.AND.V2.EQ.0.D0) GOTO 171 290 : ANG=ANG+ATAN2(V1,V2) 291 : 171 CONTINUE 292 : XV1=XV2 293 : YV1=YV2 294 : ZV1=ZV2 295 : 621 CONTINUE 296 : 62 CONTINUE 297 : IF (ABS(ANG).GT.XPI) GOTO 63 298 : GOTO 60 299 : 63 CONTINUE 300 : * Cet élément est ok, on le met de côté. 301 : DO 64 L=1,NBNN 302 : IPT4.NUM(L,J)=IPT1.NUM(L,J) 303 : 64 CONTINUE 304 : IPT4.ICOLOR(J)=IPT1.ICOLOR(J) 305 : ICO=ICO+1 306 : 60 CONTINUE 307 : 308 : NBELEM=ICO 309 : IF (NBELEM.EQ.0) GOTO 7051 310 : 311 : * sauvegarde partition ok 312 : SEGINI IPT3 313 : IPT3.ITYPEL=IPT1.ITYPEL 314 : JJ=0 315 : DO 71 J=1,IPT4.NUM(/2) 316 : IF (IPT4.NUM(1,J).EQ.0) GOTO 71 317 : JJ=JJ+1 318 : IPT3.ICOLOR(JJ)=IPT4.ICOLOR(J) 319 : DO 72 I=1,NBNN 320 : IPT3.NUM(I,JJ)=IPT4.NUM(I,J) 321 : 72 CONTINUE 322 : 71 CONTINUE 323 : SEGSUP IPT4 324 : IF (NBS.NE.0) THEN 325 : SEGDES IPT1,IPT3 326 : IOB1=IOB1+1 327 : IPT8.LISOUS(IOB1)=IPT3 328 : ENDIF 329 : 7051 CONTINUE 330 : ENDIF 331 : 332 : * Ecriture du maillage résultat et sortie 333 : * --------------------------------------- 334 : IF (((NBS.EQ.0).AND.(NBELEM.EQ.0)).OR. 335 : & ((NBS.NE.0).AND.(IOB1.EQ.0))) THEN 336 : * on n'a rien trouvé 337 : IF (IVERI.EQ.0) THEN 338 : * Tache impossible. Probablement données erronées 338 : CALL ERREUR(26) 339 : RETURN 340 : ELSE 341 : * on renvoie un maillage vide 342 : GOTO 7006 343 : ENDIF 344 : ENDIF 345 : IF (NBS.EQ.0) THEN 346 : * Si maillage simple en entrée 347 : * On écrit directement IPT3 348 : GOTO 7002 349 : ELSE 350 : * Si maillage complexe en entrée 351 : * IOB1 : nombre de partitions retenues 352 : IF (IOB1.EQ.NBS) THEN 353 : GOTO 7003 354 : ELSEIF (IOB1.EQ.1) THEN 355 : GOTO 7005 356 : ELSE 357 : GOTO 7004 358 : ENDIF 359 : ENDIF 360 : 361 : * Maillage complexe, avec moins de partitions que l'objet entré 362 : 7004 CONTINUE 363 : NBSOUS=IOB1 364 : NBREF=0 365 : NBNN=0 366 : NBELEM=0 367 : SEGINI IPT3 368 : DO IOB=1,NBSOUS 369 : IPT3.LISOUS(IOB)=IPT8.LISOUS(IOB) 370 : ENDDO 371 : SEGSUP IPT8 372 : GOTO 7002 373 : 374 : * Maillage simple, avec moins de partitions que l'objet entré 375 : 7005 CONTINUE 376 : IPT3=IPT8.LISOUS(1) 377 : SEGSUP IPT8 378 : GOTO 7002 379 : 380 : * Maillage avec autant de partitions que l'objet entré 381 : * (inutile de recréer un nouveau maillage) 382 : 7003 CONTINUE 383 : IPT3=IPT8 384 : GOTO 7002 385 : 386 : * Maillage vide 387 : 7006 CONTINUE 388 : NBSOUS=0 389 : NBREF=0 390 : NBNN=0 391 : NBELEM=0 392 : SEGINI IPT3 393 : IF (NBS.NE.0) SEGSUP IPT8 394 : 395 : 7002 CONTINUE 396 : SEGDES MELEME,IPT3 397 : SEGSUP ICPR,FER 397 : CALL ECROBJ('MAILLAGE',IPT3) 398 : SEGDES IPT3 399 : 400 : RETURN 401 : END 402 : 403 : 404 : 405 : 406 : 407 : 408 : 409 : 410 : 411 : 412 :
© Cast3M 2003 - All rights reserved.
Disclaimer