PROGRAM SETUAZ
C=======================================================================
C DATA GENERATING PROGRAM FOR BEM8LINU.FOR
C PROJECT: DOUBLE-LET PROJECT NAME: SURFACE BOUNDARY CONDITION
C DOMAIN: INFINITE
C ELEMENT: FOUR-NODED ISOPARAMETRIC ELEMENT
C DOMAIN DISCRETIZATION: UNEVEN ELEMENTS WITH VERTICAL SCAN
C EIJI FUKUMORI DECEMBER 29, 1993
C=======================================================================
IMPLICIT REAL*8 ( A-H , O-Z )
PARAMETER ( ND=2,MXE=500,MXN=MXE, MXI=100 )
PARAMETER ( NDVEC=4,MXEVEC=5000000,MXNVEC=5500000 )
C=======================================================================
DIMENSION NODEX(MXE,ND),XCOORD(MXN), YCOORD(MXN)
DIMENSION XI(MXI), YI(MXI), IELTYPE(MXE),BV(MXE,ND)
DIMENSION NODEXVEC(MXEVEC,NDVEC), XCOORDVEC(MXNVEC)
DIMENSION YCOORDVEC(MXNVEC)
C=======================================================================
CCCCCCCCCCCCCCCCCCCCCCCCC NC=16
NC=64
FMU = 1.D0
POTENTIL = 0.5D0
WRITE (*,*) 'Prescribed potential Value on The boundary=',POTENTIL
C=======================================================================
PI = 4.* DATAN( 1.D0)
RADIUS = 1.0D0
XC1 = 0.D0
WRITE (*,*) 'RADIUS = ', RADIUS
WRITE(*,*) 'TYPE IN GAP BETWEEN WIRE ELEMENTS'
READ (*,*) GAP
YC1 = RADIUS + GAP/2.D0
IF ( GAP .LE. 0.D0 ) STOP 'A VALUE OF GAP MUST BE MORE THAN ZERO'
XC2 = XC1
YC2 = -YC1
C=======================================================================
C----- NORDAL COORDINATES CREATION
SEG = NC
HSEG = SEG/2.D0
DO I = 1 , NC
XCOORD(I) = RADIUS*COS ( -I*PI/HSEG )
YCOORD(I) = RADIUS*SIN ( -I*PI/HSEG )
END DO
DO I = 1 , NC
XCOORD(I+NC) = XCOORD(I) + XC2
YCOORD(I+NC) = YCOORD(I) + YC2
XCOORD(I) = XCOORD(I) + XC1
YCOORD(I) = YCOORD(I) + YC1
END DO
NNODE = NC*2
C=======================================================================
C ELEMENT CREATION
C---------------- CIRCUMFERENCE = 2*PI*RADIUS
C MAGNETIC FLUX DENSITY ON SURFACE = ELECTRIC CURRENT/(CIRCUMFERENCE)
C-------- DIRECTIN = INTEGRAL DIRECTION ------
DIRECTIN = -1.D0
DO I = 1 , NC
J = I+1
IF ( I .EQ. NC ) J = 1
NODEX(I,1) = I
NODEX(I,2) = J
IELTYPE(I) = 1
BV(I,1) = POTENTIL
BV(I,2) = POTENTIL
END DO
DO I = 1+NC , NC*2
J = I+1
IF ( I .EQ. NC*2 ) J = 1+NC
NODEX(I,1) = I
NODEX(I,2) = J
IELTYPE(I) = 1
BV(I,1) = -POTENTIL
BV(I,2) = -POTENTIL
END DO
NE = NC*2
C=======================================================================
C--------INTERNAL POINTS
NIP = 0
NIP = NIP + 1
XI(NIP) = 0.D0
YI(NIP) = 0.0001D0
NIP = NIP + 1
XI(NIP) = 0.
YI(NIP) = (YC1-RADIUS)/4.D0
NIP = NIP + 1
XI(NIP) = 0.D0
YI(NIP) = (YC1-RADIUS)/2.D0
NIP = NIP + 1
XI(NIP) = 0.D0
YI(NIP) = (YC1-RADIUS)/4.D0*3.D0
NIP = NIP + 1
XI(NIP) = 0.D0
YI(NIP) = YC1 - RADIUS/4.D0*3.D0
NIP = NIP + 1
XI(NIP) = 0.D0
YI(NIP) = YC1 - RADIUS/4.D0
NIP = NIP + 1
XI(NIP) = 0.D0
YI(NIP) = YC1 + RADIUS/4.D0
NIP = NIP + 1
XI(NIP) = 0.
YI(NIP) = YC1 + RADIUS/4.D0*3.D0
NIP = NIP + 1
XI(NIP) = 0.D0
YI(NIP) = YC1 + RADIUS*1.D0 + RADIUS/4.D0
NIP = NIP + 1
XI(NIP) = 0.D0
YI(NIP) = YC1 + RADIUS*2.D0
NIP = NIP + 1
XI(NIP) = 0.D0
YI(NIP) = YC1 + RADIUS*3.D0
NIP = NIP + 1
XI(NIP) = 0.D0
YI(NIP) = YC1 + RADIUS*4.D0
NIP = NIP + 1
XI(NIP) = 0.D0
YI(NIP) = YC1 + RADIUS*5.D0
C
DO I = 3 , 30
NIP = NIP + 1
XI(NIP) = 0.
YI(NIP) = YC1 + RADIUS*I*2.D0
END DO
C
NNN = NIP
DO I = 1 , NNN
NIP = NIP + 1
XI(NIP) = 0.D0
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='BEM1.DAT', STATUS='UNKNOWN' )
WRITE (1,*) FMU
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
NNODE = NE
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)
C=======================================================================
C
C ELEMENT CREATION FOR VECTOR PLOT GRAPHICS
C
C=======================================================================
DL = RADIUS / 2.5D0
XLENGTH = 8*RADIUS
YLENGTH = 6*RADIUS + GAP
NEX = XLENGTH / DL
NEY = YLENGTH / DL
NEX = NEX/2*2
NEY = NEY/2*2
DX = XLENGTH / NEX
DY = YLENGTH / NEY
NDX = NEX + 1
NDY = NEY + 1
NEVEC = 0
DO I = 1 , NEY
DO J = 1 , NEX
NEVEC = NEVEC + 1
IF ( NEVEC .GT. MXEVEC ) STOP 'NEVEC > MXE'
NODEXVEC(NEVEC,1) = NDX*(I-1) + J
NODEXVEC(NEVEC,2) = NODEXVEC(NEVEC,1) + 1
NODEXVEC(NEVEC,3) = NODEXVEC(NEVEC,2) +NDX
NODEXVEC(NEVEC,4) = NODEXVEC(NEVEC,1)+NDX
END DO
END DO
C--------------- NODAL INFORMATION -------------------
NNODEVEC = 0
DO J = 1 , NDY
DO I = 1 , NDX
NNODEVEC = NNODEVEC + 1
XCOORDVEC(NNODEVEC) = (I-(NEX/2+1))*DX
YCOORDVEC(NNODEVEC) = (J-(NEY/2+1))*DY
END DO
END DO
WRITE(*,*) 'NEVEX=',NEVEC
WRITE(*,*) 'NNODEVEC=',NNODEVEC
C--------------MAKING INPUT FILE------------------------
OPEN ( 5, FILE='VECTORCG.DAT', STATUS = 'UNKNOWN' )
WRITE(5,*) NEVEC
DO I = 1 , NEVEC
WRITE (5,*) I,(NODEXVEC(I,J),J=1,NDVEC)
END DO
WRITE(5,*) NNODEVEC
DO I = 1 , NNODEVEC
WRITE(5,*) I, XCOORDVEC(I), YCOORDVEC(I)
END DO
CLOSE (5)
C=======================================================================
STOP
END