www.pudn.com > zaishenghunningtu.rar > tfour.for, change:2011-03-21,size:39282b


 
!1.NJ,NE 
!     NJ--total number of nodes 
!     NE--total number of elements 
!2.MC,NB,ND,NX,EO,VO,T     
!	MC--type of problem     
!           MC=0   plane stress 
!           MC=1   plane strain     
!     NB--total number of given displacements    
!	ND--maxium half band-width of stiffness matrix K 
!     NX--total nnumber of load groups 
!     EO--young's modulus       
!	VO--possin's ratio        
!	T--thickness of structure  
!3.NWA,NWE,NWK,NWP,NWD  
!     CONTROL PARAMETERS OF OUTPUT 
!            =1  output 
!            =0  no output 
!     NWA--output control parameter of element parameters 
!     NWE--output control parameter of element stiffness matrix 
!     NWK--output control parameter of structure stiffness matrix  
!     NWP--output control parameter of load vector 
!     NWD--output control parameter of nodal displacements 
!4.IJM(3,NE)--nodal codes of element 
!                IJM(1,I);IJM(2,I);IJM(3,I) 
!                thre  nodal codes of I-th triangle element 
!5.XY(2, NJ)--coordinates of nodes           
!6.MB(2, NB),ZB(NB)--information and values of prescribed boundary          
!                       MB(1,I)--nodal number of prescribed boundary 
!                       MB(2,I)=1--X-displicement is prescribed 
!                              =0--Y-displacement is prescribed 
!                       ZB(I)--value of prescribed boundary 
!  ID(L)----the identify of element L 
!  IDO 
!7.NF,NP 
!     NF--total number of concentrated forces 
!     NP--total number of uniform forces 
!     if NF>0 
!8.MF(2,NF),ZF(NF) 
!                       MF(1,I)--nodal number applied by concentrated forces 
!                       MF(2,I)=1--X-component is applied 
!                              =0--Y-component is applied 
!                       ZF(I)--values of forces 
!     if NP>0 
!9.MP(2,NP),ZP(NP)   
!                       MP(1,I)--the first nodal number of strip 
!                       MP(2,I)--the second nodal number of strip 
!                       ZP(I)--valueof uniform force                   
!if NX>1 repeat 7-9.(NX-1)times 
!data end of NJ=0 
!************************************************************************************* 
!     arrangement of adjustable arrays 
!       array     C(600000000)        array       IA(5000000) 
!       C(1)      XY(2,NJ)            IA(1)       IJM(3,NE) 
!       C(N1)     ZB(NB)              IA(M1)      MB(2,ZB) 
!       C(N2)     BCA(7,NE)           IA(M2)      MF(2,NB) 
!       C(N3)     SK(NT,ND)           IA(M3)      MP(2,NP) 
!       C(N4)     F(NT)               IA(MEND)    limit needed 
!       C(N5)     ZF(NF)              
!       C(N6)     ZP(NP) 
!       C(NEND)   limit needed 
!************************************************************************************* 
! error stop 
!       stop 111    exceed the limit of array C 
!       STOP 222    exceed the limit of array C/IA 
!       stop 333    element area is nonpositive 
!       stop 444    main element of matrix is nonpositve 
!*************************************************************************************   
!main program  
	use msflib 
      DIMENSION IA(5000000),C(60000000),EK(36) 
      DIMENSION emax(70000),IDO(70000),ID(70000) 
      open(1,file='结点位移.TXT',status='unknown') 
	open(2,file='刚度矩阵.TXT',status='unknown') 
	open(3,file='荷载列阵.TXT',status='unknown') 
	open(4,file='应力数据.TXT',status='unknown') 
	open(5,file='输入参数.txt',status='old') 
      open(6,file='计算结果.TXT',status='unknown') 
      open(7,file='网格数据.txt',status='old') 
      open(8,file='应变数据.TXT',status='unknown') 
	open(9,file='破坏数据.txt',status='unknown') 
      read(7,*) NJ,NE 
      if (NJ.EQ.0) STOP 
      write(6,15) 
