eq.for

静力学有限元源程序

	    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 FORM_KE'
	  END IF
        DO K = 1,TT
	    KSI = POINT3(K,1)
	    YITA = POINT3(K,2)
	    CALL CAL_DN()
          CALL CAL_BB(BB)
	    CALL CAL_DD(DD)
	    CALL MULTI(DD,BB,TEMP1,3,3,16)
	    CALL MULTI(TRANSPOSE(BB),TEMP1,TEMP2,16,3,16)
	    DO P = 1,16
	      DO Q = 1,16
              KE(P,Q) = KE(P,Q) + TEMP2(P,Q)*H3(K)*PT*DET
	      END DO
	    END DO
	  END DO
	ELSE IF (MTYPE .EQ. 8) THEN
	  CALL CAL_POINT()
	  DO I = 1,N_INT
	    DO J = 1,N_INT
	      YITA = POINT(I)
	      KSI = POINT(J)
	      CALL CAL_DN()
            CALL CAL_BB(BB)
	      CALL CAL_DD(DD)
	      CALL MULTI(DD,BB,TEMP1,3,3,16)
	      CALL MULTI(TRANSPOSE(BB),TEMP1,TEMP2,16,3,16)
	      DO P = 1,16
	        DO Q = 1,16
                KE(P,Q) = KE(P,Q) + PT*TEMP2(P,Q)*H(J)*H(I)*DET
	        END DO
	      END DO
	    END DO
	  END DO

	END IF

	END SUBROUTINE FORM_KE
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      SUBROUTINE FIND_34()
      USE DEF
	IMPLICIT NONE
	!ABOVE TO TRANSFER
	INTEGER(4)::I,J,K,Q


      DEVIDE = 0
	DO I = 1,N_ELE
	  IF ((MARK(I) .GT. 5) .AND. 
     5     (I .NE. ELE)) THEN
	    IF (MTYPE .EQ. 3) THEN
	      DO J = 1,MARK(I)
	        DO Q = 1,2
                IF (G_ELE(ELE,Q) .EQ. G_ELE(I,J)) THEN
	            DO K = 1,MARK(I)
	              IF(G_ELE(ELE,Q+1) .EQ. G_ELE(I,K)) THEN
	                DEVIDE(Q) =  1
	                IF (J .LT. K) THEN
	                  G_ELE(ELE,3+Q)=G_ELE(I,J+3)
	                ELSE IF (J .GT. K) THEN
	                  G_ELE(ELE,3+Q)=G_ELE(I,K+3)
	                END IF
	              END IF
                  END DO
                END IF
	        END DO
              IF (G_ELE(ELE,3) .EQ. G_ELE(I,J)) THEN
	          DO K = 1,MARK(I)
	            IF(G_ELE(ELE,1) .EQ. G_ELE(I,K)) THEN
	              DEVIDE(3) =  1
	              IF (J .LT. K) THEN
	                G_ELE(ELE,6)=G_ELE(I,K+3)
	              ELSE IF (J .GT. K) THEN
	                G_ELE(ELE,6)=G_ELE(I,J+3)
	              END IF
	            END IF
                END DO
              END IF
            END DO
          ELSE IF (MTYPE .EQ. 4) THEN
	      DO J = 1,MARK(I)
	        DO Q = 1,3
                IF (G_ELE(ELE,Q) .EQ. G_ELE(I,J)) THEN
	            DO K = 1,MARK(I)
	              IF(G_ELE(ELE,Q+1) .EQ. G_ELE(I,K)) THEN
	                DEVIDE(Q) =  1
	                IF (J .LT. K) THEN
	                  G_ELE(ELE,4+Q)=G_ELE(I,J+4)
	                ELSE IF (J .GT. K) THEN
	                  G_ELE(ELE,4+Q)=G_ELE(I,K+4)
	                END IF
	              END IF
                  END DO
                END IF
	        END DO
              IF (G_ELE(ELE,4) .EQ. G_ELE(I,J)) THEN
	          DO K = 1,MARK(I)
	            IF(G_ELE(ELE,1) .EQ. G_ELE(I,K)) THEN
	              DEVIDE(4) =  1
	              IF (J .LT. K) THEN
	                G_ELE(ELE,8)=G_ELE(I,K+4)
	              ELSE IF (J .GT. K) THEN
	                G_ELE(ELE,8)=G_ELE(I,J+4)
	              END IF
	            END IF
                END DO
              END IF
            END DO	  
	    END IF
	  END IF
	END DO

	END SUBROUTINE FIND_34
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      SUBROUTINE CAL_BB(BB)
      !CALCULATE THE STRAIN-DISPLACEMENT MATRIX OF ELEMENT
      USE DEF
      IMPLICIT NONE
	REAL(8)::BB(3,16)
	!ABOVE TO TRANSFER
	INTEGER(4)::I,J
	!BB(3,16): STRAIN-DISPLACEMENT MATRIX OF ELEMENT
	!  FOR TRIANGLE ELEMENT, THE LAST TWO COLUMN IS ZEROS

      BB = 0.0D0
	CALL MULTI(DN,COORD,JACOBI,2,8,2)
	CALL CAL_DET()
	CALL CAL_INVJ()
	DO I = 1,8
        BB(1,2*I-1) = JACOBI(1,1)*DN(1,I)+JACOBI(1,2)*DN(2,I)
	  BB(2,2*I) = JACOBI(2,1)*DN(1,I)+JACOBI(2,2)*DN(2,I)
	  BB(3,2*I-1) = BB(2,2*I)
	  BB(3,2*I) = BB(1,2*I-1)
	END DO

	END SUBROUTINE CAL_BB
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      SUBROUTINE CAL_DD(DD)
	!CALCULATE THE STRAIN-STESS MATRIX OF ELEMENT
	USE DEF
      IMPLICIT NONE
	REAL(8)::DD(3,3)
	!ABOVE TO TRANSFER

      DD(1,:)=(/1.0D0,PR,0.0D0/)
	DD(2,:)=(/PR,1.0D0,0.0D0/)
	DD(3,:)=(/0.0D0,0.0D0,(1.0D0-PR)/2.0D0/)
	DD = DD*PE/(1.0D0-PR*PR)
 
	END SUBROUTINE CAL_DD
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      SUBROUTINE CAL_POINT()
      !CALCULATE THE POINTS TO BE INTEGRATED
      USE DEF
	!ABOVE TO BE TRANSFER
	!INTEGER(4)::
	!REAL(4)::
      !N_INT FOR S34=3, AND S34=4 HAVE DIFFERENT MEANING
      !FOR S34=3, N_INT MEANS THE ORDERS 
	!FOR S34=4, N_INT MEANS THE POINTS	                     
      IMPLICIT NONE
      POINT = 0.0
	POINT3 = 0.0
	H=0.0
	H3=0.0
	IF ((MTYPE .EQ. 3) .OR. (MTYPE .EQ. 6)) THEN
	  IF (N_INT .EQ. 1) THEN
          POINT3(1,:) = (/1/3.0,1/3.0,1/3.0/)
	    H3(1) = 1.0/2.0
	  ELSE IF(N_INT .EQ. 2) THEN	    
          POINT3(1,:) = (/2.0/3.0,1.0/6.0,1.0/6.0/)
	    POINT3(2,:) = (/1.0/6.0,2.0/3.0,1.0/6.0/)
	    POINT3(3,:) = (/1.0/6.0,1.0/6.0,2.0/3.0/)
	    H3(1:3) = 1.0/3.0/2.0
	  ELSE IF(N_INT .EQ. 3) THEN
          POINT3(1,:) = (/1/3.0,1/3.0,1/3.0/)
	    POINT3(2,:) = (/0.6,0.2,0.2/)
	    POINT3(3,:) = (/0.2,0.6,0.2/)
	    POINT3(4,:) = (/0.2,0.2,0.6/)
	    H3(1) = -27.0/48.0/2.0 
	    H3(2:4) = 25.0/48.0/2.0
	  END IF
	ELSE IF ((MTYPE .EQ. 4) .OR. (MTYPE .EQ. 8)) THEN
	  IF (N_INT .EQ. 1) THEN
	    POINT=0.0 
	    H=2.0
	  ELSE IF (N_INT .EQ. 2) THEN
	    POINT(1) = 1.0/SQRT(3.0)
	    POINT(2) = -POINT(1)
	    H(1:2)=1.0
	  ELSE IF (N_INT .EQ. 3) THEN
	    POINT(1) =  0.0
	    POINT(2) = SQRT(15.0)/5.0
	    POINT(3) = -POINT(2)
	    H(1) = 8.0/9.0
	    H(2:3) = 5.0/9.0
	  ELSE 
	    PRINT*,'WRONG NUMBER FOR GAUSS INTEGRATE IN 
     5		SUBROUTINE CAL_POINT'
	  END IF
	ELSE
	  PRINT*,'WRONG ELEMENT TYPE IN SUBROUTINE CAL_POINT'
	END IF

      END SUBROUTINE CAL_POINT
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      SUBROUTINE BAND_K()     
	!GET THE LENGHTH OF STIFF MATRIX KK
      USE DEF
      IMPLICIT NONE
	!ABOVE TO TRANSFER
	INTEGER(4)::TN,T_NODE(8,2),I,J,K,M,GJ,GI,T,MIN
	EXTERNAL::FORM_ORDER

      !PP(1:N_NODE*2) : THE MINIMAL ROW NUMBER OF EACH COLUMN
      !DIAG(1:N_NODE*2+1) : DIAG ELEMENT IN 1D ARRAY

      TN = N_NODE*2
      DO I = 1,TN
        PP(I) = I
      END DO
      
	!TWO METHODS TO FIND THE PP(1:2*N_NODE)
      IF(1 .EQ. 1) THEN
      DO I = 1,N_ELE
	!FIND THE MINIMAL ROW NUMBER OF EACH 
	  ELE = I
	  MTYPE = MARK(I)
        CALL FORM_ORDER(T_NODE,ELE)

	  MIN = T_NODE(1,1)
	  DO J = 2,MTYPE
	    IF (T_NODE(J,1) .LT. MIN) THEN
            MIN = T_NODE(J,1)
	    END IF
	  END DO

        DO J = 1,MTYPE        
	    GI = T_NODE(J,1)
          GJ = T_NODE(J,2)
	    IF (PP(GI) .GT. MIN) THEN
	  	  PP(GI) = MIN
	    END IF
	    IF (PP(GJ) .GT. MIN) THEN
	   	  PP(GJ) = MIN
	    END IF
	  END DO
      END DO
	END IF

      IF( 1 .EQ. 0) THEN
	DO I = 1,2*N_NODE
	  DO J = 1,N_ELE
	    ELE = J
	    MTYPE = MARK(ELE)
          CALL FORM_ORDER(T_NODE,ELE)
	    DO K = 1,MTYPE
	      IF ((T_NODE(K,1) .EQ. I) .OR. (T_NODE(K,2) .EQ. I))THEN
	        DO M = 1,MTYPE
                GI = T_NODE(M,1)
			  GJ = T_NODE(M,2)              
	          IF (PP(GI) .GT. I) THEN
	            PP(GI) = I
	          END IF
	          IF (PP(GJ) .GT. I) THEN
	            PP(GJ) = I
	          END IF	          
	        END DO
	      END IF
	    END DO
	  END DO
	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,TEMP
	INTEGER(4)::GG(1:16)
      EXTERNAL::FORM_KE
       
      DO I = 1,N_ELE
        MTYPE = MARK(I)
	  ELE = I
        COORD=0.0
        DO J = 1,MTYPE
	    IF (G_ELE(I,J) .NE. 0) THEN
	      COORD(J,1) = G_NODE(G_ELE(I,J),1)
	      COORD(J,2) = G_NODE(G_ELE(I,J),2)   !RIGHT
	    END IF
	  END DO
 
        
	  CALL FORM_KE()
	  AREA(ELE) = DET
	  !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*
        GG = 0
	  DO J = 1,8
	    IF (G_ELE(I,J) .NE. 0) THEN
	      GG(2*J-1)=2*G_ELE(I,J)-1
	      GG(2*J)=2*G_ELE(I,J)
	    END IF
        END DO
	  DO J =1,16
	    DO K = 1,16
	IF ((GG(K) .NE. 0) .AND. (GG(J) .NE. 0)) THEN
	      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 IF
	    END DO
	  END DO
	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(8,2),M
	!ABOVE TO TRANSFER
	INTEGER(4)::I
      T_NODE=0
      DO I = 1,MTYPE
        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
      
	VV=0.0
      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


	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,P,Q,TT
	REAL(8)::TEMP(4),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
	MTYPE = MARK(I)
	ELE = I
	TEMP = 0.0
	COORD = 0.0
        DO J = 1,8
	    IF (G_ELE(I,J) .NE. 0) THEN
	      COORD(J,1) = G_NODE(G_ELE(I,J),1)
	      COORD(J,2) = G_NODE(G_ELE(I,J),2)
	    END IF
	  END DO
        IF (MTYPE .EQ. 3) THEN
	    CALL FIND_34()
	    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
	    CALL CAL_POINT()
          DO K = 1,TT
	      KSI = POINT3(K,1)
	      YITA = POINT3(K,2)
		  CALL CAL_N()
	      CALL  COMPUT()
	      DO J = 1,MTYPE
              IF (DEVIDE(J) .EQ. 1) THEN
                N(J) = N6(J)
                N(J) = N6(J)
                N(3+J) = N6(3+J)
                N(3+J) = N6(3+J)
              END IF
            END DO
            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)
	      DO J = 1,MTYPE
	        IF (DEVIDE(J) .EQ. 1) THEN
                VV(G_ELE(ELE,J+3)*2)=VV(G_ELE(ELE,J+3)*2)-
     5	        H3(K)*ROU*GG*PT*N(3+J)
              END IF
	      END DO
	    END DO
	  ELSE IF (MTYPE .EQ. 4) THEN
	     CALL FIND_34()
		 CALL CAL_POINT()
           DO J = 1,N_INT
	       DO K = 1,N_INT
	         YITA = POINT(J)
	         KSI = POINT(K)
	  	     CALL CAL_N()
	         CALL  COMPUT()
	         DO P = 1,MTYPE
                 IF (DEVIDE(P) .EQ. 1) THEN
                   N(P) = N8(P)
                   N(P) = N8(P)
                   N(4+P) = N8(4+P)
                   N(4+P) = N8(4+P)
                 END IF
               END DO
         		 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)
	         DO P = 1,MTYPE
	           IF (DEVIDE(P) .EQ. 1) THEN
                   TEMP(P)=TEMP(P)+H(J)*H(K)*ROU*GG*PT*N(4+J)
                 END IF
	         END DO
	       END DO