program main
C*********************************************************************72
C
cc MAIN is the main program for TRIPACK_PRB
c
c Discussion:
c
c TRIPACK_PRB tests the TRIPACK library.
c
C TRITEST: Portable Test Driver for TRIPACK
C 07/02/98
C
C
C This driver tests software package TRIPACK for construc-
C ting a constrained Delaunay triangulation of a set of
C points in the plane. All modules other than TRMSHR are
C tested unless an error is encountered, in which case the
C program terminates immediately.
C
C By default, tests are performed on a simple data set
C consisting of 12 nodes whose convex hull covers the unit
C square. The data set includes a single constraint region
C consisting of four nodes forming a smaller square at the
C center of the unit square. However, by enabling the READ
C statements below (C# in the first two columns), testing
C may be performed on an arbitrary set of up to NMAX nodes
C with up to NCMAX constraint curves. (Refer to the
C PARAMETER statements below.) A data set consists of the
C following sequence of records:
C
C N = Number of nodes (format I4) -- 3 to NMAX.
C NCC = Number of constraint curves (format I4) -- 0 to
C NCMAX.
C (LCC(I), I = 1,NCC) = Indexes of the first node in each
C constraint curve (format I4). 1 .LE. LCC(1)
C and, for I .GT. 1, LCC(I-1) + 3 .LE. LCC(I)
C .LE. N-2. (Each constraint curve has at least
C three nodes.)
C (X(I),Y(I), I = 1,N) = Nodal coordinates with non-
C constraint nodes followed by the NCC
C sequences of constraint nodes (format
C 2F13.8).
C
C The I/O units may be changed by altering LIN (input) and
C LOUT (output) in the DATA statement below.
C
C This driver must be linked to TRIPACK.
C
CHARACTER*80 TITLE
INTEGER IER, IO1, IO2, K, KSUM, LIN, LNEW, LOUT, LPLT,
. LW, LWK, N, N0, N6, NA, NB, NCC, NCMAX, NMAX,
. NN, NROW, NT, NTMX
INTEGER NEARND
LOGICAL NUMBR, PRNTX
REAL A, ARMAX, DSQ, PLTSIZ, TOL, WX1, WX2, WY1, WY2
REAL AREAP
C
PARAMETER (NMAX=100, NCMAX=5, NTMX=2*NMAX, N6=6*NMAX,
. LWK=2*NMAX, NROW=9)
C
C Array storage:
C
INTEGER LCC(NCMAX), LCT(NCMAX), LEND(NMAX), LIST(N6),
. LPTR(N6), LTRI(NROW,NTMX), NODES(LWK)
REAL DS(NMAX), X(NMAX), Y(NMAX)
C
C Tolerance for TRMTST and NEARND: upper bound on squared
C distances.
C
DATA TOL/1.E-2/
C
C Plot size for the triangulation plot.
C
DATA PLTSIZ/7.5/
C
C Default data set:
C
DATA N/12/, NCC/1/, LCC(1)/9/
DATA X(1), X(2), X(3), X(4), X(5), X(6), X(7)
. / 0., 1., .5, .15, .85, .5, 0./,
. X(8), X(9), X(10), X(11), X(12)
. / 1., .35, .65, .65, .35/,
. Y(1), Y(2), Y(3), Y(4), Y(5), Y(6), Y(7)
. / 0., 0., .15, .5, .5, .85, 1./,
. Y(8), Y(9), Y(10), Y(11), Y(12)
. / 1., .35, .35, .65, .65/
C
C Logical unit numbers for I/O:
C
DATA LIN/1/, LOUT/6/, LPLT/3/
C OPEN (LOUT,FILE='RES')
OPEN (LPLT,FILE='tripack_prb.eps')
C
C Store a plot title. It must be enclosed in parentheses.
C
TITLE = '(Triangulation created by TRITEST)'
C
C *** Read triangulation parameters -- N, NCC, LCC, X, Y.
C
C# OPEN (LIN,FILE='tritest.dat',STATUS='OLD')
C# READ (LIN,100,ERR=30) N, NCC
IF (N .LT. 3 .OR. N .GT. NMAX .OR. NCC .LT. 0
. .OR. NCC .GT. NCMAX) GO TO 31
C# IF (NCC .GT. 0) READ (LIN,100,ERR=30)
C# . (LCC(K), K = 1,NCC)
C# READ (LIN,110,ERR=30) (X(K),Y(K), K = 1,N)
C#100 FORMAT (I4)
C#110 FORMAT (2F13.8)
C
C Print a heading.
C
WRITE (LOUT,400) N
C
C *** Create the Delaunay triangulation (TRMESH), and test
C for errors (refer to TRMTST below). NODES and DS are
C used as work space.
C
CALL TRMESH (N,X,Y, LIST,LPTR,LEND,LNEW,NODES,
. NODES(N+1),DS,IER)
IF (IER .EQ. -2) THEN
WRITE (LOUT,210)
STOP
ELSEIF (IER .EQ. -4) THEN
WRITE (LOUT,212)
STOP
ELSEIF (IER .GT. 0) THEN
WRITE (LOUT,214)
STOP
ENDIF
CALL TRMTST (N,X,Y,LIST,LPTR,LEND,LNEW,TOL,
. LOUT, ARMAX,IER)
WRITE (LOUT,410) ARMAX
IF (IER .GT. 0) STOP
C
C *** Add the constraint curves (ADDCST). Note that edges
C and triangles are not removed from constraint regions.
C ADDCST forces the inclusion of triangulation edges
C connecting the sequences of constraint nodes. If it
C is necessary to alter the triangulation, the empty
C circumcircle property is no longer satisfied.
C
LW = LWK
CALL ADDCST (NCC,LCC,N,X,Y, LW,NODES,LIST,LPTR,
. LEND, IER)
IF (IER .NE. 0) THEN
WRITE (LOUT,220) IER
STOP
ENDIF
IF (LW .EQ. 0) WRITE (LOUT,430)
C
C *** Test TRPRNT, TRLIST, and TRLPRT, and TRPLOT.
C
PRNTX = .TRUE.
CALL TRPRNT (NCC,LCC,N,X,Y,LIST,LPTR,LEND,LOUT,PRNTX)
CALL TRLIST (NCC,LCC,N,LIST,LPTR,LEND,NROW, NT,LTRI,
. LCT,IER)
CALL TRLPRT (NCC,LCT,N,X,Y,NROW,NT,LTRI,LOUT,PRNTX)
C
C Set the plot window [WX1,WX2] X [WY1,WY2] to the
C smallest rectangle that contains the nodes.
C NUMBR = TRUE iff nodal indexes are to be displayed.
C
WX1 = X(1)
WX2 = WX1
WY1 = Y(1)
WY2 = WY1
DO 1 K = 2,N
IF (X(K) .LT. WX1) WX1 = X(K)
IF (X(K) .GT. WX2) WX2 = X(K)
IF (Y(K) .LT. WY1) WY1 = Y(K)
IF (Y(K) .GT. WY2) WY2 = Y(K)
1 CONTINUE
NUMBR = N .LE. 200
CALL TRPLOT (LPLT,PLTSIZ,WX1,WX2,WY1,WY2,NCC,LCC,N,
. X,Y,LIST,LPTR,LEND,TITLE,NUMBR, IER)
IF (IER .EQ. 0) THEN
WRITE (LOUT,470)
ELSE
WRITE (LOUT,270) IER
ENDIF
C
C *** Test BNODES and AREAP.
C
CALL BNODES (N,LIST,LPTR,LEND, NODES,NB,NA,NT)
A = AREAP(X,Y,NB,NODES)
WRITE (LOUT,420) NB, NA, NT, A
C
C *** Test GETNP by ordering the nodes on distance from N0
C and verifying the ordering.
C
N0 = N/2
NODES(1) = N0
DS(1) = 0.
KSUM = N0
DO 2 K = 2,N
CALL GETNP (NCC,LCC,N,X,Y,LIST,LPTR,LEND,
. K, NODES,DS, IER)
IF (IER .NE. 0 .OR. DS(K) .LT. DS(K-1)) THEN
WRITE (LOUT,230)
STOP
ENDIF
KSUM = KSUM + NODES(K)
2 CONTINUE
C
C Test for all nodal indexes included in NODES.
C
IF (KSUM .NE. (N*(N+1))/2) THEN
WRITE (LOUT,230)
STOP
ENDIF
C
C *** Test NEARND by verifying that the nearest node to K is
C node K for K = 1 to N.
C
DO 3 K = 1,N
N0 = NEARND (X(K),Y(K),1,N,X,Y,LIST,LPTR,LEND, DSQ)
IF (N0 .NE. K .OR. DSQ .GT. TOL) THEN
WRITE (LOUT,240)
STOP
ENDIF
3 CONTINUE
C
C *** Test DELARC by removing a boundary arc if possible.
C The first two nodes define a boundary arc
C in the default data set.
C
IO1 = 1
IO2 = 2
CALL DELARC (N,IO1,IO2, LIST,LPTR,LEND,LNEW, IER)
IF (IER .EQ. 1 .OR. IER .EQ. 4) THEN
WRITE (LOUT,250) IER
STOP
ENDIF
IF (IER .NE. 0) WRITE (LOUT,440)
C
C Recreate the triangulation without constraints.
C
CALL TRMESH (N,X,Y, LIST,LPTR,LEND,LNEW,NODES,
. NODES(N+1),DS,IER)
NCC = 0
C
C *** Test DELNOD by removing nodes 4 to N (in reverse
C order).
C
IF (N .LE. 3) THEN
WRITE (LOUT,450)
ELSE
NN = N
4 LW = LWK/2
CALL DELNOD (NN,NCC, LCC,NN,X,Y,LIST,LPTR,LEND,
. LNEW,LW,NODES, IER)
IF (IER .NE. 0) THEN
WRITE (LOUT,260) IER
STOP
ENDIF
IF (NN .GT. 3) GO TO 4
ENDIF
C
C Successful test.
C
WRITE (LOUT,460)
STOP
C
C Error reading the data set.
C
C# 30 WRITE (*,200)
STOP
C
C Invalid value of N or NCC.
C
31 WRITE (*,205) N, NCC
STOP
C
C Error message formats:
C
C#200 FORMAT (//5X,'*** Input data set invalid ***'/)
205 FORMAT (//5X,'*** N or NCC is outside its valid ',
. 'range: N =',I5,', NCC = ',I4,' ***'/)
210 FORMAT (//5X,'*** Error in TRMESH: the first three ',
. 'nodes are collinear ***'/)
212 FORMAT (//5X,'*** Error in TRMESH: invalid ',
. 'triangulation ***'/)
214 FORMAT (//5X,'*** Error in TRMESH: duplicate nodes ',
. 'encountered ***'/)
220 FORMAT (//5X,'*** Error in ADDCST: IER = ',I1,
. ' ***'/)
230 FORMAT (//5X,'*** Error in GETNP ***'/)
240 FORMAT (//5X,'*** Error in NEARND ***'/)
250 FORMAT (//5X,'*** Error in DELARC: IER = ',I1,
. ' ***'/)
260 FORMAT (//5X,'*** Error in DELNOD: IER = ',I1,
. ' ***'/)
270 FORMAT (//5X,'*** Error in TRPLOT: IER = ',I1,
. ' ***'/)
C
C Informative message formats:
C
400 FORMAT (///1X,21X,'TRIPACK Test: N =',I5///)
410 FORMAT (5X,'Maximum triangle aspect ratio = ',E10.3//)
420 FORMAT (/5X,'Output from BNODES and AREAP'//
. 5X,'BNODES: ',I4,' boundary nodes, ',I5,
. ' edges, ',I5,' triangles'/5X,
. 'AREAP: area of convex hull = ',E15.8//)
430 FORMAT (5X,'Subroutine EDGE not tested:'/
. 5X,' No edges were swapped by ADDCST'//)
440 FORMAT (5X,'Subroutine DELARC not tested:'/
. 5X,' Nodes 1 and 2 do not form a ',
. 'removable boundary edge.'//)
450 FORMAT (5X,'Subroutine DELNOD not tested:'/
. 5X,' N cannot be reduced below 3'//)
460 FORMAT (5X,'No triangulation errors encountered.'/)
470 FORMAT (/5X,'A triangulation plot file was ',
. 'successfully created.'/)
END
SUBROUTINE TRMTST (N,X,Y,LIST,LPTR,LEND,LNEW,TOL,
. LUN, ARMAX,IER)
INTEGER N, LIST(*), LPTR(*), LEND(N), LNEW, LUN, IER
REAL X(N), Y(N), TOL, ARMAX
C
C***********************************************************
C
C From TRPTEST
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C (817) 565-2767
C 09/01/91
C
C This subroutine tests the validity of the data structure
C representing a Delaunay triangulation created by subrou-
C tine TRMESH. The following properties are tested:
C
C 1) Each interior node has at least three neighbors, and
C each boundary node has at least two neighbors.
C
C 2) abs(LIST(LP)) is a valid nodal index in the range
C 1 to N and LIST(LP) > 0 unless LP = LEND(K) for some
C nodal index K.
C
C 3) Each pointer LEND(K) for K = 1 to N and LPTR(LP) for
C LP = 1 to LNEW-1 is a valid LIST index in the range
C 1 to LNEW-1.
C
C 4) N .GE. NB .GE. 3, NT = 2*N-NB-2, and NA = 3*N-NB-3 =
C (LNEW-1)/2, where NB, NT, and NA are the numbers of
C boundary nodes, triangles, and arcs, respectively.
C
C 5) Each circumcircle defined by the vertices of a tri-
C angle contains no nodes in its interior. This prop-
C erty distinguishes a Delaunay triangulation from an
C arbitrary triangulation of the nodes.
C
C Note that no test is made for the property that a triangu-
C lation covers the convex hull of the nodes.
C
C On input:
C
C N = Number of nodes. N .GE. 3.
C
C X,Y = Arrays of length N containing the nodal
C coordinates.
C
C LIST,LPTR,LEND = Data structure containing the tri-
C angulation. Refer to subroutine
C TRMESH.
C
C TOL = Nonnegative tolerance to allow for floating-
C point errors in the circumcircle test. An
C error situation is defined as (R**2 - D**2)/
C R**2 > TOL, where R is the radius of a circum-
C circle and D is the distance from the
C circumcenter to the nearest node. A reason-
C able value for TOL is 10*EPS, where EPS is the
C machine precision. The test is effectively
C bypassed by making TOL large. If TOL < 0, the
C tolerance is taken to be 0.
C
C LUN = Logical unit number for printing error mes-
C sages. If LUN < 0 or LUN > 99, no messages
C are printed.
C
C Input parameters are not altered by this routine.
C
C On output:
C
C ARMAX = Maximum aspect ratio (radius of inscribed
C circle divided by circumradius) of a trian-
C gle in the triangulation unless IER > 0.
C
C IER = Error indicator:
C IER = -1 if one or more null triangles (area =
C 0) are present but no (other) errors
C were encountered. A null triangle is
C an error only if it occurs in the
C the interior.
C IER = 0 if no errors or null triangles were
C encountered.
C IER = 1 if a node has too few neighbors.
C IER = 2 if a LIST entry is outside its valid
C range.
C IER = 3 if a LPTR or LEND entry is outside its
C valid range.
C IER = 4 if the triangulation parameters (N,
C NB, NT, NA, and LNEW) are inconsistent
C (or N < 3 or LNEW is invalid).
C IER = 5 if a triangle contains a node interior
C to its circumcircle.
C
C Module required by TRMTST: CIRCUM
C
C Intrinsic function called by TRMTST: MAX, ABS
C
C***********************************************************
C
INTEGER K, LP, LPL, LPMAX, LPN, N1, N2, N3, NA, NB,
. NFAIL, NN, NNB, NT, NULL
LOGICAL RATIO, RITE
REAL AR, CR, CX, CY, RS, RTOL, SA
C
C Store local variables, test for errors in input, and
C initialize counts.
C
NN = N
LPMAX = LNEW - 1
RTOL = TOL
IF (RTOL .LT. 0.) RTOL = 0.
RATIO = .TRUE.
ARMAX = 0.
RITE = LUN .GE. 0 .AND. LUN .LE. 99
IF (NN .LT. 3) GO TO 14
NB = 0
NT = 0
NULL = 0
NFAIL = 0
C
C Loop on triangles (N1,N2,N3) such that N2 and N3 index
C adjacent neighbors of N1 and are both larger than N1
C (each triangle is associated with its smallest index).
C NNB is the neighbor count for N1.
C
DO 5 N1 = 1,NN
NNB = 0
LPL = LEND(N1)
IF (LPL .LT. 1 .OR. LPL .GT. LPMAX) THEN
LP = LPL
GO TO 13
ENDIF
LP = LPL
C
C Loop on neighbors of N1.
C
1 LP = LPTR(LP)
NNB = NNB + 1
IF (LP .LT. 1 .OR. LP .GT. LPMAX) GO TO 13
N2 = LIST(LP)
IF (N2 .LT. 0) THEN
IF (LP .NE. LPL) GO TO 12
IF (N2 .EQ. 0 .OR. -N2 .GT. NN) GO TO 12
NB = NB + 1
GO TO 4
ENDIF
IF (N2 .LT. 1 .OR. N2 .GT. NN) GO TO 12
LPN = LPTR(LP)
N3 = ABS(LIST(LPN))
IF (N2 .LT. N1 .OR. N3 .LT. N1) GO TO 4
NT = NT + 1
C
C Compute the coordinates of the circumcenter of
C (N1,N2,N3).
C
CALL CIRCUM (X(N1),Y(N1),X(N2),Y(N2),X(N3),Y(N3),
. RATIO, CX,CY,CR,SA,AR)
IF (SA .EQ. 0.) THEN
NULL = NULL + 1
GO TO 4
ENDIF
ARMAX = MAX(ARMAX,AR)
C
C Test for nodes within the circumcircle.
C
RS = CR*CR*(1.-RTOL)
DO 2 K = 1,NN
IF (K .EQ. N1 .OR. K .EQ. N2 .OR.
. K .EQ. N3) GO TO 2
IF ((CX-X(K))**2 + (CY-Y(K))**2 .LT. RS) GO TO 3
2 CONTINUE
GO TO 4
C
C Node K is interior to the circumcircle of (N1,N2,N3).
C
3 NFAIL = NFAIL + 1
C
C Bottom of loop on neighbors.
C
4 IF (LP .NE. LPL) GO TO 1
IF (NNB .LT. 2 .OR. (NNB .EQ. 2 .AND.
. LIST(LPL) .GT. 0)) GO TO 11
5 CONTINUE
C
C Test parameters for consistency and check for NFAIL = 0.
C
NA = LPMAX/2
IF (NB .LT. 3 .OR. NT .NE. 2*NN-NB-2 .OR.
. NA .NE. 3*NN-NB-3) GO TO 14
IF (NFAIL .NE. 0) GO TO 15
C
C No errors were encountered.
C
IER = 0
IF (NULL .EQ. 0) RETURN
IER = -1
IF (RITE) WRITE (LUN,100) NULL
100 FORMAT (//5X,'*** TRMTST -- ',I5,' NULL TRIANGLES ',
. 'ARE PRESENT'/19X,'(NULL TRIANGLES ',
. 'ON THE BOUNDARY ARE UNAVOIDABLE) ***'//)
RETURN
C
C Node N1 has fewer than three neighbors.
C
11 IER = 1
IF (RITE) WRITE (LUN,110) N1, NNB
110 FORMAT (//5X,'*** TRMTST -- NODE ',I5,
. ' HAS ONLY ',I5,' NEIGHBORS ***'/)
RETURN
C
C N2 = LIST(LP) is outside its valid range.
C
12 IER = 2
IF (RITE) WRITE (LUN,120) N2, LP, N1
120 FORMAT (//5X,'*** TRMTST -- LIST(LP) =',I5,
. ', FOR LP =',I5,','/19X,
.'IS NOT A VALID NEIGHBOR OF ',I5,' ***'/)
RETURN
C
C LIST pointer LP is outside its valid range.
C
13 IER = 3
IF (RITE) WRITE (LUN,130) LP, LNEW, N1
130 FORMAT (//5X,'*** TRMTST -- LP =',I5,' IS NOT IN THE',
. ' RANGE 1 TO LNEW-1 FOR LNEW = ',I5/
. 19X,'LP POINTS TO A NEIGHBOR OF ',I5,
. ' ***'/)
RETURN
C
C Inconsistent triangulation parameters encountered.
C
14 IER = 4
IF (RITE) WRITE (LUN,140) N, LNEW, NB, NT, NA
140 FORMAT (//5X,'*** TRMTST -- INCONSISTENT PARAMETERS',
. ' ***'/19X,'N = ',I5,' NODES',12X,'LNEW =',I5/
. 19X,'NB = ',I5,' BOUNDARY NODES'/
. 19X,'NT = ',I5,' TRIANGLES'/
. 19X,'NA = ',I5,' ARCS'/)
RETURN
C
C Circumcircle test failure.
C
15 IER = 5
IF (RITE) WRITE (LUN,150) NFAIL
150 FORMAT (//5X,'*** TRMTST -- ',I5,' CIRCUMCIRCLES ',
. 'CONTAIN NODES IN THEIR INTERIORS ***'/)
RETURN
END