www.pudn.com > zaishenghunningtu.rar > three.for, change:2011-05-12,size:6653b


      program main 
      use msflib 
      integer,dimension(800000,3)::ijm 
 
      real,DIMENSION(600000,2):: XY 
      real temp,height,width,step,temp2,temp3 
      integer numx,numy,i,j,K,N 
      integer(4) status 
      open(6,file='骨料参数.txt',status='old') 
      OPEN(5,FILE='lll.TXT',STATUS='UNKNOWN') 
      write(*,*)'	请输入截面高度和宽度和分割步长,单位:毫米' 
      read(*,*)height,broad,step 
	print*,'请输入绘图比例2----10' 
      read*,bili 
	 
       
!计算试件的水平向点数, 
      temp=broad/step 
      if (temp-int(temp).gt.0) then 
      numx=int(temp)+2 
      else 
      numx=int(temp)+1 
      end if 
	 
!计算试件竖直向点数1 
      temp=height/step 
      if(temp-int(temp).gt.0) then 
      numy=int(temp)+2 
      else 
      numy=int(temp)+1 
      end if 
 
!计算各节点坐标,并写入lll文件 
      DO J=1,numy 
      DO I=1,numx 
      XY(I+(J-1)*numx,1)=(I-1)*step 
      XY(I+(J-1)*numx,2)=(J-1)*step 
      temp2=(I-1)*step 
	if(temp2.ge.broad) temp2=broad 
      temp3=(J-1)*step 
	if(temp3.ge.height) temp3=height 
      WRITE(5,*)temp2,temp3 
      END DO 
      END DO 
!计算各单元结点号 
      DO J=2,numy 
      Do I=2,numx 
      IJM(2*I-3+(J-2)*2*(NUMX-1),1)=I+NUMX*(J-2) 
      IJM(2*I-3+(J-2)*2*(NUMX-1),2)=I-1+NUMX*(J-1) 
      IJM(2*I-3+(J-2)*2*(NUMX-1),3)=I-1+NUMX*(J-2) 
      IJM(2*I-2+(J-2)*2*(NUMX-1),1)=I-1+NUMX*(J-1) 
      IJM(2*I-2+(J-2)*2*(NUMX-1),2)=I+NUMX*(J-2) 
      IJM(2*I-2+(J-2)*2*(NUMX-1),3)=I+NUMX*(J-1) 
      end do 
      end do 
