DELAUNAY.FOR
************************************************************************
*----------------------------MAIN PROGRAM------------------------------*
************************************************************************
CALL INPUT
CALL INITIAL
CALL NEWPOINT
END
************************************************************************
SUBROUTINE SETUP
************************************************************************
LOGICAL LI,LT,LL1,LL2,INDOMAIN
COMMON/TT/KTSTACK(5500,3),APOINT(5500,2),AINP(500,2),B(5500),
&NDOMAIN(400),RB(500)
COMMON/NN/NS,NPP,NADD,NBOUND,NPMAX,NCIRCLE,X0,Y0,R,M,N,
&INSERTNO,XSCAN,YSCAN,INDOMAIN,LT,XT,XT1,XT2,XT3,YT,YT1,YT2,YT3
************************************************************************
*
ENTRY INITIAL
NS=2
KTSTACK(1,1)=NPP-3
KTSTACK(1,2)=NPP-2
KTSTACK(1,3)=NPP-1
KTSTACK(2,1)=NPP-1
KTSTACK(2,2)=NPP
KTSTACK(2,3)=NPP-3
DO 100 INSERTNO=1,NPP-4
WRITE(*,*)'INSERTNO',INSERTNO
I=INSERTNO
CALL INSERT
100 CONTINUE
CALL MODIFY
CALL BACKGROUND
NADD=NS
CALL SORTING
CALL OUTPUT1
RETURN
*
ENTRY NEWPOINT
DO 200 KINSERT=1,NPMAX
DO 300 NVE=1,NS
NCIRCLE=NVE
CALL CIRCLE
DO 400 J=1,NPP
IF(ABS(X0-APOINT(J,1)).LT.1E-20.AND.ABS(Y0-APOINT(J,2)).LT.1E-20)
&THEN
GOTO 300
END IF
400 CONTINUE
XSCAN=X0
YSCAN=Y0
CALL SCAN
IF(INDOMAIN) THEN
GOTO 500
ELSE
WRITE(*,*)'I AM VERY SORRY!'
GOTO 300
END IF
300 CONTINUE
500 CONTINUE
NPP=NPP+1
WRITE(*,*)'I=',NPP,NS
WRITE(*,*)'DELTA=',B(NVE),B(NS),B(1)
APOINT(NPP,1)=X0
APOINT(NPP,2)=Y0
INSERTNO=NPP
CALL INSERT
CALL SORTING
IF(MOD(KINSERT,400).EQ.0.OR.KINSERT.EQ.NPMAX) CALL OUTPUT2
200 CONTINUE
RETURN
END
************************************************************************
SUBROUTINE INSERT
************************************************************************
LOGICAL LI,LT,LL1,LL2,INDOMAIN
COMMON/TT/KTSTACK(5500,3),APOINT(5500,2),AINP(500,2),B(5500),
&NDOMAIN(400),RB(500)
COMMON/NN/NS,NPP,NADD,NBOUND,NPMAX,NCIRCLE,X0,Y0,R,M,N,
&INSERTNO,XSCAN,YSCAN,INDOMAIN,LT,XT,XT1,XT2,XT3,YT,YT1,YT2,YT3
************************************************************************
DIMENSION MPOINT(5500)
*
NPOINT=0
XX=APOINT(INSERTNO,1)
YY=APOINT(INSERTNO,2)
XT=XX
YT=YY
DO 3100 J=1,NS
N1=KTSTACK(J,1)
N2=KTSTACK(J,2)
N3=KTSTACK(J,3)
XT1=APOINT(N1,1)
YT1=APOINT(N1,2)
XT2=APOINT(N2,1)
YT2=APOINT(N2,2)
XT3=APOINT(N3,1)
YT3=APOINT(N3,2)
CALL INTRIANGLE
IF(LT) THEN
MPOINT(1)=N1
MPOINT(2)=N2
MPOINT(3)=N3
NPOINT=3
DO 3200 JJ=J,NS-1
DO 3200 JM=1,3
KTSTACK(JJ,JM)=KTSTACK(JJ+1,JM)
3200 B(JJ)=B(JJ+1)
NS=NS-1
GOTO 3300
END IF
3100 CONTINUE
3300 CONTINUE
KBAK=1
3400 IF(KBAK.EQ.NPOINT) GOTO 3500
DO 3600 JK=KBAK,NPOINT
L1=MPOINT(JK)
L2=MPOINT(JK+1)
IF(JK.EQ.NPOINT) L2=MPOINT(1)
NSS=NS
IF (NSS.EQ.0) GOTO 3500
DO 3700 J=1,NSS
M1=KTSTACK(J,1)
M2=KTSTACK(J,2)
M3=KTSTACK(J,3)
LL1=.FALSE.
LL2=.FALSE.
IF(L1.EQ.M1.OR.L1.EQ.M2.OR.L1.EQ.M3) LL1=.TRUE.
IF(L2.EQ.M1.OR.L2.EQ.M2.OR.L2.EQ.M3) LL2=.TRUE.
IF(LL1.AND.LL2) THEN
IF(M1.NE.L1.AND.M1.NE.L2) MNEW=M1
IF(M2.NE.L1.AND.M2.NE.L2) MNEW=M2
IF(M3.NE.L1.AND.M3.NE.L2) MNEW=M3
NCIRCLE=J
CALL CIRCLE
LI=.FALSE.
RR=(XX-X0)**2.+(YY-Y0)**2.
IF(RR.LE.(R**2.)) LI=.TRUE.
IF(LI) THEN
JK0=JK+1
DO 3800 K1=NPOINT,JK0,-1
MPOINT(K1+1)=MPOINT(K1)
3800 CONTINUE
MPOINT(JK0)=MNEW
NPOINT=NPOINT+1
KBAK=JK
DO 3900 JJ=J,NS-1
DO 3900 JM=1,3
KTSTACK(JJ,JM)=KTSTACK(JJ+1,JM)
3900 B(JJ)=B(JJ+1)
NS=NS-1
GOTO 3400
END IF
END IF
3700 CONTINUE
3600 CONTINUE
3500 CONTINUE
NADD=0
DO 4000 J=1,NPOINT
NS=NS+1
NADD=NADD+1
KTSTACK(NS,1)=MPOINT(J)
KTSTACK(NS,2)=MPOINT(J+1)
KTSTACK(NS,3)=INSERTNO
IF(J.EQ.NPOINT) THEN
KTSTACK(NS,2)=MPOINT(1)
END IF
4000 CONTINUE
RETURN
END
************************************************************************
SUBROUTINE SUPPLY
************************************************************************
LOGICAL LI,LT,LL1,LL2,INDOMAIN
COMMON/TT/KTSTACK(5500,3),APOINT(5500,2),AINP(500,2),B(5500),
&NDOMAIN(400),RB(500)
COMMON/NN/NS,NPP,NADD,NBOUND,NPMAX,NCIRCLE,X0,Y0,R,M,N,
&INSERTNO,XSCAN,YSCAN,INDOMAIN,LT,XT,XT1,XT2,XT3,YT,YT1,YT2,YT3
************************************************************************
DIMENSION MARK(500),LINE(500,2),KISSBAK(4000,3),ID(3),LPOINT(5500)
*
ENTRY INPUT
WRITE(*,*)'HOW MANY POINTS TO INSERT?'
READ(*,*)NPMAX
OPEN(1,FILE='INPUT.DAT')
READ(1,*)NPP
DO 600 I=1,NPP
READ(1,*)APOINT(I,1),APOINT(I,2)
AINP(I,1)=APOINT(I,1)
AINP(I,2)=APOINT(I,2)
600 CONTINUE
M=NPP
READ(1,*)N
DO 700 I=1,N
READ(1,*)NDOMAIN(I)
700 CONTINUE
READ(1,*) NBOUND
CLOSE(1)
RETURN
*
ENTRY MODIFY
KBAK=1
800 CONTINUE
IF(KBAK.EQ.NS) GOTO 1000
DO 900 J=KBAK,NS
NAA=KTSTACK(J,1)
NBB=KTSTACK(J,2)
NCC=KTSTACK(J,3)
X1=APOINT(NAA,1)
X2=APOINT(NBB,1)
X3=APOINT(NCC,1)
Y1=APOINT(NAA,2)
Y2=APOINT(NBB,2)
Y3=APOINT(NCC,2)
XSCAN=(X1+X2+X3)/3.
YSCAN=(Y1+Y2+Y3)/3.
CALL SCAN
IF(.NOT.INDOMAIN) THEN
DO 1200 JJ=J,NS-1
DO 1200 JM=1,3
1200 KTSTACK(JJ,JM)=KTSTACK(JJ+1,JM)
KBAK=J
NS=NS-1
GOTO 800
END IF
900 CONTINUE
1000 CONTINUE
RETURN
*
ENTRY BACKGROUND
NB=0.
DO 1300 K=1,N
DO 1300 I=1,NDOMAIN(K)
NB=NB+1
R1=(APOINT(NB,1)-APOINT(NB+1,1))**2.
R2=(APOINT(NB,2)-APOINT(NB+1,2))**2.
IF(I.EQ.NDOMAIN(K)) THEN
R1=(APOINT(NB,1)-APOINT(NB-NDOMAIN(K)+1,1))**2.
R2=(APOINT(NB,2)-APOINT(NB-NDOMAIN(K)+1,2))**2.
END IF
RB(NB)=SQRT(R1+R2)*1.732/2.
1300 CONTINUE
DO 1400 I=1,NS
DO 1400 JM=1,3
KISSBAK(I,JM)=KTSTACK(I,JM)
1400 CONTINUE
KISS=NS
RETURN
*
ENTRY SORTING
DO 1500 I=NS-NADD+1,NS
NCIRCLE=I
CALL CIRCLE
XT=X0
YT=Y0
DO 1600 J=1,KISS
N1=KISSBAK(J,1)
N2=KISSBAK(J,2)
N3=KISSBAK(J,3)
XT1=APOINT(N1,1)
YT1=APOINT(N1,2)
XT2=APOINT(N2,1)
YT2=APOINT(N2,2)
XT3=APOINT(N3,1)
YT3=APOINT(N3,2)
CALL INTRIANGLE
IF(LT) THEN
AINP1=1./(SQRT((XT1-X0)**2.+(YT1-Y0)**2.)+1.E-30)
AINP2=1./(SQRT((XT2-X0)**2.+(YT2-Y0)**2.)+1.E-30)
AINP3=1./(SQRT((XT3-X0)**2.+(YT3-Y0)**2.)+1.E-30)
RCXY=RB(N1)*AINP1+RB(N2)*AINP2+RB(N3)*AINP3
RCXY=RCXY/(AINP1+AINP2+AINP3+1.E-30)
GOTO 1700
END IF
1600 CONTINUE
1700 B(I)=R/(RCXY+1.E-30)
1500 CONTINUE
C**************INSERT ADD TRIANGLES IN THE STACK**********
DO 1800 K=1,NADD
KK=NS-NADD+K
DO 1800 I=1,KK-1
KKK=KK-I+1
AMIN=B(KKK)
BMIN=B(KKK-1)
IF(AMIN.GT.BMIN) THEN
DO JM=1,3
ID(JM)=KTSTACK(KKK,JM)
KTSTACK(KKK,JM)=KTSTACK(KKK-1,JM)
KTSTACK(KKK-1,JM)=ID(JM)
END DO
D=B(KKK)
B(KKK)=B(KKK-1)
B(KKK-1)=D
END IF
1800 CONTINUE
RETURN
*
ENTRY CIRCLE
N1=KTSTACK(NCIRCLE,1)
N2=KTSTACK(NCIRCLE,2)
N3=KTSTACK(NCIRCLE,3)
X1=APOINT(N1,1)
Y1=APOINT(N1,2)
X2=APOINT(N2,1)
Y2=APOINT(N2,2)
X3=APOINT(N3,1)
Y3=APOINT(N3,2)
XA=0.5*(X1+X2)
YA=0.5*(Y1+Y2)
XB=0.5*(X1+X3)
YB=0.5*(Y1+Y3)
IF(Y1.EQ.Y2) GOTO 1900
IF(Y1.EQ.Y3) GOTO 2000
AKA=-1.*(X1-X2)/(Y1-Y2)
AKB=-1.*(X1-X3)/(Y1-Y3)
X0=(YB-YA+AKA*XA-AKB*XB)/(AKA-AKB)
Y0=YA+AKA*(X0-XA)
R=SQRT((X0-X1)**2.+(Y0-Y1)**2.)
RETURN
1900 CONTINUE
X0=XA
Y0=YB-(X1-X3)*(X0-XB)/(Y1-Y3)
R=SQRT((X0-X1)**2.+(Y0-Y1)**2.)
RETURN
2000 CONTINUE
X0=XB
Y0=YA-(X1-X2)*(X0-XA)/(Y1-Y2)
R=SQRT((X0-X1)**2.+(Y0-Y1)**2.)
RETURN
*
ENTRY INTRIANGLE
X123=(XT1+XT2+XT3)/3.
Y123=(YT1+YT2+YT3)/3.
Z1=(YT2-YT1)*X123-(XT2-XT1)*Y123+YT1*(XT2-XT1)-XT1*(YT2-YT1)
Z2=(YT2-YT1)*XT-(XT2-XT1)*YT+YT1*(XT2-XT1)-XT1*(YT2-YT1)
Z12=Z1*Z2
Z2=(YT3-YT2)*X123-(XT3-XT2)*Y123+YT2*(XT3-XT2)-XT2*(YT3-YT2)
Z3=(YT3-YT2)*XT-(XT3-XT2)*YT+YT2*(XT3-XT2)-XT2*(YT3-YT2)
Z23=Z2*Z3
Z1=(YT3-YT1)*X123-(XT3-XT1)*Y123+YT1*(XT3-XT1)-XT1*(YT3-YT1)
Z3=(YT3-YT1)*XT-(XT3-XT1)*YT+YT1*(XT3-XT1)-XT1*(YT3-YT1)
Z13=Z1*Z3
LT=.FALSE.
IF(Z12.GE.0..AND.Z23.GE.0..AND.Z13.GE.0.) LT=.TRUE.
RETURN
*
ENTRY SCAN
NP=0
DO 2100 K=1,N
DO 2100 I=1,NDOMAIN(K)
NP=NP+1
LINE(NP,1)=NP
LINE(NP,2)=NP+1
IF(I.EQ.NDOMAIN(K)) LINE(NP,2)=NP-NDOMAIN(K)+1
2100 CONTINUE
NMARK=0
DO 2200 I=1,NP
IF(XSCAN.LT.AINP(LINE(I,1),1).AND.XSCAN.GE.AINP(LINE(I,2),1)) THEN
NMARK=NMARK+1
MARK(NMARK)=I
ELSE
IF(XSCAN.GE.AINP(LINE(I,1),1).AND.XSCAN.LT.AINP(LINE(I,2),1)) THEN
NMARK=NMARK+1
MARK(NMARK)=I
END IF
END IF
2200 CONTINUE
MM=NMARK
2300 IF(MM.GT.0) THEN
J=MM-1
MM=0
DO 2400 I=1,J
AMIN=0.5*(AINP(LINE(MARK(I),1),2)+AINP(LINE(MARK(I),2),2))
BMIN=0.5*(AINP(LINE(MARK(I+1),1),2)+AINP(LINE(MARK(I+1),2),2))
IF(AMIN.GT.BMIN) THEN
D=MARK(I)
MARK(I)=MARK(I+1)
MARK(I+1)=D
MM=I
END IF
2400 CONTINUE
GOTO 2300
END IF
INDOMAIN=.FALSE.
DO 2500 I=1,NMARK,2
X1=AINP(LINE(MARK(I),1),1)
Y1=AINP(LINE(MARK(I),1),2)
X2=AINP(LINE(MARK(I),2),1)
Y2=AINP(LINE(MARK(I),2),2)
YDOWN=Y1+(Y1-Y2)/(X1-X2)*(XSCAN-X1)
X1=AINP(LINE(MARK(I+1),1),1)
Y1=AINP(LINE(MARK(I+1),1),2)
X2=AINP(LINE(MARK(I+1),2),1)
Y2=AINP(LINE(MARK(I+1),2),2)
YUP=Y1+(Y1-Y2)/(X1-X2)*(XSCAN-X1)
IF(YSCAN.GT.YDOWN.AND.YSCAN.LT.YUP) THEN
INDOMAIN=.TRUE.
END IF
2500 CONTINUE
RETURN
*
ENTRY OUTPUT1
NPP=NPP-4
OPEN(2,FILE='DD.DAT')
WRITE(2,*)NPP
DO 2600 I=1,NPP
WRITE(2,*)APOINT(I,1),APOINT(I,2)
2600 CONTINUE
WRITE(2,*)NS
DO 2700 I=1,NS
WRITE(2,*)KTSTACK(I,1),KTSTACK(I,2),KTSTACK(I,3)
2700 CONTINUE
DO 2800 I=1,NS
WRITE(2,*) I,B(I)
2800 CONTINUE
CLOSE(2)
RETURN
*
ENTRY OUTPUT2
WRITE(*,*)'PRINTING THE DATA!'
C WANG GE GUAN HUA
DO 555 NNN=1,3
DO I=NBOUND+1,NPP
NADD=0
DO J=1,NS
NO1=I
NUM1=KTSTACK(J,1)
NUM2=KTSTACK(J,2)
NUM3=KTSTACK(J,3)
IF (NO1.EQ.NUM1) THEN
NADD=NADD+1
LPOINT(NADD)=NUM2
NADD=NADD+1
LPOINT(NADD)=NUM3
END IF
IF (NO1.EQ.NUM2) THEN
NADD=NADD+1
LPOINT(NADD)=NUM1
NADD=NADD+1
LPOINT(NADD)=NUM3
END IF
IF (NO1.EQ.NUM3) THEN
NADD=NADD+1
LPOINT(NADD)=NUM1
NADD=NADD+1
LPOINT(NADD)=NUM2
END IF
END DO
SUM1=0
SUM2=0
DO K=1,NADD
SUM1=SUM1+APOINT(LPOINT(K),1)
SUM2=SUM2+APOINT(LPOINT(K),2)
END DO
APOINT(I,1)=SUM1/NADD
APOINT(I,2)=SUM2/NADD
END DO
555 CONTINUE
DO 222 NNN=1,1
DO I=NPP,NBOUND+1,-1
NADD=0
DO J=1,NS
NO1=I
NUM1=KTSTACK(J,1)
NUM2=KTSTACK(J,2)
NUM3=KTSTACK(J,3)
IF (NO1.EQ.NUM1) THEN
NADD=NADD+1
LPOINT(NADD)=NUM2
NADD=NADD+1
LPOINT(NADD)=NUM3
END IF
IF (NO1.EQ.NUM2) THEN
NADD=NADD+1
LPOINT(NADD)=NUM1
NADD=NADD+1
LPOINT(NADD)=NUM3
END IF
IF (NO1.EQ.NUM3) THEN
NADD=NADD+1
LPOINT(NADD)=NUM1
NADD=NADD+1
LPOINT(NADD)=NUM2
END IF
END DO
SUM1=0
SUM2=0
DO K=1,NADD
SUM1=SUM1+APOINT(LPOINT(K),1)
SUM2=SUM2+APOINT(LPOINT(K),2)
END DO
APOINT(I,1)=SUM1/NADD
APOINT(I,2)=SUM2/NADD
END DO
222 CONTINUE
OPEN(10,FILE='GRID.DAT')
WRITE(10,*)NPP
WRITE(10,*)NBOUND
DO 2900 I=1,NPP
WRITE(10,*)APOINT(I,1),APOINT(I,2)
2900 CONTINUE
WRITE(10,*)NS
DO 3000 I=1,NS
WRITE(10,*)KTSTACK(I,1),KTSTACK(I,2),KTSTACK(I,3)
3000 CONTINUE
CLOSE(10)
RETURN
END
2008DELAUNEY.FOR
************************************************************************
*----------------------------MAIN PROGRAM------------------------------*
************************************************************************
CALL INPUT
CALL INITIAL
CALL NEWPOINT
CALL TECPLOT
END
************************************************************************
SUBROUTINE SETUP
************************************************************************
LOGICAL LI,LT,LL1,LL2,INDOMAIN
COMMON/TT/KTSTACK(5500,3),APOINT(5500,2),AINP(500,2),B(5500),
&NDOMAIN(400),RB(500)
COMMON/NN/NS,NPP,NADD,NBOUND,NPMAX,NCIRCLE,X0,Y0,R,M,N,
&INSERTNO,XSCAN,YSCAN,INDOMAIN,LT,XT,XT1,XT2,XT3,YT,YT1,YT2,YT3
************************************************************************
ENTRY INITIAL
NS=2
KTSTACK(1,1)=NPP-3
KTSTACK(1,2)=NPP-2
KTSTACK(1,3)=NPP-1
KTSTACK(2,1)=NPP-1
KTSTACK(2,2)=NPP
KTSTACK(2,3)=NPP-3
DO 100 INSERTNO=1,NPP-4
WRITE(*,*)'INSERTNO',INSERTNO
I=INSERTNO
CALL INSERT
100 CONTINUE
CALL MODIFY
CALL BACKGROUND
NADD=NS
CALL SORTING
CALL OUTPUT
RETURN
ENTRY NEWPOINT
DO 200 KINSERT=1,NPMAX
DO 300 NVE=1,NS
NCIRCLE=NVE
CALL CIRCLE
DO 400 J=1,NPP
IF(ABS(X0-APOINT(J,1)).LT.1E-20.AND.ABS(Y0-APOINT(J,2)).LT.1E-20)
&THEN
GOTO 300
END IF
400 CONTINUE
XSCAN=X0
YSCAN=Y0
CALL SCAN
IF(INDOMAIN) THEN
GOTO 500
ELSE
WRITE(*,*)'MISSHAPEN TRIANGLE WARNING:
& PLEASE INCREASE BOUNDARY POINTS!'
GOTO 300
END IF
300 CONTINUE
500 CONTINUE
NPP=NPP+1
WRITE(*,*)'I=',NPP,NS
WRITE(*,*)'DELTA=',B(NVE),B(NS),B(1)
APOINT(NPP,1)=X0
APOINT(NPP,2)=Y0
INSERTNO=NPP
CALL INSERT
CALL SORTING
IF(MOD(KINSERT,400).EQ.0.OR.KINSERT.EQ.NPMAX) THEN
CALL OUTPUT
ENDIF
200 CONTINUE
RETURN
END
************************************************************************
SUBROUTINE INSERT
************************************************************************
LOGICAL LI,LT,LL1,LL2,INDOMAIN
COMMON/TT/KTSTACK(5500,3),APOINT(5500,2),AINP(500,2),B(5500),
&NDOMAIN(400),RB(500)
COMMON/NN/NS,NPP,NADD,NBOUND,NPMAX,NCIRCLE,X0,Y0,R,M,N,
&INSERTNO,XSCAN,YSCAN,INDOMAIN,LT,XT,XT1,XT2,XT3,YT,YT1,YT2,YT3
************************************************************************
DIMENSION MPOINT(5500)
NPOINT=0
XX=APOINT(INSERTNO,1)
YY=APOINT(INSERTNO,2)
XT=XX
YT=YY
DO 3100 J=1,NS
N1=KTSTACK(J,1)
N2=KTSTACK(J,2)
N3=KTSTACK(J,3)
XT1=APOINT(N1,1)
YT1=APOINT(N1,2)
XT2=APOINT(N2,1)
YT2=APOINT(N2,2)
XT3=APOINT(N3,1)
YT3=APOINT(N3,2)
CALL INTRIANGLE
IF(LT) THEN
MPOINT(1)=N1
MPOINT(2)=N2
MPOINT(3)=N3
NPOINT=3
DO 3200 JJ=J,NS-1
DO 3200 JM=1,3
KTSTACK(JJ,JM)=KTSTACK(JJ+1,JM)
3200 B(JJ)=B(JJ+1)
NS=NS-1
GOTO 3300
END IF
3100 CONTINUE
3300 CONTINUE
KBAK=1
3400 IF(KBAK.EQ.NPOINT) GOTO 3500
DO 3600 JK=KBAK,NPOINT
L1=MPOINT(JK)
L2=MPOINT(JK+1)
IF(JK.EQ.NPOINT) L2=MPOINT(1)
NSS=NS
IF (NSS.EQ.0) GOTO 3500
DO 3700 J=1,NSS
M1=KTSTACK(J,1)
M2=KTSTACK(J,2)
M3=KTSTACK(J,3)
LL1=.FALSE.
LL2=.FALSE.
IF(L1.EQ.M1.OR.L1.EQ.M2.OR.L1.EQ.M3) LL1=.TRUE.
IF(L2.EQ.M1.OR.L2.EQ.M2.OR.L2.EQ.M3) LL2=.TRUE.
IF(LL1.AND.LL2) THEN
IF(M1.NE.L1.AND.M1.NE.L2) MNEW=M1
IF(M2.NE.L1.AND.M2.NE.L2) MNEW=M2
IF(M3.NE.L1.AND.M3.NE.L2) MNEW=M3
NCIRCLE=J
CALL CIRCLE
LI=.FALSE.
RR=(XX-X0)**2.+(YY-Y0)**2.
IF(RR.LE.(R**2.)) LI=.TRUE.
IF(LI) THEN
JK0=JK+1
DO 3800 K1=NPOINT,JK0,-1
MPOINT(K1+1)=MPOINT(K1)
3800 CONTINUE
MPOINT(JK0)=MNEW
NPOINT=NPOINT+1
KBAK=JK
DO 3900 JJ=J,NS-1
DO 3900 JM=1,3
KTSTACK(JJ,JM)=KTSTACK(JJ+1,JM)
3900 B(JJ)=B(JJ+1)
NS=NS-1
GOTO 3400
END IF
END IF
3700 CONTINUE
3600 CONTINUE
3500 CONTINUE
NADD=0
DO 4000 J=1,NPOINT
NS=NS+1
NADD=NADD+1
KTSTACK(NS,1)=MPOINT(J)
KTSTACK(NS,2)=MPOINT(J+1)
KTSTACK(NS,3)=INSERTNO
IF(J.EQ.NPOINT) THEN
KTSTACK(NS,2)=MPOINT(1)
END IF
4000 CONTINUE
RETURN
END
************************************************************************
RECURSIVE SUBROUTINE SUPPLY
************************************************************************
LOGICAL LI,LT,LL1,LL2,INDOMAIN
COMMON/TT/KTSTACK(5500,3),APOINT(5500,2),AINP(500,2),B(5500),
&NDOMAIN(400),RB(500)
COMMON/NN/NS,NPP,NADD,NBOUND,NPMAX,NCIRCLE,X0,Y0,R,M,N,
&INSERTNO,XSCAN,YSCAN,INDOMAIN,LT,XT,XT1,XT2,XT3,YT,YT1,YT2,YT3
************************************************************************
DIMENSION MARK(500),LINE(500,2),KISSBAK(4000,3),ID(3),LPOINT(5500)
ENTRY INPUT
WRITE(*,*)'HOW MANY POINTS TO INSERT?'
READ(*,*)NPMAX
OPEN(1,FILE='INPUT.TXT')
READ(1,*)NPP
DO 600 I=1,NPP
READ(1,*)APOINT(I,1),APOINT(I,2)
AINP(I,1)=APOINT(I,1)
AINP(I,2)=APOINT(I,2)
600 CONTINUE
M=NPP
READ(1,*)N
DO 700 I=1,N
READ(1,*)NDOMAIN(I)
700 CONTINUE
READ(1,*) NBOUND
CLOSE(1)
RETURN
ENTRY MODIFY
C=========THIS PART IS MODIFIED=======
KBAK=1
800 CONTINUE
DO 900 J=KBAK,NS
NAA=KTSTACK(J,1)
NBB=KTSTACK(J,2)
NCC=KTSTACK(J,3)
X1=APOINT(NAA,1)
X2=APOINT(NBB,1)
X3=APOINT(NCC,1)
Y1=APOINT(NAA,2)
Y2=APOINT(NBB,2)
Y3=APOINT(NCC,2)
XSCAN=(X1+X2+X3)/3.
YSCAN=(Y1+Y2+Y3)/3.
CALL SCAN
IF(.NOT.INDOMAIN) THEN
DO 1200 JJ=J,NS-1
DO 1200 JM=1,3
1200 KTSTACK(JJ,JM)=KTSTACK(JJ+1,JM)
KBAK=J
NS=NS-1
GOTO 800
END IF
900 CONTINUE
C=========THIS PART IS MODIFIED=======
RETURN
ENTRY BACKGROUND
NB=0.
DO 1300 K=1,N
C=========THIS PART IS MODIFIED=======
NB=NB+1
R1=0
R1=R1+(APOINT(NB,1)-APOINT(NB+1,1))**2.
R1=R1+(APOINT(NB,2)-APOINT(NB+1,2))**2.
R2=0
R2=R2+(APOINT(NB,1)-APOINT(NB+NDOMAIN(K)-1,1))**2.
R2=R2+(APOINT(NB,2)-APOINT(NB+NDOMAIN(K)-1,2))**2.
RB(NB)=0.5*(SQRT(R1)+SQRT(R2))*1.732/2.
DO 1250 I=2,NDOMAIN(K)-1
NB=NB+1
R1=0
R1=R1+(APOINT(NB,1)-APOINT(NB+1,1))**2.
R1=R1+(APOINT(NB,2)-APOINT(NB+1,2))**2.
R2=0
R2=R2+(APOINT(NB,1)-APOINT(NB-1,1))**2.
R2=R2+(APOINT(NB,2)-APOINT(NB-1,2))**2.
RB(NB)=0.5*(SQRT(R1)+SQRT(R2))*1.732/2.
1250 CONTINUE
NB=NB+1
R1=0
R1=R1+(APOINT(NB,1)-APOINT(NB-NDOMAIN(K)+1,1))**2.
R1=R1+(APOINT(NB,2)-APOINT(NB-NDOMAIN(K)+1,2))**2.
R2=0
R2=R2+(APOINT(NB,1)-APOINT(NB-1,1))**2.
R2=R2+(APOINT(NB,2)-APOINT(NB-1,2))**2.
RB(NB)=0.5*(SQRT(R1)+SQRT(R2))*1.732/2.
C=========THIS PART IS MODIFIED=======
1300 CONTINUE
DO 1400 I=1,NS
DO 1400 JM=1,3
KISSBAK(I,JM)=KTSTACK(I,JM)
1400 CONTINUE
KISS=NS
RETURN
ENTRY SORTING
DO 1500 I=NS-NADD+1,NS
NCIRCLE=I
CALL CIRCLE
XT=X0
YT=Y0
DO 1600 J=1,KISS
N1=KISSBAK(J,1)
N2=KISSBAK(J,2)
N3=KISSBAK(J,3)
XT1=APOINT(N1,1)
YT1=APOINT(N1,2)
XT2=APOINT(N2,1)
YT2=APOINT(N2,2)
XT3=APOINT(N3,1)
YT3=APOINT(N3,2)
CALL INTRIANGLE
IF(LT) THEN
AINP1=1./(SQRT((XT1-X0)**2.+(YT1-Y0)**2.)+1.E-30)
AINP2=1./(SQRT((XT2-X0)**2.+(YT2-Y0)**2.)+1.E-30)
AINP3=1./(SQRT((XT3-X0)**2.+(YT3-Y0)**2.)+1.E-30)
RCXY=RB(N1)*AINP1+RB(N2)*AINP2+RB(N3)*AINP3
RCXY=RCXY/(AINP1+AINP2+AINP3+1.E-30)
GOTO 1700
END IF
1600 CONTINUE
1700 B(I)=R/(RCXY+1.E-30)
1500 CONTINUE
C**************INSERT ADD TRIANGLES IN THE STACK**********
DO 1800 K=1,NADD
KK=NS-NADD+K
DO 1800 I=1,KK-1
KKK=KK-I+1
AMIN=B(KKK)
BMIN=B(KKK-1)
IF(AMIN.GT.BMIN) THEN
DO JM=1,3
ID(JM)=KTSTACK(KKK,JM)
KTSTACK(KKK,JM)=KTSTACK(KKK-1,JM)
KTSTACK(KKK-1,JM)=ID(JM)
END DO
D=B(KKK)
B(KKK)=B(KKK-1)
B(KKK-1)=D
END IF
1800 CONTINUE
RETURN
ENTRY CIRCLE
N1=KTSTACK(NCIRCLE,1)
N2=KTSTACK(NCIRCLE,2)
N3=KTSTACK(NCIRCLE,3)
X1=APOINT(N1,1)
Y1=APOINT(N1,2)
X2=APOINT(N2,1)
Y2=APOINT(N2,2)
X3=APOINT(N3,1)
Y3=APOINT(N3,2)
XA=0.5*(X1+X2)
YA=0.5*(Y1+Y2)
XB=0.5*(X1+X3)
YB=0.5*(Y1+Y3)
IF(Y1.EQ.Y2) GOTO 1900
IF(Y1.EQ.Y3) GOTO 2000
AKA=-1.*(X1-X2)/(Y1-Y2)
AKB=-1.*(X1-X3)/(Y1-Y3)
X0=(YB-YA+AKA*XA-AKB*XB)/(AKA-AKB)
Y0=YA+AKA*(X0-XA)
R=SQRT((X0-X1)**2.+(Y0-Y1)**2.)
RETURN
1900 CONTINUE
X0=XA
Y0=YB-(X1-X3)*(X0-XB)/(Y1-Y3)
R=SQRT((X0-X1)**2.+(Y0-Y1)**2.)
RETURN
2000 CONTINUE
X0=XB
Y0=YA-(X1-X2)*(X0-XA)/(Y1-Y2)
R=SQRT((X0-X1)**2.+(Y0-Y1)**2.)
RETURN
ENTRY INTRIANGLE
X123=(XT1+XT2+XT3)/3.
Y123=(YT1+YT2+YT3)/3.
Z1=(YT2-YT1)*X123-(XT2-XT1)*Y123+YT1*(XT2-XT1)-XT1*(YT2-YT1)
Z2=(YT2-YT1)*XT-(XT2-XT1)*YT+YT1*(XT2-XT1)-XT1*(YT2-YT1)
Z12=Z1*Z2
Z2=(YT3-YT2)*X123-(XT3-XT2)*Y123+YT2*(XT3-XT2)-XT2*(YT3-YT2)
Z3=(YT3-YT2)*XT-(XT3-XT2)*YT+YT2*(XT3-XT2)-XT2*(YT3-YT2)
Z23=Z2*Z3
Z1=(YT3-YT1)*X123-(XT3-XT1)*Y123+YT1*(XT3-XT1)-XT1*(YT3-YT1)
Z3=(YT3-YT1)*XT-(XT3-XT1)*YT+YT1*(XT3-XT1)-XT1*(YT3-YT1)
Z13=Z1*Z3
LT=.FALSE.
IF(Z12.GE.0..AND.Z23.GE.0..AND.Z13.GE.0.) LT=.TRUE.
RETURN
ENTRY SCAN
NP=0
DO 2100 K=1,N
DO 2100 I=1,NDOMAIN(K)
NP=NP+1
LINE(NP,1)=NP
LINE(NP,2)=NP+1
IF(I.EQ.NDOMAIN(K)) LINE(NP,2)=NP-NDOMAIN(K)+1
2100 CONTINUE
C=========THIS PART IS MODIFIED=======
INDOMAIN=.FALSE.
NMARK=0
DO 2200 I=1,NP
IF(XSCAN.EQ.AINP(LINE(I,1),1).AND.XSCAN.EQ.AINP(LINE(I,2),1)) THEN
LL1=YSCAN.LE.AINP(LINE(I,1),2).AND.YSCAN.GE.AINP(LINE(I,2),2)
LL2=YSCAN.GE.AINP(LINE(I,1),2).AND.YSCAN.LE.AINP(LINE(I,2),2)
IF(LL1.OR.LL2) RETURN
END IF
IF(XSCAN.LT.AINP(LINE(I,1),1).AND.XSCAN.GE.AINP(LINE(I,2),1)) THEN
NMARK=NMARK+1
MARK(NMARK)=I
ELSE
IF(XSCAN.GE.AINP(LINE(I,1),1).AND.XSCAN.LT.AINP(LINE(I,2),1)) THEN
NMARK=NMARK+1
MARK(NMARK)=I
END IF
END IF
2200 CONTINUE
MM=NMARK
DO 2500 I=1,NMARK
X1=AINP(LINE(MARK(I),1),1)
Y1=AINP(LINE(MARK(I),1),2)
X2=AINP(LINE(MARK(I),2),1)
Y2=AINP(LINE(MARK(I),2),2)
YLINE=Y1+(Y1-Y2)/(X1-X2)*(XSCAN-X1)
IF(YLINE.EQ.YSCAN) THEN
INDOMAIN=.FALSE.
RETURN
ELSE IF(YLINE.GT.YSCAN) THEN
INDOMAIN=.NOT.INDOMAIN
END IF
2500 CONTINUE
C=========THIS PART IS MODIFIED=======
RETURN
ENTRY OUTPUT
OPEN(10,FILE='GRID.DAT')
WRITE(10,*)NPP
WRITE(10,*)NBOUND
DO 2900 I=1,NPP
WRITE(10,*)APOINT(I,1),APOINT(I,2)
2900 CONTINUE
WRITE(10,*)NS
DO 3000 I=1,NS
WRITE(10,*)KTSTACK(I,1),KTSTACK(I,2),KTSTACK(I,3)
3000 CONTINUE
CLOSE(10)
RETURN
ENTRY TECPLOT
OPEN(8,FILE='TECPLOT.DAT')
WRITE(8,*) 'TITLE = " Delaunay 2D Grid " '
WRITE(8,*) 'VARIABLES = "X" "Y" '
WRITE(8,*) 'ZONE T="FE-Boundary: Src= Zone 1" '
WRITE(8,*) 'N=',NPP,', E=',NS,', F=FEPOINT ET=Triangle'
WRITE(8,*) 'DT=(SINGLE SINGLE SINGLE)'
DO 3600 I=1,NPP
WRITE(8,*)APOINT(I,1),APOINT(I,2)
3600 CONTINUE
DO 3700 I=1,NS
WRITE(8,*)KTSTACK(I,1),KTSTACK(I,2),KTSTACK(I,3)
3700 CONTINUE
CLOSE(8)
RETURN
END
1723
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
