PROGRAM PDOMAIN
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 )
PARAMETER (ND3=(ND*ND-ND)/2+ND, MXEV=5000000,MXNV=6000000 )
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 NODEXVEC(MXEV,ND),XCOORDVEC(MXNV),YCOORDVEC(MXNV),
* AJVEC(MXNV)
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,MXEV,MXNV,NEVEC,NNODEVEC,NODEXVEC,XCOORDVEC,YCOORDVEC)
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,
* MXEV,MXNV,NEVEC,NNODEVEC,NODEXVEC,XCOORDVEC,YCOORDVEC,
* AJVEC )
C=======================================================================
STOP 'NORMAL TERMINATION'
END
C
C
SUBROUTINE INDCTAJ ( MXE,MXN,ND,ND3, NE,NNODE,XCOORD,YCOORD,
* FJ, SS,NODEX,SOURCE,
* MXEV,MXNV,NEVEC,NNODEVEC,NODEXVEC,XCOORDVEC,YCOORDVEC,
* AJVEC )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND), SOURCE(MXE,ND3), SS(ND,ND)
DIMENSION XCOORD(MXN), YCOORD(MXN), FJ(MXE)
DIMENSION NODEXVEC(MXEV,ND),XCOORDVEC(MXNV),YCOORDVEC(MXNV),
* AJVEC(MXNV)
C-------XCOORD AND YCOORD ARE COORDINATES OF CONDUCTOR------
C-------XCOORDVEC AND YCOORDVEC ARE COORDINATES OF observation point------
PI = 4.D0 * DATAN( 1.D0 )
DO IP = 1 , NNODEVEC
AJVEC(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 = XCOORDVEC(IP) - XCOORD(NODEJ)
DY = YCOORDVEC(IP) - YCOORD(NODEJ)
R = DSQRT ( DX*DX + DY*DY )
IF ( R.EQ. 0. ) THEN
WRITE(*,*) ' R = 0.'
R = 1.0D-10
END IF
GREEN = -DLOG(R)/(2.D0*PI)
AJVEC(IP) = AJVEC(IP) + FJ(IEL)*SS(I,J)*GREEN
END DO
END DO
END DO
END DO
C=================================================================
C23456789-123456789-123456789-123456789-123456789-123456789-123456789-12
C----------------- COMPUTATION OF B VECTOR--------------------
IUNIT = 1
OPEN ( IUNIT, FILE='VECTORB.SOL', STATUS='UNKNOWN' )
JUNIT = 2
OPEN ( JUNIT, FILE='VECTORE.SOL', STATUS='UNKNOWN' )
ALPHA = 3.2D0
DO IEL = 1 , NEVEC
I = NODEXVEC(IEL,1)
J = NODEXVEC(IEL,2)
K = NODEXVEC(IEL,3)
L = NODEXVEC(IEL,4)
XP = (XCOORDVEC(I)+XCOORDVEC(J)+XCOORDVEC(K)+XCOORDVEC(L))/4.D0
YP = (YCOORDVEC(I)+YCOORDVEC(J)+YCOORDVEC(K)+YCOORDVEC(L))/4.D0
DX = XCOORDVEC(J) - XCOORDVEC(I)
DY = YCOORDVEC(K) - YCOORDVEC(J)
DADX = ( (AJVEC(J)-AJVEC(I))/DX+(AJVEC(K)-AJVEC(L))/DX )/2.D0
DADY = ( (AJVEC(K)-AJVEC(J))/DY+(AJVEC(L)-AJVEC(I))/DY )/2.D0
BX = DADY
BY = - DADX
CALL VECPLT ( IUNIT, XP, YP, BX, BY, ALPHA )
EX = - DADX
EY = - DADY
CALL VECPLT ( JUNIT, XP, YP, EX, EY, ALPHA )
END DO
CLOSE (IUNIT)
CLOSE (JUNIT)
OPEN ( 1, FILE='CONTOUR.SOL', STATUS='UNKNOWN' )
CALL CONTOUR (MXEV,MXNV,ND,NEVEC,NNODEVEC,NODEXVEC,
* XCOORDVEC,YCOORDVEC, AJVEC )
CLOSE (1)
RETURN
END
C
C
SUBROUTINE CONTOUR (MXEV,MXNV,ND,NEVEC,NNODEVEC,NODEXVEC,
* XCOORDVEC,YCOORDVEC, P )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEXVEC(MXEV,ND), XCOORDVEC(MXNV),
* YCOORDVEC(MXNV), P(MXNV),
* S(4), B(2,4)
CCCCCCC IF ( PPMAX .EQ. PPMIN ) RETURN
NSTEP = 20
NSTEP = NSTEP/2*2 + 1
PPMAX = P(1)
PPMIN = P(1)
DO I = 2 , NNODEVEC
IF ( P(I) .GT. PPMAX ) PPMAX = P(I)
IF ( P(I) .LT. PPMIN ) PPMIN = P(I)
END DO
WRITE (*,*) PPMAX , PPMIN
DS = ( PPMAX - PPMIN ) / NSTEP
C
IC = 0
DO IEL = 1 , NE
DO I = 1 , ND
B(1,I) = XCOORDVEC(NODEXVEC(IEL,I))
B(2,I) = YCOORDVEC(NODEXVEC(IEL,I))
S(I) = P(NODEXVEC(IEL,I))
END DO
IF ( (PPMIN .LT. 0.D0) .AND. (PPMAX.GT.0.D0) ) THEN
CALL PLTSAI ( DS, NSTEP, -NSTEP*DS, B, S )
CALL PLTSAI ( DS, 1, 0.D0, B, S )
CALL PLTSAI ( DS, NSTEP, DS, B, S )
ELSE
CALL PLTSAI ( DS, NSTEP, PPMIN, B, S )
END IF
END DO
RETURN
END
C
C
SUBROUTINE VECPLT ( IUNIT, X0, Y0, U, V, FACT )
IMPLICIT REAL*8 ( A-H , O-Z )
DATA AL, BETA / 0.65D0, 0.15D0 /
DX = U * FACT
DY = V * FACT
WRITE(IUNIT,*) X0, Y0
WRITE(IUNIT,*) X0+ AL*DX , Y0+ AL*DY
RNX = DY
RNY = -DX
X = X0 + AL*DX + BETA*RNX
Y = Y0 + AL*DY + BETA*RNY
WRITE (IUNIT,*) X, Y
WRITE(IUNIT,*) X0+DX , Y0+DY
X = X0 + AL*DX - BETA*RNX
Y = Y0 + AL*DY - BETA*RNY
WRITE (IUNIT,*) X, Y
WRITE(IUNIT,*) X0+ AL*DX , Y0+ AL*DY
WRITE(IUNIT,*)
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,MXEV,MXNV,NEVEC,NNODEVEC,NODEXVEC,XCOORDVEC,YCOORDVEC)
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),FJ(MXE)
DIMENSION NODEXVEC(MXEV,ND),XCOORDVEC(MXNV),YCOORDVEC(MXNV)
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
CLOSE (1)
C========> VECTOR PLOT ELEMENTS
OPEN ( 2, FILE = 'VECTORCG.DAT', STATUS = 'UNKNOWN' )
READ(2,*) NEVEC
DO I = 1 , NEVEC
READ (2,*) IEL,(NODEXVEC(IEL,J),J=1,ND)
END DO
READ (2,*) NNODEVEC
DO I = 1 , NNODEVEC
READ (2,*) NODE, XCOORDVEC(NODE), YCOORDVEC(NODE)
END DO
CLOSE (2)
RETURN
END
C
C
SUBROUTINE PLTSAI ( DS, NSTEP, START, CRD, SS )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION CRD(2,4), SS(4), X(4), Y(4), S(4)
IF ( NSTEP .EQ. 0 ) RETURN
X(3) = ( CRD(1,1) + CRD(1,2) + CRD(1,3) + CRD(1,4) ) / 4.
Y(3) = ( CRD(2,1) + CRD(2,2) + CRD(2,3) + CRD(2,4) ) / 4.
S(3) = ( SS(1) + SS(2) + SS(3) + SS(4) ) / 4.
SMAX = DMAX1 ( SS(1), SS(2), SS(3), SS(4) )
SMIN = DMIN1 ( SS(1), SS(2), SS(3), SS(4) )
DO 52 LEVEL = 1 , NSTEP
SXY = START + (LEVEL-1) * DS
IF ( (SMAX-SXY)*(SMIN-SXY) .LT. 0 ) THEN
DO 60 IEL = 1 , 4
X(1) = CRD(1,IEL)
Y(1) = CRD(2,IEL)
S(1) = SS(IEL)
IF ( IEL .EQ. 4 ) THEN
X(2) = CRD(1,1)
Y(2) = CRD(2,1)
S(2) = SS(1)
ELSE
X(2) = CRD(1,IEL+1)
Y(2) = CRD(2,IEL+1)
S(2) = SS (IEL+1)
ENDIF
X(4) = X(1)
Y(4) = Y(1)
S(4) = S(1)
K = 0
DO 70 ISG = 1 , 3
IF ( S(ISG ) .LT. SXY ) GO TO 30
IF ( S(ISG+1) .LT. SXY ) GO TO 40
GO TO 70
30 IF ( S(ISG+1) .LT. SXY ) GO TO 70
40 T = ( SXY - S(ISG) ) / ( S(ISG+1) - S(ISG) )
X0 = X(ISG+1)*T + (1.- T)*X(ISG)
Y0 = Y(ISG+1)*T + (1.- T)*Y(ISG)
IF ( K .EQ. 0 ) GO TO 71
CALL XDRAW ( X0, Y0 )
GO TO 60
71 CALL XMOVE ( X0, Y0 )
K = 1
70 CONTINUE
60 CONTINUE
ENDIF
52 CONTINUE
RETURN
END
C
C
SUBROUTINE XMOVE ( X0, Y0 )
IMPLICIT REAL*8 ( A-H , O-Z )
WRITE (1,*)
WRITE (1,*) X0, Y0
RETURN
END
C
C
SUBROUTINE XDRAW ( X0, Y0 )
IMPLICIT REAL*8 ( A-H , O-Z )
WRITE (1,*) X0, Y0
RETURN
END