PROGRAM SETNEUM2
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: rectangular
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=======================================================================
C--------------OUTER BOUNDARY IS DEVIDED INTO 8 SEGMENTS
C ELEMENTS IN EACH SEGMENT "NEOUTER".
C HENCE, ELEMENTS ON OUTER BOUNADRY = NEOUTER * 8
NEOUTER = 4
NEINNER = 4*4
C=======================================================================
PI = 4.* DATAN( 1.D0)
RADIUSA = 6.35D0
RADIUSB = 9.5D0
RADIUSC = 12.0D0
IF ( RADIUSB .LE. RADIUSA ) STOP 'RADIUSB MUST BE > RADIUSA.'
IF ( RADIUSC .LT. RADIUSB ) STOP 'RADIUSC MUST BE >= RADIUSB.'
WRITE (*,*) 'RADIUS-A(INNER) = ', RADIUSA
WRITE (*,*) 'RADIUS-B(OUTER-VERTICA L) = ', RADIUSB
WRITE (*,*) 'RADIUS-C(OUTER-HORIZONTAL) = ', RADIUSC
C=======================================================================
C----------------- BOUNDARY CONDITIONS
C------ ON OUTER BOUNDARY: DIRICHLET ----> POTENTIAL(POUTER) = 0
C------ ON INNER BOUNDARY: NEUMANN ------> POTENTIAL DERIVATIVE = BINNER
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(*,*) 'NEUMANN BOUNDARY VALUE =', BINNER
C=======================================================================
C------------------- NORDAL COORDINATES CREATION
C------------------------- OUTER ELEMENTS -------------
C-------------- FIRST SEGMENT STARTS FROM 0 TO PI/4 ---------------
DX = RADIUSC/(2*NEOUTER)
DY = RADIUSB/(2*NEOUTER)
NODE = 0
DO I = 1 , 2*NEOUTER
NODE = NODE + 1
XCOORD (NODE) = RADIUSC
YCOORD (NODE) = DY * (I-1)
END DO
C-------------- SECOND SEGMENT STARTS FROM PI/4 TO PI/2 ---------------
DO I = 1 , 2*NEOUTER
NODE = NODE + 1
XCOORD (NODE) = RADIUSC - DX*(I-1)
YCOORD (NODE) = RADIUSB
END DO
C-------------- THIRD SEGMENT STARTS FROM PI/2 TO (3/4)PI ---------------
DO I = 1 , 2*NEOUTER
NODE = NODE + 1
XCOORD (NODE) = -DX*(I-1)
YCOORD (NODE) = RADIUSB
END DO
C-------------- 4TH SEGMENT STARTS FROM (3/4)PI TO PI ---------------
DO I = 1 , 2*NEOUTER
NODE = NODE + 1
XCOORD (NODE) = -RADIUSC
YCOORD (NODE) = RADIUSB - DY*(I-1)
END DO
C-------------- 5TH SEGMENT STARTS FROM PI TO (5/4)PI ---------------
DO I = 1 , 2*NEOUTER
NODE = NODE + 1
XCOORD (NODE) = -RADIUSC
YCOORD (NODE) = -DY*(I-1)
END DO
C-------------- 6TH SEGMENT STARTS FROM (5/4)PI TO (6/4)PI ---------------
DO I = 1 , 2*NEOUTER
NODE = NODE + 1
XCOORD (NODE) = -RADIUSC + DX*(I-1)
YCOORD (NODE) = -RADIUSB
END DO
C-------------- 7TH SEGMENT STARTS FROM (6/4)PI TO (7/4)PI ---------------
DO I = 1 , 2*NEOUTER
NODE = NODE + 1
XCOORD (NODE) = DX*(I-1)
YCOORD (NODE) = -RADIUSB
END DO
C-------------- 8TH SEGMENT STARTS FROM (6/4)PI TO (7/4)PI ---------------
DO I = 1 , 2*NEOUTER
NODE = NODE + 1
XCOORD (NODE) = RADIUSC
YCOORD (NODE) = -RADIUSB + DY*(I-1)
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
NODE = NODE + 1
XCOORD (NODE) = DREAL( BASEVEC )
YCOORD (NODE) = DIMAG( BASEVEC )
BASEVEC = BASEVEC * UNITVEC
END DO
C=======================================================================
C ELEMENT CREATION
C----------- OUTER ELENENTS -----------
C----------- DIRICHLET BOUNDARY CONDITION ---------
DO I = 1 , 8*NEOUTER
II = I*2-1
JJ = II + 1
KK = JJ + 1
IF ( I .EQ. 8*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+8*NEOUTER , NEINNER+8*NEOUTER
II = I*2-1
JJ = II + 1
KK = JJ + 1
IF ( I .EQ. NEINNER+8*NEOUTER ) KK = (1+8*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+8*NEOUTER
NNODE = NODE
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