PROGRAM SETTRS9C
C=======================================================================
C     DATA GENERATING PROGRAM FOR torsion9q.for
C                 CIRCULER DOMAIN         1/4 MODEL
C     DOMAIN SIZE: TLX BY TLY
C                         BOUNDARY CONDITIONS 
C        DIRICHLET: STRESS FUNCTION(S) ALONG THE OUTER BOUNDARY = 0
C          NEUMANN: DS/DN=0 ALONG THE SYMMETRY LINES
C     EIJI FUKUMORI NOVEMBER 28, 1995, NOV. 03, 2024
C                                  NUMBERING OF NODEAL POINTS
C  FOR CASE OF NEX=2 & NEY=2     21-----22-----23-----24-----25
C                                 |             |             |
C                                16     17     18     19     20
C                                 |             |             |
C                                11-----12-----13-----14-----15
C                                 |             |             |
C                                 6      7      8      9     10
C                                 |             |             |
C                                 1------2------3------4------5
C NOEDX OF 1ST ELEMENT = 1, 2, 3, 8, 13, 12, 11, 6, 7
C=======================================================================
      IMPLICIT REAL*8 ( A-H , O-Z )
      PARAMETER ( ND=9, MXE=31300, MXN=34410, MXB=12100 )
CCCCCCCCCCCCCCCCCCCCCCCC      PARAMETER (   TLX = 1.,   TLY = 1. )
C=======================================================================
      DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),IBNDT(MXB),
     *          BVT(MXB), ITYPET(MXB)
C=======================================================================
C         NEY: NEMBER OF VERTICAL ELEMENTS (NUMBER OF NODES: NEY+1)
C         NEX: NEMBER OF HORIZONTAL ELEMENTS (NUMBER OF NODES: NEX+1)
C=======================================================================
      TLX = 1.094693879D0
      TLY = 1.094693879D0
      NEX = 10
      NEY = 10
C
      DX = TLX / NEX
      DY = TLY / NEY
      GM = 26538461.5384615D0
      THETA = 0.0000670055862315403D0
C=======================================================================
C                    ELEMENT CREATION
      NE = 0
      DO I = 1 , NEY
      DO J = 1 , NEX
      NE = NE + 1
      IF ( NE .GT. MXE ) STOP 'NE > MXE'
      NODEX(NE,1) = (I-1)*(NEX*4+2) + (j-1)*2 + 1
      NODEX(NE,2) = NODEX(NE,1) + 1
      NODEX(NE,3) = NODEX(NE,2) + 1
      NODEX(NE,4) = NODEX(NE,3) + NEX*2 + 1
      NODEX(NE,5) = NODEX(NE,4) + NEX*2 + 1
      NODEX(NE,6) = NODEX(NE,5) - 1
      NODEX(NE,7) = NODEX(NE,6) - 1
      NODEX(NE,8) = NODEX(NE,7) - (NEX*2+1)
      NODEX(NE,9) = NODEX(NE,8) + 1
      END DO
      END DO
C=======================================================================
C                 NODAL COORDINATE CREATION
      NNODE = 0
      DO I = 1 , NEY*2+1
      DO J = 1 , NEX*2+1
      NNODE = NNODE + 1
      IF ( NNODE .GT. MXN ) STOP 'NNODE > MXN'
      XCOORD(NNODE) =  DX/2.D0*(J-1)
      YCOORD(NNODE) =  DY/2.D0*(I-1)
      END DO
      END DO
C-------------------------- Circularization ----------------------------
      DO I = 1 , NNODE
      IF ( (XCOORD(I) .NE. 0.D0) ) THEN
      ANGLE = DATAN ( YCOORD(I)/XCOORD(I) )
      R = DMAX1 (XCOORD(I), YCOORD(I))
      XCOORD(I) = R*DCOS(ANGLE)
      YCOORD(I) = R*DSIN(ANGLE)
      END IF
      END DO
C=======================================================================
      PI = 4.D0* DATAN (1.D0)
C                    BOUNDARY CONDITIONS
C--------- MOMENTUM EQUATIONS AND HEAT EQUATION
      NBT  = 0
C--------- FACE OF +Y AND -Y
C--------- -Y
C--------- +Y
      DO J = 1 , NEX*2+1
      NBT = NBT + 1
      IBNDT(NBT) = NEY*(NEX*4 + 2) + J
      ITYPET(NBT) = 1
      BVT(NBT) = 0.D0
      END DO
C--------- FACE OF +X AND -X
C           +X
      DO I = 1 , NEY*2
      NBT  = NBT  + 1
      IBNDT (NBT) = I*(2*NEX+1)
      ITYPET(NBT) = 1
      BVT(NBT ) = 0.D0
      END DO
C=======================================================================
      WRITE (*,*) ' NUMBER OF ELEMENTS (NE) = ',NE
      WRITE (*,*) ' NUMBER OF NODAL POINTS (NNODE) = ',NNODE
C=======================================================================
      OPEN ( 1, FILE='ELEMENT9.OUT', STATUS = 'UNKNOWN' )
      DO IEL = 1 , NE
      WRITE(1,*) XCOORD(NODEX(IEL,1)),YCOORD(NODEX(IEL,1))
      WRITE(1,*) XCOORD(NODEX(IEL,2)),YCOORD(NODEX(IEL,2))
      WRITE(1,*) XCOORD(NODEX(IEL,3)),YCOORD(NODEX(IEL,3))
      WRITE(1,*)
      WRITE(1,*) XCOORD(NODEX(IEL,3)),YCOORD(NODEX(IEL,3))
      WRITE(1,*) XCOORD(NODEX(IEL,4)),YCOORD(NODEX(IEL,4))
      WRITE(1,*) XCOORD(NODEX(IEL,5)),YCOORD(NODEX(IEL,5))
      WRITE(1,*)
      WRITE(1,*) XCOORD(NODEX(IEL,5)),YCOORD(NODEX(IEL,5))
      WRITE(1,*) XCOORD(NODEX(IEL,6)),YCOORD(NODEX(IEL,6))
      WRITE(1,*) XCOORD(NODEX(IEL,7)),YCOORD(NODEX(IEL,7))
      WRITE(1,*)
      WRITE(1,*) XCOORD(NODEX(IEL,7)),YCOORD(NODEX(IEL,7))
      WRITE(1,*) XCOORD(NODEX(IEL,8)),YCOORD(NODEX(IEL,8))
      WRITE(1,*) XCOORD(NODEX(IEL,1)),YCOORD(NODEX(IEL,1))
      WRITE(1,*)
      WRITE(1,*) XCOORD(NODEX(IEL,9)),YCOORD(NODEX(IEL,9))
      WRITE(1,*)
      END DO
      CLOSE (1)
C=======================================================================
C                     MAKING DATA FILES
      IW = 1
      OPEN ( IW, FILE='TRS9DATA.DAT', STATUS='UNKNOWN' )
      WRITE (IW,*) GM, THETA
      WRITE (IW,*) NE
      DO IEL = 1 , NE
        WRITE (IW,*) IEL,(NODEX(IEL,J),J=1,ND)
      END DO
      WRITE (IW,*) NNODE
      DO NODE = 1 , NNODE
        WRITE (IW,*) NODE,XCOORD(NODE),YCOORD(NODE)
      END DO
      WRITE (IW,*) NBT
      DO I = 1 , NBT
        WRITE (IW,*) IBNDT(I), ITYPET(I), BVT(I)
      END DO
      CLOSE (1)
      STOP
      END