eq.for

静力学有限元源程序

      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	!PROGRAM FOR FINITE ELEMENT METHOD
	!USED ONLY FOR PLANE PROBLEM
	!PROBLEM STYLE: PLANE STRESS(ID=1) OR PLANE STRAIN(ID=2)
	!ELEMENT TYPE: TRIANGLE AND QUADRANGLE
	!LOAD TYPE: GRAVITY, LOAD ON NODE
	!B.C.: DISPLACEMENT CONSTRAINT ONLY
      !KEY PARAMETER:
	!   ID: PROBLEM INDICATOR(1 FOR PLANE STRESS,2 FOR PLANE STRAIN)
	!   N_NODE: NUMBER OF NODE
	!   N_ELE: NUMBER OF ELEMENT
	!   N_LOAD: NUMBER OF LOAD
	!   N_BC: NUMBER OF DISPLACEMENT
	!   N_INT:NUMBERS TO INTEGRATE  ONLY FROM 1 TO 3 POINTS
	!   PE: ELASTIC MODULE
	!   PT: THICKNESS
	!   PR: POISSON RATIO
	!   ROU: ROUSITY
	!
	!KEY ARRAY:
	!   X(1:N_NODE),Y(1:NODE): X AND Y COORD OF NODE
      !   G_NODE(1:N_NODE,2): 
	!                       FIRST COLUMN REPRESENT THE X COORD OF NODE
	!                       SECOND COLUMN REPRESENT THE Y COORD OF NODE
	!   
	!   G_ELE(1:N_ELE,4): THE NODE NUMBER IN EACH ELEMENT
	!   KK(:): TOTAL STIFF MATRIX, ALLOCATABLE ONE RANK
	!   KE(8,8): STIFF MATRIX IN ELEMENT
      !   PP(1:N_NODE*2) : THE MINIMAL ROW NUMBER OF EACH COLUMN
		  !TWO METHODS TO FIND THE PP(1:2*N_NODE)
      !   DIAG(1:N_NODE*2+1) : DIAG ELEMENT IN 1D ARRAY
	!   VV(1:N_NODE*2): THE ARRAY OF DISPLACEMENT RESULT
      !   P_LOAD(1:N_LOAD,2):FIRST COLUMN REPRESENT THE X FORCE OF NODE
	!                      SECOND COLUMN REPRESENT THE Y FORCE OF NODE
	!   BC(1:N_BC,2):FIRST COLUMN REPRESENT THE X DISPLACEMENT OF NODE
	!                SECOND COLUMN REPRESENT THE Y DISPLACEMENT OF NODE
	!                1 REPRESENT THAT THE DISPLACEMENT IS FREE
	!                OTHERS REPRESENT THAT THE DISPLACEMENT IS BC(1:N_BC,1:2)
	!   LINE_LOAD(1:N_LINE_LOAD):REPRESENT THE Q IN DIFFERENT LINE
	!                WHEN Q .GT. 0, FORCE IS PRESSURE
	!                       .LT. 0, FORCE IS PULL
      !   LINE_LOAD_NODE(1:N_LINE_LOAD,2):GIVE THE NUMBER OF THE NODE BETWEEN ONE LINE
      !   X_LOAD(1:X_N_LOAD):GIVE THE Q IN DIFFERENT LINE
	!                WHEN Q .GT. 0, FORCE IS POSSITIVE IN X-COORD
	!                       .LT. 0, FORCE IS NEGATIVE IN X-COORD
	!   X_LOAD_NODE(1:X_N_LOAD,2):GIVE THE NUMBER OF THE NODE BETWEEN ONE LINE
      !   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
	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      MODULE DEF
	IMPLICIT NONE

	REAL(8)::COORD(8,2),DET,JACOBI(2,2),N(1:8),DN(2,8),KSI,YITA
!READ_DATA      
      INTEGER(4)::N_NODE,N_ELE,N_LOAD,N_BC,N_LINE_LOAD,X_N_LOAD
      INTEGER(4),ALLOCATABLE::MARK(:)
	INTEGER(4)::MTYPE,ELE      

