zonag
C ZONAG SOURCE OF166741 24/07/25 21:15:04 11950 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 elementaire massif C ip2 pointeur sur le maillage poi1 C iexx pointeur sur un segment contenant les points deja vus C (eviter comptage des points aux frontieres de 2 sous zones C et le traitement des points qui sont des noeuds de ip1) C ITR : COMPTEUR MIS A JOUR DANS ZONAG C Sorties : C iccoun pointeur sur le segment ICOUNT contenant le tableau de C correspondance suivant C IEINT(1,N) numero du point du meleme IP2 C IEINT(2,N) numero de l element du meleme IP1 le contenant C IEINT(3,N) numero de la souszone ISZ correspondant a IP1 C MAXPZ nombre de points de IP2 contenus dans la s_zone (N max) 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 iexclu INTEGER IDEJVU(nbpti) ENDSEGMENT SEGMENT icount INTEGER NBPTOT INTEGER MAXPZ(nbsz) INTEGER IEINT(3,NODES) REAL*8 QQQ(3,NODES) ENDSEGMENT SEGMENT inomil(0) segment wrk4 real*8 qsi(3),xel(3,nbnn),shpp(6,nbnn) endsegment segment wrkl real*8 shpl(6,nbnnl),xell(3,nbnnl) endsegment SEGMENT ISEG1 ENDSEGMENT SEGMENT ISEG3 INTEGER NIZO(NZO+1), NUMZO(NZO), IDEJ(NZO) ENDSEGMENT SEGMENT ISEG5 INTEGER NNMEL(ILON) ENDSEGMENT DIMENSION XPU(3) IDIMP1 = IDIM + 1 c*d write(ioimp,*) 'ZONAGE',isz,ip1,ip2 IPT2 = IP2 C* SEGACT IPT2 <- Actif en E/S et non modifie NODES = IPT2.NUM(/2) IPT1 = IP1 C* SEGACT IPT1 <- Actif en E/S et non modifie NBNN1 = IPT1.NUM(/1) NBEL1 = IPT1.NUM(/2) KELE1 = IPT1.ITYPEL c*d write(ioimp,*) 'ZONAGE',nbnn1,nbel1,kele1, nodes icount = ICCOUN C* SEGACT,icount*MOD <- Actif en E/S en etat modifiable iexclu = IEXX C* SEGACT,iexclu*MOD <- Actif en E/S en etat modifiable NBEL = NBEL1 SEGINI ISEG1 XZO = 0.D0 YZO = 0.D0 ZZO = 0.D0 XZA = XGRAND YZA = XGRAND ZZA = XGRAND XTOMI = +XGRAND XTOMA = -XGRAND YTOMI = +XGRAND YTOMA = -XGRAND ZTOMI = +XGRAND ZTOMA = -XGRAND D_X = 0.D0 D_Y = 0.D0 D_Z = 0.D0 DO I1 = 1, NBEL1 XMI = +XGRAND XMA = -XGRAND YMI = +XGRAND YMA = -XGRAND ZMI = +XGRAND ZMA = -XGRAND c* iexclu.IDEJVU(IB) = 2 <- Deja fait dans ACCRO.eso IA = (IB-1)*IDIMP1 XMI = MIN(XMI,XCOOR(IA+1)) XMA = MAX(XMA,XCOOR(IA+1)) YMI = MIN(YMI,XCOOR(IA+2)) YMA = MAX(YMA,XCOOR(IA+2)) IF (IDIM.EQ.3) THEN ZMI = MIN(ZMI,XCOOR(IA+3)) ZMA = MAX(ZMA,XCOOR(IA+3)) ENDIF ENDDO XLIM(1,I1) = XMI XLIM(2,I1) = XMA YLIM(1,I1) = YMI YLIM(2,I1) = YMA D_X = XMA - XMI XZO = MAX(XZO,D_X) XZA = MIN(XZA,D_X) XTOMI = MIN(XMI,XTOMI) XTOMA = MAX(XMA,XTOMA) D_Y = YMA - YMI YZO = MAX(YZO,D_Y) YZA = MIN(YZA,D_Y) YTOMI = MIN(YMI,YTOMI) YTOMA = MAX(YMA,YTOMA) IF (IDIM.EQ.3) THEN ZLIM(1,I1)=ZMI ZLIM(2,I1)=ZMA D_Z = ZMA - ZMI ZZO = MAX(ZZO,D_Z) ZZA = MIN(ZZA,D_Z) ZTOMI = MIN(ZMI,ZTOMI) ZTOMA = MAX(ZMA,ZTOMA) ENDIF ENDDO D_X = 1.D-4 * (XTOMA - XTOMI) XPR = MIN(XZO*1.D-2, D_X*0.5D0) XTOMI = XTOMI - D_X XTOMA = XTOMA + D_X D_X = XTOMA - XTOMI XZA = MIN(XZA*0.97D0, D_X) NXZO = D_X / XZA + 1 XZO = XZA D_Y = 1.D-4 * (YTOMA - YTOMI) YPR = MIN(YZO*1.D-2, D_Y*0.5D0) YTOMI = YTOMI - D_Y YTOMA = YTOMA + D_Y D_Y = YTOMA - YTOMI YZA = MIN(YZA*0.97D0, D_Y) NYZO = D_Y / YZA + 1 YZO = YZA c*d WRITE(6,FMT='('' NXZO NYZO'',2I7)') NXZO,NYZO NZZO = 1 IF (IDIM.EQ.3) THEN D_Z = (ZTOMA - ZTOMI) * 1.D-4 ZPR = MIN(ZZO*1.D-2, 0.5D0*D_Z) ZTOMI = ZTOMI - D_Z ZTOMA = ZTOMA + D_Z D_Z = ZTOMA - ZTOMI ZZA = MIN(ZZA*0.97D0, D_Z) IF (ZZA.EQ.0.D0) THEN NZZO = 1 ELSE NZZO = D_Z / ZZA + 1 ENDIF ZZO = ZZA c*d WRITE(6,FMT='('' zz0,zzA,ztomi,ztoma'',4e12.5)') c*d $ xzo,xza,ztomi,ztoma ENDIF c*d WRITE(6,*) XTOMI, XTOMA c*d WRITE(6,*) YTOMI ,YTOMA c*d WRITE(6,*) ZTOMI ,ZTOMA XZO = EPSI D_X = XTOMA - XTOMI ENDIF YZO = EPSI D_Y = YTOMA - YTOMI ENDIF IF (IDIM.EQ.3) THEN ZZO = EPSI D_Z = ZTOMA - ZTOMI ENDIF ENDIF C_______________________________ NXDEP = MIN(NXZO,10) NYDEP = MIN(NYZO,10) IF (IDIM.EQ.2) THEN r_z = FLOAT(NXZO)*FLOAT(NYZO) IF (r_z.GT.100000.) THEN XY=SQRT(r_z)/90 NXZO=MAX(INT(NXZO/XY),NXDEP) NYZO=MAX(INT(NYZO/XY),NYDEP) r_z = FLOAT(NXZO)*FLOAT(NYZO) IF (r_z.GT.100000.) THEN XY=SQRT(r_z)/60 NXZO=MAX(INT(NXZO/XY),NXDEP) NYZO=MAX(INT(NYZO/XY),NYDEP) ENDIF XZO= D_X /NXZO YZO= D_Y /NYZO NXZO = D_X /XZO +1 NYZO = D_Y /YZO +1 ENDIF c*d WRITE(6,FMT='('' XZO NXZO YZO NYZO '' , E12.5,I5,E12.5,I5)') c*d $ XZO ,NXZO, YZO, NYZO ELSE NZDEP=MIN(NZZO,10) r_z = FLOAT(NXZO)*FLOAT(NYZO)*FLOAT(NZZO) IF (r_z .GT. 250000.) THEN XYZ = r_z**0.33333333D0/25 NXZO=MAX(INT(NXZO/XYZ),NXDEP) NYZO=MAX(INT(NYZO/XYZ),NYDEP) NZZO=MAX(INT(NZZO/XYZ),NZDEP) r_z = FLOAT(NXZO)*FLOAT(NYZO)*FLOAT(NZZO) IF (r_z.GT.250000.) THEN XYZ = r_z**0.3333333333D0/25 NXZO=MAX(INT(NXZO/XYZ),NXDEP) NYZO=MAX(INT(NYZO/XYZ),NYDEP) NZZO=MAX(INT(NZZO/XYZ),NZDEP) c*d IF(IIMPI.NE.0)WRITE(IOIMP,FMT='('' NXZO NYZO NZZO '' c*d $ ,4I7) ') NXZO,NYZO,NZZO r_z = FLOAT(NXZO)*FLOAT(NYZO)*FLOAT(NZZO) IF (r_z.GT.250000.) THEN XYZ = r_z**0.3333333333D0/25 NXZO=MAX(INT(NXZO/XYZ),NXDEP) NYZO=MAX(INT(NYZO/XYZ),NYDEP) NZZO=MAX(INT(NZZO/XYZ),NZDEP) ENDIF ENDIF XZO=D_X / NXZO YZO=D_Y / NYZO ZZO=D_Z / NZZO NXZO= D_X / XZO +1 NYZO= D_Y / YZO +1 NZZO= D_Z / ZZO +1 ENDIF ENDIF * * ON VEUT CONSTRUIRE LA LISTE DES ELEMENTS TOUCHANT UNE ZONE * POUR CELA ON COMMENCE PAR COMPTER COMBIEN D'ELEMENTS TOUCHENT * CHAQUE ZONE ET EN MEME TEMPS ON STOCKE LES ZONES TOUCHEES * PAR CHAQUE ELEMENT ET LEUR NOMBRE NXYZO = NXZO*NYZO NZO = NXYZO*NZZO c*d WRITE(IOIMP,FMT='('' NZO NXZO NYZO NZZO '' c*d $,4I7) ') NZO,NXZO,NYZO,NZZO SEGINI ISEG3 DO I1=1,NBEL1 NIZ1X = INT( (XLIM(1,I1)-XTOMI-XPR) / XZO ) + 1 NIZ1Y = INT( (YLIM(1,I1)-YTOMI-YPR) / YZO ) + 1 NIZ2X = INT( (XLIM(2,I1)-XTOMI+XPR) / XZO ) + 1 NIZ2Y = INT( (YLIM(2,I1)-YTOMI+YPR) / YZO ) + 1 IF (IDIM.EQ.3) THEN NIZ1Z=INT((ZLIM(1,I1)-ZTOMI-ZPR)/ZZO) +1 NIZ2Z=INT((ZLIM(2,I1)-ZTOMI+ZPR)/ZZO) +1 DO L3=NIZ1Z,NIZ2Z DO L1=NIZ1Y,NIZ2Y DO L2=NIZ1X,NIZ2X NIZA = L2 + (L1-1) * NXZO + (L3-1)*NXYZO NUMZO(NIZA) = NUMZO(NIZA) +1 ENDDO ENDDO ENDDO ELSE DO L1=NIZ1Y,NIZ2Y DO L2=NIZ1X,NIZ2X NIZA = L2 + ( L1-1) * NXZO NUMZO(NIZA) = NUMZO(NIZA) +1 ENDDO ENDDO ENDIF ENDDO * * CONSTRUCTION DU TABLEAU D'ADRESSAGE DU TABLEAU DONNANT LES * ELEMENTS CONCERNES PAR UNE ZONE * ILON=0 NIZO(1)=1 DO L1=1,NZO NIZO(L1+1)=NIZO(L1)+NUMZO(L1) ILON=ILON+ NUMZO(L1) ENDDO C WRITE(6,FMT='('' ILON '',I5)') ILON C WRITE(6,109) (KKK,NUMZO(KKK),(NELZO(KI,KKK),KI=1,4),KKK=1,NBEL) C 109 FORMAT(I6,I5,4I5) C WRITE(6,110)( NIZO(KI),KI=1,NZO+1) C 110 FORMAT(16I5) SEGINI ISEG5 NIZ1X=INT((XLIM(1,I1)-XTOMI-XPR)/XZO) +1 NIZ1Y=INT((YLIM(1,I1)-YTOMI-YPR)/YZO) +1 NIZ2X=INT((XLIM(2,I1)-XTOMI+XPR)/XZO) +1 NIZ2Y=INT((YLIM(2,I1)-YTOMI+YPR)/YZO) +1 IF(IDIM.EQ.3) THEN NIZ1Z=INT((ZLIM(1,I1)-ZTOMI-ZPR)/ZZO) +1 NIZ2Z=INT((ZLIM(2,I1)-ZTOMI+ZPR)/ZZO) +1 DO 251 L3=NIZ1Z,NIZ2Z DO 2511 L1=NIZ1Y,NIZ2Y DO 25111 L2=NIZ1X,NIZ2X NIZA = L2 + (L1-1) * NXZO + (L3-1)*NXYZO IAD=NIZO(NIZA)+IDEJ(NIZA) NNMEL(IAD)=I1 IDEJ(NIZA)=IDEJ(NIZA)+1 25111 CONTINUE 2511 CONTINUE 251 CONTINUE ELSE DO 252 L1=NIZ1Y,NIZ2Y DO 2521 L2=NIZ1X,NIZ2X NIZA = L2 + ( L1-1) * NXZO IAD=NIZO(NIZA)+IDEJ(NIZA) NNMEL(IAD)=I1 IDEJ(NIZA)=IDEJ(NIZA)+1 2521 CONTINUE 252 CONTINUE ENDIF ENDDO nbnn = NBNN1 segini wrk4 segini inomil C element quaf ? si oui kell=numero lineaire correspondant C si non kell=0 IF (IERR.NE.0) RETURN IF (KELL.NE.0) THEN nsom = nbsom(KELE1) nbnnl = NBNNE(kell) IF (nsom.NE.nbnnl) THEN return endif idx = nspos(KELE1)-1 segini wrkl do i=1,nbnnl inomil(**)=i enddo ELSE kd = kdegre(KELE1) C element quadratiques if (kd.eq.3) then nso = nbsom(KELE1) idx = nspos(KELE1)-1 do 762 i=1,nbnn do 763 j=1,nso iso = ibsom(idx+j) if (i.eq.iso) goto 762 763 continue inomil(**)= i 762 continue C elements lineaires else if (kd.eq.2) then do i=1,nbnn inomil(**)=i enddo endif ENDIF c*d write(6,*) (inomil(i),i=1,inomil(/1)) C* IA1 = 0 C* IF (KELE1.EQ.14 .OR. KELE1.EQ.15) IA1 = 1 C* IF (KELE1.EQ.16 .OR. KELE1.EQ.17) IA1 = 7 C on va tourner sur les elements poi1 DO 2 K = 1, NODES IP = IPT2.NUM(1,K) C write(ioimp,*) 'on cherche l element contenant le noeud' ,ip IF (iexclu.IDEJVU(IP).NE.0) GOTO 2 IREFP = (IP-1) * IDIMP1 XPU(1) = XCOOR(IREFP+1) IF (XPU(1).LT.XTOMI.OR.XPU(1).GT.XTOMA) GOTO 2 XPU(2) = XCOOR(IREFP+2) IF (XPU(2).LT.YTOMI.OR.XPU(2).GT.YTOMA) GOTO 2 XPU(3) = 0.D0 IF (IDIM.EQ.3) THEN XPU(3) = XCOOR(IREFP+3) IF (XPU(3).LT.ZTOMI.OR.XPU(3).GT.ZTOMA) GOTO 2 ENDIF INDZO=INT((XPU(1)-XTOMI)/XZO)+ 1 +INT((XPU(2)-YTOMI)/YZO)*NXZO IF (IDIM.EQ.3) INDZO=INDZO+INT((XPU(3)-ZTOMI)/ZZO)*NXZO*NYZO IDEB = NIZO(INDZO) IFIN = NIZO(INDZO+1)-1 C write(6,fmt='('' ideb ifin'',2i5)') ideb,ifin ,indzo IF (IDEB.GT.IFIN) GO TO 2 C ITROU = 1 C------------------------------------------------- DO 11 KK = IDEB,IFIN jkk = NNMEL(KK) * write(6,*) ' ' * write(6,*) 'recherche point ',ip, 'ds elem ',jkk,'itrou',ITROU C Dans le cas "ACCRO", on cherche tout d'abord si le noeud C appartient aux deux maillages car si oui on ne le stocke pas dans C les noeuds a traiter ensuite (pas de rigidite a construire) C Deja fait auparavant ! c* DO jl = 1, NBNN1 c* IF (IP.EQ.IPT1.NUM(jl,jkk)) THEN c* ITROU = 2 c* GOTO 19 c* ENDIF c* ENDDO * cas quaf IF (kell.NE.0) THEN DO innl = 1, nbnnl jl = ibsom(idx+innl) DO iid = 1, IDIM xell(iid,innl) = xel(iid,jl) ENDDO ENDDO ELSE ENDIF IF (IRET.NE.0) GOTO 11 c*d write(6,fmt='(8E12.5)')(shpp(1,ih),ih=1,nbnn) if (kell.ne.0) then do i=1,inomil(/1) ilp= inomil(i) if (shpl(1,ilp).lt.0.D0) goto 11 enddo else do i=1,inomil(/1) ilp= inomil(i) if (shpp(1,ilp).lt.0.D0) goto 11 enddo endif GOTO 100 11 CONTINUE C on n a rien trouve pour ce point ITROU = 0 100 CONTINUE 19 CONTINUE IF (ITROU.EQ.1) THEN C on a trouve l element contenant le point K attention iexclu.IDEJVU(ip) = 1 itr = itr + 1 icount.IEINT(1,itr) = IP icount.IEINT(2,itr) = JKK icount.IEINT(3,itr) = ISZ icount.QQQ(1,itr) = qsi(1) icount.QQQ(2,itr) = qsi(2) icount.QQQ(3,itr) = qsi(3) * write(6,*) 'trouve point' ,IP, 'Element' ,JKK C write(6,2375)ip,j,(xpu(i1),i1=1,IDIM),(qsi(i2),i2=1,IDIM) 2375 format(2i4,6e12.5) C* ELSE IF (ITROU.EQ.2) THEN C* iexclu.IDEJVU(ip) = 2 C*C write(6,*) 'trouve point' ,IP, 'noeud Element' ,jll,JKK ENDIF 2 CONTINUE icount.MAXPZ(ISZ) = ITR - ICOUNT.NBPTOT icount.NBPTOT = ITR segsup wrk4,inomil IF (KELL.NE.0) SEGSUP wrkl SEGSUP,ISEG1,ISEG5,ISEG3 C C le tableau des correspondances est termine C C* SEGDES,IPT1,IPT2 <- Actifs en E/S et non modifies C* SEGDES,icount,iexclu <- Actifs en E/S et modifies RETURN END
© Cast3M 2003 - Tous droits réservés.
Mentions légales