!单元结点号调整 
      do k=1, (numx-1)*(numy-1)*2 
      if (mod(k,4).eq.3) then 
      ijm(k,2)=ijm(k,2)+1 
      ijm(K+1,2)=ijm(k+1,2)-1 
      end if 
      end do 
 
      DO K=1,(numx-1)*(numy-1)*2 
      WRITE(5,*) (IJM(K,N),N=1,3) 
      END DO 
      CLOSE(5) 
 
      NJ=numy*numx 
	WRITE(*,*)NJ 
      NE=(numy-1)*(numx-1)*2 
      read(6,*)NG 
      CALL projection(NG,NJ,NE) 
      close(6) 
      call graphicsmode() 
      status=setbkcolorrgb(#ffffff) 
      call clearscreen ($gclearscreen)
      call drawline(bili) 
      end program 
 
       
	 
	 
	subroutine projection(NG,NJ,NE) 
      dimension ijl(NE,3),xy(NJ,2),xyp(NG,2),r(NG),h(NG) 
      open(7,file='LLL.txt',status='old')
      open(8,file='网格数据.txt',status='unknown') 
!向网格数据中写入颗粒数,结点数,单元数 
      write(8,*) NJ,NE 
!读入结点坐标 
      do I=1,NJ 
      read(7,*) (xy(I,j),j=1,2) 
      write(8,*)(xy(I,j),j=1,2) 
      end do 
!读入单元编码数组 
      do I=1,NE 
      read(7,*) (ijl(I,j),j=1,3) 
      end do 
!读入骨料参数 
      do I=1,NG 
      read(6,*)(xyp(I,j),j=1,2),r(I),h(I) 
      end do 
	 
	 
!开始确定每一个单元的属性 
      do I=1,NE 
      ID=0 
      k=0 
      nx1=ijl(I,1) 
      nx2=ijl(I,2) 
      nx3=ijl(I,3) 
       
	do J=1,NG 
	d1=xy(nx1,1)-xyp(J,1) 
	d2=xy(nx1,2)-xyp(J,2) 
      if ((d1**2+d2**2).le.(r(J)-h(J))**2)   k=K+1 
	d3=xy(nx2,1)-xyp(J,1) 
	d4=xy(nx2,2)-xyp(J,2) 
      if ((d3**2+d4**2).le.(r(J)-h(J))**2)   k=K+1 
      d5=xy(nx3,1)-xyp(J,1) 
	d6=xy(nx3,2)-xyp(J,2) 
      if ((d5**2+d6**2).le.(r(J)-h(J))**2)   k=k+1 
      end do 
!该单元为骨料 
      if (k.eq.3) ID=1 
!该单元为老粘结带 
      if (k.eq.2.or.k.eq.1) ID=2 
	 
 
	k=0 
	 
	do J=1,NG 
	d1=xy(nx1,1)-xyp(J,1) 
	d2=xy(nx1,2)-xyp(J,2) 
      if ((d1**2+d2**2).le.(r(J)*r(J)))   k=K+1 
	d3=xy(nx2,1)-xyp(J,1) 
	d4=xy(nx2,2)-xyp(J,2) 
      if ((d3**2+d4**2).le.(r(J)*r(J)))   k=K+1 
      d5=xy(nx3,1)-xyp(J,1) 
	d6=xy(nx3,2)-xyp(J,2) 
      if ((d5**2+d6**2).le.(r(J)*r(J)))   k=k+1 
	end do 
!该单元老砂浆 
	if ((k.eq.3).and.(ID.eq.0)) ID=3 
      if ((k.eq.2.or.k.eq.1).and.(ID.eq.2)) ID=3 
!该单元为新粘结带 
      if ((k.eq.2.or.k.eq.1).and.(ID.eq.0)) ID=4 
      write(8,*) ijl(I,1),ijl(I,2),ijl(I,3),ID 
	end do 
      close(7) 
      close(8) 
      end subroutine 
 
	subroutine graphicsmode() 
      use msflib 
      logical  modestatus 
      integer(2)   maxx,maxy 
      type (windowconfig) myscreen 
      common  maxx,maxy 
      myscreen.numxpixels=-1 
      myscreen.numypixels=-1 
      myscreen.numtextcols=-1 
      myscreen.numtextrows=-1 
      myscreen.numcolors=-1 
      myscreen.fontsize=-1 
      myscreen.title="" 
      modestatus=setwindowconfig(myscreen) 
      modestatus=getwindowconfig(myscreen) 
      maxx=myscreen.numxpixels-1 
      maxy=myscreen.numypixels-1 
	maxx=maxy 
      end subroutine 
!newx-this function finds new x-coordinates 
      INTEGER(2) FUNCTION newx(XCOORD) 
      INTEGER(2) XCOORD,maxx,maxy 
      real(4) tempx 
      common maxx,maxy 
      tempx=maxx/2000.0 
      tempx=xcoord*tempx+0.5 
      newx=tempx 
      end  
!newy-this function finds new y-coordinates 
      INTEGER(2) FUNCTION newy(yCOORD) 
      INTEGER(2) yCOORD,maxx,maxy 
      real(4) tempy 
      common maxx,maxy 
      tempy=maxy/2000.0 
      tempy=ycoord*tempy+0.5 
      newy=tempy 
      end  
 
      subroutine drawline(bili) 
      use msflib 
      external newx,newy 
      real(8),dimension(80000,2)::xy 
      real(8),dimension(80000,3)::xyr 
      integer(2) status,newx,newy,maxx,maxy,unit(80000,4) 
      type (xycoord) poly(3) 
      common maxx,maxy 
      open(6,file='骨料参数.txt',status='old') 
      read(6,*)np 
      do i=1,np  
      read(6,*) (xyr(i,j),j=1,3) 
      c1=newx(int((xyr(i,1)-xyr(i,3))*bili)) 
	c2=newy(int((xyr(i,2)-xyr(i,3))*bili)) 
	c3=newx(int((xyr(i,1)+xyr(i,3))*bili)) 
	c4=newy(int((xyr(i,2)+xyr(i,3))*bili))  
	status=ellipse($gborder,c1,c2,c3,c4) 
      end do 
      close(6) 
       
	open(5,file='网格数据.txt',status='old') 
      read(5,*)NJ,NE 
      do i=1,NJ 
      read(5,*) (xy(i,j),j=1,2) 
      end do 
 
      do i=1,NE 
      read(5,*)(unit(i,j),j=1,4)  
      end do 
      close(5) 
      do i=1,NE 
      poly(1).xcoord=newx(int(xy(unit(i,1),1)*bili)) 
      poly(1).ycoord=newy(int(xy(unit(i,1),2)*bili)) 
      poly(2).xcoord=newx(int(xy(unit(i,2),1)*bili)) 
      poly(2).ycoord=newy(int(xy(unit(i,2),2)*bili)) 
      poly(3).xcoord=newx(int(xy(unit(i,3),1)*bili)) 
      poly(3).ycoord=newy(int(xy(unit(i,3),2)*bili)) 
!给新砂浆附绿色 
      if (unit(i,4).eq.0) status=setcolorrgb(#00FF00) 
!给骨料附深红色 
	if (unit(i,4).eq.1) status=setcolorrgb(#0000FF) 
!给老黏结带附蓝色 
      if (unit(i,4).eq.2) status=setcolorrgb(#FF0000) 
!给老砂浆附蓝色 
      if (unit(i,4).eq.3) status=setcolorrgb(#FFFF00) 
!给新黏结带附蓝色 
      if (unit(i,4).eq.4) status=setcolorrgb(#BA55D3) 
!      if (unit(i,4).eq.1) then 
!      status=setcolorrgb(#800000) 
!      status = polygon($gfillinterior, poly, int(3)) 
!      end if 
      status = polygon($gborder, poly, int(3))  ! 绘制多边形 
      end do 
      end subroutine