!CAL_DD
      INTEGER(4)::ID
	REAL(8)::PT,PE,PR

!FORM_KE
      INTEGER(4)::N_INT
	REAL(8)::KE(16,16)
	REAL(8)::POINT(1:3),H(1:3),POINT3(4,3),H3(1:4)
!BAND_K

!FORM_K
      INTEGER(4)::DOF,DEVIDE(1:4)
	REAL(8),ALLOCATABLE::AREA(:)
!FORM_LOAD
	INTEGER(4)::N_LINEAR(1:2)
	REAL(8)::ROU,Q_LINEAR
!ADD_BC
!CAL_STR
      !REAL(8)::

!PARAMETER
      REAL(8)::DN6(2,6),DN8(2,8),N6(1:6),N8(1:8)
!MAIN
	REAL(8),ALLOCATABLE::KK(:)
      INTEGER(4),ALLOCATABLE::PP(:),N_DIS(:),NP_LOAD(:),X_LOAD_NODE(:,:)
	INTEGER(4),ALLOCATABLE::G_ELE(:,:),DIAG(:),LINE_LOAD_NODE(:,:)
	REAL(8),ALLOCATABLE::G_NODE(:,:),P_LOAD(:,:),VV(:)
	REAL(8),ALLOCATABLE::SIGMA(:,:),BC(:,:)
	REAL(8),ALLOCATABLE::LINE_LOAD(:),X_LOAD(:)



      CONTAINS
	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      SUBROUTINE CAL_N()
      IMPLICIT NONE
      N = 0.0D0
	IF (MTYPE .EQ. 3) THEN
	  N(1) = 1.0-KSI-YITA
        N(2) = KSI
	  N(3) = YITA
	ELSE IF (MTYPE .EQ. 4) THEN
	  N(1) = (1.0-KSI)*(1.0-YITA)/4.0
	  N(2) = (1.0+KSI)*(1.0-YITA)/4.0
	  N(3) = (1.0+KSI)*(1.0+YITA)/4.0
	  N(4) = (1.0-KSI)*(1.0+YITA)/4.0
	ELSE IF (MTYPE .EQ. 6) THEN
        N(1) = 2.0D0*KSI*KSI-KSI
	  N(2) = 2.0D0*YITA**2-YITA
	  N(3) = (1.0D0-2.0D0*KSI-2.0D0*YITA)*(1.0D0-KSI-YITA)
	  N(4) = 4.0D0*KSI*YITA
	  N(5) = 4.0D0*YITA*(1.0D0-YITA-KSI)
	  N(6) = 4.0D0*KSI*(1.0D0-YITA-KSI)
	ELSE IF (MTYPE .EQ. 8) THEN
        N(1) = -(1.0D0-KSI)*(1.0D0-YITA)*(KSI+YITA+1.0D0)/4.0D0
        N(2) = (1.0D0+KSI)*(1.0D0-YITA)*(KSI-YITA-1.0D0)/4.0D0
        N(3) = (1.0D0+KSI)*(1.0D0+YITA)*(KSI+YITA-1.0D0)/4.0D0
        N(4) = (1.0D0-KSI)*(1.0D0+YITA)*(-KSI+YITA-1.0D0)/4.0D0
	  N(5) = (1.0D0-KSI*KSI)*(1.0D0-YITA)/2.0D0
	  N(6) = (1.0D0-YITA**2)*(1.0D0+KSI)/2.0D0
	  N(7) = (1.0D0-KSI*KSI)*(1.0D0+YITA)/2.0D0
	  N(8) = (1.0D0-YITA**2)*(1.0D0-KSI)/2.0D0
      ELSE
	  PRINT*, 'WRONG ELEMENT TYPE IN SUBROUTINE CAL_N'
	END IF
      END SUBROUTINE CAL_N
	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      SUBROUTINE CAL_DN()
      IMPLICIT NONE
      DN = 0.0D0
	IF (MTYPE .EQ. 3) THEN
        DN(1,1:3) = (/1.0D0,0.0D0,-1.0D0/)
	  DN(2,1:3) = (/0.0D0,1.0D0,-1.0D0/)
	ELSE IF (MTYPE .EQ. 4) THEN
	  DN(1,1:4)=(/YITA-1.0,1.0-YITA,1.0+YITA,-1.0-YITA/)
	  DN(2,1:4)=(/KSI-1.0,-1.0-KSI,1.0+KSI,1.0-KSI/)
	  DN = DN/4.0
	ELSE IF (MTYPE .EQ. 6) THEN
        DN(1,1:6)=(/4.0D0*KSI-1.0D0,0.0D0,4.0D0*KSI+4.0D0*YITA-3.0D0,
     5         4.0D0*YITA,-4.0D0*YITA,4.0D0*(1.0D0-YITA-2.0D0*KSI)/)
	  DN(2,1:6)=(/0.0D0,4.0D0*YITA-1.0D0,4.0D0*KSI+4.0D0*YITA-3.0D0,
     5         4.0D0*KSI,4.0D0*(1.0D0-KSI-2.0D0*YITA),-4.0D0*KSI/)
	ELSE IF (MTYPE .EQ. 8) THEN
        DN(1,1:8)=(/(1.0D0-YITA)*(2.0D0*KSI+YITA)/4.0D0,
     5	          (1.0D0-YITA)*(2.0D0*KSI-YITA)/4.0D0,
     5	          (1.0D0+YITA)*(2.0D0*KSI+YITA)/4.0D0,
     5	          (1.0D0+YITA)*(2.0D0*KSI-YITA)/4.0D0,
     5	          -(1.0D0-YITA)*KSI,(1.0D0-YITA**2)/2.0D0,
     5	          -(1.0D0+YITA)*KSI,-(1.0D0-YITA**2)/2.0D0/)
        DN(2,1:8)=(/(1.0D0-KSI)*(2.0D0*YITA+KSI)/4.0D0,
     5	          (1.0D0+KSI)*(2.0D0*YITA-KSI)/4.0D0,
     5	          (1.0D0+KSI)*(2.0D0*YITA+KSI)/4.0D0,
     5	          (1.0D0-KSI)*(2.0D0*YITA-KSI)/4.0D0,
     5	          -(1.0D0-KSI**2)/2.0D0,-(1.0D0+KSI)*YITA,
     5	          (1.0D0-KSI**2)/2.0D0,-(1.0D0-KSI)*YITA/)
	ELSE 
	  PRINT*,'WRONG ELEMENT TYPE IN SUBROUTINE CAL_DN'
	END IF
      END SUBROUTINE CAL_DN
	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	SUBROUTINE CAL_DET()
      IMPLICIT NONE
      DET = JACOBI(1,1)*JACOBI(2,2)-JACOBI(1,2)*JACOBI(2,1)
	END SUBROUTINE CAL_DET
	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	SUBROUTINE CAL_INVJ()
      IMPLICIT NONE
	REAL(8)::TEMP1

      TEMP1 = JACOBI(1,1)
      JACOBI(1,1) = JACOBI(2,2)
	JACOBI(2,2) = TEMP1
      JACOBI(2,1) = -JACOBI(2,1)
	JACOBI(1,2) = -JACOBI(1,2)
      JACOBI=JACOBI/DET

	END SUBROUTINE CAL_INVJ
	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      SUBROUTINE COMPUT()
	IMPLICIT NONE

      DN6(1,1:6)=(/4.0D0*KSI-1.0D0,0.0D0,4.0D0*KSI+4.0D0*YITA-3.0D0,
     5         4.0D0*YITA,-4.0D0*YITA,4.0D0*(1.0D0-YITA-2.0D0*KSI)/)
	DN6(2,1:6)=(/0.0D0,4.0D0*YITA-1.0D0,4.0D0*KSI+4.0D0*YITA-3.0D0,
     5         4.0D0*KSI,4.0D0*(1.0D0-KSI-2.0D0*YITA),-4.0D0*KSI/)
      N8(1) = -(1.0D0-KSI)*(1.0D0-YITA)*(KSI+YITA+1.0D0)/4.0D0
      N8(2) = (1.0D0+KSI)*(1.0D0-YITA)*(KSI-YITA-1.0D0)/4.0D0
      N8(3) = (1.0D0+KSI)*(1.0D0+YITA)*(KSI+YITA-1.0D0)/4.0D0
      N8(4) = (1.0D0-KSI)*(1.0D0+YITA)*(-KSI+YITA-1.0D0)/4.0D0
	N8(5) = (1.0D0-KSI*KSI)*(1.0D0-YITA)/2.0D0
	N8(6) = (1.0D0-YITA**2)*(1.0D0+KSI)/2.0D0
	N8(7) = (1.0D0-KSI*KSI)*(1.0D0+YITA)/2.0D0
	N8(8) = (1.0D0-YITA**2)*(1.0D0-KSI)/2.0D0
      DN8(1,1:8)=(/(1.0D0-YITA)*(2.0D0*KSI+YITA)/4.0D0,
     5	          (1.0D0-YITA)*(2.0D0*KSI-YITA)/4.0D0,
     5	          (1.0D0+YITA)*(2.0D0*KSI+YITA)/4.0D0,
     5	          (1.0D0+YITA)*(2.0D0*KSI-YITA)/4.0D0,
     5	          -(1.0D0-YITA)*KSI,(1.0D0-YITA**2)/2.0D0,
     5	          -(1.0D0+YITA)*KSI,-(1.0D0-YITA**2)/2.0D0/)
      DN8(2,1:8)=(/(1.0D0-KSI)*(2.0D0*YITA+KSI)/4.0D0,
     5	          (1.0D0+KSI)*(2.0D0*YITA-KSI)/4.0D0,
     5	          (1.0D0+KSI)*(2.0D0*YITA+KSI)/4.0D0,
     5	          (1.0D0-KSI)*(2.0D0*YITA-KSI)/4.0D0,
     5	          -(1.0D0-KSI**2)/2.0D0,-(1.0D0+KSI)*YITA,
     5	          (1.0D0-KSI**2)/2.0D0,-(1.0D0-KSI)*YITA/)
      N6(1) = 2.0D0*KSI*KSI-KSI
	N6(2) = 2.0D0*YITA**2-YITA
	N6(3) = (1.0D0-2.0D0*KSI-2.0D0*YITA)*(1.0D0-KSI-YITA)
	N6(4) = 4.0D0*KSI*YITA
	N6(5) = 4.0D0*YITA*(1.0D0-YITA-KSI)
	N6(6) = 4.0D0*KSI*(1.0D0-YITA-KSI)
      
	END SUBROUTINE COMPUT

      END MODULE DEF
	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      SUBROUTINE MULTI(A,B,C,M,N,P)
	!MULTIPLY THE TWO MATRIX AND RESULT IN C(M,P)
	INTEGER(4)::M,N,P
	REAL(8)::A(M,N),B(N,P),C(M,P),TEMP(M,P)
	!ABOVE TO TRANSFER
	INTEGER(4)::I,J

	DO I = 1,M
	  DO J = 1,P
	    TEMP(I,J) = DOT_PRODUCT(A(I,:),B(:,J))
	  END DO 
	END DO

      C = TEMP

	END SUBROUTINE MULTI
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      SUBROUTINE READ_DATA()
      !READ DATA FROM DATA.TXT
	USE DEF
      IMPLICIT NONE
	!ABOVE TO TRANSFER
	INTEGER(4)::I,J,TEMP
	INTEGER(4)::G_NODENUM
	EXTERNAL::ELE_TRANS

      !G_NODE: FIRST COLUMN REPRESENT THE NUMBER OF NODE
	!        SECOND COLUMN REPRESENT THE X COORD OF NODE
	!        THIRD COLUMN REPRESENT THE Y COORD OF NODE
	!G_ELE: THE NODE NUMBER IN EACH ELEMENT
	!P_LOAD(1:N_LOAD): VALUE OF POINT LOAD
	!NP_LOAD(1:N_LOAD):POINT NUMBER OF POINT LOAD
	OPEN(UNIT=2,FILE='D:\MY PROGRAM IN FORTRAN\A\NODE_ANSYS.TXT')
	OPEN(UNIT=3,FILE='D:\MY PROGRAM IN FORTRAN\A\ELEMENT_ANSYS.TXT')
	OPEN(UNIT=5,FILE='D:\MY PROGRAM IN FORTRAN\A\BC_ANSYS.TXT')
	OPEN(UNIT=9,FILE='D:\MY PROGRAM IN FORTRAN\A\MARK.TXT')

      DO I = 1,N_NODE
        READ(2,*) G_NODENUM,G_NODE(I,1),G_NODE(I,2)
	END DO     
      CLOSE(2)

	DO I = 1,N_ELE
	  READ(9,*) MARK(I)
	  IF ((MARK(I) .NE. 3) .AND. (MARK(I) .NE. 6) 
     4	  .AND. (MARK(I) .NE. 4) .AND. (MARK(I) .NE. 8)) THEN
	    PRINT*,'ERROR ELEMENT TYPE IN SUROUTINE READ_DATA'
	  END IF
	END DO
	CLOSE(9)

      G_ELE = 0
      DO I = 1,N_ELE
	    READ(3,*) G_ELE(I,1:MARK(I))
	END DO
	CLOSE(3)

	!DO I = 1,N_ELE
	!  IF (MARK(I) .EQ. 4) THEN
	!    CALL ELE_TRANS(G_ELE(I,1:MARK(I)))         
	!END IF                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                 
      !END DO
	!其实自己转换会遇到很多问题,ansys输出已经很好了
      
	IF (N_LOAD .NE. 0) THEN	
	  OPEN(UNIT=4,FILE='D:\MY PROGRAM IN FORTRAN\A\LOAD_ANSYS.TXT')
	  DO I = 1,N_LOAD
          READ(4,*) NP_LOAD(I),P_LOAD(I,1),P_LOAD(I,2)
        END DO
	END IF
      CLOSE(4)

      IF (X_N_LOAD .NE. 0) THEN
	  OPEN(UNIT=6,FILE='D:\MY PROGRAM IN FORTRAN\A\X_LOAD.TXT')
	  DO I = 1,X_N_LOAD
	    READ(6,*) X_LOAD(I),X_LOAD_NODE(I,1),X_LOAD_NODE(I,2)
	  END DO
	END IF
      CLOSE(6)

      IF (N_LINE_LOAD .NE. 0) THEN
	  OPEN(UNIT=7,FILE='D:\MY PROGRAM IN FORTRAN\A\LINE_LOAD.TXT')
	  DO I = 1,N_LINE_LOAD
	    READ(7,*) LINE_LOAD_NODE(I,1),LINE_LOAD(I),LINE_LOAD_NODE(I,2)
	  END DO
	END IF
      CLOSE(7)

      DO I = 1,N_BC
        READ(5,*) N_DIS(I),BC(I,1),TEMP,BC(I,2)
	END DO
	CLOSE(5)

	END SUBROUTINE READ_DATA
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      SUBROUTINE ELE_TRANS(G_E)
      !TRANSFORM THE ORDER OF NODE IN ELEMENT
      USE DEF
      IMPLICIT NONE
	INTEGER(4)::G_E(1:4)
  	!ABOVE TO TRANSFER
	INTEGER(4)::I,J,NN(1:4),T(1:4)
	REAL(8)::X(1:4),Y(1:4)
	REAL(8)::R
      EXTERNAL::Bi_Directional_Bubble_Sort

      DO I = 1,4
	  X(I) = G_NODE(G_E(I),1)
	  Y(I) = G_NODE(G_E(I),2)
      END DO
      NN(1:4)=(/1,2,3,4/)

	CALL Bi_Directional_Bubble_Sort(X,NN,4)

      IF (Y(NN(1)) .LT. Y(NN(2))) THEN
        T(1) = G_E(NN(1))
	  T(4) = G_E(NN(2))
	  IF (Y(NN(3)) .LT. Y(NN(4))) THEN
	    T(2) = G_E(NN(3))
	    T(3) = G_E(NN(4))
	  ELSE
	    T(3) = G_E(NN(3))
	    T(2) = G_E(NN(4))
	  END IF
	ELSE
        T(4) = G_E(NN(1))
	  T(1) = G_E(NN(2))
	  IF (Y(NN(3)) .LT. Y(NN(4))) THEN
	    T(2) = G_E(NN(3))
	    T(3) = G_E(NN(4))
	  ELSE
	    T(3) = G_E(NN(3))
	    T(2) = G_E(NN(4))
	  END IF
	END IF
      G_E=T

	END SUBROUTINE ELE_TRANS
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      SUBROUTINE Bi_Directional_Bubble_Sort(X,NN,NUM)
      !GIVE A ORDER OF FROM LOWER TO TOP OF FOUR NUMBERS
      REAL(8)::X(1:4)
	INTEGER(4)::NN(1:4),NUM
	!ABOVE TO TRANSFER
      REAL(8)::R
	INTEGER(4)::I,Q,LOW,UP,T

      LOW = 1
	UP = NUM
      DO WHILE (UP .GT. LOW)
        T = LOW
	  DO I = LOW,UP-1
	    IF (X(I) .GT. X(I+1)) THEN
	      R = X(I+1)
	      X(I+1) =X(I)
	      X(I) = R
	      Q = NN(I)
	      NN(I) = NN(I+1)
	      NN(I+1) = Q
	      T = I
	    END IF
        END DO
	  UP = T
	  DO I = UP,LOW+1,-1
	    IF (X(I) .LT. X(I-1)) THEN
	      R = X(I-1)
	      X(I-1) =X(I)
	      X(I) = R
	      Q = NN(I)
	      NN(I) = NN(I-1)
	      NN(I-1) = Q
	      T = I
          END IF
	  END DO
	  LOW = T
      END DO

	END SUBROUTINE Bi_Directional_Bubble_Sort
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      SUBROUTINE CAL_PE()
      ! COMPUTE THE ELASTIC MODULE AND POISSON RATIO FOR DIFFERENT PROBLEM
      USE DEF
      IMPLICIT NONE

	IF (ID .EQ. 1) THEN
	  PE = PE
	  PR = PR
	ELSE IF (ID. EQ. 2) THEN
	  PE = PE/(1.0-PR*PR)
	  PR = PR/(1.0-PR)
	END IF

	END SUBROUTINE CAL_PE

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      SUBROUTINE FORM_KE()
	!FORM THE MATRIX OF ELEMENT K
	USE DEF
      IMPLICIT NONE
	!ABOVE TO TRANSFER
	INTEGER(4)::I,J,P,Q,TT,K
	REAL(8)::DD(3,3),BB(3,16)
	REAL(8)::TEMP1(3,16),TEMP2(16,16)
	EXTERNAL::MULTI,CAL_BB,CAL_DD,CAL_POINT,FIND_34

	KE = 0.0
      IF (MTYPE .EQ. 3) THEN
        CALL FIND_34()
	  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 FORM_KE'
	  END IF
        DO K = 1,TT
	    KSI = POINT3(K,1)
	    YITA = POINT3(K,2)
	    CALL CAL_DN()
	    CALL  COMPUT()
	    DO I = 1,MTYPE
            IF (DEVIDE(I) .EQ. 1) THEN
              DN(1,I) = DN6(1,I)
              DN(2,I) = DN6(2,I)
              DN(1,3+I) = DN6(1,3+I)
              DN(2,3+I) = DN6(2,3+I)
            END IF
          END DO
          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. 4) THEN
        CALL FIND_34()
	  CALL CAL_POINT()
	  DO I = 1,N_INT
	    DO J = 1,N_INT
	      YITA = POINT(I)
	      KSI = POINT(J)
	      CALL CAL_DN()
	      CALL  COMPUT()
	      DO K = 1,MTYPE
              IF (DEVIDE(K) .EQ. 1) THEN
                DN(1,K) = DN8(1,K)
                DN(2,K) = DN8(2,K)
                DN(1,4+K) = DN8(1,4+K)
                DN(2,4+K) = DN8(2,4+K)
              END IF
            END DO 
            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

	ELSE IF (MTYPE .EQ. 6) THEN
	  CALL CAL_POINT()
	  IF (N_INT .EQ. 1) THEN