eqsolve2.for

静力学有限元源程序

	END DO
      END IF

      DIAG(1) = 1
	DO I = 2,TN+1
	  DIAG(I) = DIAG(I-1)+(I-1)-PP(I-1)+1
	END DO
      END SUBROUTINE BAND_K

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      SUBROUTINE FORM_K()
	!FORM THE TOTAL STIFF MATRIX KK
      USE DEF      
      IMPLICIT NONE
      !ABOVE TO TRANSFER
      INTEGER(4)::I,J,K,T_NODE(4,2),TEMP
	INTEGER(4),ALLOCATABLE::GG(:)
      EXTERNAL::FORM_KE,FORM_ORDER
       
      COORD=0.0
      DO I = 1,N_ELE
        IF (G_ELE(I,3) .EQ. G_ELE(I,4)) THEN
	    S34 = 3
	  ELSE
	    S34 = 4
	  END IF
	  ALLOCATE(GG(1:S34*2))
        DO J = 1,S34
	    COORD(J,1) = G_NODE(G_ELE(I,J),1)
	    COORD(J,2) = G_NODE(G_ELE(I,J),2)   !RINGT
	  END DO
 

	  CALL FORM_KE()
	  !PRINT*,'ELEMENT:',I
        !DO J = 1,S34
	  ! 	PRINT*,'LOCATION:',G_ELE(I,J)
        !END DO	 
        !DO J = 1,S34*2
        ! PRINT*,KE(J,1:S34*2)
	  !END DO
	  !PRINT*

	  DO J = 1,S34
	      GG(2*J-1)=2*G_ELE(I,J)-1
	      GG(2*J)=2*G_ELE(I,J)
        END DO
	  DO J =1,S34*2
	    DO K = 1,S34*2
	      IF (GG(K) .LE. GG(J)) THEN        !这里一定要判断一下,因为单元编号的顺序不一定从小到大
	        TEMP = DIAG(GG(J))+GG(J)-GG(K)
	        KK(TEMP) = KK(TEMP) + KE(K,J) !*2 !&&&&&&&&&&&&&&&&&&&&&&&2bei....OK
	      END IF
	    END DO
	  END DO         
        DEALLOCATE(GG)
	END DO

      END SUBROUTINE FORM_K
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      SUBROUTINE FORM_ORDER(T_NODE,M)
	!CALCULATE THE ORDER OF ELEMENT IN TOTAL STIFF MATRIX
	USE DEF
      IMPLICIT NONE
	INTEGER(4)::T_NODE(4,2),M
	!ABOVE TO TRANSFER
	INTEGER(4)::I
      T_NODE=0
      DO I = 1,4
        T_NODE(I,1) = G_ELE(M,I)*2-1
	  T_NODE(I,2) = G_ELE(M,I)*2
	END DO
	END SUBROUTINE FORM_ORDER       
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      SUBROUTINE FORM_LOAD()
      !FORM LOAD
      USE DEF
      IMPLICIT NONE
	!ABOVE TO TRANSFER
      EXTERNAL::PP_LOAD,G_LOAD,EQUAL_VERTICAL_LINE_LOAD

      IF (N_LOAD .GT. 0) THEN
      CALL PP_LOAD()
	END IF
     
	IF (ROU .GT. 0.0D0) THEN
      CALL G_LOAD()
	END IF
	
	IF (N_LINE_LOAD .GT. 0) THEN
      CALL EQUAL_VERTICAL_LINE_LOAD()
	END IF

      IF (X_N_LOAD .GT. 0) THEN
	CALL EQUAL_X_LOAD()
	END IF

      IF (Q_LINEAR .NE. 0) THEN
      CALL LINEAR_LOAD()
      END IF

	END SUBROUTINE FORM_LOAD
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      SUBROUTINE PP_LOAD()
      !FORM THE LOAD ARRAY OF POINT LOAD
      USE DEF
      IMPLICIT NONE
      !ABOVE TO TRANSFER
	INTEGER(4)::I,GI,GJ

      VV=0.0
	DO I = 1,N_LOAD
	  GI = 2*NP_LOAD(I)-1
	  GJ = 2*NP_LOAD(I)
	  VV(GI) = P_LOAD(I,1)+VV(GI)
	  VV(GJ) = P_LOAD(I,2)+VV(GJ)
	END DO
                
      END SUBROUTINE PP_LOAD
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	SUBROUTINE G_LOAD()
	!FORM THE ARRAY OF GRAVITY LOAD
      USE DEF
      IMPLICIT NONE
       !ABOVE TO TRANSFER
	INTEGER(4)::I,J,K,TT
	REAL(8)::TEMP(2,2),GG
	REAL(8)::T1,T2,T3,T4
	PARAMETER(GG=9.80)
	EXTERNAL::CAL_POINT


	DO I = 1,N_ELE
      T1=0.0
	T2=0.0
	T3=0.0
	T4=0.0
        DO J = 1,4
	    COORD(J,1) = G_NODE(G_ELE(I,J),1)
	    COORD(J,2) = G_NODE(G_ELE(I,J),2)
	  END DO
        IF (G_ELE(I,3) .EQ. G_ELE(I,4)) THEN
	     S34 = 3
	     CALL CAL_POINT()
	     IF (N_INT .EQ. 1) THEN
	       TT=1
	     ELSE IF(N_INT .EQ. 2) THEN
	       TT = 3
	     ELSE IF(N_INT .EQ. 3) THEN
	       TT = 4
	     ELSE
	       PRINT*,'ERROR N_INT IN SUBROUTINE G_LOAD'
	     END IF
           DO K = 1,TT
	       KSI = POINT3(K,2)
	       YITA = POINT3(K,3)
	       CALL CAL_N()
             VV(G_ELE(I,1)*2)=VV(G_ELE(I,1)*2)-
     5		   H3(K)*ROU*GG*PT*N(1)
             VV(G_ELE(I,2)*2)=VV(G_ELE(I,2)*2)-
     5		   H3(K)*ROU*GG*PT*N(2)
             VV(G_ELE(I,3)*2)=VV(G_ELE(I,3)*2)-
     5		   H3(K)*ROU*GG*PT*N(3)
	     END DO
	  ELSE
	     S34 = 4
		 CALL CAL_POINT()
           DO J = 1,N_INT
	       DO K = 1,N_INT
	         YITA = POINT(J)
	         KSI = POINT(K)
               CALL CAL_N()
         		 T1 = T1+H(K)*H(J)*PT*ROU*GG*N(1)
     	         T2 = T2+H(K)*H(J)*PT*ROU*GG*N(2)
     			 T3 = T3+H(K)*H(J)*PT*ROU*GG*N(3)
        		 T4 = T4+H(K)*H(J)*PT*ROU*GG*N(4)
	       END DO
	     END DO
           VV(G_ELE(I,1)*2) = VV(G_ELE(I,1)*2)-T1
           VV(G_ELE(I,2)*2) = VV(G_ELE(I,2)*2)-T2
           VV(G_ELE(I,3)*2) = VV(G_ELE(I,3)*2)-T3
           VV(G_ELE(I,4)*2) = VV(G_ELE(I,4)*2)-T4
	  END IF

      END DO

      END SUBROUTINE G_LOAD
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      SUBROUTINE EQUAL_VERTICAL_LINE_LOAD()
      !CALCULATE THE LOAD (VERTICAL LOAD ON LINE EQUALLY)
      USE DEF
      IMPLICIT NONE
	!ABOVE TO TRANSFER
	INTEGER(4)::I,J,T(1:2),GI,GJ
	REAL(8)::Q,YI,YJ,XI,XJ

      DO I = 1,N_LINE_LOAD
	  Q = LINE_LOAD(I)
	  T(1) = LINE_LOAD_NODE(I,1)
	  T(2) = LINE_LOAD_NODE(I,2)
	  XI = G_NODE(T(1),1)
	  XJ = G_NODE(T(2),1)
        YI = G_NODE(T(1),2)
	  YJ = G_NODE(T(2),2)
        DO J = 1,2
          GI = 2*T(J)-1
	    GJ = 2*T(J)
          VV(GI) = VV(GI)+Q*PT*(YI-YJ)/2.0
	    VV(GJ) = VV(GJ)+Q*PT*(XJ-XI)/2.0
	  END DO
      END DO

      END SUBROUTINE EQUAL_VERTICAL_LINE_LOAD
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      SUBROUTINE EQUAL_X_LOAD()
      !CALCULATE THE LOAD OF KIND
      USE DEF
      IMPLICIT NONE
	!ABOVE TO TRANSFER
	INTEGER(4)::T1,T2,I
	REAL(8)::XI,XJ,YI,YJ,DIST,Q

	DO I = 1,X_N_LOAD
	  Q = X_LOAD(I)
	  T1 = X_LOAD_NODE(I,1)
	  T2 = X_LOAD_NODE(I,2)
        XI = G_NODE(T1,1)
	  XJ = G_NODE(T2,1)
	  YI = G_NODE(T1,2)
	  YJ = G_NODE(T2,2)
	  DIST = SQRT((XI-XJ)**2+(YI-YJ)**2)
        VV(2*T1-1) = VV(2*T1-1)+Q*PT*DIST/2.0
	  VV(2*T2-1) = VV(2*T2-1)+Q*PT*DIST/2.0
      END DO

      END SUBROUTINE EQUAL_X_LOAD
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      SUBROUTINE LINEAR_LOAD()
      !CALCULATE THE LINEARLY LOAD ON ONE VERTICLE LINE
      USE DEF
      IMPLICIT NONE
	!ABOVE TO TRANSFER
	INTEGER(4)::I,J,NN,T(1:2),GI,GJ
	REAL(8)::Q,XI,XJ,YI,YJ
      !N_LINEAR(1:3):1.MINIMAL NUMBER OF STARTING POINT WITH THE LINEAR FORCE
	!              2.MAXIMUM NUMBER OF ENDDING POINT WITH THE LINEAR FORCE
	!Q_LINEAR:THE MAXIMUM VALUE OF THE FORCE
	!        Q_LINEAR .GT. 0: PRESSURE
	!                 .LT. 0: PULL
      
	Q = Q_LINEAR
	NN = N_LINEAR(2)-N_LINEAR(1)
      DO I = 1,NN
	  T(1) = N_LINEAR(1)+(I-1)
        T(2) = T(1)+1
        XI = G_NODE(T(1),1)
        XJ = G_NODE(T(2),1)
        YI = G_NODE(T(1),2)
        YJ = G_NODE(T(2),2)
        GI = 2*T(1)-1
        GJ = 2*T(1)
        VV(GI) = VV(GI)+Q*PT*(YI-YJ)*(1/3.0+(I-1)/NN/2.0)
	  VV(GJ) = VV(GJ)+Q*PT*(XJ-XI)*(1/3.0+(I-1)/NN/2.0)
        GI = 2*T(2)-1
        GJ = 2*T(2)
        VV(GI) = VV(GI)+Q*PT*(YI-YJ)*(1/6.0+(I-1)/NN/2.0)
	  VV(GJ) = VV(GJ)+Q*PT*(XJ-XI)*(1/6.0+(I-1)/NN/2.0)    
	END DO
      END SUBROUTINE LINEAR_LOAD
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      SUBROUTINE ADD_BC()
      !ADD THE BOUNDARY CONDITION
      USE DEF
      IMPLICIT NONE
	!ABOVE TO TRANSFER
	INTEGER(4)::I,J,GI,GJ
	REAL(8)::ALPHA

	
	ALPHA = 1.0E10

      DO I = 1,N_BC
	  GI = 2*N_DIS(I)-1
	  GJ = 2*N_DIS(I)  
	  IF (BC(I,1) .NE. 1.0) THEN
          KK(DIAG(GI)) = ALPHA*KK(DIAG(GI))
	    VV(GI) = ALPHA*BC(I,1)*KK(DIAG(GI))
	  END IF
	  IF (BC(I,2) .NE. 1.0) THEN
          KK(DIAG(GJ)) = ALPHA*KK(DIAG(GJ))
	    VV(GJ) = ALPHA*BC(I,2)*KK(DIAG(GJ))
	  END IF
      END DO

	END SUBROUTINE ADD_BC
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      SUBROUTINE SOLVE(A,P,M,D,NUM)
      !一维变带宽矩阵求解方程组
	USE DEF
      INTEGER(4)::NUM
	REAL(8)::A(DOF)
	REAL(8)::D(NUM)
	INTEGER(4)::P(NUM),M(NUM+1)
      !以上变量需要传入
      INTEGER(4)::I,J,Q
	REAL(8)::S,T

      !用一维数组A存储上三角的非零元素
      !用一维数组m储存主对角元在一维数组中的位置
      !用一维数组N储存列高
      !用一维数组p储存每一列的起始行号

      !GUASS分解a
      IF (1 .EQ. 0) THEN
	DO J=1,NUM
	  DO I = 1,J
	    S = A(M(J)+J-I)
	    DO Q = MAX(P(I),P(J)), I-1
            S=S-A(M(I)+I-Q)*A(M(J)+J-Q)/A(M(Q))
	    END DO
	    A(M(J)+J-I) = S
	  END DO
	END DO


      !求解方程
      !D是最后求出来的解
 
      DO J =1,NUM
	  S = D(J)
	  DO I = P(J),J-1
	    S = S-A(M(J)+J-I)*D(I)/A(M(I))
	  END DO
	  D(J) = S
	END DO

      D(NUM) = D(NUM)/A(M(NUM))
	DO I = NUM-1,1,-1
	  S = D(I)
	  DO J = I+1,NUM
	    IF (P(J) .LE. I) THEN
	      S = S-A(M(J)+J-I)*D(J)
	    END IF
	  END DO
	  D(I) = S/A(M(I))
	END DO

      END IF

      !LDLT
      IF (1 .EQ. 1) THEN

	DO K = 1,NUM
	  DO J = P(K),K-1
          A(M(K)) = A(M(K)) - A(M(K)+K-J)**2*A(M(J))
	  END DO
	  DO I = K+1,NUM
	  IF (P(I) .LE .K) THEN
          DO J = MAX(P(K),P(I)),K-1
            A(M(I)+I-K)=A(M(I)+I-K)-A(M(I)+I-J)*A(M(K)+K-J)*A(M(J))
	    END DO
	    A(M(I)+I-K)=A(M(I)+I-K)/A(M(K))
        END IF
	  END DO
      END DO

      DO K = 2,NUM
	  DO J = P(K),K-1
          D(K) = D(K) - A(M(K)+K-J)*D(J)
	  END DO
      END DO

      D(NUM)=D(NUM)/A(M(NUM))
	DO K = NUM-1,1,-1
	  T = 0.0D0
        DO J = K+1,NUM
	    IF (P(J) .LE. K) THEN
          T = T+A(M(J)+J-K)*D(J)
	    END IF
        END DO
	D(K) = D(K)/A(M(K))-T
      END DO
      END IF
 


      DO I = 1,NUM
	  IF (ABS(D(I)) .LT. 1.0E-16) THEN
	    D(I) = 0.0
	  END IF
	END DO

      END SUBROUTINE SOLVE
	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	SUBROUTINE CAL_STR()
      !CALCULATE THE STRESS
      USE DEF
      IMPLICIT NONE
      !ABOVE TO TRANSFER
	INTEGER(4)::I,J,K,P,Q
	REAL(8)::DD(3,3),BB(3,8),TEMP1(3,6)
	REAL(8)::TEMP2(1:8),TEMP3(1:3)
      REAL(8)::TJACO(2,2),QQ
	EXTERNAL::MULTI,CAL_BB,CAL_DD
      
	SIGMA=0.0
	DO I = 1,N_ELE
	   BB = 0.0
	   DD = 0.0
	   DO J = 1,4
	     COORD(J,1) = G_NODE(G_ELE(I,J),1)
	     COORD(J,2) = G_NODE(G_ELE(I,J),2)
         END DO 
	   IF (G_ELE(I,3) .EQ. G_ELE(I,4)) THEN
	     S34 = 3
	     KSI = 0.0
	     YITA = 0.0
	     CALL CAL_BB(BB)
	     CALL CAL_DD(DD)
	     CALL MULTI(DD,BB(:,1:6),TEMP1,3,3,6)
	     DO J = 1,S34
	       TEMP2(2*J-1) = VV(2*G_ELE(I,J)-1)
	  	   TEMP2(2*J) = VV(2*G_ELE(I,J))
		 END DO 
	     CALL MULTI(TEMP1,TEMP2(1:6),TEMP3,3,6,1)
	     CALL CAL_DN()
	     CALL MULTI(DN,COORD,JACOBI,2,4,2)
           CALL CAL_DET()
	     DO J = 1,S34
	      SIGMA(G_ELE(I,J),1)=SIGMA(G_ELE(I,J),1)+TEMP3(1)*DET/2.0
	      SIGMA(G_ELE(I,J),2)=SIGMA(G_ELE(I,J),2)+TEMP3(2)*DET/2.0
	      SIGMA(G_ELE(I,J),3)=SIGMA(G_ELE(I,J),3)+TEMP3(3)*DET/2.0
	      SIGMA(G_ELE(I,J),4)=SIGMA(G_ELE(I,J),4)+DET/2.0
	     END DO
	   ELSE 
	     S34 = 4
	     DO K =1,S34
             TEMP2(K*2-1) = VV(2*G_ELE(I,K)-1)
	       TEMP2(K*2) = VV(2*G_ELE(I,K))
	     END DO
	     DO J = 1,S34
		   CALL CAL_VEX()  
             KSI = VEX(1,J)
	       YITA = VEX(2,J)
	       CALL CAL_BB(BB)
	       CALL CAL_DD(DD)
	       CALL MULTI(DD,BB,BB,3,3,8)
	       CALL MULTI(BB,TEMP2,TEMP3,3,8,1)
		   CALL CAL_POINT()
		   QQ=0.0
		   DO Q = 1,N_INT
	         DO P = 1,N_INT
	           YITA = POINT(P)
	           KSI = POINT(Q)
	           CALL CAL_DN()
	           CALL MULTI(DN,COORD,JACOBI,2,4,2)
	           CALL CAL_DET()
	           QQ = QQ+DET*H(P)*H(Q)
	         END DO
	       END DO
	      SIGMA(G_ELE(I,J),1)=SIGMA(G_ELE(I,J),1)+TEMP3(1)*QQ
	      SIGMA(G_ELE(I,J),2)=SIGMA(G_ELE(I,J),2)+TEMP3(2)*QQ
	      SIGMA(G_ELE(I,J),3)=SIGMA(G_ELE(I,J),3)+TEMP3(3)*QQ
	      SIGMA(G_ELE(I,J),4)=SIGMA(G_ELE(I,J),4)+QQ
	     END DO
	   END IF
      END DO
	DO I = 1,N_NODE
        SIGMA(I,1) = SIGMA(I,1)/SIGMA(I,4)
        SIGMA(I,2) = SIGMA(I,2)/SIGMA(I,4)
        SIGMA(I,3) = SIGMA(I,3)/SIGMA(I,4)
      END DO

      END SUBROUTINE CAL_STR
	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      PROGRAM MAIN    

      USE DEF
      IMPLICIT NONE
      INTEGER(4)::I,J,K

      EXTERNAL::READ_DATA,CAL_PE,BAND_K,FORM_LOAD,ADD_BC,SOLVE,CAL_STR
	EXTERNAL::ELE_TRANS,FORM_K

      OPEN(UNIT=1,FILE='D:\MY PROGRAM IN FORTRAN\A\BASIC.TXT')
	OPEN(UNIT=8,FILE='D:\MY PROGRAM IN FORTRAN\A\DATA.TXT')
	!OPEN(UNIT=9,FILE='D:\MY PROGRAM IN FORTRAN\A\TEST.TXT')

	READ(1,*) ID,N_NODE,N_ELE,N_INT,PE,PT,PR,ROU
     5      ,N_BC,N_LOAD,N_LINE_LOAD,X_N_LOAD,N_LINEAR(1:2),Q_LINEAR
	

      ALLOCATE(G_ELE(1:N_ELE,4),MARK(1:N_ELE),G_NODE(1:N_NODE,2))
	ALLOCATE(LINE_LOAD_NODE(1:N_LINE_LOAD,2),P_LOAD(1:N_LOAD,2))
	ALLOCATE(X_LOAD_NODE(1:X_N_LOAD,2),X_LOAD(1:X_N_LOAD)
	ALLOCATE(NP_LOAD(1:N_LOAD),LINE_LOAD(1:N_LINE_LOAD))
	ALLOCATE(N_DIS(1:N_BC),BC(1:N_BC,2))
	ALLOCATE(VV(1:N_NODE*2),DIAG(2*N_NODE+1),PP(2*N_NODE))
	ALLOCATE(SIGMA(1:N_NODE,1:4))


      !READ DATA FROM DATA.TXT
	CALL READ_DATA()
	CALL CAL_PE()

	!FORM THE MATRIX OF ELEMENT K
      !CALL FORM_KE(S34,N_INT,PT,PR,PE,COORD,KE)

	!FORM THE MATRIX OF TOTAL K
	CALL BAND_K()
	ALLOCATE(KK(DIAG(2*N_NODE+1)-1))
      KK=0.0
	DOF = DIAG(2*N_NODE+1)-1      
	CALL FORM_K()

	!FORM THE SCALAR OF LOAD
	VV=0.0
      CALL FORM_LOAD()

	!ADD THE BOUNDARY CONDITION OF DISPLACEMENT
      CALL ADD_BC()

	!SOLVE THE EQUATION
	CALL SOLVE(KK,PP,DIAG,VV,N_NODE*2)

	!CALCULATE THE STRESS
      CALL CAL_STR()

      WRITE(8,*) 'DISPLACEMENT OF EACH NODE'
	DO I = 1,2*N_NODE,2
        WRITE(8,*) (I+1)/2,':',VV(I),VV(I+1)
      END DO
      WRITE(8,*) 
	WRITE(8,*) 'SIGMA FOR EACH NODE AFTER AVERAGING'
	DO I = 1,N_NODE
        WRITE(8,*) I,':',SIGMA(I,1),SIGMA(I,2),SIGMA(I,3)
	END DO
      WRITE(8,*) 


      CLOSE(1)
	CLOSE(6)
      CLOSE(8)

	END PROGRAM MAIN