1 : C ZONEG2 SOURCE CHAT 10/10/05 21:15:01 6772 2 : SUBROUTINE ZONEG2(IP1,IP2,ICOUNT,IFL) 3 : C==================================================================== 4 : C localisation de points dans des elements d une zone elementaire 5 : C zonage et construction d un tableau de correspondance 6 : C 7 : C Entrees : 8 : C ip1 pointeur sur le maillage massif 9 : C ip2 pointeur sur le maillage poi1 10 : C iexx pointeur sur un segment contenant les points deja vus 11 : C ( eviter comptage des point aux frontieres de 2 sous zones) 12 : C IFL = 0 noeuds communs a IP1 et IP2 non repertories( accro) 13 : C IFL = 1 "" comptabilisés ( proi) 14 : C Sorties : 15 : C iccoun pointeur sur le segment ICOUNT contenant le tableau de 16 : C correspondance voir commentaire pro2 17 : C==================================================================== 18 : IMPLICIT INTEGER(I-N) 19 : IMPLICIT REAL*8 (A-H,O-Z) 20 : *- 21 : -INC CCREEL 22 : -INC CCOPTIO 23 : -INC SMCOORD 24 : -INC SMELEME 25 : -INC CCGEOME 26 : SEGMENT /STRAV/(NP1(ITE),NP2(ITE),NP3(ITE),NPI(ITE)) 27 : dimension xpu(3) 28 : 29 : PARAMETER (us3=1.D0/3.D0) 30 : PARAMETER (Un=1.D0) 31 : 32 : segment inomil(0) 33 : segment wrk4 34 : real*8 xel(3,nbnn),shpp(6,nbnn),qsi(3) 35 : endsegment 36 : segment wrkl 37 : real*8 xell(3,nbnnl),shpl(6,nbnnl) 38 : endsegment 39 : 40 : 41 : SEGMENT ICOUNT 42 : INTEGER IEINT(2,NODES),IDEJVU(NODES),IACCRO 43 : REAL*8 QQQ(3,NODES),CKRIT(NODES) 44 : ENDSEGMENT 45 : EXTERNAL SHAPE 46 : 47 : ip=1 48 : ipo=1 49 : * WRITE(IOIMP,*) ' entree dans zoneg2' 50 : itr=0 51 : irate=0 52 : IDIM1 = IDIM + 1 53 : IPT2 = IP2 54 : ipt1=ip1 55 : 56 : * on fait le zonage sur les points 57 : SEGACT IPT2 58 : Nodes = IPT2.NUM(/2) 59 : * WRITE(IOIMP,*) ' nodes ' , nodes 60 : 61 : numnp=nodes 62 : ITE=numnp+3+1 63 : numnp3=ite-1 64 : * write(6,*) ' numnp ite numnp3 ' , numnp , ite , numnp3 65 : NF=100 66 : NMULT=100 67 : segini strav 68 : C 69 : C Cas 1D : idem qu'en 2D ou 3D mais pour eviter multiplications 70 : C -------- des tests (IF) qui alourdissent, on duplique ZONEG2 71 : C pour le 1D en 600 72 : IF (IDIM.EQ.1) GOTO 600 73 : 74 : C recherche des min et max 75 : C 76 : iref = ipt2.NUM(1,1)*idim1-idim 77 : XMIN=xcoor(iref) 78 : YMIN=xcoor(iref+1) 79 : ZMIN=xcoor(iref+2) 80 : YMAx=ymin 81 : XMAx=xmin 82 : ZMAx=zmin 83 : DO 201 Iyr = 2,nodes 84 : IB=ipt2.NUM(1,Iyr) 85 : IA=IB*IDIM1-idim 86 : IF(XCOOR(IA).LT.XMIN) XMIn=XCOOR(IA) 87 : IF(XCOOR(IA).GT.XMAx) XMAx=XCOOR(IA) 88 : IF(XCOOR(IA+1).LT.YMIn) YMIn=XCOOR(IA+1) 89 : IF(XCOOR(IA+1).GT.YMAx) YMAx=XCOOR(IA+1) 90 : IF( IDIM.EQ.3 ) THEN 91 : IF(XCOOR(IA+2).LT.ZMIn) ZMIn=XCOOR(IA+2) 92 : IF(XCOOR(IA+2).GT.ZMAx) ZMAx=XCOOR(IA+2) 93 : ENDIF 94 : 201 CONTINUE 95 : if(idim.eq.2) then 96 : ZMIN=0.D0 97 : ZMAX=0.D0 98 : ZDIST=1.D0 99 : endif 100 : XDIST=XMAx-XMIn 101 : YDIST=YMAx-YMIn 102 : xdistv=xdist 103 : ydistv=ydist 104 : zdistv=zdist 105 : IF(XDIST.EQ.0.D0) XDIST=1.D0 106 : IF(YDIST.EQ.0.D0) YDIST=1.D0 107 : * SMO=1E-5*XDIST*YDIST/NUMNP 108 : IF( IDIM.EQ.3 ) then 109 : zdist = zmax-zmin 110 : IF(ZDIST.EQ.0.D0) ZDIST=1.D0 111 : * SMO=1E-5*XDIST*YDIST*ZDIST/NUMNP 112 : endif 113 : C 114 : C*** recherche sur les elements de la plus petiete taille en x y et z 115 : C 116 : ipt1=ip1 117 : segact ipt1 118 : dminx=1.D10 119 : dminy=1.D10 120 : dminz=1.D10 121 : do 453 isous=1,ipt1.lisous(/1) 122 : meleme=ipt1.lisous(isous) 123 : segact meleme 124 : do 454 iel=1,num(/2) 125 : dmix=1.D10 126 : dmax=-1.D10 127 : dmiy=1.D10 128 : dmay=-1.D10 129 : dmiz=1.D10 130 : dmaz=-1.D10 131 : do 455 ip=1,num(/1) 132 : ire=num(ip,iel)*idim1-idim 133 : if(xcoor(ire).le.dmix) dmix=xcoor(ire) 134 : if(xcoor(ire).gt.dmax) dmax=xcoor(ire) 135 : if(xcoor(ire+1).le.dmiy) dmiy=xcoor(ire+1) 136 : if(xcoor(ire+1).gt.dmay) dmay=xcoor(ire+1) 137 : if(xcoor(ire+2).le.dmiz) dmiz=xcoor(ire+2) 138 : if(xcoor(ire+2).gt.dmaz) dmaz=xcoor(ire+2) 139 : 455 continue 140 : if(dminx.gt.(dmax-dmix) ) dminx= dmax-dmix 141 : if(dminy.gt.(dmay-dmiy) ) dminy= dmay-dmiy 142 : if(dminz.gt.(dmaz-dmiz) ) dminz= dmaz-dmiz 143 : 454 continue 144 : 453 continue 145 : * WRITE(IOIMP,*) 'dminx dminy dminz',dminx,dminy,dminz 146 : * WRITE(IOIMP,*) ' xdist,ydist,zdist',xdist,ydist,zdist 147 : 148 : if( dminx.eq.0.D0) then 149 : if(xdistv.eq.0.D0) then 150 : dminx=1.D0 151 : else 152 : dminx = xdistv/1000.D0 153 : endif 154 : endif 155 : if( dminy.eq.0.D0) then 156 : if(ydistv.eq.0.D0) then 157 : dminy=1.D0 158 : else 159 : dminy = ydistv/1000.D0 160 : endif 161 : endif 162 : if( dminz.eq.0.D0) then 163 : if(zdistv.eq.0.D0) then 164 : dminz=1.D0 165 : else 166 : dminz = zdistv/1000.D0 167 : endif 168 : endif 169 : dminx = max( dminx , xdist/10000.D0) 170 : dminy = max(dminy,ydist/10000.D0) 171 : if( idim.eq.3) dminz = max(dminz,zdist/10000.D0) 172 : ndex=xdist/dminx 173 : ndex=max(1,ndex) 174 : ndey=ydist/dminy 175 : ndey=max(1,ndey) 176 : if(idim.eq.3) then 177 : ndez=zdist/dminz 178 : ndez=max(1,ndez) 179 : else 180 : ndez=1 181 : zdist=1.D0 182 : dminz=1.01D0 183 : endif 184 : * WRITE(IOIMP,*) ' ndex ndey ndez',ndex,ndey,ndez 185 : 186 : C CALCUL DU DECOUPAGE EN ZONE 187 : NBZONE=NUMNP*NF 188 : * if( ndex*ndey*ndez.lt.nbzone) then 189 : xnprop = real(ndex)*real(ndey)*real(ndez) 190 : if(idim.eq.3) then 191 : rap = (real(nbzone)/xnprop)**0.333333333D0 192 : else 193 : rap= (real(nbzone)/xnprop)** 0.5D0 194 : endif 195 : if(rap.gt.1.D0) rap=1.d0 196 : ndec=ndex*rap 197 : mdec=ndey*rap 198 : if(idim.eq.3) then 199 : ldec=ndez*rap 200 : ldec = max(1,ldec) 201 : else 202 : ldec=1 203 : endif 204 : nbzone=ndec*mdec*ldec 205 : nf= nbzone / numnp3 +1 206 : nmult=nf 207 : * else 208 : * RAPY=YDIST/XDIST 209 : * NDEC=MAX(INT(SQRT(REAL(NBZONE)/RAPY)),1) 210 : * MDEC=NBZONE/NDEC 211 : * LDEC=1 212 : * IF(IDIM.EQ.3) THEN 213 : * NDEC=max(int((real(nbzone)*xdist*xdist/(ydist*zdist)) 214 : * $ **0.3333333D0),1) 215 : * MDEC= max(int((real(nbzone)*ydist*ydist/(xdist*zdist)) 216 : * $ **0.3333333D0),1) 217 : * LDEC=NBZONE/(MDEC*NDEC) 218 : * ENDIF 219 : NBZONE=NDEC*MDEC*LDEC 220 : * endif 221 : * WRITE (IOIMP,8933) NDEC,MDEC,ldec,NBZONE,nmult 222 : * 8933 FORMAT(' DECOUPAGE EN X ',I6,' EN Y ',I6,' EN Z ',i6, 223 : * $ ' SOIT ',I8,' ZONES ') 224 : * write(6,*) 'diviseur des zones pour aller aux super zones',nmult 225 : * write(6,*) ' nombre max de super zones ', nbzone/nmult+1 226 : XDEC=XDIST/NDEC*1.0000001D0 227 : YDEC=YDIST/MDEC*1.0000001D0 228 : ZDEC=ZDIST/LDEC*1.0000001D0 229 : * WRITE (IOIMP,8935) XMIN,XMAX,YMIN,YMAX,zmin, zmax,XDEC,YDEC,zdec 230 : * 8935 FORMAT (' XMIN,XMAX,YMIN,YMAX,XDEC,YDEC ',2(6G12.5)) 231 : NMDEC=NDEC*MDEC 232 : if(idim.eq.2) THEN 233 : DO 2 I=1,NUMNP 234 : IB=IPT2.num(1,i) 235 : IA=IB*IDIM1-idim 236 : IZONE=INT(real((XCOOR(IA)-XMIN)/XDEC))+1+ 237 : & INT(real((XCOOR(IA+1)-YMIN)/YDEC))*NDEC 238 : IPOINT=IZONE/NMULT+1 239 : * WRITE(IOIMP,*) ' point zone ipoint' ,ib,izone,ipoint 240 : np1(ipoint)=np1(ipoint)+1 241 : 2 CONTINUE 242 : ELSE 243 : DO 3 I=1,NUMNP 244 : IB=IPT2.num(1,i) 245 : IA=IB*IDIM1 -idim 246 : IZONE=INT(real((XCOOR(IA)-XMIn)/XDEC))+1+ 247 : & INT(real((XCOOR(IA+1)-YMIn)/YDEC))*NDEC + 248 : & INT(real((XCOOR(IA+2)-ZMIn)/ZDEC))*NMDEC 249 : IPOINT=IZONE/NMULT+1 250 : * WRITE(IOIMP,*) ' point zone ipoint' ,ib,izone,ipoint 251 : np1(ipoint)=np1(ipoint)+1 252 : 3 CONTINUE 253 : endif 254 : C np1= NB DE POINTS PAR SUPER-ZONE DE CLASSEMENT 255 : DO 4 I=2,ite 256 : np1(i)=np1(i-1)+np1(i) 257 : 4 CONTINUE 258 : * WRITE(IOIMP,*) ' np1 avant descente' 259 : * WRITE(IOIMP,*) ( np1( iyou),iyou=1,ite) 260 : if(idim.eq.2) then 261 : do 35 i=1,numnp 262 : ia = ipt2.num(1,i)*idim1-idim 263 : IZONE=INT(real((XCOOR(IA)-XMIN)/XDEC))+1+ 264 : & INT(real((XCOOR(IA+1)-YMIN)/YDEC))*NDEC 265 : IPOINT=IZONE/NMULT+1 266 : J=np1(ipoint) 267 : np3(j)=izone 268 : np2(j)=i 269 : np1(ipoint)=np1(ipoint)-1 270 : 35 continue 271 : else 272 : do 36 i=1,numnp 273 : ia = ipt2.num(1,i)*idim1-idim 274 : IZONE=INT(real((XCOOR(IA)-XMIN)/XDEC))+1+ 275 : & INT(real((XCOOR(IA+1)-YMIN)/YDEC))*NDEC + 276 : & INT(real((XCOOR(IA+2)-ZMIn)/ZDEC))*NMDEC 277 : IPOINT=IZONE/NMULT+1 278 : J=np1(ipoint) 279 : np3(j)=izone 280 : np2(j)=i 281 : np1(ipoint)=np1(ipoint)-1 282 : 36 continue 283 : endif 284 : C 285 : * WRITE(IOIMP,*) ' np1 ite avant permutation' ,ite 286 : * WRITE(IOIMP,*) ( np1( iyuo),iyuo=1,ite) 287 : * WRITE(IOIMP,*) ' np2 ite avant permutation' ,ite 288 : * WRITE(IOIMP,*)( np2( iyuo),iyuo=1,ite) 289 : * WRITE(IOIMP,*) ' np3 ite avant permutation' ,ite 290 : * WRITE(IOIMP,*)( np3( iyuo),iyuo=1,ite) 291 : DO 40 I=1,NUMNP3 292 : JD=NP1(I)+2 293 : JF=NP1(I+1) 294 : * write(6,*) ' i jd jf ' ,i,jd,jf 295 : IF (JD.LE.JF) THEN 296 : 45 IPERM=0 297 : DO 50 J=JD,JF 298 : IF(NP3(J-1).GT.NP3(J)) THEN 299 : L1=NP2(J-1) 300 : NP2(J-1)=NP2(J) 301 : NP2(J)=L1 302 : L1=NP3(J-1) 303 : NP3(J-1)=NP3(J) 304 : NP3(J)=L1 305 : IPERM=1 306 : ENDIF 307 : 50 CONTINUE 308 : IF (IPERM.EQ.1) GO TO 45 309 : ENDIF 310 : 40 CONTINUE 311 : 312 : * WRITE(IOIMP,*) ' dans zoneg2 apres boucle 40' 313 : C 314 : C a ce stade np1(i)=j le debut de la super zone i est en j-ieme position 315 : C +1 316 : C np2(j)=k le premier point dans la super zone i est le k-ieme 317 : C np3(j)=m sa vrai zone est la m-ieme 318 : C 319 : * 320 : * IL NE RESTE PLUS QU'A FAIRE LE TRAVAIL PROPREMENT DIT POUR CHAQUE element 321 : * on va boucler sur les points qui sont dans les zones en dessous 322 : * 323 : * WRITE(IOIMP,*) ' np1 ite apres peerm ',ite 324 : * WRITE(IOIMP,*) ( np1( iyuo),iyuo=1,ite) 325 : * WRITE(IOIMP,*) ' np2 ite apres peerm ',ite 326 : * WRITE(IOIMP,*)( np2( iyuo),iyuo=1,ite) 327 : * WRITE(IOIMP,*) ' np3 ite apres peerm ',ite 328 : * WRITE(IOIMP,*)( np3( iyuo),iyuo=1,ite) 329 : 330 : segact ipt1 331 : xpu(3)=0.d0 332 : do 200 isous=1,ipt1.lisous(/1) 333 : * WRITE(IOIMP,*) ' boucle 200 isous' , isous 334 : meleme = ipt1.lisous(isous) 335 : segact meleme 336 : KEL =itypel 337 : nbnn=num(/1) 338 : C 339 : kd = kdegre(kel) 340 : nso = nbsom(kel) 341 : iad=nspos(kel)-1 342 : segini wrk4 343 : segini inomil 344 : C element quaf ? si oui kell=numero lineaire correspondant 345 : C si non kell=0 345 : CALL NQF2NL(kel,kell) 346 : IF (IERR.NE.0) RETURN 347 : IF (KELL.NE.0) THEN 348 : nbnnl=NBNNE(KELL) 349 : segini wrkl 350 : do i=1,nbnnl 351 : inomil(**)=i 352 : enddo 353 : ELSE 354 : if(kd.eq.3) then 355 : C element quadratiques 356 : do 762 i=1,nbnn 357 : do 763 j=1,nso 358 : iso = ibsom(iad+j) 359 : if(i.eq.iso) goto 762 360 : 763 continue 361 : inomil(**)= i 362 : 762 continue 363 : elseif(kd.eq.2) then 364 : C elements lineaires 365 : do i=1,nbnn 366 : inomil(**)=i 367 : enddo 368 : endif 369 : ENDIF 370 : itest=1 371 : do 101 iel=1,num(/2) 371 : call doxe(xcoor,idim,nbnn,num,iel,xel) 372 : * cas quaf 373 : if (kell.ne.0) then 374 : nsom=nbsom(kel) 375 : if(nsom.ne.nbnnl) then 375 : call erreur(5) 376 : return 377 : endif 378 : idx=nspos(kel)-1 379 : do innl=1,nbnnl 380 : do iid=1,idim 381 : xell(iid,innl)=xel(iid,ibsom(idx+innl)) 382 : enddo 383 : enddo 384 : endif 385 : 386 : iref = num(1,iel)*idim1 - idim 387 : xmi=xcoor(iref) 388 : ymi=xcoor(iref+1) 389 : zmi=xcoor(iref+2) 390 : xma=xmi 391 : yma=ymi 392 : zma=zmi 393 : do 102 ipo=2,num(/1) 394 : iref=num(ipo,iel)*idim1-idim 395 : if(xcoor(iref).lt.xmi) xmi = xcoor(iref) 396 : if(xcoor(iref).gt.xma) xma = xcoor(iref) 397 : if(xcoor(iref+1).lt.ymi) ymi = xcoor(iref+1) 398 : if(xcoor(iref+1).gt.yma) yma = xcoor(iref+1) 399 : if(xcoor(iref+2).lt.zmi) zmi = xcoor(iref+2) 400 : if(xcoor(iref+2).gt.zma) zma = xcoor(iref+2) 401 : 102 continue 402 : * if(iel.lt.2.and.isous.eq.1) then 403 : * WRITE(IOIMP,*) ' element ',iel 404 : * WRITE(IOIMP,*) ( num(iuo,iel),iuo=1,num(/1)) 405 : * WRITE(IOIMP,*) ' xmi ymi zmi xma yma zma' 406 : * WRITE(IOIMP,*) ' boite reelle' 407 : * WRITE(IOIMP,*) xmi,ymi,zmi,xma,yma,zma 408 : * endif 409 : xdec1=(xma-xmi)/3.D0 410 : ydec1=(yma-ymi)/3.D0 411 : zdec1=(zma-zmi)/3.D0 412 : xmi=xmi-xdec1 413 : xma=xma+xdec1 414 : ymi=ymi-ydec1 415 : yma=yma+ydec1 416 : zmi=zmi-zdec1 417 : zma=zma+zdec1 418 : if(xmi.lt.xmin) xmi=xmin 419 : if(xmi.gt.xmax) xmi=xmax 420 : if(xma.lt.xmin) xma=xmin 421 : if(xma.gt.xmax) xma=xmax 422 : if(ymi.lt.ymin) ymi=ymin 423 : if(ymi.gt.ymax) ymi=ymax 424 : if(yma.lt.ymin) yma=ymin 425 : if(yma.gt.ymax) yma=ymax 426 : if(zmi.lt.zmin) zmi=zmin 427 : if(zmi.gt.zmax) zmi=zmax 428 : if(zma.lt.zmin) zma=zmin 429 : if(zma.gt.zmax) zma=zmax 430 : * if(iel.lt.2.and.isous.eq.1) then 431 : * 432 : * WRITE(IOIMP,*) ' boite apres modif' 433 : * WRITE(IOIMP,*) xmi,ymi,zmi,xma,yma,zma 434 : * endif 435 : * WRITE(IOIMP,*) ' xmi xma ymi yma zmi zma',xmi,xma,ymi,yma,zmi,zma 436 : * if(idim.eq.2) then 437 : * zmi=0.d0 438 : * zma=0.d0 439 : * endif 440 : C 441 : C 442 : C recherche des points interne à l'element 443 : C pour cela on balaye les super zones 444 : C 445 : if(idim.eq.2) then 446 : 447 : IXDEP=INT(real((Xmi-XMIN)/XDEC))+1+ 448 : & INT(real((ymi-YMIN)/YDEC))*NDEC 449 : 450 : 451 : ixpas=INT(real((Xma-XMIN)/XDEC))+1+ 452 : & INT(real((ymi-YMIN)/YDEC))*NDEC - ixdep 453 : 454 : iny=INT(real((yma-YMIn)/YDEC))-INT(real((ymi-YMIn) 455 : $ /YDEC))+1 456 : 457 : do 120 ily=1,iny 458 : ile=ixdep +( ily -1 )* ndec 459 : ilf=ile + ixpas 460 : izv= ile/nf+1 461 : * if (izv.gt.ite) WRITE(IOIMP,*) ' prob ',izv,ixpas,ixdep, 462 : * > ndec,mdec,xdec,ydec,ymi,ymin 463 : jd=np1(izv)+1 464 : do 122 j=jd,numnp3 465 : izd=np3(j) 466 : if(izd.lt.ile) go to 122 467 : if(izd.gt.ilf) go to 122 468 : ipo= np2(j) 469 : C recherche si le point est interne 470 : if(idejvu(ipo).eq.1) go to 122 471 : ip=ipt2.num(1,ipo) 472 : iref=ip*idim1-idim 473 : xpu(1)=xcoor(iref) 474 : xpu(2)=xcoor(iref+1) 475 : if(xpu(1).lt.xmi.or.xpu(1).gt.xma) go to 122 476 : if(xpu(2).lt.ymi.or.xpu(2).gt.yma) go to 122 477 : * cas quaf 478 : if (kell.ne.0) then 478 : call qsijs(xell,kell,nbnnl,idim,xpu, 479 : $ SHPL,qsi,iret) 479 : call qsijs2(xel,kel,nbnn,idim,xpu, 480 : $ SHPP,qsi,iret) 481 : else 481 : call qsijs(xel,kel,nbnn,idim,xpu,SHPP,qsi,iret) 482 : endif 483 : IF (IRET.NE.0) then 484 : irate=irate+1 485 : GOTO 122 486 : endif 487 : IF (KELL.NE.0) THEN 487 : CALL SHAPE(qsi(1),qsi(2),qsi(3),kell,shpl 488 : $ ,iret2) 489 : if (iret2.eq.0) then 490 : MOTERR(1:4)=NOMS(KELL) 490 : CALL ERREUR(68) 491 : RETURN 492 : endif 493 : ENDIF 494 : 495 : * WRITE(IOIMP,fmt='(8E12.5)')(shpp(1,ih),ih=1,nbnn) 496 : 497 : do i=1,inomil(/1) 498 : ilp= inomil(i) 499 : if (kell.ne.0) then 500 : if( shpl(1,ilp).lt.0.D0) goto 110 501 : else 502 : if( shpp(1,ilp).lt.0.D0) goto 110 503 : endif 504 : enddo 505 : C 506 : goto 100 507 : 110 xmap=1.D10 508 : do i=1,inomil(/1) 509 : ilp=inomil(i) 510 : if (kell.ne.0) then 511 : xmap = min ( xmap,shpl(1,ilp)) 512 : else 513 : xmap = min ( xmap,shpp(1,ilp)) 514 : endif 515 : enddo 516 : if( xmap . gt . ckrit(ipo) ) then 517 : do 180 jl =1,NUM(/1) 518 : if(IP.NE.NUM(jl,iel)) go to 180 519 : if(IFL.EQ.0) go to 190 520 : idejvu(ipo)=1 521 : itr=itr+1 522 : go to 181 523 : 180 continue 524 : 181 continue 525 : ckrit(ipo)=xmap 526 : ieint(1,ipo)=isous 527 : ieint(2,ipo)=iel 528 : QQQ(1,ipo) = qsi(1) 529 : QQQ(2,ipo) = qsi(2) 530 : endif 531 : 190 continue 532 : go to 122 533 : 100 CONTINUE 534 : 535 : C on a trouve l element contenant le point ipo attention 536 : C si le point appartient aux deux maillages on ne le stocke pas 537 : C pour accro mais on le stocke pour proi 538 : do 18 jl =1,NUM(/1) 539 : if((IP.EQ.NUM(jl,iel)).AND.(IFL.EQ.0)) GOTO 540 : $ 19 541 : 18 continue 542 : if(idejvu(ipo).EQ.0) then 543 : idejvu(ipo) = 1 544 : itr = itr + 1 545 : IEINT(1,ipo) = isous 546 : IEINT(2,ipo) = iel 547 : QQQ(1,ipo) = qsi(1) 548 : QQQ(2,ipo) = qsi(2) 549 : * WRITE(IOIMP,*) 'trouve point' ,IP, 'Element' ,J 550 : * WRITE(IOIMP,2375)ip,j,(xpu(i1),i1=1,idim),(qsi(i2),i2=1,idim) 551 : endif 552 : 19 continue 553 : * 2375 format(2i4,6e12.5) 554 : 122 continue 555 : 120 continue 556 : endif 557 : if(idim.eq.3) then 558 : 559 : 560 : 561 : IXDEP=INT(real((Xmi-XMIN)/XDEC))+1+ 562 : & INT(real((ymi-YMIN)/YDEC))*NDEC + 563 : & INT(real((zmi-zmin)/ZDEC))*NMDEC 564 : 565 : 566 : ixpas=INT(real((Xma-XMIN)/XDEC))+1+ 567 : & INT(real((ymi-YMIN)/YDEC))*NDEC + 568 : & INT(real((zmi-zmin)/ZDEC))*NMDEC - ixdep 569 : 570 : 571 : iny=INT(real((yma-YMIn)/YDEC))-INT(real((ymi-YMIn) 572 : $ /YDEC))+1 573 : inz= INT(real((zma-zmin)/ZDEC))-INT(real((zmi-zmin) 574 : $ /ZDEC))+1 575 : 576 : * if(iel.lt.2.and.isous.eq.1) WRITE(IOIMP,*) ' inz iny',inz,iny 577 : do 240 itz=1,inz 578 : idexx= ixdep + (itz-1)*nmdec 579 : do 220 ily=1,iny 580 : ile=idexx +( ily -1 )* ndec 581 : ilf=ile + ixpas 582 : * if(iel.lt.2.and.isous.eq.1)WRITE(IOIMP,*)' ile,ilf' , ile,ilf 583 : izv= ile/nf + 1 584 : jd=np1(izv)+1 585 : * if(iel.lt.2.and.isous.eq.1) WRITE(IOIMP,*) ' jd,jf ',jd,jf 586 : do 222 j=jd,numnp3 587 : izd=np3(j) 588 : * WRITE(IOIMP,*) ' izd ', izd 589 : if(izd.lt.ile) go to 222 590 : if(izd.gt.ilf) go to 223 591 : ipo= np2(j) 592 : C recherche si le point est interne 593 : if(idejvu(ipo).eq.1) go to 222 594 : ip=ipt2.num(1,ipo) 595 : * if(iel.lt.2.and.isous.eq.1)WRITE(IOIMP,*) ' ipo ip ',ipo,ip 596 : iref=ip*idim1-idim 597 : xpu(1)=xcoor(iref) 598 : xpu(2)=xcoor(iref+1) 599 : xpu(3)=xcoor(iref+2) 600 : 601 : * if(iel.lt.2.and.isous.eq.1) WRITE(IOIMP,*) 'cooor ' 602 : * $ ,xpu(1),xpu(2),xpu(3) 603 : if(xpu(1).lt.xmi.or.xpu(1).gt.xma) go to 222 604 : if(xpu(2).lt.ymi.or.xpu(2).gt.yma) go to 222 605 : if(xpu(3).lt.zmi.or.xpu(3).gt.zma) go to 222 605 : call qsijs(xel,kel,nbnn,idim,xpu,SHPP,qsi 606 : $ ,iret) 607 : IF (IRET.NE.0) THEN 608 : irate=irate+1 609 : if(irate.eq.itest) then 610 : * WRITE(IOIMP,*) 'solution pas trouvee', irate 611 : itest=itest *2 612 : endif 613 : GOTO 222 614 : endif 615 : IF (KELL.NE.0) THEN 615 : CALL SHAPE(qsi(1),qsi(2),qsi(3),kell,shpl 616 : $ ,iret2) 617 : ENDIF 618 : C 619 : * WRITE(IOIMP,fmt='(8E12.5)')(shpp(1,ih),ih=1,nbnn) 620 : do i=1,inomil(/1) 621 : ilp= inomil(i) 622 : if (kell.ne.0) then 623 : if( shpl(1,ilp).lt.0.D0) goto 210 624 : else 625 : if( shpp(1,ilp).lt.0.D0) goto 210 626 : endif 627 : enddo 628 : C 629 : goto 300 630 : 210 xmap=1.D10 631 : do i=1,inomil(/1) 632 : ilp=inomil(i) 633 : if (kell.ne.0) then 634 : xmap = min ( xmap,shpl(1,ilp)) 635 : else 636 : xmap = min ( xmap,shpp(1,ilp)) 637 : endif 638 : enddo 639 : if( xmap . gt . ckrit(ipo) ) then 640 : do 280 jl =1,NUM(/1) 641 : if(IP.NE.NUM(jl,iel)) go to 280 642 : if(IFL.EQ.0) go to 290 643 : idejvu(ipo)=1 644 : itr=itr+1 645 : * if(iel.lt.3.and.isous.eq.1) WRITE(IOIMP,*)'egalite de point' 646 : C 647 : C 648 : C 649 : go to 281 650 : 280 continue 651 : 281 continue 652 : ckrit(ipo)=xmap 653 : ieint(1,ipo)=isous 654 : ieint(2,ipo)=iel 655 : QQQ(1,ipo) = qsi(1) 656 : QQQ(2,ipo) = qsi(2) 657 : QQQ(3,ipo) = qsi(3) 658 : endif 659 : 290 continue 660 : go to 222 661 : 300 CONTINUE 662 : 663 : C on a trouve l element contenant le point ipo attention 664 : C si le point appartient aux deux maillages on ne le stocke pas 665 : C pour accro mais on le stocke pour proi 666 : do 28 jl =1,NUM(/1) 667 : if((IP.EQ.NUM(jl,iel)).AND.(IFL.EQ.0)) 668 : $ GOTO 29 669 : 28 continue 670 : if(idejvu(ipo).EQ.0) then 671 : 672 : idejvu(ipo) = 1 673 : * if(iel.lt.3.and.isous.eq.1) WRITE(IOIMP,*)'point interne' 674 : itr = itr + 1 675 : IEINT(1,ipo) = isous 676 : IEINT(2,ipo) = iel 677 : QQQ(1,ipo) = qsi(1) 678 : QQQ(2,ipo) = qsi(2) 679 : QQQ(3,ipo) = qsi(3) 680 : * if(iel.lt.3.and.isous.eq.1) then 681 : * WRITE(IOIMP,*) 'trouve point' ,IP, 'Element' ,iel 682 : * WRITE(IOIMP,2375)ip,iel,(xpu(i1),i1=1,idim),(qsi(i2),i2=1,idim 683 : C ) 684 : * endif 685 : endif 686 : 29 continue 687 : 222 continue 688 : 223 continue 689 : 220 continue 690 : 240 continue 691 : endif 692 : 101 continue 693 : segsup inomil,wrk4 694 : IF (KELL.NE.0) THEN 695 : SEGSUP wrkl 696 : ENDIF 697 : segdes meleme 698 : 200 continue 699 : segdes ipt1 700 : GOTO 800 701 : 702 : C========================== 703 : C= CAS 1D - DIMENSION 1 = 704 : C========================== 705 : 600 CONTINUE 706 : 707 : C Recherche des min et max 708 : XMIN=XGRAND 709 : XMAX=-XGRAND 710 : DO i=1,nodes 711 : j=IPT2.NUM(1,i)*IDIM1-IDIM 712 : xx=XCOOR(j) 713 : XMIN=MIN(xx,XMIN) 714 : XMAX=MAX(xx,XMAX) 715 : ENDDO 716 : XDIST=XMAX-XMIN 717 : xdistv=XDIST 718 : IF (XDIST.EQ.0.D0) XDIST=1.D0 719 : 720 : C Recherche sur les elements de la plus petite taille en x 721 : SEGACT,ipt1 722 : nbsous=ipt1.lisous(/1) 723 : dminx=XGRAND 724 : DO isous=1,nbsous 725 : meleme=ipt1.lisous(isous) 726 : SEGACT,meleme 727 : nbnn=num(/1) 728 : DO iel=1,num(/2) 729 : dmix=XGRAND 730 : dmax=-XGRAND 731 : DO i=1,nbnn 732 : j=num(i,iel)*IDIM1-IDIM 733 : xx=XCOOR(j) 734 : dmix=MIN(dmix,xx) 735 : dmax=MAX(dmax,xx) 736 : ENDDO 737 : xx=dmax-dmix 738 : dminx=MIN(dminx,xx) 739 : ENDDO 740 : ENDDO 741 : C- WRITE(IOIMP,*) 'dminx',dminx,' xdist',xdist 742 : IF (dminx.EQ.0.D0) THEN 743 : IF (xdistv.EQ.0.D0) THEN 744 : dminx=1.D0 745 : ELSE 746 : dminx=xdistv/1000.D0 747 : ENDIF 748 : ENDIF 749 : dminx=MAX(dminx,xdist/10000.D0) 750 : 751 : ndex=xdist/dminx 752 : ndex=MAX(1,ndex) 753 : C- WRITE(IOIMP,*) ' ndex=',ndex 754 : 755 : C Calcul du decoupage en zones 756 : nbzone=NUMNP*NF 757 : rap=REAL(nbzone)/ndex 758 : rap=MIN(Un,rap) 759 : ndec=ndex*rap 760 : nbzone=ndec 761 : XDEC=xdist/ndec*1.0000001D0 762 : nf=(nbzone/numnp3)+1 763 : nmult=nf 764 : C- WRITE(IOIMP,601) ndec,nbzone 765 : C- 601 FORMAT(' DECOUPAGE EN X ',I6,' SOIT ',I8,' ZONES ') 766 : C- WRITE(IOIMP,602) XMIN,XMAX,XDEC 767 : C- 602 FORMAT(' XMIN,XMAX,XDEC ',3G12.5) 768 : 769 : DO i=1,numnp 770 : j=IPT2.num(1,i)*IDIM1-IDIM 771 : izone=INT((XCOOR(j)-XMIN)/XDEC)+1 772 : ipoint=izone/nmult+1 773 : np1(ipoint)=np1(ipoint)+1 774 : C- WRITE(IOIMP,*) 'point,izone,ipoint =',i,j,izone,ipoint 775 : ENDDO 776 : C np1 = nb de points par super-zone de classement 777 : DO i=2,NUMNP3 778 : np1(i)=np1(i)+np1(i-1) 779 : ENDDO 780 : C- WRITE(IOIMP,*) ' np1 avant descente' 781 : C- WRITE(IOIMP,*) (np1(i),i=1,ite) 782 : DO i=1,numnp 783 : j=IPT2.num(1,i)*IDIM1-IDIM 784 : izone=INT((XCOOR(j)-XMIN)/XDEC)+1 785 : ipoint=izone/nmult+1 786 : j=np1(ipoint) 787 : np3(j)=izone 788 : np2(j)=i 789 : np1(ipoint)=j-1 790 : ENDDO 791 : 792 : C- WRITE(IOIMP,*) ' np1 ite avant permutation' ,ite 793 : C- WRITE(IOIMP,*) (np1(i),i=1,ite) 794 : C- WRITE(IOIMP,*) (np2(i),i=1,ite) 795 : C- WRITE(IOIMP,*) (np3(i),i=1,ite) 796 : DO i=1,NUMNP 797 : jD=np1(i)+2 798 : jF=np1(i+1) 799 : IF (jD.LE.jF) THEN 800 : 610 IPERM=0 801 : DO j=jD,jF 802 : k=j-1 803 : IF (np3(k).GT.np3(j)) THEN 804 : L1=np2(k) 805 : np2(k)=np2(j) 806 : np2(j)=L1 807 : L1=np3(k) 808 : np3(k)=np3(j) 809 : np3(j)=L1 810 : IPERM=1 811 : ENDIF 812 : ENDDO 813 : IF (IPERM.EQ.1) GOTO 610 814 : ENDIF 815 : ENDDO 816 : C- WRITE(IOIMP,*) ' dans zoneg2 apres permutation' 817 : C- WRITE(IOIMP,*) ' np1 ite apres perm ',ite 818 : C- WRITE(IOIMP,*) (np1(i),i=1,ite) 819 : C- WRITE(IOIMP,*) (np2(i),i=1,ite) 820 : C- WRITE(IOIMP,*) (np3(i),i=1,ite) 821 : 822 : C A ce stade : 823 : C - np1(i)=j le debut de la super-zone i est en j-eme position + 1 824 : C - np2(j)=k le premier point dans la super zone i est le k-eme 825 : C - np3(j)=m sa vrai zone est la m-eme 826 : C Il ne reste plus qu'a faire le travail pour chaque element. 827 : C On va boucler sur les points qui sont dans les zones en dessous. 828 : 829 : xpu(2)=0.D0 830 : xpu(3)=0.D0 831 : C SEGACT,ipt1 832 : DO isous=1,nbsous 833 : C- WRITE(IOIMP,*) ' boucle isous=',isous 834 : meleme=ipt1.lisous(isous) 835 : C SEGACT,meleme 836 : nbnn=num(/1) 837 : 838 : SEGINI,wrk4 839 : SEGINI,inomil 840 : 841 : kel=ITYPEL 842 : kd=KDEGRE(kel) 843 : C Elements quadratiques 844 : IF (kd.EQ.3) then 845 : jd=NSPOS(kel) 846 : jf=jd+NBSOM(kel)-1 847 : DO i=1,nbnn 848 : DO j=jd,jf 849 : IF (i.EQ.ibsom(j)) GOTO 620 850 : ENDDO 851 : inomil(**)=i 852 : 620 CONTINUE 853 : ENDDO 854 : C Elements lineaires 855 : ELSE IF (kd.EQ.2) THEN 856 : DO i=1,nbnn 857 : inomil(**)=i 858 : ENDDO 859 : ENDIF 860 : 861 : DO iel=1,num(/2) 861 : CALL DOXE(XCOOR,IDIM,nbnn,num,iel,xel) 862 : j=num(1,iel)*IDIM1-IDIM 863 : xmi=XCOOR(j) 864 : xma=xmi 865 : DO i=2,nbnn 866 : j=num(i,iel)*IDIM1-IDIM 867 : xx=XCOOR(j) 868 : xmi=MIN(xmi,xx) 869 : xma=MAX(xma,xx) 870 : ENDDO 871 : C- IF (iel.EQ.1.AND.isous.EQ.1) THEN 872 : C- WRITE(IOIMP,*) ' element ',iel 873 : C- WRITE(IOIMP,*) (num(i,iel),i=1,nbnn) 874 : C- WRITE(IOIMP,*) ' boite reelle : xmi xma ',xmi,xma 875 : C- ENDIF 876 : xdec1=us3*(xma-xmi) 877 : xmi=xmi-xdec1 878 : xma=xma+xdec1 879 : xmi=MAX(xmi,xmin) 880 : xmi=MIN(xmax,xmi) 881 : xma=MAX(xma,xmin) 882 : xma=MIN(xma,xmax) 883 : C- IF (iel.EQ.1.AND.isous.EQ.1) THEN 884 : C- WRITE(IOIMP,*) ' element ',iel 885 : C- WRITE(IOIMP,*) ' boite apres modif : xmi xma ',xmi,xma 886 : C- ENDIF 887 : C- WRITE(IOIMP,*) ' xmi xma :',xmi,xma 888 : 889 : C Recherche des points internes a l'element par balayage des super zones 890 : ixdep=INT((xmi-xmin)/xdec)+1 891 : ixpas=INT((xma-xmin)/xdec)+1 892 : ile=ixdep 893 : ilf=ile+ixpas 894 : izv=ile/nf+1 895 : C- IF (izv.GT.ite) WRITE(IOIMP,*) ' prob ',izv,ixpas,ixdep,ndec,xdec 896 : jd=np1(izv)+1 897 : DO 630 j=jd,numnp3 898 : izd=np3(j) 899 : IF (izd.LT.ile.OR.izd.GT.ilf) GOTO 630 900 : C Recherche si le point est interne 901 : ipo=np2(j) 902 : IF (idejvu(ipo).EQ.1) GOTO 630 903 : ip=ipt2.num(1,ipo) 904 : xpu(1)=XCOOR(ip*IDIM1-IDIM) 905 : IF (xpu(1).LT.xmi.OR.xpu(1).GT.xma) GOTO 630 905 : CALL QSIJS(xel,kel,nbnn,IDIM,xpu,shpp,qsi,IRET) 906 : C- WRITE(IOIMP,FMT='(8E12.5)') (shpp(1,i),i=1,nbnn) 907 : IF (IRET.NE.0) then 908 : irate=irate+1 909 : GOTO 630 910 : ENDIF 911 : 912 : do i=1,inomil(/1) 913 : ilp= inomil(i) 914 : if( shpp(1,ilp).lt.0.D0) goto 1110 915 : enddo 916 : C 917 : goto 1100 918 : 1110 xmap=1.D10 919 : do i=1,inomil(/1) 920 : ilp=inomil(i) 921 : xmap = min ( xmap,shpp(1,ilp)) 922 : enddo 923 : if( xmap . gt . ckrit(ipo) ) then 924 : do 1180 jl =1,NUM(/1) 925 : if(IP.NE.NUM(jl,iel)) go to 1180 926 : if(IFL.EQ.0) go to 1190 927 : idejvu(ipo)=1 928 : itr=itr+1 929 : go to 1181 930 : 1180 continue 931 : 1181 continue 932 : ckrit(ipo)=xmap 933 : ieint(1,ipo)=isous 934 : ieint(2,ipo)=iel 935 : QQQ(1,ipo) = qsi(1) 936 : endif 937 : 1190 continue 938 : go to 630 939 : 1100 CONTINUE 940 : 941 : 942 : 943 : 944 : C on a trouve l element contenant le point ipo attention 945 : C si le point appartient aux deux maillages on ne le stocke pas 946 : C pour accro mais on le stocke pour proi 947 : DO i=1,nbnn 948 : IF ((ip.EQ.NUM(i,iel)).AND.(IFL.EQ.0)) GOTO 630 949 : ENDDO 950 : IF (idejvu(ipo).EQ.0) THEN 951 : idejvu(ipo)=1 952 : itr=itr+1 953 : IEINT(1,ipo)=isous 954 : IEINT(2,ipo)=iel 955 : QQQ(1,ipo)=qsi(1) 956 : C WRITE(IOIMP,*) 'Trouve point de numero global ',IP 957 : C WRITE(IOIMP,*) 'dans le ',iel,'eme élément ', 958 : C $ 'de la ',isous,'eme sous-zone' 959 : C WRITE(IOIMP,*) 'Coordonnée globale = ',xpu(1) 960 : C WRITE(IOIMP,*) 'Coordonnée locale = ',qsi(1) 961 : CC- WRITE(IOIMP,*) 'Trouve point' ,IP, 'Element' ,J 962 : CC- WRITE(IOIMP,603) ip,j,xpu(1),qsi(1) 963 : C 603 FORMAT(2I4,2E12.5) 964 : ENDIF 965 : 630 CONTINUE 966 : ENDDO 967 : SEGSUP,inomil,wrk4 968 : SEGDES,meleme 969 : ENDDO 970 : SEGDES IPT1 971 : 800 CONTINUE 972 : IACCRO=IACCRO+itr 973 : C- WRITE(IOIMP,*) ' nb points non converges dans qsijs ',irate 974 : SEGDES,IPT2 975 : SEGSUP STRAV 976 : RETURN 977 : END 978 : 979 : 980 : 981 : 982 : 983 : 984 : 985 : 986 : 987 : 988 : 989 :
© Cast3M 2003 - All rights reserved.
Disclaimer