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(*,*)'	请输入截面高度和宽度和分割步长，单位：毫米'
print*,'请输入绘图比例2----10'

!计算试件的水平向点数,
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
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
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
write(8,*)(xy(I,j),j=1,2)
end do
!读入单元编码数组
do I=1,NE
end do
!读入骨料参数
do I=1,NG
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')
do i=1,np
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')
do i=1,NJ
end do

do i=1,NE
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

```