PROGRAM SET3
C=======================================================================
C      DATA GENERATING PROGRAM FOR DOMAIN.FOR
C   PROJECT: DOUBLE-LET   PROJECT NAME: DOMAIN
C    DOMAIN: INFINITE
C    ELEMENT: FOUR-NODED ISOPARAMETRIC ELEMENT
C    DOMAIN DISCRETIZATION: UNEVEN ELEMENTS WITH VERTICAL SCAN
C     EIJI FUKUMORI DECEMBER 29, 1993 
C=======================================================================
      PARAMETER ( ND=4,MXE=500,MXN=520, MXI=100 )
      PARAMETER (MXD=10)
C=======================================================================
      DIMENSION NODEX(MXE,ND),XCOORD(MXN), YCOORD(MXN),FJ(MXE),ICIR(16)
      DIMENSION XI(MXI), YI(MXI)
C=======================================================================
      DATA NC, ICIR /16, 20,25,24,23,22,21,16,11,6,1,2,3,4,5,10,15/
C=======================================================================
      PI = 4.* ATAN( 1.)
      RADIUS = 1.
      XC1 = 0.
      YC1 = 2.
      WRITE(*,*) 'TYPE IN HALF OF D.  I.E. (D/2)'
      READ (*,*) YC1
      IF ( YC1 .LT. 2*RADIUS ) STOP 'YC1>2.'
      XC2 = XC1
      YC2 = -YC1
C=======================================================================
C-----     RADIUS OF WELL = 1. I.E., A=1
      DO I = 1 , NC
      XCOORD(ICIR(I)) = COS ( I*PI/8. )
      YCOORD(ICIR(I)) = SIN ( I*PI/8. )
      END DO
      XCOORD(7) = -0.5
      XCOORD(8) =  0.
      XCOORD(9) =  0.5
      YCOORD(7) = -0.5
      YCOORD(8) = -0.5
      YCOORD(9) = -0.5
C
      XCOORD(12) = -0.5
      XCOORD(13) =  0.
      XCOORD(14) =  0.5
      YCOORD(12) =  0.
      YCOORD(13) =  0.
      YCOORD(14) =  0.
C
      XCOORD(17) = -0.5
      XCOORD(18) =  0.
      XCOORD(19) =  0.5
      YCOORD(17) =  0.5
      YCOORD(18) =  0.5
      YCOORD(19) =  0.5
      DO I = 1 , 25
      XCOORD(I+25) = XCOORD(I) + XC2
      YCOORD(I+25) = YCOORD(I) + YC2
      XCOORD(I) = XCOORD(I) + XC1
      YCOORD(I) = YCOORD(I) + YC1
      END DO
      NNODE = 25*2
C=======================================================================
C                    ELEMENT CREATION
      NEX = 4
      NEY = 4
      NDX = NEX + 1
      NDY = NEY + 1
      NE = 0
      DO I = 1 , NEY
      DO J = 1 , NEX
      NE = NE + 1
      IF ( NE .GT. MXE ) STOP 'NE > MXE'
      NODEX(NE,1) = NDX*(I-1) + J
      NODEX(NE,2) = NODEX(NE,1) + 1
      NODEX(NE,3) = NODEX(NE,2) +NDX
      NODEX(NE,4) = NODEX(NE,1)+NDX
      FJ(NE) = 1./PI
      END DO
      END DO
C
      DO I = 1 , NEY
      DO J = 1 , NEX
      NE = NE + 1
      IF ( NE .GT. MXE ) STOP 'NE > MXE'
      NODEX(NE,1) = NDX*(I-1) + J + NDX*NDY
      NODEX(NE,2) = NODEX(NE,1) + 1
      NODEX(NE,3) = NODEX(NE,2) +NDX
      NODEX(NE,4) = NODEX(NE,1)+NDX
      FJ(NE) = -1./PI
      END DO
      END DO
C=======================================================================
      NIP = 0
      NIP = NIP + 1
      XI(NIP) = 0.
      YI(NIP) = 0.0001
      NIP = NIP + 1
      XI(NIP) = 0.
      YI(NIP) = (YC1-RADIUS)/4.
      NIP = NIP + 1
      XI(NIP) = 0.
      YI(NIP) = (YC1-RADIUS)/2.
      NIP = NIP + 1
      XI(NIP) = 0.
      YI(NIP) = (YC1-RADIUS)/4.*3.
      NIP = NIP + 1
      XI(NIP) = 0.
      YI(NIP) = YC1 - RADIUS/4.*3.
      NIP = NIP + 1
      XI(NIP) = 0.
      YI(NIP) = YC1 - RADIUS/4.
      NIP = NIP + 1
      XI(NIP) = 0.
      YI(NIP) = YC1 + RADIUS/4.
      NIP = NIP + 1
      XI(NIP) = 0.
      YI(NIP) = YC1 + RADIUS/4.*3.
      NIP = NIP + 1
      XI(NIP) = 0.
      YI(NIP) = YC1 + RADIUS*1. + RADIUS/4.
      NIP = NIP + 1
      XI(NIP) = 0.
      YI(NIP) = YC1 + RADIUS*2.
      NIP = NIP + 1
      XI(NIP) = 0.
      YI(NIP) = YC1 + RADIUS*3.
      NIP = NIP + 1
      XI(NIP) = 0.
      YI(NIP) = YC1 + RADIUS*4.
      NIP = NIP + 1
      XI(NIP) = 0.
      YI(NIP) = YC1 + RADIUS*5.
C
      DO I = 3 , 30
      NIP = NIP + 1
      XI(NIP) = 0.
      YI(NIP) = YC1 + RADIUS*I*2
      END DO
C
      NNN = NIP
      DO I = 1 , NNN
      NIP = NIP + 1
      XI(NIP) = 0.
      YI(NIP) = -YI(I)
      END DO
C
C
      DO I = 1 , NIP-1
      YMIN = YI(I)
      K = I
      DO J = I+1 , NIP
      IF ( YI(J) .LT. YMIN ) THEN
      YMIN = YI(J)
      K = J
      END IF
      END DO
      YI(K) = YI(I)
      YI(I) = YMIN
      END DO
C=======================================================================
C                     MAKING DATA FILES
      OPEN ( 1, FILE='DOMAIN.DAT', STATUS = 'UNKNOWN' )
      WRITE(1,100) NE
C------ ELEMENTS CONECTIVITY INFORMATION
      DO I = 1 , NE
      WRITE (1,100) I,(NODEX(I,J),J=1,ND), FJ(I)
      END DO
C------ NODAL COORDINATES
      WRITE(1,100) NNODE
      DO I = 1 , NNODE
      WRITE(1,200) I, XCOORD(I), YCOORD(I)
      END DO
C------ INTERNAL POINT
      WRITE(1,*) NIP
      IF ( NIP .GE. 1 ) THEN
      DO I = 1 , NIP
      WRITE(1,*) XI(I), YI(I)
      END DO
      END IF
      CLOSE (1)
  100 FORMAT ( 5I7 ,  E20.10      )
  200 FORMAT ( I7  , 2E20.10 )
  300 FORMAT ( 2I7 , 2E16.5 )
  400 FORMAT ( I7  , 2E16.5 )
      WRITE(*,*) 'NE=',NE
      WRITE(*,*) 'NNODE=',NNODE
      STOP
      END