15    FORMAT(//,5X,'**********plane problem**********',/) 
      read(5,*)MC,NX,NF,NP,EO,VO,T 
      NF=int(sqrt(NJ+0.05))         !作用集中荷载的结点个数 
	NB=NF+1                       !给定位移个数 
      ND=2*NF+4                     !结构刚度阵的半带宽 
	NX=1                      
      print*,'请输入荷载值及步长' 
      read(*,*)force,step 
      print*,'请输入显示因子' 
      read(*,*)yz 
	print*,'请输入试件的高和宽' 
	read(*,*) height,broad 
      write(6,20)NJ,NE,MC,NX,NB,ND,EO,VO,T 
20    FORMAT(//,1X,'NJ=',I5,3X,'NE=',I5,3X,'MC=',I2,3X,'NX=',I2,3X,'NB=' 
     &,I4,3X,'ND=',I5,3X,'EO=',E10.4,3X,'VO=',F7.2,3X,'T=',F7.2) 
       
	read(5,*)ft0,ft1,ft2,ft3,ft4    !砂浆、黏结带、骨料的抗拉强度 
      read(5,*)e0,e1,e2,e3,e4        !砂浆、黏结带、骨料的初始弹模 
	read(5,*) NWA,NWE,NWK,NWP,NWD 
      write(6,25)NWA,NWE,NWK,NWP,NWD 
25    format(/10X,'NWA=',I2,3x,'NWE=',I2,3x,'NWK=',I2,3x,'NWP=',I2,3x 
     &,'NWD=',I2) 
      
       
	NT=2*NJ            !结点位移个数 
!!!!!!!!!!!!!!!!!!!!!!! calculate the limit of adjustable arrays!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
      M1=3*NE+1 
      M2=M1+2*NB 
      N1=2*NJ+1 
      N2=N1+NB 
      N3=7*NE+N2 
      N4=N3+NT*ND 
      N5=N4+NT 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!check the limit of array C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!    
      NEND=N5 
      IF (NEND.LE.60000000) GOTO 35 
      WRITE(*,*)' exceed the limit of array C' 
      WRITE(*,30) NEND 
30    FORMAT(/,'***********NEND**********',I6,1X,'>10000000') 
      STOP 111 
 
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!input data!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
35    CALL input(NE,NJ,NB,IA(1),C(1),IA(M1),C(N1),ID,IDO,NF) 
      WRITE(*,40) 
40    FORMAT(/10X,'########input passed######') 
       
!!!!!!!!!!!!!!!!!!!!!!!!!! calculating element parameters!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
      CALL ABC(NE,NJ,NWA,IA(1),C(1),C(N2))				     
      WRITE(*,55) 
55    FORMAT(/10X,'########ABC passed########') 
       
 
 
	ntempt0=0         !超过弹性应变单元数的控制变量 
      ntempt1=0         !超过残余应变单元数的控制变量 
	ntempt2=0         !超过极限应变单元数的控制变量 
       
600   continue 
       
	DO I=N3,N4-1       !把数组SK(NT,ND)中的数附零 
      C(I)=0.0 
      END DO 
       
 
 
      DO 65 K=1,NE 
      IO=K 
      if (ID(io).eq.14) then 
      vo=0.2 
      eo=0.1*ft1/emax(io) 
      if (eo.lt.e1/100) eo=e1/100 
      if (eo.gt.e1/50) eo=e1/50 
      end if 
       
      if (ID(io).eq.13) then 
      vo=0.2 
      eo=0.1*ft4/emax(io) 
      if (eo.lt.e4/100) eo=e4/100 
	if (eo.gt.e4/30) eo=e4/30 
      end if 
 
      if (ID(io).eq.12) then 
      vo=0.2 
      eo=0.1*ft2/emax(io) 
      if (eo.lt.e2/100) eo=e2/100 
	if (eo.gt.e2/30) eo=e2/30 
      end if 
 
      if (ID(io).eq.11) then 
      vo=0.2 
      eo=0.1*ft3/emax(io) 
      if (eo.lt.e3/100) eo=e3/100 
	if (eo.gt.e3/40) eo=e4/40 
      end if 
 
      if (ID(io).eq.10) then 
      vo=0.2 
      eo=0.1*ft0/emax(io) 
      if (eo.lt.e0/100) eo=e0/100 
      if (eo.gt.e0/40) eo=e0/40 
      end if 
 
      if (ID(io).eq.9) then 
      vo=0.2 
      eo=e1*(4.9/4.0*ft1/e1/emax(io)-0.9/4.0 ) 
      if (eo.lt.e1/50) eo=e1/50 
      if (eo.gt.e1) eo=e1 
      end if 
 
	if (ID(io).eq.8) then 
      vo=0.15 
      eo=e4*(2.9/4.0*ft4/e4/emax(io)-0.9/2.0 ) 
      if (eo.lt.e4/30) eo=e4/30 
      if (eo.gt.e4) eo=e4 
      end if 
 
	if (ID(io).eq.7) then 
      vo=0.2 
      eo=e2*(2.9/4.0*ft2/e2/emax(io)-0.9/2.0 ) 
      if (eo.lt.e2/30) eo=e2/30 
      if (eo.gt.e2) eo=e2 
      end if 
 
	if (ID(io).eq.6) then 
      vo=0.2 
      eo=e3*(3.9/2.0*ft3/e3/emax(io)-0.9/3.0 ) 
      if (eo.lt.e3/40) eo=e3/40 
      if (eo.gt.e3) eo=e3 
      end if 
  
      if (ID(io).eq.5) then 
      vo=0.2 
      eo=e0*(3.9/2.0*ft0/e0/emax(io)-0.9/3.0 ) 
      if (eo.lt.e0/40) eo=e0/40 
      if (eo.gt.e0) eo=e0 
      end if 
 
      if (ID(io).eq.1) then    !给骨料附泊松比、弹模 
      vo=0.2 
      eo=e1 
      end if 
 
      if (ID(io).eq.2) then   !给老黏结带附泊松比、弹模 
      vo=0.2 
      eo=e2 
      end if 
 
	if (ID(io).eq.3) then   !给老砂浆附泊松比、弹模 
      vo=0.2 
      eo=e3 
      end if 
 
      if (ID(io).eq.4) then   !给新黏结带附泊松比、弹模 
      vo=0.2 
      eo=e4 
      end if 
 
	if (ID(io).eq.0) then   !给新砂浆附泊松比、弹模 
      vo=0.2 
      eo=e0 
      end if 
       
	if (ID(io).eq.15) then 
      vo=0.2 
      eo=1 
      end if 
       
	 
	IF (MC.EQ.0) GOTO 45 
! plane strain problem 
      E=EO/(1.0-VO*VO) 
      V=VO/(1.0-VO) 
      GOTO 50  
! plane stress problem 
45    E=EO 
      V=VO 
 
50    NX1=NX               !荷载组数为1 
      A1=E/(1.0-V*V)/4.0   !材料系数 
      A2=0.5*(1.0-V)       !材料系数 
       
       
       
!!!!!!!!!!!!!!!!!!!!!!!!!! form matrix K!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
      CALL KE(IO,NE,NWE,T,A1,A2,V,EK,C(N2)) 
      CALL SUMK(IO,NE,ND,NT,IA(1),C(N3),EK) 
65    continue 
	 
	WRITE(*,70) 
70    FORMAT(/10x,'#######sumK passed########') 
      write(2,*)(C(I),I=N3,N4-1) 
 
	CALL CHECK(NT,ND,NWK,C(N3)) 
      WRITE(*,75) 
75    FORMAT(/10x,'#######check passed#######') 
 
 
 
 
80    continue 
!!!!!!!!!!!!!!!!!!! calculate the limit of adjustable-size arrays!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
!!!!!!!!!!!!!!!!!!!!!!check the limit of array C and array IA!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
      M3=2*NF+M2 
      N6=N5+NF 
      NEND=N6+NP-1 
      MEND=M3+2*NP-1 
      NM=0 
      IF(NEND.LE.60000000) GOTO 85 
      WRITE(*,*)' exceed the limit of array C' 
      WRITE(*,30) NEND 
      NM=1 
85    IF(MEND.LE.5000000) GOTO 95 
      WRITE(*,*)' exceed the limit of array IA' 
      WRITE(*,90) MEND 
90    FORMAT(/,'##########MEND#######',I6,1X,'>5000000') 
      STOP 222 
95    IF(NM.EQ.1) STOP 222 
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!form vector P!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
      DO 100 I=N4,N5-1 
	C(I)=0.0 
100   CONTINUE  
	IF(NF.GT.0) CALL PF(NF,NP,NT,NWP,C(N4),IA(M2),C(N5),force) 
      WRITE(*,105) 
105   FORMAT(/10x,'########PF passed#########') 
      IF (NP.GT.0) CALL PP(NP,NT,NJ,NWP,C(1),C(N4),IA(M3),C(N6)) 
      WRITE(*,110) 
110   FORMAT(/10x,'########PP passed#########') 
   
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!introduce given displacements!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
      CALL DBC(NT,ND,NB,C(N3),C(N4),IA(M1),C(N1)) 
      WRITE(*,115) 
115   FORMAT(/10x,'########DBC passed########') 
  
 
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!solve linear equations KA=P!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
      call GAUSS(NT,ND,NWD,NX,NX1,C(N3),C(N4)) 
      WRITE(*,120) 
120   FORMAT(/10x,'#######GAUSS passed#######') 
   
 
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!calculate element stresses!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
      CALL STRESS(NJ,NE,NT,A1,A2,V,IA(1),C(N2),C(N4),ID,MC,nrh,ncy,yl, 
     & yb,IDO,nph,emax,force,ft0,ft1,ft2,ft3,ft4,T,NF,e0,e1,e2,e3,e4, 
     & height,broad) 
 	   
	WRITE(*,125) 
125   FORMAT(/10x,'######STRESS passed#######') 
       
	 
	 
	NX1=NX1-1 
      IF (NX1.GT.0) GOTO 80  
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!drawing pictures!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!	 
	if (ntempt0.eq.0) call graphicsmode() 
      status=settextcolorrgb(#000000) 
      if ((nrh.ne.ntempt0).or.(ncy.ne.ntempt1).or.(nph.ne.ntempt2)) then 
      ntempt0=nrh 
      ntempt1=ncy  
      ntempt2=nph 
      end if 
	WRITE(*,130) 
130   FORMAT(/10x,'#####graphicsmode passed#####') 
	 
	 
	call graphic(NJ,NE,ID,C(1),IA(1),IDO,step,yz) 
      write(*,150)nrh,ncy,nph 
150   format(1x,'软化',I10,'残余',I10,'破坏',I10)    
       
 
       
	         
	 
!      if (force.ge.0) write(8,*)yb*1000000 
!                      write(4,*)yl 
	if (force.le.0) write(8,*)yb*(-1000000) 
	                write(4,*)yl*(-1) 
      print*,yb,yl,force 
	if (step.eq.0) stop 
!      if (force.ge.0) force=force+step 
      if (force.lE.0) force=force-step 
      print*,force 
	 
      print*,'程序计算中,请勿动,谢谢!' 
       
	GOTO 600 
 
      end program  
 
     
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
!   input informations of nodes and elements boundary conditions   !  
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
      subroutine input(NE,NJ,NB,IJM,XY,MB,ZB,ID,IDO,NF) 
      DIMENSION IJM(3,NE),XY(2,NJ),MB(2,NB),ZB(NB),ID(NE),IDO(NE) 
       
	do J=1,NJ 
      read(7,*)(XY(I,J),I=1,2) 
      end do 
    
      do L=1,NE 
      read(7,*)(IJM(I,L),I=1,3),ID(L) 
      end do 
      IDO=ID 
 
       
 
      do n=1,NF 
      MB(1,n)=n        !给定下边界的结点号 
      MB(2,n)=0        !给定了Y方向的位移 
      ZB(n)=0          !试件底面各给定结点位移为零 
      end do 
       
	MB(1,NB)=NB/2     !下边界的中间点号 
      MB(2,NB)=1             !给定了X向的位移 
      ZB(NB)=0               !中间点位移为零 
 
      write(6,10) 
10    format(/10X,'单元结点IJM(3,NE)',/) 
      write(6,20)((IJM(M,I),M=1,3),I=1,NE) 
20    format(1X,3(I5,4X)) 
      write(6,30) 
30    format(/10X,'结点坐标XY(2,NJ)',/) 
      write(6,40)((XY(M,I),M=1,2),I=1,NJ) 
40    format(1X,E12.5,3X,E12.5) 
      write(6,50) 
50    format(/10X,'给定的位移信息和数值',20X,'MB',15X,'ZB',/) 
      DO I=1,NB 
      write(6,60)(MB(K,I),K=1,2),ZB(I) 
60    format(10X,2I5,5X,E14.6) 
      end do 
       
	return 
      end subroutine   
 
       
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
!                     calculate element parameters	             ! 
!                          input: IJM, XY                          ! 
!                           output: BAC                            ! 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 	 
!X(2,5)--current calculating element's coordinates of nodes   
!        [xi,xj,xm,xi,xj 
!         yi,yj,ym,yi,yj] 
!B(7)--current calculating element's parameters 
!       B(1)-B(3)  bi,bj,bm 
!       B(4)-B(6)  ci,cj,cm 
!       B(7)     element's area 
	SUBROUTINE ABC(NE,NJ,NWA,IJM,XY,BCA) 
      DIMENSION IJM(3,NE), XY(2,NJ),BCA(7,NE),X(2,5),B(7) 
      IF (NWA.EQ.1) WRITE(6,5) 
5     FORMAT(/10X,'parameters of elements BCA(7,NE)'/) 
       
	DO 80 I=1,NE 
      DO 10 K=1,3 
      K1=IJM(K,I) 
      DO 10 J=1,2 
      X(J,K)=XY(J,K1)    !单元I上结点K的坐标 
10    continue 
 
      DO 20 J=1,2 
      X(J,4)=X(J,1) 
      X(J,5)=X(J,2) 
20    continue        
   
      DO 30 K=1,3 
	B(K)=X(2,K+1)-X(2,K+2)           
      B(K+3)=X(1,K+2)-X(1,K+1) 
30    continue 	 
       
	B(7)=(B(1)*B(5)-B(4)*B(2))*0.5      
 
	IF(NWA.EQ.1) WRITE(6,40)I,B 
40    FORMAT(1X,'NE=',I3,/3X,7E10.4) 
	 
      IF(B(7).LE.0.0) GOTO 60 
 
	DO J=1,7 
	BCA(J,I)=B(J) 
	END DO 
      GOTO 80 
 
60    WRITE(6,70)I,(IJM(J,I),J=1,3) 
70    FORMAT(/5X,'element', I5, 5X, 'area is nonpositive',5X,'IJM',3I5) 
	STOP 333 
	 
80    continue 
	return 
	END SUBROUTINE 
 
	 
	 
	 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
!         calculate stiffness matrix of element                 ! 
!               input: BCA                                      ! 
!               output: EK                                      ! 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 	 
!BCA(7,NE)--array of element parameter 
!EK(6,6)--stiffness matrix of element  
	SUBROUTINE KE(IO,NE,NWE,T,A1,A2,V,EK,BCA) 
      dimension B(7),BCA(7,NE),EK(6,6) 
	EK=0.0 
	DO I=1,7 
	B(I)=BCA(I,IO) 
	END DO 
	A=A1/B(7)*T     !单刚前面的系数 
	DO I=1,3 
	 DO J=1,3 
	 I1=2*I 
	 J1=2*J 
	 EK(I1-1,J1-1)=A*(B(I)*B(J)+A2*B(I+3)*B(J+3)) 
	 EK(I1-1,J1)=A*(V*B(I)*B(J+3)+A2*B(I+3)*B(J)) 
	 EK(I1,J1-1)=A*(V*B(I+3)*B(J)+A2*B(I)*B(J+3)) 
	 EK(I1,J1)=A*(B(I+3)*B(J+3)+A2*B(I)*B(J)) 
	 END DO 
	END DO 
 
	DO I=3,6 
	 DO J=1,I 
	 EK(I,J)=EK(J,I) 
	 end do 
	end do 
 
	IF (NWE.EQ.0) GOTO 60 
	WRITE(6,40) IO 
40    FORMAT(/1X,'EK     NE=',I5) 
	WRITE(6,50) EK 
50    FORMAT(1X,6E11.4) 
60    return 
	END SUBROUTINE 
 
	 
	 
	 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
! form the structure stiffness matrix of the maximum half band-width ! 
!         input:  EK , IJM                                           ! 
!         output:  SK                                                ! 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!	 
!IO--the number of element, coming from main program 
!EK(6,6)--stiffness matrix of element 
!IJM(3,NE)--nodal codes of element 
!SK(NT,ND)--the structure stiffness matrix of the maximum half band-width 
	SUBROUTINE SUMK(IO,NE,ND,NT,IJM,SK,EK) 
	DIMENSION IJ(3),IJM(3,NE),SK(NT,ND),EK(6,6) 
	DO I=1,3 
	IJ(I)=IJM(I,IO) 
	end do 
	 
	DO 20 I=1,3 
	DO 20 J=1,3 
	    IF (IJ(I).GT.IJ(J)) GOTO 20    !下三角元素,舍去 
		M=2*IJ(I)-1 
		N=2*(IJ(J)-IJ(I))+1 
		MO=2*I-1 
		NO=2*J-1 
		SK(M,N)=SK(M,N)+EK(MO,NO) 
		SK(M,N+1)=SK(M,N+1)+EK(MO,NO+1) 
		SK(M+1,N)=SK(M+1,N)+EK(MO+1,NO+1) 
		IF (IJ(I).EQ.IJ(J)) GOTO 20 
		SK(M+1,N-1)=SK(M+1,N-1)+EK(MO+1,NO) 
20    continue 
	return 
	END SUBROUTINE 
 
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
!check whether the main elements of stiffness matrixis K are positive   ! 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
      SUBROUTINE CHECK(NT,ND,NWK,SK) 
	dimension SK(NT,ND) 
      IF (NWK.EQ.0) GOTO 30 
	WRITE(6,10) 
10    FORMAT(/,10X,'globle stiffness matrix',/) 
	WRITE(6,20)((SK(I,J),I=1,NT),J=1,ND) 
20    FORMAT(1X,5E13.6) 
	 
30    M=0 
	DO I=1,NT 
	IF (SK(I,1).GT.1E-10) GOTO 50				  
	WRITE(6,40)I,SK(I,1) 
40    FORMAT(/10X,'main element is nonpositive NT=',I6,5X,E12.6) 
	M=M+1 
50    END DO 
	IF (M.GT.0) GOTO 60 
	GOTO 70 
60    STOP 444 
70    RETURN  
	END SUBROUTINE 
 
	 
	 
	 
	 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
!      form the vector of nodal loads produced by concentrated forces      ! 
!		input:  NF,MF,ZF,F			                                     ! 
!		 output:   F				                                     ! 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
!F--the equivalent vector of nodal loads 
	SUBROUTINE PF(NF,NP,NT,NWP,F,MF,ZF,force) 
      DIMENSION MF(2,NF),ZF(NF),F(NT) 
       
	do I=1,NF           	 
	MF(1,I)=NT/2-I+1   !加荷载的结点码 
	MF(2,I)=0          !y方向加荷载 
	ZF(I)=force 
	end do					 
      ZF(1)=force/2 
      ZF(NF)=ZF(1) 
 
	 
	DO I=1,NF 
	N=2*MF(1,I)-MF(2,I) 
	F(N)=ZF(I) 
	END DO 
	   
	 
	WRITE(3,50) 
50    FORMAT(/,10X,'load vector'/ ) 
	WRITE(3,60)F 
60    FORMAT(1X,E14.6) 
 
	END SUBROUTINE  
 
	 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
!  form the vector of nodal force produced by uniform force              !                                       
!       input:    F,XY,NP,MP,ZP		                                   ! 
!       output:      F                                                   ! 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
!F--the equivalent vector of nodal loads 
	SUBROUTINE PP(NP,NT,NJ,NWP,XY,F,MP,ZP) 
	DIMENSION MP(2,NP),	ZP(NP),XY(2,NJ),F(NT) 
	READ(5,*)((MP(I,L),I=1,2),L=1,NP),(ZP(L),L=1,NP) 
	WRITE(6,10) 
10    FORMAT(/, 10X,'uniform pressure MP(2,NP) ZP(NP)',/) 
	DO L=1,NP 
      WRITE(6,20)(MP(I,L),I=1,2),ZP(L) 
20    FORMAT(1X,2I10,10X,E12.6) 
	END DO 
 
	DO I=1,NP 
	N1=MP(1,I) 
	N2=MP(2,I) 
	PX=XY(2,N1)-XY(2,N2) 
	PY=XY(1,N2)-XY(1,N1) 
	PX=.5*ZP(I)*PX 
	PY=.5*ZP(I)*PY 
	F(2*N1-1)=F(2*N1-1)+PX 
	F(2*N1)=F(2*N1)+PY 
  	F(2*N2-1)=F(2*N2-1)+PX 
  	F(2*N2)=F(2*N2)+PY 
	END DO 
 
	IF (NWP.EQ.0)	GOTO 70 
	WRITE(6,50) 
50    FORMAT(/,10X,'uniform loads'/) 
	WRITE(6,60)F 
60    FORMAT(1X,5E14.6) 
70    RETURN  
	END SUBROUTINE  
 
 
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
!     introduce given nodal displacement                                    ! 
!       input:A,B,NB,MB,ZB                                                  ! 
!       output: A,B                                                         ! 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!      
!A(NT,ND)--coefficient matrix 
!B(NT)--the equivalent vector of nodal loads 
!NB--total number of given displacements 
      subroutine DBC(NT,ND,NB,A,B,MB,ZB) 
      dimension MB(2,NB),ZB(NB),A(NT,ND),B(NT) 
      NX=1 
      NX1=1 
	 
      do 60 I=1,NB 
      N=2*MB(1,I)-MB(2,I) 
      Z=ZB(I) 
      
	IF (ABS(Z).LT.1E-10) GOTO 20 
      IF (NX.NE.NX1) GOTO 10 
      write(*,*)nx1 
	A(N,1)=A(N,1)*1E+15 
10    B(N)=A(N,1)*Z 
      GOTO 60 
20    IF (NX.NE.NX1) GOTO 50 
      A(N,1)=1.0 
      DO 30 J=2,ND 
      A(N,J)=0.0 
30    CONTINUE 
      DO 40 K=2,ND 
      IF (N.LT.K) GOTO 50 
      M=N-K+1 
      A(M,K)=0.0 
40    CONTINUE 
 
50    B(N)=0.0 
60    CONTINUE 
       
	RETURN 
      END SUBROUTINE 
 
 
 
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
!                    gauss elimination                                ! 
!                      input:A,B                                      ! 
!                      output:B                                       ! 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
!B(NT)--the equivalent vector of nodal loads before calculation 
!     --nodal displacement after calculation 
 
      subroutine GAUSS(NT,ND,NWD,NX,NX1,A,B) 
      dimension A(NT,ND),B(NT) 
      NX=1														  
      NX1=1 
      N=NT-1 
	 
      if (NX.eq.NX1) goto 10 
      assign 50 to M 
      goto 20 
10    assign 30 to M 
20    do 60 k=1,N 
      M1=ND-1 
      if ((M1).gt.(NT-K)) M1=NT-K 
 
      do 60 L=1,M1 
      C=A(k,L+1)/A(k,1) 
      if (ABS(c).lt.1e-18) goto 60 
      goto M,(30,50) 
30    M2=ND-L 
       
      DO 40 J=1,M2 
      A(K+L,J)=A(K+L,J)-C*A(K,J+L) 
40    CONTINUE 
 
50    B(K+L)=B(K+L)-C*B(K) 
60    CONTINUE 
 
      B(NT)=B(NT)/A(NT,1) 
 
      DO 80 K=1,N 
      I=NT-K 
      M1=ND 
      C=B(I) 
      IF ((K+1).LT.(ND)) M1=K+1 
 
      DO 70 J=2,M1 
      L=I+J-1 
      C=C-A(I,J)*B(L) 
70    CONTINUE 
      B(I)=C/A(I,1) 
80    CONTINUE 
 
!output nodal displacements 
       
      WRITE(1,90) 
      WRITE(1,100) 
90    FORMAT(/,10X,'displacements of nodes',/) 
100   FORMAT(1X,2(1X,'NG',9X,'U',12X,'V',9X)) 
      goto 150 
      N=NT/2 
      N11=N/2 
      IF (FLOAT(N11)+.3-FLOAT(N)/2.0.GT.1E-7) N=N-1 
      DO 140 I=1,N,2 
      NT2=NT/4 
      NT3=NT/2 
      FL1=FLOAT(NT2)+.3-FLOAT(NT3)/2 
      IF ((FL1.LT.0).AND.(I.EQ.N)) GOTO 120 
      J=2*I-1 
      K=2*I 
      I1=I+1 
      J1=J+2 
      K1=K+2 
      WRITE(6,110)I,B(J),B(K),I1,B(J1),B(K1) 
110   FORMAT(1X,2(I5,3X,E11.5,2x,e11.5,5X)) 
      GOTO 140 
120   J=2*I-1 
      K=2*I 
      WRITE(6,130)I,B(J),B(K) 
130   FORMAT(1X,I5,3X,E11.5,2X,E11.5) 
140   CONTINUE 
150   RETURN 
      END SUBROUTINE	    
 
	 
	 
	 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
!		       calculate the element stresses and strain             ! 
!                   input: F,IJM,BCA,A1,A2,V  	                     ! 
!                    output: S-X,S-Y,S-XY,S1,S2,AG 	                 ! 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!	 
 
      SUBROUTINE STRESS(NJ,NE,NT,A1,A2,V,IJM,BCA,F,ID,MC,nrh,ncy,yl, 
     & yb,IDO,nph,emax,force,ft0,ft1,ft2,ft3,ft4,T,NF,e0,e1,e2,e3,e4, 
     & height,broad) 
       
	DIMENSION IJM(3,NE),ID(NE),BCA(7,NE) 
	DIMENSION F(NT),B(7),R(6),IDO(NE),emax(NE) 
	 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
!            calculate the element strain                                   ! 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
!F(NT)--nodal displacement 
!nrh--the number of element whose strain excess elastic strain 
!ncy--the number of element whose strain excess residul strain 
!nph--the number of element whose strain excess ultimate strain  
!ybx--liner strain of X 
!yby--liner strain of Y 
!ybxy--sheering strain 
!ybmax--the maximum principle strain	 
      nrh=0	 
      ncy=0 
      nph=0 
       
	 
      do 10 I=1,NE 
      if ((ID(I).ge.5).and.(ID(I).le.9))	nrh=nrh+1 
      if ((ID(I).ge.10).and.(ID(I).le.14))	ncy=ncy+1 
      if (ID(I).eq.15) nph=nph+1 
10    continue     
       
      DO 20 I=1,NE 
      IO=I 
       
 
      DO J=1,7 
	B(J)=BCA(J,IO) 
	end do 
 
	 
	DO J=1,3 
	N=IJM(J,IO)*2 
      R(2*J-1)=F(N-1) 
      R(2*J)=F(N) 
	END DO 
	     
	   
      ybx=(B(1)*R(1)+B(2)*R(3)+B(3)*R(5))/2/B(7)              !X向应变 
       
	yby=(B(4)*R(2)+B(5)*R(4)+B(6)*R(6))/2/B(7)              !Y向应变 
       
	y0=(B(4)*R(1)+B(1)*R(2)+B(5)*R(3)+B(2)*R(4)+B(6)*R(5)+B(3)*R(6)) 
      ybxy=y0/2/B(7)                                          !剪应变 
       
	ybmax=(ybx+yby+sqrt((ybx-yby)*(ybx-yby)+ybxy*ybxy))/2 
 
      if (ID(IO).eq.14) then 
      vo=0.2 
      eo=0.1*ft1/emax(io) 
      if (eo.lt.e1/100) eo=e1/100 
      if (eo.gt.e1/50) eo=e1/50 
      end if 
 
	if (ID(IO).eq.13) then 
      vo=0.2 
      eo=0.1*ft4/emax(io) 
      if (eo.lt.e4/100) eo=e4/100 
      if (eo.gt.e4/30) eo=e4/30 
      end if 
 
	if (ID(IO).eq.12) then 
      vo=0.2 
      eo=0.1*ft2/emax(io) 
      if (eo.lt.e2/100) eo=e2/100 
      if (eo.gt.e2/30) eo=e2/30 
      end if 
 
      if (ID(io).eq.11) then 
      vo=0.2 
      eo=0.1*ft3/emax(io) 
      if (eo.lt.e3/100) eo=e3/100 
	if (eo.gt.e3/40) eo=e4/40 
      end if 
 
      if (ID(IO).eq.10) then 
      vo=0.2 
      eo=0.1*ft0/emax(io) 
      if (eo.lt.e0/100) eo=e0/100 
	if (eo.gt.e0/40) eo=e0/40 
      end if 
 
      if (ID(io).eq.9) then 
      vo=0.2 
      eo=e1*(4.9/4.0*ft1/e1/emax(io)-0.9/4.0 ) 
      if (eo.lt.e1/50) eo=e1/50 
      if (eo.gt.e1) eo=e1 
      end if 
 
	if (ID(io).eq.8) then 
      vo=0.2 
      eo=e4*(2.9/4.0*ft4/e4/emax(io)-0.9/2.0 ) 
      if (eo.lt.e4/30) eo=e4/30 
      if (eo.gt.e4) eo=e4 
      end if 
  
      if (ID(io).eq.7) then 
      vo=0.2 
      eo=e2*(2.9/2.0*ft2/e2/emax(io)-0.9/2.0 ) 
      if (eo.lt.e2/30) eo=e2/30 
      if (eo.gt.e2) eo=e2 
      end if 
 
	if (ID(io).eq.6) then 
      vo=0.2 
      eo=e3*(3.9/2.0*ft3/e3/emax(io)-0.9/3.0 ) 
      if (eo.lt.e3/40) eo=e3/40 
      if (eo.gt.e3) eo=e3 
      end if 
 
      if (ID(io).eq.5) then 
      vo=0.2 
      eo=e0*(3.9/3.0*ft0/e0/emax(io)-0.9/3.0) 
      if (eo.lt.e0/40) eo=e0/40 
      if (eo.gt.e0) eo=e0  
      end if 
 
      if (ID(io).eq.1) then    !给骨料附泊松比、弹模 
      vo=0.2 
      eo=e1 
      end if 
 
      if (ID(io).eq.2) then    !给老黏结带附泊松比、弹模 
      vo=0.22 
      eo=e2 
      end if 
 
	if (ID(io).eq.3) then   !给老砂浆附泊松比、弹模 
      vo=0.2 
      eo=e3 
      end if 
 
	if (ID(io).eq.4) then   !给新黏结带附泊松比、弹模 
      vo=0.2 
      eo=e4 
      end if 
 
      if (ID(io).eq.0) then   !给新砂浆附泊松比、弹模 
      vo=0.22 
      eo=e0 
      end if 
      if (ID(io).eq.15) then 
      vo=0.2 
      eo=1 
      end if 
       
	if (emax(IO).lt.ybmax) emax(IO)=ybmax 
	IF (MC.EQ.0) GOTO 45 
       
!平面应变问题 
 
      E=EO/(1.0-VO*VO) 
      V=VO/(1.0-VO) 
 
      GOTO 50  
 
 !平面应力问题 
45    E=EO 
      V=VO 
50    continue 
      A1=E/(1.0-V*V)/4.0 
      A2=0.5*(1.0-V) 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
!                calculate the element stresses                           ! 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
!S1--normal stress of X 
!S2--normal stress of Y 
!S3--sheering stress 
!X1--the maximun principle stress 
!X2--the minmum principle stress 
!CTA--angle between principal stress S1 and X axial 
       
      S1=0.0 
	S2=0.0 
	S3=0.0 
	    
	DO J=1,7 
	B(J)=BCA(J,IO) 
	end do 
 
	A=2*A1/B(7)         !应力矩阵前的系数 
	      
	DO J=1,3 
	N=IJM(J,IO)*2 
      R(2*J-1)=F(N-1) 
      R(2*J)=F(N) 
	END DO 
	     
	   DO J=1,3 
	   K=2*J 
	   S1=S1+A*(B(J)*R(K-1)+V*B(J+3)*R(K))          !x向的应力 
   	   S2=S2+A*(V*B(J)*R(K-1)+B(J+3)*R(K))          !y向的应力 
	   S3=S3+A*A2*(B(J+3)*R(K-1)+B(J)*R(K))         !剪应力 
	   END DO 
	   P=0.5*(S1+S2) 
	   Q=0.5*(S1-S2) 
	   X1=P+SQRT(Q*Q+S3*S3)                         !最大主应力 
	   X2=2*P-X1                                    !最小主应力 
	   CTA=0 
	  
	IF (S3.GT.0) CTA=ATAN((X1-S1)/S3)               !主平面与x轴的夹角 
 
 
	 
!判断超过弹性单元 
      if (ID(IO).le.4) then	  
       
	if (ID(IO).eq.0) then           !如果是新砂浆 
      if ((x1.gt.ft0).or.(x2.gt.ft0).or.(x1.lt.-10*ft0).or.(x2.lt. 
     & -10*ft0).or.(abs(s3).gt.ft0) )then 
!	if (ybmax.gt.(ft0/e0)) then     !如果砂浆单元超过弹性阶段 
!      IF(S2.GT.ft0)THEN 
	ID(IO)=5                        !给超过弹性段的砂浆单元附4 
	nrh=nrh+1 
	write(9,*)IO,ID(IO),force 
	end if 
      end if	  
       
	if (ID(IO).eq.3) then           !如果是老砂浆 
      if ((x1.gt.ft3).or.(x2.gt.ft3).or.(x1.lt.-10*ft3).or.(x2.lt. 
     & -10*ft3).or.(abs(s3).gt.ft3) )then 
!	if (ybmax.gt.(ft3/e3)) then     !如果砂浆单元超过弹性阶段 
!      IF(S2.GT.ft3)THEN 
	ID(IO)=6                        !给超过弹性段的砂浆单元附4 
	nrh=nrh+1 
	write(9,*)IO,ID(IO),force 
	end if 
      end if 
  
  	if (ID(IO).eq.2)  then           !如果是老黏结带 
      if ((x1.gt.ft2).or.(x2.gt.ft2).or.(x1.lt.-10*ft2).or.(x2.lt. 
     &-10*ft2).or.(abs(s3).gt.ft2) )then 
!	if (ybmax.gt.(ft2/e2)) then      !如果黏结带单元超过弹性阶段 
!      IF(S2.GT.ft2)THEN 
      ID(IO)=7                         !给超过弹性段的黏结带单元附5 
	nrh=nrh+1 
	write(9,*)IO,ID(IO),force 
	end if 
      end if 
 
	if (ID(IO).eq.4)  then           !如果是新黏结带 
      if ((x1.gt.ft4).or.(x2.gt.ft4).or.(x1.lt.-10*ft4).or.(x2.lt. 
     &-10*ft4).or.(abs(s3).gt.ft4) )then 
!	if (ybmax.gt.(ft4/e4)) then      !如果黏结带单元超过弹性阶段 
!      IF(S2.GT.ft4)THEN 
      ID(IO)=8                         !给超过弹性段的黏结带单元附5 
	nrh=nrh+1 
	write(9,*)IO,ID(IO),force 
	end if 
      end if 
 
  	if (ID(IO).eq.1)  then           !如果是骨料 
      if ((x1.gt.ft1).or.(x2.gt.ft1).or.(x1.lt.-10*ft1).or.(x2.lt. 
     & -10*ft1).or.(abs(s3).gt.ft1) )then 
!	if (ybmax.gt.(ft1/e1)) then      !如果骨料单元超过弹性阶段 
!      IF(S2.GT.ft1)THEN 
      ID(IO)=9                         !给超过弹性段的骨料单元附6 
	nrh=nrh+1 
      write(9,*)IO,ID(IO),force 
	end if 
      end if 
 
      end if	   
 
 
 
!判断超过残余应变单元 
      if ((ID(IO).ge.5).and.(ID(IO).le.9)) then 
       
	if (IDO(IO).eq.0) then             !如果是新砂浆 
             
	if (ybmax.gt.4*ft0/e0)  then       !砂浆单元的应变超过残余应变 
      ID(IO)=10                           !给超过残余应变的砂浆单元附7 
	nrh=nrh-1 
	ncy=ncy+1 
	write(9,*)IO,ID(IO),force 
	end if 
      end if  
 
	if (IDO(IO).eq.3) then             !如果是老砂浆 
             
	if (ybmax.gt.4*ft3/e3)  then       !砂浆单元的应变超过残余应变 
      ID(IO)=11                           !给超过残余应变的砂浆单元附7 
	nrh=nrh-1 
	ncy=ncy+1 
	write(9,*)IO,ID(IO),force 
	end if 
      end if  
  
      if (IDO(IO).eq.2) then              !如果是老黏结带  
      if (ybmax.gt.3*ft2/e2)  then        !黏结带单元的应变超过残余应变 
	ID(IO)=12                            !给超过残余应变的黏结带单元附8 
      nrh=nrh-1  
      ncy=ncy+1 
	write(9,*)IO,ID(IO),force 
      end if 
      end if  
 
	if (IDO(IO).eq.4) then              !如果是新黏结带  
      if (ybmax.gt.3*ft4/e4)  then        !黏结带单元的应变超过残余应变 
	ID(IO)=13                            !给超过残余应变的黏结带单元附8 
      nrh=nrh-1  
      ncy=ncy+1 
	write(9,*)IO,ID(IO),force 
      end if 
      end if  
 
      if (IDO(IO).eq.1) then              !如果是骨料 
      if (ybmax.gt.5*ft1/e1)  then        !骨料单元的应变超过残余应变 
	ID(IO)=14                            !给超过残余应变的黏结带单元附9 
	nrh=nrh-1 
	ncy=ncy+1 
	write(9,*)IO,ID(IO),force  
      end if 
      end if  
 
      end if 
 
!  判断破坏单元  
      if ((ID(IO).ge.10).and.(ID(IO).le.14)) then 
 
      if (IDO(IO).eq.0) then               !如果是新砂浆 
      if (ybmax.gt.10*ft0/e0)  then       !砂浆单元的应变超过极限应变 
	ID(IO)=15                           !给超过极限应变的砂浆单元附15 
	ncy=ncy-1 
	nph=nph+1 
	write(9,*)IO,ID(IO),force 
	end if 
      end if  
 
	if (IDO(IO).eq.3) then               !如果是老砂浆 
      if (ybmax.gt.10*ft3/e3)  then       !砂浆单元的应变超过极限应变 
	ID(IO)=15                           !给超过极限应变的砂浆单元附15 
	ncy=ncy-1 
	nph=nph+1 
	write(9,*)IO,ID(IO),force 
	end if 
      end if  
  
      if (IDO(IO).eq.2) then              !如果是老黏结带  
      if (ybmax.gt.10*ft2/e2) then       !黏结带单元的应变超过极限应变 
      ID(IO)=15                          !给超过极限应变的黏结带单元附15 
      ncy=ncy-1 
	nph=nph+1 
      write(9,*)IO,ID(IO),force 
      end if 
      end if  
 
	if (IDO(IO).eq.4) then              !如果是新黏结带  
      if (ybmax.gt.10*ft4/e4) then       !黏结带单元的应变超过极限应变 
      ID(IO)=15                          !给超过极限应变的黏结带单元附15 
      ncy=ncy-1 
	nph=nph+1 
      write(9,*)IO,ID(IO),force 
      end if 
      end if  
 
      if (IDO(IO).eq.1) then              !如果是骨料      
      if (ybmax.gt.10*ft1/e1) then       !骨料单元的应变超过极限应变 
      ID(IO)=15                          !给超过极限应变的骨料单元附15 
	ncy=ncy-1 
	nph=nph+1 
      write(9,*)IO,ID(IO),force  
      end if 
      end if  
 
      end if	  
 
20    continue 
       
	WRITE(*,200) 
200   FORMAT(/1x,'输入截面的宽度和高度') 
 
      yl=force*(NF-1)/broad           !求得是顶面的平均应力 
         
	total=0  
      do m=1,NF  
      total=total+F(NT-2*m+2)     !叠加顶面各接点的竖向位移 
      end do 
      yb=total/NF/height           !求得是顶面的平均应变 
       
       
	end subroutine 
 
	    
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
!                      setting the window's mode                                ! 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
      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 
      maxy=maxx 
      end subroutine 
 
       
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
!                  draw responding picture                                      ! 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!	 
	INTEGER(2) FUNCTION NEWX(XCOORD) 
      INTEGER(2) XCOORD,MAXX,MAXY 
      real(4) tempx 
      common maxx,maxy 
      tempx=maxx/1500.0 
      tempx=xcoord*tempx+0.5 
      newx=tempx 
      end  
 
      INTEGER(2) FUNCTION NEWy(yCOORD) 
      INTEGER(2) yCOORD,MAXX,MAXY 
      real(4) tempy 
      common maxx,maxy 
      tempy=maxy/1500.0 
      tempy=ycoord*tempy+0.5 
      newy=tempy 
      end  
 
 
 
      subroutine graphic(NJ,NE,ID,XY,IJM,IDO,step,yz) 
      use msflib 
      external newx,newy 
      dimension XY(2,NJ) 
      integer(2) status,newx,newy,maxx,maxy 
      dimension IJM(3,NE)	,ID(NE) ,IDO(NE) 
      type (xycoord) poly(3) 
      common maxx,maxy 
      status=setbkcolorrgb(#ffffff)                        !给界面附纯白色 
      CALL CLEARSCREEN($gclearscreen)                   !调用内部函数图形子程序         
 
      do I=1,NE 
      poly(1).xcoord=newx(int(xy(1,ijm(1,i))*yz))+400 
      poly(1).ycoord=newy(int(xy(2,ijm(1,i))*yz))+80 
      poly(2).xcoord=newx(int(xy(1,ijm(2,i))*yz))+400  
      poly(2).ycoord=newy(int(xy(2,ijm(2,i))*yz))+80 
      poly(3).xcoord=newx(int(xy(1,ijm(3,i))*yz))+400 
      poly(3).ycoord=newy(int(xy(2,ijm(3,i))*yz))+80 
 
      if (IDO(i).eq.0) then   
	status=setcolorrgb(#00FF00)    !给新砂浆附绿色 
      status = polygon($gborder, poly, int(3)) 
      end if 
 
	if (IDO(i).eq.3) then   
	status=setcolorrgb(#FFFF00)    !给老砂浆附绿色 
      status = polygon($gborder, poly, int(3)) 
      end if 
  
 
      if (IDO(i).eq.1) then 
	status=setcolorrgb(#0000ff)    !给骨料附红色 
	status = polygon($gborder, poly, int(3))   
      end if 
	 
	if (IDO(i).eq.2) then           
	status=setcolorrgb(#ff0000)    !给老黏结带附蓝 
      status = polygon($gborder, poly, int(3)) 
      end if 
 
	if (IDO(i).eq.4) then           
	status=setcolorrgb(#BA55D3)    !给新黏结带附蓝 
      status = polygon($gborder, poly, int(3)) 
      end if 
       
       
	if (ID(i).ge.5)  then 
	status=setcolorrgb(#00ff00)      !给超过弹性单元的新砂浆附绿色 
	if (step.eq.0) status=setcolorrgb(#ff0000)  
      status = POLYGON($gfillinterior, poly, int(3)) 
      status = POLYGON($gborder, poly, int(3)) 
      end if 
 
      if (ID(i).ge.6)  then 
	status=setcolorrgb(#ffff00)      !给超过弹性单元的老砂浆附青色 
	if (step.eq.0) status=setcolorrgb(#ff0000)  
      status = POLYGON($gfillinterior, poly, int(3)) 
      status = POLYGON($gborder, poly, int(3)) 
      end if 
   
      if ((ID(i).ge.8).or.(ID(i).ge.8))  then 
	status=setcolorrgb(#BA55D3)      !给超过弹性单元的新黏结带附洋红 
	if (step.eq.0) status=setcolorrgb(#ff0000)  
      status = POLYGON($gfillinterior, poly, int(3)) 
      status = POLYGON($gborder, poly, int(3)) 
      end if 
 
	if ((ID(i).ge.7).or.(ID(i).ge.7))  then 
	status=setcolorrgb(#ff0000)      !给超过弹性单元的老黏结带附蓝色 
	if (step.eq.0) status=setcolorrgb(#ff0000)  
      status = POLYGON($gfillinterior, poly, int(3)) 
      status = POLYGON($gborder, poly, int(3)) 
      end if 
	 
	 
	if (ID(i).ge.9)  then 
	status=setcolorrgb(#0000ff)      !给超过弹性单元的骨料附红色 
	if (step.eq.0) status=setcolorrgb(#ff0000)  
       
      status = POLYGON($gborder, poly, int(3)) 
      end if 
	 
	if (ID(i).ge.10)  then             
	status=setcolorrgb(#00ff00)      !给超过残余应变单元的新砂浆附绿色 
      if (step.eq.0) status=setcolorrgb(#0000ff)  
      status = POLYGON($gfillinterior, poly, int(3)) 
      status = POLYGON($gborder, poly, int(3)) 
      end if 
 
	if (ID(i).ge.11)  then             
	status=setcolorrgb(#ffff00)      !给超过残余应变单元的老砂浆附绿色 
      if (step.eq.0) status=setcolorrgb(#0000ff)  
      status = POLYGON($gfillinterior, poly, int(3)) 
      status = POLYGON($gborder, poly, int(3)) 
      end if 
       
	if ((ID(i).ge.13).or.(ID(i).ge.13))  then             
	status=setcolorrgb(#BA55D3)      !给超过残余应变单元的新黏结带附深灰色 
      if (step.eq.0) status=setcolorrgb(#0000ff)  
      status = POLYGON($gfillinterior, poly, int(3)) 
      status = POLYGON($gborder, poly, int(3)) 
      end if 
 
	if ((ID(i).ge.12).or.(ID(i).ge.12))  then             
	status=setcolorrgb(#ff0000)      !给超过残余应变单元的老黏结带附深灰色 
      if (step.eq.0) status=setcolorrgb(#0000ff)  
      status = POLYGON($gfillinterior, poly, int(3)) 
      status = POLYGON($gborder, poly, int(3)) 
      end if 
 
	if (ID(i).ge.14)  then             
	status=setcolorrgb(#0000ff)      !给超过残余应变单元的骨料附暗洋红色 
      if (step.eq.0) status=setcolorrgb(#0000ff)  
      status = POLYGON($gfillinterior, poly, int(3)) 
      status = POLYGON($gborder, poly, int(3)) 
      end if 
 
	if (ID(i).ge.15)  then          !给超过极限应变单元的砂浆,骨料、黏结带附黑 
	status=setcolorrgb(#000000) 
      status = POLYGON($gfillinterior, poly, int(3)) 
      end if 
 
      end do 
 
      end subroutine