PROGRAM SETNEUM
C=======================================================================
C              DATA GENERATING PROGRAM FOR BEM8QUDQ.FOR
C                   PROJECT:COAXIAL CABLE ANALYSIS
C  PROJECT NAME: COAXIAL NEUMANN BOUNDARY ON INNER DIRICHLET ON OUTER
C             INNER SHAPE: CIRCLE   OUTER SHAPE: CIRCLE
C                            DOMAIN: FINITE
C                 ELEMENT: 3-NODED PARAMETRIC ELEMENT
C                     EIJI FUKUMORI FEB 08, 2022
C=======================================================================
      IMPLICIT REAL*8 ( A-H , O-Z )
      PARAMETER ( ND=3,MXE=500,MXN=2*MXE, MXI=100 )
C=======================================================================
      DIMENSION NODEX(MXE,ND),XCOORD(MXN), YCOORD(MXN)
      DIMENSION XI(MXI), YI(MXI), IELTYPE(MXE),BV(MXE,ND)
      COMPLEX*16 BASEVEC , UNITVEC 
C=======================================================================
      NEOUTER = 4*2
      NEINNER = 4*1
C=======================================================================
      PI = 4.* DATAN( 1.D0)
      RADIUSA = 0.7D0
      RADIUSB = 2.4D0
      IF ( RADIUSB .LE. RADIUSA ) STOP 'RADIUSB MUST BE > RADIUSA.'
      WRITE (*,*) 'RADIUS-A(INNER) = ', RADIUSA
      WRITE (*,*) 'RADIUS-B(OUTER) = ', RADIUSB
C=======================================================================
C----------------- BOUNDARY CONDITIONS
C---------- ON OUTER BOUNDARY: DIRICHLET ----> POTENTIAL = 0
C---------- ON INNER BOUNDARY: NEUMANN ------> POTENTIAL DERIVATIVE = B
C---------- B INTO DOMAIN     -----> THE VALUE IS NEGATIVE
C---------- B OUT FROM DOMAIN -----> THE VALUE IS POSITIVE
       Q = 1.D0
       POUTER =  0.D0
       BINNER = -Q/(2.D0*PI*RADIUSA)
      WRITE (*,*) 'CHARGE = ', Q
      WRITE (*,*) 'DIRICHLET: POTENTIAL ON OUTER = ', POUTER
      WRITE (*,*) 'NEUMANN: DPDN ON INNER = ', BINNER
C=======================================================================
C-----     NORDAL COORDINATES CREATION
C------------ OUTER ELEMENTS -------------
      BASEVEC = DCMPLX( RADIUSB, 0.D0 )
      UNITANG = 2.D0*PI / (2*NEOUTER)
      UNITVEC =  DCMPLX( DCOS(UNITANG), DSIN(UNITANG) )
      DO I = 1 , 2*NEOUTER
      XCOORD (I) =  DREAL( BASEVEC )
      YCOORD (I) =  DIMAG( BASEVEC )
      BASEVEC = BASEVEC * UNITVEC 
      END DO
C------------ INNER ELEMENTS -------------
      BASEVEC = DCMPLX( RADIUSA, 0.D0 )
      UNITANG = -2.D0*PI / (2*NEINNER)
      UNITVEC =  DCMPLX( DCOS(UNITANG), DSIN(UNITANG) )
      DO I = 1 , 2*NEINNER
      J = 2*NEOUTER + I
      XCOORD (J) =  DREAL( BASEVEC )
      YCOORD (J) =  DIMAG( BASEVEC )
      BASEVEC = BASEVEC * UNITVEC 
      END DO
C=======================================================================
C                    ELEMENT CREATION
C----------- OUTER ELENENTS -----------
C----------- DIRICHLET BOUNDARY CONDITION ---------
      DO I = 1 , NEOUTER
      II = I*2-1
      JJ = II + 1
      KK = JJ + 1
      IF ( I .EQ. NEOUTER ) KK = 1
      NODEX(I,1) = II
      NODEX(I,2) = JJ
      NODEX(I,3) = KK
      BV(I,1) = POUTER
      BV(I,2) = POUTER
      BV(I,3) = POUTER
      IELTYPE(I) = 1
      END DO
C---------- INNER ELEMENTS ---------------
C---------- NEUMANN BOUNDARY CONDITION ------------
      DO I = 1+NEOUTER , NEINNER+NEOUTER
      II = I*2-1
      JJ = II + 1
      KK = JJ + 1
      IF ( I .EQ. NEINNER+NEOUTER ) KK = (1+NEOUTER)*2-1
      NODEX(I,1) = II
      NODEX(I,2) = JJ
      NODEX(I,3) = KK
      BV(I,1) = BINNER
      BV(I,2) = BINNER
      BV(I,3) = BINNER
      IELTYPE(I) = 2
      END DO
C=======================================================================
      NE = NEINNER+NEOUTER
      NNODE = 2*(NEINNER+NEOUTER)
C=======================================================================
C--------INTERNAL POINTS
      NIP = 20
      DR = (RADIUSB - RADIUSA) / (NIP+1)
      DO I = 1 , NIP
      YI(I) = 0.D0
      XI(I) = RADIUSA + DR*I
      END DO
C=======================================================================
C                     MAKING DATA FILES
      OPEN ( 1, FILE='BEM2.DAT', STATUS='UNKNOWN' )
      WRITE (1,*) NE
      DO I = 1 , NE
      WRITE (1,*) I,(NODEX(I,J),J=1,ND),IELTYPE(I), (BV(I,J),J=1,ND)
      END DO
      WRITE (1,*) NNODE
      DO I = 1 , NNODE
      WRITE (1,*) I, XCOORD(I), YCOORD(I)
      END DO
C------ INTERNAL POINT
      WRITE(1,*) NIP
      IF ( NIP .GE. 1 ) THEN
      DO I = 1 , NIP
      WRITE(1,*) I, XI(I), YI(I)
      END DO
      END IF
      CLOSE (1)
      WRITE(*,*) 'NE=',NE
      WRITE(*,*) 'NNODE=',NNODE
C------------- BOUNDARY configuration
      OPEN ( 2, FILE='BOUNDARY.DAT', STATUS='UNKNOWN' )
      DO IEL = 1 , NE
      WRITE (2,*) XCOORD(NODEX(IEL,1)), YCOORD(NODEX(IEL,1))
      WRITE (2,*) XCOORD(NODEX(IEL,2)), YCOORD(NODEX(IEL,2))
      WRITE (2,*) XCOORD(NODEX(IEL,3)), YCOORD(NODEX(IEL,3))
      WRITE (2,*)
      END DO
      CLOSE (2)
      STOP
      END