C ZONEG2 SOURCE OF166741 25/07/28 21:15:16 12336 SUBROUTINE ZONEG2(IP1,IP2,ICOUNT,IFL) C==================================================================== C localisation de points dans des elements d une zone elementaire C zonage et construction d un tableau de correspondance C C Entrees : C ip1 pointeur sur le maillage massif C ip2 pointeur sur le maillage poi1 C iexx pointeur sur un segment contenant les points deja vus C ( eviter comptage des point aux frontieres de 2 sous zones) C IFL = 0 noeuds communs a IP1 et IP2 non repertories( accro) C IFL = 1 "" comptabilisés ( proi) C Sorties : C iccoun pointeur sur le segment ICOUNT contenant le tableau de C correspondance voir commentaire pro2 C==================================================================== IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) -INC PPARAM -INC CCOPTIO -INC CCREEL -INC CCGEOME -INC SMCOORD -INC SMELEME SEGMENT /STRAV/(NP1(ITE),NP2(ITE),NP3(ITE),NPI(ITE)) dimension xpu(3) PARAMETER (us3=1.D0/3.D0) PARAMETER (Un=1.D0) segment inomil(0) segment wrk4 real*8 xel(3,nbnn),shpp(6,nbnn),qsi(3) endsegment segment wrkl real*8 xell(3,nbnnl),shpl(6,nbnnl) endsegment SEGMENT ICOUNT INTEGER IEINT(2,NODES),IDEJVU(NODES),IACCRO REAL*8 QQQ(3,NODES),CKRIT(NODES) ENDSEGMENT EXTERNAL SHAPE REAL*8 ZDIST ZDIST=0.D0 ip=1 ipo=1 * WRITE(IOIMP,*) ' entree dans zoneg2' itr=0 irate=0 IDIM1 = IDIM + 1 IPT2 = IP2 ipt1=ip1 * on fait le zonage sur les points SEGACT IPT2 Nodes = IPT2.NUM(/2) * WRITE(IOIMP,*) ' nodes ' , nodes numnp=nodes ITE=numnp+3+1 numnp3=ite-1 * write(6,*) ' numnp ite numnp3 ' , numnp , ite , numnp3 NF=100 NMULT=100 segini strav C C Cas 1D : idem qu'en 2D ou 3D mais pour eviter multiplications C -------- des tests (IF) qui alourdissent, on duplique ZONEG2 C pour le 1D en 600 IF (IDIM.EQ.1) GOTO 600 C recherche des min et max C iref = ipt2.NUM(1,1)*idim1-idim XMIN=xcoor(iref) YMIN=xcoor(iref+1) ZMIN=xcoor(iref+2) YMAx=ymin XMAx=xmin ZMAx=zmin DO 201 Iyr = 2,nodes IB=ipt2.NUM(1,Iyr) IA=IB*IDIM1-idim IF(XCOOR(IA).LT.XMIN) XMIn=XCOOR(IA) IF(XCOOR(IA).GT.XMAx) XMAx=XCOOR(IA) IF(XCOOR(IA+1).LT.YMIn) YMIn=XCOOR(IA+1) IF(XCOOR(IA+1).GT.YMAx) YMAx=XCOOR(IA+1) IF( IDIM.EQ.3 ) THEN IF(XCOOR(IA+2).LT.ZMIn) ZMIn=XCOOR(IA+2) IF(XCOOR(IA+2).GT.ZMAx) ZMAx=XCOOR(IA+2) ENDIF 201 CONTINUE if(idim.eq.2) then ZMIN=0.D0 ZMAX=0.D0 ZDIST=1.D0 ZDISTV=1.D0 endif XDIST=XMAx-XMIn YDIST=YMAx-YMIn xdistv=xdist ydistv=ydist IF(XDIST.EQ.0.D0) THEN NDEX=1 XDIST=1.D0 ELSE NDEX=0 ENDIF IF(YDIST.EQ.0.D0) THEN NDEY=1 YDIST=1.D0 ELSE NDEY=0 ENDIF * SMO=1E-5*XDIST*YDIST/NUMNP IF( IDIM.EQ.3 ) then zdist = zmax-zmin zdistv=zdist IF(ZDIST.EQ.0.D0) THEN NDEZ=1 ZDIST=1.D0 * SMO=1E-5*XDIST*YDIST*ZDIST/NUMNP ELSE NDEZ=0 ENDIF endif * write(ioimp,*) 'Maillage POI1' * write(ioimp,*) 'xmin,xmax,ymin,ymax,zmin,zmax=',xmin,xmax,ymin * $ ,ymax,zmin,zmax * write(ioimp,*) 'xdistv,ydistv,zdistv=',xdistv,ydistv,zdistv * write(ioimp,*) 'xdist,ydist,zdist=',xdist,ydist,zdist C C*** recherche sur les elements de la plus petite taille en x y et z C ipt1=ip1 segact ipt1 dminx=1.D10 dminy=1.D10 dminz=1.D10 do 453 isous=1,ipt1.lisous(/1) meleme=ipt1.lisous(isous) segact meleme do 454 iel=1,num(/2) dmix=1.D10 dmax=-1.D10 dmiy=1.D10 dmay=-1.D10 dmiz=1.D10 dmaz=-1.D10 do 455 ip=1,num(/1) ire=num(ip,iel)*idim1-idim if(xcoor(ire).le.dmix) dmix=xcoor(ire) if(xcoor(ire).gt.dmax) dmax=xcoor(ire) if(xcoor(ire+1).le.dmiy) dmiy=xcoor(ire+1) if(xcoor(ire+1).gt.dmay) dmay=xcoor(ire+1) if(xcoor(ire+2).le.dmiz) dmiz=xcoor(ire+2) if(xcoor(ire+2).gt.dmaz) dmaz=xcoor(ire+2) 455 continue if(dminx.gt.(dmax-dmix) ) dminx= dmax-dmix if(dminy.gt.(dmay-dmiy) ) dminy= dmay-dmiy if(dminz.gt.(dmaz-dmiz) ) dminz= dmaz-dmiz 454 continue 453 continue * write(ioimp,*) 'Maillage Massif' * WRITE(IOIMP,*) 'dminx,dminy,dminz=',dminx,dminy,dminz * WRITE(IOIMP,*) ' xdist,ydist,zdist',xdist,ydist,zdist if( dminx.eq.0.D0) then if(xdistv.eq.0.D0) then dminx=1.D0 else dminx = xdistv/1000.D0 endif endif if( dminy.eq.0.D0) then if(ydistv.eq.0.D0) then dminy=1.D0 else dminy = ydistv/1000.D0 endif endif if (idim.eq.3) then if( dminz.eq.0.D0) then if(zdistv.eq.0.D0) then dminz=1.D0 else dminz = zdistv/1000.D0 endif endif endif * * write(ioimp,*) 'Apres correction 1' * WRITE(IOIMP,*) 'dminx,dminy,dminz=',dminx,dminy,dminz dminx = max( dminx , xdist/10000.D0) dminy = max(dminy,ydist/10000.D0) if( idim.eq.3) dminz = max(dminz,zdist/10000.D0) * write(ioimp,*) 'Apres correction 2' * WRITE(IOIMP,*) 'dminx,dminy,dminz=',dminx,dminy,dminz if (ndex.eq.0) ndex=xdist/dminx ndex=max(1,ndex) if (ndey.eq.0) ndey=ydist/dminy ndey=max(1,ndey) if(idim.eq.3) then if (ndez.eq.0) ndez=zdist/dminz ndez=max(1,ndez) else ndez=1 zdist=1.D0 dminz=1.01D0 endif * WRITE(IOIMP,*) ' ndex ndey ndez',ndex,ndey,ndez C CALCUL DU DECOUPAGE EN ZONE NBZONE=NUMNP*NF * if( ndex*ndey*ndez.lt.nbzone) then * WRITE(IOIMP,*) ' ndex ndey ndez nbzone=',ndex,ndey,ndez,nbzone xnprop = real(ndex)*real(ndey)*real(ndez) * dimension effective : là où les découpages sont plus grands que 1 idimf=0 if (ndex.gt.1) idimf=idimf+1 if (ndey.gt.1) idimf=idimf+1 if (ndez.gt.1) idimf=idimf+1 * write(ioimp,*) ' idimf=',idimf * if(idim.eq.3) then * rap = (real(nbzone)/xnprop)**0.333333333D0 * else * rap= (real(nbzone)/xnprop)** 0.5D0 * endif if (idimf.gt.0) then rap = (real(nbzone)/xnprop)**(1.D0/real(idimf)) * write(ioimp,*) ' rap=',rap if(rap.gt.1.D0) rap=1.d0 else rap=1.d0 endif * write(ioimp,*) ' rap=',rap if (ndex.gt.1) then ndec=ndex*rap ndec=max(1,ndec) else ndec=1 endif * if (ndey.gt.1) then mdec=ndey*rap mdec = max(1,mdec) else mdec=1 endif * if(idim.eq.3) then if (ndez.gt.1) then ldec=ndez*rap ldec = max(1,ldec) else ldec=1 endif else ldec=1 endif * nbzone=ndec*mdec*ldec nf= nbzone / numnp3 +1 nmult=nf * write(ioimp,*) 'Apres correction 3' * WRITE(IOIMP,*) ' ndec mdec ldec nbzone=',ndec,mdec,ldec,nbzone * else * RAPY=YDIST/XDIST * NDEC=MAX(INT(SQRT(REAL(NBZONE)/RAPY)),1) * MDEC=NBZONE/NDEC * LDEC=1 * IF(IDIM.EQ.3) THEN * NDEC=max(int((real(nbzone)*xdist*xdist/(ydist*zdist)) * $ **0.3333333D0),1) * MDEC= max(int((real(nbzone)*ydist*ydist/(xdist*zdist)) * $ **0.3333333D0),1) * LDEC=NBZONE/(MDEC*NDEC) * ENDIF NBZONE=NDEC*MDEC*LDEC * endif * WRITE (IOIMP,8933) NDEC,MDEC,ldec,NBZONE,nmult * 8933 FORMAT(' DECOUPAGE EN X ',I6,' EN Y ',I6,' EN Z ',i6, * $ ' SOIT ',I8,' ZONES ') * write(6,*) 'diviseur des zones pour aller aux super zones',nmult * write(6,*) ' nombre max de super zones ', nbzone/nmult+1 XDEC=XDIST/NDEC*1.0000001D0 YDEC=YDIST/MDEC*1.0000001D0 ZDEC=ZDIST/LDEC*1.0000001D0 * WRITE (IOIMP,8935) XMIN,XMAX,YMIN,YMAX,zmin, zmax,XDEC,YDEC,zdec * 8935 FORMAT (' XMIN,XMAX,YMIN,YMAX,XDEC,YDEC ',2(6G12.5)) NMDEC=NDEC*MDEC if(idim.eq.2) THEN DO 2 I=1,NUMNP IB=IPT2.num(1,i) IA=IB*IDIM1-idim IZONE=INT(real((XCOOR(IA)-XMIN)/XDEC))+1+ & INT(real((XCOOR(IA+1)-YMIN)/YDEC))*NDEC IPOINT=IZONE/NMULT+1 * WRITE(IOIMP,*) ' point zone ipoint' ,ib,izone,ipoint np1(ipoint)=np1(ipoint)+1 2 CONTINUE ELSE DO 3 I=1,NUMNP IB=IPT2.num(1,i) IA=IB*IDIM1 -idim IZONE=INT(real((XCOOR(IA)-XMIn)/XDEC))+1+ & INT(real((XCOOR(IA+1)-YMIn)/YDEC))*NDEC + & INT(real((XCOOR(IA+2)-ZMIn)/ZDEC))*NMDEC IPOINT=IZONE/NMULT+1 * WRITE(IOIMP,*) ' point zone ipoint' ,ib,izone,ipoint np1(ipoint)=np1(ipoint)+1 3 CONTINUE endif C np1= NB DE POINTS PAR SUPER-ZONE DE CLASSEMENT DO 4 I=2,ite np1(i)=np1(i-1)+np1(i) 4 CONTINUE * WRITE(IOIMP,*) ' np1 avant descente' * WRITE(IOIMP,*) ( np1( iyou),iyou=1,ite) if(idim.eq.2) then do 35 i=1,numnp ia = ipt2.num(1,i)*idim1-idim IZONE=INT(real((XCOOR(IA)-XMIN)/XDEC))+1+ & INT(real((XCOOR(IA+1)-YMIN)/YDEC))*NDEC IPOINT=IZONE/NMULT+1 J=np1(ipoint) np3(j)=izone np2(j)=i np1(ipoint)=np1(ipoint)-1 35 continue else do 36 i=1,numnp ia = ipt2.num(1,i)*idim1-idim IZONE=INT(real((XCOOR(IA)-XMIN)/XDEC))+1+ & INT(real((XCOOR(IA+1)-YMIN)/YDEC))*NDEC + & INT(real((XCOOR(IA+2)-ZMIn)/ZDEC))*NMDEC IPOINT=IZONE/NMULT+1 J=np1(ipoint) np3(j)=izone np2(j)=i np1(ipoint)=np1(ipoint)-1 36 continue endif C * WRITE(IOIMP,*) ' np1 ite avant permutation' ,ite * WRITE(IOIMP,*) ( np1( iyuo),iyuo=1,ite) * WRITE(IOIMP,*) ' np2 ite avant permutation' ,ite * WRITE(IOIMP,*)( np2( iyuo),iyuo=1,ite) * WRITE(IOIMP,*) ' np3 ite avant permutation' ,ite * WRITE(IOIMP,*)( np3( iyuo),iyuo=1,ite) DO 40 I=1,NUMNP3 JD=NP1(I)+2 JF=NP1(I+1) * write(6,*) ' i jd jf ' ,i,jd,jf IF (JD.LE.JF) THEN 45 IPERM=0 DO 50 J=JD,JF IF(NP3(J-1).GT.NP3(J)) THEN L1=NP2(J-1) NP2(J-1)=NP2(J) NP2(J)=L1 L1=NP3(J-1) NP3(J-1)=NP3(J) NP3(J)=L1 IPERM=1 ENDIF 50 CONTINUE IF (IPERM.EQ.1) GO TO 45 ENDIF 40 CONTINUE * WRITE(IOIMP,*) ' dans zoneg2 apres boucle 40' C C a ce stade np1(i)=j le debut de la super zone i est en j-ieme position C +1 C np2(j)=k le premier point dans la super zone i est le k-ieme C np3(j)=m sa vrai zone est la m-ieme C * * IL NE RESTE PLUS QU'A FAIRE LE TRAVAIL PROPREMENT DIT POUR CHAQUE element * on va boucler sur les points qui sont dans les zones en dessous * * WRITE(IOIMP,*) ' np1 ite apres peerm ',ite * WRITE(IOIMP,*) ( np1( iyuo),iyuo=1,ite) * WRITE(IOIMP,*) ' np2 ite apres peerm ',ite * WRITE(IOIMP,*)( np2( iyuo),iyuo=1,ite) * WRITE(IOIMP,*) ' np3 ite apres peerm ',ite * WRITE(IOIMP,*)( np3( iyuo),iyuo=1,ite) segact ipt1 xpu(3)=0.d0 do 200 isous=1,ipt1.lisous(/1) * WRITE(IOIMP,*) ' boucle 200 isous' , isous meleme = ipt1.lisous(isous) segact meleme KEL =itypel nbnn=num(/1) C kd = kdegre(kel) nso = nbsom(kel) iad=nspos(kel)-1 segini wrk4 segini inomil C element quaf ? si oui kell=numero lineaire correspondant C si non kell=0 CALL NQF2NL(kel,kell) IF (IERR.NE.0) RETURN IF (KELL.NE.0) THEN nbnnl=NBNNE(KELL) segini wrkl do i=1,nbnnl inomil(**)=i enddo ELSE if(kd.eq.3) then C element quadratiques do 762 i=1,nbnn do 763 j=1,nso iso = ibsom(iad+j) if(i.eq.iso) goto 762 763 continue inomil(**)= i 762 continue elseif(kd.eq.2) then C elements lineaires do i=1,nbnn inomil(**)=i enddo endif ENDIF itest=1 do 101 iel=1,num(/2) call doxe(xcoor,idim,nbnn,num,iel,xel) * cas quaf if (kell.ne.0) then nsom=nbsom(kel) if(nsom.ne.nbnnl) then call erreur(5) return endif idx=nspos(kel)-1 do innl=1,nbnnl do iid=1,idim xell(iid,innl)=xel(iid,ibsom(idx+innl)) enddo enddo endif iref = num(1,iel)*idim1 - idim xmi=xcoor(iref) ymi=xcoor(iref+1) zmi=xcoor(iref+2) xma=xmi yma=ymi zma=zmi do 102 ipo=2,num(/1) iref=num(ipo,iel)*idim1-idim if(xcoor(iref).lt.xmi) xmi = xcoor(iref) if(xcoor(iref).gt.xma) xma = xcoor(iref) if(xcoor(iref+1).lt.ymi) ymi = xcoor(iref+1) if(xcoor(iref+1).gt.yma) yma = xcoor(iref+1) if(xcoor(iref+2).lt.zmi) zmi = xcoor(iref+2) if(xcoor(iref+2).gt.zma) zma = xcoor(iref+2) 102 continue * if(iel.lt.2.and.isous.eq.1) then * WRITE(IOIMP,*) ' element ',iel * WRITE(IOIMP,*) ( num(iuo,iel),iuo=1,num(/1)) * WRITE(IOIMP,*) ' xmi ymi zmi xma yma zma' * WRITE(IOIMP,*) ' boite reelle' * WRITE(IOIMP,*) xmi,ymi,zmi,xma,yma,zma * endif xdec1=(xma-xmi)/3.D0 ydec1=(yma-ymi)/3.D0 zdec1=(zma-zmi)/3.D0 xmi=xmi-xdec1 xma=xma+xdec1 ymi=ymi-ydec1 yma=yma+ydec1 zmi=zmi-zdec1 zma=zma+zdec1 if(xmi.lt.xmin) xmi=xmin if(xmi.gt.xmax) xmi=xmax if(xma.lt.xmin) xma=xmin if(xma.gt.xmax) xma=xmax if(ymi.lt.ymin) ymi=ymin if(ymi.gt.ymax) ymi=ymax if(yma.lt.ymin) yma=ymin if(yma.gt.ymax) yma=ymax if(zmi.lt.zmin) zmi=zmin if(zmi.gt.zmax) zmi=zmax if(zma.lt.zmin) zma=zmin if(zma.gt.zmax) zma=zmax * if(iel.lt.2.and.isous.eq.1) then * * WRITE(IOIMP,*) ' boite apres modif' * WRITE(IOIMP,*) xmi,ymi,zmi,xma,yma,zma * endif * WRITE(IOIMP,*) ' xmi xma ymi yma zmi zma',xmi,xma,ymi,yma,zmi,zma * if(idim.eq.2) then * zmi=0.d0 * zma=0.d0 * endif C C C recherche des points interne à l'element C pour cela on balaye les super zones C if(idim.eq.2) then IXDEP=INT(real((Xmi-XMIN)/XDEC))+1+ & INT(real((ymi-YMIN)/YDEC))*NDEC ixpas=INT(real((Xma-XMIN)/XDEC))+1+ & INT(real((ymi-YMIN)/YDEC))*NDEC - ixdep iny=INT(real((yma-YMIn)/YDEC))-INT(real((ymi-YMIn) $ /YDEC))+1 do 120 ily=1,iny ile=ixdep +( ily -1 )* ndec ilf=ile + ixpas izv= ile/nf+1 * if (izv.gt.ite) WRITE(IOIMP,*) ' prob ',izv,ixpas,ixdep, * > ndec,mdec,xdec,ydec,ymi,ymin jd=np1(izv)+1 do 122 j=jd,numnp3 izd=np3(j) if(izd.lt.ile) go to 122 if(izd.gt.ilf) go to 122 ipo= np2(j) C recherche si le point est interne if(idejvu(ipo).eq.1) go to 122 ip=ipt2.num(1,ipo) iref=ip*idim1-idim xpu(1)=xcoor(iref) xpu(2)=xcoor(iref+1) if(xpu(1).lt.xmi.or.xpu(1).gt.xma) go to 122 if(xpu(2).lt.ymi.or.xpu(2).gt.yma) go to 122 * cas quaf qsi(1) = 0.D0 qsi(2) = 0.D0 qsi(3) = 0.D0 if (kell.ne.0) then call qsijs(xell,kell,nbnnl,idim,xpu, $ SHPL,qsi,iret) call qsijs(xel ,kel ,nbnn ,idim,xpu, $ SHPP,qsi,iret) else call qsijs(xel,kel,nbnn,idim,xpu,SHPP,qsi,iret) endif IF (IRET.NE.0) then irate=irate+1 GOTO 122 endif IF (KELL.NE.0) THEN CALL SHAPE(qsi(1),qsi(2),qsi(3),kell,shpl $ ,iret2) if (iret2.eq.0) then MOTERR(1:4)=NOMS(KELL) CALL ERREUR(68) RETURN endif ENDIF * WRITE(IOIMP,fmt='(8E12.5)')(shpp(1,ih),ih=1,nbnn) do i=1,inomil(/1) ilp= inomil(i) if (kell.ne.0) then if( shpl(1,ilp).lt.0.D0) goto 110 else if( shpp(1,ilp).lt.0.D0) goto 110 endif enddo C goto 100 110 xmap=1.D10 do i=1,inomil(/1) ilp=inomil(i) if (kell.ne.0) then xmap = min ( xmap,shpl(1,ilp)) else xmap = min ( xmap,shpp(1,ilp)) endif enddo if( xmap . gt . ckrit(ipo) ) then do 180 jl =1,NUM(/1) if(IP.NE.NUM(jl,iel)) go to 180 if(IFL.EQ.0) go to 190 idejvu(ipo)=1 itr=itr+1 go to 181 180 continue 181 continue ckrit(ipo)=xmap ieint(1,ipo)=isous ieint(2,ipo)=iel QQQ(1,ipo) = qsi(1) QQQ(2,ipo) = qsi(2) endif 190 continue go to 122 100 CONTINUE C on a trouve l element contenant le point ipo attention C si le point appartient aux deux maillages on ne le stocke pas C pour accro mais on le stocke pour proi do 18 jl =1,NUM(/1) if((IP.EQ.NUM(jl,iel)).AND.(IFL.EQ.0)) GOTO $ 19 18 continue if(idejvu(ipo).EQ.0) then idejvu(ipo) = 1 itr = itr + 1 IEINT(1,ipo) = isous IEINT(2,ipo) = iel QQQ(1,ipo) = qsi(1) QQQ(2,ipo) = qsi(2) * WRITE(IOIMP,*) 'trouve point' ,IP, 'Element' ,J * WRITE(IOIMP,2375)ip,j,(xpu(i1),i1=1,idim),(qsi(i2),i2=1,idim) endif 19 continue * 2375 format(2i4,6e12.5) 122 continue 120 continue endif if(idim.eq.3) then IXDEP=INT(real((Xmi-XMIN)/XDEC))+1+ & INT(real((ymi-YMIN)/YDEC))*NDEC + & INT(real((zmi-zmin)/ZDEC))*NMDEC ixpas=INT(real((Xma-XMIN)/XDEC))+1+ & INT(real((ymi-YMIN)/YDEC))*NDEC + & INT(real((zmi-zmin)/ZDEC))*NMDEC - ixdep iny=INT(real((yma-YMIn)/YDEC))-INT(real((ymi-YMIn) $ /YDEC))+1 inz= INT(real((zma-zmin)/ZDEC))-INT(real((zmi-zmin) $ /ZDEC))+1 * if(iel.lt.2.and.isous.eq.1) WRITE(IOIMP,*) ' inz iny',inz,iny do 240 itz=1,inz idexx= ixdep + (itz-1)*nmdec do 220 ily=1,iny ile=idexx +( ily -1 )* ndec ilf=ile + ixpas * if(iel.lt.2.and.isous.eq.1)WRITE(IOIMP,*)' ile,ilf' , ile,ilf izv= ile/nf + 1 jd=np1(izv)+1 * if(iel.lt.2.and.isous.eq.1) WRITE(IOIMP,*) ' jd,jf ',jd,jf do 222 j=jd,numnp3 izd=np3(j) * WRITE(IOIMP,*) ' izd ', izd if(izd.lt.ile) go to 222 if(izd.gt.ilf) go to 223 ipo= np2(j) C recherche si le point est interne if(idejvu(ipo).eq.1) go to 222 ip=ipt2.num(1,ipo) * if(iel.lt.2.and.isous.eq.1)WRITE(IOIMP,*) ' ipo ip ',ipo,ip iref=ip*idim1-idim xpu(1)=xcoor(iref) xpu(2)=xcoor(iref+1) xpu(3)=xcoor(iref+2) * if(iel.lt.2.and.isous.eq.1) WRITE(IOIMP,*) 'cooor ' * $ ,xpu(1),xpu(2),xpu(3) if(xpu(1).lt.xmi.or.xpu(1).gt.xma) go to 222 if(xpu(2).lt.ymi.or.xpu(2).gt.yma) go to 222 if(xpu(3).lt.zmi.or.xpu(3).gt.zma) go to 222 qsi(1) = 0.D0 qsi(2) = 0.D0 qsi(3) = 0.D0 call qsijs(xel,kel,nbnn,idim,xpu,SHPP,qsi,iret) IF (IRET.NE.0) THEN irate=irate+1 if(irate.eq.itest) then * WRITE(IOIMP,*) 'solution pas trouvee', irate itest=itest *2 endif GOTO 222 endif IF (KELL.NE.0) THEN CALL SHAPE(qsi(1),qsi(2),qsi(3),kell,shpl $ ,iret2) ENDIF C * WRITE(IOIMP,fmt='(8E12.5)')(shpp(1,ih),ih=1,nbnn) do i=1,inomil(/1) ilp= inomil(i) if (kell.ne.0) then if( shpl(1,ilp).lt.0.D0) goto 210 else if( shpp(1,ilp).lt.0.D0) goto 210 endif enddo C goto 300 210 xmap=1.D10 do i=1,inomil(/1) ilp=inomil(i) if (kell.ne.0) then xmap = min ( xmap,shpl(1,ilp)) else xmap = min ( xmap,shpp(1,ilp)) endif enddo if( xmap . gt . ckrit(ipo) ) then do 280 jl =1,NUM(/1) if(IP.NE.NUM(jl,iel)) go to 280 if(IFL.EQ.0) go to 290 idejvu(ipo)=1 itr=itr+1 * if(iel.lt.3.and.isous.eq.1) WRITE(IOIMP,*)'egalite de point' go to 281 280 continue 281 continue ckrit(ipo)=xmap ieint(1,ipo)=isous ieint(2,ipo)=iel QQQ(1,ipo) = qsi(1) QQQ(2,ipo) = qsi(2) QQQ(3,ipo) = qsi(3) endif 290 continue go to 222 300 CONTINUE C on a trouve l element contenant le point ipo attention C si le point appartient aux deux maillages on ne le stocke pas C pour accro mais on le stocke pour proi do 28 jl =1,NUM(/1) if((IP.EQ.NUM(jl,iel)).AND.(IFL.EQ.0)) $ GOTO 29 28 continue if(idejvu(ipo).EQ.0) then idejvu(ipo) = 1 * if(iel.lt.3.and.isous.eq.1) WRITE(IOIMP,*)'point interne' itr = itr + 1 IEINT(1,ipo) = isous IEINT(2,ipo) = iel QQQ(1,ipo) = qsi(1) QQQ(2,ipo) = qsi(2) QQQ(3,ipo) = qsi(3) * if(iel.lt.3.and.isous.eq.1) then * WRITE(IOIMP,*) 'trouve point' ,IP, 'Element' ,iel * WRITE(IOIMP,2375)ip,iel,(xpu(i1),i1=1,idim),(qsi(i2),i2=1,idim * endif endif 29 continue 222 continue 223 continue 220 continue 240 continue endif 101 continue segsup inomil,wrk4 IF (KELL.NE.0) THEN SEGSUP wrkl ENDIF segdes meleme 200 continue segdes ipt1 GOTO 800 C========================== C= CAS 1D - DIMENSION 1 = C========================== 600 CONTINUE C Recherche des min et max XMIN=XGRAND XMAX=-XGRAND DO i=1,nodes j=IPT2.NUM(1,i)*IDIM1-IDIM xx=XCOOR(j) XMIN=MIN(xx,XMIN) XMAX=MAX(xx,XMAX) ENDDO XDIST=XMAX-XMIN xdistv=XDIST IF (XDIST.EQ.0.D0) XDIST=1.D0 C Recherche sur les elements de la plus petite taille en x SEGACT,ipt1 nbsous=ipt1.lisous(/1) dminx=XGRAND DO isous=1,nbsous meleme=ipt1.lisous(isous) SEGACT,meleme nbnn=num(/1) DO iel=1,num(/2) dmix=XGRAND dmax=-XGRAND DO i=1,nbnn j=num(i,iel)*IDIM1-IDIM xx=XCOOR(j) dmix=MIN(dmix,xx) dmax=MAX(dmax,xx) ENDDO xx=dmax-dmix dminx=MIN(dminx,xx) ENDDO ENDDO C- WRITE(IOIMP,*) 'dminx',dminx,' xdist',xdist IF (dminx.EQ.0.D0) THEN IF (xdistv.EQ.0.D0) THEN dminx=1.D0 ELSE dminx=xdistv/1000.D0 ENDIF ENDIF dminx=MAX(dminx,xdist/10000.D0) ndex=xdist/dminx ndex=MAX(1,ndex) C- WRITE(IOIMP,*) ' ndex=',ndex C Calcul du decoupage en zones nbzone=NUMNP*NF rap=REAL(nbzone)/ndex rap=MIN(Un,rap) ndec=ndex*rap ndec = max(1,ndec) nbzone=ndec XDEC=xdist/ndec*1.0000001D0 nf=(nbzone/numnp3)+1 nmult=nf C- WRITE(IOIMP,601) ndec,nbzone C- 601 FORMAT(' DECOUPAGE EN X ',I6,' SOIT ',I8,' ZONES ') C- WRITE(IOIMP,602) XMIN,XMAX,XDEC C- 602 FORMAT(' XMIN,XMAX,XDEC ',3G12.5) DO i=1,numnp j=IPT2.num(1,i)*IDIM1-IDIM izone=INT((XCOOR(j)-XMIN)/XDEC)+1 ipoint=izone/nmult+1 np1(ipoint)=np1(ipoint)+1 C- WRITE(IOIMP,*) 'point,izone,ipoint =',i,j,izone,ipoint ENDDO C np1 = nb de points par super-zone de classement DO i=2,NUMNP3 np1(i)=np1(i)+np1(i-1) ENDDO C- WRITE(IOIMP,*) ' np1 avant descente' C- WRITE(IOIMP,*) (np1(i),i=1,ite) DO i=1,numnp j=IPT2.num(1,i)*IDIM1-IDIM izone=INT((XCOOR(j)-XMIN)/XDEC)+1 ipoint=izone/nmult+1 j=np1(ipoint) np3(j)=izone np2(j)=i np1(ipoint)=j-1 ENDDO C- WRITE(IOIMP,*) ' np1 ite avant permutation' ,ite C- WRITE(IOIMP,*) (np1(i),i=1,ite) C- WRITE(IOIMP,*) (np2(i),i=1,ite) C- WRITE(IOIMP,*) (np3(i),i=1,ite) DO i=1,NUMNP jD=np1(i)+2 jF=np1(i+1) IF (jD.LE.jF) THEN 610 IPERM=0 DO j=jD,jF k=j-1 IF (np3(k).GT.np3(j)) THEN L1=np2(k) np2(k)=np2(j) np2(j)=L1 L1=np3(k) np3(k)=np3(j) np3(j)=L1 IPERM=1 ENDIF ENDDO IF (IPERM.EQ.1) GOTO 610 ENDIF ENDDO C- WRITE(IOIMP,*) ' dans zoneg2 apres permutation' C- WRITE(IOIMP,*) ' np1 ite apres perm ',ite C- WRITE(IOIMP,*) (np1(i),i=1,ite) C- WRITE(IOIMP,*) (np2(i),i=1,ite) C- WRITE(IOIMP,*) (np3(i),i=1,ite) C A ce stade : C - np1(i)=j le debut de la super-zone i est en j-eme position + 1 C - np2(j)=k le premier point dans la super zone i est le k-eme C - np3(j)=m sa vrai zone est la m-eme C Il ne reste plus qu'a faire le travail pour chaque element. C On va boucler sur les points qui sont dans les zones en dessous. xpu(2)=0.D0 xpu(3)=0.D0 C SEGACT,ipt1 DO isous=1,nbsous C- WRITE(IOIMP,*) ' boucle isous=',isous meleme=ipt1.lisous(isous) C SEGACT,meleme nbnn=num(/1) SEGINI,wrk4 SEGINI,inomil kel=ITYPEL kd=KDEGRE(kel) C Elements quadratiques IF (kd.EQ.3) then jd=NSPOS(kel) jf=jd+NBSOM(kel)-1 DO i=1,nbnn DO j=jd,jf IF (i.EQ.ibsom(j)) GOTO 620 ENDDO inomil(**)=i 620 CONTINUE ENDDO C Elements lineaires ELSE IF (kd.EQ.2) THEN DO i=1,nbnn inomil(**)=i ENDDO ENDIF DO iel=1,num(/2) CALL DOXE(XCOOR,IDIM,nbnn,num,iel,xel) j=num(1,iel)*IDIM1-IDIM xmi=XCOOR(j) xma=xmi DO i=2,nbnn j=num(i,iel)*IDIM1-IDIM xx=XCOOR(j) xmi=MIN(xmi,xx) xma=MAX(xma,xx) ENDDO C- IF (iel.EQ.1.AND.isous.EQ.1) THEN C- WRITE(IOIMP,*) ' element ',iel C- WRITE(IOIMP,*) (num(i,iel),i=1,nbnn) C- WRITE(IOIMP,*) ' boite reelle : xmi xma ',xmi,xma C- ENDIF xdec1=us3*(xma-xmi) xmi=xmi-xdec1 xma=xma+xdec1 xmi=MAX(xmi,xmin) xmi=MIN(xmax,xmi) xma=MAX(xma,xmin) xma=MIN(xma,xmax) C- IF (iel.EQ.1.AND.isous.EQ.1) THEN C- WRITE(IOIMP,*) ' element ',iel C- WRITE(IOIMP,*) ' boite apres modif : xmi xma ',xmi,xma C- ENDIF C- WRITE(IOIMP,*) ' xmi xma :',xmi,xma C Recherche des points internes a l'element par balayage des super zones ixdep=INT((xmi-xmin)/xdec)+1 ixpas=INT((xma-xmin)/xdec)+1 ile=ixdep ilf=ile+ixpas izv=ile/nf+1 C- IF (izv.GT.ite) WRITE(IOIMP,*) ' prob ',izv,ixpas,ixdep,ndec,xdec jd=np1(izv)+1 DO 630 j=jd,numnp3 izd=np3(j) IF (izd.LT.ile.OR.izd.GT.ilf) GOTO 630 C Recherche si le point est interne ipo=np2(j) IF (idejvu(ipo).EQ.1) GOTO 630 ip=ipt2.num(1,ipo) xpu(1)=XCOOR(ip*IDIM1-IDIM) IF (xpu(1).LT.xmi.OR.xpu(1).GT.xma) GOTO 630 qsi(1) = 0.D0 qsi(2) = 0.D0 qsi(3) = 0.D0 CALL QSIJS(xel,kel,nbnn,IDIM,xpu,shpp,qsi,IRET) C- WRITE(IOIMP,FMT='(8E12.5)') (shpp(1,i),i=1,nbnn) IF (IRET.NE.0) then irate=irate+1 GOTO 630 ENDIF do i=1,inomil(/1) ilp= inomil(i) if( shpp(1,ilp).lt.0.D0) goto 1110 enddo goto 1100 1110 xmap=1.D10 do i=1,inomil(/1) ilp=inomil(i) xmap = min ( xmap,shpp(1,ilp)) enddo if( xmap . gt . ckrit(ipo) ) then do 1180 jl =1,NUM(/1) if(IP.NE.NUM(jl,iel)) go to 1180 if(IFL.EQ.0) go to 1190 idejvu(ipo)=1 itr=itr+1 go to 1181 1180 continue 1181 continue ckrit(ipo)=xmap ieint(1,ipo)=isous ieint(2,ipo)=iel QQQ(1,ipo) = qsi(1) endif 1190 continue go to 630 1100 CONTINUE C on a trouve l element contenant le point ipo attention C si le point appartient aux deux maillages on ne le stocke pas C pour accro mais on le stocke pour proi DO i=1,nbnn IF ((ip.EQ.NUM(i,iel)).AND.(IFL.EQ.0)) GOTO 630 ENDDO IF (idejvu(ipo).EQ.0) THEN idejvu(ipo)=1 itr=itr+1 IEINT(1,ipo)=isous IEINT(2,ipo)=iel QQQ(1,ipo)=qsi(1) C WRITE(IOIMP,*) 'Trouve point de numero global ',IP C WRITE(IOIMP,*) 'dans le ',iel,'eme élément ', C $ 'de la ',isous,'eme sous-zone' C WRITE(IOIMP,*) 'Coordonnée globale = ',xpu(1) C WRITE(IOIMP,*) 'Coordonnée locale = ',qsi(1) CC- WRITE(IOIMP,*) 'Trouve point' ,IP, 'Element' ,J CC- WRITE(IOIMP,603) ip,j,xpu(1),qsi(1) C 603 FORMAT(2I4,2E12.5) ENDIF 630 CONTINUE ENDDO SEGSUP,inomil,wrk4 SEGDES,meleme ENDDO SEGDES IPT1 800 CONTINUE IACCRO=IACCRO+itr C- WRITE(IOIMP,*) ' nb points non converges dans qsijs ',irate SEGDES,IPT2 SEGSUP STRAV RETURN END