1 : C PROOB4 SOURCE CHAT 05/01/13 02:34:16 5004 2 : C CE PROGRAMME PERMET LA PROJECTION CONIQUE D'UN MAILLAGE 3D SUR 3 : C UNE SUFACE CARACTERISEE PAR UN AUTRE MAILLAGE 3D. 4 : 5 : C ENTREES/ MAI2 : UN MAILLAGE DE TYPE MELEME A PROJETER. 6 : C MAI1 : UN MAILLAGE REPRESENTANT LA SURFACE SUR LAQUELLE MAI2 7 : C EST PROJETE. 8 : C IP1 : UN POINT DEFINISSANT L'OEIL DE LA PROJECTION. 9 : 10 : C DATE: SEPTEMBRE 1996. 11 : 12 : C ********************************************************************** 13 : Subroutine proob4(mai2,mai1,ip1) 14 : IMPLICIT INTEGER(I-N) 15 : IMPLICIT REAL*8(A-H,O-Z) 16 : 17 : -INC CCOPTIO 18 : -INC SMELEME 19 : -INC SMCOORD 20 : -INC CCREEL 21 : *- 22 : SEGMENT ISEG1 23 : REAL*8 XLIM(2,NBEL),YLIM(2,NBEL),ZLIM(2,NBEL) 24 : ENDSEGMENT 25 : 26 : SEGMENT ISEG3 27 : INTEGER NIZO(NZO+1) 28 : ENDSEGMENT 29 : 30 : SEGMENT ISEG4 31 : INTEGER NUMZO(NZO) 32 : ENDSEGMENT 33 : 34 : SEGMENT ISEG5 35 : INTEGER NNMEL(ILON),IDEJ(NZO) 36 : ENDSEGMENT 37 : 38 : SEGMENT ISEG6 39 : REAL*8 AM(4,4),AL(4),A(4),B(4),C(4) 40 : REAL*8 XPU(2),P(3) 41 : ENDSEGMENT 42 : 43 : SEGMENT MCOOR1 44 : REAL*8 XCOOR1(XCOOR(/1)) 45 : ENDSEGMENT 46 : C LA DERNIERE COORDONNEE D'UN POINT N'EST PAS UTILISEE. 47 : 48 : SEGMENT MCOOR2 49 : REAL*8 XCOOR2(XCOOR(/1)) 50 : ENDSEGMENT 51 : C LA DERNIERE COORDONNEE D'UN POINT N'EST PAS UTILISEE. 52 : 53 : SEGMENT MCOOR3 54 : REAL*8 XCOOR3(XCOOR(/1)) 55 : ENDSEGMENT 56 : C LA DERNIERE COORDONNEE D'UN POINT N'EST PAS UTILISEE. 57 : 58 : SEGMENT MCOOR4 59 : REAL*8 XCOOR4(XCOOR(/1)) 60 : ENDSEGMENT 61 : C LA DERNIERE COORDONNEE D'UN POINT N'EST PAS UTILISEE. 62 : 63 : SEGMENT MCOR 64 : INTEGER ICOR(IONMAX+1) 65 : ENDSEGMENT 66 : 67 : SEGMENT MCORRE 68 : INTEGER ICORRE((XCOOR(/1)/(IDIM+1))) 69 : ENDSEGMENT 70 : 71 : SEGMENT ICP1 72 : INTEGER ICPR1(XCOOR(/1)/(IDIM+1)) 73 : ENDSEGMENT 74 : 75 : SEGMENT ICP2 76 : INTEGER ICPR2(XCOOR(/1)/(IDIM+1)) 77 : ENDSEGMENT 78 : 79 : SEGMENT MATRIC 80 : REAL*8 TABMAT(3,3),V(3),W(3),Q(3),invma(3,3) 81 : ENDSEGMENT 82 : 83 : SEGMENT CRIT 84 : INTEGER ELECRI(NBELEM),TABCRIT(NBELEM+1) 85 : ENDSEGMENT 86 : 87 : 88 : segini matric 89 : nbpts=xcoor(/1)/(idim+1) 90 : zchoix=100. 91 : w(1)=xcoor((ip1-1)*(idim+1)+1) 92 : w(2)=xcoor((ip1-1)*(idim+1)+2) 93 : w(3)=xcoor((ip1-1)*(idim+1)+3) 94 : 95 : C PASSAGE DU MAILLAGE MAI1 EN TRI3 95 : call change(mai1,4) 96 : 97 : C LISTAGE DES POINTS APPARTENANT A MAI1 98 : segini icp1 99 : nbpts = xcoor(/1)/( idim+1) 100 : meleme=mai1 101 : segact meleme 102 : nbelem=num(/2) 103 : nbnn = num(/1) 104 : Do 10 i=1,nbelem 105 : Do 20 j=1,nbnn 106 : ia=num(j,i) 107 : icpr1(ia)=1 108 : 20 continue 109 : 10 continue 110 : 111 : C LISTAGE DES POINTS APPARTENANT A MAI2 112 : ipt1=mai2 113 : segact ipt1 114 : segini icp2 115 : nbelem= ipt1.num(/2) 116 : nbnn = ipt1.num(/1) 117 : Do 30 i=1,nbelem 118 : Do 40 j=1,nbnn 119 : ia=ipt1.num(j,i) 120 : If (icpr1(ia).eq.1) then 121 : icpr2(ia)=2 122 : else 123 : icpr2(ia)=1 124 : Endif 125 : 40 continue 126 : 30 continue 127 : 128 : C POSITION DU CENTRE DE GRAVITE DE MAI2 129 : sx=0. 130 : sy=0. 131 : sz=0. 132 : isd=0 133 : Do 2875 ia=1,nbpts 134 : If (icpr2(ia).ne.0) then 135 : sx=sx+xcoor((ia-1)*(idim+1)+1) 136 : sy=sy+xcoor((ia-1)*(idim+1)+2) 137 : sz=sz+xcoor((ia-1)*(idim+1)+3) 138 : isd = isd+1 139 : Endif 140 : 2875 continue 141 : 142 : C CREATION DE LA MATRICE DE CHANGEMENT DE BASE 143 : v(1)=(sx/isd)-w(1) 144 : v(2)=(sy/isd)-w(2) 145 : v(3)=(sz/isd)-w(3) 146 : If ((abs(v(1))).lt.(1.E-3) .and. (abs(v(2))).lt.(1.E-3) .and. 147 : $ (abs(v(3))).lt.(1.E-3)) then 148 : v(1)=v(1)+0.1 149 : v(2)=v(2)+0.1 150 : v(3)=v(3)+0.1 151 : Endif 152 : segini mcoor3 153 : segini mcoor1 154 : segini mcoor4 155 : segini mcoor2 156 : r1=sqrt(v(1)*v(1)+v(2)*v(2)+v(3)*v(3)) 157 : v(1) = V(1) / r1 158 : v(2) = v(2) / r1 159 : v(3) = v(3) / r1 160 : r2=sqrt(v(1)*v(1)+v(2)*v(2)) 161 : if ( r2 . gt. 1.e-9) then 162 : tabmat(1,1)=-v(2)/r2 163 : tabmat(1,2)=v(1)/r2 164 : tabmat(1,3)=0. 165 : else 166 : tabmat(1,1)=1. 167 : tabmat(1,2)=0. 168 : tabmat(1,3)=0. 169 : endif 170 : tabmat(2,1)=-v(3)*tabmat(1,2) 171 : tabmat(2,2)=v(3)*tabmat(1,1) 172 : tabmat(2,3)=v(1)*tabmat(1,2)-v(2)*tabmat(1,1) 173 : tabmat(3,1)=v(1) 174 : tabmat(3,2)=v(2) 175 : tabmat(3,3)=v(3) 176 : * pas besoin de diviser par le determinant car egal à 1. 177 : INVMA(1,1)= TABMAT(2,2)*TABMAT(3,3)-TABMAT(3,2)*TABMAT(2,3) 178 : INVMA(1,2)=-TABMAT(1,2)*TABMAT(3,3)+TABMAT(3,2)*TABMAT(1,3) 179 : INVMA(1,3)= TABMAT(1,2)*TABMAT(2,3)-TABMAT(2,2)*TABMAT(1,3) 180 : INVMA(2,1)=-TABMAT(2,1)*TABMAT(3,3)+TABMAT(3,1)*TABMAT(2,3) 181 : INVMA(2,2)= TABMAT(1,1)*TABMAT(3,3)-TABMAT(3,1)*TABMAT(1,3) 182 : INVMA(2,3)=-TABMAT(1,1)*TABMAT(2,3)+TABMAT(2,1)*TABMAT(1,3) 183 : INVMA(3,1)= TABMAT(2,1)*TABMAT(3,2)-TABMAT(2,2)*TABMAT(3,1) 184 : INVMA(3,2)=-TABMAT(1,1)*TABMAT(3,2)+TABMAT(3,1)*TABMAT(1,2) 185 : INVMA(3,3)= TABMAT(1,1)*TABMAT(2,2)-TABMAT(2,1)*TABMAT(1,2) 186 : C CALCUL DES COORDONNEES DANS LE NOUVEAU REPERE POUR MELEME 187 : nbpts=xcoor(/1)/(idim+1) 188 : nbelem=num(/2) 189 : nbnn=num(/1) 190 : segini crit 191 : Do 31 i1=1,nbelem 192 : Do 38 k=1,nbnn 193 : ia=num(k,i1) 194 : xia1=xcoor((ia-1)*(idim+1)+1)-w(1) 195 : xia2=xcoor((ia-1)*(idim+1)+2)-w(2) 196 : xia3=xcoor((ia-1)*(idim+1)+3)-w(3) 197 : Do 32 j=1,3 198 : xcoor3((ia-1)*(idim+1)+j)=((xia1*tabmat(j,1))+ 199 : $(xia2*tabmat(j,2))+(xia3*tabmat(j,3))) 200 : 32 continue 201 : 38 continue 202 : 31 continue 203 : 204 : C CALCUL DE LA TAILLE LA PLUS GRANDE DES ELEMENTS SUIVANT OZ 205 : zlong=-xgrand 206 : Do 1534 i1=1,nbelem 207 : zmin=xgrand 208 : zmax=-xgrand 209 : Do 1233 i2=1,nbnn 210 : ib=num(i2,i1) 211 : ia=(ib-1)*4 212 : If (xcoor3(ia+3).lt.zmin) zmin=xcoor3(ia+3) 213 : If (xcoor3(ia+3).gt.zmax) zmax=xcoor3(ia+3) 214 : 1233 continue 215 : zlong=max(zlong,zmax-zmin) 216 : 1534 continue 217 : If (zlong.lt.(5.E-2)) zlong=(5.E-2) 218 : 219 : C DETERMINATION DES ELEMENTS TOUCHANT LA TRANCHE CRITIQUE 220 : Do 317 i1=1,nbelem 221 : ime=0 222 : Do 387 k=1,nbnn 223 : ia=num(k,i1) 224 : q(1)=(xcoor3((ia-1)*(idim+1)+1)) 225 : q(2)=(xcoor3((ia-1)*(idim+1)+2)) 226 : q(3)=(xcoor3((ia-1)*(idim+1)+3)) 227 : If (abs(xcoor3((ia-1)*(idim+1)+3)).gt.(0.55*zlong)) then 228 : xcoor1((ia-1)*(idim+1)+1)=(zchoix/(abs(q(3))))*q(1) 229 : xcoor1((ia-1)*(idim+1)+2)=(zchoix/(abs(q(3))))*q(2) 230 : xcoor1((ia-1)*(idim+1)+3)=zchoix 231 : else 232 : ime=ime+1 233 : Endif 234 : 387 continue 235 : If (ime.ne.0) then 236 : tabcrit(1)=tabcrit(1)+1 237 : tabcrit(tabcrit(1)+1)=i1 238 : elecri(i1)=1 239 : Endif 240 : 317 continue 241 : 242 : C CHANGEMENT DE REPERE POUR IPT1 243 : Do 33 ia=1,nbpts 244 : If (icpr2(ia).eq.1) then 245 : xia1=xcoor((ia-1)*(idim+1)+1)-w(1) 246 : xia2=xcoor((ia-1)*(idim+1)+2)-w(2) 247 : xia3=xcoor((ia-1)*(idim+1)+3)-w(3) 248 : Do 34 j=1,3 249 : xcoor4((ia-1)*(idim+1)+j)=((xia1*tabmat(j,1))+ 250 : $(xia2*tabmat(j,2))+(xia3*tabmat(j,3))) 251 : 34 continue 252 : If (abs(xcoor4((ia-1)*(idim+1)+3)).lt.(0.55*zlong)) then 253 : icpr2(ia)=3 254 : else 255 : q(1)=(xcoor4((ia-1)*(idim+1)+1)) 256 : q(2)=(xcoor4((ia-1)*(idim+1)+2)) 257 : q(3)=(xcoor4((ia-1)*(idim+1)+3)) 258 : xcoor2((ia-1)*(idim+1)+1)=(zchoix/(abs(q(3))))*q(1) 259 : xcoor2((ia-1)*(idim+1)+2)=(zchoix/(abs(q(3))))*q(2) 260 : xcoor2((ia-1)*(idim+1)+3)=zchoix 261 : Endif 262 : Endif 263 : If (icpr2(ia).eq.2) then 264 : xia1=xcoor((ia-1)*(idim+1)+1)-w(1) 265 : xia2=xcoor((ia-1)*(idim+1)+2)-w(2) 266 : xia3=xcoor((ia-1)*(idim+1)+3)-w(3) 267 : Do 35 j=1,3 268 : xcoor2((ia-1)*(idim+1)+j)=((xia1*tabmat(j,1))+ 269 : $(xia2*tabmat(j,2))+(xia3*tabmat(j,3))) 270 : 35 continue 271 : q(1)=(xcoor4((ia-1)*(idim+1)+1)) 272 : q(2)=(xcoor4((ia-1)*(idim+1)+2)) 273 : q(3)=(xcoor4((ia-1)*(idim+1)+3)) 274 : xcoor2((ia-1)*(idim+1)+1)=(zchoix/(abs(q(3))))*q(1) 275 : xcoor2((ia-1)*(idim+1)+2)=(zchoix/(abs(q(3))))*q(2) 276 : xcoor2((ia-1)*(idim+1)+3)=zchoix 277 : Endif 278 : 33 continue 279 : 280 : C******************** DEBUT DU ZONAGE DE MAI1 ********************** 281 : 282 : C ON CALCULE LA TAILLE MAXI D'UN ELEMENT DANS TOUTES LES DIRECTIONS 283 : C AFIN DE CREER UN ZONAGE DE L'ESPACE. EN MEME TEMPS ON CALCULE 284 : C LA DIMENSION HORS TOUT DU MAILLAGE 285 : C 286 : IDIM1=4 287 : NBEL = NUM(/2) 288 : NBNN=NUM(/1) 289 : SEGINI ISEG1 290 : ILOC=0 291 : XZO=0.D0 292 : YZO=0.D0 293 : ZZO=0.D0 294 : XZA=XGRAND 295 : YZA=XGRAND 296 : ZZA=XGRAND 297 : XTOMI=XGRAND 298 : XTOMA=-XGRAND 299 : YTOMI=XGRAND 300 : YTOMA=-XGRAND 301 : ZTOMI=XGRAND 302 : ZTOMA=-XGRAND 303 : DO 1 I1=1,NBEL 304 : If (elecri(i1).eq.0) then 305 : XMI=XGRAND 306 : YMI=XGRAND 307 : ZMI=XGRAND 308 : YMA=-XGRAND 309 : XMA=-XGRAND 310 : ZMA=-XGRAND 311 : DO 2 I2 = 1,NBNN 312 : IB=NUM(I2,I1) 313 : IA=(IB-1)*IDIM1 314 : IF(XCOOR1(IA+1).LT.XMI) XMI=XCOOR1(IA+1) 315 : IF(XCOOR1(IA+1).GT.XMA) XMA=XCOOR1(IA+1) 316 : IF(XCOOR1(IA+2).LT.YMI) YMI=XCOOR1(IA+2) 317 : IF(XCOOR1(IA+2).GT.YMA) YMA=XCOOR1(IA+2) 318 : 2 CONTINUE 319 : XLIM(1,I1)=XMI 320 : XLIM(2,I1)=XMA 321 : YLIM(1,I1)=YMI 322 : YLIM(2,I1)=YMA 323 : XZO=MAX (XZO,XMA-XMI) 324 : YZO=MAX (YZO,YMA-YMI) 325 : XZA=MIN(XZA,XMA-XMI) 326 : YZA=MIN(YZA,YMA-YMI) 327 : IF(XMI.LT.XTOMI) XTOMI=XMI 328 : IF(XMA.GT.XTOMA) XTOMA=XMA 329 : IF(YMI.LT.YTOMI) YTOMI=YMI 330 : IF(YMA.GT.YTOMA) YTOMA=YMA 331 : Endif 332 : 1 CONTINUE 333 : XPR=MIN(XZO*1.D-2,(XTOMA-XTOMI)/2.D+4) 334 : YPR=MIN(YZO*1.D-2,(YTOMA-YTOMI)/2.D+4) 335 : XZA=XZA*0.97 336 : YZA=YZA*0.97 337 : XTOMI= XTOMI - (XTOMA-XTOMI)/1.D+4 338 : XTOMA= XTOMA + (XTOMA-XTOMI)/1.D+4 339 : YTOMI= YTOMI - (YTOMA-YTOMI)/1.D+4 340 : YTOMA= YTOMA + (YTOMA-YTOMI)/1.D+4 341 : XZA=MIN ( XZA, XTOMA-XTOMI) 342 : YZA=MIN ( YZA, YTOMA-YTOMI) 343 : XZA = max ( xza , (XTOMA-XTOMI) / 50) 344 : YZA = max ( yza , (YTOMA-YTOMI) / 50) 345 : NXZO=(XTOMA-XTOMI)/XZA + 1 346 : NYZO=(YTOMA-YTOMI)/YZA + 1 347 : XZO=XZA 348 : YZO=YZA 349 : NZZO=1 350 : NXDEP=MIN(NXZO,10) 351 : NYDEP=MIN(NYZO,10) 352 : IF(FLOAT(NXZO)*FLOAT(NYZO).GT.10000.) THEN 353 : XY=SQRT(FLOAT(NXZO)*FLOAT(NYZO))/90 354 : NXZO=MAX(INT(NXZO/XY),NXDEP) 355 : NYZO=MAX(INT(NYZO/XY),NYDEP) 356 : IF(FLOAT(NXZO)*FLOAT(NYZO).GT.10000.) THEN 357 : XY=SQRT(FLOAT(NXZO)*FLOAT(NYZO))/60 358 : NXZO=MAX(INT(NXZO/XY),NXDEP) 359 : NYZO=MAX(INT(NYZO/XY),NYDEP) 360 : ENDIF 361 : XZO=(XTOMA-XTOMI)/NXZO 362 : YZO=(YTOMA-YTOMI)/NYZO 363 : NXZO=(XTOMA-XTOMI)/XZO +1 364 : NYZO=(YTOMA-YTOMI)/YZO +1 365 : ENDIF 366 : 367 : C 368 : C ON VEUT CONSTRUIRE LA LISTE DES ELEMENTS TOUCHANT UNE ZONE 369 : C POUR CELA ON COMMENCE PAR COMPTER COMBIEN D'ELEMENT TOUCHENT 370 : C CHAQUE ZONE ET EN MEME TEMPS ON STOCKE LES ZONES TOUCHEES 371 : C PAR CHAQUE ELEMENT ET LEUR NOMBRE 372 : C 373 : 374 : NZO=NXZO*NYZO 375 : IF(IIMPI.NE.0)WRITE(IOIMP,FMT='('' NZO NXZO NYZO NZZO '' 376 : $,4I7) ') NZO,NXZO,NYZO,NZZO 377 : NXYZO=NXZO*NYZO 378 : SEGINI ISEG3 379 : SEGINI ISEG4 380 : DO 3 I1=1,NBEL 381 : If (elecri(i1).eq.0) then 382 : NIZ1X=INT((XLIM(1,I1)-XTOMI-XPR)/XZO) +1 383 : NIZ1Y=INT((YLIM(1,I1)-YTOMI-YPR)/YZO) +1 384 : NIZ2X=INT((XLIM(2,I1)-XTOMI+XPR)/XZO) +1 385 : NIZ2Y=INT((YLIM(2,I1)-YTOMI+YPR)/YZO) +1 386 : DO 201 L1=NIZ1Y,NIZ2Y 387 : DO 201 L2=NIZ1X,NIZ2X 388 : NIZA = L2 + ( L1-1) * NXZO 389 : NUMZO(NIZA) = NUMZO(NIZA) +1 390 : 201 CONTINUE 391 : Endif 392 : 3 CONTINUE 393 : C 394 : C CONSTRUCTION DU TABLEAU D'ADRESSAGE DU TABLEAU DONNANT LES 395 : C ELEMENTS CONCERNEES PAR UNE ZONE 396 : C 397 : ILON=0 398 : NIZO(1)=1 399 : DO 202 L1=1,NZO 400 : NIZO(L1+1)=NIZO(L1)+NUMZO(L1) 401 : ILON=ILON+ NUMZO(L1) 402 : 202 CONTINUE 403 : 110 FORMAT(16I5) 404 : SEGINI ISEG5 405 : DO 5 I1=1,NBEL 406 : If (elecri(i1).eq.0) then 407 : NIZ1X=INT((XLIM(1,I1)-XTOMI-XPR)/XZO) +1 408 : NIZ1Y=INT((YLIM(1,I1)-YTOMI-YPR)/YZO) +1 409 : NIZ2X=INT((XLIM(2,I1)-XTOMI+XPR)/XZO) +1 410 : NIZ2Y=INT((YLIM(2,I1)-YTOMI+YPR)/YZO) +1 411 : DO 203 L1=NIZ1Y,NIZ2Y 412 : DO 203 L2=NIZ1X,NIZ2X 413 : NIZA = L2 + ( L1-1) * NXZO 414 : IAD=NIZO(NIZA)+IDEJ(NIZA) 415 : NNMEL(IAD)=I1 416 : IDEJ(NIZA)=IDEJ(NIZA)+1 417 : 203 CONTINUE 418 : Endif 419 : 5 CONTINUE 420 : C *********************** FIN DU ZONAGE DE MAI1 ************************** 421 : prec = 1.E-5 422 : prec2 = -prec 423 : prec3 = -5.E-3 424 : C ******RECHERCHE DU NOMBRE D'ELEMENT MAX CONTENU DANS UNE ZONE*********** 425 : la=0 426 : ls=0 427 : Do 420 i=1,nzo 428 : ls=ls+numzo(i) 429 : If (numzo(i).gt.la) then 430 : la=numzo(i) 431 : Endif 432 : 420 continue 433 : ionmax=la 434 : 435 : C ****************** TRAITEMENT SUR LES POINTS DE MAI2 ******************* 436 : nbelem=num(/2) 437 : nbnn=num(/1) 438 : nbpts=xcoor(/1)/(idim+1) 439 : nbpts2=nbpts 440 : segini mcor,mcorre,iseg6 441 : Do 405 ia=1,nbpts2 442 : If (icpr2(ia).eq.2) then 443 : icorre(ia)=ia 444 : Endif 445 : If (icpr2(ia).eq.1) then 446 : xpu(1)=xcoor2((ia-1)*(idim+1)+1) 447 : xpu(2)=xcoor2((ia-1)*(idim+1)+2) 448 : If (tabcrit(1).eq.0) then 449 : if (((xpu(1)-xtomi).lt.prec.) .or. 450 : $ ((xpu(2)-ytomi).lt.prec.) .or. 451 : $ ((xpu(1)-xtoma).gt.prec2.) .or. 452 : $ ((xpu(2)-ytoma).gt.prec2.)) then 452 : call erreur(21) 453 : return 454 : endif 455 : Endif 456 : C RECHERCHE DU NUMERO DE LA ZONE CORRESPONDANTE 457 : k2=int(((xcoor2((ia-1)*(idim+1)+1))-xtomi-xpr)/xzo)+1 458 : k1=int(((xcoor2((ia-1)*(idim+1)+2))-ytomi-ypr)/yzo)+1 459 : If ((k2.gt.0.) .and. (k2.lt.(nxzo+1)) .and. 460 : $ (k1.gt.0.) .and. (k1.lt.(nyzo+1))) then 461 : niza=k2+(k1-1)*nxzo 462 : C INVENTAIRE DES ELEMENTS CONCERNES PAR LE POINT IA 463 : iad=nizo(niza) 464 : Do 406 i=1,numzo(niza) 465 : i1=nnmel(iad) 466 : C CALCUL DES COORDONNEES BARYCENTRIQUES DE IA DANS L'ELEMENT I1 467 : j1=num(1,i1) 468 : j2=num(2,i1) 469 : j3=num(3,i1) 470 : j1idim=(j1-1)*(idim+1) 471 : j2idim=(j2-1)*(idim+1) 472 : j3idim=(j3-1)*(idim+1) 473 : Do 408 l=1,2 474 : p(l+1)=xpu(l) 475 : am(l+1,1)=xcoor1(j1idim+l) 476 : am(l+1,2)=xcoor1(j2idim+l) 477 : am(l+1,3)=xcoor1(j3idim+l) 478 : 408 continue 479 : x1=am(2,1) 480 : x2=am(2,2) 481 : x3=am(2,3) 482 : y1=am(3,1) 483 : y2=am(3,2) 484 : y3=am(3,3) 485 : x=p(2) 486 : y=p(3) 487 : detam=x1*y2+x2*y3+x3*y1-y1*x2-y2*x3-y3*x1 488 : a(1)=x2*y3-y2*x3 489 : a(2)=x3*y1-x1*y3 490 : a(3)=x1*y2-x2*y1 491 : b(1)=y2-y3 492 : b(2)=y3-y1 493 : b(3)=y1-y2 494 : c(1)=x3-x2 495 : c(2)=x1-x3 496 : c(3)=x2-x1 497 : Do 409 k=1,nbnn 498 : al(k)=(a(k)+b(k)*x+c(k)*y)/detam 499 : 409 continue 500 : C TEST DE POSITION SUR LES COORDONNEES BARYCENTRIQUES 501 : If (al(1).gt.prec3 . and. al(2).gt.prec3 . and. 502 : $ al(3).gt.prec3) then 503 : C LE POINT EST A L'INTERIEUR DE L'ELEMENT I1 504 : icor(1)=icor(1)+1 505 : icor(icor(1)+1)=i1 506 : Endif 507 : iad=iad+1 508 : 509 : 406 continue 510 : Endif 511 : C PASSAGE TRANCHE CRITIQUE POUR UN POINT DU ZONAGE 512 : If (k2.lt.2 .or. k2.gt.(nxzo-1)) then 513 : Do 1462 iz=1,tabcrit(1) 514 : i1=tabcrit(iz+1) 515 : m1=num(1,i1) 516 : m2=num(2,i1) 517 : m3=num(3,i1) 518 : C CALCUL DES COORDONNEES DU PROJETE DE IA SUR L'ELEMENT ICOR(2) 519 : xm1=xcoor3((m1-1)*(idim+1)+1) 520 : ym1=xcoor3((m1-1)*(idim+1)+2) 521 : zm1=xcoor3((m1-1)*(idim+1)+3) 522 : xm2=xcoor3((m2-1)*(idim+1)+1) 523 : ym2=xcoor3((m2-1)*(idim+1)+2) 524 : zm2=xcoor3((m2-1)*(idim+1)+3) 525 : xm3=xcoor3((m3-1)*(idim+1)+1) 526 : ym3=xcoor3((m3-1)*(idim+1)+2) 527 : zm3=xcoor3((m3-1)*(idim+1)+3) 528 : q(1)=(xcoor4((ia-1)*(idim+1)+1)) 529 : q(2)=(xcoor4((ia-1)*(idim+1)+2)) 530 : q(3)=(xcoor4((ia-1)*(idim+1)+3)) 531 : C Calcul du zib 532 : o1=(ym1-ym2)*(zm1-zm3)-(ym1-ym3)*(zm1-zm2) 533 : o2=(zm1-zm2)*(xm1-xm3)-(zm1-zm3)*(xm1-xm2) 534 : o3=(xm1-xm2)*(ym1-ym3)-(xm1-xm3)*(ym1-ym2) 535 : o4=-(o1*xm1+o2*ym1+o3*zm1) 536 : zk=-o4/(q(1)*o1+q(2)*o2+q(3)*o3) 537 : xnew=q(1)*zk 538 : ynew=q(2)*zk 539 : znew=q(3)*zk 540 : C PASSAGE DANS UN REPERE DONT L'AXE Z EST ORTHO AU PLAN DE I1 541 : r1=sqrt(o1*o1+o2*o2+o3*o3) 542 : if( r1 . EQ. 0. ) then 542 : call erreur ( 21 ) 543 : endif 544 : o1 = o1/r1 545 : o2 = o2/r1 546 : o3 = o3/r1 547 : r2=sqrt(o1*o1+o2*o2) 548 : if ( r2 . gt. 1.e-5) then 549 : tabmat(1,1)=-o2/r2 550 : tabmat(1,2)=o1/r2 551 : tabmat(1,3)=0. 552 : else 553 : tabmat(1,1)=1. 554 : tabmat(1,2)=0. 555 : tabmat(1,3)=0. 556 : endif 557 : tabmat(2,1)=-o3*tabmat(1,2) 558 : tabmat(2,2)=o3*tabmat(1,1) 559 : tabmat(2,3)=o1*tabmat(1,2)-o2*tabmat(1,1) 560 : tabmat(3,1)=o1 561 : tabmat(3,2)=o2 562 : tabmat(3,3)=o3 563 : x=((xnew*tabmat(1,1))+(ynew*tabmat(1,2))+(znew*tabmat(1,3))) 564 : y=((xnew*tabmat(2,1))+(ynew*tabmat(2,2))+(znew*tabmat(3,3))) 565 : x1=((xm1*tabmat(1,1))+(ym1*tabmat(1,2))+(zm1*tabmat(1,3))) 566 : y1=((xm1*tabmat(2,1))+(ym1*tabmat(2,2))+(zm1*tabmat(3,3))) 567 : x2=((xm2*tabmat(1,1))+(ym2*tabmat(1,2))+(zm2*tabmat(1,3))) 568 : y2=((xm2*tabmat(2,1))+(ym2*tabmat(2,2))+(zm2*tabmat(3,3))) 569 : x3=((xm3*tabmat(1,1))+(ym3*tabmat(1,2))+(zm3*tabmat(1,3))) 570 : y3=((xm3*tabmat(2,1))+(ym3*tabmat(2,2))+(zm3*tabmat(3,3))) 571 : C CALCUL DES COORDONNEES BARYCENTRIQUES 572 : detam=x1*y2+x2*y3+x3*y1-y1*x2-y2*x3-y3*x1 573 : a(1)=x2*y3-y2*x3 574 : a(2)=x3*y1-x1*y3 575 : a(3)=x1*y2-x2*y1 576 : b(1)=y2-y3 577 : b(2)=y3-y1 578 : b(3)=y1-y2 579 : c(1)=x3-x2 580 : c(2)=x1-x3 581 : c(3)=x2-x1 582 : Do 798 k=1,nbnn 583 : al(k)=(a(k)+b(k)*x+c(k)*y)/detam 584 : 798 continue 585 : C TEST DE POSITION SUR LES COORDONNEES BARYCENTRIQUES 586 : If (al(1).gt.prec3 . and. al(2).gt.prec3 . and. 587 : $ al(3).gt.prec3) then 588 : C LE POINT EST A L'INTERIEUR DE L'ELEMENT I1 589 : icor(1)=icor(1)+1 590 : icor(icor(1)+1)=i1 591 : Endif 592 : 1462 continue 593 : Endif 594 : 595 : If (icor(1).eq.0) then 596 : interr(1) = ia 596 : call erreur(782) 597 : return 598 : Endif 599 : 600 : Endif 601 : C FIN DU IF ICPR2(IA).EQ.1 602 : 603 : If (icpr2(ia).eq.3) then 604 : 605 : C PASSAGE SUR LES ELEMENTS DE LA TRANCHE CRITIQUE 606 : Do 1463 iz=1,tabcrit(1) 607 : i1=tabcrit(iz+1) 608 : m1=num(1,i1) 609 : m2=num(2,i1) 610 : m3=num(3,i1) 611 : C CALCUL DES COORDONNEES DU PROJETE DE IA SUR L'ELEMENT ICOR(2) 612 : xm1=xcoor3((m1-1)*(idim+1)+1) 613 : ym1=xcoor3((m1-1)*(idim+1)+2) 614 : zm1=xcoor3((m1-1)*(idim+1)+3) 615 : xm2=xcoor3((m2-1)*(idim+1)+1) 616 : ym2=xcoor3((m2-1)*(idim+1)+2) 617 : zm2=xcoor3((m2-1)*(idim+1)+3) 618 : xm3=xcoor3((m3-1)*(idim+1)+1) 619 : ym3=xcoor3((m3-1)*(idim+1)+2) 620 : zm3=xcoor3((m3-1)*(idim+1)+3) 621 : q(1)=(xcoor4((ia-1)*(idim+1)+1)) 622 : q(2)=(xcoor4((ia-1)*(idim+1)+2)) 623 : q(3)=(xcoor4((ia-1)*(idim+1)+3)) 624 : C Calcul du zib 625 : o1=(ym1-ym2)*(zm1-zm3)-(ym1-ym3)*(zm1-zm2) 626 : o2=(zm1-zm2)*(xm1-xm3)-(zm1-zm3)*(xm1-xm2) 627 : o3=(xm1-xm2)*(ym1-ym3)-(xm1-xm3)*(ym1-ym2) 628 : o4=-(o1*xm1+o2*ym1+o3*zm1) 629 : zk=-o4/(q(1)*o1+q(2)*o2+q(3)*o3) 630 : xnew=q(1)*zk 631 : ynew=q(2)*zk 632 : znew=q(3)*zk 633 : C PASSAGE DANS UN REPERE DONT L'AXE Z EST ORTHO AU PLAN DE I1 634 : r1=sqrt(o1*o1+o2*o2+o3*o3) 635 : if( r1 . EQ. 0. ) then 635 : call erreur ( 21 ) 636 : endif 637 : o1 = o1/r1 638 : o2 = o2/r1 639 : o3 = o3/r1 640 : r2=sqrt(o1*o1+o2*o2) 641 : if ( r2 . gt. 1.e-5) then 642 : tabmat(1,1)=-o2/r2 643 : tabmat(1,2)=o1/r2 644 : tabmat(1,3)=0. 645 : else 646 : tabmat(1,1)=1. 647 : tabmat(1,2)=0. 648 : tabmat(1,3)=0. 649 : endif 650 : tabmat(2,1)=-o3*tabmat(1,2) 651 : tabmat(2,2)=o3*tabmat(1,1) 652 : tabmat(2,3)=o1*tabmat(1,2)-o2*tabmat(1,1) 653 : tabmat(3,1)=o1 654 : tabmat(3,2)=o2 655 : tabmat(3,3)=o3 656 : x=((xnew*tabmat(1,1))+(ynew*tabmat(1,2))+(znew*tabmat(1,3))) 657 : y=((xnew*tabmat(2,1))+(ynew*tabmat(2,2))+(znew*tabmat(3,3))) 658 : x1=((xm1*tabmat(1,1))+(ym1*tabmat(1,2))+(zm1*tabmat(1,3))) 659 : y1=((xm1*tabmat(2,1))+(ym1*tabmat(2,2))+(zm1*tabmat(3,3))) 660 : x2=((xm2*tabmat(1,1))+(ym2*tabmat(1,2))+(zm2*tabmat(1,3))) 661 : y2=((xm2*tabmat(2,1))+(ym2*tabmat(2,2))+(zm2*tabmat(3,3))) 662 : x3=((xm3*tabmat(1,1))+(ym3*tabmat(1,2))+(zm3*tabmat(1,3))) 663 : y3=((xm3*tabmat(2,1))+(ym3*tabmat(2,2))+(zm3*tabmat(3,3))) 664 : C CALCUL DES COORDONNEES BARYCENTRIQUES 665 : detam=x1*y2+x2*y3+x3*y1-y1*x2-y2*x3-y3*x1 666 : a(1)=x2*y3-y2*x3 667 : a(2)=x3*y1-x1*y3 668 : a(3)=x1*y2-x2*y1 669 : b(1)=y2-y3 670 : b(2)=y3-y1 671 : b(3)=y1-y2 672 : c(1)=x3-x2 673 : c(2)=x1-x3 674 : c(3)=x2-x1 675 : Do 799 k=1,nbnn 676 : al(k)=(a(k)+b(k)*x+c(k)*y)/detam 677 : 799 continue 678 : C TEST DE POSITION SUR LES COORDONNEES BARYCENTRIQUES 679 : If (al(1).gt.prec3 . and. al(2).gt.prec3 . and. 680 : $al(3).gt.prec3) then 681 : C LE POINT EST A L'INTERIEUR DE L'ELEMENT I1 682 : icor(1)=icor(1)+1 683 : icor(icor(1)+1)=i1 684 : Endif 685 : 1463 continue 686 : If (icor(1).eq.0) then 687 : interr(1)=ia 687 : call erreur(782) 688 : return 689 : Endif 690 : Endif 691 : C FIN DE L'INVENTAIRE 692 : 693 : C CALCUL DES PROJETES DE MAI2 SUR MAI1 694 : If ((icpr2(ia).eq.1) .or. (icpr2(ia).eq.3)) then 695 : 696 : dmin=1.e30 697 : nbpts=nbpts+1 698 : segadj mcoord 699 : n=0 700 : Do 500 m=1,icor(1) 701 : m1=num(1,icor(m+1)) 702 : m2=num(2,icor(m+1)) 703 : m3=num(3,icor(m+1)) 704 : C Calcul des coordonnees du projete 705 : xm1=xcoor3((m1-1)*(idim+1)+1) 706 : ym1=xcoor3((m1-1)*(idim+1)+2) 707 : zm1=xcoor3((m1-1)*(idim+1)+3) 708 : xm2=xcoor3((m2-1)*(idim+1)+1) 709 : ym2=xcoor3((m2-1)*(idim+1)+2) 710 : zm2=xcoor3((m2-1)*(idim+1)+3) 711 : xm3=xcoor3((m3-1)*(idim+1)+1) 712 : ym3=xcoor3((m3-1)*(idim+1)+2) 713 : zm3=xcoor3((m3-1)*(idim+1)+3) 714 : q(1)=(xcoor4((ia-1)*(idim+1)+1)) 715 : q(2)=(xcoor4((ia-1)*(idim+1)+2)) 716 : q(3)=(xcoor4((ia-1)*(idim+1)+3)) 717 : C Calcul du zib 718 : o1=(ym1-ym2)*(zm1-zm3)-(ym1-ym3)*(zm1-zm2) 719 : o2=(zm1-zm2)*(xm1-xm3)-(zm1-zm3)*(xm1-xm2) 720 : o3=(xm1-xm2)*(ym1-ym3)-(xm1-xm3)*(ym1-ym2) 721 : o4=-(o1*xm1+o2*ym1+o3*zm1) 722 : zk=-o4/(q(1)*o1+q(2)*o2+q(3)*o3) 723 : xnew=q(1)*zk 724 : ynew=q(2)*zk 725 : znew=q(3)*zk 726 : delta1=sqrt(xnew*xnew+ynew*ynew+znew*znew) 727 : delta2=sqrt(q(1)*q(1)+q(2)* q(2)+q(3)*q(3)) 728 : pds=(xnew*q(1))+(ynew*q(2))+(znew*q(3)) 729 : If (pds.gt.0.) then 730 : If ((delta1-delta2).gt.0.) then 731 : If ((delta1-delta2).lt.dmin) then 732 : icorre(ia) = nbpts 733 : dmin=(delta1-delta2) 734 : xcoor((nbpts-1)*(idim+1)+1)=xnew 735 : xcoor((nbpts-1)*(idim+1)+2)=ynew 736 : xcoor((nbpts-1)*(idim+1)+3)=znew 737 : Endif 738 : n=n+1 739 : Endif 740 : Endif 741 : 500 continue 742 : If (n.eq.0) then 742 : call erreur(21) 743 : return 744 : Endif 745 : 746 : C *********************RETOUR A L'ANCIENNE BASE********************* 747 : ib=nbpts 748 : xa1=xcoor((ib-1)*(idim+1)+1) 749 : xa2=xcoor((ib-1)*(idim+1)+2) 750 : xa3=xcoor((ib-1)*(idim+1)+3) 751 : do 3200 j=1,3 752 : xcoor((ib-1)*(idim+1)+j)=(xa1*invma(j,1))+ 753 : $(xa2*invma(j,2))+(xa3*invma(j,3))+w(J) 754 : 3200 continue 755 : icor(1)=0 756 : n=0 757 : Endif 758 : 405 continue 759 : segsup iseg1,iseg3,iseg4,iseg5,iseg6,mcoor1,mcoor2 760 : segsup icp1,icp2,mcoor3,mcoor4 761 : C CREATION DU MAILLAGE PROJETE MAI3 762 : segini,ipt3=ipt1 763 : ipt4 = IPT3 764 : do 2456 IU = 1, max( 1,IPT3.LISOUS(/1)) 765 : If( IPT3.LISOUS(/1).NE.0) THEN 766 : IPT5=IPT3.LISOUS(IU) 767 : SEGINI,IPT4=IPT5 768 : NBREF = 0 769 : NBSOUS=0 770 : SEGADJ IPT4 771 : ENDIF 772 : nbel1=ipt4.num(/2) 773 : nbnn1=ipt4.num(/1) 774 : Do 600 i=1,nbel1 775 : Do 601 j=1,nbnn1 776 : ia=ipt4.num(j,i) 777 : ib=icorre(ia) 778 : ipt4.num(j,i)=ib 779 : 601 continue 780 : 600 continue 781 : 2456 continue 782 : do 2457 IU = 1, max( 0,IPT3.LISREF(/1)) 783 : If( IPT3.LISREF(/1).NE.0) THEN 784 : IPT5 = IPT3.LISREF(IU) 785 : SEGINI,IPT4=IPT5 786 : ENDIF 787 : nbel1=ipt4.num(/2) 788 : nbnn1=ipt4.num(/1) 789 : Do 602 i=1,nbel1 790 : Do 603 j=1,nbnn1 791 : ia=ipt4.num(j,i) 792 : ib=icorre(ia) 793 : ipt4.num(j,i)=ib 794 : 603 continue 795 : 602 continue 796 : SEGDES IPT4 797 : 2457 continue 797 : call ecrobj('MAILLAGE',ipt3) 798 : segdes ipt3 799 : segsup mcor,mcorre,matric 800 : Return 801 : End 802 : 803 : 804 : 805 : 806 : 807 : 808 : 809 : 810 :
© Cast3M 2003 - All rights reserved.
Disclaimer