PROGRAM DOMAIN
C=======================================================================
C 2-DIM FEM PROGRAM FOR SOLVING POISSON EQUATION WITH
C INFINITE DOMAIN AND EXTERNAL SOURCE TERM
C USING UPPER HALF BANDED MATRIX
C EQUATION: DP/DX + DP/DY = F(X,Y)
C ELEMENT : 4-NODED ISO-PARAMETERIC
C NUMBERING ORDER: (-1,-1),(+1,-1),(+1,+1),(-1,+1)
C ORIGINAL:1984 EIJI FUKUMORI BUFFALO NY & REVISED: 1994 ACHI
C=======================================================================
IMPLICIT REAL*8 ( A-H , O-Z )
PARAMETER (ND=4,MXE=50000,MXN=52000,INTEPT=2 , MXNC=10000)
PARAMETER ( MXI=10000, ND3=(ND*ND-ND)/2+ND )
DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN), SK(ND,ND),
* X(ND), Y(ND),BX(2,ND),SAI(INTEPT), W(INTEPT),
* BPP(2,ND,INTEPT,INTEPT),
* SOURCE(MXE,ND3), SS(ND,ND), SF(ND,INTEPT,INTEPT)
DIMENSION F1(ND), F2(ND), E1(ND), E2(ND), BP(2,ND,ND), FJ(MXE)
DIMENSION XI(MXI), YI(MXI), AJ(MXI), ICIR(MXNC)
C=======================================================================
CALL GRULE ( INTEPT, SAI, W )
CALL DERIV ( ND, INTEPT, X, Y, SAI, BPP )
CALL SHAPEF( ND, INTEPT, X, SAI, SF )
C=======================================================================
CALL INPUT (MXNC, ND,MXE,MXN,NE,NNODE,NODEX,XCOORD,YCOORD,
* FJ,MXI,NIP,XI,YI, NC, ICIR)
C=======================================================================
CALL GDM ( MXE,MXN,INTEPT,ND,ND3,BPP,W,NE,SF,
* XCOORD,YCOORD,NODEX,SS, SOURCE)
C=======================================================================
CALL INDCTAJ ( MXE,MXN,ND,ND3, NE,NNODE,XCOORD,YCOORD,
* FJ, SS,NODEX,SOURCE,MXI,NIP,XI,YI, AJ, AJMIN, AJMAX )
C=======================================================================
CALL WIREAJ ( MXE,MXN,ND,ND3, NE,NNODE,XCOORD,YCOORD,
* FJ, SS,NODEX,SOURCE,MXI,NC,XI,YI, AJ, AJAVE, AJMINS )
CALL AREACMP ( MXE,MXN,ND,ND3, NE,NNODE,XCOORD,YCOORD,
* SS,NODEX,SOURCE )
C=======================================================================
DELTAAJ = AJMAX - AJMIN
OPEN ( 3, FILE='DIFFEREN.CEP', STATUS='UNKNOWN' )
WRITE (3,*) 'NUMBER OF ELEMENTS=', NE
WRITE (3,*) 'POTENTIAL-MIN(WITHIN CONDUCTOR)=', AJMIN
WRITE (3,*) 'POTENTIAL-MAX(WITHIN CONDUCTOR)=', AJMAX
WRITE (3,*) 'DIFFERENCE IN MIN AND MAX =', DELTAAJ
WRITE (3,*) 'AVERAGE VALUE OF POTENTIAL ON SURFACE*2 =', AJAVE*2
WRITE (3,*) 'MINIMUM VALUE OF POTENTIAL ON SURFACE*2 =', AJMINS*2
CLOSE (3)
STOP 'NORMAL TERMINATION'
END
C
C
SUBROUTINE AREACMP ( MXE,MXN,ND,ND3, NE,NNODE,XCOORD,YCOORD,
* SS,NODEX,SOURCE )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND), SOURCE(MXE,ND3), SS(ND,ND)
DIMENSION XCOORD(MXN), YCOORD(MXN)
OPEN ( 1, FILE='AREA.SOL', STATUS='UNKNOWN' )
AREA = 0.D0
DO IEL = 1 ,NE
M = 0
DO K = 1 , ND
DO L = K , ND
M = M + 1
SS(K,L) = SOURCE(IEL,M)
SS(L,K) = SS(K,L)
END DO
END DO
DO I = 1 , ND
DO J = 1 , ND
NODEJ = NODEX(IEL,J)
AREA = AREA + SS(I,J)
END DO
END DO
END DO
WRITE (1,*) 'CROSS-SECTIONAL-AREA= ', AREA
CLOSE (1)
RETURN
END
C
C
SUBROUTINE WIREAJ ( MXE,MXN,ND,ND3, NE,NNODE,XCOORD,YCOORD,
* FJ, SS,NODEX,SOURCE,MXI,NC,XI,YI, AJ, AJAVE, AJMINS )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND), SOURCE(MXE,ND3), SS(ND,ND), AJ(MXI)
DIMENSION XCOORD(MXN), YCOORD(MXN), FJ(MXE), XI(MXI), YI(MXI)
C
OPEN ( 4, FILE='DPDN.DAT', STATUS='UNKNOWN' )
READ (4,*) NCNC, XC1, YC1
DO I = 1 , NCNC
READ (4,*) XI(I), YI(I)
END DO
CLOSE (4)
OPEN ( 1, FILE='DPDN.SOL', STATUS='UNKNOWN' )
WRITE(1,*)'X-COORDINATE Y-COORDINATE ANGLE(RAD) DP/DN P1 P2'
PI = 4.D0 * DATAN( 1.D0 )
DO IP = 1 , NCNC
AJ(IP) = 0.D0
DO IEL = 1 ,NE
M = 0
DO K = 1 , ND
DO L = K , ND
M = M + 1
SS(K,L) = SOURCE(IEL,M)
SS(L,K) = SS(K,L)
END DO
END DO
DO I = 1 , ND
DO J = 1 , ND
NODEJ = NODEX(IEL,J)
DX = XI(IP) - XCOORD(NODEJ)
DY = YI(IP) - YCOORD(NODEJ)
R = DSQRT ( DX*DX + DY*DY )
IF ( R.EQ. 0. ) THEN
WRITE(*,*) ' R = 0.'
R = 1.0D-13
END IF
GREEN = -DLOG(R)/(2.D0*PI)
AJ(IP) = AJ(IP) + FJ(IEL)*SS(I,J)*GREEN
END DO
END DO
END DO
END DO
AJAVE = 0.D0
AJMINS = 1.D10
NC = NCNC/2
DO I = 1 , NC
I2 = I*2
I1 = I2 - 1
DX = XI(I2)-XI(I1)
DY = YI(I2)-YI(I1)
DN = DSQRT ( DX**2 + DY**2 )
DADN = - (AJ(I2)-AJ(I1))/DN
WRITE (1,*) XI(I1),YI(I1),ANGLE(PI,XI(I1)-XC1,YI(I1)-YC1),
* DADN, AJ(I1), AJ(I2)
AJAVE = AJAVE + AJ(I1)
IF ( AJ(I1) .LT. AJMINS ) AJMINS = AJ(I1)
END DO
AJAVE = AJAVE / NC
CLOSE (1)
RETURN
END
C
C
FUNCTION ANGLE (PI,X,Y)
IMPLICIT REAL*8 ( A-H , O-Z )
IF ( Y .EQ. 0.D0 .AND. X .GT. 0.D0 ) THEN
ANGLE=0.D0
RETURN
END IF
IF ( Y .GT. 0.D0 .AND. X .EQ. 0.D0 ) THEN
ANGLE=0.5D0*PI
RETURN
END IF
IF ( Y .EQ. 0.D0 .AND. X .LT. 0.D0 ) THEN
ANGLE=PI
RETURN
END IF
IF ( Y .LT. 0.D0 .AND. X .EQ. 0.D0 ) THEN
ANGLE=1.5D0*PI
RETURN
END IF
IF ( X .NE. 0.D0 ) ANG = DATAN (DABS(Y)/DABS(X))
IF ( Y .GT. 0.D0 .AND. X .GT. 0.D0 ) THEN
ANGLE=ANG
RETURN
END IF
IF ( Y .GT. 0.D0 .AND. X .LT. 0.D0 ) THEN
ANGLE=PI-ANG
RETURN
END IF
IF ( Y .LT. 0.D0 .AND. X .LT. 0.D0 ) THEN
ANGLE=PI+ANG
RETURN
END IF
IF ( Y .LT. 0.D0 .AND. X .GT. 0.D0 ) THEN
ANGLE=2.D0*PI-ANG
RETURN
END IF
RETURN
END
C
C
SUBROUTINE INDCTAJ ( MXE,MXN,ND,ND3, NE,NNODE,XCOORD,YCOORD,
* FJ, SS,NODEX,SOURCE,MXI,NIP,XI,YI, AJ, AJMIN, AJMAX )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND), SOURCE(MXE,ND3), SS(ND,ND), AJ(MXI)
DIMENSION XCOORD(MXN), YCOORD(MXN), FJ(MXE), XI(MXI), YI(MXI)
C
OPEN ( 1, FILE='INTERNAL.SOL', STATUS='UNKNOWN' )
WRITE(1,*)' X-COORDINATE Y-COORDINATE POTENTIAL'
PI = 4.D0 * DATAN( 1.D0 )
DO IP = 1 , NIP
AJ(IP) = 0.D0
DO IEL = 1 ,NE
M = 0
DO K = 1 , ND
DO L = K , ND
M = M + 1
SS(K,L) = SOURCE(IEL,M)
SS(L,K) = SS(K,L)
END DO
END DO
DO I = 1 , ND
DO J = 1 , ND
NODEJ = NODEX(IEL,J)
DX = XI(IP) - XCOORD(NODEJ)
DY = YI(IP) - YCOORD(NODEJ)
R = DSQRT ( DX*DX + DY*DY )
IF ( R.EQ. 0. ) THEN
WRITE(*,*) ' R = 0.'
R = 1.0D-13
END IF
GREEN = -DLOG(R)/(2.D0*PI)
AJ(IP) = AJ(IP) + FJ(IEL)*SS(I,J)*GREEN
END DO
END DO
END DO
WRITE(1,*) XI(IP), YI(IP), AJ(IP)
END DO
CLOSE (1)
C========= search FOR MAX AND MIN IN AJ(I)
AJMAX = AJ(1)
AJMIN = AJ(1)
DO I = 2, NIP
IF ( AJ(I) .GT. AJMAX ) AJMAX = AJ(I)
IF ( AJ(I) .LT. AJMIN ) AJMIN = AJ(I)
END DO
RETURN
END
C
C
SUBROUTINE GRULE ( N , SAI , W )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION SAI(N) , W(N)
IF ( N .LT. 2 ) STOP'N<2'
IF ( N .GT. 6 ) STOP'N>6'
GO TO ( 99, 20, 30, 40, 50, 60 ) , N
99 STOP
20 SAI(1) = DSQRT(3.D0)/3.D0
W(1) = 1.D0
GO TO 88
30 SAI(1) = DSQRT(15.D0)/5.D0
SAI(2) = 0.D0
W(1) = 5.D0/ 9.D0
W(2) = 8.D0/ 9.D0
GO TO 88
40 SAI(1) = 0.33998104358485D0
SAI(2) = 0.86113631159405D0
W(1) = 0.65214515486254D0
W(2) = 0.34785484513745D0
GO TO 88
50 SAI(1) = 0.90617984593866D0
SAI(2) = 0.53846931010568D0
SAI(3) = 0.D0
W(1) = 0.23692688505619D0
W(2) = 0.47862867049937D0
W(3) = 5.12D0 / 9.D0
GO TO 88
60 SAI(1) = 0.23861918608320D0
SAI(2) = 0.66120938646626D0
SAI(3) = 0.93246951420315D0
W(1) = 0.46791393457269D0
W(2) = 0.36076157304814D0
W(3) = 0.17132449237917D0
88 NN = N / 2
DO 11 I = 1 , NN
J = N - I + 1
SAI(J) = - SAI(I)
W(J) = W(I)
11 CONTINUE
RETURN
END
C
C
SUBROUTINE ISOPARA ( ND, E1 , E2 , F )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION F(ND)
F(1) = 0.25*(1.- E1)*(1.- E2)
F(2) = 0.25*(1.+ E1)*(1.- E2)
F(3) = 0.25*(1.+ E1)*(1.+ E2)
F(4) = 0.25*(1.- E1)*(1.+ E2)
RETURN
END
C
C
SUBROUTINE GDM ( MXE,MXN,INTEPT,ND,ND3,BPP,W,NE,SF,
* XCOORD,YCOORD,NODEX,SS, SOURCE )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),SOURCE(MXE,ND3)
DIMENSION BPP(2,ND,INTEPT,INTEPT),W(INTEPT),SF(ND,INTEPT,INTEPT)
DIMENSION SS(ND,ND)
C
DO IEL = 1 ,NE
DO I = 1 , ND
DO J = I , ND
SS(I,J) = 0.D0
END DO
END DO
C
DO K = 1 , INTEPT
DO L = 1 , INTEPT
WEIGHT = W(K) * W(L)
YAC11 = 0.
YAC12 = 0.
YAC21 = 0.
YAC22 = 0.
DO I = 1 , ND
YAC11 = YAC11 + BPP(1,I,K,L) * XCOORD(NODEX(IEL,I))
YAC12 = YAC12 + BPP(1,I,K,L) * YCOORD(NODEX(IEL,I))
YAC21 = YAC21 + BPP(2,I,K,L) * XCOORD(NODEX(IEL,I))
YAC22 = YAC22 + BPP(2,I,K,L) * YCOORD(NODEX(IEL,I))
END DO
DETJ = YAC11 * YAC22 - YAC12 * YAC21
BETA = WEIGHT * DETJ
C
DO I = 1 , ND
DO J = I , ND
SS(I,J) = SS(I,J) + SF(I,K,L) * SF(J,K,L) * BETA
END DO
END DO
C
END DO
END DO
C
M = 0
DO K = 1 , ND
DO L = K , ND
M = M + 1
SOURCE(IEL,M) = SS(K,L)
END DO
END DO
C
END DO
RETURN
END
C
C
SUBROUTINE SHAPEF ( ND, INTEPT, F, SAI, SF )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION F(ND) , SAI(INTEPT) , SF(ND,INTEPT,INTEPT)
DO K = 1 , INTEPT
E1 = SAI (K)
DO L = 1 , INTEPT
E2 = SAI( L )
CALL ISOPARA ( ND, E1 , E2 , F )
DO I = 1 , ND
SF(I,K,L) = F(I)
END DO
END DO
END DO
RETURN
END
C
C
SUBROUTINE DERIV ( ND, INTEPT, F0, F1, SAI, BPP )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION SAI(INTEPT),BPP(2,ND,INTEPT,INTEPT),F0(ND),F1(ND)
DO 40 K = 1 , INTEPT
E1 = SAI (K)
DO 30 L = 1 , INTEPT
E2 = SAI (L)
CALL ISOPARA ( ND, E1+0.5D0 , E2 , F1 )
CALL ISOPARA ( ND, E1-0.5D0 , E2 , F0 )
DO I = 1 , ND
BPP(1,I,K,L) = F1(I) - F0(I)
END DO
CALL ISOPARA ( ND, E1 , E2+0.5D0 , F1 )
CALL ISOPARA ( ND, E1 , E2-0.5D0 , F0 )
DO I = 1 , ND
BPP(2,I,K,L) = F1(I) - F0(I)
END DO
30 CONTINUE
40 CONTINUE
RETURN
END
C
C
SUBROUTINE INPUT( MXNC,ND,MXE,MXN, NE,NNODE,NODEX,XCOORD,YCOORD,
* FJ,MXI,NIP,XI,YI, NC, ICIR)
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),FJ(MXE)
DIMENSION XI(MXI), YI(MXI), ICIR(MXNC)
C========> FILE OPEN
OPEN ( 1, FILE = 'DOMAIN.DAT', STATUS = 'UNKNOWN' )
C========> ELEMENTS
READ (1,*) NE
WRITE(*,*)'NUMBER OF ELEMENTS =', NE
IF ( NE .GT. MXE ) STOP 'NE > MXE'
DO I = 1 , NE
READ (1,*) IEL,(NODEX(IEL,J),J=1,ND), FJ(IEL)
END DO
C========> NODAL COORDINATES
READ (1,*) NNODE
WRITE(*,*)'NUMBER OF NODAL POINTS =', NNODE
IF ( NNODE .GT. MXN ) STOP 'NNODE > MXN'
DO I = 1 , NNODE
READ (1,*) NODE,XCOORD(NODE),YCOORD(NODE)
END DO
C========> INTERNAL POINTS
READ(1,*) NIP
IF ( NIP .GE. 1 ) THEN
DO I = 1 , NIP
READ (1,*) XI(I), YI(I)
END DO
END IF
READ (1,*) NC
DO I = 1 , NC
READ (1,*) ICIR(I)
END DO
CLOSE (1)
RETURN
END