eqsolve2.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
	INTEGER(4)::S34
	REAL(8)::COORD(4,2),DET,JACOBI(2,2),N(1:4),DN(2,4),KSI,YITA
!READ_DATA      
      INTEGER(4)::N_NODE,N_ELE,N_LOAD,N_BC,N_LINE_LOAD,X_N_LOAD
      INTEGER(4)::MARK(:)      

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

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

!FORM_K
      INTEGER(4)::DOF
!FORM_LOAD
	INTEGER(4)::N_LINEAR(1:2)
	REAL(8)::ROU,Q_LINEAR
!ADD_BC
!CAL_STR
      REAL(8)::VEX(2,4)
!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 (S34 .EQ. 3) THEN
	  N(1) = 1.0-KSI-YITA
        N(2) = KSI
	  N(3) = YITA
	ELSE
	  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
	END IF
      END SUBROUTINE CAL_N
	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      SUBROUTINE CAL_DN()
      IMPLICIT NONE
      DN = 0.0
	IF (S34 .EQ. 3) THEN
        DN(1,:) = (/-1.0D0,1.0D0,0.0D0,0.0D0/)
	  DN(2,:) = (/-1.0D0,0.0D0,1.0D0,0.0D0/)
	ELSE
	  DN(1,:)=(/YITA-1.0,1.0-YITA,1.0+YITA,-1.0-YITA/)
	  DN(2,:)=(/KSI-1.0,-1.0-KSI,1.0+KSI,1.0-KSI/)
	  DN = DN/4.0
	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 CAL_VEX()
      IMPLICIT NONE
      VEX(:,1)=(/-1.0D0,-1.0D0/)
	VEX(:,2)=(/1.0D0,-1.0D0/)
	VEX(:,3)=(/1.0D0,1.0D0/)
	VEX(:,4)=(/-1.0D0,1.0D0/)
      END SUBROUTINE CAL_VEX

      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
	INTEGER(4)::NE_ANSYS(1:N_ELE,14),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')

      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(3,*) NE_ANSYS(I,1:14)
	END DO
	CLOSE(3)

      DO I = 1,N_ELE
	  DO J = 1,4
	    G_ELE(I,J) = NE_ANSYS(I,J)
	  END DO
	END DO

	DO I = 1,N_ELE
	  IF (G_ELE(I,3) .NE. G_ELE(I,4)) THEN
	    CALL ELE_TRANS(G_ELE(I,:))         
	  END IF 
	END DO
      
	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),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
	REAL(8)::DD(3,3),BB(3,8),TEMP1(3,6),TEMP2(3,6),TEMP3(6,6)
	REAL(8)::TEMP4(3,8),TEMP5(8,8)
	EXTERNAL::MULTI,CAL_BB,CAL_DD,CAL_POINT

	KE = 0.0
      IF (S34 .EQ. 3) THEN
	  KSI = 0.0
	  YITA = 0.0
        CALL CAL_BB(BB)
	  DO I = 1,3
	    TEMP1(I,:) = BB(I,1:6)
	  END DO
	  CALL CAL_DD(DD)
	  CALL MULTI(DD,TEMP1,TEMP2,3,3,6)
	  CALL MULTI(TRANSPOSE(TEMP1),TEMP2,TEMP3,6,3,6)
        TEMP3 = PT*DET*TEMP3/2.0      !REMEMBER THE AREA
	  DO I = 1,6
	    KE(I,1:6) = TEMP3(I,:)
	  END DO
	ELSE IF (S34 .EQ. 4) THEN
	  CALL CAL_POINT()
	  DO I = 1,N_INT
	    DO J = 1,N_INT
	      YITA = POINT(I)
	      KSI = POINT(J)
            CALL CAL_BB(BB)
	      CALL CAL_DD(DD)
	      CALL MULTI(DD,BB,TEMP4,3,3,8)
	      CALL MULTI(TRANSPOSE(BB),TEMP4,TEMP5,8,3,8)
	      DO P = 1,S34*2
	        DO Q = 1,S34*2
                KE(P,Q) = KE(P,Q) + TEMP5(P,Q)*H(J)*H(I)
	        END DO
	      END DO
	    END DO
	  END DO
        KE = PT*DET*KE
	END IF

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

      BB = 0.0D0
	CALL CAL_DN()
	CALL MULTI(DN,COORD,JACOBI,2,4,2)
	CALL CAL_DET()
	CALL CAL_INVJ()
	DO I = 1,S34
        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 (S34 .EQ. 3) 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 (S34 .EQ. 4) 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(4,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 
        CALL FORM_ORDER(T_NODE,I)

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

	  IF (G_ELE(I,3) .EQ. G_ELE(I,4)) THEN
	    T = 6
	  ELSE 
	    T = 8
	  END IF
        DO J = 1,T/2        
	    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
          CALL FORM_ORDER(T_NODE,J)
	    IF (G_ELE(J,3) .EQ. G_ELE(J,4)) THEN
	      T=3
	    ELSE
	      T=4
	    END IF
	    DO K = 1,T
	      IF ((T_NODE(K,1) .EQ. I) .OR. (T_NODE(K,2) .EQ. I))THEN
	        DO M = 1,T
                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