SUBROUTINE ADDCST (NCC,LCC,N,X,Y, LWK,IWK,LIST,LPTR,
. LEND, IER)
INTEGER NCC, LCC(*), N, LWK, IWK(LWK), LIST(*),
. LPTR(*), LEND(N), IER
REAL X(N), Y(N)
C
C***********************************************************
C
C From TRIPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C renka@cs.unt.edu
C 11/12/94
C
C This subroutine provides for creation of a constrained
C Delaunay triangulation which, in some sense, covers an
C arbitrary connected region R rather than the convex hull
C of the nodes. This is achieved simply by forcing the
C presence of certain adjacencies (triangulation arcs) cor-
C responding to constraint curves. The union of triangles
C coincides with the convex hull of the nodes, but triangles
C in R can be distinguished from those outside of R. The
C only modification required to generalize the definition of
C the Delaunay triangulation is replacement of property 5
C (refer to TRMESH) by the following:
C
C 5') If a node is contained in the interior of the cir-
C cumcircle of a triangle, then every interior point
C of the triangle is separated from the node by a
C constraint arc.
C
C In order to be explicit, we make the following defini-
C tions. A constraint region is the open interior of a
C simple closed positively oriented polygonal curve defined
C by an ordered sequence of three or more distinct nodes
C (constraint nodes) P(1),P(2),...,P(K), such that P(I) is
C adjacent to P(I+1) for I = 1,...,K with P(K+1) = P(1).
C Thus, the constraint region is on the left (and may have
C nonfinite area) as the sequence of constraint nodes is
C traversed in the specified order. The constraint regions
C must not contain nodes and must not overlap. The region
C R is the convex hull of the nodes with constraint regions
C excluded.
C
C Note that the terms boundary node and boundary arc are
C reserved for nodes and arcs on the boundary of the convex
C hull of the nodes.
C
C The algorithm is as follows: given a triangulation
C which includes one or more sets of constraint nodes, the
C corresponding adjacencies (constraint arcs) are forced to
C be present (Subroutine EDGE). Any additional new arcs
C required are chosen to be locally optimal (satisfy the
C modified circumcircle property).
C
C
C On input:
C
C NCC = Number of constraint curves (constraint re-
C gions). NCC .GE. 0.
C
C LCC = Array of length NCC (or dummy array of length
C 1 if NCC = 0) containing the index (for X, Y,
C and LEND) of the first node of constraint I in
C LCC(I) for I = 1 to NCC. Thus, constraint I
C contains K = LCC(I+1) - LCC(I) nodes, K .GE.
C 3, stored in (X,Y) locations LCC(I), ...,
C LCC(I+1)-1, where LCC(NCC+1) = N+1.
C
C N = Number of nodes in the triangulation, including
C constraint nodes. N .GE. 3.
C
C X,Y = Arrays of length N containing the coordinates
C of the nodes with non-constraint nodes in the
C first LCC(1)-1 locations, followed by NCC se-
C quences of constraint nodes. Only one of
C these sequences may be specified in clockwise
C order to represent an exterior constraint
C curve (a constraint region with nonfinite
C area).
C
C The above parameters are not altered by this routine.
C
C LWK = Length of IWK. This must be at least 2*NI
C where NI is the maximum number of arcs which
C intersect a constraint arc to be added. NI
C is bounded by N-3.
C
C IWK = Integer work array of length LWK (used by
C Subroutine EDGE to add constraint arcs).
C
C LIST,LPTR,LEND = Data structure defining the trian-
C gulation. Refer to Subroutine
C TRMESH.
C
C On output:
C
C LWK = Required length of IWK unless IER = 1 or IER =
C 3. In the case of IER = 1, LWK is not altered
C from its input value.
C
C IWK = Array containing the endpoint indexes of the
C new arcs which were swapped in by the last
C call to Subroutine EDGE.
C
C LIST,LPTR,LEND = Triangulation data structure with
C all constraint arcs present unless
C IER .NE. 0. These arrays are not
C altered if IER = 1.
C
C IER = Error indicator:
C IER = 0 if no errors were encountered.
C IER = 1 if NCC, N, or an LCC entry is outside
C its valid range, or LWK .LT. 0 on
C input.
C IER = 2 if more space is required in IWK.
C IER = 3 if the triangulation data structure is
C invalid, or failure (in EDGE or OPTIM)
C was caused by collinear nodes on the
C convex hull boundary. An error mes-
C sage is written to logical unit 6 in
C this case.
C IER = 4 if intersecting constraint arcs were
C encountered.
C IER = 5 if a constraint region contains a
C node.
C
C Modules required by ADDCST: EDGE, LEFT, LSTPTR, OPTIM,
C SWAP, SWPTST
C
C Intrinsic functions called by ADDCST: ABS, MAX
C
C***********************************************************
C
INTEGER I, IFRST, ILAST, K, KBAK, KFOR, KN, LCCIP1,
. LP, LPB, LPF, LPL, LW, LWD2, N1, N2
LWD2 = LWK/2
C
C Test for errors in input parameters.
C
IER = 1
IF (NCC .LT. 0 .OR. LWK .LT. 0) RETURN
IF (NCC .EQ. 0) THEN
IF (N .LT. 3) RETURN
LWK = 0
GO TO 9
ELSE
LCCIP1 = N+1
DO 1 I = NCC,1,-1
IF (LCCIP1 - LCC(I) .LT. 3) RETURN
LCCIP1 = LCC(I)
1 CONTINUE
IF (LCCIP1 .LT. 1) RETURN
ENDIF
C
C Force the presence of constraint arcs. The outer loop is
C on constraints in reverse order. IFRST and ILAST are
C the first and last nodes of constraint I.
C
LWK = 0
IFRST = N+1
DO 3 I = NCC,1,-1
ILAST = IFRST - 1
IFRST = LCC(I)
C
C Inner loop on constraint arcs N1-N2 in constraint I.
C
N1 = ILAST
DO 2 N2 = IFRST,ILAST
LW = LWD2
CALL EDGE (N1,N2,X,Y, LW,IWK,LIST,LPTR,LEND, IER)
LWK = MAX(LWK,2*LW)
IF (IER .EQ. 4) IER = 3
IF (IER .NE. 0) RETURN
N1 = N2
2 CONTINUE
3 CONTINUE
C
C Test for errors. The outer loop is on constraint I with
C first and last nodes IFRST and ILAST, and the inner loop
C is on constraint nodes K with (KBAK,K,KFOR) a subse-
C quence of constraint I.
C
IER = 4
IFRST = N+1
DO 8 I = NCC,1,-1
ILAST = IFRST - 1
IFRST = LCC(I)
KBAK = ILAST
DO 7 K = IFRST,ILAST
KFOR = K + 1
IF (K .EQ. ILAST) KFOR = IFRST
C
C Find the LIST pointers LPF and LPB of KFOR and KBAK as
C neighbors of K.
C
LPF = 0
LPB = 0
LPL = LEND(K)
LP = LPL
C
4 LP = LPTR(LP)
KN = ABS(LIST(LP))
IF (KN .EQ. KFOR) LPF = LP
IF (KN .EQ. KBAK) LPB = LP
IF (LP .NE. LPL) GO TO 4
C
C A pair of intersecting constraint arcs was encountered
C if and only if a constraint arc is missing (introduc-
C tion of the second caused the first to be swapped out).
C
IF (LPF .EQ. 0 .OR. LPB .EQ. 0) RETURN
C
C Loop on neighbors KN of node K which follow KFOR and
C precede KBAK. The constraint region contains no nodes
C if and only if all such nodes KN are in constraint I.
C
LP = LPF
5 LP = LPTR(LP)
IF (LP .EQ. LPB) GO TO 6
KN = ABS(LIST(LP))
IF (KN .LT. IFRST .OR. KN .GT. ILAST) GO TO 10
GO TO 5
C
C Bottom of loop.
C
6 KBAK = K
7 CONTINUE
8 CONTINUE
C
C No errors encountered.
C
9 IER = 0
RETURN
C
C A constraint region contains a node.
C
10 IER = 5
RETURN
END
SUBROUTINE ADDNOD (K,XK,YK,IST,NCC, LCC,N,X,Y,LIST,
. LPTR,LEND,LNEW, IER)
INTEGER K, IST, NCC, LCC(*), N, LIST(*), LPTR(*),
. LEND(*), LNEW, IER
REAL XK, YK, X(*), Y(*)
C
C***********************************************************
C
C From TRIPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C renka@cs.unt.edu
C 06/27/98
C
C Given a triangulation of N nodes in the plane created by
C Subroutine TRMESH or TRMSHR, this subroutine updates the
C data structure with the addition of a new node in position
C K. If node K is inserted into X and Y (K .LE. N) rather
C than appended (K = N+1), then a corresponding insertion
C must be performed in any additional arrays associated
C with the nodes. For example, an array of data values Z
C must be shifted down to open up position K for the new
C value: set Z(I+1) to Z(I) for I = N,N-1,...,K. For
C optimal efficiency, new nodes should be appended whenever
C possible. Insertion is necessary, however, to add a non-
C constraint node when constraints are present (refer to
C Subroutine ADDCST).
C
C Note that a constraint node cannot be added by this
C routine. In order to insert a constraint node, it is
C necessary to add the node with no constraints present
C (call this routine with NCC = 0), update LCC by increment-
C ing the appropriate entries, and then create (or restore)
C the constraints by a call to ADDCST.
C
C The algorithm consists of the following steps: node K
C is located relative to the triangulation (TRFIND), its
C index is added to the data structure (INTADD or BDYADD),
C and a sequence of swaps (SWPTST and SWAP) are applied to
C the arcs opposite K so that all arcs incident on node K
C and opposite node K (excluding constraint arcs) are local-
C ly optimal (satisfy the circumcircle test). Thus, if a
C (constrained) Delaunay triangulation is input, a (con-
C strained) Delaunay triangulation will result. All indexes
C are incremented as necessary for an insertion.
C
C
C On input:
C
C K = Nodal index (index for X, Y, and LEND) of the
C new node to be added. 1 .LE. K .LE. LCC(1).
C (K .LE. N+1 if NCC=0).
C
C XK,YK = Cartesian coordinates of the new node (to be
C stored in X(K) and Y(K)). The node must not
C lie in a constraint region.
C
C IST = Index of a node at which TRFIND begins the
C search. Search time depends on the proximity
C of this node to node K. 1 .LE. IST .LE. N.
C
C NCC = Number of constraint curves. NCC .GE. 0.
C
C The above parameters are not altered by this routine.
C
C LCC = List of constraint curve starting indexes (or
C dummy array of length 1 if NCC = 0). Refer to
C Subroutine ADDCST.
C
C N = Number of nodes in the triangulation before K is
C added. N .GE. 3. Note that N will be incre-
C mented following the addition of node K.
C
C X,Y = Arrays of length at least N+1 containing the
C Cartesian coordinates of the nodes in the
C first N positions with non-constraint nodes
C in the first LCC(1)-1 locations if NCC > 0.
C
C LIST,LPTR,LEND,LNEW = Data structure associated with
C the triangulation of nodes 1
C to N. The arrays must have
C sufficient length for N+1
C nodes. Refer to TRMESH.
C
C On output:
C
C LCC = List of constraint curve starting indexes in-
C cremented by 1 to reflect the insertion of K
C unless NCC = 0 or (IER .NE. 0 and IER .NE.
C -4).
C
C N = Number of nodes in the triangulation including K
C unless IER .NE. 0 and IER .NE. -4. Note that
C all comments refer to the input value of N.
C
C X,Y = Arrays updated with the insertion of XK and YK
C in the K-th positions (node I+1 was node I be-
C fore the insertion for I = K to N if K .LE. N)
C unless IER .NE. 0 and IER .NE. -4.
C
C LIST,LPTR,LEND,LNEW = Data structure updated with
C the addition of node K unless
C IER .NE. 0 and IER .NE. -4.
C
C IER = Error indicator:
C IER = 0 if no errors were encountered.
C IER = -1 if K, IST, NCC, N, or an LCC entry is
C outside its valid range on input.
C IER = -2 if all nodes (including K) are col-
C linear.
C IER = L if nodes L and K coincide for some L.
C IER = -3 if K lies in a constraint region.
C IER = -4 if an error flag is returned by SWAP
C implying that the triangulation
C (geometry) was bad on input.
C
C The errors conditions are tested in the order
C specified.
C
C Modules required by ADDNOD: BDYADD, CRTRI, INDXCC,
C INSERT, INTADD, JRAND,
C LEFT, LSTPTR, SWAP,
C SWPTST, TRFIND
C
C Intrinsic function called by ADDNOD: ABS
C
C***********************************************************
C
INTEGER INDXCC, LSTPTR
INTEGER I, I1, I2, I3, IBK, IO1, IO2, IN1, KK, L,
. LCCIP1, LP, LPF, LPO1, NM1
LOGICAL CRTRI, SWPTST
KK = K
C
C Test for an invalid input parameter.
C
IF (KK .LT. 1 .OR. IST .LT. 1 .OR. IST .GT. N
. .OR. NCC .LT. 0 .OR. N .LT. 3) GO TO 7
LCCIP1 = N+1
DO 1 I = NCC,1,-1
IF (LCCIP1-LCC(I) .LT. 3) GO TO 7
LCCIP1 = LCC(I)
1 CONTINUE
IF (KK .GT. LCCIP1) GO TO 7
C
C Find a triangle (I1,I2,I3) containing K or the rightmost
C (I1) and leftmost (I2) visible boundary nodes as viewed
C from node K.
C
CALL TRFIND (IST,XK,YK,N,X,Y,LIST,LPTR,LEND, I1,I2,I3)
C
C Test for collinear nodes, duplicate nodes, and K lying in
C a constraint region.
C
IF (I1 .EQ. 0) GO TO 8
IF (I3 .NE. 0) THEN
L = I1
IF (XK .EQ. X(L) .AND. YK .EQ. Y(L)) GO TO 9
L = I2
IF (XK .EQ. X(L) .AND. YK .EQ. Y(L)) GO TO 9
L = I3
IF (XK .EQ. X(L) .AND. YK .EQ. Y(L)) GO TO 9
IF (NCC .GT. 0 .AND. CRTRI(NCC,LCC,I1,I2,I3) )
. GO TO 10
ELSE
C
C K is outside the convex hull of the nodes and lies in a
C constraint region iff an exterior constraint curve is
C present.
C
IF (NCC .GT. 0 .AND. INDXCC(NCC,LCC,N,LIST,LEND)
. .NE. 0) GO TO 10
ENDIF
C
C No errors encountered.
C
IER = 0
NM1 = N
N = N + 1
IF (KK .LT. N) THEN
C
C Open a slot for K in X, Y, and LEND, and increment all
C nodal indexes which are greater than or equal to K.
C Note that LIST, LPTR, and LNEW are not yet updated with
C either the neighbors of K or the edges terminating on K.
C
DO 2 IBK = NM1,KK,-1
X(IBK+1) = X(IBK)
Y(IBK+1) = Y(IBK)
LEND(IBK+1) = LEND(IBK)
2 CONTINUE
DO 3 I = 1,NCC
LCC(I) = LCC(I) + 1
3 CONTINUE
L = LNEW - 1
DO 4 I = 1,L
IF (LIST(I) .GE. KK) LIST(I) = LIST(I) + 1
IF (LIST(I) .LE. -KK) LIST(I) = LIST(I) - 1
4 CONTINUE
IF (I1 .GE. KK) I1 = I1 + 1
IF (I2 .GE. KK) I2 = I2 + 1
IF (I3 .GE. KK) I3 = I3 + 1
ENDIF
C
C Insert K into X and Y, and update LIST, LPTR, LEND, and
C LNEW with the arcs containing node K.
C
X(KK) = XK
Y(KK) = YK
IF (I3 .EQ. 0) THEN
CALL BDYADD (KK,I1,I2, LIST,LPTR,LEND,LNEW )
ELSE
CALL INTADD (KK,I1,I2,I3, LIST,LPTR,LEND,LNEW )
ENDIF
C
C Initialize variables for optimization of the triangula-
C tion.
C
LP = LEND(KK)
LPF = LPTR(LP)
IO2 = LIST(LPF)
LPO1 = LPTR(LPF)
IO1 = ABS(LIST(LPO1))
C
C Begin loop: find the node opposite K.
C
5 LP = LSTPTR(LEND(IO1),IO2,LIST,LPTR)
IF (LIST(LP) .LT. 0) GO TO 6
LP = LPTR(LP)
IN1 = ABS(LIST(LP))
IF ( CRTRI(NCC,LCC,IO1,IO2,IN1) ) GO TO 6
C
C Swap test: if a swap occurs, two new arcs are
C opposite K and must be tested.
C
IF ( .NOT. SWPTST(IN1,KK,IO1,IO2,X,Y) ) GO TO 6
CALL SWAP (IN1,KK,IO1,IO2, LIST,LPTR,LEND, LPO1)
IF (LPO1 .EQ. 0) GO TO 11
IO1 = IN1
GO TO 5
C
C No swap occurred. Test for termination and reset
C IO2 and IO1.
C
6 IF (LPO1 .EQ. LPF .OR. LIST(LPO1) .LT. 0) RETURN
IO2 = IO1
LPO1 = LPTR(LPO1)
IO1 = ABS(LIST(LPO1))
GO TO 5
C
C A parameter is outside its valid range on input.
C
7 IER = -1
RETURN
C
C All nodes are collinear.
C
8 IER = -2
RETURN
C
C Nodes L and K coincide.
C
9 IER = L
RETURN
C
C Node K lies in a constraint region.
C
10 IER = -3
RETURN
C
C Zero pointer returned by SWAP.
C
11 IER = -4
RETURN
END
REAL FUNCTION AREAP (X,Y,NB,NODES)
INTEGER NB, NODES(NB)
REAL X(*), Y(*)
C
C***********************************************************
C
C From TRIPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C renka@cs.unt.edu
C 09/21/90
C
C Given a sequence of NB points in the plane, this func-
C tion computes the signed area bounded by the closed poly-
C gonal curve which passes through the points in the
C specified order. Each simple closed curve is positively
C oriented (bounds positive area) if and only if the points
C are specified in counterclockwise order. The last point
C of the curve is taken to be the first point specified, and
C this point should therefore not be specified twice.
C
C The area of a triangulation may be computed by calling
C AREAP with values of NB and NODES determined by Subroutine
C BNODES.
C
C
C On input:
C
C X,Y = Arrays of length N containing the Cartesian
C coordinates of a set of points in the plane
C for some N .GE. NB.
C
C NB = Length of NODES.
C
C NODES = Array of length NB containing the ordered
C sequence of nodal indexes (in the range
C 1 to N) which define the polygonal curve.
C
C Input parameters are not altered by this function.
C
C On output:
C
C AREAP = Signed area bounded by the polygonal curve,
C or zero if NB < 3.
C
C Modules required by AREAP: None
C
C***********************************************************
C
INTEGER I, ND1, ND2, NNB
REAL A
C
C Local parameters:
C
C A = Partial sum of signed (and doubled) trapezoid
C areas
C I = DO-loop and NODES index
C ND1,ND2 = Elements of NODES
C NNB = Local copy of NB
C
NNB = NB
A = 0.
IF (NNB .LT. 3) GO TO 2
ND2 = NODES(NNB)
C
C Loop on line segments NODES(I-1) -> NODES(I), where
C NODES(0) = NODES(NB), adding twice the signed trapezoid
C areas (integrals of the linear interpolants) to A.
C
DO 1 I = 1,NNB
ND1 = ND2
ND2 = NODES(I)
A = A + (X(ND2)-X(ND1))*(Y(ND1)+Y(ND2))
1 CONTINUE
C
C A contains twice the negative signed area of the region.
C
2 AREAP = -A/2.
RETURN
END
SUBROUTINE BDYADD (KK,I1,I2, LIST,LPTR,LEND,LNEW )
INTEGER KK, I1, I2, LIST(*), LPTR(*), LEND(*), LNEW
C
C***********************************************************
C
C From TRIPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C renka@cs.unt.edu
C 02/22/91
C
C This subroutine adds a boundary node to a triangulation
C of a set of points in the plane. The data structure is
C updated with the insertion of node KK, but no optimization
C is performed.
C
C
C On input:
C
C KK = Index of a node to be connected to the sequence
C of all visible boundary nodes. KK .GE. 1 and
C KK must not be equal to I1 or I2.
C
C I1 = First (rightmost as viewed from KK) boundary
C node in the triangulation which is visible from
C node KK (the line segment KK-I1 intersects no
C arcs.
C
C I2 = Last (leftmost) boundary node which is visible
C from node KK. I1 and I2 may be determined by
C Subroutine TRFIND.
C
C The above parameters are not altered by this routine.
C
C LIST,LPTR,LEND,LNEW = Triangulation data structure
C created by TRMESH or TRMSHR.
C Nodes I1 and I2 must be in-
C cluded in the triangulation.
C
C On output:
C
C LIST,LPTR,LEND,LNEW = Data structure updated with
C the addition of node KK. Node
C KK is connected to I1, I2, and
C all boundary nodes in between.
C
C Module required by BDYADD: INSERT
C
C***********************************************************
C
INTEGER K, LP, LSAV, N1, N2, NEXT, NSAV
K = KK
N1 = I1
N2 = I2
C
C Add K as the last neighbor of N1.
C
LP = LEND(N1)
LSAV = LPTR(LP)
LPTR(LP) = LNEW
LIST(LNEW) = -K
LPTR(LNEW) = LSAV
LEND(N1) = LNEW
LNEW = LNEW + 1
NEXT = -LIST(LP)
LIST(LP) = NEXT
NSAV = NEXT
C
C Loop on the remaining boundary nodes between N1 and N2,
C adding K as the first neighbor.
C
1 LP = LEND(NEXT)
CALL INSERT (K,LP,LIST,LPTR,LNEW)
IF (NEXT .EQ. N2) GO TO 2
NEXT = -LIST(LP)
LIST(LP) = NEXT
GO TO 1
C
C Add the boundary nodes between N1 and N2 as neighbors
C of node K.
C
2 LSAV = LNEW
LIST(LNEW) = N1
LPTR(LNEW) = LNEW + 1
LNEW = LNEW + 1
NEXT = NSAV
C
3 IF (NEXT .EQ. N2) GO TO 4
LIST(LNEW) = NEXT
LPTR(LNEW) = LNEW + 1
LNEW = LNEW + 1
LP = LEND(NEXT)
NEXT = LIST(LP)
GO TO 3
C
4 LIST(LNEW) = -N2
LPTR(LNEW) = LSAV
LEND(K) = LNEW
LNEW = LNEW + 1
RETURN
END
SUBROUTINE BNODES (N,LIST,LPTR,LEND, NODES,NB,NA,NT)
INTEGER N, LIST(*), LPTR(*), LEND(N), NODES(*), NB,
. NA, NT
C
C***********************************************************
C
C From TRIPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C renka@cs.unt.edu
C 09/01/88
C
C Given a triangulation of N points in the plane, this
C subroutine returns an array containing the indexes, in
C counterclockwise order, of the nodes on the boundary of
C the convex hull of the set of points.
C
C
C On input:
C
C N = Number of nodes in the triangulation. N .GE. 3.
C
C LIST,LPTR,LEND = Data structure defining the trian-
C gulation. Refer to Subroutine
C TRMESH.
C
C The above parameters are not altered by this routine.
C
C NODES = Integer array of length at least NB
C (NB .LE. N).
C
C On output:
C
C NODES = Ordered sequence of boundary node indexes
C in the range 1 to N.
C
C NB = Number of boundary nodes.
C
C NA,NT = Number of arcs and triangles, respectively,
C in the triangulation.
C
C Modules required by BNODES: None
C
C***********************************************************
C
INTEGER K, LP, N0, NST
C
C Set NST to the first boundary node encountered.
C
NST = 1
1 LP = LEND(NST)
IF (LIST(LP) .LT. 0) GO TO 2
NST = NST + 1
GO TO 1
C
C Initialization.
C
2 NODES(1) = NST
K = 1
N0 = NST
C
C Traverse the boundary in counterclockwise order.
C
3 LP = LEND(N0)
LP = LPTR(LP)
N0 = LIST(LP)
IF (N0 .EQ. NST) GO TO 4
K = K + 1
NODES(K) = N0
GO TO 3
C
C Termination.
C
4 NB = K
NT = 2*N - NB - 2
NA = NT + N - 1
RETURN
END
SUBROUTINE CIRCUM (X1,Y1,X2,Y2,X3,Y3,RATIO, XC,YC,CR,
. SA,AR)
LOGICAL RATIO
REAL X1, Y1, X2, Y2, X3, Y3, XC, YC, CR, SA, AR
C
C***********************************************************
C
C From TRIPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C renka@cs.unt.edu
C 12/10/96
C
C Given three vertices defining a triangle, this subrou-
C tine returns the circumcenter, circumradius, signed
C triangle area, and, optionally, the aspect ratio of the
C triangle.
C
C
C On input:
C
C X1,...,Y3 = Cartesian coordinates of the vertices.
C
C RATIO = Logical variable with value TRUE if and only
C if the aspect ratio is to be computed.
C
C Input parameters are not altered by this routine.
C
C On output:
C
C XC,YC = Cartesian coordinates of the circumcenter
C (center of the circle defined by the three
C points) unless SA = 0, in which XC and YC
C are not altered.
C
C CR = Circumradius (radius of the circle defined by
C the three points) unless SA = 0 (infinite
C radius), in which case CR is not altered.
C
C SA = Signed triangle area with positive value if
C and only if the vertices are specified in
C counterclockwise order: (X3,Y3) is strictly
C to the left of the directed line from (X1,Y1)
C toward (X2,Y2).
C
C AR = Aspect ratio r/CR, where r is the radius of the
C inscribed circle, unless RATIO = FALSE, in
C which case AR is not altered. AR is in the
C range 0 to .5, with value 0 iff SA = 0 and
C value .5 iff the vertices define an equilateral
C triangle.
C
C Modules required by CIRCUM: None
C
C Intrinsic functions called by CIRCUM: ABS, SQRT
C
C***********************************************************
C
INTEGER I
REAL DS(3), FX, FY, U(3), V(3)
C
C Set U(K) and V(K) to the x and y components, respectively,
C of the directed edge opposite vertex K.
C
U(1) = X3 - X2
U(2) = X1 - X3
U(3) = X2 - X1
V(1) = Y3 - Y2
V(2) = Y1 - Y3
V(3) = Y2 - Y1
C
C Set SA to the signed triangle area.
C
SA = (U(1)*V(2) - U(2)*V(1))/2.
IF (SA .EQ. 0.) THEN
IF (RATIO) AR = 0.
RETURN
ENDIF
C
C Set DS(K) to the squared distance from the origin to
C vertex K.
C
DS(1) = X1*X1 + Y1*Y1
DS(2) = X2*X2 + Y2*Y2
DS(3) = X3*X3 + Y3*Y3
C
C Compute factors of XC and YC.
C
FX = 0.
FY = 0.
DO 1 I = 1,3
FX = FX - DS(I)*V(I)
FY = FY + DS(I)*U(I)
1 CONTINUE
XC = FX/(4.*SA)
YC = FY/(4.*SA)
CR = SQRT( (XC-X1)**2 + (YC-Y1)**2 )
IF (.NOT. RATIO) RETURN
C
C Compute the squared edge lengths and aspect ratio.
C
DO 2 I = 1,3
DS(I) = U(I)*U(I) + V(I)*V(I)
2 CONTINUE
AR = 2.*ABS(SA)/
. ( (SQRT(DS(1)) + SQRT(DS(2)) + SQRT(DS(3)))*CR )
RETURN
END
LOGICAL FUNCTION CRTRI (NCC,LCC,I1,I2,I3)
INTEGER NCC, LCC(*), I1, I2, I3
C
C***********************************************************
C
C From TRIPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C renka@cs.unt.edu
C 08/14/91
C
C This function returns TRUE if and only if triangle (I1,
C I2,I3) lies in a constraint region.
C
C
C On input:
C
C NCC,LCC = Constraint data structure. Refer to Sub-
C routine ADDCST.
C
C I1,I2,I3 = Nodal indexes of the counterclockwise-
C ordered vertices of a triangle.
C
C Input parameters are altered by this function.
C
C CRTRI = TRUE iff (I1,I2,I3) is a constraint region
C triangle.
C
C Note that input parameters are not tested for validity.
C
C Modules required by CRTRI: None
C
C Intrinsic functions called by CRTRI: MAX, MIN
C
C***********************************************************
C
INTEGER I, IMAX, IMIN
IMAX = MAX(I1,I2,I3)
C
C Find the index I of the constraint containing IMAX.
C
I = NCC + 1
1 I = I - 1
IF (I .LE. 0) GO TO 2
IF (IMAX .LT. LCC(I)) GO TO 1
IMIN = MIN(I1,I2,I3)
C
C P lies in a constraint region iff I1, I2, and I3 are nodes
C of the same constraint (IMIN >= LCC(I)), and (IMIN,IMAX)
C is (I1,I3), (I2,I1), or (I3,I2).
C
CRTRI = IMIN .GE. LCC(I) .AND. ((IMIN .EQ. I1 .AND.
. IMAX .EQ. I3) .OR. (IMIN .EQ. I2 .AND.
. IMAX .EQ. I1) .OR. (IMIN .EQ. I3 .AND.
. IMAX .EQ. I2))
RETURN
C
C NCC .LE. 0 or all vertices are non-constraint nodes.
C
2 CRTRI = .FALSE.
RETURN
END
SUBROUTINE DELARC (N,IO1,IO2, LIST,LPTR,LEND,
. LNEW, IER)
INTEGER N, IO1, IO2, LIST(*), LPTR(*), LEND(N), LNEW,
. IER
C
C***********************************************************
C
C From TRIPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C renka@cs.unt.edu
C 11/12/94
C
C This subroutine deletes a boundary arc from a triangula-
C tion. It may be used to remove a null triangle from the
C convex hull boundary. Note, however, that if the union of
C triangles is rendered nonconvex, Subroutines DELNOD, EDGE,
C and TRFIND may fail. Thus, Subroutines ADDCST, ADDNOD,
C DELNOD, EDGE, and NEARND should not be called following
C an arc deletion.
C
C
C On input:
C
C N = Number of nodes in the triangulation. N .GE. 4.
C
C IO1,IO2 = Indexes (in the range 1 to N) of a pair of
C adjacent boundary nodes defining the arc
C to be removed.
C
C The above parameters are not altered by this routine.
C
C LIST,LPTR,LEND,LNEW = Triangulation data structure
C created by TRMESH or TRMSHR.
C
C On output:
C
C LIST,LPTR,LEND,LNEW = Data structure updated with
C the removal of arc IO1-IO2
C unless IER > 0.
C
C IER = Error indicator:
C IER = 0 if no errors were encountered.
C IER = 1 if N, IO1, or IO2 is outside its valid
C range, or IO1 = IO2.
C IER = 2 if IO1-IO2 is not a boundary arc.
C IER = 3 if the node opposite IO1-IO2 is al-
C ready a boundary node, and thus IO1
C or IO2 has only two neighbors or a
C deletion would result in two triangu-
C lations sharing a single node.
C IER = 4 if one of the nodes is a neighbor of
C the other, but not vice versa, imply-
C ing an invalid triangulation data
C structure.
C
C Modules required by DELARC: DELNB, LSTPTR
C
C Intrinsic function called by DELARC: ABS
C
C***********************************************************
C
INTEGER LSTPTR
INTEGER LP, LPH, LPL, N1, N2, N3
N1 = IO1
N2 = IO2
C
C Test for errors, and set N1->N2 to the directed boundary
C edge associated with IO1-IO2: (N1,N2,N3) is a triangle
C for some N3.
C
IF (N .LT. 4 .OR. N1 .LT. 1 .OR. N1 .GT. N .OR.
. N2 .LT. 1 .OR. N2 .GT. N .OR. N1 .EQ. N2) THEN
IER = 1
RETURN
ENDIF
C
LPL = LEND(N2)
IF (-LIST(LPL) .NE. N1) THEN
N1 = N2
N2 = IO1
LPL = LEND(N2)
IF (-LIST(LPL) .NE. N1) THEN
IER = 2
RETURN
ENDIF
ENDIF
C
C Set N3 to the node opposite N1->N2 (the second neighbor
C of N1), and test for error 3 (N3 already a boundary
C node).
C
LPL = LEND(N1)
LP = LPTR(LPL)
LP = LPTR(LP)
N3 = ABS(LIST(LP))
LPL = LEND(N3)
IF (LIST(LPL) .LE. 0) THEN
IER = 3
RETURN
ENDIF
C
C Delete N2 as a neighbor of N1, making N3 the first
C neighbor, and test for error 4 (N2 not a neighbor
C of N1). Note that previously computed pointers may
C no longer be valid following the call to DELNB.
C
CALL DELNB (N1,N2,N, LIST,LPTR,LEND,LNEW, LPH)
IF (LPH .LT. 0) THEN
IER = 4
RETURN
ENDIF
C
C Delete N1 as a neighbor of N2, making N3 the new last
C neighbor.
C
CALL DELNB (N2,N1,N, LIST,LPTR,LEND,LNEW, LPH)
C
C Make N3 a boundary node with first neighbor N2 and last
C neighbor N1.
C
LP = LSTPTR(LEND(N3),N1,LIST,LPTR)
LEND(N3) = LP
LIST(LP) = -N1
C
C No errors encountered.
C
IER = 0
RETURN
END
SUBROUTINE DELNB (N0,NB,N, LIST,LPTR,LEND,LNEW, LPH)
INTEGER N0, NB, N, LIST(*), LPTR(*), LEND(N), LNEW,
. LPH
C
C***********************************************************
C
C From TRIPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C renka@cs.unt.edu
C 07/30/98
C
C This subroutine deletes a neighbor NB from the adjacency
C list of node N0 (but N0 is not deleted from the adjacency
C list of NB) and, if NB is a boundary node, makes N0 a
C boundary node. For pointer (LIST index) LPH to NB as a
C neighbor of N0, the empty LIST,LPTR location LPH is filled
C in with the values at LNEW-1, pointer LNEW-1 (in LPTR and
C possibly in LEND) is changed to LPH, and LNEW is decremen-
C ted. This requires a search of LEND and LPTR entailing an
C expected operation count of O(N).
C
C
C On input:
C
C N0,NB = Indexes, in the range 1 to N, of a pair of
C nodes such that NB is a neighbor of N0.
C (N0 need not be a neighbor of NB.)
C
C N = Number of nodes in the triangulation. N .GE. 3.
C
C The above parameters are not altered by this routine.
C
C LIST,LPTR,LEND,LNEW = Data structure defining the
C triangulation.
C
C On output:
C
C LIST,LPTR,LEND,LNEW = Data structure updated with
C the removal of NB from the ad-
C jacency list of N0 unless
C LPH < 0.
C
C LPH = List pointer to the hole (NB as a neighbor of
C N0) filled in by the values at LNEW-1 or error
C indicator:
C LPH > 0 if no errors were encountered.
C LPH = -1 if N0, NB, or N is outside its valid
C range.
C LPH = -2 if NB is not a neighbor of N0.
C
C Modules required by DELNB: None
C
C Intrinsic function called by DELNB: ABS
C
C***********************************************************
C
INTEGER I, LNW, LP, LPB, LPL, LPP, NN
C
C Local parameters:
C
C I = DO-loop index
C LNW = LNEW-1 (output value of LNEW)
C LP = LIST pointer of the last neighbor of NB
C LPB = Pointer to NB as a neighbor of N0
C LPL = Pointer to the last neighbor of N0
C LPP = Pointer to the neighbor of N0 that precedes NB
C NN = Local copy of N
C
NN = N
C
C Test for error 1.
C
IF (N0 .LT. 1 .OR. N0 .GT. NN .OR. NB .LT. 1 .OR.
. NB .GT. NN .OR. NN .LT. 3) THEN
LPH = -1
RETURN
ENDIF
C
C Find pointers to neighbors of N0:
C
C LPL points to the last neighbor,
C LPP points to the neighbor NP preceding NB, and
C LPB points to NB.
C
LPL = LEND(N0)
LPP = LPL
LPB = LPTR(LPP)
1 IF (LIST(LPB) .EQ. NB) GO TO 2
LPP = LPB
LPB = LPTR(LPP)
IF (LPB .NE. LPL) GO TO 1
C
C Test for error 2 (NB not found).
C
IF (ABS(LIST(LPB)) .NE. NB) THEN
LPH = -2
RETURN
ENDIF
C
C NB is the last neighbor of N0. Make NP the new last
C neighbor and, if NB is a boundary node, then make N0
C a boundary node.
C
LEND(N0) = LPP
LP = LEND(NB)
IF (LIST(LP) .LT. 0) LIST(LPP) = -LIST(LPP)
GO TO 3
C
C NB is not the last neighbor of N0. If NB is a boundary
C node and N0 is not, then make N0 a boundary node with
C last neighbor NP.
C
2 LP = LEND(NB)
IF (LIST(LP) .LT. 0 .AND. LIST(LPL) .GT. 0) THEN
LEND(N0) = LPP
LIST(LPP) = -LIST(LPP)
ENDIF
C
C Update LPTR so that the neighbor following NB now fol-
C lows NP, and fill in the hole at location LPB.
C
3 LPTR(LPP) = LPTR(LPB)
LNW = LNEW-1
LIST(LPB) = LIST(LNW)
LPTR(LPB) = LPTR(LNW)
DO 4 I = NN,1,-1
IF (LEND(I) .EQ. LNW) THEN
LEND(I) = LPB
GO TO 5
ENDIF
4 CONTINUE
C
5 DO 6 I = 1,LNW-1
IF (LPTR(I) .EQ. LNW) THEN
LPTR(I) = LPB
ENDIF
6 CONTINUE
C
C No errors encountered.
C
LNEW = LNW
LPH = LPB
RETURN
END
SUBROUTINE DELNOD (K,NCC, LCC,N,X,Y,LIST,LPTR,LEND,
. LNEW,LWK,IWK, IER)
INTEGER K, NCC, LCC(*), N, LIST(*), LPTR(*),
. LEND(*), LNEW, LWK, IWK(2,*), IER
REAL X(*), Y(*)
C
C***********************************************************
C
C From TRIPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C renka@cs.unt.edu
C 06/28/98
C
C This subroutine deletes node K (along with all arcs
C incident on node K) from a triangulation of N nodes in the
C plane, and inserts arcs as necessary to produce a triangu-
C lation of the remaining N-1 nodes. If a Delaunay triangu-
C lation is input, a Delaunay triangulation will result, and
C thus, DELNOD reverses the effect of a call to Subroutine
C ADDNOD.
C
C Note that a constraint node cannot be deleted by this
C routine. In order to delete a constraint node, it is
C necessary to call this routine with NCC = 0, decrement the
C appropriate LCC entries (LCC(I) such that LCC(I) > K), and
C then create (or restore) the constraints by a call to Sub-
C routine ADDCST.
C
C
C On input:
C
C K = Index (for X and Y) of the node to be deleted.
C 1 .LE. K .LT. LCC(1). (K .LE. N if NCC=0).
C
C NCC = Number of constraint curves. NCC .GE. 0.
C
C The above parameters are not altered by this routine.
C
C LCC = List of constraint curve starting indexes (or
C dummy array of length 1 if NCC = 0). Refer to
C Subroutine ADDCST.
C
C N = Number of nodes in the triangulation on input.
C N .GE. 4. Note that N will be decremented
C following the deletion.
C
C X,Y = Arrays of length N containing the coordinates
C of the nodes with non-constraint nodes in the
C first LCC(1)-1 locations if NCC > 0.
C
C LIST,LPTR,LEND,LNEW = Data structure defining the
C triangulation. Refer to Sub-
C routine TRMESH.
C
C LWK = Number of columns reserved for IWK. LWK must
C be at least NNB-3, where NNB is the number of
C neighbors of node K, including an extra
C pseudo-node if K is a boundary node.
C
C IWK = Integer work array dimensioned 2 by LWK (or
C array of length .GE. 2*LWK).
C
C On output:
C
C LCC = List of constraint curve starting indexes de-
C cremented by 1 to reflect the deletion of K
C unless NCC = 0 or 1 .LE. IER .LE. 4.
C
C N = New number of nodes (input value minus one) un-
C less 1 .LE. IER .LE. 4.
C
C X,Y = Updated arrays of length N-1 containing nodal
C coordinates (with elements K+1,...,N shifted
C up a position and thus overwriting element K)
C unless 1 .LE. IER .LE. 4. (N here denotes the
C input value.)
C
C LIST,LPTR,LEND,LNEW = Updated triangulation data
C structure reflecting the dele-
C tion unless IER .NE. 0. Note
C that the data structure may
C have been altered if IER .GE.
C 3.
C
C LWK = Number of IWK columns required unless IER = 1
C or IER = 3.
C
C IWK = Indexes of the endpoints of the new arcs added
C unless LWK = 0 or 1 .LE. IER .LE. 4. (Arcs
C are associated with columns, or pairs of
C adjacent elements if IWK is declared as a
C singly-subscripted array.)
C
C IER = Error indicator:
C IER = 0 if no errors were encountered.
C IER = 1 if K, NCC, N, or an LCC entry is out-
C side its valid range or LWK < 0 on
C input.
C IER = 2 if more space is required in IWK.
C Refer to LWK.
C IER = 3 if the triangulation data structure is
C invalid on input.
C IER = 4 if K is an interior node with 4 or
C more neighbors, and the number of
C neighbors could not be reduced to 3
C by swaps. This could be caused by
C floating point errors with collinear
C nodes or by an invalid data structure.
C IER = 5 if an error flag was returned by
C OPTIM. An error message is written
C to the standard output unit in this
C event.
C
C Note that the deletion may result in all remaining nodes
C being collinear. This situation is not flagged.
C
C Modules required by DELNOD: DELNB, LEFT, LSTPTR, NBCNT,
C OPTIM, SWAP, SWPTST
C
C Intrinsic function called by DELNOD: ABS
C
C***********************************************************
C
INTEGER LSTPTR, NBCNT
LOGICAL LEFT
INTEGER I, IERR, IWL, J, LCCIP1, LNW, LP, LP21, LPF,
. LPH, LPL, LPL2, LPN, LWKL, N1, N2, NFRST, NIT,
. NL, NN, NNB, NR
LOGICAL BDRY
REAL X1, X2, XL, XR, Y1, Y2, YL, YR
C
C Set N1 to K and NNB to the number of neighbors of N1 (plus
C one if N1 is a boundary node), and test for errors. LPF
C and LPL are LIST indexes of the first and last neighbors
C of N1, IWL is the number of IWK columns containing arcs,
C and BDRY is TRUE iff N1 is a boundary node.
C
N1 = K
NN = N
IF (NCC .LT. 0 .OR. N1 .LT. 1 .OR. NN .LT. 4 .OR.
. LWK .LT. 0) GO TO 21
LCCIP1 = NN+1
DO 1 I = NCC,1,-1
IF (LCCIP1-LCC(I) .LT. 3) GO TO 21
LCCIP1 = LCC(I)
1 CONTINUE
IF (N1 .GE. LCCIP1) GO TO 21
LPL = LEND(N1)
LPF = LPTR(LPL)
NNB = NBCNT(LPL,LPTR)
BDRY = LIST(LPL) .LT. 0
IF (BDRY) NNB = NNB + 1
IF (NNB .LT. 3) GO TO 23
LWKL = LWK
LWK = NNB - 3
IF (LWKL .LT. LWK) GO TO 22
IWL = 0
IF (NNB .EQ. 3) GO TO 5
C
C Initialize for loop on arcs N1-N2 for neighbors N2 of N1,
C beginning with the second neighbor. NR and NL are the
C neighbors preceding and following N2, respectively, and
C LP indexes NL. The loop is exited when all possible
C swaps have been applied to arcs incident on N1. If N1
C is interior, the number of neighbors will be reduced
C to 3.
C
X1 = X(N1)
Y1 = Y(N1)
NFRST = LIST(LPF)
NR = NFRST
XR = X(NR)
YR = Y(NR)
LP = LPTR(LPF)
N2 = LIST(LP)
X2 = X(N2)
Y2 = Y(N2)
LP = LPTR(LP)
C
C Top of loop: set NL to the neighbor following N2.
C
2 NL = ABS(LIST(LP))
IF (NL .EQ. NFRST .AND. BDRY) GO TO 5
XL = X(NL)
YL = Y(NL)
C
C Test for a convex quadrilateral. To avoid an incorrect
C test caused by collinearity, use the fact that if N1
C is a boundary node, then N1 LEFT NR->NL and if N2 is
C a boundary node, then N2 LEFT NL->NR.
C
LPL2 = LEND(N2)
IF ( (BDRY .OR. LEFT(XR,YR,XL,YL,X1,Y1)) .AND.
. (LIST(LPL2) .LT. 0 .OR.
. LEFT(XL,YL,XR,YR,X2,Y2)) ) GO TO 3
C
C Nonconvex quadrilateral -- no swap is possible.
C
NR = N2
XR = X2
YR = Y2
GO TO 4
C
C The quadrilateral defined by adjacent triangles
C (N1,N2,NL) and (N2,N1,NR) is convex. Swap in
C NL-NR and store it in IWK. Indexes larger than N1
C must be decremented since N1 will be deleted from
C X and Y.
C
3 CALL SWAP (NL,NR,N1,N2, LIST,LPTR,LEND, LP21)
IWL = IWL + 1
IF (NL .LE. N1) THEN
IWK(1,IWL) = NL
ELSE
IWK(1,IWL) = NL - 1
ENDIF
IF (NR .LE. N1) THEN
IWK(2,IWL) = NR
ELSE
IWK(2,IWL) = NR - 1
ENDIF
C
C Recompute the LIST indexes LPL,LP and decrement NNB.
C
LPL = LEND(N1)
NNB = NNB - 1
IF (NNB .EQ. 3) GO TO 5
LP = LSTPTR(LPL,NL,LIST,LPTR)
IF (NR .EQ. NFRST) GO TO 4
C
C NR is not the first neighbor of N1.
C Back up and test N1-NR for a swap again: Set N2 to
C NR and NR to the previous neighbor of N1 -- the
C neighbor of NR which follows N1. LP21 points to NL
C as a neighbor of NR.
C
N2 = NR
X2 = XR
Y2 = YR
LP21 = LPTR(LP21)
LP21 = LPTR(LP21)
NR = ABS(LIST(LP21))
XR = X(NR)
YR = Y(NR)
GO TO 2
C
C Bottom of loop -- test for invalid termination.
C
4 IF (N2 .EQ. NFRST) GO TO 24
N2 = NL
X2 = XL
Y2 = YL
LP = LPTR(LP)
GO TO 2
C
C Delete N1 from the adjacency list of N2 for all neighbors
C N2 of N1. LPL points to the last neighbor of N1.
C LNEW is stored in local variable LNW.
C
5 LP = LPL
LNW = LNEW
C
C Loop on neighbors N2 of N1, beginning with the first.
C
6 LP = LPTR(LP)
N2 = ABS(LIST(LP))
CALL DELNB (N2,N1,N, LIST,LPTR,LEND,LNW, LPH)
IF (LPH .LT. 0) GO TO 23
C
C LP and LPL may require alteration.
C
IF (LPL .EQ. LNW) LPL = LPH
IF (LP .EQ. LNW) LP = LPH
IF (LP .NE. LPL) GO TO 6
C
C Delete N1 from X, Y, and LEND, and remove its adjacency
C list from LIST and LPTR. LIST entries (nodal indexes)
C which are larger than N1 must be decremented.
C
NN = NN - 1
IF (N1 .GT. NN) GO TO 9
DO 7 I = N1,NN
X(I) = X(I+1)
Y(I) = Y(I+1)
LEND(I) = LEND(I+1)
7 CONTINUE
C
DO 8 I = 1,LNW-1
IF (LIST(I) .GT. N1) LIST(I) = LIST(I) - 1
IF (LIST(I) .LT. -N1) LIST(I) = LIST(I) + 1
8 CONTINUE
C
C For LPN = first to last neighbors of N1, delete the
C preceding neighbor (indexed by LP).
C
C Each empty LIST,LPTR location LP is filled in with the
C values at LNW-1, and LNW is decremented. All pointers
C (including those in LPTR and LEND) with value LNW-1
C must be changed to LP.
C
C LPL points to the last neighbor of N1.
C
9 IF (BDRY) NNB = NNB - 1
LPN = LPL
DO 13 J = 1,NNB
LNW = LNW - 1
LP = LPN
LPN = LPTR(LP)
LIST(LP) = LIST(LNW)
LPTR(LP) = LPTR(LNW)
IF (LPTR(LPN) .EQ. LNW) LPTR(LPN) = LP
IF (LPN .EQ. LNW) LPN = LP
DO 10 I = NN,1,-1
IF (LEND(I) .EQ. LNW) THEN
LEND(I) = LP
GO TO 11
ENDIF
10 CONTINUE
C
11 DO 12 I = LNW-1,1,-1
IF (LPTR(I) .EQ. LNW) LPTR(I) = LP
12 CONTINUE
13 CONTINUE
C
C Decrement LCC entries.
C
DO 14 I = 1,NCC
LCC(I) = LCC(I) - 1
14 CONTINUE
C
C Update N and LNEW, and optimize the patch of triangles
C containing K (on input) by applying swaps to the arcs
C in IWK.
C
N = NN
LNEW = LNW
IF (IWL .GT. 0) THEN
NIT = 4*IWL
CALL OPTIM (X,Y,IWL, LIST,LPTR,LEND,NIT,IWK, IERR)
IF (IERR .NE. 0) GO TO 25
ENDIF
C
C Successful termination.
C
IER = 0
RETURN
C
C Invalid input parameter.
C
21 IER = 1
RETURN
C
C Insufficient space reserved for IWK.
C
22 IER = 2
RETURN
C
C Invalid triangulation data structure. NNB < 3 on input or
C N2 is a neighbor of N1 but N1 is not a neighbor of N2.
C
23 IER = 3
RETURN
C
C K is an interior node with 4 or more neighbors, but the
C number of neighbors could not be reduced.
C
24 IER = 4
RETURN
C
C Error flag returned by OPTIM.
C
25 IER = 5
WRITE (*,100) NIT, IERR
RETURN
100 FORMAT (//5X,'*** Error in OPTIM: NIT = ',I4,
. ', IER = ',I1,' ***'/)
END
SUBROUTINE EDGE (IN1,IN2,X,Y, LWK,IWK,LIST,LPTR,
. LEND, IER)
INTEGER IN1, IN2, LWK, IWK(2,*), LIST(*), LPTR(*),
. LEND(*), IER
REAL X(*), Y(*)
C
C***********************************************************
C
C From TRIPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C renka@cs.unt.edu
C 06/23/98
C
C Given a triangulation of N nodes and a pair of nodal
C indexes IN1 and IN2, this routine swaps arcs as necessary
C to force IN1 and IN2 to be adjacent. Only arcs which
C intersect IN1-IN2 are swapped out. If a Delaunay triangu-
C lation is input, the resulting triangulation is as close
C as possible to a Delaunay triangulation in the sense that
C all arcs other than IN1-IN2 are locally optimal.
C
C A sequence of calls to EDGE may be used to force the
C presence of a set of edges defining the boundary of a non-
C convex and/or multiply connected region (refer to Subrou-
C tine ADDCST), or to introduce barriers into the triangula-
C tion. Note that Subroutine GETNP will not necessarily
C return closest nodes if the triangulation has been con-
C strained by a call to EDGE. However, this is appropriate
C in some applications, such as triangle-based interpolation
C on a nonconvex domain.
C
C
C On input:
C
C IN1,IN2 = Indexes (of X and Y) in the range 1 to N
C defining a pair of nodes to be connected
C by an arc.
C
C X,Y = Arrays of length N containing the Cartesian
C coordinates of the nodes.
C
C The above parameters are not altered by this routine.
C
C LWK = Number of columns reserved for IWK. This must
C be at least NI -- the number of arcs which
C intersect IN1-IN2. (NI is bounded by N-3.)
C
C IWK = Integer work array of length at least 2*LWK.
C
C LIST,LPTR,LEND = Data structure defining the trian-
C gulation. Refer to Subroutine
C TRMESH.
C
C On output:
C
C LWK = Number of arcs which intersect IN1-IN2 (but
C not more than the input value of LWK) unless
C IER = 1 or IER = 3. LWK = 0 if and only if
C IN1 and IN2 were adjacent (or LWK=0) on input.
C
C IWK = Array containing the indexes of the endpoints
C of the new arcs other than IN1-IN2 unless IER
C .GT. 0 or LWK = 0. New arcs to the left of
C IN2-IN1 are stored in the first K-1 columns
C (left portion of IWK), column K contains
C zeros, and new arcs to the right of IN2-IN1
C occupy columns K+1,...,LWK. (K can be deter-
C mined by searching IWK for the zeros.)
C
C LIST,LPTR,LEND = Data structure updated if necessary
C to reflect the presence of an arc
C connecting IN1 and IN2 unless IER
C .NE. 0. The data structure has
C been altered if IER = 4.
C
C IER = Error indicator:
C IER = 0 if no errors were encountered.
C IER = 1 if IN1 .LT. 1, IN2 .LT. 1, IN1 = IN2,
C or LWK .LT. 0 on input.
C IER = 2 if more space is required in IWK.
C IER = 3 if IN1 and IN2 could not be connected
C due to either an invalid data struc-
C ture or collinear nodes (and floating
C point error).
C IER = 4 if an error flag was returned by
C OPTIM.
C
C An error message is written to the standard output unit
C in the case of IER = 3 or IER = 4.
C
C Modules required by EDGE: LEFT, LSTPTR, OPTIM, SWAP,
C SWPTST
C
C Intrinsic function called by EDGE: ABS
C
C***********************************************************
C
LOGICAL LEFT
INTEGER I, IERR, IWC, IWCP1, IWEND, IWF, IWL, LFT, LP,
. LPL, LP21, NEXT, NIT, NL, NR, N0, N1, N2,
. N1FRST, N1LST
REAL DX, DY, X0, Y0, X1, Y1, X2, Y2
C
C Local parameters:
C
C DX,DY = Components of arc N1-N2
C I = DO-loop index and column index for IWK
C IERR = Error flag returned by Subroutine OPTIM
C IWC = IWK index between IWF and IWL -- NL->NR is
C stored in IWK(1,IWC)->IWK(2,IWC)
C IWCP1 = IWC + 1
C IWEND = Input or output value of LWK
C IWF = IWK (column) index of the first (leftmost) arc
C which intersects IN1->IN2
C IWL = IWK (column) index of the last (rightmost) are
C which intersects IN1->IN2
C LFT = Flag used to determine if a swap results in the
C new arc intersecting IN1-IN2 -- LFT = 0 iff
C N0 = IN1, LFT = -1 implies N0 LEFT IN1->IN2,
C and LFT = 1 implies N0 LEFT IN2->IN1
C LP21 = Unused parameter returned by SWAP
C LP = List pointer (index) for LIST and LPTR
C LPL = Pointer to the last neighbor of IN1 or NL
C N0 = Neighbor of N1 or node opposite NR->NL
C N1,N2 = Local copies of IN1 and IN2
C N1FRST = First neighbor of IN1
C N1LST = (Signed) last neighbor of IN1
C NEXT = Node opposite NL->NR
C NIT = Flag or number of iterations employed by OPTIM
C NL,NR = Endpoints of an arc which intersects IN1-IN2
C with NL LEFT IN1->IN2
C X0,Y0 = Coordinates of N0
C X1,Y1 = Coordinates of IN1
C X2,Y2 = Coordinates of IN2
C
C
C Store IN1, IN2, and LWK in local variables and test for
C errors.
C
N1 = IN1
N2 = IN2
IWEND = LWK
IF (N1 .LT. 1 .OR. N2 .LT. 1 .OR. N1 .EQ. N2 .OR.
. IWEND .LT. 0) GO TO 31
C
C Test for N2 as a neighbor of N1. LPL points to the last
C neighbor of N1.
C
LPL = LEND(N1)
N0 = ABS(LIST(LPL))
LP = LPL
1 IF (N0 .EQ. N2) GO TO 30
LP = LPTR(LP)
N0 = LIST(LP)
IF (LP .NE. LPL) GO TO 1
C
C Initialize parameters.
C
IWL = 0
NIT = 0
C
C Store the coordinates of N1 and N2.
C
2 X1 = X(N1)
Y1 = Y(N1)
X2 = X(N2)
Y2 = Y(N2)
C
C Set NR and NL to adjacent neighbors of N1 such that
C NR LEFT N2->N1 and NL LEFT N1->N2,
C (NR Forward N1->N2 or NL Forward N1->N2), and
C (NR Forward N2->N1 or NL Forward N2->N1).
C
C Initialization: Set N1FRST and N1LST to the first and
C (signed) last neighbors of N1, respectively, and
C initialize NL to N1FRST.
C
LPL = LEND(N1)
N1LST = LIST(LPL)
LP = LPTR(LPL)
N1FRST = LIST(LP)
NL = N1FRST
IF (N1LST .LT. 0) GO TO 4
C
C N1 is an interior node. Set NL to the first candidate
C for NR (NL LEFT N2->N1).
C
3 IF ( LEFT(X2,Y2,X1,Y1,X(NL),Y(NL)) ) GO TO 4
LP = LPTR(LP)
NL = LIST(LP)
IF (NL .NE. N1FRST) GO TO 3
C
C All neighbors of N1 are strictly left of N1->N2.
C
GO TO 5
C
C NL = LIST(LP) LEFT N2->N1. Set NR to NL and NL to the
C following neighbor of N1.
C
4 NR = NL
LP = LPTR(LP)
NL = ABS(LIST(LP))
IF ( LEFT(X1,Y1,X2,Y2,X(NL),Y(NL)) ) THEN
C
C NL LEFT N1->N2 and NR LEFT N2->N1. The Forward tests
C are employed to avoid an error associated with
C collinear nodes.
C
DX = X2-X1
DY = Y2-Y1
IF ((DX*(X(NL)-X1)+DY*(Y(NL)-Y1) .GE. 0. .OR.
. DX*(X(NR)-X1)+DY*(Y(NR)-Y1) .GE. 0.) .AND.
. (DX*(X(NL)-X2)+DY*(Y(NL)-Y2) .LE. 0. .OR.
. DX*(X(NR)-X2)+DY*(Y(NR)-Y2) .LE. 0.)) GO TO 6
C
C NL-NR does not intersect N1-N2. However, there is
C another candidate for the first arc if NL lies on
C the line N1-N2.
C
IF ( .NOT. LEFT(X2,Y2,X1,Y1,X(NL),Y(NL)) ) GO TO 5
ENDIF
C
C Bottom of loop.
C
IF (NL .NE. N1FRST) GO TO 4
C
C Either the triangulation is invalid or N1-N2 lies on the
C convex hull boundary and an edge NR->NL (opposite N1 and
C intersecting N1-N2) was not found due to floating point
C error. Try interchanging N1 and N2 -- NIT > 0 iff this
C has already been done.
C
5 IF (NIT .GT. 0) GO TO 33
NIT = 1
N1 = N2
N2 = IN1
GO TO 2
C
C Store the ordered sequence of intersecting edges NL->NR in
C IWK(1,IWL)->IWK(2,IWL).
C
6 IWL = IWL + 1
IF (IWL .GT. IWEND) GO TO 32
IWK(1,IWL) = NL
IWK(2,IWL) = NR
C
C Set NEXT to the neighbor of NL which follows NR.
C
LPL = LEND(NL)
LP = LPTR(LPL)
C
C Find NR as a neighbor of NL. The search begins with
C the first neighbor.
C
7 IF (LIST(LP) .EQ. NR) GO TO 8
LP = LPTR(LP)
IF (LP .NE. LPL) GO TO 7
C
C NR must be the last neighbor, and NL->NR cannot be a
C boundary edge.
C
IF (LIST(LP) .NE. NR) GO TO 33
C
C Set NEXT to the neighbor following NR, and test for
C termination of the store loop.
C
8 LP = LPTR(LP)
NEXT = ABS(LIST(LP))
IF (NEXT .EQ. N2) GO TO 9
C
C Set NL or NR to NEXT.
C
IF ( LEFT(X1,Y1,X2,Y2,X(NEXT),Y(NEXT)) ) THEN
NL = NEXT
ELSE
NR = NEXT
ENDIF
GO TO 6
C
C IWL is the number of arcs which intersect N1-N2.
C Store LWK.
C
9 LWK = IWL
IWEND = IWL
C
C Initialize for edge swapping loop -- all possible swaps
C are applied (even if the new arc again intersects
C N1-N2), arcs to the left of N1->N2 are stored in the
C left portion of IWK, and arcs to the right are stored in
C the right portion. IWF and IWL index the first and last
C intersecting arcs.
C
IWF = 1
C
C Top of loop -- set N0 to N1 and NL->NR to the first edge.
C IWC points to the arc currently being processed. LFT
C .LE. 0 iff N0 LEFT N1->N2.
C
10 LFT = 0
N0 = N1
X0 = X1
Y0 = Y1
NL = IWK(1,IWF)
NR = IWK(2,IWF)
IWC = IWF
C
C Set NEXT to the node opposite NL->NR unless IWC is the
C last arc.
C
11 IF (IWC .EQ. IWL) GO TO 21
IWCP1 = IWC + 1
NEXT = IWK(1,IWCP1)
IF (NEXT .NE. NL) GO TO 16
NEXT = IWK(2,IWCP1)
C
C NEXT RIGHT N1->N2 and IWC .LT. IWL. Test for a possible
C swap.
C
IF ( .NOT. LEFT(X0,Y0,X(NR),Y(NR),X(NEXT),Y(NEXT)) )
. GO TO 14
IF (LFT .GE. 0) GO TO 12
IF ( .NOT. LEFT(X(NL),Y(NL),X0,Y0,X(NEXT),Y(NEXT)) )
. GO TO 14
C
C Replace NL->NR with N0->NEXT.
C
CALL SWAP (NEXT,N0,NL,NR, LIST,LPTR,LEND, LP21)
IWK(1,IWC) = N0
IWK(2,IWC) = NEXT
GO TO 15
C
C Swap NL-NR for N0-NEXT, shift columns IWC+1,...,IWL to
C the left, and store N0-NEXT in the right portion of
C IWK.
C
12 CALL SWAP (NEXT,N0,NL,NR, LIST,LPTR,LEND, LP21)
DO 13 I = IWCP1,IWL
IWK(1,I-1) = IWK(1,I)
IWK(2,I-1) = IWK(2,I)
13 CONTINUE
IWK(1,IWL) = N0
IWK(2,IWL) = NEXT
IWL = IWL - 1
NR = NEXT
GO TO 11
C
C A swap is not possible. Set N0 to NR.
C
14 N0 = NR
X0 = X(N0)
Y0 = Y(N0)
LFT = 1
C
C Advance to the next arc.
C
15 NR = NEXT
IWC = IWC + 1
GO TO 11
C
C NEXT LEFT N1->N2, NEXT .NE. N2, and IWC .LT. IWL.
C Test for a possible swap.
C
16 IF ( .NOT. LEFT(X(NL),Y(NL),X0,Y0,X(NEXT),Y(NEXT)) )
. GO TO 19
IF (LFT .LE. 0) GO TO 17
IF ( .NOT. LEFT(X0,Y0,X(NR),Y(NR),X(NEXT),Y(NEXT)) )
. GO TO 19
C
C Replace NL->NR with NEXT->N0.
C
CALL SWAP (NEXT,N0,NL,NR, LIST,LPTR,LEND, LP21)
IWK(1,IWC) = NEXT
IWK(2,IWC) = N0
GO TO 20
C
C Swap NL-NR for N0-NEXT, shift columns IWF,...,IWC-1 to
C the right, and store N0-NEXT in the left portion of
C IWK.
C
17 CALL SWAP (NEXT,N0,NL,NR, LIST,LPTR,LEND, LP21)
DO 18 I = IWC-1,IWF,-1
IWK(1,I+1) = IWK(1,I)
IWK(2,I+1) = IWK(2,I)
18 CONTINUE
IWK(1,IWF) = N0
IWK(2,IWF) = NEXT
IWF = IWF + 1
GO TO 20
C
C A swap is not possible. Set N0 to NL.
C
19 N0 = NL
X0 = X(N0)
Y0 = Y(N0)
LFT = -1
C
C Advance to the next arc.
C
20 NL = NEXT
IWC = IWC + 1
GO TO 11
C
C N2 is opposite NL->NR (IWC = IWL).
C
21 IF (N0 .EQ. N1) GO TO 24
IF (LFT .LT. 0) GO TO 22
C
C N0 RIGHT N1->N2. Test for a possible swap.
C
IF ( .NOT. LEFT(X0,Y0,X(NR),Y(NR),X2,Y2) ) GO TO 10
C
C Swap NL-NR for N0-N2 and store N0-N2 in the right
C portion of IWK.
C
CALL SWAP (N2,N0,NL,NR, LIST,LPTR,LEND, LP21)
IWK(1,IWL) = N0
IWK(2,IWL) = N2
IWL = IWL - 1
GO TO 10
C
C N0 LEFT N1->N2. Test for a possible swap.
C
22 IF ( .NOT. LEFT(X(NL),Y(NL),X0,Y0,X2,Y2) ) GO TO 10
C
C Swap NL-NR for N0-N2, shift columns IWF,...,IWL-1 to the
C right, and store N0-N2 in the left portion of IWK.
C
CALL SWAP (N2,N0,NL,NR, LIST,LPTR,LEND, LP21)
I = IWL
23 IWK(1,I) = IWK(1,I-1)
IWK(2,I) = IWK(2,I-1)
I = I - 1
IF (I .GT. IWF) GO TO 23
IWK(1,IWF) = N0
IWK(2,IWF) = N2
IWF = IWF + 1
GO TO 10
C
C IWF = IWC = IWL. Swap out the last arc for N1-N2 and
C store zeros in IWK.
C
24 CALL SWAP (N2,N1,NL,NR, LIST,LPTR,LEND, LP21)
IWK(1,IWC) = 0
IWK(2,IWC) = 0
C
C Optimization procedure --
C
IF (IWC .GT. 1) THEN
C
C Optimize the set of new arcs to the left of IN1->IN2.
C
NIT = 3*(IWC-1)
CALL OPTIM (X,Y,IWC-1, LIST,LPTR,LEND,NIT,IWK, IERR)
IF (IERR .NE. 0) GO TO 34
ENDIF
IF (IWC .LT. IWEND) THEN
C
C Optimize the set of new arcs to the right of IN1->IN2.
C
NIT = 3*(IWEND-IWC)
CALL OPTIM (X,Y,IWEND-IWC, LIST,LPTR,LEND,NIT,
. IWK(1,IWC+1), IERR)
IF (IERR .NE. 0) GO TO 34
ENDIF
C
C Successful termination.
C
IER = 0
RETURN
C
C IN1 and IN2 were adjacent on input.
C
30 IER = 0
RETURN
C
C Invalid input parameter.
C
31 IER = 1
RETURN
C
C Insufficient space reserved for IWK.
C
32 IER = 2
RETURN
C
C Invalid triangulation data structure or collinear nodes
C on convex hull boundary.
C
33 IER = 3
WRITE (*,130) IN1, IN2
130 FORMAT (//5X,'*** Error in EDGE: Invalid triangula',
. 'tion or null triangles on boundary'/
. 9X,'IN1 =',I4,', IN2=',I4/)
RETURN
C
C Error flag returned by OPTIM.
C
34 IER = 4
WRITE (*,140) NIT, IERR
140 FORMAT (//5X,'*** Error in OPTIM: NIT = ',I4,
. ', IER = ',I1,' ***'/)
RETURN
END
SUBROUTINE GETNP (NCC,LCC,N,X,Y,LIST,LPTR,LEND,
. L, NPTS,DS, IER)
INTEGER NCC, LCC(*), N, LIST(*), LPTR(*), LEND(N),
. L, NPTS(L), IER
REAL X(N), Y(N), DS(L)
C
C***********************************************************
C
C From TRIPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C renka@cs.unt.edu
C 11/12/94
C
C Given a triangulation of N nodes and an array NPTS con-
C taining the indexes of L-1 nodes ordered by distance from
C NPTS(1), this subroutine sets NPTS(L) to the index of the
C next node in the sequence -- the node, other than NPTS(1),
C ...,NPTS(L-1), which is closest to NPTS(1). Thus, the
C ordered sequence of K closest nodes to N1 (including N1)
C may be determined by K-1 calls to GETNP with NPTS(1) = N1
C and L = 2,3,...,K for K .GE. 2. Note that NPTS must in-
C clude constraint nodes as well as non-constraint nodes.
C Thus, a sequence of K1 closest non-constraint nodes to N1
C must be obtained as a subset of the closest K2 nodes to N1
C for some K2 .GE. K1.
C
C The terms closest and distance have special definitions
C when constraint nodes are present in the triangulation.
C Nodes N1 and N2 are said to be visible from each other if
C and only if the line segment N1-N2 intersects no con-
C straint arc (except possibly itself) and is not an interi-
C or constraint arc (arc whose interior lies in a constraint
C region). A path from N1 to N2 is an ordered sequence of
C nodes, with N1 first and N2 last, such that adjacent path
C elements are visible from each other. The path length is
C the sum of the Euclidean distances between adjacent path
C nodes. Finally, the distance from N1 to N2 is defined to
C be the length of the shortest path from N1 to N2.
C
C The algorithm uses the property of a Delaunay triangula-
C tion that the K-th closest node to N1 is a neighbor of one
C of the K-1 closest nodes to N1. With the definition of
C distance used here, this property holds when constraints
C are present as long as non-constraint arcs are locally
C optimal.
C
C
C On input:
C
C NCC = Number of constraints. NCC .GE. 0.
C
C LCC = List of constraint curve starting indexes (or
C dummy array of length 1 if NCC = 0). Refer to
C Subroutine ADDCST.
C
C N = Number of nodes in the triangulation. N .GE. 3.
C
C X,Y = Arrays of length N containing the coordinates
C of the nodes with non-constraint nodes in the
C first LCC(1)-1 locations if NCC > 0.
C
C LIST,LPTR,LEND = Triangulation data structure. Re-
C fer to Subroutine TRMESH.
C
C L = Number of nodes in the sequence on output. 2
C .LE. L .LE. N.
C
C NPTS = Array of length .GE. L containing the indexes
C of the L-1 closest nodes to NPTS(1) in the
C first L-1 locations.
C
C DS = Array of length .GE. L containing the distance
C (defined above) between NPTS(1) and NPTS(I) in
C the I-th position for I = 1,...,L-1. Thus,
C DS(1) = 0.
C
C Input parameters other than NPTS(L) and DS(L) are not
C altered by this routine.
C
C On output:
C
C NPTS = Array updated with the index of the L-th
C closest node to NPTS(1) in position L unless
C IER .NE. 0.
C
C DS = Array updated with the distance between NPTS(1)
C and NPTS(L) in position L unless IER .NE. 0.
C
C IER = Error indicator:
C IER = 0 if no errors were encountered.
C IER = -1 if NCC, N, L, or an LCC entry is
C outside its valid range on input.
C IER = K if NPTS(K) is not a valid index in
C the range 1 to N.
C
C Module required by GETNP: INTSEC
C
C Intrinsic functions called by GETNP: ABS, MIN, SQRT
C
C***********************************************************
C
LOGICAL INTSEC
INTEGER I, IFRST, ILAST, J, K, KM1, LCC1, LM1, LP,
. LPCL, LPK, LPKL, N1, NC, NF1, NF2, NJ, NK,
. NKBAK, NKFOR, NL, NN
LOGICAL ISW, VIS, NCF, NJF, SKIP, SKSAV, LFT1, LFT2,
. LFT12
REAL DC, DL, X1, XC, XJ, XK, Y1, YC, YJ, YK
C
C Store parameters in local variables and test for errors.
C LCC1 indexes the first constraint node.
C
IER = -1
NN = N
LCC1 = NN+1
LM1 = L-1
IF (NCC .LT. 0 .OR. LM1 .LT. 1 .OR. LM1 .GE. NN)
. RETURN
IF (NCC .EQ. 0) THEN
IF (NN .LT. 3) RETURN
ELSE
DO 1 I = NCC,1,-1
IF (LCC1 - LCC(I) .LT. 3) RETURN
LCC1 = LCC(I)
1 CONTINUE
IF (LCC1 .LT. 1) RETURN
ENDIF
C
C Test for an invalid index in NPTS.
C
DO 2 K = 1,LM1
NK = NPTS(K)
IF (NK .LT. 1 .OR. NK .GT. NN) THEN
IER = K
RETURN
ENDIF
2 CONTINUE
C
C Store N1 = NPTS(1) and mark the elements of NPTS.
C
N1 = NPTS(1)
X1 = X(N1)
Y1 = Y(N1)
DO 3 K = 1,LM1
NK = NPTS(K)
LEND(NK) = -LEND(NK)
3 CONTINUE
C
C Candidates NC for NL = NPTS(L) are the unmarked visible
C neighbors of nodes NK in NPTS. ISW is an initialization
C switch set to .TRUE. when NL and its distance DL from N1
C have been initialized with the first candidate encount-
C ered.
C
ISW = .FALSE.
DL = 0.
C
C Loop on marked nodes NK = NPTS(K). LPKL indexes the last
C neighbor of NK in LIST.
C
DO 16 K = 1,LM1
KM1 = K - 1
NK = NPTS(K)
XK = X(NK)
YK = Y(NK)
LPKL = -LEND(NK)
NKFOR = 0
NKBAK = 0
VIS = .TRUE.
IF (NK .GE. LCC1) THEN
C
C NK is a constraint node. Set NKFOR and NKBAK to the
C constraint nodes which follow and precede NK. IFRST
C and ILAST are set to the first and last nodes in the
C constraint containing NK.
C
IFRST = NN + 1
DO 4 I = NCC,1,-1
ILAST = IFRST - 1
IFRST = LCC(I)
IF (NK .GE. IFRST) GO TO 5
4 CONTINUE
C
5 IF (NK .LT. ILAST) THEN
NKFOR = NK + 1
ELSE
NKFOR = IFRST
ENDIF
IF (NK .GT. IFRST) THEN
NKBAK = NK - 1
ELSE
NKBAK = ILAST
ENDIF
C
C Initialize VIS to TRUE iff NKFOR precedes NKBAK in the
C adjacency list for NK -- the first neighbor is visi-
C ble and is not NKBAK.
C
LPK = LPKL
6 LPK = LPTR(LPK)
NC = ABS(LIST(LPK))
IF (NC .NE. NKFOR .AND. NC .NE. NKBAK) GO TO 6
VIS = NC .EQ. NKFOR
ENDIF
C
C Loop on neighbors NC of NK, bypassing marked and nonvis-
C ible neighbors.
C
LPK = LPKL
7 LPK = LPTR(LPK)
NC = ABS(LIST(LPK))
IF (NC .EQ. NKBAK) VIS = .TRUE.
C
C VIS = .FALSE. iff NK-NC is an interior constraint arc
C (NK is a constraint node and NC lies strictly between
C NKFOR and NKBAK).
C
IF (.NOT. VIS) GO TO 15
IF (NC .EQ. NKFOR) VIS = .FALSE.
IF (LEND(NC) .LT. 0) GO TO 15
C
C Initialize distance DC between N1 and NC to Euclidean
C distance.
C
XC = X(NC)
YC = Y(NC)
DC = SQRT((XC-X1)*(XC-X1) + (YC-Y1)*(YC-Y1))
IF (ISW .AND. DC .GE. DL) GO TO 15
IF (K .EQ. 1) GO TO 14
C
C K .GE. 2. Store the pointer LPCL to the last neighbor
C of NC.
C
LPCL = LEND(NC)
C
C Set DC to the length of the shortest path from N1 to NC
C which has not previously been encountered and which is
C a viable candidate for the shortest path from N1 to NL.
C This is Euclidean distance iff NC is visible from N1.
C Since the shortest path from N1 to NL contains only ele-
C ments of NPTS which are constraint nodes (in addition to
C N1 and NL), only these need be considered for the path
C from N1 to NC. Thus, for distance function D(A,B) and
C J = 1,...,K, DC = min(D(N1,NJ) + D(NJ,NC)) over con-
C straint nodes NJ = NPTS(J) which are visible from NC.
C
DO 13 J = 1,KM1
NJ = NPTS(J)
IF (J .GT. 1 .AND. NJ .LT. LCC1) GO TO 13
C
C If NC is a visible neighbor of NJ, a path from N1 to NC
C containing NJ has already been considered. Thus, NJ may
C be bypassed if it is adjacent to NC.
C
LP = LPCL
8 LP = LPTR(LP)
IF ( NJ .EQ. ABS(LIST(LP)) ) GO TO 12
IF (LP .NE. LPCL) GO TO 8
C
C NJ is a constraint node (unless J=1) not adjacent to NC,
C and is visible from NC iff NJ-NC is not intersected by
C a constraint arc. Loop on constraints I in reverse
C order --
C
XJ = X(NJ)
YJ = Y(NJ)
IFRST = NN+1
DO 11 I = NCC,1,-1
ILAST = IFRST - 1
IFRST = LCC(I)
NF1 = ILAST
NCF = NF1 .EQ. NC
NJF = NF1 .EQ. NJ
SKIP = NCF .OR. NJF
C
C Loop on boundary constraint arcs NF1-NF2 which contain
C neither NC nor NJ. NCF and NJF are TRUE iff NC (or NJ)
C has been encountered in the constraint, and SKIP =
C .TRUE. iff NF1 = NC or NF1 = NJ.
C
DO 10 NF2 = IFRST,ILAST
IF (NF2 .EQ. NC) NCF = .TRUE.
IF (NF2 .EQ. NJ) NJF = .TRUE.
SKSAV = SKIP
SKIP = NF2 .EQ. NC .OR. NF2 .EQ. NJ
C
C The last constraint arc in the constraint need not be
C tested if none of the arcs have been skipped.
C
IF ( SKSAV .OR. SKIP .OR.
. (NF2 .EQ. ILAST .AND.
. .NOT. NCF .AND. .NOT. NJF) ) GO TO 9
IF ( INTSEC(X(NF1),Y(NF1),X(NF2),Y(NF2),
. XC,YC,XJ,YJ) ) GO TO 12
9 NF1 = NF2
10 CONTINUE
IF (.NOT. NCF .OR. .NOT. NJF) GO TO 11
C
C NC and NJ are constraint nodes in the same constraint.
C NC-NJ is intersected by an interior constraint arc iff
C 1) NC LEFT NF2->NF1 and (NJ LEFT NF1->NC and NJ LEFT
C NC->NF2) or
C 2) NC .NOT. LEFT NF2->NF1 and (NJ LEFT NF1->NC or
C NJ LEFT NC->NF2),
C where NF1, NC, NF2 are consecutive constraint nodes.
C
IF (NC .NE. IFRST) THEN
NF1 = NC - 1
ELSE
NF1 = ILAST
ENDIF
IF (NC .NE. ILAST) THEN
NF2 = NC + 1
ELSE
NF2 = IFRST
ENDIF
LFT1 = (XC-X(NF1))*(YJ-Y(NF1)) .GE.
. (XJ-X(NF1))*(YC-Y(NF1))
LFT2 = (X(NF2)-XC)*(YJ-YC) .GE.
. (XJ-XC)*(Y(NF2)-YC)
LFT12 = (X(NF1)-X(NF2))*(YC-Y(NF2)) .GE.
. (XC-X(NF2))*(Y(NF1)-Y(NF2))
IF ( (LFT1 .AND. LFT2) .OR. (.NOT. LFT12
. .AND. (LFT1 .OR. LFT2)) ) GO TO 12
11 CONTINUE
C
C NJ is visible from NC. Exit the loop with DC = Euclidean
C distance if J = 1.
C
IF (J .EQ. 1) GO TO 14
DC = MIN(DC,DS(J) + SQRT((XC-XJ)*(XC-XJ) +
. (YC-YJ)*(YC-YJ)))
GO TO 13
C
C NJ is not visible from NC or is adjacent to NC. Initial-
C ize DC with D(N1,NK) + D(NK,NC) if J = 1.
C
12 IF (J .EQ. 1) DC = DS(K) + SQRT((XC-XK)*(XC-XK)
. + (YC-YK)*(YC-YK))
13 CONTINUE
C
C Compare DC with DL.
C
IF (ISW .AND. DC .GE. DL) GO TO 15
C
C The first (or a closer) candidate for NL has been
C encountered.
C
14 NL = NC
DL = DC
ISW = .TRUE.
15 IF (LPK .NE. LPKL) GO TO 7
16 CONTINUE
C
C Unmark the elements of NPTS and store NL and DL.
C
DO 17 K = 1,LM1
NK = NPTS(K)
LEND(NK) = -LEND(NK)
17 CONTINUE
NPTS(L) = NL
DS(L) = DL
IER = 0
RETURN
END
INTEGER FUNCTION INDXCC (NCC,LCC,N,LIST,LEND)
INTEGER NCC, LCC(*), N, LIST(*), LEND(N)
C
C***********************************************************
C
C From TRIPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C renka@cs.unt.edu
C 08/25/91
C
C Given a constrained Delaunay triangulation, this func-
C tion returns the index, if any, of an exterior constraint
C curve (an unbounded constraint region). An exterior con-
C straint curve is assumed to be present if and only if the
C clockwise-ordered sequence of boundary nodes is a subse-
C quence of a constraint node sequence. The triangulation
C adjacencies corresponding to constraint edges may or may
C not have been forced by a call to ADDCST, and the con-
C straint region may or may not be valid (contain no nodes).
C
C
C On input:
C
C NCC = Number of constraints. NCC .GE. 0.
C
C LCC = List of constraint curve starting indexes (or
C dummy array of length 1 if NCC = 0). Refer to
C Subroutine ADDCST.
C
C N = Number of nodes in the triangulation. N .GE. 3.
C
C LIST,LEND = Data structure defining the triangula-
C tion. Refer to Subroutine TRMESH.
C
C Input parameters are not altered by this function. Note
C that the parameters are not tested for validity.
C
C On output:
C
C INDXCC = Index of the exterior constraint curve, if
C present, or 0 otherwise.
C
C Modules required by INDXCC: None
C
C***********************************************************
C
INTEGER I, IFRST, ILAST, LP, N0, NST, NXT
INDXCC = 0
IF (NCC .LT. 1) RETURN
C
C Set N0 to the boundary node with smallest index.
C
N0 = 0
1 N0 = N0 + 1
LP = LEND(N0)
IF (LIST(LP) .GT. 0) GO TO 1
C
C Search in reverse order for the constraint I, if any, that
C contains N0. IFRST and ILAST index the first and last
C nodes in constraint I.
C
I = NCC
ILAST = N
2 IFRST = LCC(I)
IF (N0 .GE. IFRST) GO TO 3
IF (I .EQ. 1) RETURN
I = I - 1
ILAST = IFRST - 1
GO TO 2
C
C N0 is in constraint I which indexes an exterior constraint
C curve iff the clockwise-ordered sequence of boundary
C node indexes beginning with N0 is increasing and bounded
C above by ILAST.
C
3 NST = N0
C
4 NXT = -LIST(LP)
IF (NXT .EQ. NST) GO TO 5
IF (NXT .LE. N0 .OR. NXT .GT. ILAST) RETURN
N0 = NXT
LP = LEND(N0)
GO TO 4
C
C Constraint I contains the boundary node sequence as a
C subset.
C
5 INDXCC = I
RETURN
END
SUBROUTINE INSERT (K,LP, LIST,LPTR,LNEW )
INTEGER K, LP, LIST(*), LPTR(*), LNEW
C
C***********************************************************
C
C From TRIPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C renka@cs.unt.edu
C 09/01/88
C
C This subroutine inserts K as a neighbor of N1 following
C N2, where LP is the LIST pointer of N2 as a neighbor of
C N1. Note that, if N2 is the last neighbor of N1, K will
C become the first neighbor (even if N1 is a boundary node).
C
C
C On input:
C
C K = Index of the node to be inserted.
C
C LP = LIST pointer of N2 as a neighbor of N1.
C
C The above parameters are not altered by this routine.
C
C LIST,LPTR,LNEW = Data structure defining the trian-
C gulation. Refer to Subroutine
C TRMESH.
C
C On output:
C
C LIST,LPTR,LNEW = Data structure updated with the
C addition of node K.
C
C Modules required by INSERT: None
C
C***********************************************************
C
INTEGER LSAV
C
LSAV = LPTR(LP)
LPTR(LP) = LNEW
LIST(LNEW) = K
LPTR(LNEW) = LSAV
LNEW = LNEW + 1
RETURN
END
SUBROUTINE INTADD (KK,I1,I2,I3, LIST,LPTR,LEND,LNEW )
INTEGER KK, I1, I2, I3, LIST(*), LPTR(*), LEND(*),
. LNEW
C
C***********************************************************
C
C From TRIPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C renka@cs.unt.edu
C 02/22/91
C
C This subroutine adds an interior node to a triangulation
C of a set of points in the plane. The data structure is
C updated with the insertion of node KK into the triangle
C whose vertices are I1, I2, and I3. No optimization of the
C triangulation is performed.
C
C
C On input:
C
C KK = Index of the node to be inserted. KK .GE. 1
C and KK must not be equal to I1, I2, or I3.
C
C I1,I2,I3 = Indexes of the counterclockwise-ordered
C sequence of vertices of a triangle which
C contains node KK.
C
C The above parameters are not altered by this routine.
C
C LIST,LPTR,LEND,LNEW = Data structure defining the
C triangulation. Refer to Sub-
C routine TRMESH. Triangle
C (I1,I2,I3) must be included
C in the triangulation.
C
C On output:
C
C LIST,LPTR,LEND,LNEW = Data structure updated with
C the addition of node KK. KK
C will be connected to nodes I1,
C I2, and I3.
C
C Modules required by INTADD: INSERT, LSTPTR
C
C***********************************************************
C
INTEGER LSTPTR
INTEGER K, LP, N1, N2, N3
K = KK
C
C Initialization.
C
N1 = I1
N2 = I2
N3 = I3
C
C Add K as a neighbor of I1, I2, and I3.
C
LP = LSTPTR(LEND(N1),N2,LIST,LPTR)
CALL INSERT (K,LP,LIST,LPTR,LNEW)
LP = LSTPTR(LEND(N2),N3,LIST,LPTR)
CALL INSERT (K,LP,LIST,LPTR,LNEW)
LP = LSTPTR(LEND(N3),N1,LIST,LPTR)
CALL INSERT (K,LP,LIST,LPTR,LNEW)
C
C Add I1, I2, and I3 as neighbors of K.
C
LIST(LNEW) = N1
LIST(LNEW+1) = N2
LIST(LNEW+2) = N3
LPTR(LNEW) = LNEW + 1
LPTR(LNEW+1) = LNEW + 2
LPTR(LNEW+2) = LNEW
LEND(K) = LNEW + 2
LNEW = LNEW + 3
RETURN
END
LOGICAL FUNCTION INTSEC (X1,Y1,X2,Y2,X3,Y3,X4,Y4)
REAL X1, Y1, X2, Y2, X3, Y3, X4, Y4
C
C***********************************************************
C
C From TRIPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C renka@cs.unt.edu
C 09/01/88
C
C Given a pair of line segments P1-P2 and P3-P4, this
C function returns the value .TRUE. if and only if P1-P2
C shares one or more points with P3-P4. The line segments
C include their endpoints, and the four points need not be
C distinct. Thus, either line segment may consist of a
C single point, and the segments may meet in a V (which is
C treated as an intersection). Note that an incorrect
C decision may result from floating point error if the four
C endpoints are nearly collinear.
C
C
C On input:
C
C X1,Y1 = Coordinates of P1.
C
C X2,Y2 = Coordinates of P2.
C
C X3,Y3 = Coordinates of P3.
C
C X4,Y4 = Coordinates of P4.
C
C Input parameters are not altered by this function.
C
C On output:
C
C INTSEC = Logical value defined above.
C
C Modules required by INTSEC: None
C
C***********************************************************
C
REAL A, B, D, DX12, DX31, DX34, DY12, DY31, DY34
C
C Test for overlap between the smallest rectangles that
C contain the line segments and have sides parallel to
C the axes.
C
IF ((X1 .LT. X3 .AND. X1 .LT. X4 .AND. X2 .LT. X3
. .AND. X2 .LT. X4) .OR.
. (X1 .GT. X3 .AND. X1 .GT. X4 .AND. X2 .GT. X3
. .AND. X2 .GT. X4) .OR.
. (Y1 .LT. Y3 .AND. Y1 .LT. Y4 .AND. Y2 .LT. Y3
. .AND. Y2 .LT. Y4) .OR.
. (Y1 .GT. Y3 .AND. Y1 .GT. Y4 .AND. Y2 .GT. Y3
. .AND. Y2 .GT. Y4)) THEN
INTSEC = .FALSE.
RETURN
ENDIF
C
C Compute A = P4-P3 X P1-P3, B = P2-P1 X P1-P3, and
C D = P2-P1 X P4-P3 (Z components).
C
DX12 = X2 - X1
DY12 = Y2 - Y1
DX34 = X4 - X3
DY34 = Y4 - Y3
DX31 = X1 - X3
DY31 = Y1 - Y3
A = DX34*DY31 - DX31*DY34
B = DX12*DY31 - DX31*DY12
D = DX12*DY34 - DX34*DY12
IF (D .EQ. 0.) GO TO 1
C
C D .NE. 0 and the point of intersection of the lines de-
C fined by the line segments is P = P1 + (A/D)*(P2-P1) =
C P3 + (B/D)*(P4-P3).
C
INTSEC = A/D .GE. 0. .AND. A/D .LE. 1. .AND.
. B/D .GE. 0. .AND. B/D .LE. 1.
RETURN
C
C D .EQ. 0 and thus either the line segments are parallel,
C or one (or both) of them is a single point.
C
1 INTSEC = A .EQ. 0. .AND. B .EQ. 0.
RETURN
END
INTEGER FUNCTION JRAND (N, IX,IY,IZ )
INTEGER N, IX, IY, IZ
C
C***********************************************************
C
C From STRIPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C renka@cs.unt.edu
C 07/28/98
C
C This function returns a uniformly distributed pseudo-
C random integer in the range 1 to N.
C
C
C On input:
C
C N = Maximum value to be returned.
C
C N is not altered by this function.
C
C IX,IY,IZ = Integer seeds initialized to values in
C the range 1 to 30,000 before the first
C call to JRAND, and not altered between
C subsequent calls (unless a sequence of
C random numbers is to be repeated by
C reinitializing the seeds).
C
C On output:
C
C IX,IY,IZ = Updated integer seeds.
C
C JRAND = Random integer in the range 1 to N.
C
C Reference: B. A. Wichmann and I. D. Hill, "An Efficient
C and Portable Pseudo-random Number Generator",
C Applied Statistics, Vol. 31, No. 2, 1982,
C pp. 188-190.
C
C Modules required by JRAND: None
C
C Intrinsic functions called by JRAND: INT, MOD, REAL
C
C***********************************************************
C
REAL U, X
C
C Local parameters:
C
C U = Pseudo-random number uniformly distributed in the
C interval (0,1).
C X = Pseudo-random number in the range 0 to 3 whose frac-
C tional part is U.
C
IX = MOD(171*IX,30269)
IY = MOD(172*IY,30307)
IZ = MOD(170*IZ,30323)
X = (REAL(IX)/30269.) + (REAL(IY)/30307.) +
. (REAL(IZ)/30323.)
U = X - INT(X)
JRAND = REAL(N)*U + 1.
RETURN
END
LOGICAL FUNCTION LEFT (X1,Y1,X2,Y2,X0,Y0)
REAL X1, Y1, X2, Y2, X0, Y0
C
C***********************************************************
C
C From TRIPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C renka@cs.unt.edu
C 09/01/88
C
C This function determines whether node N0 is to the left
C or to the right of the line through N1-N2 as viewed by an
C observer at N1 facing N2.
C
C
C On input:
C
C X1,Y1 = Coordinates of N1.
C
C X2,Y2 = Coordinates of N2.
C
C X0,Y0 = Coordinates of N0.
C
C Input parameters are not altered by this function.
C
C On output:
C
C LEFT = .TRUE. if and only if (X0,Y0) is on or to the
C left of the directed line N1->N2.
C
C Modules required by LEFT: None
C
C***********************************************************
C
REAL DX1, DY1, DX2, DY2
C
C Local parameters:
C
C DX1,DY1 = X,Y components of the vector N1->N2
C DX2,DY2 = X,Y components of the vector N1->N0
C
DX1 = X2-X1
DY1 = Y2-Y1
DX2 = X0-X1
DY2 = Y0-Y1
C
C If the sign of the vector cross product of N1->N2 and
C N1->N0 is positive, then sin(A) > 0, where A is the
C angle between the vectors, and thus A is in the range
C (0,180) degrees.
C
LEFT = DX1*DY2 .GE. DX2*DY1
RETURN
END
INTEGER FUNCTION LSTPTR (LPL,NB,LIST,LPTR)
INTEGER LPL, NB, LIST(*), LPTR(*)
C
C***********************************************************
C
C From TRIPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C renka@cs.unt.edu
C 09/01/88
C
C This function returns the index (LIST pointer) of NB in
C the adjacency list for N0, where LPL = LEND(N0).
C
C
C On input:
C
C LPL = LEND(N0)
C
C NB = Index of the node whose pointer is to be re-
C turned. NB must be connected to N0.
C
C LIST,LPTR = Data structure defining the triangula-
C tion. Refer to Subroutine TRMESH.
C
C Input parameters are not altered by this function.
C
C On output:
C
C LSTPTR = Pointer such that LIST(LSTPTR) = NB or
C LIST(LSTPTR) = -NB, unless NB is not a
C neighbor of N0, in which case LSTPTR = LPL.
C
C Modules required by LSTPTR: None
C
C***********************************************************
C
INTEGER LP, ND
C
LP = LPTR(LPL)
1 ND = LIST(LP)
IF (ND .EQ. NB) GO TO 2
LP = LPTR(LP)
IF (LP .NE. LPL) GO TO 1
C
2 LSTPTR = LP
RETURN
END
INTEGER FUNCTION NBCNT (LPL,LPTR)
INTEGER LPL, LPTR(*)
C
C***********************************************************
C
C From TRIPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C renka@cs.unt.edu
C 09/01/88
C
C This function returns the number of neighbors of a node
C N0 in a triangulation created by Subroutine TRMESH (or
C TRMSHR).
C
C
C On input:
C
C LPL = LIST pointer to the last neighbor of N0 --
C LPL = LEND(N0).
C
C LPTR = Array of pointers associated with LIST.
C
C Input parameters are not altered by this function.
C
C On output:
C
C NBCNT = Number of neighbors of N0.
C
C Modules required by NBCNT: None
C
C***********************************************************
C
INTEGER K, LP
C
LP = LPL
K = 1
C
1 LP = LPTR(LP)
IF (LP .EQ. LPL) GO TO 2
K = K + 1
GO TO 1
C
2 NBCNT = K
RETURN
END
INTEGER FUNCTION NEARND (XP,YP,IST,N,X,Y,LIST,LPTR,
. LEND, DSQ)
INTEGER IST, N, LIST(*), LPTR(*), LEND(N)
REAL XP, YP, X(N), Y(N), DSQ
C
C***********************************************************
C
C From TRIPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C renka@cs.unt.edu
C 06/27/98
C
C Given a point P in the plane and a Delaunay triangula-
C tion created by Subroutine TRMESH or TRMSHR, this function
C returns the index of the nearest triangulation node to P.
C
C The algorithm consists of implicitly adding P to the
C triangulation, finding the nearest neighbor to P, and
C implicitly deleting P from the triangulation. Thus, it
C is based on the fact that, if P is a node in a Delaunay
C triangulation, the nearest node to P is a neighbor of P.
C
C
C On input:
C
C XP,YP = Cartesian coordinates of the point P to be
C located relative to the triangulation.
C
C IST = Index of a node at which TRFIND begins the
C search. Search time depends on the proximity
C of this node to P.
C
C N = Number of nodes in the triangulation. N .GE. 3.
C
C X,Y = Arrays of length N containing the Cartesian
C coordinates of the nodes.
C
C LIST,LPTR,LEND = Data structure defining the trian-
C gulation. Refer to TRMESH.
C
C Input parameters are not altered by this function.
C
C On output:
C
C NEARND = Nodal index of the nearest node to P, or 0
C if N < 3 or the triangulation data struc-
C ture is invalid.
C
C DSQ = Squared distance between P and NEARND unless
C NEARND = 0.
C
C Note that the number of candidates for NEARND
C (neighbors of P) is limited to LMAX defined in
C the PARAMETER statement below.
C
C Modules required by NEARND: JRAND, LEFT, LSTPTR, TRFIND
C
C Intrinsic function called by NEARND: ABS
C
C***********************************************************
C
INTEGER LSTPTR
INTEGER LMAX
PARAMETER (LMAX=25)
INTEGER I1, I2, I3, L, LISTP(LMAX), LP, LP1, LP2,
. LPL, LPTRP(LMAX), N1, N2, N3, NR, NST
REAL COS1, COS2, DS1, DSR, DX11, DX12, DX21,
. DX22, DY11, DY12, DY21, DY22, SIN1, SIN2
C
C Store local parameters and test for N invalid.
C
IF (N .LT. 3) GO TO 7
NST = IST
IF (NST .LT. 1 .OR. NST .GT. N) NST = 1
C
C Find a triangle (I1,I2,I3) containing P, or the rightmost
C (I1) and leftmost (I2) visible boundary nodes as viewed
C from P.
C
CALL TRFIND (NST,XP,YP,N,X,Y,LIST,LPTR,LEND, I1,I2,I3)
C
C Test for collinear nodes.
C
IF (I1 .EQ. 0) GO TO 7
C
C Store the linked list of 'neighbors' of P in LISTP and
C LPTRP. I1 is the first neighbor, and 0 is stored as
C the last neighbor if P is not contained in a triangle.
C L is the length of LISTP and LPTRP, and is limited to
C LMAX.
C
IF (I3 .NE. 0) THEN
LISTP(1) = I1
LPTRP(1) = 2
LISTP(2) = I2
LPTRP(2) = 3
LISTP(3) = I3
LPTRP(3) = 1
L = 3
ELSE
N1 = I1
L = 1
LP1 = 2
LISTP(L) = N1
LPTRP(L) = LP1
C
C Loop on the ordered sequence of visible boundary nodes
C N1 from I1 to I2.
C
1 LPL = LEND(N1)
N1 = -LIST(LPL)
L = LP1
LP1 = L+1
LISTP(L) = N1
LPTRP(L) = LP1
IF (N1 .NE. I2 .AND. LP1 .LT. LMAX) GO TO 1
L = LP1
LISTP(L) = 0
LPTRP(L) = 1
ENDIF
C
C Initialize variables for a loop on arcs N1-N2 opposite P
C in which new 'neighbors' are 'swapped' in. N1 follows
C N2 as a neighbor of P, and LP1 and LP2 are the LISTP
C indexes of N1 and N2.
C
LP2 = 1
N2 = I1
LP1 = LPTRP(1)
N1 = LISTP(LP1)
C
C Begin loop: find the node N3 opposite N1->N2.
C
2 LP = LSTPTR(LEND(N1),N2,LIST,LPTR)
IF (LIST(LP) .LT. 0) GO TO 4
LP = LPTR(LP)
N3 = ABS(LIST(LP))
C
C Swap test: Exit the loop if L = LMAX.
C
IF (L .EQ. LMAX) GO TO 5
DX11 = X(N1) - X(N3)
DX12 = X(N2) - X(N3)
DX22 = X(N2) - XP
DX21 = X(N1) - XP
C
DY11 = Y(N1) - Y(N3)
DY12 = Y(N2) - Y(N3)
DY22 = Y(N2) - YP
DY21 = Y(N1) - YP
C
COS1 = DX11*DX12 + DY11*DY12
COS2 = DX22*DX21 + DY22*DY21
IF (COS1 .GE. 0. .AND. COS2 .GE. 0.) GO TO 4
IF (COS1 .LT. 0. .AND. COS2 .LT. 0.) GO TO 3
C
SIN1 = DX11*DY12 - DX12*DY11
SIN2 = DX22*DY21 - DX21*DY22
IF (SIN1*COS2 + COS1*SIN2 .GE. 0.) GO TO 4
C
C Swap: Insert N3 following N2 in the adjacency list for P.
C The two new arcs opposite P must be tested.
C
3 L = L+1
LPTRP(LP2) = L
LISTP(L) = N3
LPTRP(L) = LP1
LP1 = L
N1 = N3
GO TO 2
C
C No swap: Advance to the next arc and test for termination
C on N1 = I1 (LP1 = 1) or N1 followed by 0.
C
4 IF (LP1 .EQ. 1) GO TO 5
LP2 = LP1
N2 = N1
LP1 = LPTRP(LP1)
N1 = LISTP(LP1)
IF (N1 .EQ. 0) GO TO 5
GO TO 2
C
C Set NR and DSR to the index of the nearest node to P and
C its squared distance from P, respectively.
C
5 NR = I1
DSR = (X(NR)-XP)**2 + (Y(NR)-YP)**2
DO 6 LP = 2,L
N1 = LISTP(LP)
IF (N1 .EQ. 0) GO TO 6
DS1 = (X(N1)-XP)**2 + (Y(N1)-YP)**2
IF (DS1 .LT. DSR) THEN
NR = N1
DSR = DS1
ENDIF
6 CONTINUE
DSQ = DSR
NEARND = NR
RETURN
C
C Invalid input.
C
7 NEARND = 0
RETURN
END
SUBROUTINE OPTIM (X,Y,NA, LIST,LPTR,LEND,NIT,IWK, IER)
INTEGER NA, LIST(*), LPTR(*), LEND(*), NIT, IWK(2,NA),
. IER
REAL X(*), Y(*)
C
C***********************************************************
C
C From TRIPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C renka@cs.unt.edu
C 06/27/98
C
C Given a set of NA triangulation arcs, this subroutine
C optimizes the portion of the triangulation consisting of
C the quadrilaterals (pairs of adjacent triangles) which
C have the arcs as diagonals by applying the circumcircle
C test and appropriate swaps to the arcs.
C
C An iteration consists of applying the swap test and
C swaps to all NA arcs in the order in which they are
C stored. The iteration is repeated until no swap occurs
C or NIT iterations have been performed. The bound on the
C number of iterations may be necessary to prevent an
C infinite loop caused by cycling (reversing the effect of a
C previous swap) due to floating point inaccuracy when four
C or more nodes are nearly cocircular.
C
C
C On input:
C
C X,Y = Arrays containing the nodal coordinates.
C
C NA = Number of arcs in the set. NA .GE. 0.
C
C The above parameters are not altered by this routine.
C
C LIST,LPTR,LEND = Data structure defining the trian-
C gulation. Refer to Subroutine
C TRMESH.
C
C NIT = Maximum number of iterations to be performed.
C A reasonable value is 3*NA. NIT .GE. 1.
C
C IWK = Integer array dimensioned 2 by NA containing
C the nodal indexes of the arc endpoints (pairs
C of endpoints are stored in columns).
C
C On output:
C
C LIST,LPTR,LEND = Updated triangulation data struc-
C ture reflecting the swaps.
C
C NIT = Number of iterations performed.
C
C IWK = Endpoint indexes of the new set of arcs
C reflecting the swaps.
C
C IER = Error indicator:
C IER = 0 if no errors were encountered.
C IER = 1 if a swap occurred on the last of
C MAXIT iterations, where MAXIT is the
C value of NIT on input. The new set
C of arcs in not necessarily optimal
C in this case.
C IER = 2 if NA < 0 or NIT < 1 on input.
C IER = 3 if IWK(2,I) is not a neighbor of
C IWK(1,I) for some I in the range 1
C to NA. A swap may have occurred in
C this case.
C IER = 4 if a zero pointer was returned by
C Subroutine SWAP.
C
C Modules required by OPTIM: LSTPTR, SWAP, SWPTST
C
C Intrinsic function called by OPTIM: ABS
C
C***********************************************************
C
LOGICAL SWPTST
INTEGER I, IO1, IO2, ITER, LP, LP21, LPL, LPP, MAXIT,
. N1, N2, NNA
LOGICAL SWP
C
C Local parameters:
C
C I = Column index for IWK
C IO1,IO2 = Nodal indexes of the endpoints of an arc in IWK
C ITER = Iteration count
C LP = LIST pointer
C LP21 = Parameter returned by SWAP (not used)
C LPL = Pointer to the last neighbor of IO1
C LPP = Pointer to the node preceding IO2 as a neighbor
C of IO1
C MAXIT = Input value of NIT
C N1,N2 = Nodes opposite IO1->IO2 and IO2->IO1,
C respectively
C NNA = Local copy of NA
C SWP = Flag set to TRUE iff a swap occurs in the
C optimization loop
C
NNA = NA
MAXIT = NIT
IF (NNA .LT. 0 .OR. MAXIT .LT. 1) GO TO 7
C
C Initialize iteration count ITER and test for NA = 0.
C
ITER = 0
IF (NNA .EQ. 0) GO TO 5
C
C Top of loop --
C SWP = TRUE iff a swap occurred in the current iteration.
C
1 IF (ITER .EQ. MAXIT) GO TO 6
ITER = ITER + 1
SWP = .FALSE.
C
C Inner loop on arcs IO1-IO2 --
C
DO 4 I = 1,NNA
IO1 = IWK(1,I)
IO2 = IWK(2,I)
C
C Set N1 and N2 to the nodes opposite IO1->IO2 and
C IO2->IO1, respectively. Determine the following:
C
C LPL = pointer to the last neighbor of IO1,
C LP = pointer to IO2 as a neighbor of IO1, and
C LPP = pointer to the node N2 preceding IO2.
C
LPL = LEND(IO1)
LPP = LPL
LP = LPTR(LPP)
2 IF (LIST(LP) .EQ. IO2) GO TO 3
LPP = LP
LP = LPTR(LPP)
IF (LP .NE. LPL) GO TO 2
C
C IO2 should be the last neighbor of IO1. Test for no
C arc and bypass the swap test if IO1 is a boundary
C node.
C
IF (ABS(LIST(LP)) .NE. IO2) GO TO 8
IF (LIST(LP) .LT. 0) GO TO 4
C
C Store N1 and N2, or bypass the swap test if IO1 is a
C boundary node and IO2 is its first neighbor.
C
3 N2 = LIST(LPP)
IF (N2 .LT. 0) GO TO 4
LP = LPTR(LP)
N1 = ABS(LIST(LP))
C
C Test IO1-IO2 for a swap, and update IWK if necessary.
C
IF ( .NOT. SWPTST(N1,N2,IO1,IO2,X,Y) ) GO TO 4
CALL SWAP (N1,N2,IO1,IO2, LIST,LPTR,LEND, LP21)
IF (LP21 .EQ. 0) GO TO 9
SWP = .TRUE.
IWK(1,I) = N1
IWK(2,I) = N2
4 CONTINUE
IF (SWP) GO TO 1
C
C Successful termination.
C
5 NIT = ITER
IER = 0
RETURN
C
C MAXIT iterations performed without convergence.
C
6 NIT = MAXIT
IER = 1
RETURN
C
C Invalid input parameter.
C
7 NIT = 0
IER = 2
RETURN
C
C IO2 is not a neighbor of IO1.
C
8 NIT = ITER
IER = 3
RETURN
C
C Zero pointer returned by SWAP.
C
9 NIT = ITER
IER = 4
RETURN
END
REAL FUNCTION STORE (X)
REAL X
C
C***********************************************************
C
C From TRIPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C renka@cs.unt.edu
C 03/18/90
C
C This function forces its argument X to be stored in a
C memory location, thus providing a means of determining
C floating point number characteristics (such as the machine
C precision) when it is necessary to avoid computation in
C high precision registers.
C
C
C On input:
C
C X = Value to be stored.
C
C X is not altered by this function.
C
C On output:
C
C STORE = Value of X after it has been stored and
C possibly truncated or rounded to the single
C precision word length.
C
C Modules required by STORE: None
C
C***********************************************************
C
REAL Y
COMMON/STCOM/Y
C
Y = X
STORE = Y
RETURN
END
SUBROUTINE SWAP (IN1,IN2,IO1,IO2, LIST,LPTR,
. LEND, LP21)
INTEGER IN1, IN2, IO1, IO2, LIST(*), LPTR(*), LEND(*),
. LP21
C
C***********************************************************
C
C From TRIPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C renka@cs.unt.edu
C 06/22/98
C
C Given a triangulation of a set of points on the unit
C sphere, this subroutine replaces a diagonal arc in a
C strictly convex quadrilateral (defined by a pair of adja-
C cent triangles) with the other diagonal. Equivalently, a
C pair of adjacent triangles is replaced by another pair
C having the same union.
C
C
C On input:
C
C IN1,IN2,IO1,IO2 = Nodal indexes of the vertices of
C the quadrilateral. IO1-IO2 is re-
C placed by IN1-IN2. (IO1,IO2,IN1)
C and (IO2,IO1,IN2) must be trian-
C gles on input.
C
C The above parameters are not altered by this routine.
C
C LIST,LPTR,LEND = Data structure defining the trian-
C gulation. Refer to Subroutine
C TRMESH.
C
C On output:
C
C LIST,LPTR,LEND = Data structure updated with the
C swap -- triangles (IO1,IO2,IN1) and
C (IO2,IO1,IN2) are replaced by
C (IN1,IN2,IO2) and (IN2,IN1,IO1)
C unless LP21 = 0.
C
C LP21 = Index of IN1 as a neighbor of IN2 after the
C swap is performed unless IN1 and IN2 are
C adjacent on input, in which case LP21 = 0.
C
C Module required by SWAP: LSTPTR
C
C Intrinsic function called by SWAP: ABS
C
C***********************************************************
C
INTEGER LSTPTR
INTEGER LP, LPH, LPSAV
C
C Local parameters:
C
C LP,LPH,LPSAV = LIST pointers
C
C
C Test for IN1 and IN2 adjacent.
C
LP = LSTPTR(LEND(IN1),IN2,LIST,LPTR)
IF (ABS(LIST(LP)) .EQ. IN2) THEN
LP21 = 0
RETURN
ENDIF
C
C Delete IO2 as a neighbor of IO1.
C
LP = LSTPTR(LEND(IO1),IN2,LIST,LPTR)
LPH = LPTR(LP)
LPTR(LP) = LPTR(LPH)
C
C If IO2 is the last neighbor of IO1, make IN2 the
C last neighbor.
C
IF (LEND(IO1) .EQ. LPH) LEND(IO1) = LP
C
C Insert IN2 as a neighbor of IN1 following IO1
C using the hole created above.
C
LP = LSTPTR(LEND(IN1),IO1,LIST,LPTR)
LPSAV = LPTR(LP)
LPTR(LP) = LPH
LIST(LPH) = IN2
LPTR(LPH) = LPSAV
C
C Delete IO1 as a neighbor of IO2.
C
LP = LSTPTR(LEND(IO2),IN1,LIST,LPTR)
LPH = LPTR(LP)
LPTR(LP) = LPTR(LPH)
C
C If IO1 is the last neighbor of IO2, make IN1 the
C last neighbor.
C
IF (LEND(IO2) .EQ. LPH) LEND(IO2) = LP
C
C Insert IN1 as a neighbor of IN2 following IO2.
C
LP = LSTPTR(LEND(IN2),IO2,LIST,LPTR)
LPSAV = LPTR(LP)
LPTR(LP) = LPH
LIST(LPH) = IN1
LPTR(LPH) = LPSAV
LP21 = LPH
RETURN
END
LOGICAL FUNCTION SWPTST (IN1,IN2,IO1,IO2,X,Y)
INTEGER IN1, IN2, IO1, IO2
REAL X(*), Y(*)
C
C***********************************************************
C
C From TRIPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C renka@cs.unt.edu
C 09/01/88
C
C This function applies the circumcircle test to a quadri-
C lateral defined by a pair of adjacent triangles. The
C diagonal arc (shared triangle side) should be swapped for
C the other diagonl if and only if the fourth vertex is
C strictly interior to the circumcircle of one of the
C triangles (the decision is independent of the choice of
C triangle). Equivalently, the diagonal is chosen to maxi-
C mize the smallest of the six interior angles over the two
C pairs of possible triangles (the decision is for no swap
C if the quadrilateral is not strictly convex).
C
C When the four vertices are nearly cocircular (the
C neutral case), the preferred decision is no swap -- in
C order to avoid unnecessary swaps and, more important, to
C avoid cycling in Subroutine OPTIM which is called by
C DELNOD and EDGE. Thus, a tolerance SWTOL (stored in
C SWPCOM by TRMESH or TRMSHR) is used to define 'nearness'
C to the neutral case.
C
C
C On input:
C
C IN1,IN2,IO1,IO2 = Nodal indexes of the vertices of
C the quadrilateral. IO1-IO2 is the
C triangulation arc (shared triangle
C side) to be replaced by IN1-IN2 if
C the decision is to swap. The
C triples (IO1,IO2,IN1) and (IO2,
C IO1,IN2) must define triangles (be
C in counterclockwise order) on in-
C put.
C
C X,Y = Arrays containing the nodal coordinates.
C
C Input parameters are not altered by this routine.
C
C On output:
C
C SWPTST = .TRUE. if and only if the arc connecting
C IO1 and IO2 is to be replaced.
C
C Modules required by SWPTST: None
C
C***********************************************************
C
REAL DX11, DX12, DX22, DX21, DY11, DY12, DY22, DY21,
. SIN1, SIN2, COS1, COS2, SIN12, SWTOL
C
C Tolerance stored by TRMESH or TRMSHR.
C
COMMON/SWPCOM/SWTOL
C
C Local parameters:
C
C DX11,DY11 = X,Y components of the vector IN1->IO1
C DX12,DY12 = X,Y components of the vector IN1->IO2
C DX22,DY22 = X,Y components of the vector IN2->IO2
C DX21,DY21 = X,Y components of the vector IN2->IO1
C SIN1 = Cross product of the vectors IN1->IO1 and
C IN1->IO2 -- proportional to sin(T1), where
C T1 is the angle at IN1 formed by the vectors
C COS1 = Inner product of the vectors IN1->IO1 and
C IN1->IO2 -- proportional to cos(T1)
C SIN2 = Cross product of the vectors IN2->IO2 and
C IN2->IO1 -- proportional to sin(T2), where
C T2 is the angle at IN2 formed by the vectors
C COS2 = Inner product of the vectors IN2->IO2 and
C IN2->IO1 -- proportional to cos(T2)
C SIN12 = SIN1*COS2 + COS1*SIN2 -- proportional to
C sin(T1+T2)
C
C
C Compute the vectors containing the angles T1 and T2.
C
DX11 = X(IO1) - X(IN1)
DX12 = X(IO2) - X(IN1)
DX22 = X(IO2) - X(IN2)
DX21 = X(IO1) - X(IN2)
C
DY11 = Y(IO1) - Y(IN1)
DY12 = Y(IO2) - Y(IN1)
DY22 = Y(IO2) - Y(IN2)
DY21 = Y(IO1) - Y(IN2)
C
C Compute inner products.
C
COS1 = DX11*DX12 + DY11*DY12
COS2 = DX22*DX21 + DY22*DY21
C
C The diagonals should be swapped iff (T1+T2) > 180
C degrees. The following two tests ensure numerical
C stability: the decision must be FALSE when both
C angles are close to 0, and TRUE when both angles
C are close to 180 degrees.
C
IF (COS1 .GE. 0. .AND. COS2 .GE. 0.) GO TO 2
IF (COS1 .LT. 0. .AND. COS2 .LT. 0.) GO TO 1
C
C Compute vector cross products (Z-components).
C
SIN1 = DX11*DY12 - DX12*DY11
SIN2 = DX22*DY21 - DX21*DY22
SIN12 = SIN1*COS2 + COS1*SIN2
IF (SIN12 .GE. -SWTOL) GO TO 2
C
C Swap.
C
1 SWPTST = .TRUE.
RETURN
C
C No swap.
C
2 SWPTST = .FALSE.
RETURN
END
SUBROUTINE TRFIND (NST,PX,PY,N,X,Y,LIST,LPTR,LEND, I1,
. I2,I3)
INTEGER NST, N, LIST(*), LPTR(*), LEND(N), I1, I2, I3
REAL PX, PY, X(N), Y(N)
C
C***********************************************************
C
C From TRIPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C renka@cs.unt.edu
C 07/28/98
C
C This subroutine locates a point P relative to a triangu-
C lation created by Subroutine TRMESH or TRMSHR. If P is
C contained in a triangle, the three vertex indexes are
C returned. Otherwise, the indexes of the rightmost and
C leftmost visible boundary nodes are returned.
C
C
C On input:
C
C NST = Index of a node at which TRFIND begins the
C search. Search time depends on the proximity
C of this node to P.
C
C PX,PY = X and y coordinates of the point P to be
C located.
C
C N = Number of nodes in the triangulation. N .GE. 3.
C
C X,Y = Arrays of length N containing the coordinates
C of the nodes in the triangulation.
C
C LIST,LPTR,LEND = Data structure defining the trian-
C gulation. Refer to Subroutine
C TRMESH.
C
C Input parameters are not altered by this routine.
C
C On output:
C
C I1,I2,I3 = Nodal indexes, in counterclockwise order,
C of the vertices of a triangle containing
C P if P is contained in a triangle. If P
C is not in the convex hull of the nodes,
C I1 indexes the rightmost visible boundary
C node, I2 indexes the leftmost visible
C boundary node, and I3 = 0. Rightmost and
C leftmost are defined from the perspective
C of P, and a pair of points are visible
C from each other if and only if the line
C segment joining them intersects no trian-
C gulation arc. If P and all of the nodes
C lie on a common line, then I1 = I2 = I3 =
C 0 on output.
C
C Modules required by TRFIND: JRAND, LEFT, LSTPTR, STORE
C
C Intrinsic function called by TRFIND: ABS
C
C***********************************************************
C
INTEGER JRAND, LSTPTR
LOGICAL LEFT
REAL STORE
INTEGER IX, IY, IZ, LP, N0, N1, N1S, N2, N2S, N3, N4,
. NB, NF, NL, NP, NPP
LOGICAL FRWRD
REAL B1, B2, XA, XB, XC, XP, YA, YB, YC, YP
C
SAVE IX, IY, IZ
DATA IX/1/, IY/2/, IZ/3/
C
C Local parameters:
C
C B1,B2 = Unnormalized barycentric coordinates of P with
C respect to (N1,N2,N3)
C IX,IY,IZ = Integer seeds for JRAND
C LP = LIST pointer
C N0,N1,N2 = Nodes in counterclockwise order defining a
C cone (with vertex N0) containing P
C N1S,N2S = Saved values of N1 and N2
C N3,N4 = Nodes opposite N1->N2 and N2->N1, respectively
C NB = Index of a boundary node -- first neighbor of
C NF or last neighbor of NL in the boundary
C traversal loops
C NF,NL = First and last neighbors of N0, or first
C (rightmost) and last (leftmost) nodes
C visible from P when P is exterior to the
C triangulation
C NP,NPP = Indexes of boundary nodes used in the boundary
C traversal loops
C XA,XB,XC = Dummy arguments for FRWRD
C YA,YB,YC = Dummy arguments for FRWRD
C XP,YP = Local variables containing the components of P
C
C Statement function:
C
C FRWRD = TRUE iff C is forward of A->B
C iff B,A->C> .GE. 0.
C
FRWRD(XA,YA,XB,YB,XC,YC) = (XB-XA)*(XC-XA) +
. (YB-YA)*(YC-YA) .GE. 0.
C
C Initialize variables.
C
XP = PX
YP = PY
N0 = NST
IF (N0 .LT. 1 .OR. N0 .GT. N)
. N0 = JRAND(N, IX,IY,IZ )
C
C Set NF and NL to the first and last neighbors of N0, and
C initialize N1 = NF.
C
1 LP = LEND(N0)
NL = LIST(LP)
LP = LPTR(LP)
NF = LIST(LP)
N1 = NF
C
C Find a pair of adjacent neighbors N1,N2 of N0 that define
C a wedge containing P: P LEFT N0->N1 and P RIGHT N0->N2.
C
IF (NL .GT. 0) GO TO 2
C
C N0 is a boundary node. Test for P exterior.
C
NL = -NL
IF ( .NOT. LEFT(X(N0),Y(N0),X(NF),Y(NF),XP,YP) ) THEN
NL = N0
GO TO 9
ENDIF
IF ( .NOT. LEFT(X(NL),Y(NL),X(N0),Y(N0),XP,YP) ) THEN
NB = NF
NF = N0
NP = NL
NPP = N0
GO TO 11
ENDIF
GO TO 3
C
C N0 is an interior node. Find N1.
C
2 IF ( LEFT(X(N0),Y(N0),X(N1),Y(N1),XP,YP) ) GO TO 3
LP = LPTR(LP)
N1 = LIST(LP)
IF (N1 .EQ. NL) GO TO 6
GO TO 2
C
C P is to the left of edge N0->N1. Initialize N2 to the
C next neighbor of N0.
C
3 LP = LPTR(LP)
N2 = ABS(LIST(LP))
IF ( .NOT. LEFT(X(N0),Y(N0),X(N2),Y(N2),XP,YP) )
. GO TO 7
N1 = N2
IF (N1 .NE. NL) GO TO 3
IF ( .NOT. LEFT(X(N0),Y(N0),X(NF),Y(NF),XP,YP) )
. GO TO 6
IF (XP .EQ. X(N0) .AND. YP .EQ. Y(N0)) GO TO 5
C
C P is left of or on edges N0->NB for all neighbors NB
C of N0.
C All points are collinear iff P is left of NB->N0 for
C all neighbors NB of N0. Search the neighbors of N0.
C NOTE -- N1 = NL and LP points to NL.
C
4 IF ( .NOT. LEFT(X(N1),Y(N1),X(N0),Y(N0),XP,YP) )
. GO TO 5
LP = LPTR(LP)
N1 = ABS(LIST(LP))
IF (N1 .EQ. NL) GO TO 17
GO TO 4
C
C P is to the right of N1->N0, or P=N0. Set N0 to N1 and
C start over.
C
5 N0 = N1
GO TO 1
C
C P is between edges N0->N1 and N0->NF.
C
6 N2 = NF
C
C P is contained in the wedge defined by line segments
C N0->N1 and N0->N2, where N1 is adjacent to N2. Set
C N3 to the node opposite N1->N2, and save N1 and N2 to
C test for cycling.
C
7 N3 = N0
N1S = N1
N2S = N2
C
C Top of edge hopping loop. Test for termination.
C
8 IF ( LEFT(X(N1),Y(N1),X(N2),Y(N2),XP,YP) ) THEN
C
C P LEFT N1->N2 and hence P is in (N1,N2,N3) unless an
C error resulted from floating point inaccuracy and
C collinearity. Compute the unnormalized barycentric
C coordinates of P with respect to (N1,N2,N3).
C
B1 = (X(N3)-X(N2))*(YP-Y(N2)) -
. (XP-X(N2))*(Y(N3)-Y(N2))
B2 = (X(N1)-X(N3))*(YP-Y(N3)) -
. (XP-X(N3))*(Y(N1)-Y(N3))
IF (STORE(B1+1.) .GE. 1. .AND.
. STORE(B2+1.) .GE. 1.) GO TO 16
C
C Restart with N0 randomly selected.
C
N0 = JRAND(N, IX,IY,IZ )
GO TO 1
ENDIF
C
C Set N4 to the neighbor of N2 which follows N1 (node
C opposite N2->N1) unless N1->N2 is a boundary edge.
C
LP = LSTPTR(LEND(N2),N1,LIST,LPTR)
IF (LIST(LP) .LT. 0) THEN
NF = N2
NL = N1
GO TO 9
ENDIF
LP = LPTR(LP)
N4 = ABS(LIST(LP))
C
C Select the new edge N1->N2 which intersects the line
C segment N0-P, and set N3 to the node opposite N1->N2.
C
IF ( LEFT(X(N0),Y(N0),X(N4),Y(N4),XP,YP) ) THEN
N3 = N1
N1 = N4
N2S = N2
IF (N1 .NE. N1S .AND. N1 .NE. N0) GO TO 8
ELSE
N3 = N2
N2 = N4
N1S = N1
IF (N2 .NE. N2S .AND. N2 .NE. N0) GO TO 8
ENDIF
C
C The starting node N0 or edge N1-N2 was encountered
C again, implying a cycle (infinite loop). Restart
C with N0 randomly selected.
C
N0 = JRAND(N, IX,IY,IZ )
GO TO 1
C
C Boundary traversal loops. NL->NF is a boundary edge and
C P RIGHT NL->NF. Save NL and NF.
9 NP = NL
NPP = NF
C
C Find the first (rightmost) visible boundary node NF. NB
C is set to the first neighbor of NF, and NP is the last
C neighbor.
C
10 LP = LEND(NF)
LP = LPTR(LP)
NB = LIST(LP)
IF ( .NOT. LEFT(X(NF),Y(NF),X(NB),Y(NB),XP,YP) )
. GO TO 12
C
C P LEFT NF->NB and thus NB is not visible unless an error
C resulted from floating point inaccuracy and collinear-
C ity of the 4 points NP, NF, NB, and P.
C
11 IF ( FRWRD(X(NF),Y(NF),X(NP),Y(NP),XP,YP) .OR.
. FRWRD(X(NF),Y(NF),X(NP),Y(NP),X(NB),Y(NB)) ) THEN
I1 = NF
GO TO 13
ENDIF
C
C Bottom of loop.
C
12 NP = NF
NF = NB
GO TO 10
C
C Find the last (leftmost) visible boundary node NL. NB
C is set to the last neighbor of NL, and NPP is the first
C neighbor.
C
13 LP = LEND(NL)
NB = -LIST(LP)
IF ( .NOT. LEFT(X(NB),Y(NB),X(NL),Y(NL),XP,YP) )
. GO TO 14
C
C P LEFT NB->NL and thus NB is not visible unless an error
C resulted from floating point inaccuracy and collinear-
C ity of the 4 points P, NB, NL, and NPP.
C
IF ( FRWRD(X(NL),Y(NL),X(NPP),Y(NPP),XP,YP) .OR.
. FRWRD(X(NL),Y(NL),X(NPP),Y(NPP),X(NB),Y(NB)) )
. GO TO 15
C
C Bottom of loop.
C
14 NPP = NL
NL = NB
GO TO 13
C
C NL is the leftmost visible boundary node.
C
15 I2 = NL
I3 = 0
RETURN
C
C P is in the triangle (N1,N2,N3).
C
16 I1 = N1
I2 = N2
I3 = N3
RETURN
C
C All points are collinear.
C
17 I1 = 0
I2 = 0
I3 = 0
RETURN
END
SUBROUTINE TRLIST (NCC,LCC,N,LIST,LPTR,LEND,NROW, NT,
. LTRI,LCT,IER)
INTEGER NCC, LCC(*), N, LIST(*), LPTR(*), LEND(N),
. NROW, NT, LTRI(NROW,*), LCT(*), IER
C
C***********************************************************
C
C From TRIPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C renka@cs.unt.edu
C 03/22/97
C
C This subroutine converts a triangulation data structure
C from the linked list created by Subroutine TRMESH or
C TRMSHR to a triangle list.
C
C On input:
C
C NCC = Number of constraints. NCC .GE. 0.
C
C LCC = List of constraint curve starting indexes (or
C dummy array of length 1 if NCC = 0). Refer to
C Subroutine ADDCST.
C
C N = Number of nodes in the triangulation. N .GE. 3.
C
C LIST,LPTR,LEND = Linked list data structure defin-
C ing the triangulation. Refer to
C Subroutine TRMESH.
C
C NROW = Number of rows (entries per triangle) re-
C served for the triangle list LTRI. The value
C must be 6 if only the vertex indexes and
C neighboring triangle indexes are to be
C stored, or 9 if arc indexes are also to be
C assigned and stored. Refer to LTRI.
C
C The above parameters are not altered by this routine.
C
C LTRI = Integer array of length at least NROW*NT,
C where NT is at most 2N-5. (A sufficient
C length is 12N if NROW=6 or 18N if NROW=9.)
C
C LCT = Integer array of length NCC or dummy array of
C length 1 if NCC = 0.
C
C On output:
C
C NT = Number of triangles in the triangulation unless
C IER .NE. 0, in which case NT = 0. NT = 2N - NB
C - 2, where NB is the number of boundary nodes.
C
C LTRI = NROW by NT array whose J-th column contains
C the vertex nodal indexes (first three rows),
C neighboring triangle indexes (second three
C rows), and, if NROW = 9, arc indexes (last
C three rows) associated with triangle J for
C J = 1,...,NT. The vertices are ordered
C counterclockwise with the first vertex taken
C to be the one with smallest index. Thus,
C LTRI(2,J) and LTRI(3,J) are larger than
C LTRI(1,J) and index adjacent neighbors of
C node LTRI(1,J). For I = 1,2,3, LTRI(I+3,J)
C and LTRI(I+6,J) index the triangle and arc,
C respectively, which are opposite (not shared
C by) node LTRI(I,J), with LTRI(I+3,J) = 0 if
C LTRI(I+6,J) indexes a boundary arc. Vertex
C indexes range from 1 to N, triangle indexes
C from 0 to NT, and, if included, arc indexes
C from 1 to NA = NT+N-1. The triangles are or-
C dered on first (smallest) vertex indexes,
C except that the sets of constraint triangles
C (triangles contained in the closure of a con-
C straint region) follow the non-constraint
C triangles.
C
C LCT = Array of length NCC containing the triangle
C index of the first triangle of constraint J in
C LCT(J). Thus, the number of non-constraint
C triangles is LCT(1)-1, and constraint J con-
C tains LCT(J+1)-LCT(J) triangles, where
C LCT(NCC+1) = NT+1.
C
C IER = Error indicator.
C IER = 0 if no errors were encountered.
C IER = 1 if NCC, N, NROW, or an LCC entry is
C outside its valid range on input.
C IER = 2 if the triangulation data structure
C (LIST,LPTR,LEND) is invalid. Note,
C however, that these arrays are not
C completely tested for validity.
C
C Modules required by TRLIST: None
C
C Intrinsic function called by TRLIST: ABS
C
C***********************************************************
C
INTEGER I, I1, I2, I3, ISV, J, JLAST, KA, KN, KT, L,
. LCC1, LP, LP2, LPL, LPLN1, N1, N1ST, N2, N3,
. NM2, NN
LOGICAL ARCS, CSTRI, PASS2
C
C Test for invalid input parameters and store the index
C LCC1 of the first constraint node (if any).
C
NN = N
IF (NCC .LT. 0 .OR. (NROW .NE. 6 .AND.
. NROW .NE. 9)) GO TO 12
LCC1 = NN+1
IF (NCC .EQ. 0) THEN
IF (NN .LT. 3) GO TO 12
ELSE
DO 1 I = NCC,1,-1
IF (LCC1-LCC(I) .LT. 3) GO TO 12
LCC1 = LCC(I)
1 CONTINUE
IF (LCC1 .LT. 1) GO TO 12
ENDIF
C
C Initialize parameters for loop on triangles KT = (N1,N2,
C N3), where N1 < N2 and N1 < N3. This requires two
C passes through the nodes with all non-constraint
C triangles stored on the first pass, and the constraint
C triangles stored on the second.
C
C ARCS = TRUE iff arc indexes are to be stored.
C KA,KT = Numbers of currently stored arcs and triangles.
C N1ST = Starting index for the loop on nodes (N1ST = 1 on
C pass 1, and N1ST = LCC1 on pass 2).
C NM2 = Upper bound on candidates for N1.
C PASS2 = TRUE iff constraint triangles are to be stored.
C
ARCS = NROW .EQ. 9
KA = 0
KT = 0
N1ST = 1
NM2 = NN-2
PASS2 = .FALSE.
C
C Loop on nodes N1: J = constraint containing N1,
C JLAST = last node in constraint J.
C
2 J = 0
JLAST = LCC1 - 1
DO 11 N1 = N1ST,NM2
IF (N1 .GT. JLAST) THEN
C
C N1 is the first node in constraint J+1. Update J and
C JLAST, and store the first constraint triangle index
C if in pass 2.
C
J = J + 1
IF (J .LT. NCC) THEN
JLAST = LCC(J+1) - 1
ELSE
JLAST = NN
ENDIF
IF (PASS2) LCT(J) = KT + 1
ENDIF
C
C Loop on pairs of adjacent neighbors (N2,N3). LPLN1 points
C to the last neighbor of N1, and LP2 points to N2.
C
LPLN1 = LEND(N1)
LP2 = LPLN1
3 LP2 = LPTR(LP2)
N2 = LIST(LP2)
LP = LPTR(LP2)
N3 = ABS(LIST(LP))
IF (N2 .LT. N1 .OR. N3 .LT. N1) GO TO 10
C
C (N1,N2,N3) is a constraint triangle iff the three nodes
C are in the same constraint and N2 < N3. Bypass con-
C straint triangles on pass 1 and non-constraint triangles
C on pass 2.
C
CSTRI = N1 .GE. LCC1 .AND. N2 .LT. N3 .AND.
. N3 .LE. JLAST
IF ((CSTRI .AND. .NOT. PASS2) .OR.
. (.NOT. CSTRI .AND. PASS2)) GO TO 10
C
C Add a new triangle KT = (N1,N2,N3).
C
KT = KT + 1
LTRI(1,KT) = N1
LTRI(2,KT) = N2
LTRI(3,KT) = N3
C
C Loop on triangle sides (I1,I2) with neighboring triangles
C KN = (I1,I2,I3).
C
DO 9 I = 1,3
IF (I .EQ. 1) THEN
I1 = N3
I2 = N2
ELSEIF (I .EQ. 2) THEN
I1 = N1
I2 = N3
ELSE
I1 = N2
I2 = N1
ENDIF
C
C Set I3 to the neighbor of I1 which follows I2 unless
C I2->I1 is a boundary arc.
C
LPL = LEND(I1)
LP = LPTR(LPL)
4 IF (LIST(LP) .EQ. I2) GO TO 5
LP = LPTR(LP)
IF (LP .NE. LPL) GO TO 4
C
C I2 is the last neighbor of I1 unless the data structure
C is invalid. Bypass the search for a neighboring
C triangle if I2->I1 is a boundary arc.
C
IF (ABS(LIST(LP)) .NE. I2) GO TO 13
KN = 0
IF (LIST(LP) .LT. 0) GO TO 8
C
C I2->I1 is not a boundary arc, and LP points to I2 as
C a neighbor of I1.
C
5 LP = LPTR(LP)
I3 = ABS(LIST(LP))
C
C Find L such that LTRI(L,KN) = I3 (not used if KN > KT),
C and permute the vertex indexes of KN so that I1 is
C smallest.
C
IF (I1 .LT. I2 .AND. I1 .LT. I3) THEN
L = 3
ELSEIF (I2 .LT. I3) THEN
L = 2
ISV = I1
I1 = I2
I2 = I3
I3 = ISV
ELSE
L = 1
ISV = I1
I1 = I3
I3 = I2
I2 = ISV
ENDIF
C
C Test for KN > KT (triangle index not yet assigned).
C
IF (I1 .GT. N1 .AND. .NOT. PASS2) GO TO 9
C
C Find KN, if it exists, by searching the triangle list in
C reverse order.
C
DO 6 KN = KT-1,1,-1
IF (LTRI(1,KN) .EQ. I1 .AND. LTRI(2,KN) .EQ.
. I2 .AND. LTRI(3,KN) .EQ. I3) GO TO 7
6 CONTINUE
GO TO 9
C
C Store KT as a neighbor of KN.
C
7 LTRI(L+3,KN) = KT
C
C Store KN as a neighbor of KT, and add a new arc KA.
C
8 LTRI(I+3,KT) = KN
IF (ARCS) THEN
KA = KA + 1
LTRI(I+6,KT) = KA
IF (KN .NE. 0) LTRI(L+6,KN) = KA
ENDIF
9 CONTINUE
C
C Bottom of loop on triangles.
C
10 IF (LP2 .NE. LPLN1) GO TO 3
11 CONTINUE
C
C Bottom of loop on nodes.
C
IF (.NOT. PASS2 .AND. NCC .GT. 0) THEN
PASS2 = .TRUE.
N1ST = LCC1
GO TO 2
ENDIF
C
C No errors encountered.
C
NT = KT
IER = 0
RETURN
C
C Invalid input parameter.
C
12 NT = 0
IER = 1
RETURN
C
C Invalid triangulation data structure: I1 is a neighbor of
C I2, but I2 is not a neighbor of I1.
C
13 NT = 0
IER = 2
RETURN
END
SUBROUTINE TRLPRT (NCC,LCT,N,X,Y,NROW,NT,LTRI,LOUT,
. PRNTX)
INTEGER NCC, LCT(*), N, NROW, NT, LTRI(NROW,NT),
. LOUT
LOGICAL PRNTX
REAL X(N), Y(N)
C
C***********************************************************
C
C From TRLPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C renka@cs.unt.edu
C 07/02/98
C
C Given a triangulation of a set of points in the plane,
C this subroutine prints the triangle list created by
C Subroutine TRLIST and, optionally, the nodal coordinates
C on logical unit LOUT. The numbers of boundary nodes,
C triangles, and arcs, and the constraint region triangle
C indexes, if any, are also printed.
C
C All parameters other than LOUT and PRNTX should be
C unaltered from their values on output from TRLIST.
C
C
C On input:
C
C NCC = Number of constraints.
C
C LCT = List of constraint triangle starting indexes
C (or dummy array of length 1 if NCC = 0).
C
C N = Number of nodes in the triangulation.
C 3 .LE. N .LE. 9999.
C
C X,Y = Arrays of length N containing the coordinates
C of the nodes in the triangulation -- not used
C unless PRNTX = TRUE.
C
C NROW = Number of rows (entries per triangle) re-
C served for the triangle list LTRI. The value
C must be 6 if only the vertex indexes and
C neighboring triangle indexes are stored, or 9
C if arc indexes are also stored.
C
C NT = Number of triangles in the triangulation.
C 1 .LE. NT .LE. 9999.
C
C LTRI = NROW by NT array whose J-th column contains
C the vertex nodal indexes (first three rows),
C neighboring triangle indexes (second three
C rows), and, if NROW = 9, arc indexes (last
C three rows) associated with triangle J for
C J = 1,...,NT.
C
C LOUT = Logical unit number for output. 0 .LE. LOUT
C .LE. 99. Output is printed on unit 6 if LOUT
C is outside its valid range on input.
C
C PRNTX = Logical variable with value TRUE if and only
C if X and Y are to be printed (to 6 decimal
C places).
C
C None of the parameters are altered by this routine.
C
C Modules required by TRLPRT: None
C
C***********************************************************
C
INTEGER I, K, LUN, NA, NB, NL, NLMAX, NMAX
DATA NMAX/9999/, NLMAX/60/
C
C Local parameters:
C
C I = DO-loop, nodal index, and row index for LTRI
C K = DO-loop and triangle index
C LUN = Logical unit number for output
C NA = Number of triangulation arcs
C NB = Number of boundary nodes
C NL = Number of lines printed on the current page
C NLMAX = Maximum number of print lines per page
C NMAX = Maximum value of N and NT (4-digit format)
C
LUN = LOUT
IF (LUN .LT. 0 .OR. LUN .GT. 99) LUN = 6
C
C Print a heading and test for invalid input.
C
WRITE (LUN,100)
NL = 1
IF (N .LT. 3 .OR. N .GT. NMAX .OR.
. (NROW .NE. 6 .AND. NROW .NE. 9) .OR.
. NT .LT. 1 .OR. NT .GT. NMAX) THEN
C
C Print an error message and bypass the loops.
C
WRITE (LUN,110) N, NROW, NT
GO TO 3
ENDIF
IF (PRNTX) THEN
C
C Print X and Y.
C
WRITE (LUN,101)
NL = 6
DO 1 I = 1,N
IF (NL .GE. NLMAX) THEN
WRITE (LUN,106)
NL = 0
ENDIF
WRITE (LUN,102) I, X(I), Y(I)
NL = NL + 1
1 CONTINUE
ENDIF
C
C Print the triangulation LTRI.
C
IF (NL .GT. NLMAX/2) THEN
WRITE (LUN,106)
NL = 0
ENDIF
IF (NROW .EQ. 6) THEN
WRITE (LUN,103)
ELSE
WRITE (LUN,104)
ENDIF
NL = NL + 5
DO 2 K = 1,NT
IF (NL .GE. NLMAX) THEN
WRITE (LUN,106)
NL = 0
ENDIF
WRITE (LUN,105) K, (LTRI(I,K), I = 1,NROW)
NL = NL + 1
2 CONTINUE
C
C Print NB, NA, and NT (boundary nodes, arcs, and
C triangles).
C
NB = 2*N - NT - 2
NA = NT + N - 1
IF (NL .GT. NLMAX-6) WRITE (LUN,106)
WRITE (LUN,107) NB, NA, NT
C
C Print NCC and LCT.
C
3 WRITE (LUN,108) NCC
IF (NCC .GT. 0) WRITE (LUN,109) (LCT(I), I = 1,NCC)
RETURN
C
C Print formats:
C
100 FORMAT (///,24X,'TRIPACK (TRLIST) Output')
101 FORMAT (//16X,'Node',7X,'X(Node)',10X,'Y(Node)'//)
102 FORMAT (16X,I4,2E17.6)
103 FORMAT (//1X,'Triangle',8X,'Vertices',12X,'Neighbors'/
. 4X,'KT',7X,'N1',5X,'N2',5X,'N3',4X,'KT1',4X,
. 'KT2',4X,'KT3'/)
104 FORMAT (//1X,'Triangle',8X,'Vertices',12X,'Neighbors',
. 14X,'Arcs'/
. 4X,'KT',7X,'N1',5X,'N2',5X,'N3',4X,'KT1',4X,
. 'KT2',4X,'KT3',4X,'KA1',4X,'KA2',4X,'KA3'/)
105 FORMAT (2X,I4,2X,6(3X,I4),3(2X,I5))
106 FORMAT (///)
107 FORMAT (/1X,'NB = ',I4,' Boundary Nodes',5X,
. 'NA = ',I5,' Arcs',5X,'NT = ',I5,
. ' Triangles')
108 FORMAT (/1X,'NCC =',I3,' Constraint Curves')
109 FORMAT (1X,9X,14I5)
110 FORMAT (//1X,10X,'*** Invalid Parameter: N =',I5,
. ', NROW =',I5,', NT =',I5,' ***')
END
SUBROUTINE TRMESH (N,X,Y, LIST,LPTR,LEND,LNEW,NEAR,
. NEXT,DIST,IER)
INTEGER N, LIST(*), LPTR(*), LEND(N), LNEW, NEAR(N),
. NEXT(N), IER
REAL X(N), Y(N), DIST(N)
C
C***********************************************************
C
C From TRIPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C renka@cs.unt.edu
C 06/28/98
C
C This subroutine creates a Delaunay triangulation of a
C set of N arbitrarily distributed points in the plane re-
C ferred to as nodes. The Delaunay triangulation is defined
C as a set of triangles with the following five properties:
C
C 1) The triangle vertices are nodes.
C 2) No triangle contains a node other than its vertices.
C 3) The interiors of the triangles are pairwise disjoint.
C 4) The union of triangles is the convex hull of the set
C of nodes (the smallest convex set which contains
C the nodes).
C 5) The interior of the circumcircle of each triangle
C contains no node.
C
C The first four properties define a triangulation, and the
C last property results in a triangulation which is as close
C as possible to equiangular in a certain sense and which is
C uniquely defined unless four or more nodes lie on a common
C circle. This property makes the triangulation well-suited
C for solving closest point problems and for triangle-based
C interpolation.
C
C The triangulation can be generalized to a constrained
C Delaunay triangulation by a call to Subroutine ADDCST.
C This allows for user-specified boundaries defining a non-
C convex and/or multiply connected region.
C
C The algorithm for constructing the triangulation has
C expected time complexity O(N*log(N)) for most nodal dis-
C tributions. Also, since the algorithm proceeds by adding
C nodes incrementally, the triangulation may be updated with
C the addition (or deletion) of a node very efficiently.
C The adjacency information representing the triangulation
C is stored as a linked list requiring approximately 13N
C storage locations.
C
C
C The following is a list of the software package modules
C which a user may wish to call directly:
C
C ADDCST - Generalizes the Delaunay triangulation to allow
C for user-specified constraints.
C
C ADDNOD - Updates the triangulation by appending or
C inserting a new node.
C
C AREAP - Computes the area bounded by a closed polygonal
C curve such as the boundary of the triangula-
C tion or of a constraint region.
C
C BNODES - Returns an array containing the indexes of the
C boundary nodes in counterclockwise order.
C Counts of boundary nodes, triangles, and arcs
C are also returned.
C
C CIRCUM - Computes the area, circumcenter, circumradius,
C and, optionally, the aspect ratio of a trian-
C gle defined by user-specified vertices.
C
C DELARC - Deletes a boundary arc from the triangulation.
C
C DELNOD - Updates the triangulation with the deletion of a
C node.
C
C EDGE - Forces a pair of nodes to be connected by an arc
C in the triangulation.
C
C GETNP - Determines the ordered sequence of L closest
C nodes to a given node, along with the associ-
C ated distances. The distance between nodes is
C taken to be the length of the shortest connec-
C ting path which intersects no constraint
C region.
C
C INTSEC - Determines whether or not an arbitrary pair of
C line segments share a common point.
C
C JRAND - Generates a uniformly distributed pseudo-random
C integer.
C
C LEFT - Locates a point relative to a line.
C
C NEARND - Returns the index of the nearest node to an
C arbitrary point, along with its squared
C distance.
C
C STORE - Forces a value to be stored in main memory so
C that the precision of floating point numbers
C in memory locations rather than registers is
C computed.
C
C TRLIST - Converts the triangulation data structure to a
C triangle list more suitable for use in a fin-
C ite element code.
C
C TRLPRT - Prints the triangle list created by Subroutine
C TRLIST.
C
C TRMESH - Creates a Delaunay triangulation of a set of
C nodes.
C
C TRMSHR - Creates a Delaunay triangulation (more effici-
C ently than TRMESH) of a set of nodes lying at
C the vertices of a (possibly skewed) rectangu-
C lar grid.
C
C TRPLOT - Creates a level-2 Encapsulated Postscript (EPS)
C file containing a triangulation plot.
C
C TRPRNT - Prints the triangulation data structure and,
C optionally, the nodal coordinates.
C
C
C On input:
C
C N = Number of nodes in the triangulation. N .GE. 3.
C
C X,Y = Arrays of length N containing the Cartesian
C coordinates of the nodes. (X(K),Y(K)) is re-
C ferred to as node K, and K is referred to as
C a nodal index. The first three nodes must not
C be collinear.
C
C The above parameters are not altered by this routine.
C
C LIST,LPTR = Arrays of length at least 6N-12.
C
C LEND = Array of length at least N.
C
C NEAR,NEXT,DIST = Work space arrays of length at
C least N. The space is used to
C efficiently determine the nearest
C triangulation node to each un-
C processed node for use by ADDNOD.
C
C On output:
C
C LIST = Set of nodal indexes which, along with LPTR,
C LEND, and LNEW, define the triangulation as a
C set of N adjacency lists -- counterclockwise-
C ordered sequences of neighboring nodes such
C that the first and last neighbors of a bound-
C ary node are boundary nodes (the first neigh-
C bor of an interior node is arbitrary). In
C order to distinguish between interior and
C boundary nodes, the last neighbor of each
C boundary node is represented by the negative
C of its index.
C
C LPTR = Set of pointers (LIST indexes) in one-to-one
C correspondence with the elements of LIST.
C LIST(LPTR(I)) indexes the node which follows
C LIST(I) in cyclical counterclockwise order
C (the first neighbor follows the last neigh-
C bor).
C
C LEND = Set of pointers to adjacency lists. LEND(K)
C points to the last neighbor of node K for
C K = 1,...,N. Thus, LIST(LEND(K)) < 0 if and
C only if K is a boundary node.
C
C LNEW = Pointer to the first empty location in LIST
C and LPTR (list length plus one). LIST, LPTR,
C LEND, and LNEW are not altered if IER < 0,
C and are incomplete if IER > 0.
C
C NEAR,NEXT,DIST = Garbage.
C
C IER = Error indicator:
C IER = 0 if no errors were encountered.
C IER = -1 if N < 3 on input.
C IER = -2 if the first three nodes are
C collinear.
C IER = -4 if an error flag was returned by a
C call to SWAP in ADDNOD. This is an
C internal error and should be reported
C to the programmer.
C IER = L if nodes L and M coincide for some
C M > L. The linked list represents
C a triangulation of nodes 1 to M-1
C in this case.
C
C Modules required by TRMESH: ADDNOD, BDYADD, INSERT,
C INTADD, JRAND, LEFT,
C LSTPTR, STORE, SWAP,
C SWPTST, TRFIND
C
C Intrinsic function called by TRMESH: ABS
C
C***********************************************************
C
LOGICAL LEFT
REAL STORE
INTEGER I, I0, J, K, KM1, LCC(1), LP, LPL, NCC, NEXTI,
. NN
REAL D, D1, D2, D3, EPS, SWTOL
COMMON/SWPCOM/SWTOL
C
C Local parameters:
C
C D = Squared distance from node K to node I
C D1,D2,D3 = Squared distances from node K to nodes 1, 2,
C and 3, respectively
C EPS = Half the machine precision
C I,J = Nodal indexes
C I0 = Index of the node preceding I in a sequence of
C unprocessed nodes: I = NEXT(I0)
C K = Index of node to be added and DO-loop index:
C K > 3
C KM1 = K-1
C LCC(1) = Dummy array
C LP = LIST index (pointer) of a neighbor of K
C LPL = Pointer to the last neighbor of K
C NCC = Number of constraint curves
C NEXTI = NEXT(I)
C NN = Local copy of N
C SWTOL = Tolerance for function SWPTST
C
NN = N
IF (NN .LT. 3) THEN
IER = -1
RETURN
ENDIF
C
C Compute a tolerance for function SWPTST: SWTOL = 10*
C (machine precision)
C
EPS = 1.
1 EPS = EPS/2.
SWTOL = STORE(EPS + 1.)
IF (SWTOL .GT. 1.) GO TO 1
SWTOL = EPS*20.
C
C Store the first triangle in the linked list.
C
IF ( .NOT. LEFT(X(1),Y(1),X(2),Y(2),X(3),Y(3)) ) THEN
C
C The initial triangle is (3,2,1) = (2,1,3) = (1,3,2).
C
LIST(1) = 3
LPTR(1) = 2
LIST(2) = -2
LPTR(2) = 1
LEND(1) = 2
C
LIST(3) = 1
LPTR(3) = 4
LIST(4) = -3
LPTR(4) = 3
LEND(2) = 4
C
LIST(5) = 2
LPTR(5) = 6
LIST(6) = -1
LPTR(6) = 5
LEND(3) = 6
C
ELSEIF ( .NOT. LEFT(X(2),Y(2),X(1),Y(1),X(3),Y(3)) )
. THEN
C
C The initial triangle is (1,2,3).
C
LIST(1) = 2
LPTR(1) = 2
LIST(2) = -3
LPTR(2) = 1
LEND(1) = 2
C
LIST(3) = 3
LPTR(3) = 4
LIST(4) = -1
LPTR(4) = 3
LEND(2) = 4
C
LIST(5) = 1
LPTR(5) = 6
LIST(6) = -2
LPTR(6) = 5
LEND(3) = 6
C
ELSE
C
C The first three nodes are collinear.
C
IER = -2
RETURN
ENDIF
C
C Initialize LNEW and test for N = 3.
C
LNEW = 7
IF (NN .EQ. 3) THEN
IER = 0
RETURN
ENDIF
C
C A nearest-node data structure (NEAR, NEXT, and DIST) is
C used to obtain an expected-time (N*log(N)) incremental
C algorithm by enabling constant search time for locating
C each new node in the triangulation.
C
C For each unprocessed node K, NEAR(K) is the index of the
C triangulation node closest to K (used as the starting
C point for the search in Subroutine TRFIND) and DIST(K)
C is an increasing function of the distance between nodes
C K and NEAR(K).
C
C Since it is necessary to efficiently find the subset of
C unprocessed nodes associated with each triangulation
C node J (those that have J as their NEAR entries), the
C subsets are stored in NEAR and NEXT as follows: for
C each node J in the triangulation, I = NEAR(J) is the
C first unprocessed node in J's set (with I = 0 if the
C set is empty), L = NEXT(I) (if I > 0) is the second,
C NEXT(L) (if L > 0) is the third, etc. The nodes in each
C set are initially ordered by increasing indexes (which
C maximizes efficiency) but that ordering is not main-
C tained as the data structure is updated.
C
C Initialize the data structure for the single triangle.
C
NEAR(1) = 0
NEAR(2) = 0
NEAR(3) = 0
DO 2 K = NN,4,-1
D1 = (X(K)-X(1))**2 + (Y(K)-Y(1))**2
D2 = (X(K)-X(2))**2 + (Y(K)-Y(2))**2
D3 = (X(K)-X(3))**2 + (Y(K)-Y(3))**2
IF (D1 .LE. D2 .AND. D1 .LE. D3) THEN
NEAR(K) = 1
DIST(K) = D1
NEXT(K) = NEAR(1)
NEAR(1) = K
ELSEIF (D2 .LE. D1 .AND. D2 .LE. D3) THEN
NEAR(K) = 2
DIST(K) = D2
NEXT(K) = NEAR(2)
NEAR(2) = K
ELSE
NEAR(K) = 3
DIST(K) = D3
NEXT(K) = NEAR(3)
NEAR(3) = K
ENDIF
2 CONTINUE
C
C Add the remaining nodes. Parameters for ADDNOD are as
C follows:
C
C K = Index of the node to be added.
C NEAR(K) = Index of the starting node for the search in
C TRFIND.
C NCC = Number of constraint curves.
C LCC = Dummy array (since NCC = 0).
C KM1 = Number of nodes in the triangulation.
C
NCC = 0
DO 7 K = 4,NN
KM1 = K-1
CALL ADDNOD (K,X(K),Y(K),NEAR(K),NCC, LCC,KM1,X,Y,
. LIST,LPTR,LEND,LNEW, IER)
IF (IER .NE. 0) RETURN
C
C Remove K from the set of unprocessed nodes associated
C with NEAR(K).
C
I = NEAR(K)
IF (NEAR(I) .EQ. K) THEN
NEAR(I) = NEXT(K)
ELSE
I = NEAR(I)
3 I0 = I
I = NEXT(I0)
IF (I .NE. K) GO TO 3
NEXT(I0) = NEXT(K)
ENDIF
NEAR(K) = 0
C
C Loop on neighbors J of node K.
C
LPL = LEND(K)
LP = LPL
4 LP = LPTR(LP)
J = ABS(LIST(LP))
C
C Loop on elements I in the sequence of unprocessed nodes
C associated with J: K is a candidate for replacing J
C as the nearest triangulation node to I. The next value
C of I in the sequence, NEXT(I), must be saved before I
C is moved because it is altered by adding I to K's set.
C
I = NEAR(J)
5 IF (I .EQ. 0) GO TO 6
NEXTI = NEXT(I)
C
C Test for the distance from I to K less than the distance
C from I to J.
C
D = (X(K)-X(I))**2 + (Y(K)-Y(I))**2
IF (D .LT. DIST(I)) THEN
C
C Replace J by K as the nearest triangulation node to I:
C update NEAR(I) and DIST(I), and remove I from J's set
C of unprocessed nodes and add it to K's set.
C
NEAR(I) = K
DIST(I) = D
IF (I .EQ. NEAR(J)) THEN
NEAR(J) = NEXTI
ELSE
NEXT(I0) = NEXTI
ENDIF
NEXT(I) = NEAR(K)
NEAR(K) = I
ELSE
I0 = I
ENDIF
C
C Bottom of loop on I.
C
I = NEXTI
GO TO 5
C
C Bottom of loop on neighbors J.
C
6 IF (LP .NE. LPL) GO TO 4
7 CONTINUE
RETURN
END
SUBROUTINE TRMSHR (N,NX,X,Y, NIT, LIST,LPTR,LEND,LNEW,
. IER)
INTEGER N, NX, NIT, LIST(*), LPTR(*), LEND(N), LNEW,
. IER
REAL X(N), Y(N)
C
C***********************************************************
C
C From TRIPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C renka@cs.unt.edu
C 06/27/98
C
C This subroutine creates a Delaunay triangulation of a
C set of N nodes in the plane, where the nodes are the vert-
C ices of an NX by NY skewed rectangular grid with the
C natural ordering. Thus, N = NX*NY, and the nodes are
C ordered from left to right beginning at the top row so
C that adjacent nodes have indexes which differ by 1 in the
C x-direction and by NX in the y-direction. A skewed rec-
C tangular grid is defined as one in which each grid cell is
C a strictly convex quadrilateral (and is thus the convex
C hull of its four vertices). Equivalently, any transfor-
C mation from a rectangle to a grid cell which is bilinear
C in both components has an invertible Jacobian.
C
C If the nodes are not distributed and ordered as defined
C above, Subroutine TRMESH must be called in place of this
C routine. Refer to Subroutine ADDCST for the treatment of
C constraints.
C
C The first phase of the algorithm consists of construc-
C ting a triangulation by choosing a diagonal arc in each
C grid cell. If NIT = 0, all diagonals connect lower left
C to upper right corners and no error checking or additional
C computation is performed. Otherwise, each diagonal arc is
C chosen to be locally optimal, and boundary arcs are added
C where necessary in order to cover the convex hull of the
C nodes. (This is the first iteration.) If NIT > 1 and no
C error was detected, the triangulation is then optimized by
C a sequence of up to NIT-1 iterations in which interior
C arcs of the triangulation are tested and swapped if appro-
C priate. The algorithm terminates when an iteration
C results in no swaps and/or when the allowable number of
C iterations has been performed. NIT = 0 is sufficient to
C produce a Delaunay triangulation if the original grid is
C actually rectangular, and NIT = 1 is sufficient if it is
C close to rectangular. Note, however, that the ordering
C and distribution of nodes is not checked for validity in
C the case NIT = 0, and the triangulation will not be valid
C unless the rectangular grid covers the convex hull of the
C nodes.
C
C
C On input:
C
C N = Number of nodes in the grid. N = NX*NY for some
C NY .GE. 2.
C
C NX = Number of grid points in the x-direction. NX
C .GE. 2.
C
C X,Y = Arrays of length N containing coordinates of
C the nodes with the ordering and distribution
C defined in the header comments above.
C (X(K),Y(K)) is referred to as node K.
C
C The above parameters are not altered by this routine.
C
C NIT = Nonnegative integer specifying the maximum
C number of iterations to be employed. Refer
C to the header comments above.
C
C LIST,LPTR = Arrays of length at least 6N-12.
C
C LEND = Array of length at least N.
C
C On output:
C
C NIT = Number of iterations employed.
C
C LIST,LPTR,LEND,LNEW = Data structure defining the
C triangulation. Refer to Sub-
C routine TRMESH.
C
C IER = Error indicator:
C IER = 0 if no errors were encountered.
C IER = K if the grid element with upper left
C corner at node K is not a strictly
C convex quadrilateral. The algorithm
C is terminated when the first such
C occurrence is detected. Note that
C this test is not performed if NIT = 0
C on input.
C IER = -1 if N, NX, or NIT is outside its valid
C range on input.
C IER = -2 if NIT > 1 on input, and the optimi-
C zation loop failed to converge within
C the allowable number of iterations.
C The triangulation is valid but not
C optimal in this case.
C
C Modules required by TRMSHR: INSERT, LEFT, LSTPTR, NBCNT,
C STORE, SWAP, SWPTST
C
C Intrinsic function called by TRMSHR: ABS
C
C***********************************************************
C
INTEGER LSTPTR, NBCNT
LOGICAL LEFT, SWPTST
REAL STORE
INTEGER I, ITER, J, K, KP1, LP, LPF, LPK, LPL, LPP,
. M1, M2, M3, M4, MAXIT, N0, N1, N2, N3, N4, NI,
. NJ, NM1, NN, NNB
LOGICAL TST
REAL EPS, SWTOL
COMMON/SWPCOM/SWTOL
C
C Store local variables and test for errors in input
C parameters.
C
NI = NX
NJ = N/NI
NN = NI*NJ
MAXIT = NIT
NIT = 0
IF (N .NE. NN .OR. NJ .LT. 2 .OR. NI .LT. 2 .OR.
. MAXIT .LT. 0) THEN
IER = -1
RETURN
ENDIF
IER = 0
C
C Compute a tolerance for function SWPTST: SWTOL = 10*
C (machine precision)
C
EPS = 1.
1 EPS = EPS/2.
SWTOL = STORE(EPS + 1.)
IF (SWTOL .GT. 1.) GO TO 1
SWTOL = EPS*20.
C
C Loop on grid points (I,J) corresponding to nodes K =
C (J-1)*NI + I. TST = TRUE iff diagonals are to be
C chosen by the swap test. M1, M2, M3, and M4 are the
C slopes (-1, 0, or 1) of the diagonals in quadrants 1
C to 4 (counterclockwise beginning with the upper right)
C for a coordinate system with origin at node K.
C
TST = MAXIT .GT. 0
M1 = 0
M4 = 0
LP = 0
KP1 = 1
DO 6 J = 1,NJ
DO 5 I = 1,NI
M2 = M1
M3 = M4
K = KP1
KP1 = K + 1
LPF = LP + 1
IF (J .EQ. NJ .AND. I .NE. NI) GO TO 2
IF (I .NE. 1) THEN
IF (J .NE. 1) THEN
C
C K is not in the top row, leftmost column, or bottom row
C (unless K is the lower right corner). Take the first
C neighbor to be the node above K.
C
LP = LP + 1
LIST(LP) = K - NI
LPTR(LP) = LP + 1
IF (M2 .LE. 0) THEN
LP = LP + 1
LIST(LP) = K - 1 - NI
LPTR(LP) = LP + 1
ENDIF
ENDIF
C
C K is not in the leftmost column. The next (or first)
C neighbor is to the left of K.
C
LP = LP + 1
LIST(LP) = K - 1
LPTR(LP) = LP + 1
IF (J .EQ. NJ) GO TO 3
IF (M3 .GE. 0) THEN
LP = LP + 1
LIST(LP) = K - 1 + NI
LPTR(LP) = LP + 1
ENDIF
ENDIF
C
C K is not in the bottom row. The next (or first)
C neighbor is below K.
C
LP = LP + 1
LIST(LP) = K + NI
LPTR(LP) = LP + 1
C
C Test for a negative diagonal in quadrant 4 unless K is
C in the rightmost column. The quadrilateral associated
C with the quadrant is tested for strict convexity un-
C less NIT = 0 on input.
C
IF (I .EQ. NI) GO TO 3
M4 = 1
IF (.NOT. TST) GO TO 2
IF ( LEFT(X(KP1),Y(KP1),X(K+NI),Y(K+NI),X(K),Y(K))
. .OR. LEFT(X(K),Y(K),X(KP1+NI),Y(KP1+NI),
. X(K+NI),Y(K+NI))
. .OR. LEFT(X(K+NI),Y(K+NI),X(KP1),Y(KP1),
. X(KP1+NI),Y(KP1+NI))
. .OR. LEFT(X(KP1+NI),Y(KP1+NI),X(K),Y(K),
. X(KP1),Y(KP1)) ) GO TO 12
IF ( SWPTST(KP1,K+NI,K,KP1+NI,X,Y) ) GO TO 2
M4 = -1
LP = LP + 1
LIST(LP) = KP1 + NI
LPTR(LP) = LP + 1
C
C The next (or first) neighbor is to the right of K.
C
2 LP = LP + 1
LIST(LP) = KP1
LPTR(LP) = LP + 1
C
C Test for a positive diagonal in quadrant 1 (the neighbor
C of K-NI which follows K is not K+1) unless K is in the
C top row.
C
IF (J .EQ. 1) GO TO 3
IF (TST) THEN
M1 = -1
LPK = LSTPTR(LEND(K-NI),K,LIST,LPTR)
LPK = LPTR(LPK)
IF (LIST(LPK) .NE. KP1) THEN
M1 = 1
LP = LP + 1
LIST(LP) = KP1 - NI
LPTR(LP) = LP + 1
ENDIF
ENDIF
C
C If K is in the leftmost column (and not the top row) or
C in the bottom row (and not the rightmost column), then
C the next neighbor is the node above K.
C
IF (I .NE. 1 .AND. J .NE. NJ) GO TO 4
LP = LP + 1
LIST(LP) = K - NI
LPTR(LP) = LP + 1
IF (I .EQ. 1) GO TO 3
C
C K is on the bottom row (and not the leftmost or right-
C most column).
C
IF (M2 .LE. 0) THEN
LP = LP + 1
LIST(LP) = K - 1 - NI
LPTR(LP) = LP + 1
ENDIF
LP = LP + 1
LIST(LP) = K - 1
LPTR(LP) = LP + 1
C
C K is a boundary node.
C
3 LIST(LP) = -LIST(LP)
C
C Bottom of loop. Store LEND and correct LPTR(LP).
C LPF and LP point to the first and last neighbors
C of K.
C
4 LEND(K) = LP
LPTR(LP) = LPF
5 CONTINUE
6 CONTINUE
C
C Store LNEW, and terminate the algorithm if NIT = 0 on
C input.
C
LNEW = LP + 1
IF (MAXIT .EQ. 0) RETURN
C
C Add boundary arcs where necessary in order to cover the
C convex hull of the nodes. N1, N2, and N3 are consecu-
C tive boundary nodes in counterclockwise order, and N0
C is the starting point for each loop around the boundary.
C
N0 = 1
N1 = N0
N2 = NI + 1
C
C TST is set to TRUE if an arc is added. The boundary
C loop is repeated until a traversal results in no
C added arcs.
C
7 TST = .FALSE.
C
C Top of boundary loop. Set N3 to the first neighbor of
C N2, and test for N3 LEFT N1 -> N2.
C
8 LPL = LEND(N2)
LP = LPTR(LPL)
N3 = LIST(LP)
IF ( LEFT(X(N1),Y(N1),X(N2),Y(N2),X(N3),Y(N3)) )
. N1 = N2
IF (N1 .NE. N2) THEN
C
C Add the boundary arc N1-N3. If N0 = N2, the starting
C point is changed to N3, since N2 will be removed from
C the boundary. N3 is inserted as the first neighbor of
C N1, N2 is changed to an interior node, and N1 is
C inserted as the last neighbor of N3.
C
TST = .TRUE.
IF (N2 .EQ. N0) N0 = N3
LP = LEND(N1)
CALL INSERT (N3,LP, LIST,LPTR,LNEW )
LIST(LPL) = -LIST(LPL)
LP = LEND(N3)
LIST(LP) = N2
CALL INSERT (-N1,LP, LIST,LPTR,LNEW )
LEND(N3) = LNEW - 1
ENDIF
C
C Bottom of loops. Test for termination.
C
N2 = N3
IF (N1 .NE. N0) GO TO 8
IF (TST) GO TO 7
C
C Terminate the algorithm if NIT = 1 on input.
C
NIT = 1
IF (MAXIT .EQ. 1) RETURN
C
C Optimize the triangulation by applying the swap test and
C appropriate swaps to the interior arcs. The loop is
C repeated until no swaps are performed or MAXIT itera-
C tions have been applied. ITER is the current iteration,
C and TST is set to TRUE if a swap occurs.
C
ITER = 1
NM1 = NN - 1
9 ITER = ITER + 1
TST = .FALSE.
C
C Loop on interior arcs N1-N2, where N2 > N1 and
C (N1,N2,N3) and (N2,N1,N4) are adjacent triangles.
C
C Top of loop on nodes N1.
C
DO 11 N1 = 1,NM1
LPL = LEND(N1)
N4 = LIST(LPL)
LPF = LPTR(LPL)
N2 = LIST(LPF)
LP = LPTR(LPF)
N3 = LIST(LP)
NNB = NBCNT(LPL,LPTR)
C
C Top of loop on neighbors N2 of N1. NNB is the number of
C neighbors of N1.
C
DO 10 I = 1,NNB
C
C Bypass the swap test if N1 is a boundary node and N2 is
C the first neighbor (N4 < 0), N2 < N1, or N1-N2 is a
C diagonal arc (already locally optimal) when ITER = 2.
C
IF ( N4 .GT. 0 .AND. N2 .GT. N1 .AND.
. (ITER .NE. 2 .OR. ABS(N1+NI-N2) .NE. 1) )
. THEN
IF (SWPTST(N3,N4,N1,N2,X,Y) ) THEN
C
C Swap diagonal N1-N2 for N3-N4, set TST to TRUE, and set
C N2 to N4 (the neighbor preceding N3).
C
CALL SWAP (N3,N4,N1,N2, LIST,LPTR,LEND, LPP)
IF (LPP .NE. 0) THEN
TST = .TRUE.
N2 = N4
ENDIF
ENDIF
ENDIF
C
C Bottom of neighbor loop.
C
IF (LIST(LPL) .EQ. -N3) GO TO 11
N4 = N2
N2 = N3
LP = LSTPTR(LPL,N2,LIST,LPTR)
LP = LPTR(LP)
N3 = ABS(LIST(LP))
10 CONTINUE
11 CONTINUE
C
C Test for termination.
C
IF (TST .AND. ITER .LT. MAXIT) GO TO 9
NIT = ITER
IF (TST) IER = -2
RETURN
C
C Invalid grid cell encountered.
C
12 IER = K
RETURN
END
SUBROUTINE TRPLOT (LUN,PLTSIZ,WX1,WX2,WY1,WY2,NCC,LCC,
. N,X,Y,LIST,LPTR,LEND,TITLE,
. NUMBR, IER)
CHARACTER*(*) TITLE
INTEGER LUN, NCC, LCC(*), N, LIST(*), LPTR(*),
. LEND(N), IER
LOGICAL NUMBR
REAL PLTSIZ, WX1, WX2, WY1, WY2, X(N), Y(N)
C
C***********************************************************
C
C From TRIPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C renka@cs.unt.edu
C 07/15/98
C
C This subroutine creates a level-2 Encapsulated Post-
C script (EPS) file containing a triangulation plot.
C
C
C On input:
C
C LUN = Logical unit number in the range 0 to 99.
C The unit should be opened with an appropriate
C file name before the call to this routine.
C
C PLTSIZ = Plot size in inches. The window is mapped,
C with aspect ratio preserved, to a rectangu-
C lar viewport with maximum side-length equal
C to .88*PLTSIZ (leaving room for labels out-
C side the viewport). The viewport is
C centered on the 8.5 by 11 inch page, and
C its boundary is drawn. 1.0 .LE. PLTSIZ
C .LE. 8.5.
C
C WX1,WX2,WY1,WY2 = Parameters defining a rectangular
C window against which the triangu-
C lation is clipped. (Only the
C portion of the triangulation that
C lies in the window is drawn.)
C (WX1,WY1) and (WX2,WY2) are the
C lower left and upper right cor-
C ners, respectively. WX1 < WX2 and
C WY1 < WY2.
C
C NCC = Number of constraint curves. Refer to Subrou-
C tine ADDCST. NCC .GE. 0.
C
C LCC = Array of length NCC (or dummy parameter if
C NCC = 0) containing the index of the first
C node of constraint I in LCC(I). For I = 1 to
C NCC, LCC(I+1)-LCC(I) .GE. 3, where LCC(NCC+1)
C = N+1.
C
C N = Number of nodes in the triangulation. N .GE. 3.
C
C X,Y = Arrays of length N containing the coordinates
C of the nodes with non-constraint nodes in the
C first LCC(1)-1 locations.
C
C LIST,LPTR,LEND = Data structure defining the trian-
C gulation. Refer to Subroutine
C TRMESH.
C
C TITLE = Type CHARACTER variable or constant contain-
C ing a string to be centered above the plot.
C The string must be enclosed in parentheses;
C i.e., the first and last characters must be
C '(' and ')', respectively, but these are not
C displayed. TITLE may have at most 80 char-
C acters including the parentheses.
C
C NUMBR = Option indicator: If NUMBR = TRUE, the
C nodal indexes are plotted next to the nodes.
C
C Input parameters are not altered by this routine.
C
C On output:
C
C IER = Error indicator:
C IER = 0 if no errors were encountered.
C IER = 1 if LUN, PLTSIZ, NCC, or N is outside
C its valid range. LCC is not tested
C for validity.
C IER = 2 if WX1 >= WX2 or WY1 >= WY2.
C IER = 3 if an error was encountered in writing
C to unit LUN.
C
C Various plotting options can be controlled by altering
C the data statement below.
C
C Modules required by TRPLOT: None
C
C Intrinsic functions called by TRPLOT: ABS, CHAR, NINT,
C REAL
C
C***********************************************************
C
INTEGER I, IFRST, IH, ILAST, IPX1, IPX2, IPY1, IPY2,
. IW, LP, LPL, N0, N0BAK, N0FOR, N1, NLS
LOGICAL ANNOT, CNSTR, PASS1
REAL DASHL, DX, DY, FSIZN, FSIZT, R, SFX, SFY, T,
. TX, TY, X0, Y0
C
DATA ANNOT/.TRUE./, DASHL/4.0/, FSIZN/10.0/,
. FSIZT/16.0/
C
C Local parameters:
C
C ANNOT = Logical variable with value TRUE iff the plot
C is to be annotated with the values of WX1,
C WX2, WY1, and WY2
C CNSTR Logical variable used to flag constraint arcs:
C TRUE iff N0-N1 lies in a constraint region
C DASHL = Length (in points, at 72 points per inch) of
C dashes and spaces in a dashed line pattern
C used for drawing constraint arcs
C DX = Window width WX2-WX1
C DY = Window height WY2-WY1
C FSIZN = Font size in points for labeling nodes with
C their indexes if NUMBR = TRUE
C FSIZT = Font size in points for the title (and
C annotation if ANNOT = TRUE)
C I = Constraint index (1 to NCC)
C IFRST = Index of the first node in constraint I
C IH = Height of the viewport in points
C ILAST = Index of the last node in constraint I
C IPX1,IPY1 = X and y coordinates (in points) of the lower
C left corner of the bounding box or viewport
C IPX2,IPY2 = X and y coordinates (in points) of the upper
C right corner of the bounding box or viewport
C IW = Width of the viewport in points
C LP = LIST index (pointer)
C LPL = Pointer to the last neighbor of N0
C N0 = Nodal index and DO-loop index
C N0BAK = Predecessor of N0 in a constraint curve
C (sequence of adjacent constraint nodes)
C N0FOR = Successor to N0 in a constraint curve
C N1 = Index of a neighbor of N0
C NLS = Index of the last non-constraint node
C PASS1 = Logical variable used to flag the first pass
C through the constraint nodes
C R = Aspect ratio DX/DY
C SFX,SFY = Scale factors for mapping world coordinates
C (window coordinates in [WX1,WX2] X [WY1,WY2])
C to viewport coordinates in [IPX1,IPX2] X
C [IPY1,IPY2]
C T = Temporary variable
C TX,TY = Translation vector for mapping world coordi-
C nates to viewport coordinates
C X0,Y0 = X(N0),Y(N0) or label location
C
C
C Test for error 1, and set NLS to the last non-constraint
C node.
C
IF (LUN .LT. 0 .OR. LUN .GT. 99 .OR.
. PLTSIZ .LT. 1.0 .OR. PLTSIZ .GT. 8.5 .OR.
. NCC .LT. 0 .OR. N .LT. 3) GO TO 11
NLS = N
IF (NCC .GT. 0) NLS = LCC(1)-1
C
C Compute the aspect ratio of the window.
C
DX = WX2 - WX1
DY = WY2 - WY1
IF (DX .LE. 0.0 .OR. DY .LE. 0.0) GO TO 12
R = DX/DY
C
C Compute the lower left (IPX1,IPY1) and upper right
C (IPX2,IPY2) corner coordinates of the bounding box.
C The coordinates, specified in default user space units
C (points, at 72 points/inch with origin at the lower
C left corner of the page), are chosen to preserve the
C aspect ratio R, and to center the plot on the 8.5 by 11
C inch page. The center of the page is (306,396), and
C T = PLTSIZ/2 in points.
C
T = 36.0*PLTSIZ
IF (R .GE. 1.0) THEN
IPX1 = 306 - NINT(T)
IPX2 = 306 + NINT(T)
IPY1 = 396 - NINT(T/R)
IPY2 = 396 + NINT(T/R)
ELSE
IPX1 = 306 - NINT(T*R)
IPX2 = 306 + NINT(T*R)
IPY1 = 396 - NINT(T)
IPY2 = 396 + NINT(T)
ENDIF
C
C Output header comments.
C
WRITE (LUN,100,ERR=13) IPX1, IPY1, IPX2, IPY2
100 FORMAT ('%!PS-Adobe-3.0 EPSF-3.0'/
. '%%BoundingBox:',4I4/
. '%%Title: Triangulation'/
. '%%Creator: TRIPACK'/
. '%%EndComments')
C
C Set (IPX1,IPY1) and (IPX2,IPY2) to the corner coordinates
C of a viewport obtained by shrinking the bounding box by
C 12% in each dimension.
C
IW = NINT(0.88*REAL(IPX2-IPX1))
IH = NINT(0.88*REAL(IPY2-IPY1))
IPX1 = 306 - IW/2
IPX2 = 306 + IW/2
IPY1 = 396 - IH/2
IPY2 = 396 + IH/2
C
C Set the line thickness to 2 points, and draw the
C viewport boundary.
C
T = 2.0
WRITE (LUN,110,ERR=13) T
WRITE (LUN,120,ERR=13) IPX1, IPY1
WRITE (LUN,130,ERR=13) IPX1, IPY2
WRITE (LUN,130,ERR=13) IPX2, IPY2
WRITE (LUN,130,ERR=13) IPX2, IPY1
WRITE (LUN,140,ERR=13)
WRITE (LUN,150,ERR=13)
110 FORMAT (F12.6,' setlinewidth')
120 FORMAT (2I4,' moveto')
130 FORMAT (2I4,' lineto')
140 FORMAT ('closepath')
150 FORMAT ('stroke')
C
C Set up a mapping from the window to the viewport.
C
SFX = REAL(IW)/DX
SFY = REAL(IH)/DY
TX = IPX1 - SFX*WX1
TY = IPY1 - SFY*WY1
WRITE (LUN,160,ERR=13) TX, TY, SFX, SFY
160 FORMAT (2F12.6,' translate'/
. 2F12.6,' scale')
C
C The line thickness (believe it or fucking not) must be
C changed to reflect the new scaling which is applied to
C all subsequent output. Set it to 1.0 point.
C
T = 2.0/(SFX+SFY)
WRITE (LUN,110,ERR=13) T
C
C Save the current graphics state, and set the clip path to
C the boundary of the window.
C
WRITE (LUN,170,ERR=13)
WRITE (LUN,180,ERR=13) WX1, WY1
WRITE (LUN,190,ERR=13) WX2, WY1
WRITE (LUN,190,ERR=13) WX2, WY2
WRITE (LUN,190,ERR=13) WX1, WY2
WRITE (LUN,200,ERR=13)
170 FORMAT ('gsave')
180 FORMAT (2F12.6,' moveto')
190 FORMAT (2F12.6,' lineto')
200 FORMAT ('closepath clip newpath')
C
C Draw the edges N0->N1, where N1 > N0, beginning with a
C loop on non-constraint nodes N0. LPL points to the
C last neighbor of N0.
C
DO 3 N0 = 1,NLS
X0 = X(N0)
Y0 = Y(N0)
LPL = LEND(N0)
LP = LPL
C
C Loop on neighbors N1 of N0.
C
2 LP = LPTR(LP)
N1 = ABS(LIST(LP))
IF (N1 .GT. N0) THEN
C
C Add the edge to the path.
C
WRITE (LUN,210,ERR=13) X0, Y0, X(N1), Y(N1)
210 FORMAT (2F12.6,' moveto',2F12.6,' lineto')
ENDIF
IF (LP .NE. LPL) GO TO 2
3 CONTINUE
C
C Loop through the constraint nodes twice. The non-
C constraint arcs incident on constraint nodes are
C drawn (with solid lines) on the first pass, and the
C constraint arcs (both boundary and interior, if any)
C are drawn (with dashed lines) on the second pass.
C
PASS1 = .TRUE.
C
C Loop on constraint nodes N0 with (N0BAK,N0,N0FOR) a sub-
C sequence of constraint I. The outer loop is on
C constraints I with first and last nodes IFRST and ILAST.
C
4 IFRST = N+1
DO 8 I = NCC,1,-1
ILAST = IFRST - 1
IFRST = LCC(I)
N0BAK = ILAST
DO 7 N0 = IFRST,ILAST
N0FOR = N0 + 1
IF (N0 .EQ. ILAST) N0FOR = IFRST
LPL = LEND(N0)
X0 = X(N0)
Y0 = Y(N0)
LP = LPL
C
C Loop on neighbors N1 of N0. CNSTR = TRUE iff N0-N1 is a
C constraint arc.
C
C Initialize CNSTR to TRUE iff the first neighbor of N0
C strictly follows N0FOR and precedes or coincides with
C N0BAK (in counterclockwise order).
C
5 LP = LPTR(LP)
N1 = ABS(LIST(LP))
IF (N1 .NE. N0FOR .AND. N1 .NE. N0BAK) GO TO 5
CNSTR = N1 .EQ. N0BAK
LP = LPL
C
C Loop on neighbors N1 of N0. Update CNSTR and test for
C N1 > N0.
C
6 LP = LPTR(LP)
N1 = ABS(LIST(LP))
IF (N1 .EQ. N0FOR) CNSTR = .TRUE.
IF (N1 .GT. N0) THEN
C
C Draw the edge iff (PASS1=TRUE and CNSTR=FALSE) or
C (PASS1=FALSE and CNSTR=TRUE); i.e., CNSTR and PASS1
C have opposite values.
C
IF (CNSTR .NEQV. PASS1)
. WRITE (LUN,210,ERR=13) X0, Y0, X(N1), Y(N1)
ENDIF
IF (N1 .EQ. N0BAK) CNSTR = .FALSE.
C
C Bottom of loops.
C
IF (LP .NE. LPL) GO TO 6
N0BAK = N0
7 CONTINUE
8 CONTINUE
IF (PASS1) THEN
C
C End of first pass: paint the path and change to dashed
C lines for subsequent drawing. Since the scale factors
C are applied to everything, the dash length must be
C specified in world coordinates.
C
PASS1 = .FALSE.
WRITE (LUN,150,ERR=13)
T = DASHL*2.0/(SFX+SFY)
WRITE (LUN,220,ERR=13) T
220 FORMAT ('[',F12.6,'] 0 setdash')
GO TO 4
ENDIF
C
C Paint the path and restore the saved graphics state (with
C no clip path).
C
WRITE (LUN,150,ERR=13)
WRITE (LUN,230,ERR=13)
230 FORMAT ('grestore')
IF (NUMBR) THEN
C
C Nodes in the window are to be labeled with their indexes.
C Convert FSIZN from points to world coordinates, and
C output the commands to select a font and scale it.
C
T = FSIZN*2.0/(SFX+SFY)
WRITE (LUN,240,ERR=13) T
240 FORMAT ('/Helvetica findfont'/
. F12.6,' scalefont setfont')
C
C Loop on nodes N0 with coordinates (X0,Y0).
C
DO 9 N0 = 1,N
X0 = X(N0)
Y0 = Y(N0)
IF (X0 .LT. WX1 .OR. X0 .GT. WX2 .OR.
. Y0 .LT. WY1 .OR. Y0 .GT. WY2) GO TO 9
C
C Move to (X0,Y0), and draw the label N0. The first char-
C acter will have its lower left corner about one
C character width to the right of the nodal position.
C
WRITE (LUN,180,ERR=13) X0, Y0
WRITE (LUN,250,ERR=13) N0
250 FORMAT ('(',I3,') show')
9 CONTINUE
ENDIF
C
C Convert FSIZT from points to world coordinates, and output
C the commands to select a font and scale it.
C
T = FSIZT*2.0/(SFX+SFY)
WRITE (LUN,240,ERR=13) T
C
C Display TITLE centered above the plot:
C
Y0 = WY2 + 3.0*T
WRITE (LUN,260,ERR=13) TITLE, (WX1+WX2)/2.0, Y0
260 FORMAT (A80/' stringwidth pop 2 div neg ',F12.6,
. ' add ',F12.6,' moveto')
WRITE (LUN,270,ERR=13) TITLE
270 FORMAT (A80/' show')
IF (ANNOT) THEN
C
C Display the window extrema below the plot.
C
X0 = WX1
Y0 = WY1 - 100.0/(SFX+SFY)
WRITE (LUN,180,ERR=13) X0, Y0
WRITE (LUN,280,ERR=13) WX1, WX2
Y0 = Y0 - 2.0*T
WRITE (LUN,290,ERR=13) X0, Y0, WY1, WY2
280 FORMAT ('(Window: WX1 = ',E9.3,', WX2 = ',E9.3,
. ') show')
290 FORMAT ('(Window: ) stringwidth pop ',F12.6,' add',
. F12.6,' moveto'/
. '( WY1 = ',E9.3,', WY2 = ',E9.3,') show')
ENDIF
C
C Paint the path and output the showpage command and
C end-of-file indicator.
C
WRITE (LUN,300,ERR=13)
300 FORMAT ('stroke'/
. 'showpage'/
. '%%EOF')
C
C HP's interpreters require a one-byte End-of-PostScript-Job
C indicator (to eliminate a timeout error message):
C ASCII 4.
C
WRITE (LUN,310,ERR=13) CHAR(4)
310 FORMAT (A1)
C
C No error encountered.
C
IER = 0
RETURN
C
C Invalid input parameter.
C
11 IER = 1
RETURN
C
C DX or DY is not positive.
C
12 IER = 2
RETURN
C
C Error writing to unit LUN.
C
13 IER = 3
RETURN
END
SUBROUTINE TRPRNT (NCC,LCC,N,X,Y,LIST,LPTR,LEND,LOUT,
. PRNTX)
INTEGER NCC, LCC(*), N, LIST(*), LPTR(*), LEND(N),
. LOUT
LOGICAL PRNTX
REAL X(N), Y(N)
C
C***********************************************************
C
C From TRIPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C renka@cs.unt.edu
C 07/30/98
C
C Given a triangulation of a set of points in the plane,
C this subroutine prints the adjacency lists and, option-
C ally, the nodal coordinates on logical unit LOUT. The
C list of neighbors of a boundary node is followed by index
C 0. The numbers of boundary nodes, triangles, and arcs,
C and the constraint curve starting indexes, if any, are
C also printed.
C
C
C On input:
C
C NCC = Number of constraints.
C
C LCC = List of constraint curve starting indexes (or
C dummy array of length 1 if NCC = 0).
C
C N = Number of nodes in the triangulation.
C 3 .LE. N .LE. 9999.
C
C X,Y = Arrays of length N containing the coordinates
C of the nodes in the triangulation -- not used
C unless PRNTX = TRUE.
C
C LIST,LPTR,LEND = Data structure defining the trian-
C gulation. Refer to Subroutine
C TRMESH.
C
C LOUT = Logical unit number for output. 0 .LE. LOUT
C .LE. 99. Output is printed on unit 6 if LOUT
C is outside its valid range on input.
C
C PRNTX = Logical variable with value TRUE if and only
C if X and Y are to be printed (to 6 decimal
C places).
C
C None of the parameters are altered by this routine.
C
C Modules required by TRPRNT: None
C
C***********************************************************
C
INTEGER I, INC, K, LP, LPL, LUN, NA, NABOR(100), NB,
. ND, NL, NLMAX, NMAX, NODE, NN, NT
DATA NMAX/9999/, NLMAX/60/
C
NN = N
LUN = LOUT
IF (LUN .LT. 0 .OR. LUN .GT. 99) LUN = 6
C
C Print a heading and test the range of N.
C
WRITE (LUN,100) NN
IF (NN .LT. 3 .OR. NN .GT. NMAX) THEN
C
C N is outside its valid range.
C
WRITE (LUN,110)
GO TO 5
ENDIF
C
C Initialize NL (the number of lines printed on the current
C page) and NB (the number of boundary nodes encountered).
C
NL = 6
NB = 0
IF (.NOT. PRNTX) THEN
C
C Print LIST only. K is the number of neighbors of NODE
C which are stored in NABOR.
C
WRITE (LUN,101)
DO 2 NODE = 1,NN
LPL = LEND(NODE)
LP = LPL
K = 0
C
1 K = K + 1
LP = LPTR(LP)
ND = LIST(LP)
NABOR(K) = ND
IF (LP .NE. LPL) GO TO 1
IF (ND .LE. 0) THEN
C
C NODE is a boundary node. Correct the sign of the last
C neighbor, add 0 to the end of the list, and increment
C NB.
C
NABOR(K) = -ND
K = K + 1
NABOR(K) = 0
NB = NB + 1
ENDIF
C
C Increment NL and print the list of neighbors.
C
INC = (K-1)/14 + 2
NL = NL + INC
IF (NL .GT. NLMAX) THEN
WRITE (LUN,106)
NL = INC
ENDIF
WRITE (LUN,103) NODE, (NABOR(I), I = 1,K)
IF (K .NE. 14) WRITE (LUN,105)
2 CONTINUE
ELSE
C
C Print X, Y, and LIST.
C
WRITE (LUN,102)
DO 4 NODE = 1,NN
LPL = LEND(NODE)
LP = LPL
K = 0
3 K = K + 1
LP = LPTR(LP)
ND = LIST(LP)
NABOR(K) = ND
IF (LP .NE. LPL) GO TO 3
IF (ND .LE. 0) THEN
C
C NODE is a boundary node.
C
NABOR(K) = -ND
K = K + 1
NABOR(K) = 0
NB = NB + 1
ENDIF
C
C Increment NL and print X, Y, and NABOR.
C
INC = (K-1)/8 + 2
NL = NL + INC
IF (NL .GT. NLMAX) THEN
WRITE (LUN,106)
NL = INC
ENDIF
WRITE (LUN,104) NODE, X(NODE), Y(NODE),
. (NABOR(I), I = 1,K)
IF (K .NE. 8) WRITE (LUN,105)
4 CONTINUE
ENDIF
C
C Print NB, NA, and NT (boundary nodes, arcs, and
C triangles).
C
NT = 2*NN - NB - 2
NA = NT + NN - 1
IF (NL .GT. NLMAX-6) WRITE (LUN,106)
WRITE (LUN,107) NB, NA, NT
C
C Print NCC and LCC.
C
5 WRITE (LUN,108) NCC
IF (NCC .GT. 0) WRITE (LUN,109) (LCC(I), I = 1,NCC)
RETURN
C
C Print formats:
C
100 FORMAT (///,26X,'Adjacency Sets, N = ',I5//)
101 FORMAT (1X,'Node',32X,'Neighbors of Node'//)
102 FORMAT (1X,'Node',5X,'X(Node)',8X,'Y(Node)',
. 20X,'Neighbors of Node'//)
103 FORMAT (1X,I4,5X,14I5/(1X,9X,14I5))
104 FORMAT (1X,I4,2E15.6,5X,8I5/(1X,39X,8I5))
105 FORMAT (1X)
106 FORMAT (///)
107 FORMAT (/1X,'NB = ',I4,' Boundary Nodes',5X,
. 'NA = ',I5,' Arcs',5X,'NT = ',I5,
. ' Triangles')
108 FORMAT (/1X,'NCC =',I3,' Constraint Curves')
109 FORMAT (1X,9X,14I5)
110 FORMAT (1X,10X,'*** N is outside its valid',
. ' range ***')
END