1 : C KALPRE SOURCE CHAT 11/03/16 21:25:56 6902 2 : SUBROUTINE KALPRE (MYMOD,INFOEL,SKFACE,XDEC,NELD) 3 : IMPLICIT INTEGER(I-N) 4 : IMPLICIT REAL*8 (A-H,O-Z) 5 : C---------------------------------------------------------------------- 6 : C Facteurs de forme OPTION CACHE 3D 7 : C Calcul du decoupage des faces et initialisations 8 : C SP appelle par facgen 9 : C entree: 10 : C MYMOD : pointeur de l'objet modèle 11 : C INFOEL : informations concernant le type des éléments des maillages 12 : C SKFACE : pointeur sur l'objet 'faces' pour le calcul des FF 13 : C SHC3D : pointeur sur la surface de projection 14 : C XDEC : paramŠtre pour le decoupage des faces 15 : C NELD : nombre d'éléments virtuels 16 : C (un élément COQ -> 2 facteurs de forme <=> 2 éléments virtuels) 17 : C sortie: 18 : C SKFACE : pointeur sur l'objet 'faces' pour le calcul des FF 19 : C SHC3D : pointeur sur la surface de projection 20 : C---------------------------------------------------------------------- 21 : C 22 : -INC CCOPTIO 23 : -INC SMCOORD 24 : -INC SMMODEL 25 : -INC SMELEME 26 : -INC TFFOR3D 27 : C 28 : C ----------------------------------- 29 : SEGMENT SDECOU 30 : INTEGER KDECOU(MFACE) 31 : ENDSEGMENT 32 : POINTEUR MYMOD.MMODEL 33 : C ----------------------------------- 34 : C Stockage d'informations concernant le type des éléments des maillages 35 : SEGMENT ,INFOEL 36 : LOGICAL KCOQ(N1),KQUAD(N1) 37 : ENDSEGMENT 38 : C ----------------------------------- 39 : DIMENSION A(3,4),UU(4),IN(3) 40 : DIMENSION A1(3,3),A2(3,3),U1(4),U2(4),G1(3),G2(3) 41 : DIMENSION GG(3,400),SS(400) 42 : LOGICAL ICOQ 43 : C 44 : C-------------------------------------------------------------------- 45 : C 46 : C 47 : DMIN=1.E-5 48 : DMAX=1./DMIN 49 : C... nombre max de triangles de subdivision 50 : NPM=20 51 : C... quadratiques 52 : NSPA=1 53 : 54 : IF(INFOEL.EQ.0) THEN 55 : ICOQ=.FALSE. 56 : ELSE 57 : ICOQ=.TRUE. 58 : SEGACT INFOEL 59 : ENDIF 60 : C 61 : C>>> CREATION DE L'OBJET FACE 62 : C 63 : NFACE = 0 64 : NELD = 0 65 : C 66 : SEGACT MYMOD 67 : MTYP = MYMOD.KMODEL(/1) 68 : SEGACT INFOEL 69 : C 70 : C On va compter le nombre d'éléments virtuels (égal au nombre de 71 : C facteurs de forme par ligne) : NELD , et le nombre de faces : NFACE . 72 : DO 10 ITYP=1,MTYP 73 : IMODE1 = MYMOD.KMODEL(ITYP) 74 : SEGACT IMODE1 75 : IPT1 = IMODE1.IMAMOD 76 : SEGDES IMODE1 77 : SEGACT IPT1 78 : NEL = IPT1.NUM(/2) 79 : NPGEO = IPT1.NUM(/1) 80 : SEGDES IPT1 81 : C 82 : C un élément COQ -> 2 facteurs de formes 83 : IF (KCOQ(ITYP)) NEL=2*NEL 84 : NELD = NELD + NEL 85 : C 86 : C un élément carré -> 2 éléments triangles 87 : IF(NPGEO.EQ.4.OR.NPGEO.EQ.8) NEL = 2*NEL 88 : NFACE=NFACE+NEL 89 : 10 CONTINUE 90 : C 91 : IF (KIMP.GE.1) THEN 92 : WRITE( 6,*) ' NOMBRE TOTAL D ELEMENTS ',NELD,' DE FACES ',NFACE 93 : ENDIF 94 : C 95 : C 96 : MFACE = NFACE 97 : MS = 3 98 : SEGINI SKFACE 99 : SEGINI SDECOU 100 : C 101 : C KEL : INDICE ELEMENT DECOUPAGE (-> numéro global de la facette) 102 : C IEL : INDICE ELEMENT D'ORIGINE (-> numero global de l'élément virtuel) 103 : C Ces indices sont différents car on découpe les quadrangles en triangles 104 : C 105 : KEL = 1 106 : IEL = 0 107 : C 108 : DO 100 ITYP = 1, MTYP 109 : IMODE1 = MYMOD.KMODEL(ITYP) 110 : SEGACT IMODE1 111 : IPT1 = IMODE1.IMAMOD 112 : SEGDES IMODE1 113 : SEGACT IPT1 114 : 115 : NSGEO = IPT1.NUM(/1) 116 : IF(ICOQ.AND.KQUAD(ITYP)) THEN 117 : NS = NSGEO/2 118 : ELSE 119 : NS=NSGEO 120 : ENDIF 121 : IF(ICOQ.AND.KQUAD(ITYP)) NSPA=2 122 : NEL = IPT1.NUM(/2) 123 : C 124 : C... TRI3 ou TRI6 125 : IF (NS.EQ.3) THEN 126 : *//////// 127 : * WRITE (6,*) 'Le maillage ',IPT1, 128 : * # 'est constitué d éléments triangles' 129 : *//////// 130 : DO 110 I = 1,NEL 131 : IF (KCOQ(ITYP)) THEN 132 : IEL1 = IEL + (2*I-1) 133 : ELSE 134 : IEL1 = IEL + I 135 : ENDIF 136 : KCORR(KEL) = IEL1 137 : 138 : DO 111 IS = 1,NSGEO,NSPA 139 : LS=IS 140 : IF(ICOQ.AND.KQUAD(ITYP)) LS=(IS+1)/2 141 : C** NUK(IS,KEL) = IPT1.NUM(IS,I) 142 : NUK(LS,KEL) = IPT1.NUM(IS,I) 143 : IREF = (IDIM+1)*(NUK(LS,KEL)-1) 144 : DO 112 K = 1,KES 145 : C** A(K,LS) = XCOOR(IREF+K) 146 : A(K,LS) = XCOOR(IREF+K) 147 : 112 CONTINUE 148 : *//////// 149 : * WRITE (6,*) 'noeud numero',IPT1.NUM(IS,I) 150 : * WRITE (6,*) 'Coordonnées :',(A(K,IS),K=1,KES) 151 : *//////// 152 : 111 CONTINUE 153 : IF (KIMP.GE.3) THEN 154 : WRITE (6,*) 'La facette ',KEL,'repose sur les points', 155 : # (IPT1.NUM(IS,I),IS=1,NS) 156 : ENDIF 156 : CALL KNORM(KES,A,NS,UU,S(KEL)) 157 : DO 113 L = 1,4 158 : U(L,KEL) = UU(L) 159 : 113 CONTINUE 160 : C WRITE(6,*) ' KEL S U ',KEL,S(KEL),(U(K,NEL),K=1,4) 161 : 162 : KEL = KEL + 1 163 : *//////// 164 : IF (KCOQ(ITYP)) THEN 165 : C On remplit KCOOR , NUK , U et S pour l'élément inverse 166 : KCORR(KEL) = IEL1 + 1 167 : NUK(1,KEL) = NUK(3,KEL-1) 168 : NUK(2,KEL) = NUK(2,KEL-1) 169 : NUK(3,KEL) = NUK(1,KEL-1) 170 : DO 114 L = 1,4 171 : U(L,KEL) = - UU(L) 172 : 114 CONTINUE 173 : S(KEL) = S(KEL-1) 174 : KEL = KEL + 1 175 : ENDIF 176 : C 177 : 110 CONTINUE 178 : C 179 : IF (KCOQ(ITYP)) THEN 180 : IEL = IEL + NEL*2 181 : ELSE 182 : IEL = IEL + NEL 183 : ENDIF 184 : C 185 : C... QUA4 ou QUA8 186 : ELSE 187 : *//////// 188 : * WRITE (6,*) 'Le maillage ',IPT1, 189 : * # 'est constitué d éléments carrés' 190 : *//////// 191 : C 192 : DO 120 I = 1,NEL 193 : C 194 : IF (KCOQ(ITYP)) THEN 195 : IEL1 = IEL + (2*I-1) 196 : ELSE 197 : IEL1 = IEL + I 198 : ENDIF 199 : DO 121 IS = 1,NSGEO,NSPA 200 : LS=IS 201 : IF(ICOQ.AND.KQUAD(ITYP)) LS=(IS+1)/2 202 : IREF = (IDIM+1)*(IPT1.NUM(IS,I)-1) 203 : *//////// 204 : * WRITE (6,*) 'noeud numero',IPT1.NUM(IS,I) 205 : * WRITE (6,*) 'IREF =',IREF 206 : *//////// 207 : DO 122 K = 1,KES 208 : C** A(K,IS) = XCOOR(IREF+K) 209 : A(K,LS) = XCOOR(IREF+K) 210 : 122 CONTINUE 211 : *//////// 212 : * WRITE (6,*) 'noeud numero',IPT1.NUM(IS,I) 213 : * WRITE (6,*) 'Coordonnées :',(A(K,IS),K=1,KES) 214 : *//////// 215 : 121 CONTINUE 216 : C 217 : C On regarde comment on va découper le quadrangle 218 : DIA1 = 0. 219 : DO 123 K = 1,KES 220 : DIA1 = DIA1 + (A(K,1)-A(K,3))**2 221 : 123 CONTINUE 222 : DIA2 = 0. 223 : DO 124 K = 1,KES 224 : DIA2 = DIA2 + (A(K,2)-A(K,4))**2 225 : 124 CONTINUE 226 : *////// 227 : * WRITE(6,*) ' DIA1 DIA2 ',DIA1,DIA2 228 : *////// 229 : 230 : C stockage des connectivités des traiangles 231 : 232 : IF(ICOQ.AND.KQUAD(ITYP)) THEN 233 : 234 : * cas des qua8 235 : 236 : IF (DIA1.LE.DIA2) THEN 237 : NUK(1,KEL) = IPT1.NUM(1,I) 238 : NUK(2,KEL) = IPT1.NUM(3,I) 239 : NUK(3,KEL) = IPT1.NUM(5,I) 240 : KEL = KEL + 1 241 : IF (KCOQ(ITYP)) THEN 242 : NUK(1,KEL) = IPT1.NUM(5,I) 243 : NUK(2,KEL) = IPT1.NUM(3,I) 244 : NUK(3,KEL) = IPT1.NUM(1,I) 245 : KEL = KEL + 1 246 : ENDIF 247 : NUK(1,KEL) = IPT1.NUM(5,I) 248 : NUK(2,KEL) = IPT1.NUM(7,I) 249 : NUK(3,KEL) = IPT1.NUM(1,I) 250 : KEL = KEL + 1 251 : IF (KCOQ(ITYP)) THEN 252 : NUK(1,KEL) = IPT1.NUM(1,I) 253 : NUK(2,KEL) = IPT1.NUM(7,I) 254 : NUK(3,KEL) = IPT1.NUM(5,I) 255 : KEL = KEL + 1 256 : ENDIF 257 : ELSE 258 : NUK(1,KEL) = IPT1.NUM(3,I) 259 : NUK(2,KEL) = IPT1.NUM(5,I) 260 : NUK(3,KEL) = IPT1.NUM(7,I) 261 : KEL = KEL + 1 262 : IF (KCOQ(ITYP)) THEN 263 : NUK(1,KEL) = IPT1.NUM(7,I) 264 : NUK(2,KEL) = IPT1.NUM(5,I) 265 : NUK(3,KEL) = IPT1.NUM(3,I) 266 : KEL = KEL + 1 267 : ENDIF 268 : NUK(1,KEL) = IPT1.NUM(7,I) 269 : NUK(2,KEL) = IPT1.NUM(1,I) 270 : NUK(3,KEL) = IPT1.NUM(3,I) 271 : KEL = KEL + 1 272 : IF (KCOQ(ITYP)) THEN 273 : NUK(1,KEL) = IPT1.NUM(3,I) 274 : NUK(2,KEL) = IPT1.NUM(1,I) 275 : NUK(3,KEL) = IPT1.NUM(7,I) 276 : KEL = KEL + 1 277 : ENDIF 278 : ENDIF 279 : 280 : ELSE 281 : 282 : * cas des qua4 283 : 284 : IF (DIA1.LE.DIA2) THEN 285 : NUK(1,KEL) = IPT1.NUM(1,I) 286 : NUK(2,KEL) = IPT1.NUM(2,I) 287 : NUK(3,KEL) = IPT1.NUM(3,I) 288 : KEL = KEL + 1 289 : IF (KCOQ(ITYP)) THEN 290 : NUK(1,KEL) = IPT1.NUM(3,I) 291 : NUK(2,KEL) = IPT1.NUM(2,I) 292 : NUK(3,KEL) = IPT1.NUM(1,I) 293 : KEL = KEL + 1 294 : ENDIF 295 : NUK(1,KEL) = IPT1.NUM(3,I) 296 : NUK(2,KEL) = IPT1.NUM(4,I) 297 : NUK(3,KEL) = IPT1.NUM(1,I) 298 : KEL = KEL + 1 299 : IF (KCOQ(ITYP)) THEN 300 : NUK(1,KEL) = IPT1.NUM(1,I) 301 : NUK(2,KEL) = IPT1.NUM(4,I) 302 : NUK(3,KEL) = IPT1.NUM(3,I) 303 : KEL = KEL + 1 304 : ENDIF 305 : ELSE 306 : NUK(1,KEL) = IPT1.NUM(2,I) 307 : NUK(2,KEL) = IPT1.NUM(3,I) 308 : NUK(3,KEL) = IPT1.NUM(4,I) 309 : KEL = KEL + 1 310 : IF (KCOQ(ITYP)) THEN 311 : NUK(1,KEL) = IPT1.NUM(4,I) 312 : NUK(2,KEL) = IPT1.NUM(3,I) 313 : NUK(3,KEL) = IPT1.NUM(2,I) 314 : KEL = KEL + 1 315 : ENDIF 316 : NUK(1,KEL) = IPT1.NUM(4,I) 317 : NUK(2,KEL) = IPT1.NUM(1,I) 318 : NUK(3,KEL) = IPT1.NUM(2,I) 319 : KEL = KEL + 1 320 : IF (KCOQ(ITYP)) THEN 321 : NUK(1,KEL) = IPT1.NUM(2,I) 322 : NUK(2,KEL) = IPT1.NUM(1,I) 323 : NUK(3,KEL) = IPT1.NUM(4,I) 324 : KEL = KEL + 1 325 : ENDIF 326 : ENDIF 327 : 328 : ENDIF 329 : 330 : C 331 : DO 125 IL = 1,2 332 : C 333 : IF (KCOQ(ITYP)) THEN 334 : INDICE = KEL - IL*2 335 : ELSE 336 : INDICE = KEL - IL 337 : ENDIF 338 : KCORR(INDICE) = IEL1 339 : C 340 : DO 126 IS = 1,3 341 : IREF = (IDIM+1)*(NUK(IS,INDICE)-1) 342 : DO 127 K = 1,KES 343 : A(K,IS) = XCOOR(IREF+K) 344 : C WRITE(6,*) 'A ',K,IS,A(K,IS) 345 : 127 CONTINUE 346 : 126 CONTINUE 346 : CALL KNORM(KES,A,3,UU,S(INDICE)) 347 : DO 128 L = 1,4 348 : U(L,INDICE) = UU(L) 349 : 128 CONTINUE 350 : CC WRITE(6,*) ' KEL S U ',INDICE,S(INDICE),(U(K,INDICE),K=1,4) 351 : 352 : 125 CONTINUE 353 : C 354 : IF (KCOQ(ITYP)) THEN 355 : C On remplit KCOOR , U et S pour les éléments inverses 356 : DO 135 IL = 1,2 357 : INDICE = KEL - IL*2 +1 358 : KCORR(INDICE) = IEL1 +1 359 : DO 138 L = 1,4 360 : U(L,INDICE) = - U(L,INDICE-1) 361 : 138 CONTINUE 362 : S(INDICE)=S(INDICE-1) 363 : 135 CONTINUE 364 : ENDIF 365 : C 366 : 120 CONTINUE 367 : C 368 : IF (KCOQ(ITYP)) THEN 369 : IEL = IEL + NEL*2 370 : ELSE 371 : IEL = IEL + NEL 372 : ENDIF 373 : C 374 : C 375 : ENDIF 376 : C 377 : SEGDES IPT1 378 : C 379 : 100 CONTINUE 380 : C 381 : C On a pris toute l'information nécessaire concernant le maillage 382 : SEGDES MYMOD 383 : SEGDES INFOEL 384 : *//////// 385 : IF(KIMP.GE.2) THEN 386 : WRITE (6,*) 'Subdivision' 387 : ENDIF 388 : *//////// 389 : C----------------------------------------------------------- 390 : 391 : DO 200 K1=1,MFACE 392 : 393 : DO 202 K = 1,KES 394 : G1(K)=0. 395 : DO 201 ISS = 1,MS 396 : IREF = (IDIM+1)*(NUK(ISS,K1)-1) 397 : A1(K,ISS) = XCOOR(IREF+K) 398 : G1(K) = G1(K) + A1(K,ISS) 399 : 201 CONTINUE 400 : G1(K)=G1(K)/3. 401 : 202 CONTINUE 402 : DO 203 K=1,KES+1 403 : U1(K) = U(K,K1) 404 : 203 CONTINUE 405 : C 406 : IF (XDEC.GE.0.01) THEN 407 : DK1 = DMAX 408 : DO 400 K2 = 1,MFACE 409 : C 410 : DO 412 K=1,KES+1 411 : U2(K) = U(K,K2) 412 : 412 CONTINUE 413 : 414 : DO 212 K = 1,KES 415 : G2(K)=0. 416 : DO 211 ISS = 1,MS 417 : IREF = (IDIM+1)*(NUK(ISS,K2)-1) 418 : A2(K,ISS) = XCOOR(IREF+K) 419 : G2(K)=G2(K)+A2(K,ISS) 420 : 211 CONTINUE 421 : G2(K)=G2(K)/3. 422 : 212 CONTINUE 423 : 424 : C>>> VISIBILITE A PRIORI 425 : C 426 : 426 : CALL KPRIOR(KES,MS,MS,A1,A2,U1,U2,KVU) 427 : IF (KVU.NE.0) THEN 428 : D1=G1(1)-G2(1) 429 : D2=G1(2)-G2(2) 430 : D3=G1(3)-G2(3) 431 : DKK1 = DSQRT(D1*D1+D2*D2+D3*D3) 432 : C WRITE(6,*) ' DKK1 CK1 ',DKK1,CK1 433 : IF(DKK1.GE.DMIN) THEN 434 : CK1 = DABS(U1(1)*D1+ U1(2)*D2 + U1(3)*D3)/DKK1 435 : IF(CK1.GE.0.01) THEN 436 : DK1 = DMIN1(DK1,DKK1) 437 : ENDIF 438 : ENDIF 439 : ENDIF 440 : 400 CONTINUE 441 : DR = DK1/XDEC 442 : ELSE 443 : DR = DMAX 444 : ENDIF 445 : 446 : C 447 : 901 CONTINUE 447 : CALL KCREPA(DR,A1,S(K1),KES,MS,NPAT,GG,SS) 448 : 449 : IF(NPAT.GT.NPM) THEN 450 : DR=DR*2. 451 : 452 : IF (KIMP.GE.2) THEN 453 : WRITE(6,900) K1,NPAT 454 : 900 FORMAT(1X,' FACE ',I4,' : SUBDIVISION EXCESSIVE ',I4) 455 : ENDIF 456 : GOTO 901 457 : ENDIF 458 : NPATCH = NPAT 459 : KDECOU(K1)=NPAT 460 : MSP = KES 461 : SEGINI IPATCH 462 : 463 : DO 130 IP = 1,NPAT 464 : SP(IP) = SS(IP)/S(K1) 465 : DO 129 K=1,KES 466 : GP(K,IP) = GG(K,IP) 467 : 129 CONTINUE 468 : C WRITE(6,*) ' SP GP ',SP(IP),GP(1,IP),GP(2,IP),GP(3,IP) 469 : 130 CONTINUE 470 : 471 : KPATCH(K1) = IPATCH 472 : SEGDES IPATCH 473 : 474 : 200 CONTINUE 475 : C----------------------------------------------------------- 476 : 477 : IF(KIMP.GE.2) THEN 478 : WRITE(6,*) 479 : WRITE(6,*) 'NOMBRE DE SUBDIVISIONS PAR FACE ' 480 : WRITE(6,1000) (KDECOU(I),I=1,MFACE) 481 : 1000 FORMAT(1X,10(I4)) 482 : ENDIF 483 : C 484 : SEGDES SKFACE 485 : SEGSUP SDECOU 486 : C 487 : RETURN 488 : END 489 : 490 : 491 : 492 : 493 : 494 : 495 : 496 :
© Cast3M 2003 - All rights reserved.
Disclaimer