PROGRAM SETNEUM0
C=======================================================================
C ********************** CONSTANT ELEMENT ************************
C DATA GENERATING PROGRAM FOR BEM8QUDQ.FOR
C PROJECT:COAXIAL CABLE ANALYSIS
C COAXIAL NEUMANN BOUNDARY ON INNER ALONG WITH DIRICHLET ON OUTER
C INNER SHAPE: CIRCLE OUTER SHAPE: RECTANGULAR--
C************************** QUARTER DOMAIN *****************************
C DOMAIN: FINITE
C ELEMENT: 2-NODED CONSTANT
C EIJI FUKUMORI MARCH 19, 2022
C=======================================================================
IMPLICIT REAL*8 ( A-H , O-Z )
PARAMETER ( ND=2,MXE=50000,MXN=MXE, MXI=10000 )
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".
NEOUTER = 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)
BVN = 0.D0
WRITE(*,*) 'NEUMANN BOUNDARY VALUE =', BINNER
C=======================================================================
C------------------- NORDAL COORDINATES CREATION
C------------------------- OUTER ELEMENTS -------------
DX = RADIUSC/NEOUTER
DXX = (RADIUSC-RADIUSA)/NEOUTER
DY = RADIUSB/NEOUTER
DYY = (RADIUSB-RADIUSA)/NEOUTER
C----------- FIRST SEGMENT STARTS FROM (RADIUSC,0) TO (RADIUSC,RADIUSB)
NODE = 0
DO I = 1 , NEOUTER
NODE = NODE + 1
XCOORD (NODE) = RADIUSC
YCOORD (NODE) = DY * (I-1)
END DO
C---------- SECOND SEGMENT STARTS FROM (RADIUSC,RADIUSB) TO (0,RADIUSB)
DO I = 1 , NEOUTER
NODE = NODE + 1
XCOORD (NODE) = RADIUSC - DX*(I-1)
YCOORD (NODE) = RADIUSB
END DO
C----------- THIRD SEGMENT STARTS FROM (0,RADIUSB) TO (0,0) ----------
DO I = 1 , NEOUTER
NODE = NODE + 1
XCOORD (NODE) = 0.D0
YCOORD (NODE) = RADIUSB - DYY*(I-1)
END DO
C------------ 4TH SEGMENT ALONG INNER CONDUCTOR SURFACE -------------
BASEVEC = DCMPLX( 0.D0, RADIUSA )
UNITANG = -(PI/2.D0) / NEOUTER
UNITVEC = DCMPLX( DCOS(UNITANG), DSIN(UNITANG) )
DO I = 1 , NEOUTER
NODE = NODE + 1
XCOORD (NODE) = DREAL( BASEVEC )
YCOORD (NODE) = DIMAG( BASEVEC )
BASEVEC = BASEVEC * UNITVEC
END DO
C--------- 5TH SEGMENT STARTS FROM (RADIUSA,0) TO (RADIUSC,0) ----------
DO I = 1 , NEOUTER
NODE = NODE + 1
XCOORD (NODE) = RADIUSA+DXX*(I-1)
YCOORD (NODE) = 0.D0
END DO
C=======================================================================
C ELEMENT CREATION
C------- SEGMENT 1 AND 2: OUTER ELENENTS: TWO SEGMENTS ARE DIRICLET-
DO I = 1 , 2*NEOUTER
II = I
JJ = II + 1
NODEX(I,1) = II
NODEX(I,2) = JJ
BV(I,1) = POUTER
BV(I,2) = POUTER
IELTYPE(I) = 1
END DO
C-------------- SEGMENT 3: NATURAL BOUNDARY: VERTICAL SECTION
DO I = 2*NEOUTER+1 , 3*NEOUTER
II = I
JJ = II + 1
NODEX(I,1) = II
NODEX(I,2) = JJ
BV(I,1) = BVN
BV(I,2) = BVN
IELTYPE(I) = 2
END DO
C---------------- SEGMENT 4: INNER CONDUCTOR SURFACE SECTION IS NEUMANN------
DO I = 3*NEOUTER+1 , 4*NEOUTER
II = I
JJ = II + 1
NODEX(I,1) = II
NODEX(I,2) = JJ
BV(I,1) = BINNER
BV(I,2) = BINNER
IELTYPE(I) = 2
END DO
C-------------- SEGMENT 5: NATURAL BOUNDARY: HORIZONTAL SECTION
DO I = 4*NEOUTER+1 , 5*NEOUTER
II = I
JJ = II + 1
IF ( I .EQ. 5*NEOUTER ) JJ = 1
NODEX(I,1) = II
NODEX(I,2) = JJ
BV(I,1) = BVN
BV(I,2) = BVN
IELTYPE(I) = 2
END DO
C=======================================================================
NE = 5*NEOUTER
NNODE = NODE
C=======================================================================
C--------INTERNAL POINTS
NIP = 20
REND = DSQRT(RADIUSB**2+RADIUSC**2)
SINANG = RADIUSB/REND
COSANG = RADIUSC/REND
DR = (REND - RADIUSA) / (NIP+1)
DO I = 1 , NIP
R = RADIUSA + I*DR
YI(I) = R*SINANG
XI(I) = R*COSANG
END DO
C=======================================================================
C MAKING DATA FILES
OPEN ( 1, FILE='BEM0.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
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,*)
END DO
CLOSE (2)
OPEN ( 2, FILE='OBPOINT.DAT', STATUS='UNKNOWN' )
DO IEL = 1 , NE
XM = ( XCOORD(NODEX(IEL,1)) + XCOORD(NODEX(IEL,2)) ) / 2.D0
YM = ( YCOORD(NODEX(IEL,1)) + YCOORD(NODEX(IEL,2)) ) / 2.D0
WRITE (2,*) XM, YM
WRITE (2,*)
END DO
CLOSE (2)
STOP
END