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