C ALGORITHM 751, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 22, NO. 1, March, 1996, P. 1--8. C ####### With remark from renka (to appear) 4/dec/1998 #! /bin/sh # This is a shell archive, meaning: # 1. Remove everything above the #! /bin/sh line. # 2. Save the resulting text in a file. # 3. Execute the file with /bin/sh (not csh) to create the files: # Drivers/ # Drivers/Fortran77/ # Drivers/Fortran77/Sp/ # Drivers/Fortran77/Sp/RES # Drivers/Fortran77/Sp/RES.eps # Drivers/Fortran77/Sp/driver.f # Src/ # Src/Fortran77/ # Src/Fortran77/Sp/ # Src/Fortran77/Sp/src.f # This archive created: Fri Dec 4 14:11:31 1998 export PATH; PATH=/bin:$PATH if test ! -d 'Drivers' then mkdir 'Drivers' fi cd 'Drivers' if test ! -d 'Fortran77' then mkdir 'Fortran77' fi cd 'Fortran77' if test ! -d 'Sp' then mkdir 'Sp' fi cd 'Sp' if test -f 'RES' then echo shar: will not over-write existing file "'RES'" else cat << SHAR_EOF > 'RES' TRIPACK Test: N = 12 Maximum triangle aspect ratio = 0.480E+00 Adjacency Sets, N = 12 Node X(Node) Y(Node) Neighbors of Node 1 0.000000E+00 0.000000E+00 2 3 9 4 7 0 2 0.100000E+01 0.000000E+00 8 5 10 3 1 0 3 0.500000E+00 0.150000E+00 10 9 1 2 4 0.150000E+00 0.500000E+00 7 1 9 12 5 0.850000E+00 0.500000E+00 8 11 10 2 6 0.500000E+00 0.850000E+00 8 7 12 11 7 0.000000E+00 0.100000E+01 1 4 12 6 8 0 8 0.100000E+01 0.100000E+01 7 6 11 5 2 0 9 0.350000E+00 0.350000E+00 10 11 12 4 1 3 10 0.650000E+00 0.350000E+00 2 5 11 9 3 11 0.650000E+00 0.650000E+00 8 6 12 9 10 5 12 0.350000E+00 0.650000E+00 9 11 6 7 4 NB = 4 Boundary Nodes NA = 29 Arcs NT = 18 Triangles NCC = 1 Constraint Curves 9 TRIPACK (TRLIST) Output Node X(Node) Y(Node) 1 0.000000E+00 0.000000E+00 2 0.100000E+01 0.000000E+00 3 0.500000E+00 0.150000E+00 4 0.150000E+00 0.500000E+00 5 0.850000E+00 0.500000E+00 6 0.500000E+00 0.850000E+00 7 0.000000E+00 0.100000E+01 8 0.100000E+01 0.100000E+01 9 0.350000E+00 0.350000E+00 10 0.650000E+00 0.350000E+00 11 0.650000E+00 0.650000E+00 12 0.350000E+00 0.650000E+00 Triangle Vertices Neighbors Arcs KT N1 N2 N3 KT1 KT2 KT3 KA1 KA2 KA3 1 1 2 3 7 2 0 8 2 1 2 1 3 9 8 3 1 10 3 2 3 1 9 4 9 4 2 12 5 3 4 1 4 7 10 0 3 13 4 5 5 2 8 5 11 6 0 15 7 6 6 2 5 10 12 7 5 16 9 7 7 2 10 3 8 1 6 11 8 9 8 3 10 9 17 2 7 26 10 11 9 4 9 12 18 10 3 28 14 12 10 4 12 7 14 4 9 19 13 14 11 5 8 11 16 12 5 22 17 15 12 5 11 10 17 6 11 25 16 17 13 6 8 7 0 14 16 18 20 23 14 6 7 12 10 15 13 19 21 20 15 6 12 11 18 16 14 27 24 21 16 6 11 8 11 13 15 22 23 24 17 9 10 11 12 18 8 25 29 26 18 9 11 12 15 9 17 27 28 29 NB = 4 Boundary Nodes NA = 29 Arcs NT = 18 Triangles NCC = 1 Constraint Curves 17 A triangulation plot file was successfully created. Output from BNODES and AREAP BNODES: 4 boundary nodes, 29 edges, 18 triangles AREAP: area of convex hull = 0.10000000E+01 No triangulation errors encountered. SHAR_EOF fi # end of overwriting check if test -f 'RES.eps' then echo shar: will not over-write existing file "'RES.eps'" else cat << SHAR_EOF > 'RES.eps' %!PS-Adobe-3.0 EPSF-3.0 %%BoundingBox: 36 126 576 666 %%Title: Triangulation %%Creator: TRIPACK %%EndComments 2.000000 setlinewidth 69 159 moveto 69 633 lineto 543 633 lineto 543 159 lineto closepath stroke 69.000000 159.000000 translate 475.000000 475.000000 scale 0.002105 setlinewidth gsave 0.000000 0.000000 moveto 1.000000 0.000000 lineto 1.000000 1.000000 lineto 0.000000 1.000000 lineto closepath clip newpath 0.000000 0.000000 moveto 1.000000 0.000000 lineto 0.000000 0.000000 moveto 0.500000 0.150000 lineto 0.000000 0.000000 moveto 0.350000 0.350000 lineto 0.000000 0.000000 moveto 0.150000 0.500000 lineto 0.000000 0.000000 moveto 0.000000 1.000000 lineto 1.000000 0.000000 moveto 1.000000 1.000000 lineto 1.000000 0.000000 moveto 0.850000 0.500000 lineto 1.000000 0.000000 moveto 0.650000 0.350000 lineto 1.000000 0.000000 moveto 0.500000 0.150000 lineto 0.500000 0.150000 moveto 0.650000 0.350000 lineto 0.500000 0.150000 moveto 0.350000 0.350000 lineto 0.150000 0.500000 moveto 0.000000 1.000000 lineto 0.150000 0.500000 moveto 0.350000 0.350000 lineto 0.150000 0.500000 moveto 0.350000 0.650000 lineto 0.850000 0.500000 moveto 1.000000 1.000000 lineto 0.850000 0.500000 moveto 0.650000 0.650000 lineto 0.850000 0.500000 moveto 0.650000 0.350000 lineto 0.500000 0.850000 moveto 1.000000 1.000000 lineto 0.500000 0.850000 moveto 0.000000 1.000000 lineto 0.500000 0.850000 moveto 0.350000 0.650000 lineto 0.500000 0.850000 moveto 0.650000 0.650000 lineto 0.000000 1.000000 moveto 0.350000 0.650000 lineto 0.000000 1.000000 moveto 1.000000 1.000000 lineto 1.000000 1.000000 moveto 0.650000 0.650000 lineto stroke [ 0.008421] 0 setdash 0.350000 0.350000 moveto 0.650000 0.350000 lineto 0.350000 0.350000 moveto 0.650000 0.650000 lineto 0.350000 0.350000 moveto 0.350000 0.650000 lineto 0.650000 0.350000 moveto 0.650000 0.650000 lineto 0.650000 0.650000 moveto 0.350000 0.650000 lineto stroke grestore /Helvetica findfont 0.021053 scalefont setfont 0.000000 0.000000 moveto ( 1) show 1.000000 0.000000 moveto ( 2) show 0.500000 0.150000 moveto ( 3) show 0.150000 0.500000 moveto ( 4) show 0.850000 0.500000 moveto ( 5) show 0.500000 0.850000 moveto ( 6) show 0.000000 1.000000 moveto ( 7) show 1.000000 1.000000 moveto ( 8) show 0.350000 0.350000 moveto ( 9) show 0.650000 0.350000 moveto ( 10) show 0.650000 0.650000 moveto ( 11) show 0.350000 0.650000 moveto ( 12) show /Helvetica findfont 0.033684 scalefont setfont (Triangulation created by TRITEST) stringwidth pop 2 div neg 0.500000 add 1.101053 moveto (Triangulation created by TRITEST) show 0.000000 -0.105263 moveto (Window: WX1 = 0.000E+00, WX2 = 0.100E+01) show (Window: ) stringwidth pop 0.000000 add -0.172632 moveto ( WY1 = 0.000E+00, WY2 = 0.100E+01) show stroke showpage %%EOF  SHAR_EOF fi # end of overwriting check if test -f 'driver.f' then echo shar: will not over-write existing file "'driver.f'" else cat << SHAR_EOF > 'driver.f' C C C TRITEST: Portable Test Driver for TRIPACK C 07/02/98 C C C This driver tests software package TRIPACK for construc- C ting a constrained Delaunay triangulation of a set of C points in the plane. All modules other than TRMSHR are C tested unless an error is encountered, in which case the C program terminates immediately. C C By default, tests are performed on a simple data set C consisting of 12 nodes whose convex hull covers the unit C square. The data set includes a single constraint region C consisting of four nodes forming a smaller square at the C center of the unit square. However, by enabling the READ C statements below (C# in the first two columns), testing C may be performed on an arbitrary set of up to NMAX nodes C with up to NCMAX constraint curves. (Refer to the C PARAMETER statements below.) A data set consists of the C following sequence of records: C C N = Number of nodes (format I4) -- 3 to NMAX. C NCC = Number of constraint curves (format I4) -- 0 to C NCMAX. C (LCC(I), I = 1,NCC) = Indexes of the first node in each C constraint curve (format I4). 1 .LE. LCC(1) C and, for I .GT. 1, LCC(I-1) + 3 .LE. LCC(I) C .LE. N-2. (Each constraint curve has at least C three nodes.) C (X(I),Y(I), I = 1,N) = Nodal coordinates with non- C constraint nodes followed by the NCC C sequences of constraint nodes (format C 2F13.8). C C The I/O units may be changed by altering LIN (input) and C LOUT (output) in the DATA statement below. C C This driver must be linked to TRIPACK. C CHARACTER*80 TITLE INTEGER IER, IO1, IO2, K, KSUM, LIN, LNEW, LOUT, LPLT, . LW, LWK, N, N0, N6, NA, NB, NCC, NCMAX, NMAX, . NN, NROW, NT, NTMX INTEGER NEARND LOGICAL NUMBR, PRNTX REAL A, ARMAX, DSQ, PLTSIZ, TOL, WX1, WX2, WY1, WY2 REAL AREAP C PARAMETER (NMAX=100, NCMAX=5, NTMX=2*NMAX, N6=6*NMAX, . LWK=2*NMAX, NROW=9) C C Array storage: C INTEGER LCC(NCMAX), LCT(NCMAX), LEND(NMAX), LIST(N6), . LPTR(N6), LTRI(NROW,NTMX), NODES(LWK) REAL DS(NMAX), X(NMAX), Y(NMAX) C C Tolerance for TRMTST and NEARND: upper bound on squared C distances. C DATA TOL/1.E-2/ C C Plot size for the triangulation plot. C DATA PLTSIZ/7.5/ C C Default data set: C DATA N/12/, NCC/1/, LCC(1)/9/ DATA X(1), X(2), X(3), X(4), X(5), X(6), X(7) . / 0., 1., .5, .15, .85, .5, 0./, . X(8), X(9), X(10), X(11), X(12) . / 1., .35, .65, .65, .35/, . Y(1), Y(2), Y(3), Y(4), Y(5), Y(6), Y(7) . / 0., 0., .15, .5, .5, .85, 1./, . Y(8), Y(9), Y(10), Y(11), Y(12) . / 1., .35, .35, .65, .65/ C C Logical unit numbers for I/O: C DATA LIN/1/, LOUT/2/, LPLT/3/ OPEN (LOUT,FILE='RES') OPEN (LPLT,FILE='RES.eps') C C Store a plot title. It must be enclosed in parentheses. C TITLE = '(Triangulation created by TRITEST)' C C *** Read triangulation parameters -- N, NCC, LCC, X, Y. C C# OPEN (LIN,FILE='tritest.dat',STATUS='OLD') C# READ (LIN,100,ERR=30) N, NCC IF (N .LT. 3 .OR. N .GT. NMAX .OR. NCC .LT. 0 . .OR. NCC .GT. NCMAX) GO TO 31 C# IF (NCC .GT. 0) READ (LIN,100,ERR=30) C# . (LCC(K), K = 1,NCC) C# READ (LIN,110,ERR=30) (X(K),Y(K), K = 1,N) C#100 FORMAT (I4) C#110 FORMAT (2F13.8) C C Print a heading. C WRITE (LOUT,400) N C C *** Create the Delaunay triangulation (TRMESH), and test C for errors (refer to TRMTST below). NODES and DS are C used as work space. C CALL TRMESH (N,X,Y, LIST,LPTR,LEND,LNEW,NODES, . NODES(N+1),DS,IER) IF (IER .EQ. -2) THEN WRITE (LOUT,210) STOP ELSEIF (IER .EQ. -4) THEN WRITE (LOUT,212) STOP ELSEIF (IER .GT. 0) THEN WRITE (LOUT,214) STOP ENDIF CALL TRMTST (N,X,Y,LIST,LPTR,LEND,LNEW,TOL, . LOUT, ARMAX,IER) WRITE (LOUT,410) ARMAX IF (IER .GT. 0) STOP C C *** Add the constraint curves (ADDCST). Note that edges C and triangles are not removed from constraint regions. C ADDCST forces the inclusion of triangulation edges C connecting the sequences of constraint nodes. If it C is necessary to alter the triangulation, the empty C circumcircle property is no longer satisfied. C LW = LWK CALL ADDCST (NCC,LCC,N,X,Y, LW,NODES,LIST,LPTR, . LEND, IER) IF (IER .NE. 0) THEN WRITE (LOUT,220) IER STOP ENDIF IF (LW .EQ. 0) WRITE (LOUT,430) C C *** Test TRPRNT, TRLIST, and TRLPRT, and TRPLOT. C PRNTX = .TRUE. CALL TRPRNT (NCC,LCC,N,X,Y,LIST,LPTR,LEND,LOUT,PRNTX) CALL TRLIST (NCC,LCC,N,LIST,LPTR,LEND,NROW, NT,LTRI, . LCT,IER) CALL TRLPRT (NCC,LCT,N,X,Y,NROW,NT,LTRI,LOUT,PRNTX) C C Set the plot window [WX1,WX2] X [WY1,WY2] to the C smallest rectangle that contains the nodes. C NUMBR = TRUE iff nodal indexes are to be displayed. C WX1 = X(1) WX2 = WX1 WY1 = Y(1) WY2 = WY1 DO 1 K = 2,N IF (X(K) .LT. WX1) WX1 = X(K) IF (X(K) .GT. WX2) WX2 = X(K) IF (Y(K) .LT. WY1) WY1 = Y(K) IF (Y(K) .GT. WY2) WY2 = Y(K) 1 CONTINUE NUMBR = N .LE. 200 CALL TRPLOT (LPLT,PLTSIZ,WX1,WX2,WY1,WY2,NCC,LCC,N, . X,Y,LIST,LPTR,LEND,TITLE,NUMBR, IER) IF (IER .EQ. 0) THEN WRITE (LOUT,470) ELSE WRITE (LOUT,270) IER ENDIF C C *** Test BNODES and AREAP. C CALL BNODES (N,LIST,LPTR,LEND, NODES,NB,NA,NT) A = AREAP(X,Y,NB,NODES) WRITE (LOUT,420) NB, NA, NT, A C C *** Test GETNP by ordering the nodes on distance from N0 C and verifying the ordering. C N0 = N/2 NODES(1) = N0 DS(1) = 0. KSUM = N0 DO 2 K = 2,N CALL GETNP (NCC,LCC,N,X,Y,LIST,LPTR,LEND, . K, NODES,DS, IER) IF (IER .NE. 0 .OR. DS(K) .LT. DS(K-1)) THEN WRITE (LOUT,230) STOP ENDIF KSUM = KSUM + NODES(K) 2 CONTINUE C C Test for all nodal indexes included in NODES. C IF (KSUM .NE. (N*(N+1))/2) THEN WRITE (LOUT,230) STOP ENDIF C C *** Test NEARND by verifying that the nearest node to K is C node K for K = 1 to N. C DO 3 K = 1,N N0 = NEARND (X(K),Y(K),1,N,X,Y,LIST,LPTR,LEND, DSQ) IF (N0 .NE. K .OR. DSQ .GT. TOL) THEN WRITE (LOUT,240) STOP ENDIF 3 CONTINUE C C *** Test DELARC by removing a boundary arc if possible. C The first two nodes define a boundary arc C in the default data set. C IO1 = 1 IO2 = 2 CALL DELARC (N,IO1,IO2, LIST,LPTR,LEND,LNEW, IER) IF (IER .EQ. 1 .OR. IER .EQ. 4) THEN WRITE (LOUT,250) IER STOP ENDIF IF (IER .NE. 0) WRITE (LOUT,440) C C Recreate the triangulation without constraints. C CALL TRMESH (N,X,Y, LIST,LPTR,LEND,LNEW,NODES, . NODES(N+1),DS,IER) NCC = 0 C C *** Test DELNOD by removing nodes 4 to N (in reverse C order). C IF (N .LE. 3) THEN WRITE (LOUT,450) ELSE NN = N 4 LW = LWK/2 CALL DELNOD (NN,NCC, LCC,NN,X,Y,LIST,LPTR,LEND, . LNEW,LW,NODES, IER) IF (IER .NE. 0) THEN WRITE (LOUT,260) IER STOP ENDIF IF (NN .GT. 3) GO TO 4 ENDIF C C Successful test. C WRITE (LOUT,460) STOP C C Error reading the data set. C C# 30 WRITE (*,200) STOP C C Invalid value of N or NCC. C 31 WRITE (*,205) N, NCC STOP C C Error message formats: C C#200 FORMAT (//5X,'*** Input data set invalid ***'/) 205 FORMAT (//5X,'*** N or NCC is outside its valid ', . 'range: N =',I5,', NCC = ',I4,' ***'/) 210 FORMAT (//5X,'*** Error in TRMESH: the first three ', . 'nodes are collinear ***'/) 212 FORMAT (//5X,'*** Error in TRMESH: invalid ', . 'triangulation ***'/) 214 FORMAT (//5X,'*** Error in TRMESH: duplicate nodes ', . 'encountered ***'/) 220 FORMAT (//5X,'*** Error in ADDCST: IER = ',I1, . ' ***'/) 230 FORMAT (//5X,'*** Error in GETNP ***'/) 240 FORMAT (//5X,'*** Error in NEARND ***'/) 250 FORMAT (//5X,'*** Error in DELARC: IER = ',I1, . ' ***'/) 260 FORMAT (//5X,'*** Error in DELNOD: IER = ',I1, . ' ***'/) 270 FORMAT (//5X,'*** Error in TRPLOT: IER = ',I1, . ' ***'/) C C Informative message formats: C 400 FORMAT (///1X,21X,'TRIPACK Test: N =',I5///) 410 FORMAT (5X,'Maximum triangle aspect ratio = ',E10.3//) 420 FORMAT (/5X,'Output from BNODES and AREAP'// . 5X,'BNODES: ',I4,' boundary nodes, ',I5, . ' edges, ',I5,' triangles'/5X, . 'AREAP: area of convex hull = ',E15.8//) 430 FORMAT (5X,'Subroutine EDGE not tested:'/ . 5X,' No edges were swapped by ADDCST'//) 440 FORMAT (5X,'Subroutine DELARC not tested:'/ . 5X,' Nodes 1 and 2 do not form a ', . 'removable boundary edge.'//) 450 FORMAT (5X,'Subroutine DELNOD not tested:'/ . 5X,' N cannot be reduced below 3'//) 460 FORMAT (5X,'No triangulation errors encountered.'/) 470 FORMAT (/5X,'A triangulation plot file was ', . 'successfully created.'/) END SUBROUTINE TRMTST (N,X,Y,LIST,LPTR,LEND,LNEW,TOL, . LUN, ARMAX,IER) INTEGER N, LIST(*), LPTR(*), LEND(N), LNEW, LUN, IER REAL X(N), Y(N), TOL, ARMAX C C*********************************************************** C C From TRPTEST C Robert J. Renka C Dept. of Computer Science C Univ. of North Texas C (817) 565-2767 C 09/01/91 C C This subroutine tests the validity of the data structure C representing a Delaunay triangulation created by subrou- C tine TRMESH. The following properties are tested: C C 1) Each interior node has at least three neighbors, and C each boundary node has at least two neighbors. C C 2) abs(LIST(LP)) is a valid nodal index in the range C 1 to N and LIST(LP) > 0 unless LP = LEND(K) for some C nodal index K. C C 3) Each pointer LEND(K) for K = 1 to N and LPTR(LP) for C LP = 1 to LNEW-1 is a valid LIST index in the range C 1 to LNEW-1. C C 4) N .GE. NB .GE. 3, NT = 2*N-NB-2, and NA = 3*N-NB-3 = C (LNEW-1)/2, where NB, NT, and NA are the numbers of C boundary nodes, triangles, and arcs, respectively. C C 5) Each circumcircle defined by the vertices of a tri- C angle contains no nodes in its interior. This prop- C erty distinguishes a Delaunay triangulation from an C arbitrary triangulation of the nodes. C C Note that no test is made for the property that a triangu- C lation covers the convex hull of the nodes. C C On input: C C N = Number of nodes. N .GE. 3. C C X,Y = Arrays of length N containing the nodal C coordinates. C C LIST,LPTR,LEND = Data structure containing the tri- C angulation. Refer to subroutine C TRMESH. C C TOL = Nonnegative tolerance to allow for floating- C point errors in the circumcircle test. An C error situation is defined as (R**2 - D**2)/ C R**2 > TOL, where R is the radius of a circum- C circle and D is the distance from the C circumcenter to the nearest node. A reason- C able value for TOL is 10*EPS, where EPS is the C machine precision. The test is effectively C bypassed by making TOL large. If TOL < 0, the C tolerance is taken to be 0. C C LUN = Logical unit number for printing error mes- C sages. If LUN < 0 or LUN > 99, no messages C are printed. C C Input parameters are not altered by this routine. C C On output: C C ARMAX = Maximum aspect ratio (radius of inscribed C circle divided by circumradius) of a trian- C gle in the triangulation unless IER > 0. C C IER = Error indicator: C IER = -1 if one or more null triangles (area = C 0) are present but no (other) errors C were encountered. A null triangle is C an error only if it occurs in the C the interior. C IER = 0 if no errors or null triangles were C encountered. C IER = 1 if a node has too few neighbors. C IER = 2 if a LIST entry is outside its valid C range. C IER = 3 if a LPTR or LEND entry is outside its C valid range. C IER = 4 if the triangulation parameters (N, C NB, NT, NA, and LNEW) are inconsistent C (or N < 3 or LNEW is invalid). C IER = 5 if a triangle contains a node interior C to its circumcircle. C C Module required by TRMTST: CIRCUM C C Intrinsic function called by TRMTST: MAX, ABS C C*********************************************************** C INTEGER K, LP, LPL, LPMAX, LPN, N1, N2, N3, NA, NB, . NFAIL, NN, NNB, NT, NULL LOGICAL RATIO, RITE REAL AR, CR, CX, CY, RS, RTOL, SA C C Store local variables, test for errors in input, and C initialize counts. C NN = N LPMAX = LNEW - 1 RTOL = TOL IF (RTOL .LT. 0.) RTOL = 0. RATIO = .TRUE. ARMAX = 0. RITE = LUN .GE. 0 .AND. LUN .LE. 99 IF (NN .LT. 3) GO TO 14 NB = 0 NT = 0 NULL = 0 NFAIL = 0 C C Loop on triangles (N1,N2,N3) such that N2 and N3 index C adjacent neighbors of N1 and are both larger than N1 C (each triangle is associated with its smallest index). C NNB is the neighbor count for N1. C DO 5 N1 = 1,NN NNB = 0 LPL = LEND(N1) IF (LPL .LT. 1 .OR. LPL .GT. LPMAX) THEN LP = LPL GO TO 13 ENDIF LP = LPL C C Loop on neighbors of N1. C 1 LP = LPTR(LP) NNB = NNB + 1 IF (LP .LT. 1 .OR. LP .GT. LPMAX) GO TO 13 N2 = LIST(LP) IF (N2 .LT. 0) THEN IF (LP .NE. LPL) GO TO 12 IF (N2 .EQ. 0 .OR. -N2 .GT. NN) GO TO 12 NB = NB + 1 GO TO 4 ENDIF IF (N2 .LT. 1 .OR. N2 .GT. NN) GO TO 12 LPN = LPTR(LP) N3 = ABS(LIST(LPN)) IF (N2 .LT. N1 .OR. N3 .LT. N1) GO TO 4 NT = NT + 1 C C Compute the coordinates of the circumcenter of C (N1,N2,N3). C CALL CIRCUM (X(N1),Y(N1),X(N2),Y(N2),X(N3),Y(N3), . RATIO, CX,CY,CR,SA,AR) IF (SA .EQ. 0.) THEN NULL = NULL + 1 GO TO 4 ENDIF ARMAX = MAX(ARMAX,AR) C C Test for nodes within the circumcircle. C RS = CR*CR*(1.-RTOL) DO 2 K = 1,NN IF (K .EQ. N1 .OR. K .EQ. N2 .OR. . K .EQ. N3) GO TO 2 IF ((CX-X(K))**2 + (CY-Y(K))**2 .LT. RS) GO TO 3 2 CONTINUE GO TO 4 C C Node K is interior to the circumcircle of (N1,N2,N3). C 3 NFAIL = NFAIL + 1 C C Bottom of loop on neighbors. C 4 IF (LP .NE. LPL) GO TO 1 IF (NNB .LT. 2 .OR. (NNB .EQ. 2 .AND. . LIST(LPL) .GT. 0)) GO TO 11 5 CONTINUE C C Test parameters for consistency and check for NFAIL = 0. C NA = LPMAX/2 IF (NB .LT. 3 .OR. NT .NE. 2*NN-NB-2 .OR. . NA .NE. 3*NN-NB-3) GO TO 14 IF (NFAIL .NE. 0) GO TO 15 C C No errors were encountered. C IER = 0 IF (NULL .EQ. 0) RETURN IER = -1 IF (RITE) WRITE (LUN,100) NULL 100 FORMAT (//5X,'*** TRMTST -- ',I5,' NULL TRIANGLES ', . 'ARE PRESENT'/19X,'(NULL TRIANGLES ', . 'ON THE BOUNDARY ARE UNAVOIDABLE) ***'//) RETURN C C Node N1 has fewer than three neighbors. C 11 IER = 1 IF (RITE) WRITE (LUN,110) N1, NNB 110 FORMAT (//5X,'*** TRMTST -- NODE ',I5, . ' HAS ONLY ',I5,' NEIGHBORS ***'/) RETURN C C N2 = LIST(LP) is outside its valid range. C 12 IER = 2 IF (RITE) WRITE (LUN,120) N2, LP, N1 120 FORMAT (//5X,'*** TRMTST -- LIST(LP) =',I5, . ', FOR LP =',I5,','/19X, .'IS NOT A VALID NEIGHBOR OF ',I5,' ***'/) RETURN C C LIST pointer LP is outside its valid range. C 13 IER = 3 IF (RITE) WRITE (LUN,130) LP, LNEW, N1 130 FORMAT (//5X,'*** TRMTST -- LP =',I5,' IS NOT IN THE', . ' RANGE 1 TO LNEW-1 FOR LNEW = ',I5/ . 19X,'LP POINTS TO A NEIGHBOR OF ',I5, . ' ***'/) RETURN C C Inconsistent triangulation parameters encountered. C 14 IER = 4 IF (RITE) WRITE (LUN,140) N, LNEW, NB, NT, NA 140 FORMAT (//5X,'*** TRMTST -- INCONSISTENT PARAMETERS', . ' ***'/19X,'N = ',I5,' NODES',12X,'LNEW =',I5/ . 19X,'NB = ',I5,' BOUNDARY NODES'/ . 19X,'NT = ',I5,' TRIANGLES'/ . 19X,'NA = ',I5,' ARCS'/) RETURN C C Circumcircle test failure. C 15 IER = 5 IF (RITE) WRITE (LUN,150) NFAIL 150 FORMAT (//5X,'*** TRMTST -- ',I5,' CIRCUMCIRCLES ', . 'CONTAIN NODES IN THEIR INTERIORS ***'/) RETURN END SHAR_EOF fi # end of overwriting check cd .. cd .. cd .. if test ! -d 'Src' then mkdir 'Src' fi cd 'Src' if test ! -d 'Fortran77' then mkdir 'Fortran77' fi cd 'Fortran77' if test ! -d 'Sp' then mkdir 'Sp' fi cd 'Sp' if test -f 'src.f' then echo shar: will not over-write existing file "'src.f'" else cat << SHAR_EOF > 'src.f' 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 SHAR_EOF fi # end of overwriting check cd .. cd .. cd .. # End of shell archive exit 0