C ALGORITHM 624 COLLECTED ALGORITHMS FROM ACM. C ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.10, NO. 4, C DEC., 1984, P. 453. SUBROUTINE ADNODE (KK,X,Y, IADJ,IEND, IER) INTEGER KK, IADJ(1), IEND(KK), IER REAL X(KK), Y(KK) LOGICAL SWPTST EXTERNAL INDEX C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C THIS ROUTINE ADDS NODE KK TO A TRIANGULATION OF A SET C OF POINTS IN THE PLANE PRODUCING A NEW TRIANGULATION. A C SEQUENCE OF EDGE SWAPS IS THEN APPLIED TO THE MESH, C RESULTING IN AN OPTIMAL TRIANGULATION. ADNODE IS PART C OF AN INTERPOLATION PACKAGE WHICH ALSO PROVIDES ROUTINES C TO INITIALIZE THE DATA STRUCTURE, PLOT THE MESH, AND C DELETE ARCS. C C INPUT PARAMETERS - KK - INDEX OF THE NODE TO BE ADDED C TO THE MESH. KK .GE. 4. C C X,Y - VECTORS OF COORDINATES OF THE C NODES IN THE MESH. (X(I),Y(I)) C DEFINES NODE I FOR I = 1,..,KK. C C IADJ - SET OF ADJACENCY LISTS OF NODES C 1,..,KK-1. C C IEND - POINTERS TO THE ENDS OF C ADJACENCY LISTS IN IADJ FOR C EACH NODE IN THE MESH. C C IADJ AND IEND MAY BE CREATED BY TRMESH. C C KK, X, AND Y ARE NOT ALTERED BY THIS ROUTINE. C C OUTPUT PARAMETERS - IADJ,IEND - UPDATED WITH THE ADDITION C OF NODE KK AS THE LAST C ENTRY. C C IER - ERROR INDICATOR C IER = 0 IF NO ERRORS C WERE ENCOUNTERED. C IER = 1 IF ALL NODES C (INCLUDING KK) ARE C COLLINEAR. C C MODULES REFERENCED BY ADNODE - TRFIND, INTADD, BDYADD, C SHIFTD, INDEX, SWPTST, C SWAP C C*********************************************************** C INTEGER K, KM1, I1, I2, I3, INDKF, INDKL, NABOR1, . IO1, IO2, IN1, INDK1, IND2F, IND21 REAL XK, YK C C LOCAL PARAMETERS - C C K = LOCAL COPY OF KK C KM1 = K - 1 C I1,I2,I3 = VERTICES OF A TRIANGLE CONTAINING K C INDKF = IADJ INDEX OF THE FIRST NEIGHBOR OF K C INDKL = IADJ INDEX OF THE LAST NEIGHBOR OF K C NABOR1 = FIRST NEIGHBOR OF K BEFORE ANY SWAPS OCCUR C IO1,IO2 = ADJACENT NEIGHBORS OF K DEFINING AN ARC TO C BE TESTED FOR A SWAP C IN1 = VERTEX OPPOSITE K -- FIRST NEIGHBOR OF IO2 C WHICH PRECEDES IO1. IN1,IO1,IO2 ARE IN C COUNTERCLOCKWISE ORDER. C INDK1 = INDEX OF IO1 IN THE ADJACENCY LIST FOR K C IND2F = INDEX OF THE FIRST NEIGHBOR OF IO2 C IND21 = INDEX OF IO1 IN THE ADJACENCY LIST FOR IO2 C XK,YK = X(K), Y(K) C IER = 0 K = KK C C INITIALIZATION C KM1 = K - 1 XK = X(K) YK = Y(K) C C ADD NODE K TO THE MESH C CALL TRFIND(KM1,XK,YK,X,Y,IADJ,IEND, I1,I2,I3) IF (I1 .EQ. 0) GO TO 5 IF (I3 .EQ. 0) CALL BDYADD(K,I1,I2, IADJ,IEND ) IF (I3 .NE. 0) CALL INTADD(K,I1,I2,I3, IADJ,IEND ) C C INITIALIZE VARIABLES FOR OPTIMIZATION OF THE MESH C INDKF = IEND(KM1) + 1 INDKL = IEND(K) NABOR1 = IADJ(INDKF) IO2 = NABOR1 INDK1 = INDKF + 1 IO1 = IADJ(INDK1) C C BEGIN LOOP -- FIND THE VERTEX OPPOSITE K C 1 IND2F = 1 IF (IO2 .NE. 1) IND2F = IEND(IO2-1) + 1 IND21 = INDEX(IO2,IO1,IADJ,IEND) IF (IND2F .EQ. IND21) GO TO 2 IN1 = IADJ(IND21-1) GO TO 3 C C IN1 IS THE LAST NEIGHBOR OF IO2 C 2 IND21 = IEND(IO2) IN1 = IADJ(IND21) IF (IN1 .EQ. 0) GO TO 4 C C SWAP TEST -- IF A SWAP OCCURS, TWO NEW ARCS ARE OPPOSITE K C AND MUST BE TESTED. INDK1 AND INDKF MUST BE C DECREMENTED. C 3 IF ( .NOT. SWPTST(IN1,K,IO1,IO2,X,Y) ) GO TO 4 CALL SWAP(IN1,K,IO1,IO2, IADJ,IEND ) IO1 = IN1 INDK1 = INDK1 - 1 INDKF = INDKF - 1 GO TO 1 C C NO SWAP OCCURRED. RESET IO2 AND IO1, AND TEST FOR C TERMINATION. C 4 IF (IO1 .EQ. NABOR1) RETURN IO2 = IO1 INDK1 = INDK1 + 1 IF (INDK1 .GT. INDKL) INDK1 = INDKF IO1 = IADJ(INDK1) IF (IO1 .NE. 0) GO TO 1 RETURN C C ALL NODES ARE COLLINEAR C 5 IER = 1 RETURN END FUNCTION AREA (X,Y,NB,NODES) INTEGER NB, NODES(NB) REAL X(1), Y(1) C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C GIVEN A SEQUENCE OF NB POINTS IN THE PLANE, THIS C FUNCTION COMPUTES THE 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 THUS THIS POINT NEED NOT BE SPECIFIED TWICE. HOWEVER, ANY C POINT MAY BE SPECIFIED MORE THAN ONCE IN ORDER TO DEFINE A C MULTIPLY CONNECTED DOMAIN. C THE AREA OF A TRIANGULATION MAY BE COMPUTED BY CALLING C AREA WITH VALUES OF NB AND NODES DETERMINED BY SUBROUTINE C BNODES. C C INPUT PARAMETERS - X,Y - N-VECTORS OF COORDINATES OF C POINTS IN THE PLANE FOR N .GE. C NB. NODE I HAS COORDINATES C (X(I),Y(I)) FOR I = 1, 2, ..., C N. C C NB - LENGTH OF NODES. C C NODES - VECTOR OF NODE INDICES IN THE C RANGE 1 TO N DEFINING THE C POLYGONAL CURVE. C C INPUT PARAMETERS ARE NOT ALTERED BY THIS FUNCTION. C C OUTPUT PARAMETER - AREA - SIGNED AREA BOUNDED BY THE C POLYGONAL CURVE DEFINED C ABOVE. C C MODULES REFERENCED BY AREA - NONE C C*********************************************************** C INTEGER NNB, ND, I REAL A, X0, Y0, DX1, DY1, DX2, DY2 C C LOCAL PARAMETERS - C C NNB = LOCAL COPY OF NB C ND = ELEMENT OF NODES C I = DO-LOOP AND NODES INDEX C A = PARTIAL SUM OF SIGNED (AND DOUBLED) TRIANGLE C AREAS C X0,Y0 = X(NODES(1)), Y(NODES(1)) C DX1,DY1 = COMPONENTS OF THE VECTOR NODES(1)-NODES(I) FOR C I = 2, 3, ..., NB-1 C DX2,DY2 = COMPONENTS OF THE VECTOR NODES(1)-NODES(I) FOR C I = 3, 4, ..., NB C NNB = NB A = 0. IF (NNB .LT. 3) GO TO 2 C C INITIALIZATION C ND = NODES(1) X0 = X(ND) Y0 = Y(ND) ND = NODES(2) DX1 = X(ND) - X0 DY1 = Y(ND) - Y0 C C LOOP ON TRIANGLES (NODES(1),NODES(I),NODES(I+1)), C I = 2, 3, ..., NB-1, ADDING TWICE THEIR SIGNED C AREAS TO A C DO 1 I = 3,NNB ND = NODES(I) DX2 = X(ND) - X0 DY2 = Y(ND) - Y0 A = A + DX1*DY2 - DX2*DY1 DX1 = DX2 DY1 = DY2 1 CONTINUE C C A CONTAINS TWICE THE SIGNED AREA OF THE REGION C 2 AREA = A/2. RETURN END SUBROUTINE BDYADD (KK,I1,I2, IADJ,IEND ) INTEGER KK, I1, I2, IADJ(1), IEND(KK) C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C THIS ROUTINE ADDS A BOUNDARY NODE TO A TRIANGULATION C OF A SET OF KK-1 POINTS IN THE PLANE. IADJ AND IEND ARE C UPDATED WITH THE INSERTION OF NODE KK. C C INPUT PARAMETERS - KK - INDEX OF AN EXTERIOR NODE TO BE C ADDED. KK .GE. 4. C C I1 - FIRST (RIGHTMOST AS VIEWED FROM C KK) BOUNDARY NODE IN THE MESH C WHICH IS VISIBLE FROM KK - THE C LINE SEGMENT KK-I1 INTERSECTS C NO ARCS. C C I2 - LAST (LEFTMOST) BOUNDARY NODE C WHICH IS VISIBLE FROM KK. C C IADJ - SET OF ADJACENCY LISTS OF NODES C IN THE MESH. C C IEND - POINTERS TO THE ENDS OF C ADJACENCY LISTS IN IADJ FOR C EACH NODE IN THE MESH. C C IADJ AND IEND MAY BE CREATED BY TRMESH AND MUST CONTAIN C THE VERTICES I1 AND I2. I1 AND I2 MAY BE DETERMINED BY C TRFIND. C C KK, I1, AND I2 ARE NOT ALTERED BY THIS ROUTINE. C C OUTPUT PARAMETERS - IADJ,IEND - UPDATED WITH THE ADDITION C OF NODE KK AS THE LAST C ENTRY. NODE KK WILL BE C CONNECTED TO I1, I2, AND C ALL BOUNDARY NODES BETWEEN C THEM. NO OPTIMIZATION OF C THE MESH IS PERFORMED. C C MODULE REFERENCED BY BDYADD - SHIFTD C C INTRINSIC FUNCTIONS CALLED BY BDYADD - MIN0, MAX0 C C*********************************************************** C INTEGER K, KM1, NRIGHT, NLEFT, NF, NL, N1, N2, I, . IMIN, IMAX, KEND, NEXT, INDX C C LOCAL PARAMETERS - C C K = LOCAL COPY OF KK C KM1 = K - 1 C NRIGHT,NLEFT = LOCAL COPIES OF I1, I2 C NF,NL = INDICES OF IADJ BOUNDING THE PORTION OF THE C ARRAY TO BE SHIFTED C N1 = IADJ INDEX OF THE FIRST NEIGHBOR OF NLEFT C N2 = IADJ INDEX OF THE LAST NEIGHBOR OF NRIGHT C I = DO-LOOP INDEX C IMIN,IMAX = BOUNDS ON DO-LOOP INDEX -- FIRST AND LAST C ELEMENTS OF IEND TO BE INCREMENTED C KEND = POINTER TO THE LAST NEIGHBOR OF K IN IADJ C NEXT = NEXT BOUNDARY NODE TO BE CONNECTED TO KK C INDX = INDEX FOR IADJ C K = KK KM1 = K - 1 NRIGHT = I1 NLEFT = I2 C C INITIALIZE VARIABLES C NL = IEND(KM1) N1 = 1 IF (NLEFT .NE. 1) N1 = IEND(NLEFT-1) + 1 N2 = IEND(NRIGHT) NF = MAX0(N1,N2) C C INSERT K AS A NEIGHBOR OF MAX(NRIGHT,NLEFT) C CALL SHIFTD(NF,NL,2, IADJ ) IADJ(NF+1) = K IMIN = MAX0(NRIGHT,NLEFT) DO 1 I = IMIN,KM1 IEND(I) = IEND(I) + 2 1 CONTINUE C C INITIALIZE KEND AND INSERT K AS A NEIGHBOR OF C MIN(NRIGHT,NLEFT) C KEND = NL + 3 NL = NF - 1 NF = MIN0(N1,N2) CALL SHIFTD(NF,NL,1, IADJ ) IADJ(NF) = K IMAX = IMIN - 1 IMIN = MIN0(NRIGHT,NLEFT) DO 2 I = IMIN,IMAX IEND(I) = IEND(I) + 1 2 CONTINUE C C INSERT NRIGHT AS THE FIRST NEIGHBOR OF K C IADJ(KEND) = NRIGHT C C INITIALIZE INDX FOR LOOP ON BOUNDARY NODES BETWEEN NRIGHT C AND NLEFT C INDX = IEND(NRIGHT) - 2 3 NEXT = IADJ(INDX) IF (NEXT .EQ. NLEFT) GO TO 4 C C CONNECT NEXT AND K C KEND = KEND + 1 IADJ(KEND) = NEXT INDX = IEND(NEXT) IADJ(INDX) = K INDX = INDX - 1 GO TO 3 C C INSERT NLEFT AND 0 AS THE LAST NEIGHBORS OF K C 4 IADJ(KEND+1) = NLEFT KEND = KEND + 2 IADJ(KEND) = 0 IEND(K) = KEND RETURN END SUBROUTINE BNODES (N,IADJ,IEND, NB,NA,NT,NODES) INTEGER N, IADJ(1), IEND(N), NB, NA, NT, NODES(1) C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C GIVEN A TRIANGULATION OF N POINTS IN THE PLANE, THIS C ROUTINE RETURNS A VECTOR CONTAINING THE INDICES, IN C COUNTERCLOCKWISE ORDER, OF THE NODES ON THE BOUNDARY OF C THE CONVEX HULL OF THE SET OF POINTS. C C INPUT PARAMETERS - N - NUMBER OF NODES IN THE MESH. C C IADJ - SET OF ADJACENCY LISTS OF C NODES IN THE MESH. C C IEND - POINTERS TO THE ENDS OF C ADJACENCY LISTS IN IADJ FOR C EACH NODE IN THE MESH. C C NODES - VECTOR OF LENGTH .GE. NB. C (NB .LE. N). C C IADJ AND IEND MAY BE CREATED BY TRMESH AND ARE NOT C ALTERED BY THIS ROUTINE. C C OUTPUT PARAMETERS - NB - NUMBER OF BOUNDARY NODES. C C NA,NT - NUMBER OF ARCS AND TRIANGLES, C RESPECTIVELY, IN THE MESH. C C NODES - VECTOR OF NB BOUNDARY NODE C INDICES RANGING FROM 1 TO N. C C MODULES REFERENCED BY BNODES - NONE C C*********************************************************** C INTEGER NST, INDL, K, N0, INDF C C LOCAL PARAMETERS - C C NST = FIRST ELEMENT OF NODES -- ARBITRARILY CHOSEN C INDL = IADJ INDEX OF THE LAST NEIGHBOR OF NST C K = NODES INDEX C N0 = BOUNDARY NODE TO BE ADDED TO NODES C INDF = IADJ INDEX OF THE FIRST NEIGHBOR OF N0 C C SET NST TO THE FIRST BOUNDARY NODE ENCOUNTERED C NST = 1 1 INDL = IEND(NST) IF (IADJ(INDL) .EQ. 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 INDF = 1 IF (N0 .GT. 1) INDF = IEND(N0-1) + 1 N0 = IADJ(INDF) 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 DELETE (NN,NOUT1,NOUT2, IADJ,IEND, IER) INTEGER NN, NOUT1, NOUT2, IADJ(1), IEND(NN), IER EXTERNAL INDEX C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C THIS ROUTINE DELETES A BOUNDARY EDGE FROM A TRIANGU- C LATION OF A SET OF POINTS IN THE PLANE. IT MAY BE NEC- C ESSARY TO FORCE CERTAIN EDGES TO BE PRESENT BEFORE CALL- C ING DELETE (SEE SUBROUTINE EDGE). NOTE THAT SUBROUTINES C EDGE, TRFIND, AND THE ROUTINES WHICH CALL TRFIND (ADNODE, C UNIF, INTRC1, AND INTRC0) SHOULD NOT BE CALLED FOLLOWING C A DELETION. C C INPUT PARAMETERS - NN - NUMBER OF NODES IN THE TRIAN- C GULATION. C C NOUT1,NOUT2 - PAIR OF ADJACENT NODES ON THE C BOUNDARY DEFINING THE ARC TO C BE REMOVED. NOUT2 MUST BE THE C LAST NONZERO NEIGHBOR OF NOUT1. C C THE ABOVE PARAMETERS ARE NOT ALTERED BY THIS ROUTINE. C C IADJ,IEND - DATA STRUCTURE DEFINING THE C TRIANGULATION (SEE SUBROUTINE C TRMESH). C C OUTPUT PARAMETERS - IADJ,IEND - UPDATED WITH THE REMOVAL C OF THE ARC NOUT1-NOUT2 C IF IER .EQ. 0. C C IER - ERROR INDICATOR C IER = 0 IF NO ERRORS WERE C ENCOUNTERED. C IER = 1 IF NOUT1 OR NOUT2 C IS NOT ON THE C BOUNDARY. C IER = 2 IF NOUT1 OR NOUT2 C HAS ONLY 2 NONZERO C NEIGHBORS. C IER = 3 IF NOUT2 IS NOT C THE LAST NEIGHBOR C OF NOUT1. C IER = 4 IF A DELETION C WOULD DIVIDE THE C MESH INTO TWO C REGIONS. C C MODULES REFERENCED BY DELETE - SHIFTD, INDEX C C*********************************************************** C INTEGER N, IOUT1, IOUT2, IO1, IO2, IND12, IND21, . ITEMP, IND1F, IND1L, IND2F, IND2L, NEWBD, . INDNF, INDNL, INDN0, INDFP2, INDLM3, NF, NL, . I, IMAX C C LOCAL PARAMETERS - C C N = LOCAL COPY OF NN C IOUT1,IOUT2 = LOCAL COPIES OF NOUT1 AND NOUT2 C IO1,IO2 = NOUT1,NOUT2 IN ORDER OF INCREASING MAGNITUDE C IND12 = INDEX OF IO2 IN THE ADJACENCY LIST FOR IO1 C IND21 = INDEX OF IO1 IN THE ADJACENCY LIST FOR IO2 C ITEMP = TEMPORARY STORAGE LOCATION FOR PERMUTATIONS C IND1F = IADJ INDEX OF THE FIRST NEIGHBOR OF IO1 C IND1L = IADJ INDEX OF THE LAST NEIGHBOR OF IO1 C IND2F = IADJ INDEX OF THE FIRST NEIGHBOR OF IO2 C IND2L = IADJ INDEX OF THE LAST NEIGHBOR OF IO2 C NEWBD = THE NEIGHBOR COMMON TO NOUT1 AND NOUT2 C INDNF = IADJ INDEX OF THE FIRST NEIGHBOR OF NEWBD C INDNL = IADJ INDEX OF THE LAST NEIGHBOR OF NEWBD C INDN0 = INDEX OF 0 IN THE ADJACENCY LIST FOR NEWBD C BEFORE PERMUTING THE NEIGHBORS C INDFP2 = INDNF + 2 C INDLM3 = INDNL - 3 C NF,NL = BOUNDS ON THE PORTION OF IADJ TO BE SHIFTED C I = DO-LOOP INDEX C IMAX = UPPER BOUND ON DO-LOOP FOR SHIFTING IEND C N = NN IOUT1 = NOUT1 IOUT2 = NOUT2 C C INITIALIZE INDICES C IND1F = 1 IF (IOUT1 .GT. 1) IND1F = IEND(IOUT1-1) + 1 IND1L = IEND(IOUT1) IND2F = 1 IF (IOUT2 .GT. 1) IND2F = IEND(IOUT2-1) + 1 IND2L = IEND(IOUT2) NEWBD = IADJ(IND1L-2) INDN0 = INDEX(NEWBD,IOUT2,IADJ,IEND) INDNL = IEND(NEWBD) C C ORDER VERTICES SUCH THAT THE NEIGHBORS OF IO1 PRECEDE C THOSE OF IO2 C IF (IOUT1 .GT. IOUT2) GO TO 1 IO1 = IOUT1 IO2 = IOUT2 IND12 = IND1L - 1 IND21 = IND2F GO TO 2 1 IO1 = IOUT2 IO2 = IOUT1 IND12 = IND2F IND21 = IND1L - 1 C C CHECK FOR ERRORS C 2 IF ( (IADJ(IND1L) .NE. 0) .OR. (IADJ(IND2L) .NE. 0) ) . GO TO 21 IF ( (IND1L-IND1F .LE. 2) .OR. (IND2L-IND2F .LE. 2) ) . GO TO 22 IF (IADJ(IND1L-1) .NE. IOUT2) GO TO 23 IF (IADJ(INDNL) .EQ. 0) GO TO 24 C C DELETE THE EDGE IO1-IO2 AND MAKE NEWBD A BOUNDARY NODE C IF (NEWBD .LT. IO1) GO TO 8 IF (NEWBD .LT. IO2) GO TO 6 C C THE VERTICES ARE ORDERED IO1, IO2, NEWBD. C DELETE IO2 AS A NEIGHBOR OF IO1. C NF = IND12 + 1 NL = IND21 - 1 CALL SHIFTD(NF,NL,-1, IADJ ) IMAX = IO2 - 1 DO 3 I = IO1,IMAX IEND(I) = IEND(I) - 1 3 CONTINUE C C DELETE IO1 AS A NEIGHBOR OF IO2 C NF = NL + 2 NL = INDN0 CALL SHIFTD(NF,NL,-2, IADJ ) IMAX = NEWBD - 1 DO 4 I = IO2,IMAX IEND(I) = IEND(I) - 2 4 CONTINUE C C SHIFT THE BOTTOM OF IADJ UP 1 LEAVING ROOM FOR 0 AS A C NEIGHBOR OF NEWBD C INDN0 = INDN0 - 1 NF = NL + 1 NL = IEND(N) IF (NF .LE. NL) CALL SHIFTD(NF,NL,-1, IADJ ) DO 5 I = NEWBD,N IEND(I) = IEND(I) - 1 5 CONTINUE GO TO 12 C C THE VERTICES ARE ORDERED IO1, NEWBD, IO2. C DELETE IO2 AS A NEIGHBOR OF IO1 LEAVING ROOM FOR 0 AS A C NEIGHBOR OF NEWBD. C 6 NF = IND12 + 1 NL = INDN0 CALL SHIFTD(NF,NL,-1, IADJ ) IMAX = NEWBD - 1 DO 7 I = IO1,IMAX IEND(I) = IEND(I) - 1 7 CONTINUE GO TO 10 C C THE VERTICES ARE ORDERED NEWBD, IO1, IO2. C DELETE IO2 AS A NEIGHBOR OF IO1 LEAVING ROOM FOR 0 AS A C NEIGHBOR OF NEWBD. C 8 INDN0 = INDN0 + 1 NF = INDN0 NL = IND12 - 1 IF (NF .LE. NL) CALL SHIFTD(NF,NL,1, IADJ ) IMAX = IO1 - 1 DO 9 I = NEWBD,IMAX IEND(I) = IEND(I) + 1 9 CONTINUE C C DELETE IO1 AS A NEIGHBOR OF IO2 C 10 NF = IND21 + 1 NL = IEND(N) CALL SHIFTD(NF,NL,-1, IADJ ) DO 11 I = IO2,N IEND(I) = IEND(I) - 1 11 CONTINUE C C PERMUTE THE NEIGHBORS OF NEWBD WITH END-AROUND SHIFTS SO C THAT 0 IS THE LAST NEIGHBOR C 12 INDNF = 1 IF (NEWBD .GT. 1) INDNF = IEND(NEWBD-1) + 1 INDNL = IEND(NEWBD) IF (INDN0-INDNF .GE. INDNL-INDN0) GO TO 16 C C SHIFT UPWARD C IF (INDN0 .GT. INDNF) GO TO 13 CALL SHIFTD(INDNF+1,INDNL,-1, IADJ ) GO TO 20 13 INDFP2 = INDNF + 2 IF (INDN0 .LT. INDFP2) GO TO 15 DO 14 I = INDFP2,INDN0 ITEMP = IADJ(INDNF) CALL SHIFTD(INDNF+1,INDNL,-1, IADJ ) IADJ(INDNL) = ITEMP 14 CONTINUE C C THE LAST SHIFT IS BY 2 C 15 ITEMP = IADJ(INDNF) CALL SHIFTD(INDFP2,INDNL,-2, IADJ ) IADJ(INDNL-1) = ITEMP GO TO 20 C C SHIFT DOWNWARD C 16 IF (INDN0 .EQ. INDNL) GO TO 20 IF (INDN0 .LT. INDNL-1) GO TO 17 CALL SHIFTD(INDNF,INDNL-2,1, IADJ ) IADJ(INDNF) = IADJ(INDNL) GO TO 20 17 INDLM3 = INDNL - 3 IF (INDN0 .GT. INDLM3) GO TO 19 DO 18 I = INDN0,INDLM3 ITEMP = IADJ(INDNL) CALL SHIFTD(INDNF,INDNL-1,1, IADJ ) IADJ(INDNF) = ITEMP 18 CONTINUE C C THE LAST SHIFT IS BY 2 C 19 ITEMP = IADJ(INDNL-1) CALL SHIFTD(INDNF,INDLM3,2, IADJ ) IADJ(INDNF+1) = IADJ(INDNL) IADJ(INDNF) = ITEMP C C INSERT 0 AS THE LAST NEIGHBOR OF NEWBD C 20 IADJ(INDNL) = 0 IER = 0 RETURN C C ONE OF THE VERTICES IS NOT ON THE BOUNDARY C 21 IER = 1 RETURN C C ONE OF THE VERTICES HAS ONLY TWO NONZERO NEIGHBORS. THE C TRIANGULATION WOULD BE DESTROYED BY A DELETION C 22 IER = 2 RETURN C C NOUT2 IS NOT THE LAST NONZERO NEIGHBOR OF NOUT1 C 23 IER = 3 RETURN C C A DELETION WOULD DIVIDE THE MESH INTO TWO REGIONS C CONNECTED AT A SINGLE NODE C 24 IER = 4 RETURN END SUBROUTINE EDGE (IN1,IN2,X,Y, LWK,IWK,IADJ,IEND, IER) LOGICAL SWPTST INTEGER IN1, IN2, LWK, IWK(2,1), IADJ(1), IEND(1), IER REAL X(1), Y(1) C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C GIVEN A TRIANGULATION OF N NODES AND A PAIR OF NODAL C INDICES 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 THIESSEN TRIANGU- C LATION IS INPUT, THE RESULTING TRIANGULATION IS AS CLOSE C AS POSSIBLE TO A THIESSEN TRIANGULATION IN THE SENSE THAT C ALL ARCS OTHER THAN IN1-IN2 ARE LOCALLY OPTIMAL. 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 REGION. SUBSEQUENT DELETION OF EDGES OUTSIDE THIS C REGION (BY SUBROUTINE DELETE) RESULTS IN A NONCONVEX TRI- C ANGULATION WHICH MAY SERVE AS A FINITE ELEMENT GRID. C (EDGE SHOULD NOT BE CALLED AFTER A CALL TO DELETE.) IF, C ON THE OTHER HAND, INTERPOLATION IS TO BE PERFORMED IN THE C NONCONVEX REGION, EDGES MUST NOT BE DELETED, BUT IT IS C STILL ADVANTAGEOUS TO HAVE THE NONCONVEX BOUNDARY PRESENT C IF IT IS DESIRABLE THAT INTERPOLATED VALUES BE INFLUENCED C BY THE GEOMETRY. NOTE THAT SUBROUTINE GETNP WHICH IS USED C TO SELECT THE NODES ENTERING INTO LOCAL DERIVATIVE ESTI- C MATES WILL NOT NECESSARILY RETURN CLOSEST NODES IF THE C TRIANGULATION HAS BEEN RENDERED NONOPTIMAL BY A CALL TO C EDGE. HOWEVER, THE EFFECT WILL BE MERELY TO FURTHER EN- C HANCE THE INFLUENCE OF THE NONCONVEX GEOMETRY ON INTERPO- C LATED VALUES. C C INPUT PARAMETERS - IN1,IN2 - INDICES (OF X AND Y) IN THE C RANGE 1,...,N DEFINING A PAIR C OF NODES TO BE CONNECTED BY C AN ARC. C C X,Y - N-VECTORS CONTAINING CARTE- C SIAN COORDINATES OF THE C NODES. C C THE ABOVE PARAMETERS ARE NOT ALTERED BY THIS ROUTINE. C C LWK - NUMBER OF COLUMNS RESERVED C FOR IWK. THIS MUST BE AT C LEAST NI -- THE NUMBER OF C ARCS WHICH INTERSECT IN1-IN2. C (NI IS BOUNDED BY N-3). C C IWK - INTEGER WORK ARRAY DIMENSION- C ED 2 BY LWK (OR VECTOR OF C LENGTH .GE. 2*LWK). C C IADJ,IEND - DATA STRUCTURE DEFINING THE C TRIANGULATION. SEE SUBROU- C TINE TRMESH. C C OUTPUT PARAMETERS - LWK - NUMBER OF IWK COLUMNS REQUIRED C IF IER = 0 OR IER = 2. LWK = 0 C IFF IN1 AND IN2 WERE ADJACENT C ON INPUT. C C IWK - CONTAINS THE INDICES OF THE END- C POINTS OF THE NEW ARCS OTHER C THAN IN1-IN2 UNLESS IER .GT. 0 C OR LWK = 0. NEW ARCS TO THE C LEFT OF IN1->IN2 ARE STORED IN C THE FIRST K-1 COLUMNS (LEFT POR- C TION OF IWK), COLUMN K CONTAINS C ZEROS, AND NEW ARCS TO THE RIGHT C OF IN1->IN2 OCCUPY COLUMNS K+1, C ...,LWK. (K CAN BE DETERMINED C BY SEARCHING IWK FOR THE ZEROS.) C C IADJ,IEND - UPDATED IF NECESSARY TO REFLECT C THE PRESENCE OF AN ARC CONNECT- C ING IN1 AND IN2, UNALTERED IF C IER .NE. 0. C C IER - ERROR INDICATOR C IER = 0 IF NO ERRORS WERE EN- C COUNTERED. C IER = 1 IF IN1 .LT. 1, IN2 .LT. C 1, IN1 = IN2, OR LWK C .LT. 0 ON INPUT. C IER = 2 IF MORE SPACE IS REQUIR- C ED IN IWK. SEE LWK. C IER = 3 IF IN1 AND IN2 COULD NOT C BE CONNECTED DUE TO AN C INVALID DATA STRUCTURE. C C MODULES REFERENCED BY EDGE - SWAP, INDEX, SHIFTD, SWPTST C C*********************************************************** C INTEGER N1, N2, IWEND, IWL, INDF, INDX, N1LST, NL, NR, . NEXT, IWF, LFT, N0, IWC, IWCP1, IWCM1, I, IO1, . IO2, INDL REAL X1, Y1, X2, Y2, X0, Y0 LOGICAL SWP, LEFT C C LOCAL PARAMETERS - C C N1,N2 = LOCAL COPIES OF IN1 AND IN2 OR NODES OPPOSITE AN C ARC IO1-IO2 TO BE TESTED FOR A SWAP IN THE C OPTIMIZATION LOOPS C IWEND = INPUT OR OUTPUT VALUE OF LWK C IWL = IWK (COLUMN) INDEX OF THE LAST (RIGHTMOST) ARC C WHICH INTERSECTS IN1->IN2 C INDF = IADJ INDEX OF THE FIRST NEIGHBOR OF IN1 OR IO1 C INDX = IADJ INDEX OF A NEIGHBOR OF IN1, NL, OR IO1 C N1LST = LAST NEIGHBOR OF IN1 C NL,NR = ENDPOINTS OF AN ARC WHICH INTERSECTS IN1-IN2 C WITH NL LEFT IN1->IN2 C NEXT = NODE OPPOSITE NL->NR C IWF = IWK (COLUMN) INDEX OF THE FIRST (LEFTMOST) ARC 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 N0 = NODE OPPOSITE NR->NL 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 IWCM1 = IWC - 1 C I = DO-LOOP INDEX AND COLUMN INDEX FOR IWK C IO1,IO2 = ENDPOINTS OF AN ARC TO BE TESTED FOR A SWAP IN C THE OPTIMIZATION LOOPS C INDL = IADJ INDEX OF THE LAST NEIGHBOR OF IO1 C X1,Y1 = COORDINATES OF IN1 C X2,Y2 = COORDINATES OF IN2 C X0,Y0 = COORDINATES OF N0 C SWP = FLAG SET TO .TRUE. IFF A SWAP OCCURS IN AN OPTI- C MIZATION LOOP C LEFT = STATEMENT FUNCTION WHICH RETURNS THE VALUE C .TRUE. IFF (XP,YP) IS ON OR TO THE LEFT OF THE C VECTOR (XA,YA)->(XB,YB) C LEFT(XA,YA,XB,YB,XP,YP) = (XB-XA)*(YP-YA) .GE. . (XP-XA)*(YB-YA) C C STORE IN1, IN2, AND LWK IN LOCAL VARIABLES AND CHECK 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 35 C C STORE THE COORDINATES OF N1 AND N2 AND INITIALIZE IWL. C X1 = X(N1) Y1 = Y(N1) X2 = X(N2) Y2 = Y(N2) IWL = 0 C C SET NR AND NL TO ADJACENT NEIGHBORS OF N1 SUCH THAT C NR LEFT N2->N1 AND NL LEFT N1->N2. C C SET INDF AND INDX TO THE INDICES OF THE FIRST AND LAST C NEIGHBORS OF N1 AND SET N1LST TO THE LAST NEIGHBOR. C INDF = 1 IF (N1 .GT. 1) INDF = IEND(N1-1) + 1 INDX = IEND(N1) N1LST = IADJ(INDX) IF (N1LST .EQ. 0) INDX = INDX - 1 IF (N1LST .EQ. 0) GO TO 2 C C N1 IS AN INTERIOR NODE. LOOP THROUGH THE NEIGHBORS NL C IN REVERSE ORDER UNTIL NL LEFT N1->N2. C NL = N1LST 1 IF ( LEFT(X1,Y1,X2,Y2,X(NL),Y(NL)) ) GO TO 2 INDX = INDX - 1 NL = IADJ(INDX) IF (INDX .GT. INDF) GO TO 1 C C NL IS THE FIRST NEIGHBOR OF N1. SET NR TO THE LAST C NEIGHBOR AND TEST FOR AN ARC N1-N2. C NR = N1LST IF (NL .EQ. N2) GO TO 34 GO TO 4 C C NL = IADJ(INDX) LEFT N1->N2 AND INDX .GT. INDF. SET C NR TO THE PRECEDING NEIGHBOR OF N1. C 2 INDX = INDX - 1 NR = IADJ(INDX) IF ( LEFT(X2,Y2,X1,Y1,X(NR),Y(NR)) ) GO TO 3 IF (INDX .GT. INDF) GO TO 2 C C SET NL AND NR TO THE FIRST AND LAST NEIGHBORS OF N1 AND C TEST FOR AN INVALID DATA STRUCTURE (N1 CANNOT BE A C BOUNDARY NODE AND CANNOT BE ADJACENT TO N2). C NL = NR NR = N1LST IF (NR .EQ. 0 .OR. NR .EQ. N2) GO TO 37 GO TO 4 C C SET NL TO THE NEIGHBOR FOLLOWING NR AND TEST FOR AN ARC C N1-N2. C 3 NL = IADJ(INDX+1) IF (NL .EQ. N2 .OR. NR .EQ. N2) GO TO 34 C C STORE THE ORDERED SEQUENCE OF INTERSECTING EDGES NL->NR IN C IWK(1,IWL)->IWK(2,IWL). C 4 IWL = IWL + 1 IF (IWL .LE. IWEND) IWK(1,IWL) = NL IF (IWL .LE. IWEND) IWK(2,IWL) = NR C C SET NEXT TO THE NEIGHBOR OF NL WHICH FOLLOWS NR. C INDX = IEND(NL) IF (IADJ(INDX) .NE. NR) GO TO 5 C C NR IS THE LAST NEIGHBOR OF NL. SET NEXT TO THE FIRST C NEIGHBOR. C INDX = 0 IF (NL .NE. 1) INDX = IEND(NL-1) GO TO 6 C C NR IS NOT THE LAST NEIGHBOR OF NL. LOOP THROUGH THE C NEIGHBORS IN REVERSE ORDER. C 5 INDX = INDX - 1 IF (IADJ(INDX) .NE. NR) GO TO 5 C C STORE NEXT, TEST FOR AN INVALID TRIANGULATION (NL->NR C CANNOT BE A BOUNDARY EDGE), AND TEST FOR TERMINATION C OF THE LOOP. C 6 NEXT = IADJ(INDX+1) IF (NEXT .EQ. 0) GO TO 37 IF (NEXT .EQ. N2) GO TO 8 C C SET NL OR NR TO NEXT. C IF ( LEFT(X1,Y1,X2,Y2,X(NEXT),Y(NEXT)) ) GO TO 7 NR = NEXT GO TO 4 7 NL = NEXT GO TO 4 C C IWL IS THE NUMBER OF ARCS WHICH INTERSECT N1-N2. STORE C LWK AND TEST FOR SUFFICIENT SPACE. C 8 LWK = IWL IF (IWL .GT. IWEND) GO TO 36 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 IER = 0 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 9 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 10 IF (IWC .EQ. IWL) GO TO 21 IWCP1 = IWC + 1 NEXT = IWK(1,IWCP1) IF (NEXT .NE. NL) GO TO 15 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 13 IF (LFT .GE. 0) GO TO 11 IF ( .NOT. LEFT(X(NL),Y(NL),X0,Y0,X(NEXT),Y(NEXT)) ) . GO TO 13 C C REPLACE NL->NR WITH N0->NEXT. C CALL SWAP(NEXT,N0,NL,NR, IADJ,IEND ) IWK(1,IWC) = N0 IWK(2,IWC) = NEXT GO TO 14 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 11 CALL SWAP(NEXT,N0,NL,NR, IADJ,IEND ) DO 12 I = IWCP1,IWL IWK(1,I-1) = IWK(1,I) 12 IWK(2,I-1) = IWK(2,I) IWK(1,IWL) = N0 IWK(2,IWL) = NEXT IWL = IWL - 1 NR = NEXT GO TO 10 C C A SWAP IS NOT POSSIBLE. SET N0 TO NR. C 13 N0 = NR X0 = X(N0) Y0 = Y(N0) LFT = 1 C C ADVANCE TO THE NEXT ARC. C 14 NR = NEXT IWC = IWC + 1 GO TO 10 C C NEXT LEFT N1->N2, NEXT .NE. N2, AND IWC .LT. IWL. C TEST FOR A POSSIBLE SWAP. C 15 IF ( .NOT. LEFT(X(NL),Y(NL),X0,Y0,X(NEXT),Y(NEXT)) ) . GO TO 19 IF (LFT .LE. 0) GO TO 16 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, IADJ,IEND ) 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 16 CALL SWAP(NEXT,N0,NL,NR, IADJ,IEND ) I = IWC 17 IF (I .EQ. IWF) GO TO 18 IWK(1,I) = IWK(1,I-1) IWK(2,I) = IWK(2,I-1) I = I - 1 GO TO 17 18 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 10 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 9 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, IADJ,IEND ) IWK(1,IWL) = N0 IWK(2,IWL) = N2 IWL = IWL - 1 GO TO 9 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 9 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, IADJ,IEND ) 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 9 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, IADJ,IEND ) IWK(1,IWC) = 0 IWK(2,IWC) = 0 IF (IWC .EQ. 1) GO TO 29 C C OPTIMIZATION LOOPS -- OPTIMIZE THE SET OF NEW ARCS TO THE C LEFT OF IN1->IN2. THE LOOP IS REPEATED UNTIL NO SWAPS C ARE PERFORMED. C IWCM1 = IWC - 1 25 SWP = .FALSE. DO 28 I = 1,IWCM1 IO1 = IWK(1,I) IO2 = IWK(2,I) C C SET N1 TO THE NEIGHBOR OF IO1 WHICH FOLLOWS IO2 AND SET C N2 TO THE NEIGHBOR OF IO1 WHICH PRECEDES IO2. C INDF = 1 IF (IO1 .GT. 1) INDF = IEND(IO1-1) + 1 INDL = IEND(IO1) INDX = INDL IF (IADJ(INDX) .NE. IO2) GO TO 26 C C IO2 IS THE LAST NEIGHBOR OF IO1. C N1 = IADJ(INDF) N2 = IADJ(INDX-1) GO TO 27 C C IO2 IS NOT THE LAST NEIGHBOR OF IO1. LOOP THROUGH THE C NEIGHBORS IN REVERSE ORDER. C 26 INDX = INDX - 1 IF (IADJ(INDX) .NE. IO2) GO TO 26 N1 = IADJ(INDX+1) IF (INDX .NE. INDF) N2 = IADJ(INDX-1) IF (INDX .EQ. INDF) N2 = IADJ(INDL) C C TEST IO1-IO2 FOR A SWAP. C 27 IF ( .NOT. SWPTST(N1,N2,IO1,IO2,X,Y) ) GO TO 28 SWP = .TRUE. CALL SWAP(N1,N2,IO1,IO2, IADJ,IEND ) IWK(1,I) = N1 IWK(2,I) = N2 28 CONTINUE IF (SWP) GO TO 25 C C TEST FOR TERMINATION. C 29 IF (IWC .EQ. IWEND) RETURN IWCP1 = IWC + 1 C C OPTIMIZE THE SET OF NEW ARCS TO THE RIGHT OF IN1->IN2. C 30 SWP = .FALSE. DO 33 I = IWCP1,IWEND 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. C INDF = 1 IF (IO1 .GT. 1) INDF = IEND(IO1-1) + 1 INDL = IEND(IO1) INDX = INDL IF (IADJ(INDX) .NE. IO2) GO TO 31 C N1 = IADJ(INDF) N2 = IADJ(INDX-1) GO TO 32 C 31 INDX = INDX - 1 IF (IADJ(INDX) .NE. IO2) GO TO 31 N1 = IADJ(INDX+1) IF (INDX .NE. INDF) N2 = IADJ(INDX-1) IF (INDX .EQ. INDF) N2 = IADJ(INDL) C 32 IF ( .NOT. SWPTST(N1,N2,IO1,IO2,X,Y) ) GO TO 33 SWP = .TRUE. CALL SWAP(N1,N2,IO1,IO2, IADJ,IEND ) IWK(1,I) = N1 IWK(2,I) = N2 33 CONTINUE IF (SWP) GO TO 30 RETURN C C IN1 AND IN2 WERE ADJACENT ON INPUT. C 34 IER = 0 LWK = 0 RETURN C C PARAMETER OUT OF RANGE C 35 IER = 1 RETURN C C INSUFFICIENT SPACE IN IWK C 36 IER = 2 RETURN C C INVALID TRIANGULATION DATA STRUCTURE C 37 IER = 3 RETURN END SUBROUTINE GETNP (X,Y,IADJ,IEND,L, NPTS, DS,IER) INTEGER IADJ(1), IEND(1), L, NPTS(L), IER REAL X(1), Y(1), DS C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C GIVEN A THIESSEN TRIANGULATION OF N NODES AND AN ARRAY C NPTS CONTAINING THE INDICES OF L-1 NODES ORDERED BY C EUCLIDEAN DISTANCE FROM NPTS(1), THIS SUBROUTINE SETS C NPTS(L) TO THE INDEX OF THE NEXT NODE IN THE SEQUENCE -- C THE NODE, OTHER THAN NPTS(1),...,NPTS(L-1), WHICH IS C CLOSEST TO NPTS(1). THUS, THE ORDERED SEQUENCE OF K C CLOSEST NODES TO N1 (INCLUDING N1) MAY BE DETERMINED BY C K-1 CALLS TO GETNP WITH NPTS(1) = N1 AND L = 2,3,...,K C FOR K .GE. 2. C THE ALGORITHM USES THE FACT THAT, IN A THIESSEN TRIAN- C GULATION, THE K-TH CLOSEST NODE TO A GIVEN NODE N1 IS A C NEIGHBOR OF ONE OF THE K-1 CLOSEST NODES TO N1. C C INPUT PARAMETERS - X,Y - VECTORS OF LENGTH N CONTAINING C THE CARTESIAN COORDINATES OF THE C NODES. C C IADJ - SET OF ADJACENCY LISTS OF NODES C IN THE TRIANGULATION. C C IEND - POINTERS TO THE ENDS OF ADJACENCY C LISTS FOR EACH NODE IN THE TRI- C ANGULATION. C C L - NUMBER OF NODES IN THE SEQUENCE C ON OUTPUT. 2 .LE. L .LE. N. C C NPTS - ARRAY OF LENGTH .GE. L CONTAIN- C ING THE INDICES OF THE L-1 CLOS- C EST NODES TO NPTS(1) IN THE FIRST C L-1 LOCATIONS. C C IADJ AND IEND MAY BE CREATED BY SUBROUTINE TRMESH. C C INPUT PARAMETERS OTHER THAN NPTS ARE NOT ALTERED BY THIS C ROUTINE. C C OUTPUT PARAMETERS - NPTS - UPDATED WITH THE INDEX OF THE C L-TH CLOSEST NODE TO NPTS(1) IN C POSITION L UNLESS IER = 1. C C DS - SQUARED EUCLIDEAN DISTANCE BE- C TWEEN NPTS(1) AND NPTS(L) C UNLESS IER = 1. C C IER - ERROR INDICATOR C IER = 0 IF NO ERRORS WERE EN- C COUNTERED. C IER = 1 IF L IS OUT OF RANGE. C C MODULES REFERENCED BY GETNP - NONE C C INTRINSIC FUNCTION CALLED BY GETNP - IABS C C*********************************************************** C INTEGER LM1, N1, I, NI, NP, INDF, INDL, INDX, NB REAL X1, Y1, DNP, DNB C C LOCAL PARAMETERS - C C LM1 = L - 1 C N1 = NPTS(1) C I = NPTS INDEX AND DO-LOOP INDEX C NI = NPTS(I) C NP = CANDIDATE FOR NPTS(L) C INDF = IADJ INDEX OF THE FIRST NEIGHBOR OF NI C INDL = IADJ INDEX OF THE LAST NEIGHBOR OF NI C INDX = IADJ INDEX IN THE RANGE INDF,...,INDL C NB = NEIGHBOR OF NI AND CANDIDATE FOR NP C X1,Y1 = COORDINATES OF N1 C DNP,DNB = SQUARED DISTANCES FROM N1 TO NP AND NB, C RESPECTIVELY C LM1 = L - 1 IF (LM1 .LT. 1) GO TO 4 IER = 0 N1 = NPTS(1) X1 = X(N1) Y1 = Y(N1) C C MARK THE ELEMENTS OF NPTS C DO 1 I = 1,LM1 NI = NPTS(I) IEND(NI) = -IEND(NI) 1 CONTINUE C C CANDIDATES FOR NP = NPTS(L) ARE THE UNMARKED NEIGHBORS C OF NODES IN NPTS. NP=0 IS A FLAG TO SET NP TO THE C FIRST CANDIDATE ENCOUNTERED. C NP = 0 DNP = 0. C C LOOP ON NODES NI IN NPTS C DO 2 I = 1,LM1 NI = NPTS(I) INDF = 1 IF (NI .GT. 1) INDF = IABS(IEND(NI-1)) + 1 INDL = -IEND(NI) C C LOOP ON NEIGHBORS NB OF NI C DO 2 INDX = INDF,INDL NB = IADJ(INDX) IF (NB .EQ. 0 .OR. IEND(NB) .LT. 0) GO TO 2 C C NB IS AN UNMARKED NEIGHBOR OF NI. REPLACE NP IF NB IS C CLOSER TO N1 OR IS THE FIRST CANDIDATE ENCOUNTERED. C DNB = (X(NB)-X1)**2 + (Y(NB)-Y1)**2 IF (NP .NE. 0 .AND. DNB .GE. DNP) GO TO 2 NP = NB DNP = DNB 2 CONTINUE NPTS(L) = NP DS = DNP C C UNMARK THE ELEMENTS OF NPTS C DO 3 I = 1,LM1 NI = NPTS(I) IEND(NI) = -IEND(NI) 3 CONTINUE RETURN C C L IS OUT OF RANGE C 4 IER = 1 RETURN END INTEGER FUNCTION INDEX (NVERTX,NABOR,IADJ,IEND) INTEGER NVERTX, NABOR, IADJ(1), IEND(1) C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C THIS FUNCTION RETURNS THE INDEX OF NABOR IN THE C ADJACENCY LIST FOR NVERTX. C C INPUT PARAMETERS - NVERTX - NODE WHOSE ADJACENCY LIST IS C TO BE SEARCHED. C C NABOR - NODE WHOSE INDEX IS TO BE C RETURNED. NABOR MUST BE C CONNECTED TO NVERTX. C C IADJ - SET OF ADJACENCY LISTS. C C IEND - POINTERS TO THE ENDS OF C ADJACENCY LISTS IN IADJ. C C INPUT PARAMETERS ARE NOT ALTERED BY THIS FUNCTION. C C OUTPUT PARAMETER - INDEX - IADJ(INDEX) = NABOR. C C MODULES REFERENCED BY INDEX - NONE C C*********************************************************** C INTEGER NB, INDX C C LOCAL PARAMETERS - C C NB = LOCAL COPY OF NABOR C INDX = INDEX FOR IADJ C NB = NABOR C C INITIALIZATION C INDX = IEND(NVERTX) + 1 C C SEARCH THE LIST OF NVERTX NEIGHBORS FOR NB C 1 INDX = INDX - 1 IF (IADJ(INDX) .NE. NB) GO TO 1 C INDEX = INDX RETURN END SUBROUTINE INTADD (KK,I1,I2,I3, IADJ,IEND ) INTEGER KK, I1, I2, I3, IADJ(1), IEND(KK) C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C THIS ROUTINE ADDS AN INTERIOR NODE TO A TRIANGULATION C OF A SET OF KK-1 POINTS IN THE PLANE. IADJ AND IEND ARE C UPDATED WITH THE INSERTION OF NODE KK IN THE TRIANGLE C WHOSE VERTICES ARE I1, I2, AND I3. C C INPUT PARAMETERS - KK - INDEX OF NODE TO BE C INSERTED. KK .GE. 4. C C I1,I2,I3 - INDICES OF THE VERTICES OF C A TRIANGLE CONTAINING NODE C KK -- IN COUNTERCLOCKWISE C ORDER. C C IADJ - SET OF ADJACENCY LISTS C OF NODES IN THE MESH. C C IEND - POINTERS TO THE ENDS OF C ADJACENCY LISTS IN IADJ FOR C EACH NODE IN THE MESH. C C IADJ AND IEND MAY BE CREATED BY TRMESH AND MUST CONTAIN C THE VERTICES I1, I2, AND I3. I1,I2,I3 MAY BE DETERMINED C BY TRFIND. C C KK, I1, I2, AND I3 ARE NOT ALTERED BY THIS ROUTINE. C C OUTPUT PARAMETERS - IADJ,IEND - UPDATED WITH THE ADDITION C OF NODE KK AS THE LAST C ENTRY. NODE KK WILL BE C CONNECTED TO NODES I1, I2, C AND I3. NO OPTIMIZATION C OF THE MESH IS PERFORMED. C C MODULE REFERENCED BY INTADD - SHIFTD C C INTRINSIC FUNCTION CALLED BY INTADD - MOD C C*********************************************************** C INTEGER K, KM1, N(3), NFT(3), IP1, IP2, IP3, INDX, NF, . NL, N1, N2, IMIN, IMAX, I, ITEMP C C LOCAL PARAMETERS - C C K = LOCAL COPY OF KK C KM1 = K - 1 C N = VECTOR CONTAINING I1, I2, I3 C NFT = POINTERS TO THE TOPS OF THE 3 SETS OF IADJ C ELEMENTS TO BE SHIFTED DOWNWARD C IP1,IP2,IP3 = PERMUTATION INDICES FOR N AND NFT C INDX = INDEX FOR IADJ AND N C NF,NL = INDICES OF FIRST AND LAST ENTRIES IN IADJ C TO BE SHIFTED DOWN C N1,N2 = FIRST 2 VERTICES OF A NEW TRIANGLE -- C (N1,N2,KK) C IMIN,IMAX = BOUNDS ON DO-LOOP INDEX -- FIRST AND LAST C ELEMENTS OF IEND TO BE INCREMENTED C I = DO-LOOP INDEX C ITEMP = TEMPORARY STORAGE LOCATION C K = KK C C INITIALIZATION C N(1) = I1 N(2) = I2 N(3) = I3 C C SET UP NFT C DO 2 I = 1,3 N1 = N(I) INDX = MOD(I,3) + 1 N2 = N(INDX) INDX = IEND(N1) + 1 C C FIND THE INDEX OF N2 AS A NEIGHBOR OF N1 C 1 INDX = INDX - 1 IF (IADJ(INDX) .NE. N2) GO TO 1 NFT(I) = INDX + 1 2 CONTINUE C C ORDER THE VERTICES BY DECREASING MAGNITUDE. C N(IP(I+1)) PRECEDES N(IP(I)) IN IEND FOR C I = 1,2. C IP1 = 1 IP2 = 2 IP3 = 3 IF ( N(2) .LE. N(1) ) GO TO 3 IP1 = 2 IP2 = 1 3 IF ( N(3) .LE. N(IP1) ) GO TO 4 IP3 = IP1 IP1 = 3 4 IF ( N(IP3) .LE. N(IP2) ) GO TO 5 ITEMP = IP2 IP2 = IP3 IP3 = ITEMP C C ADD NODE K TO THE ADJACENCY LISTS OF EACH VERTEX AND C UPDATE IEND. FOR EACH VERTEX, A SET OF IADJ ELEMENTS C IS SHIFTED DOWNWARD AND K IS INSERTED. SHIFTING STARTS C AT THE END OF THE ARRAY. C 5 KM1 = K - 1 NL = IEND(KM1) NF = NFT(IP1) IF (NF .LE. NL) CALL SHIFTD(NF,NL,3, IADJ ) IADJ(NF+2) = K IMIN = N(IP1) IMAX = KM1 DO 6 I = IMIN,IMAX IEND(I) = IEND(I) + 3 6 CONTINUE C NL = NF - 1 NF = NFT(IP2) CALL SHIFTD(NF,NL,2, IADJ ) IADJ(NF+1) = K IMAX = IMIN - 1 IMIN = N(IP2) DO 7 I = IMIN,IMAX IEND(I) = IEND(I) + 2 7 CONTINUE C NL = NF - 1 NF = NFT(IP3) CALL SHIFTD(NF,NL,1, IADJ ) IADJ(NF) = K IMAX = IMIN - 1 IMIN = N(IP3) DO 8 I = IMIN,IMAX IEND(I) = IEND(I) + 1 8 CONTINUE C C ADD NODE K TO IEND AND ITS NEIGHBORS TO IADJ C INDX = IEND(KM1) IEND(K) = INDX + 3 DO 9 I = 1,3 INDX = INDX + 1 IADJ(INDX) = N(I) 9 CONTINUE RETURN END SUBROUTINE PERMUT (NN,IP, A ) INTEGER NN, IP(NN) REAL A(NN) C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C THIS ROUTINE APPLIES A SET OF PERMUTATIONS TO A VECTOR. C C INPUT PARAMETERS - NN - LENGTH OF A AND IP. C C IP - VECTOR CONTAINING THE SEQUENCE OF C INTEGERS 1,...,NN PERMUTED IN THE C SAME FASHION THAT A IS TO BE PER- C MUTED. C C A - VECTOR TO BE PERMUTED. C C NN AND IP ARE NOT ALTERED BY THIS ROUTINE. C C OUTPUT PARAMETERS - A - REORDERED VECTOR REFLECTING THE C PERMUTATIONS DEFINED BY IP. C C MODULES REFERENCED BY PERMUT - NONE C C*********************************************************** C INTEGER N, K, J, IPJ REAL TEMP C C LOCAL PARAMETERS - C C N = LOCAL COPY OF NN C K = INDEX FOR IP AND FOR THE FIRST ELEMENT OF A IN A C PERMUTATION C J = INDEX FOR IP AND A, J .GE. K C IPJ = IP(J) C TEMP = TEMPORARY STORAGE FOR A(K) C N = NN IF (N .LT. 2) RETURN K = 1 C C LOOP ON PERMUTATIONS C 1 J = K TEMP = A(K) C C APPLY PERMUTATION TO A. IP(J) IS MARKED (MADE NEGATIVE) C AS BEING INCLUDED IN THE PERMUTATION. C 2 IPJ = IP(J) IP(J) = -IPJ IF (IPJ .EQ. K) GO TO 3 A(J) = A(IPJ) J = IPJ GO TO 2 3 A(J) = TEMP C C SEARCH FOR AN UNMARKED ELEMENT OF IP C 4 K = K + 1 IF (K .GT. N) GO TO 5 IF (IP(K) .GT. 0) GO TO 1 GO TO 4 C C ALL PERMUTATIONS HAVE BEEN APPLIED. UNMARK IP. C 5 DO 6 K = 1,N IP(K) = -IP(K) 6 CONTINUE RETURN END SUBROUTINE QSORT (N,X, IND) INTEGER N, IND(N) REAL X(N) C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C THIS SUBROUTINE USES AN ORDER N*LOG(N) QUICK SORT TO C SORT THE REAL ARRAY X INTO INCREASING ORDER. THE ALGOR- C ITHM IS AS FOLLOWS. IND IS INITIALIZED TO THE ORDERED C SEQUENCE OF INDICES 1,...,N, AND ALL INTERCHANGES ARE C APPLIED TO IND. X IS DIVIDED INTO TWO PORTIONS BY PICKING C A CENTRAL ELEMENT T. THE FIRST AND LAST ELEMENTS ARE COM- C PARED WITH T, AND INTERCHANGES ARE APPLIED AS NECESSARY SO C THAT THE THREE VALUES ARE IN ASCENDING ORDER. INTER- C CHANGES ARE THEN APPLIED SO THAT ALL ELEMENTS GREATER THAN C T ARE IN THE UPPER PORTION OF THE ARRAY AND ALL ELEMENTS C LESS THAN T ARE IN THE LOWER PORTION. THE UPPER AND LOWER C INDICES OF ONE OF THE PORTIONS ARE SAVED IN LOCAL ARRAYS, C AND THE PROCESS IS REPEATED ITERATIVELY ON THE OTHER C PORTION. WHEN A PORTION IS COMPLETELY SORTED, THE PROCESS C BEGINS AGAIN BY RETRIEVING THE INDICES BOUNDING ANOTHER C UNSORTED PORTION. C C INPUT PARAMETERS - N - LENGTH OF THE ARRAY X. C C X - VECTOR OF LENGTH N TO BE SORTED. C C IND - VECTOR OF LENGTH .GE. N. C C N AND X ARE NOT ALTERED BY THIS ROUTINE. C C OUTPUT PARAMETER - IND - SEQUENCE OF INDICES 1,...,N C PERMUTED IN THE SAME FASHION AS X C WOULD BE. THUS, THE ORDERING ON C X IS DEFINED BY Y(I) = X(IND(I)). C C MODULES REFERENCED BY QSORT - NONE C C INTRINSIC FUNCTIONS CALLED BY QSORT - IFIX, FLOAT C C*********************************************************** C C NOTE -- IU AND IL MUST BE DIMENSIONED .GE. LOG(N) WHERE C LOG HAS BASE 2. C C*********************************************************** C INTEGER IU(21), IL(21) INTEGER M, I, J, K, L, IJ, IT, ITT, INDX REAL R, T C C LOCAL PARAMETERS - C C IU,IL = TEMPORARY STORAGE FOR THE UPPER AND LOWER C INDICES OF PORTIONS OF THE ARRAY X C M = INDEX FOR IU AND IL C I,J = LOWER AND UPPER INDICES OF A PORTION OF X C K,L = INDICES IN THE RANGE I,...,J C IJ = RANDOMLY CHOSEN INDEX BETWEEN I AND J C IT,ITT = TEMPORARY STORAGE FOR INTERCHANGES IN IND C INDX = TEMPORARY INDEX FOR X C R = PSEUDO RANDOM NUMBER FOR GENERATING IJ C T = CENTRAL ELEMENT OF X C IF (N .LE. 0) RETURN C C INITIALIZE IND, M, I, J, AND R C DO 1 I = 1,N IND(I) = I 1 CONTINUE M = 1 I = 1 J = N R = .375 C C TOP OF LOOP C 2 IF (I .GE. J) GO TO 10 IF (R .GT. .5898437) GO TO 3 R = R + .0390625 GO TO 4 3 R = R - .21875 C C INITIALIZE K C 4 K = I C C SELECT A CENTRAL ELEMENT OF X AND SAVE IT IN T C IJ = I + IFIX(R*FLOAT(J-I)) IT = IND(IJ) T = X(IT) C C IF THE FIRST ELEMENT OF THE ARRAY IS GREATER THAN T, C INTERCHANGE IT WITH T C INDX = IND(I) IF (X(INDX) .LE. T) GO TO 5 IND(IJ) = INDX IND(I) = IT IT = INDX T = X(IT) C C INITIALIZE L C 5 L = J C C IF THE LAST ELEMENT OF THE ARRAY IS LESS THAN T, C INTERCHANGE IT WITH T C INDX = IND(J) IF (X(INDX) .GE. T) GO TO 7 IND(IJ) = INDX IND(J) = IT IT = INDX T = X(IT) C C IF THE FIRST ELEMENT OF THE ARRAY IS GREATER THAN T, C INTERCHANGE IT WITH T C INDX = IND(I) IF (X(INDX) .LE. T) GO TO 7 IND(IJ) = INDX IND(I) = IT IT = INDX T = X(IT) GO TO 7 C C INTERCHANGE ELEMENTS K AND L C 6 ITT = IND(L) IND(L) = IND(K) IND(K) = ITT C C FIND AN ELEMENT IN THE UPPER PART OF THE ARRAY WHICH IS C NOT LARGER THAN T C 7 L = L - 1 INDX = IND(L) IF (X(INDX) .GT. T) GO TO 7 C C FIND AN ELEMENT IN THE LOWER PART OF THE ARRAY WHCIH IS C NOT SMALLER THAN T C 8 K = K + 1 INDX = IND(K) IF (X(INDX) .LT. T) GO TO 8 C C IF K .LE. L, INTERCHANGE ELEMENTS K AND L C IF (K .LE. L) GO TO 6 C C SAVE THE UPPER AND LOWER SUBSCRIPTS OF THE PORTION OF THE C ARRAY YET TO BE SORTED C IF (L-I .LE. J-K) GO TO 9 IL(M) = I IU(M) = L I = K M = M + 1 GO TO 11 C 9 IL(M) = K IU(M) = J J = L M = M + 1 GO TO 11 C C BEGIN AGAIN ON ANOTHER UNSORTED PORTION OF THE ARRAY C 10 M = M - 1 IF (M .EQ. 0) RETURN I = IL(M) J = IU(M) C 11 IF (J-I .GE. 11) GO TO 4 IF (I .EQ. 1) GO TO 2 I = I - 1 C C SORT ELEMENTS I+1,...,J. NOTE THAT 1 .LE. I .LT. J AND C J-I .LT. 11. C 12 I = I + 1 IF (I .EQ. J) GO TO 10 INDX = IND(I+1) T = X(INDX) IT = INDX INDX = IND(I) IF (X(INDX) .LE. T) GO TO 12 K = I C 13 IND(K+1) = IND(K) K = K - 1 INDX = IND(K) IF (T .LT. X(INDX)) GO TO 13 IND(K+1) = IT GO TO 12 END SUBROUTINE REORDR (N,IFLAG, A,B,C, IND) INTEGER N, IFLAG, IND(N) REAL A(N), B(N), C(N) C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C THIS SUBROUTINE USES AN ORDER N*LOG(N) QUICK SORT TO C REORDER THE REAL ARRAY A INTO INCREASING ORDER. A RECORD C OF THE PERMUTATIONS APPLIED TO A IS STORED IN IND, AND C THESE PERMUTATIONS MAY BE APPLIED TO ONE OR TWO ADDITIONAL C VECTORS BY THIS ROUTINE. ANY OTHER VECTOR V MAY BE PER- C MUTED IN THE SAME FASHION BY CALLING SUBROUTINE PERMUT C WITH N, IND, AND V AS PARAMETERS. C A SET OF NODES (X(I),Y(I)) AND DATA VALUES Z(I) MAY BE C PREPROCESSED BY REORDR FOR INCREASED EFFICIENCY IN THE C TRIANGULATION ROUTINE TRMESH. EFFICIENCY IS INCREASED BY C A FACTOR OF APPROXIMATELY SQRT(N)/6 FOR RANDOMLY DISTRIB- C UTED NODES, AND THE PREPROCESSING IS ALSO USEFUL FOR C DETECTING DUPLICATE NODES. EITHER X OR Y MAY BE USED AS C THE SORT KEY (ASSOCIATED WITH A). C C INPUT PARAMETERS - N - NUMBER OF NODES. C C IFLAG - NUMBER OF VECTORS TO BE PERMUTED. C IFLAG .LE. 0 IF A, B, AND C ARE TO C REMAIN UNALTERED. C IFLAG .EQ. 1 IF ONLY A IS TO BE C PERMUTED. C IFLAG .EQ. 2 IF A AND B ARE TO BE C PERMUTED. C IFLAG .GE. 3 IF A, B, AND C ARE TO C BE PERMUTED. C C A,B,C - VECTORS OF LENGTH N TO BE SORTED C (ON THE COMPONENTS OF A), OR DUMMY C PARAMETERS, DEPENDING ON IFLAG. C C IND - VECTOR OF LENGTH .GE. N. C C N, IFLAG, AND ANY DUMMY PARAMETERS ARE NOT ALTERED BY THIS C ROUTINE. C C OUTPUT PARAMETERS - A,B,C - SORTED OR UNALTERED VECTORS. C C IND - SEQUENCE OF INDICES 1,...,N C PERMUTED IN THE SAME FASHION C AS THE REAL VECTORS. THUS, C THE ORDERING MAY BE APPLIED TO C A REAL VECTOR V AND STORED IN C W BY SETTING W(I) = V(IND(I)), C OR V MAY BE OVERWRITTEN WITH C THE ORDERING BY A CALL TO PER- C MUT. C C MODULES REFERENCED BY REORDR - QSORT, PERMUT C C*********************************************************** C INTEGER NN, NV C C LOCAL PARAMETERS - C C NN = LOCAL COPY OF N C NV = LOCAL COPY OF IFLAG C NN = N NV = IFLAG CALL QSORT(NN,A, IND) IF (NV .LE. 0) RETURN CALL PERMUT(NN,IND, A ) IF (NV .EQ. 1) RETURN CALL PERMUT(NN,IND, B ) IF (NV .EQ. 2) RETURN CALL PERMUT(NN,IND, C ) RETURN END SUBROUTINE SHIFTD (NFRST,NLAST,KK, IARR ) INTEGER NFRST, NLAST, KK, IARR(1) C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C THIS ROUTINE SHIFTS A SET OF CONTIGUOUS ELEMENTS OF AN C INTEGER ARRAY KK POSITIONS DOWNWARD (UPWARD IF KK .LT. 0). C THE LOOPS ARE UNROLLED IN ORDER TO INCREASE EFFICIENCY. C C INPUT PARAMETERS - NFRST,NLAST - BOUNDS ON THE PORTION OF C IARR TO BE SHIFTED. ALL C ELEMENTS BETWEEN AND C INCLUDING THE BOUNDS ARE C SHIFTED UNLESS NFRST .GT. C NLAST, IN WHICH CASE NO C SHIFT OCCURS. C C KK - NUMBER OF POSITIONS EACH C ELEMENT IS TO BE SHIFTED. C IF KK .LT. 0 SHIFT UP. C IF KK .GT. 0 SHIFT DOWN. C C IARR - INTEGER ARRAY OF LENGTH C .GE. NLAST + MAX(KK,0). C C NFRST, NLAST, AND KK ARE NOT ALTERED BY THIS ROUTINE. C C OUTPUT PARAMETER - IARR - SHIFTED ARRAY. C C MODULES REFERENCED BY SHIFTD - NONE C C*********************************************************** C INTEGER INC, K, NF, NL, NLP1, NS, NSL, I, IBAK, INDX, . IMAX DATA INC/5/ C C LOCAL PARAMETERS - C C INC = DO-LOOP INCREMENT (UNROLLING FACTOR) -- IF INC IS C CHANGED, STATEMENTS MUST BE ADDED TO OR DELETED C FROM THE DO-LOOPS C K = LOCAL COPY OF KK C NF = LOCAL COPY OF NFRST C NL = LOCAL COPY OF NLAST C NLP1 = NL + 1 C NS = NUMBER OF SHIFTS C NSL = NUMBER OF SHIFTS DONE IN UNROLLED DO-LOOP (MULTIPLE C OF INC) C I = DO-LOOP INDEX AND INDEX FOR IARR C IBAK = INDEX FOR DOWNWARD SHIFT OF IARR C INDX = INDEX FOR IARR C IMAX = BOUND ON DO-LOOP INDEX C K = KK NF = NFRST NL = NLAST IF (NF .GT. NL .OR. K .EQ. 0) RETURN NLP1 = NL + 1 NS = NLP1 - NF NSL = INC*(NS/INC) IF ( K .LT. 0) GO TO 4 C C SHIFT DOWNWARD STARTING FROM THE BOTTOM C IF (NSL .LE. 0) GO TO 2 DO 1 I = 1,NSL,INC IBAK = NLP1 - I INDX = IBAK + K IARR(INDX) = IARR(IBAK) IARR(INDX-1) = IARR(IBAK-1) IARR(INDX-2) = IARR(IBAK-2) IARR(INDX-3) = IARR(IBAK-3) IARR(INDX-4) = IARR(IBAK-4) 1 CONTINUE C C PERFORM THE REMAINING NS-NSL SHIFTS ONE AT A TIME C 2 IBAK = NLP1 - NSL 3 IF (IBAK .LE. NF) RETURN IBAK = IBAK - 1 INDX = IBAK + K IARR(INDX) = IARR(IBAK) GO TO 3 C C SHIFT UPWARD STARTING FROM THE TOP C 4 IF (NSL .LE. 0) GO TO 6 IMAX = NLP1 - INC DO 5 I = NF,IMAX,INC INDX = I + K IARR(INDX) = IARR(I) IARR(INDX+1) = IARR(I+1) IARR(INDX+2) = IARR(I+2) IARR(INDX+3) = IARR(I+3) IARR(INDX+4) = IARR(I+4) 5 CONTINUE C C PERFORM THE REMAINING NS-NSL SHIFTS ONE AT A TIME C 6 I = NSL + NF 7 IF (I .GT. NL) RETURN INDX = I + K IARR(INDX) = IARR(I) I = I + 1 GO TO 7 END SUBROUTINE SWAP (NIN1,NIN2,NOUT1,NOUT2, IADJ,IEND ) INTEGER NIN1, NIN2, NOUT1, NOUT2, IADJ(1), IEND(1) EXTERNAL INDEX C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C THIS SUBROUTINE SWAPS THE DIAGONALS IN A CONVEX QUADRI- C LATERAL. C C INPUT PARAMETERS - NIN1,NIN2,NOUT1,NOUT2 - NODAL INDICES C OF A PAIR OF ADJACENT TRIANGLES C WHICH FORM A CONVEX QUADRILAT- C ERAL. NOUT1 AND NOUT2 ARE CON- C NECTED BY AN ARC WHICH IS TO BE C REPLACED BY THE ARC NIN1-NIN2. C (NIN1,NOUT1,NOUT2) MUST BE TRI- C ANGLE VERTICES IN COUNTERCLOCK- C WISE ORDER. C C THE ABOVE PARAMETERS ARE NOT ALTERED BY THIS ROUTINE. C C IADJ,IEND - TRIANGULATION DATA STRUCTURE C (SEE SUBROUTINE TRMESH). C C OUTPUT PARAMETERS - IADJ,IEND - UPDATED WITH THE ARC C REPLACEMENT. C C MODULES REFERENCED BY SWAP - INDEX, SHIFTD C C*********************************************************** C INTEGER IN(2), IO(2), IP1, IP2, J, K, NF, NL, I, . IMIN, IMAX C C LOCAL PARAMETERS - C C IN = NIN1 AND NIN2 ORDERED BY INCREASING MAGNITUDE C (THE NEIGHBORS OF IN(1) PRECEDE THOSE OF C IN(2) IN IADJ) C IO = NOUT1 AND NOUT2 IN INCREASING ORDER C IP1,IP2 = PERMUTATION OF (1,2) SUCH THAT IO(IP1) C PRECEDES IO(IP2) AS A NEIGHBOR OF IN(1) C J,K = PERMUTATION OF (1,2) USED AS INDICES OF IN C AND IO C NF,NL = IADJ INDICES BOUNDARY A PORTION OF THE ARRAY C TO BE SHIFTED C I = IEND INDEX C IMIN,IMAX = BOUNDS ON THE PORTION OF IEND TO BE INCRE- C MENTED OR DECREMENTED C IN(1) = NIN1 IN(2) = NIN2 IO(1) = NOUT1 IO(2) = NOUT2 IP1 = 1 C C ORDER THE INDICES SO THAT IN(1) .LT. IN(2) AND IO(1) .LT. C IO(2), AND CHOOSE IP1 AND IP2 SUCH THAT (IN(1),IO(IP1), C IO(IP2)) FORMS A TRIANGLE. C IF (IN(1) .LT. IN(2)) GO TO 1 IN(1) = IN(2) IN(2) = NIN1 IP1 = 2 1 IF (IO(1) .LT. IO(2)) GO TO 2 IO(1) = IO(2) IO(2) = NOUT1 IP1 = 3 - IP1 2 IP2 = 3 - IP1 IF (IO(2) .LT. IN(1)) GO TO 8 IF (IN(2) .LT. IO(1)) GO TO 12 C C IN(1) AND IO(1) PRECEDE IN(2) AND IO(2). FOR (J,K) = C (1,2) AND (2,1), DELETE IO(K) AS A NEIGHBOR OF IO(J) C BY SHIFTING A PORTION OF IADJ EITHER UP OR DOWN AND C AND INSERT IN(K) AS A NEIGHBOR OF IN(J). C DO 7 J = 1,2 K = 3 - J IF (IN(J) .GT. IO(J)) GO TO 4 C C THE NEIGHBORS OF IN(J) PRECEDE THOSE OF IO(J) -- SHIFT C DOWN BY 1 C NF = 1 + INDEX(IN(J),IO(IP1),IADJ,IEND) NL = -1 + INDEX(IO(J),IO(K),IADJ,IEND) IF (NF .LE. NL) CALL SHIFTD(NF,NL,1, IADJ ) IADJ(NF) = IN(K) IMIN = IN(J) IMAX = IO(J)-1 DO 3 I = IMIN,IMAX 3 IEND(I) = IEND(I) + 1 GO TO 6 C C THE NEIGHBORS OF IO(J) PRECEDE THOSE OF IN(J) -- SHIFT C UP BY 1 C 4 NF = 1 + INDEX(IO(J),IO(K),IADJ,IEND) NL = -1 + INDEX(IN(J),IO(IP2),IADJ,IEND) IF (NF .LE. NL) CALL SHIFTD(NF,NL,-1, IADJ ) IADJ(NL) = IN(K) IMIN = IO(J) IMAX = IN(J) - 1 DO 5 I = IMIN,IMAX 5 IEND(I) = IEND(I) - 1 C C REVERSE (IP1,IP2) FOR (J,K) = (2,1) C 6 IP1 = IP2 IP2 = 3 - IP1 7 CONTINUE RETURN C C THE VERTICES ARE ORDERED (IO(1),IO(2),IN(1),IN(2)). C DELETE IO(2) BY SHIFTING UP BY 1 C 8 NF = 1 + INDEX(IO(1),IO(2),IADJ,IEND) NL = -1 + INDEX(IO(2),IO(1),IADJ,IEND) IF (NF .LE. NL) CALL SHIFTD(NF,NL,-1, IADJ ) IMIN = IO(1) IMAX = IO(2)-1 DO 9 I = IMIN,IMAX 9 IEND(I) = IEND(I) - 1 C C DELETE IO(1) BY SHIFTING UP BY 2 AND INSERT IN(2) C NF = NL + 2 NL = -1 + INDEX(IN(1),IO(IP2),IADJ,IEND) IF (NF .LE. NL) CALL SHIFTD(NF,NL,-2, IADJ ) IADJ(NL-1) = IN(2) IMIN = IO(2) IMAX = IN(1)-1 DO 10 I = IMIN,IMAX 10 IEND(I) = IEND(I) - 2 C C SHIFT UP BY 1 AND INSERT IN(1) C NF = NL + 1 NL = -1 + INDEX(IN(2),IO(IP1),IADJ,IEND) CALL SHIFTD(NF,NL,-1, IADJ ) IADJ(NL) = IN(1) IMIN = IN(1) IMAX = IN(2)-1 DO 11 I = IMIN,IMAX 11 IEND(I) = IEND(I) - 1 RETURN C C THE VERTICES ARE ORDERED (IN(1),IN(2),IO(1),IO(2)). C DELETE IO(1) BY SHIFTING DOWN BY 1 C 12 NF = 1 + INDEX(IO(1),IO(2),IADJ,IEND) NL = -1 + INDEX(IO(2),IO(1),IADJ,IEND) IF (NF .LE. NL) CALL SHIFTD(NF,NL,1, IADJ ) IMIN = IO(1) IMAX = IO(2) - 1 DO 13 I = IMIN,IMAX 13 IEND(I) = IEND(I) + 1 C C DELETE IO(2) BY SHIFTING DOWN BY 2 AND INSERT IN(1) C NL = NF - 2 NF = 1 + INDEX(IN(2),IO(IP2),IADJ,IEND) IF (NF .LE. NL) CALL SHIFTD(NF,NL,2, IADJ ) IADJ(NF+1) = IN(1) IMIN = IN(2) IMAX = IO(1) - 1 DO 14 I = IMIN,IMAX 14 IEND(I) = IEND(I) + 2 C C SHIFT DOWN BY 1 AND INSERT IN(2) C NL = NF - 1 NF = 1 + INDEX(IN(1),IO(IP1),IADJ,IEND) CALL SHIFTD(NF,NL,1, IADJ ) IADJ(NF) = IN(2) IMIN = IN(1) IMAX = IN(2) - 1 DO 15 I = IMIN,IMAX 15 IEND(I) = IEND(I) + 1 RETURN END LOGICAL FUNCTION SWPTST (IN1,IN2,IO1,IO2,X,Y) INTEGER IN1, IN2, IO1, IO2 REAL X(1), Y(1) C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C THIS FUNCTION DECIDES WHETHER OR NOT TO REPLACE A C DIAGONAL ARC IN A QUADRILATERAL WITH THE OTHER DIAGONAL. C THE DETERMINATION IS BASED ON THE SIZES OF THE ANGLES C CONTAINED IN THE 2 TRIANGLES DEFINED BY THE DIAGONAL. C THE DIAGONAL IS CHOSEN TO MAXIMIZE THE SMALLEST OF THE C SIX ANGLES OVER THE TWO PAIRS OF TRIANGLES. C C INPUT PARAMETERS - IN1,IN2,IO1,IO2 - NODE INDICES OF THE C FOUR POINTS DEFINING THE C QUADRILATERAL. IO1 AND IO2 C ARE CURRENTLY CONNECTED BY A C DIAGONAL ARC. THIS ARC C SHOULD BE REPLACED BY AN ARC C CONNECTING IN1, IN2 IF THE C DECISION IS MADE TO SWAP. C IN1,IO1,IO2 MUST BE IN C COUNTERCLOCKWISE ORDER. C C X,Y - VECTORS OF NODAL COORDINATES. C (X(I),Y(I)) ARE THE COORD- C INATES OF NODE I FOR I = IN1, C IN2, IO1, OR IO2. C C NONE OF THE INPUT PARAMETERS ARE ALTERED BY THIS ROUTINE. C C OUTPUT PARAMETER - SWPTST - .TRUE. IFF THE ARC CONNECTING C IO1 AND IO2 IS TO BE REPLACED C C MODULES REFERENCED BY SWPTST - NONE C C*********************************************************** C REAL DX11, DX12, DX22, DX21, DY11, DY12, DY22, DY21, . SIN1, SIN2, COS1, COS2, SIN12 C C LOCAL PARAMETERS - C C DX11,DY11 = X,Y COORDINATES OF THE VECTOR IN1-IO1 C DX12,DY12 = X,Y COORDINATES OF THE VECTOR IN1-IO2 C DX22,DY22 = X,Y COORDINATES OF THE VECTOR IN2-IO2 C DX21,DY21 = X,Y COORDINATES OF THE VECTOR IN2-IO1 C SIN1 = CROSS PRODUCT OF THE VECTORS IN1-IO1 AND C IN1-IO2 -- PROPORTIONAL TO SIN(T1) WHERE T1 C IS THE ANGLE AT IN1 FORMED BY THE VECTORS C COS1 = INNER PRODUCT OF THE VECTORS IN1-IO1 AND C IN1-IO2 -- PROPORTIONAL TO COS(T1) C SIN2 = CROSS PRODUCT OF THE VECTORS IN2-IO2 AND C IN2-IO1 -- PROPORTIONAL TO SIN(T2) WHERE T2 C IS THE ANGLE AT IN2 FORMED BY THE VECTORS C COS2 = INNER PRODUCT OF THE VECTORS IN2-IO2 AND C IN2-IO1 -- PROPORTIONAL TO COS(T2) C SIN12 = SIN1*COS2 + COS1*SIN2 -- PROPORTIONAL TO C SIN(T1+T2) C SWPTST = .FALSE. C C COMPUTE THE VECTORS CONTAINING THE ANGLES T1, T2 C DX11 = X(IO1) - X(IN1) DX12 = X(IO2) - X(IN1) DX22 = X(IO2) - X(IN2) DX21 = X(IO1) - X(IN2) C DY11 = Y(IO1) - Y(IN1) DY12 = Y(IO2) - Y(IN1) DY22 = Y(IO2) - Y(IN2) DY21 = Y(IO1) - Y(IN2) C C COMPUTE INNER PRODUCTS C COS1 = DX11*DX12 + DY11*DY12 COS2 = DX22*DX21 + DY22*DY21 C C THE DIAGONALS SHOULD BE SWAPPED IFF (T1+T2) .GT. 180 C DEGREES. THE FOLLOWING TWO TESTS INSURE NUMERICAL C STABILITY. C IF (COS1 .GE. 0. .AND. COS2 .GE. 0.) RETURN IF (COS1 .LT. 0. .AND. COS2 .LT. 0.) GO TO 1 C C COMPUTE VECTOR CROSS PRODUCTS C SIN1 = DX11*DY12 - DX12*DY11 SIN2 = DX22*DY21 - DX21*DY22 SIN12 = SIN1*COS2 + COS1*SIN2 IF (SIN12 .GE. 0.) RETURN 1 SWPTST = .TRUE. RETURN END SUBROUTINE TRFIND (NST,PX,PY,X,Y,IADJ,IEND, I1,I2,I3) INTEGER NST, IADJ(1), IEND(1), I1, I2, I3 REAL PX, PY, X(1), Y(1) C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C THIS ROUTINE LOCATES A POINT P IN A THIESSEN TRIANGU- C LATION, RETURNING THE VERTEX INDICES OF A TRIANGLE WHICH C CONTAINS P. TRFIND IS PART OF AN INTERPOLATION PACKAGE C WHICH PROVIDES SUBROUTINES FOR CREATING THE MESH. C C INPUT PARAMETERS - NST - INDEX OF NODE AT WHICH TRFIND C BEGINS SEARCH. SEARCH TIME C DEPENDS ON THE PROXIMITY OF C NST TO P. C C PX,PY - X AND Y-COORDINATES OF THE C POINT TO BE LOCATED. C C X,Y - VECTORS OF COORDINATES OF C NODES IN THE MESH. (X(I),Y(I)) C DEFINES NODE I FOR I = 1,...,N C WHERE N .GE. 3. C C IADJ - SET OF ADJACENCY LISTS OF C NODES IN THE MESH. C C IEND - POINTERS TO THE ENDS OF C ADJACENCY LISTS IN IADJ FOR C EACH NODE IN THE MESH. C C IADJ AND IEND MAY BE CREATED BY TRMESH. C C INPUT PARAMETERS ARE NOT ALTERED BY THIS ROUTINE. C C OUTPUT PARAMETERS - I1,I2,I3 - VERTEX INDICES IN COUNTER- C CLOCKWISE ORDER - VERTICES C OF A TRIANGLE CONTAINING P C IF P IS AN INTERIOR NODE. C IF P IS OUTSIDE OF THE C BOUNDARY OF THE MESH, I1 C AND I2 ARE THE FIRST (RIGHT C -MOST) AND LAST (LEFTMOST) C NODES WHICH ARE VISIBLE C FROM P, AND I3 = 0. IF P C AND ALL OF THE NODES LIE ON C A SINGLE LINE THEN I1 = I2 C = I3 = 0. C C MODULES REFERENCED BY TRFIND - NONE C C INTRINSIC FUNCTION CALLED BY TRFIND - MAX0 C C*********************************************************** C INTEGER N0, N1, N2, N3, N4, INDX, IND, NF, . NL, NEXT REAL XP, YP LOGICAL LEFT C C LOCAL PARAMETERS - C C XP,YP = LOCAL VARIABLES CONTAINING PX AND PY C N0,N1,N2 = NODES IN COUNTERCLOCKWISE ORDER DEFINING A C CONE (WITH VERTEX N0) CONTAINING P C N3,N4 = NODES OPPOSITE N1-N2 AND N2-N1, RESPECTIVELY C INDX,IND = INDICES FOR IADJ C NF,NL = FIRST AND LAST NEIGHBORS OF N0 IN IADJ, OR C FIRST (RIGHTMOST) AND LAST (LEFTMOST) NODES C VISIBLE FROM P WHEN P IS OUTSIDE THE C BOUNDARY C NEXT = CANDIDATE FOR I1 OR I2 WHEN P IS OUTSIDE OF C THE BOUNDARY C LEFT = STATEMENT FUNCTION WHICH COMPUTES THE SIGN OF C A CROSS PRODUCT (Z-COMPONENT). LEFT(X1,..., C Y0) = .TRUE. IFF (X0,Y0) IS ON OR TO THE C LEFT OF THE VECTOR FROM (X1,Y1) TO (X2,Y2). C LEFT(X1,Y1,X2,Y2,X0,Y0) = (X2-X1)*(Y0-Y1) .GE. . (X0-X1)*(Y2-Y1) XP = PX YP = PY C C INITIALIZE VARIABLES AND FIND A CONE CONTAINING P C N0 = MAX0(NST,1) 1 INDX = IEND(N0) NL = IADJ(INDX) INDX = 1 IF (N0 .NE. 1) INDX = IEND(N0-1) + 1 NF = IADJ(INDX) N1 = NF IF (NL .NE. 0) GO TO 3 C C N0 IS A BOUNDARY NODE. SET NL TO THE LAST NONZERO C NEIGHBOR OF N0. C IND = IEND(N0) - 1 NL = IADJ(IND) IF ( LEFT(X(N0),Y(N0),X(NF),Y(NF),XP,YP) ) GO TO 2 C C P IS OUTSIDE THE BOUNDARY C NL = N0 GO TO 16 2 IF ( LEFT(X(NL),Y(NL),X(N0),Y(N0),XP,YP) ) GO TO 4 C C P IS OUTSIDE THE BOUNDARY AND N0 IS THE RIGHTMOST C VISIBLE BOUNDARY NODE C I1 = N0 GO TO 18 C C N0 IS AN INTERIOR NODE. FIND N1. C 3 IF ( LEFT(X(N0),Y(N0),X(N1),Y(N1),XP,YP) ) GO TO 4 INDX = INDX + 1 N1 = IADJ(INDX) IF (N1 .EQ. NL) GO TO 7 GO TO 3 C C P IS TO THE LEFT OF ARC N0-N1. INITIALIZE N2 TO THE NEXT C NEIGHBOR OF N0. C 4 INDX = INDX + 1 N2 = IADJ(INDX) IF ( .NOT. LEFT(X(N0),Y(N0),X(N2),Y(N2),XP,YP) ) . GO TO 8 N1 = N2 IF (N1 .NE. NL) GO TO 4 IF ( .NOT. LEFT(X(N0),Y(N0),X(NF),Y(NF),XP,YP) ) . GO TO 7 IF (XP .EQ. X(N0) .AND. YP .EQ. Y(N0)) GO TO 6 C C P IS LEFT OF OR ON ARCS N0-NB FOR ALL NEIGHBORS NB C OF N0. C ALL POINTS ARE COLLINEAR IFF P IS LEFT OF NB-N0 FOR C ALL NEIGHBORS NB OF N0. SEARCH THE NEIGHBORS OF N0 C IN REVERSE ORDER. NOTE -- N1 = NL AND INDX POINTS TO C NL. C 5 IF ( .NOT. LEFT(X(N1),Y(N1),X(N0),Y(N0),XP,YP) ) . GO TO 6 IF (N1 .EQ. NF) GO TO 20 INDX = INDX - 1 N1 = IADJ(INDX) GO TO 5 C C P IS TO THE RIGHT OF N1-N0, OR P=N0. SET N0 TO N1 AND C START OVER. C 6 N0 = N1 GO TO 1 C C P IS BETWEEN ARCS N0-N1 AND N0-NF C 7 N2 = NF C C P IS CONTAINED IN A CONE DEFINED BY LINE SEGMENTS N0-N1 C AND N0-N2 WHERE N1 IS ADJACENT TO N2 C 8 N3 = N0 9 IF ( LEFT(X(N1),Y(N1),X(N2),Y(N2),XP,YP) ) GO TO 13 C C SET N4 TO THE FIRST NEIGHBOR OF N2 FOLLOWING N1 C INDX = IEND(N2) IF (IADJ(INDX) .NE. N1) GO TO 10 C C N1 IS THE LAST NEIGHBOR OF N2. C SET N4 TO THE FIRST NEIGHBOR. C INDX = 1 IF (N2 .NE. 1) INDX = IEND(N2-1) + 1 N4 = IADJ(INDX) GO TO 11 C C N1 IS NOT THE LAST NEIGHBOR OF N2 C 10 INDX = INDX-1 IF (IADJ(INDX) .NE. N1) GO TO 10 N4 = IADJ(INDX+1) IF (N4 .NE. 0) GO TO 11 C C P IS OUTSIDE THE BOUNDARY C NF = N2 NL = N1 GO TO 16 C C DEFINE A NEW ARC N1-N2 WHICH INTERSECTS THE LINE C SEGMENT N0-P C 11 IF ( LEFT(X(N0),Y(N0),X(N4),Y(N4),XP,YP) ) GO TO 12 N3 = N2 N2 = N4 GO TO 9 12 N3 = N1 N1 = N4 GO TO 9 C C P IS IN THE TRIANGLE (N1,N2,N3) AND NOT ON N2-N3. IF C N3-N1 OR N1-N2 IS A BOUNDARY ARC CONTAINING P, TREAT P C AS EXTERIOR. C 13 INDX = IEND(N1) IF (IADJ(INDX) .NE. 0) GO TO 15 C C N1 IS A BOUNDARY NODE. N3-N1 IS A BOUNDARY ARC IFF N3 C IS THE LAST NONZERO NEIGHBOR OF N1. C IF (N3 .NE. IADJ(INDX-1)) GO TO 14 C C N3-N1 IS A BOUNDARY ARC C IF ( .NOT. LEFT(X(N1),Y(N1),X(N3),Y(N3),XP,YP) ) . GO TO 14 C C P LIES ON N1-N3 C I1 = N1 I2 = N3 I3 = 0 RETURN C C N3-N1 IS NOT A BOUNDARY ARC CONTAINING P. N1-N2 IS A C BOUNDARY ARC IFF N2 IS THE FIRST NEIGHBOR OF N1. C 14 INDX = 1 IF (N1 .NE. 1) INDX = IEND(N1-1) + 1 IF (N2 .NE. IADJ(INDX)) GO TO 15 C C N1-N2 IS A BOUNDARY ARC C IF ( .NOT. LEFT(X(N2),Y(N2),X(N1),Y(N1),XP,YP) ) . GO TO 15 C C P LIES ON N1-N2 C I1 = N2 I2 = N1 I3 = 0 RETURN C C P DOES NOT LIE ON A BOUNDARY ARC. C 15 I1 = N1 I2 = N2 I3 = N3 RETURN C C NF AND NL ARE ADJACENT BOUNDARY NODES WHICH ARE VISIBLE C FROM P. FIND THE FIRST VISIBLE BOUNDARY NODE. C SET NEXT TO THE FIRST NEIGHBOR OF NF. C 16 INDX = 1 IF (NF .NE. 1) INDX = IEND(NF-1) + 1 NEXT = IADJ(INDX) IF ( LEFT(X(NF),Y(NF),X(NEXT),Y(NEXT),XP,YP) ) . GO TO 17 NF = NEXT GO TO 16 C C NF IS THE FIRST (RIGHTMOST) VISIBLE BOUNDARY NODE C 17 I1 = NF C C FIND THE LAST VISIBLE BOUNDARY NODE. NL IS THE FIRST C CANDIDATE FOR I2. C SET NEXT TO THE LAST NEIGHBOR OF NL. C 18 INDX = IEND(NL) - 1 NEXT = IADJ(INDX) IF ( LEFT(X(NEXT),Y(NEXT),X(NL),Y(NL),XP,YP) ) . GO TO 19 NL = NEXT GO TO 18 C C NL IS THE LAST (LEFTMOST) VISIBLE BOUNDARY NODE C 19 I2 = NL I3 = 0 RETURN C C ALL POINTS ARE COLLINEAR C 20 I1 = 0 I2 = 0 I3 = 0 RETURN END SUBROUTINE TRMESH (N,X,Y, IADJ,IEND,IER) INTEGER N, IADJ(1), IEND(N), IER REAL X(N), Y(N) C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C THIS ROUTINE CREATES A THIESSEN TRIANGULATION OF N C ARBITRARILY SPACED POINTS IN THE PLANE REFERRED TO AS C NODES. THE TRIANGULATION IS OPTIMAL IN THE SENSE THAT IT C IS AS NEARLY EQUIANGULAR AS POSSIBLE. TRMESH IS PART OF C AN INTERPOLATION PACKAGE WHICH ALSO PROVIDES SUBROUTINES C TO REORDER THE NODES, ADD A NEW NODE, DELETE AN ARC, PLOT C THE MESH, AND PRINT THE DATA STRUCTURE. C UNLESS THE NODES ARE ALREADY ORDERED IN SOME REASONABLE C FASHION, THEY SHOULD BE REORDERED BY SUBROUTINE REORDR FOR C INCREASED EFFICIENCY BEFORE CALLING TRMESH. C C INPUT PARAMETERS - N - NUMBER OF NODES IN THE MESH. C N .GE. 3. C C X,Y - N-VECTORS OF COORDINATES. C (X(I),Y(I)) DEFINES NODE I. C C IADJ - VECTOR OF LENGTH .GE. 6*N-9. C C IEND - VECTOR OF LENGTH .GE. N. C C N, X, AND Y ARE NOT ALTERED BY THIS ROUTINE. C C OUTPUT PARAMETERS - IADJ - ADJACENCY LISTS OF NEIGHBORS IN C COUNTERCLOCKWISE ORDER. THE C LIST FOR NODE I+1 FOLLOWS THAT C FOR NODE I WHERE X AND Y DEFINE C THE ORDER. THE VALUE 0 DENOTES C THE BOUNDARY (OR A PSEUDO-NODE C AT INFINITY) AND IS ALWAYS THE C LAST NEIGHBOR OF A BOUNDARY C NODE. IADJ IS UNCHANGED IF IER C .NE. 0. C C IEND - POINTERS TO THE ENDS OF C ADJACENCY LISTS (SETS OF C NEIGHBORS) IN IADJ. THE C NEIGHBORS OF NODE 1 BEGIN IN C IADJ(1). FOR K .GT. 1, THE C NEIGHBORS OF NODE K BEGIN IN C IADJ(IEND(K-1)+1) AND K HAS C IEND(K) - IEND(K-1) NEIGHBORS C INCLUDING (POSSIBLY) THE C BOUNDARY. IADJ(IEND(K)) .EQ. 0 C IFF NODE K IS ON THE BOUNDARY. C IEND IS UNCHANGED IF IER = 1. C IF IER = 2 IEND CONTAINS THE C INDICES OF A SEQUENCE OF N C NODES ORDERED FROM LEFT TO C RIGHT WHERE LEFT AND RIGHT ARE C DEFINED BY ASSUMING NODE 1 IS C TO THE LEFT OF NODE 2. C C IER - ERROR INDICATOR C IER = 0 IF NO ERRORS WERE C ENCOUNTERED. C IER = 1 IF N .LT. 3. C IER = 2 IF N .GE. 3 AND ALL C NODES ARE COLLINEAR. C C MODULES REFERENCED BY TRMESH - SHIFTD, ADNODE, TRFIND, C INTADD, BDYADD, SWPTST, C SWAP, INDEX C C*********************************************************** C INTEGER NN, K, KM1, NL, NR, IND, INDX, N0, ITEMP, . IERR, KM1D2, KMI, I, KMIN REAL XL, YL, XR, YR, DXR, DYR, XK, YK, DXK, DYK, . CPROD, SPROD C C LOCAL PARAMETERS - C C NN = LOCAL COPY OF N C K = NODE (INDEX) TO BE INSERTED INTO IEND C KM1 = K-1 - (VARIABLE) LENGTH OF IEND C NL,NR = IEND(1), IEND(KM1) -- LEFTMOST AND RIGHTMOST C NODES IN IEND AS VIEWED FROM THE RIGHT OF C 1-2 WHEN IEND CONTAINS THE INITIAL ORDERED C SET OF NODAL INDICES C XL,YL,XR,YR = X AND Y COORDINATES OF NL AND NR C DXR,DYR = XR-XL, YR-YL C XK,YK = X AND Y COORDINATES OF NODE K C DXK,DYK = XK-XL, YK-YL C CPROD = VECTOR CROSS PRODUCT OF NL-NR AND NL-K -- C USED TO DETERMINE THE POSITION OF NODE K C WITH RESPECT TO THE LINE DEFINED BY THE C NODES IN IEND C SPROD = SCALAR PRODUCT USED TO DETERMINE THE C INTERVAL CONTAINING NODE K WHEN K IS ON C THE LINE DEFINED BY THE NODES IN IEND C IND,INDX = INDICES FOR IEND AND IADJ, RESPECTIVELY C N0,ITEMP = TEMPORARY NODES (INDICES) C IERR = DUMMY PARAMETER FOR CALL TO ADNODE C KM1D2,KMI,I = KM1/2, K-I, DO-LOOP INDEX -- USED IN IEND C REORDERING LOOP C KMIN = FIRST NODE INDEX SENT TO ADNODE C NN = N IER = 1 IF (NN .LT. 3) RETURN IER = 0 C C INITIALIZE IEND, NL, NR, AND K C IEND(1) = 1 IEND(2) = 2 XL = X(1) YL = Y(1) XR = X(2) YR = Y(2) K = 2 C C BEGIN LOOP ON NODES 3,4,... C 1 DXR = XR-XL DYR = YR-YL C C NEXT LOOP BEGINS HERE IF NL AND NR ARE UNCHANGED C 2 IF (K .EQ. NN) GO TO 13 KM1 = K K = KM1 + 1 XK = X(K) YK = Y(K) DXK = XK-XL DYK = YK-YL CPROD = DXR*DYK - DXK*DYR IF (CPROD .GT. 0.) GO TO 6 IF (CPROD .LT. 0.) GO TO 8 C C NODE K LIES ON THE LINE CONTAINING NODES 1,2,...,K-1. C SET SPROD TO (NL-NR,NL-K). C SPROD = DXR*DXK + DYR*DYK IF (SPROD .GT. 0.) GO TO 3 C C NODE K IS TO THE LEFT OF NL. INSERT K AS THE FIRST C (LEFTMOST) NODE IN IEND AND SET NL TO K. C CALL SHIFTD(1,KM1,1, IEND ) IEND(1) = K XL = XK YL = YK GO TO 1 C C NODE K IS TO THE RIGHT OF NL. FIND THE LEFTMOST NODE C N0 WHICH LIES TO THE RIGHT OF K. C SET SPROD TO (N0-NL,N0-K). C 3 DO 4 IND = 2,KM1 N0 = IEND(IND) SPROD = (XL-X(N0))*(XK-X(N0)) + . (YL-Y(N0))*(YK-Y(N0)) IF (SPROD .GE. 0.) GO TO 5 4 CONTINUE C C NODE K IS TO THE RIGHT OF NR. INSERT K AS THE LAST C (RIGHTMOST) NODE IN IEND AND SET NR TO K. C IEND(K) = K XR = XK YR = YK GO TO 1 C C NODE K LIES BETWEEN IEND(IND-1) AND IEND(IND). INSERT K C IN IEND. C 5 CALL SHIFTD(IND,KM1,1, IEND ) IEND(IND) = K GO TO 2 C C NODE K IS TO THE LEFT OF NL-NR. REORDER IEND SO THAT NL C IS THE LEFTMOST NODE AS VIEWED FROM K. C 6 KM1D2 = KM1/2 DO 7 I = 1,KM1D2 KMI = K-I ITEMP = IEND(I) IEND(I) = IEND(KMI) IEND(KMI) = ITEMP 7 CONTINUE C C NODE K IS TO THE RIGHT OF NL-NR. CREATE A TRIANGULATION C CONSISTING OF NODES 1,2,...,K. C 8 NL = IEND(1) NR = IEND(KM1) C C CREATE THE ADJACENCY LISTS FOR THE FIRST K-1 NODES. C INSERT NEIGHBORS IN REVERSE ORDER. EACH NODE HAS FOUR C NEIGHBORS EXCEPT NL AND NR WHICH HAVE THREE. C DO 9 IND = 1,KM1 N0 = IEND(IND) INDX = 4*N0 IF (N0 .GE. NL) INDX = INDX-1 IF (N0 .GE. NR) INDX = INDX-1 IADJ(INDX) = 0 INDX = INDX-1 IF (IND .LT. KM1) IADJ(INDX) = IEND(IND+1) IF (IND .LT. KM1) INDX = INDX-1 IADJ(INDX) = K IF (IND .EQ. 1) GO TO 9 IADJ(INDX-1) = IEND(IND-1) 9 CONTINUE C C CREATE THE ADJACENCY LIST FOR NODE K C INDX = 5*KM1 - 1 IADJ(INDX) = 0 DO 10 IND = 1,KM1 INDX = INDX-1 IADJ(INDX) = IEND(IND) 10 CONTINUE C C REPLACE IEND ELEMENTS WITH POINTERS TO IADJ C INDX = 0 DO 11 IND = 1,KM1 INDX = INDX + 4 IF (IND .EQ. NL .OR. IND .EQ. NR) INDX = INDX-1 IEND(IND) = INDX 11 CONTINUE INDX = INDX + K IEND(K) = INDX C C ADD THE REMAINING NODES TO THE TRIANGULATION C IF (K .EQ. NN) RETURN KMIN = K+1 DO 12 K = KMIN,NN CALL ADNODE(K,X,Y, IADJ,IEND, IERR) 12 CONTINUE RETURN C C ALL NODES ARE COLLINEAR C 13 IER = 2 RETURN END SUBROUTINE TRPLOT (N,X,Y,IADJ,IEND,ITITLE,NC, . NUMBR, IER) INTEGER N, IADJ(1), IEND(N), ITITLE(1), NUMBR, IER REAL X(N), Y(N) C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C THIS SUBROUTINE PLOTS THE ARCS OF A TRIANGULATION OF N C NODES IN THE PLANE. CARDS WITH C* IN THE FIRST TWO COL- C UMNS MUST BE REPLACED WITH CALLS TO USER-SUPPLIED GRAPHICS C SUBROUTINES IN ORDER TO GET ANY USE OUT OF THIS ROUTINE. C C INPUT PARAMETERS - N - NUMBER OF NODES IN THE TRIANGULA- C TION. N .GE. 3. C C X,Y - CARTESIAN COORDINATES OF THE NODES. C C IADJ,IEND - TRIANGULATION DATA STRUCTURE (SEE C SUBROUTINE TRMESH). C C ITITLE - INTEGER ARRAY CONTAINING A LINE OF C TEXT TO BE CENTERED ABOVE THE PLOT C IF NC .GT. 0. ITITLE MUST BE INIT- C IALIZED WITH HOLLERITH CONSTANTS OR C READ WITH AN A-FORMAT. ITS DIMEN- C SION DEPENDS ON NC AND THE NUMBER C OF CHARACTERS STORED IN A COMPUTER C WORD. C C NC - NUMBER OF CHARACTERS IN ITITLE. C 0 .LE. NC .LE. 40. NO TITLE IS C PLOTTED IF NC = 0. C C NUMBR - OPTION INDICATOR. IF NUMBR .NE. 0, C THE NODAL INDICES ARE PLOTTTED NEXT C TO THE NODES. C C INPUT PARAMETERS ARE NOT ALTERED BY THIS ROUTINE. C C OUTPUT PARAMETER - IER - ERROR INDICATOR C IER = 0 IF NO ERRORS WERE ENCOUN- C TERED. C IER = 1 IF N OR NC IS OUT OF C RANGE. C IER = 2 IF THE NODES LIE ON A C HORIZONTAL OR VERTICAL C LINE. C C MODULES REFERENCED BY TRPLOT - NONE C C INTRINSIC FUNCTIONS CALLED BY TRPLOT - AMIN1, AMAX1, IABS C C*********************************************************** C INTEGER NN, I, N0, INDF, INDL, INDX, N1, NEW REAL XMIN, XMAX, YMIN, YMAX, DX, DY, X0, Y0 LOGICAL ST0 C C LOCAL PARAMETERS - C C NN = LOCAL COPY OF N C I = DO-LOOP INDEX C N0 = NODE WHICH IS TO BE CONNECTED TO ITS NEIGHBORS C WITH LINE SEGMENTS C INDF = IADJ INDEX OF THE FIRST NEIGHBOR OF N0 C INDL = IADJ INDEX OF THE LAST NEIGHBOR OF N0 C INDX = IADJ INDEX IN THE RANGE INDF,...,INDL C N1 = NEIGHBOR OF N0 C NEW = NEIGHBOR OF N0 AND CANDIDATE FOR NEXT VALUE C OF N0 C XMIN,XMAX = MINIMUM AND MAXIMUM X COORDINATES C YMIN,YMAX = MINIMUM AND MAXIMUM Y COORDINATES C DX,DY = XMAX-XMIN AND YMAX-YMIN (DATA SPACE DIMEN- C SIONS) C X0,Y0 = COORDINATES OF N0 C ST0 = SWITCH USED TO ALTERNATE DIRECTION OF PEN C MOVEMENT -- TRUE IFF PEN STARTS AT N0 C NN = N C C CHECK FOR INVALID PARAMETERS C IF (NN .LT. 3 .OR. NC .LT. 0 .OR. NC .GT. 40) . GO TO 13 IER = 0 C C COMPUTE DATA SPACE DIMENSIONS C XMIN = X(1) YMIN = Y(1) XMAX = XMIN YMAX = YMIN DO 1 I = 1,NN XMIN = AMIN1(X(I),XMIN) YMIN = AMIN1(Y(I),YMIN) XMAX = AMAX1(X(I),XMAX) 1 YMAX = AMAX1(Y(I),YMAX) DX = XMAX - XMIN DY = YMAX - YMIN IF (DX .EQ. 0. .OR. DY .EQ. 0.) GO TO 14 C* COMMANDS WHICH PERFORM THE FOLLOWING FUNCTIONS SHOULD C* BE INSERTED HERE -- C* INITIALIZE THE PLOTTING ENVIRONMENT IF NECESSARY, C* COMPUTE PLOTTER SPACE DIMENSIONS, C* ESTABLISH A LINEAR MAPPING FROM THE DATA SPACE TO C* THE PLOTTER SPACE, C* OPTIONALLY, DRAW AND LABEL AXES, AND C* PLOT THE TITLE IF NC .NE. 0. C C INITIALIZE FOR LOOP ON NODES. EACH NODE N0 IS TO BE CON- C NECTED TO ITS NEIGHBORS WHICH HAVE LARGER INDICES. N0 C IS THEN MARKED BY MAKING THE CORRESPONDING IEND ENTRY C NEGATIVE, AND THE SEARCH FOR THE NEXT UNMARKED NODE BE- C GINS WITH THE NEIGHBORS OF N0. C N0 = 1 INDF = 1 C C TOP OF LOOP -- SET INDL, X0, AND Y0, AND INITIALIZE ST0 C AND INDX C 2 INDL = IEND(N0) X0 = X(N0) Y0 = Y(N0) ST0 = .TRUE. INDX = INDL IF (IADJ(INDX) .EQ. 0) INDX = INDX - 1 C C LOOP ON NEIGHBORS OF N0 IN REVERSE ORDER C 3 N1 = IADJ(INDX) IF (N1 .LT. N0) GO TO 4 C C CONNECT N0 AND N1 -- THE DIRECTION OF PEN MOVEMENT C ALTERNATES BETWEEN AWAY FROM N0 AND TOWARD N0 FOR C REDUCED PEN-UP TIME. C C* IF (ST0) CALL LINE (X0,Y0,X(N1),Y(N1)) C* IF (.NOT. ST0) CALL LINE (X(N1),Y(N1),X0,Y0) ST0 = .NOT. ST0 C C TEST FOR TERMINATION OF LOOP ON NEIGHBORS C 4 IF (INDX .EQ. INDF) GO TO 5 INDX = INDX - 1 GO TO 3 C C MARK N0 AS HAVING BEEN PROCESSED C 5 IEND(N0) = -INDL C C SEARCH THE NEIGHBORS OF N0 FOR AN UNMARKED NODE C STARTING WITH IADJ(INDX) = IADJ(INDF) C 6 NEW = IADJ(INDX) IF (NEW .EQ. 0) GO TO 8 IF (IEND(NEW) .LT. 0) GO TO 7 N0 = NEW INDF = IABS(IEND(N0-1)) + 1 GO TO 2 C C TEST FOR TERMINATION C 7 IF (INDX .EQ. INDL) GO TO 8 INDX = INDX + 1 GO TO 6 C C ALL NEIGHBORS OF N0 ARE MARKED. SEARCH IEND FOR AN C UNMARKED NODE. C 8 DO 9 N0 = 2,NN IF (IEND(N0) .LT. 0) GO TO 9 INDF = -IEND(N0-1) + 1 GO TO 2 9 CONTINUE C C ALL NODES HAVE BEEN MARKED. RESTORE IEND. C DO 10 N0 = 1,NN 10 IEND(N0) = -IEND(N0) IF (NUMBR .EQ. 0) GO TO 12 C C PLOT THE NODAL INDICES C C* THE NUMBR OPTION MAY BE IMPLEMENTED BY INSERTING A C* LOOP ON THE NODAL COORDINATES WHICH WRITES INDICES C* NEXT TO THE NODES. C C TERMINATE PLOTTING -- MOVE TO A NEW FRAME. C 12 CONTINUE C* CALL FRAME RETURN C C N OR NC IS OUT OF RANGE C 13 IER = 1 RETURN C C NODES ARE COLLINEAR C 14 IER = 2 RETURN END SUBROUTINE TRPRNT (N,LUNIT,X,Y,IADJ,IEND,IFLAG) INTEGER N, LUNIT, IADJ(1), IEND(N), IFLAG REAL X(N), Y(N) C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C GIVEN A TRIANGULATION OF A SET OF POINTS IN THE PLANE, C THIS ROUTINE PRINTS THE ADJACENCY LISTS AND, OPTIONALLY, C THE NODAL COORDINATES. THE NUMBERS OF BOUNDARY NODES, C TRIANGLES, AND ARCS ARE ALSO PRINTED. C C INPUT PARAMETERS - N - NUMBER OF NODES IN THE MESH. C 3 .LE. N .LE. 9999. C C LUNIT - LOGICAL UNIT FOR OUTPUT. 1 C .LE. LUNIT .LE. 99. OUTPUT IS C PRINTED ON UNIT 6 IF LUNIT IS C OUT OF RANGE. C C X,Y - VECTORS OF COORDINATES OF THE C NODES IN THE MESH. NOT USED C UNLESS IFLAG = 0. C C IADJ - SET OF ADJACENCY LISTS OF NODES C IN THE MESH. C C IEND - POINTERS TO THE ENDS OF C ADJACENCY LISTS IN IADJ FOR C EACH NODE IN THE MESH. C C IFLAG - OPTION INDICATOR C IFLAG = 0 IF X AND Y ARE TO BE C PRINTED (TO 6 DECIMAL C PLACES). C IFLAG = 1 IF X AND Y ARE NOT C TO BE PRINTED. C C IADJ AND IEND MAY BE CREATED BY TRMESH. C C NONE OF THE PARAMETERS ARE ALTERED BY THIS ROUTINE. C C MODULES REFERENCED BY TRPRNT - NONE C C*********************************************************** C INTEGER NN, NMAX, LUN, NODE, INDF, INDL, NL, NLMAX, . INC, I, NB, NT, NA DATA NMAX/9999/, NLMAX/60/ C C LOCAL PARAMETERS - C C NN = LOCAL COPY OF N C NMAX = UPPER BOUND ON N C LUN = LOCAL COPY OF LUNIT C NODE = INDEX OF A NODE C INDF,INDL = IADJ INDICES OF THE FIRST AND LAST NEIGHBORS C OF NODE C NL = NUMBER OF LINES PRINTED ON A PAGE C NLMAX = MAXIMUM NUMBER OF PRINT LINES PER PAGE EXCEPT C FOR THE LAST PAGE WHICH HAS 3 ADDITIONAL C LINES C INC = INCREMENT FOR NL C I = IADJ INDEX FOR IMPLIED DO-LOOP C NB = NUMBER OF BOUNDARY NODES C NT = NUMBER OF TRIANGLES C NA = NUMBER OF ARCS (UNDIRECTED EDGES) C NN = N LUN = LUNIT IF (LUN .LT. 1 .OR. LUN .GT. 99) LUN = 6 C C PRINT HEADING AND INITIALIZE NL C WRITE (LUN,100) NN IF (NN .LT. 3 .OR. NN .GT. NMAX) GO TO 5 NL = 6 IF (IFLAG .EQ. 0) GO TO 2 C C PRINT IADJ ONLY C WRITE (LUN,101) NB = 0 INDF = 1 DO 1 NODE = 1,NN INDL = IEND(NODE) IF (IADJ(INDL) .EQ. 0) NB = NB + 1 INC = (INDL - INDF)/14 + 2 NL = NL + INC IF (NL .GT. NLMAX) WRITE (LUN,106) IF (NL .GT. NLMAX) NL = INC WRITE (LUN,103) NODE, (IADJ(I), I = INDF,INDL) IF (INDL-INDF .NE. 13) WRITE (LUN,105) INDF = INDL + 1 1 CONTINUE GO TO 4 C C PRINT X, Y, AND IADJ C 2 WRITE (LUN,102) NB = 0 INDF = 1 DO 3 NODE = 1,NN INDL = IEND(NODE) IF (IADJ(INDL) .EQ. 0) NB = NB + 1 INC = (INDL - INDF)/8 + 2 NL = NL + INC IF (NL .GT. NLMAX) WRITE (LUN,106) IF (NL .GT. NLMAX) NL = INC WRITE (LUN,104) NODE, X(NODE), Y(NODE), . (IADJ(I), I = INDF,INDL) IF (INDL-INDF .NE. 7) WRITE (LUN,105) INDF = INDL + 1 3 CONTINUE C C PRINT NB, NA, AND NT C 4 NT = 2*NN - NB - 2 NA = NT + NN - 1 WRITE (LUN,107) NB, NA, NT RETURN C C N IS OUT OF RANGE C 5 WRITE (LUN,108) RETURN C C PRINT FORMATS C 100 FORMAT (1H1,26X,23HADJACENCY SETS, N = ,I5//) 101 FORMAT (1H ,4HNODE,32X,17HNEIGHBORS OF NODE//) 102 FORMAT (1H ,4HNODE,5X,7HX(NODE),8X,7HY(NODE), . 20X,17HNEIGHBORS OF NODE//) 103 FORMAT (1H ,I4,5X,14I5/(1H ,9X,14I5)) 104 FORMAT (1H ,I4,2E15.6,5X,8I5/(1H ,39X,8I5)) 105 FORMAT (1H ) 106 FORMAT (1H1) 107 FORMAT (/1H ,5HNB = ,I4,15H BOUNDARY NODES,10X, . 5HNA = ,I5,5H ARCS,10X,5HNT = ,I5, . 10H TRIANGLES) 108 FORMAT (1H ,10X,25H*** N IS OUT OF RANGE ***) END SUBROUTINE COORDS (X,Y,X1,X2,X3,Y1,Y2,Y3, R,IER) INTEGER IER REAL X, Y, X1, X2, X3, Y1, Y2, Y3, R(3) C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C THIS ROUTINE COMPUTES THE THREE BARYCENTRIC COORDINATES C OF A POINT IN THE PLANE FOR A GIVEN TRIANGLE. C C INPUT PARAMETERS - X,Y - X AND Y COORDINATES OF THE POINT C WHOSE BARYCENTRIC COORDINATES ARE C DESIRED. C C X1,X2,X3,Y1,Y2,Y3 - COORDINATES OF THE VERTICES OF C THE TRIANGLE. C C INPUT PARAMETERS ARE NOT ALTERED BY THIS ROUTINE. C C OUTPUT PARAMETERS - R - 3-VECTOR OF BARYCENTRIC COORDI- C NATES UNLESS IER = 1. NOTE THAT C R(I) .LT. 0. IFF (X,Y) IS TO THE C RIGHT OF THE VECTOR FROM VERTEX C I+1 TO VERTEX I+2 (CYCLICAL C ARITHMETIC). C C IER - ERROR INDICATOR C IER = 0 IF NO ERRORS WERE C ENCOUNTERED. C IER = 1 IF THE VERTICES OF THE C TRIANGLE ARE COLLINEAR. C C MODULES REFERENCED BY COORDS - NONE C C*********************************************************** C REAL U(3), V(3), AREA, XP, YP C C LOCAL PARAMETERS - C C U(K),V(K) = X AND Y COMPONENTS OF THE VECTOR REPRESENTING C THE SIDE OPPOSITE VERTEX K FOR K = 1,2,3. C AREA = TWICE THE AREA OF THE TRIANGLE. C XP,YP = X-X1, Y-Y1 C U(1) = X3-X2 U(2) = X1-X3 U(3) = X2-X1 C V(1) = Y3-Y2 V(2) = Y1-Y3 V(3) = Y2-Y1 C C AREA = 3-1 X 3-2 C AREA = U(1)*V(2) - U(2)*V(1) IF (AREA .EQ. 0.) GO TO 1 C C R(1) = (2-3 X 2-(X,Y))/AREA, R(2) = (1-(X,Y) X 1-3)/AREA, C R(3) = (1-2 X 1-(X,Y))/AREA C R(1) = (U(1)*(Y-Y2) - V(1)*(X-X2))/AREA XP = X - X1 YP = Y - Y1 R(2) = (U(2)*YP - V(2)*XP)/AREA R(3) = (U(3)*YP - V(3)*XP)/AREA IER = 0 RETURN C C VERTICES ARE COLLINEAR C 1 IER = 1 RETURN END SUBROUTINE GIVENS ( A,B, C,S) REAL A, B, C, S C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C THIS ROUTINE CONSTRUCTS THE GIVENS PLANE ROTATION -- C ( C S) C G = ( ) WHERE C*C + S*S = 1 -- WHICH ZEROS THE SECOND C (-S C) C ENTRY OF THE 2-VECTOR (A B)-TRANSPOSE. A CALL TO GIVENS C IS NORMALLY FOLLOWED BY A CALL TO ROTATE WHICH APPLIES C THE TRANSFORMATION TO A 2 BY N MATRIX. THIS ROUTINE WAS C TAKEN FROM LINPACK. C C INPUT PARAMETERS - A,B - COMPONENTS OF THE 2-VECTOR TO BE C ROTATED. C C OUTPUT PARAMETERS - A - OVERWRITTEN BY R = +/-SQRT(A*A C + B*B) C C B - OVERWRITTEN BY A VALUE Z WHICH C ALLOWS C AND S TO BE RECOVERED C AS FOLLOWS - C C = SQRT(1-Z*Z), S=Z IF ABS(Z) C .LE. 1. C C = 1/Z, S = SQRT(1-C*C) IF C ABS(Z) .GT. 1. C C C - +/-(A/R) C C S - +/-(B/R) C C MODULES REFERENCED BY GIVENS - NONE C C INTRINSIC FUNCTIONS CALLED BY GIVENS - ABS, SQRT C C*********************************************************** C REAL AA, BB, R, U, V C C LOCAL PARAMETERS - C C AA,BB = LOCAL COPIES OF A AND B C R = C*A + S*B = +/-SQRT(A*A+B*B) C U,V = VARIABLES USED TO SCALE A AND B FOR COMPUTING R C AA = A BB = B IF (ABS(AA) .LE. ABS(BB)) GO TO 1 C C ABS(A) .GT. ABS(B) C U = AA + AA V = BB/U R = SQRT(.25 + V*V) * U C = AA/R S = V * (C + C) C C NOTE THAT R HAS THE SIGN OF A, C .GT. 0, AND S HAS C SIGN(A)*SIGN(B) C B = S A = R RETURN C C ABS(A) .LE. ABS(B) C 1 IF (BB .EQ. 0.) GO TO 2 U = BB + BB V = AA/U C C STORE R IN A C A = SQRT(.25 + V*V) * U S = BB/A C = V * (S + S) C C NOTE THAT R HAS THE SIGN OF B, S .GT. 0, AND C HAS C SIGN(A)*SIGN(B) C B = 1. IF (C .NE. 0.) B = 1./C RETURN C C A = B = 0. C 2 C = 1. S = 0. RETURN END SUBROUTINE GRADG (N,X,Y,Z,IADJ,IEND,EPS, NIT, . ZXZY, IER) INTEGER N, IADJ(1), IEND(N), NIT, IER REAL X(N), Y(N), Z(N), EPS, ZXZY(2,N) C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C GIVEN A TRIANGULATION OF N NODES IN THE PLANE WITH C ASSOCIATED DATA VALUES, THIS ROUTINE USES A GLOBAL METHOD C TO COMPUTE ESTIMATED GRADIENTS AT THE NODES. THE METHOD C CONSISTS OF MINIMIZING A QUADRATIC FUNCTIONAL Q(G) OVER C THE N-VECTOR G OF GRADIENTS WHERE Q APPROXIMATES THE LIN- C EARIZED CURVATURE OF AN INTERPOLANT F OVER THE TRIANGULA- C TION. THE RESTRICTION OF F TO AN ARC OF THE TRIANGULATION C IS TAKEN TO BE THE HERMITE CUBIC INTERPOLANT OF THE DATA C VALUES AND TANGENTIAL GRADIENT COMPONENTS AT THE END- C POINTS OF THE ARC, AND Q IS THE SUM OF THE LINEARIZED C CURVATURES OF F ALONG THE ARCS -- THE INTEGRALS OVER THE C ARCS OF D2F(T)**2 WHERE D2F(T) IS THE SECOND DERIVATIVE C OF F WITH RESPECT TO DISTANCE T ALONG THE ARC. THIS MIN- C IMIZATION PROBLEM CORRESPONDS TO AN ORDER 2N SYMMETRIC C POSITIVE-DEFINITE SPARSE LINEAR SYSTEM WHICH IS SOLVED FOR C THE X AND Y PARTIAL DERIVATIVES BY THE BLOCK GAUSS-SEIDEL C METHOD WITH 2 BY 2 BLOCKS. C AN ALTERNATIVE METHOD, SUBROUTINE GRADL, COMPUTES A C LOCAL APPROXIMATION TO THE PARTIALS AT A SINGLE NODE AND C MAY BE MORE ACCURATE, DEPENDING ON THE DATA VALUES AND C DISTRIBUTION OF NODES (NEITHER METHOD EMERGED AS SUPERIOR C IN TESTS FOR ACCURACY). HOWEVER, IN TESTS RUN ON AN IBM C 370, GRADG WAS FOUND TO BE ABOUT 3.6 TIMES AS FAST FOR C NIT = 4. C C INPUT PARAMETERS - N - NUMBER OF NODES. N .GE. 3. C C X,Y - CARTESIAN COORDINATES OF THE NODES. C C Z - DATA VALUES AT THE NODES. Z(I) IS C ASSOCIATED WITH (X(I),Y(I)). C C IADJ,IEND - DATA STRUCTURE DEFINING THE TRIAN- C GULATION. SEE SUBROUTINE TRMESH. C C EPS - NONNEGATIVE CONVERGENCE CRITERION. C THE METHOD IS TERMINATED WHEN THE C MAXIMUM CHANGE IN A GRADIENT COMPO- C NENT BETWEEN ITERATIONS IS AT MOST C EPS. EPS = 1.E-2 IS SUFFICIENT FOR C EFFECTIVE CONVERGENCE. C C THE ABOVE PARAMETERS ARE NOT ALTERED BY THIS ROUTINE. C C NIT - MAXIMUM NUMBER OF GAUSS-SEIDEL C ITERATIONS TO BE APPLIED. THIS C MAXIMUM WILL LIKELY BE ACHIEVED IF C EPS IS SMALLER THAN THE MACHINE C PRECISION. OPTIMAL EFFICIENCY WAS C ACHIEVED IN TESTING WITH EPS = 0 C AND NIT = 3 OR 4. C C ZXZY - 2 BY N ARRAY WHOSE COLUMNS CONTAIN C INITIAL ESTIMATES OF THE PARTIAL C DERIVATIVES (ZERO VECTORS ARE C SUFFICIENT). C C OUTPUT PARAMETERS - NIT - NUMBER OF GAUSS-SEIDEL ITERA- C TIONS EMPLOYED. C C ZXZY - ESTIMATED X AND Y PARTIAL DERIV- C ATIVES AT THE NODES WITH X PAR- C TIALS IN THE FIRST ROW. ZXZY IS C NOT CHANGED IF IER = 2. C C IER - ERROR INDICATOR C IER = 0 IF THE CONVERGENCE CRI- C TERION WAS ACHIEVED. C IER = 1 IF CONVERGENCE WAS NOT C ACHIEVED WITHIN NIT C ITERATIONS. C IER = 2 IF N OR EPS IS OUT OF C RANGE OR NIT .LT. 0 ON C INPUT. C C MODULES REFERENCED BY GRADG - NONE C C INTRINSIC FUNCTIONS CALLED BY GRADG - SQRT, AMAX1, ABS C C*********************************************************** C INTEGER NN, MAXIT, ITER, K, INDF, INDL, INDX, NB REAL TOL, DGMAX, XK, YK, ZK, ZXK, ZYK, A11, A12, . A22, R1, R2, DELX, DELY, DELXS, DELYS, DSQ, . DCUB, T, DZX, DZY C C LOCAL PARAMETERS - C C NN = LOCAL COPY OF N C MAXIT = INPUT VALUE OF NIT C ITER = NUMBER OF ITERATIONS USED C K = DO-LOOP AND NODE INDEX C INDF,INDL = IADJ INDICES OF THE FIRST AND LAST NEIGHBORS C OF K C INDX = IADJ INDEX IN THE RANGE INDF,...,INDL C NB = NEIGHBOR OF K C TOL = LOCAL COPY OF EPS C DGMAX = MAXIMUM CHANGE IN A GRADIENT COMPONENT BE- C TWEEN ITERATIONS C XK,YK,ZK = X(K), Y(K), Z(K) C ZXK,ZYK = INITIAL VALUES OF ZXZY(1,K) AND ZXZY(2,K) C A11,A12,A22 = MATRIX COMPONENTS OF THE 2 BY 2 BLOCK A*DG C = R WHERE A IS SYMMETRIC, DG = (DZX,DZY) C IS THE CHANGE IN THE GRADIENT AT K, AND R C IS THE RESIDUAL C R1,R2 = COMPONENTS OF THE RESIDUAL -- DERIVATIVES OF C Q WITH RESPECT TO THE COMPONENTS OF THE C GRADIENT AT NODE K C DELX,DELY = COMPONENTS OF THE ARC NB-K C DELXS,DELYS = DELX**2, DELY**2 C DSQ = SQUARE OF THE DISTANCE D BETWEEN K AND NB C DCUB = D**3 C T = FACTOR OF R1 AND R2 C DZX,DZY = SOLUTION OF THE 2 BY 2 SYSTEM -- CHANGE IN C DERIVATIVES AT K FROM THE PREVIOUS ITERATE C NN = N TOL = EPS MAXIT = NIT C C ERROR CHECKS AND INITIALIZATION C IF (NN .LT. 3 .OR. TOL .LT. 0. .OR. MAXIT .LT. 0) . GO TO 5 ITER = 0 C C TOP OF ITERATION LOOP C 1 IF (ITER .EQ. MAXIT) GO TO 4 DGMAX = 0. INDL = 0 DO 3 K = 1,NN XK = X(K) YK = Y(K) ZK = Z(K) ZXK = ZXZY(1,K) ZYK = ZXZY(2,K) C C INITIALIZE COMPONENTS OF THE 2 BY 2 SYSTEM C A11 = 0. A12 = 0. A22 = 0. R1 = 0. R2 = 0. C C LOOP ON NEIGHBORS NB OF K C INDF = INDL + 1 INDL = IEND(K) DO 2 INDX = INDF,INDL NB = IADJ(INDX) IF (NB .EQ. 0) GO TO 2 C C COMPUTE THE COMPONENTS OF ARC NB-K C DELX = X(NB) - XK DELY = Y(NB) - YK DELXS = DELX*DELX DELYS = DELY*DELY DSQ = DELXS + DELYS DCUB = DSQ*SQRT(DSQ) C C UPDATE THE SYSTEM COMPONENTS FOR NODE NB C A11 = A11 + DELXS/DCUB A12 = A12 + DELX*DELY/DCUB A22 = A22 + DELYS/DCUB T = ( 1.5*(Z(NB)-ZK) - ((ZXZY(1,NB)/2.+ZXK)*DELX + . (ZXZY(2,NB)/2.+ZYK)*DELY) )/DCUB R1 = R1 + T*DELX R2 = R2 + T*DELY 2 CONTINUE C C SOLVE THE 2 BY 2 SYSTEM AND UPDATE DGMAX C DZY = (A11*R2 - A12*R1)/(A11*A22 - A12*A12) DZX = (R1 - A12*DZY)/A11 DGMAX = AMAX1(DGMAX,ABS(DZX),ABS(DZY)) C C UPDATE THE PARTIALS AT NODE K C ZXZY(1,K) = ZXK + DZX 3 ZXZY(2,K) = ZYK + DZY C C INCREMENT ITER AND TEST FOR CONVERGENCE C ITER = ITER + 1 IF (DGMAX .GT. TOL) GO TO 1 C C METHOD CONVERGED C NIT = ITER IER = 0 RETURN C C METHOD FAILED TO CONVERGE WITHIN NIT ITERATIONS C 4 IER = 1 RETURN C C PARAMETER OUT OF RANGE C 5 NIT = 0 IER = 2 RETURN END SUBROUTINE GRADL (N,K,X,Y,Z,IADJ,IEND, DX,DY,IER) INTEGER N, K, IADJ(1), IEND(N), IER REAL X(N), Y(N), Z(N), DX, DY C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C GIVEN A THIESSEN TRIANGULATION OF N POINTS IN THE PLANE C WITH ASSOCIATED DATA VALUES Z, THIS SUBROUTINE ESTIMATES C X AND Y PARTIAL DERIVATIVES AT NODE K. THE DERIVATIVES C ARE TAKEN TO BE THE PARTIALS AT K OF A QUADRATIC FUNCTION C WHICH INTERPOLATES Z(K) AND FITS THE DATA VALUES AT A SET C OF NEARBY NODES IN A WEIGHTED LEAST SQUARES SENSE. A MAR- C QUARDT STABILIZATION FACTOR IS USED IF NECESSARY TO ENSURE C A WELL-CONDITIONED SYSTEM AND A LINEAR FITTING FUNCTION IS C USED IF N .LT. 6. THUS, A UNIQUE SOLUTION EXISTS UNLESS C THE NODES ARE COLLINEAR. C AN ALTERNATIVE ROUTINE, GRADG, EMPLOYS A GLOBAL METHOD C TO COMPUTE THE PARTIAL DERIVATIVES AT ALL OF THE NODES AT C ONCE. THAT METHOD IS MORE EFFICIENT (WHEN ALL PARTIALS C ARE NEEDED) AND MAY BE MORE ACCURATE, DEPENDING ON THE C DATA. C C INPUT PARAMETERS - N - NUMBER OF NODES IN THE TRIANGULA- C TION. N .GE. 3. C C K - NODE AT WHICH DERIVATIVES ARE C SOUGHT. 1 .LE. K .LE. N. C C X,Y - N-VECTORS CONTAINING THE CARTESIAN C COORDINATES OF THE NODES. C C Z - N-VECTOR CONTAINING THE DATA VALUES C ASSOCIATED WITH THE NODES. C C IADJ - SET OF ADJACENCY LISTS. C C IEND - POINTERS TO THE ENDS OF ADJACENCY C LISTS FOR EACH NODE. C C IADJ AND IEND MAY BE CREATED BY SUBROUTINE TRMESH. C C INPUT PARAMETERS ARE NOT ALTERED BY THIS ROUTINE. C C OUTPUT PARAMETERS - DX,DY - ESTIMATED PARTIAL DERIVATIVES C AT NODE K UNLESS IER .LT. 0. C C IER - ERROR INDICATOR C IER .GT. 0 IF NO ERRORS WERE C ENCOUNTERED. IER C CONTAINS THE NUMBER C OF NODES (INCLUDING C K) USED IN THE FIT. C IER = 3, 4, OR 5 IM- C PLIES A LINEAR FIT. C IER = -1 IF N OR K IS OUT OF C RANGE. C IER = -2 IF ALL NODES ARE C COLLINEAR. C C MODULES REFERENCED BY GRADL - GETNP, SETUP, GIVENS, C ROTATE C C INTRINSIC FUNCTIONS CALLED BY GRADL - MIN0, FLOAT, SQRT, C AMIN1, ABS C C*********************************************************** C INTEGER NN, KK, LMN, LMX, LMIN, LMAX, LM1, LNP, . NPTS(30), IERR, NP, I, J, IM1, JP1, IP1, L REAL SUM, DS, R, RS, RTOL, AVSQ, AV, XK, YK, ZK, . A(6,6), C, S, DMIN, DTOL, SF DATA LMN/10/ DATA LMX/30/, RTOL/1.E-5/, DTOL/.01/, SF/1./ C C LOCAL PARAMETERS - C C NN,KK = LOCAL COPIES OF N AND K C LMN,LMX = MINIMUM AND MAXIMUM VALUES OF LNP FOR N C SUFFICIENTLY LARGE. IN MOST CASES LMN-1 C NODES ARE USED IN THE FIT. 4 .LE. LMN .LE. C LMX. C LMIN,LMAX = MIN(LMN,N), MIN(LMX,N) C LM1 = LMIN-1 OR LNP-1 C LNP = LENGTH OF NPTS C NPTS = ARRAY CONTAINING THE INDICES OF A SEQUENCE OF C NODES ORDERED BY DISTANCE FROM K. NPTS(1)=K C AND THE FIRST LNP-1 ELEMENTS OF NPTS ARE C USED IN THE LEAST SQUARES FIT. UNLESS LNP C EXCEEDS LMAX, NPTS(LNP) DETERMINES R. C IERR = ERROR FLAG FOR CALLS TO GETNP (NOT CHECKED) C NP = ELEMENT OF NPTS TO BE ADDED TO THE SYSTEM C I,J = DO-LOOP INDICES C IM1,JP1 = I-1, J+1 C IP1 = I+1 C L = NUMBER OF COLUMNS OF A**T TO WHICH A ROTATION C IS APPLIED C SUM = SUM OF SQUARED EUCLIDEAN DISTANCES BETWEEN C NODE K AND THE NODES USED IN THE LEAST C SQUARES FIT C DS = SQUARED DISTANCE BETWEEN NODE K AND AN ELE- C MENT OF NPTS C R = DISTANCE BETWEEN NODE K AND NPTS(LNP) OR SOME C POINT FURTHER FROM K THAN NPTS(LMAX) IF C NPTS(LMAX) IS USED IN THE FIT. R IS A C RADIUS OF INFLUENCE WHICH ENTERS INTO THE C WEIGHTS (SEE SUBROUTINE SETUP). C RS = R*R C RTOL = TOLERANCE FOR DETERMINING R. IF THE RELATIVE C CHANGE IN DS BETWEEN TWO ELEMENTS OF NPTS IS C NOT GREATER THAN RTOL THEY ARE TREATED AS C BEING THE SAME DISTANCE FROM NODE K C AVSQ = AV*AV C AV = ROOT-MEAN-SQUARE DISTANCE BETWEEN K AND THE C NODES (OTHER THAN K) IN THE LEAST SQUARES C FIT. THE FIRST 3 COLUMNS OF THE SYSTEM ARE C SCALED BY 1/AVSQ, THE NEXT 2 BY 1/AV. C XK,YK,ZK = COORDINATES AND DATA VALUE ASSOCIATED WITH K C A = TRANSPOSE OF THE AUGMENTED REGRESSION MATRIX C C,S = COMPONENTS OF THE PLANE ROTATION DETERMINED C BY SUBROUTINE GIVENS C DMIN = MINIMUM OF THE MAGNITUDES OF THE DIAGONAL C ELEMENTS OF THE REGRESSION MATRIX AFTER C ZEROS ARE INTRODUCED BELOW THE DIAGONAL C DTOL = TOLERANCE FOR DETECTING AN ILL-CONDITIONED C SYSTEM. THE SYSTEM IS ACCEPTED WHEN DMIN C .GE. DTOL C SF = MARQUARDT STABILIZATION FACTOR USED TO DAMP C OUT THE FIRST 3 SOLUTION COMPONENTS (SECOND C PARTIALS OF THE QUADRATIC) WHEN THE SYSTEM C IS ILL-CONDITIONED. AS SF INCREASES, THE C FITTING FUNCTION APPROACHES A LINEAR C NN = N KK = K C C CHECK FOR ERRORS AND INITIALIZE LMIN, LMAX C IF (NN .LT. 3 .OR. KK .LT. 1 .OR. KK .GT. NN) . GO TO 16 LMIN = MIN0(LMN,NN) LMAX = MIN0(LMX,NN) C C COMPUTE NPTS, LNP, AVSQ, AV, AND R. C SET NPTS TO THE CLOSEST LMIN-1 NODES TO K. C SUM = 0. NPTS(1) = KK LM1 = LMIN - 1 DO 1 LNP = 2,LM1 CALL GETNP (X,Y,IADJ,IEND,LNP, NPTS, DS,IERR) 1 SUM = SUM + DS C C ADD ADDITIONAL NODES TO NPTS UNTIL THE RELATIVE INCREASE C IN DS IS AT LEAST RTOL. C DO 2 LNP = LMIN,LMAX CALL GETNP (X,Y,IADJ,IEND,LNP, NPTS, RS,IERR) IF ((RS-DS)/DS .LE. RTOL) GO TO 2 IF (LNP .GT. 6) GO TO 3 2 SUM = SUM + RS C C USE ALL LMAX NODES IN THE LEAST SQUARES FIT. RS IS C ARBITRARILY INCREASED BY 10 PER CENT. C RS = 1.1*RS LNP = LMAX + 1 C C THERE ARE LNP-2 EQUATIONS CORRESPONDING TO NODES NPTS(2), C ...,NPTS(LNP-1). C 3 AVSQ = SUM/FLOAT(LNP-2) AV = SQRT(AVSQ) R = SQRT(RS) XK = X(KK) YK = Y(KK) ZK = Z(KK) IF (LNP .LT. 7) GO TO 12 C C SET UP THE FIRST 5 EQUATIONS OF THE AUGMENTED REGRESSION C MATRIX (TRANSPOSED) AS THE COLUMNS OF A, AND ZERO OUT C THE LOWER TRIANGLE (UPPER TRIANGLE OF A) WITH GIVENS C ROTATIONS C DO 5 I = 1,5 NP = NPTS(I+1) CALL SETUP (XK,YK,ZK,X(NP),Y(NP),Z(NP),AV,AVSQ, . R, A(1,I)) IF (I .EQ. 1) GO TO 5 IM1 = I - 1 DO 4 J = 1,IM1 JP1 = J + 1 L = 6 - J CALL GIVENS (A(J,J),A(J,I),C,S) 4 CALL ROTATE (L,C,S,A(JP1,J),A(JP1,I)) 5 CONTINUE C C ADD THE ADDITIONAL EQUATIONS TO THE SYSTEM USING C THE LAST COLUMN OF A -- I .LE. LNP. C I = 7 6 IF (I .EQ. LNP) GO TO 8 NP = NPTS(I) CALL SETUP (XK,YK,ZK,X(NP),Y(NP),Z(NP),AV,AVSQ, . R, A(1,6)) DO 7 J = 1,5 JP1 = J + 1 L = 6 - J CALL GIVENS (A(J,J),A(J,6),C,S) 7 CALL ROTATE (L,C,S,A(JP1,J),A(JP1,6)) I = I + 1 GO TO 6 C C TEST THE SYSTEM FOR ILL-CONDITIONING C 8 DMIN = AMIN1( ABS(A(1,1)),ABS(A(2,2)),ABS(A(3,3)), . ABS(A(4,4)),ABS(A(5,5)) ) IF (DMIN .GE. DTOL) GO TO 15 IF (LNP .GT. LMAX) GO TO 9 C C ADD ANOTHER NODE TO THE SYSTEM AND INCREASE R -- C I .EQ. LNP C LNP = LNP + 1 IF (LNP .LE. LMAX) CALL GETNP (X,Y,IADJ,IEND,LNP, . NPTS, RS,IERR) R = SQRT(1.1*RS) GO TO 6 C C STABILIZE THE SYSTEM BY DAMPING SECOND PARTIALS --ADD C MULTIPLES OF THE FIRST THREE UNIT VECTORS TO THE FIRST C THREE EQUATIONS. C 9 DO 11 I = 1,3 A(I,6) = SF IP1 = I + 1 DO 10 J = IP1,6 10 A(J,6) = 0. DO 11 J = I,5 JP1 = J + 1 L = 6 - J CALL GIVENS (A(J,J),A(J,6),C,S) 11 CALL ROTATE (L,C,S,A(JP1,J),A(JP1,6)) GO TO 14 C C 4 .LE. LNP .LE. 6 (2, 3, OR 4 EQUATIONS) -- FIT A PLANE TO C THE DATA USING THE LAST 3 COLUMNS OF A. C 12 NP = NPTS(2) CALL SETUP (XK,YK,ZK,X(NP),Y(NP),Z(NP),AV,AVSQ, . R, A(1,4)) NP = NPTS(3) CALL SETUP (XK,YK,ZK,X(NP),Y(NP),Z(NP),AV,AVSQ, . R, A(1,5)) CALL GIVENS (A(4,4),A(4,5),C,S) CALL ROTATE (2,C,S,A(5,4),A(5,5)) IF (LNP .EQ. 4) GO TO 14 C LM1 = LNP - 1 DO 13 I = 4,LM1 NP = NPTS(I) CALL SETUP (XK,YK,ZK,X(NP),Y(NP),Z(NP),AV,AVSQ, . R, A(1,6)) CALL GIVENS (A(4,4),A(4,6),C,S) CALL ROTATE (2,C,S,A(5,4),A(5,6)) CALL GIVENS (A(5,5),A(5,6),C,S) 13 CALL ROTATE (1,C,S,A(6,5),A(6,6)) C C TEST THE LINEAR FIT FOR ILL-CONDITIONING C 14 DMIN = AMIN1( ABS(A(4,4)),ABS(A(5,5)) ) IF (DMIN .LT. DTOL) GO TO 17 C C SOLVE THE 2 BY 2 TRIANGULAR SYSTEM FOR THE DERIVATIVES C 15 DY = A(6,5)/A(5,5) DX = (A(6,4) - A(5,4)*DY)/A(4,4)/AV DY = DY/AV IER = LNP - 1 RETURN C C N OR K IS OUT OF RANGE C 16 IER = -1 RETURN C C NO UNIQUE SOLUTION DUE TO COLLINEAR NODES C 17 IER = -2 RETURN END SUBROUTINE INTRC0 (N,PX,PY,X,Y,Z,IADJ,IEND, IST, PZ, . IER) INTEGER N, IADJ(1), IEND(N), IST, IER REAL PX, PY, X(N), Y(N), Z(N), PZ C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C GIVEN A TRIANGULATION OF A SET OF POINTS IN THE PLANE, C THIS ROUTINE COMPUTES THE VALUE AT (PX,PY) OF A PIECEWISE C LINEAR SURFACE WHICH INTERPOLATES DATA VALUES AT THE C VERTICES OF THE TRIANGLES. THE SURFACE IS EXTENDED IN A C CONTINUOUS FASHION BEYOND THE BOUNDARY OF THE TRIANGULAR C MESH, ALLOWING EXTRAPOLATION. INTRC0 IS PART OF AN C INTERPOLATION PACKAGE WHICH PROVIDES ROUTINES TO GENERATE, C UPDATE, AND PLOT THE MESH. C C INPUT PARAMETERS - N - NUMBER OF NODES IN THE MESH. C N .GE. 3. C C PX,PY - POINT AT WHICH THE INTERPOLATED C VALUE IS DESIRED. C C X,Y - VECTORS OF COORDINATES OF THE C NODES IN THE MESH. C C Z - VECTOR OF DATA VALUES AT THE C NODES. C C IADJ - SET OF ADJACENCY LISTS OF NODES C IN THE MESH. C C IEND - POINTERS TO THE ENDS OF C ADJACENCY LISTS IN IADJ FOR C EACH NODE IN THE MESH. C C IST - INDEX OF THE STARTING NODE IN C THE SEARCH FOR A TRIANGLE CON- C TAINING (PX,PY). 1 .LE. IST C .LE. N. THE OUTPUT VALUE OF C IST FROM A PREVIOUS CALL MAY C BE A GOOD CHOICE. C C IADJ AND IEND MAY BE CREATED BY TRMESH. C C INPUT PARAMETERS OTHER THAN IST ARE NOT ALTERED BY THIS C ROUTINE. C C OUTPUT PARAMETERS - IST - INDEX OF ONE OF THE VERTICES OF C THE TRIANGLE CONTAINING (PX,PY) C UNLESS IER .LT. 0. C C PZ - VALUE OF THE INTERPOLATORY C SURFACE AT (PX,PY) OR ZERO C IF IER .LT. 0. C C IER - ERROR INDICATOR C IER = 0 IF NO ERRORS WERE C ENCOUNTERED. C IER = 1 IF NO ERRORS WERE EN- C COUNTERED AND EXTRAPO- C LATION WAS PERFORMED. C IER = -1 IF N OR IST IS OUT OF C RANGE. C IER = -2 IF THE NODES ARE COL- C LINEAR. C C MODULES REFERENCED BY INTRC0 - TRFIND, COORDS C C*********************************************************** C INTEGER I1, I2, I3, N1, N2, INDX REAL XP, YP, R(3), X1, Y1, X2, Y2, DP C C LOCAL PARAMETERS - C C I1,I2,I3 = VERTEX INDICES RETURNED BY TRFIND C N1,N2 = ENDPOINTS OF THE CLOSEST BOUNDARY EDGE TO P C WHEN P IS OUTSIDE OF THE MESH BOUNDARY C INDX = IADJ INDEX OF N1 AS A NEIGHBOR OF N2 C XP,YP = LOCAL COPIES OF THE COORDINATES OF P=(PX,PY) C R = BARYCENTRIC COORDINATES C X1,Y1 = X,Y COORDINATES OF N1 C X2,Y2 = X,Y COORDINATES OF N2 C DP = INNER PRODUCT OF N1-N2 AND P-N2 C IF (N .LT. 3 .OR. IST .LT. 1 .OR. IST .GT. N) . GO TO 5 XP = PX YP = PY C C FIND A TRIANGLE CONTAINING P IF P IS WITHIN THE MESH C BOUNDARY C CALL TRFIND(IST,XP,YP,X,Y,IADJ,IEND, I1,I2,I3) IF (I1 .EQ. 0) GO TO 6 IST = I1 IF (I3 .EQ. 0) GO TO 1 C C COMPUTE BARYCENTRIC COORDINATES C CALL COORDS(XP,YP,X(I1),X(I2),X(I3),Y(I1),Y(I2), . Y(I3), R,IER) IF (IER .NE. 0) GO TO 6 PZ = R(1)*Z(I1) + R(2)*Z(I2) + R(3)*Z(I3) RETURN C C P IS OUTSIDE OF THE MESH BOUNDARY. EXTRAPOLATE TO P BY C EXTENDING THE INTERPOLATORY SURFACE AS A CONSTANT C BEYOND THE BOUNDARY. THUS PZ IS THE SURFACE FUNCTION C VALUE AT Q WHERE Q IS THE CLOSEST BOUNDARY POINT TO P. C C DETERMINE Q BY TRAVERSING THE BOUNDARY STARTING FROM THE C RIGHTMOST VISIBLE NODE I1. C 1 N2 = I1 C C SET N1 TO THE LAST NONZERO NEIGHBOR OF N2 AND COMPUTE DP C 2 INDX = IEND(N2) - 1 N1 = IADJ(INDX) X1 = X(N1) Y1 = Y(N1) X2 = X(N2) Y2 = Y(N2) DP = (X1-X2)*(XP-X2) + (Y1-Y2)*(YP-Y2) IF (DP .LE. 0.) GO TO 3 IF ((XP-X1)*(X2-X1) + (YP-Y1)*(Y2-Y1) .GT. 0.) GO TO 4 N2 = N1 GO TO 2 C C N2 IS THE CLOSEST BOUNDARY POINT TO P C 3 PZ = Z(N2) IER = 1 RETURN C C THE CLOSEST BOUNDARY POINT TO P LIES ON N2-N1. COMPUTE C ITS COORDINATES WITH RESPECT TO N2-N1. C 4 R(1) = DP/( (X2-X1)**2 + (Y2-Y1)**2 ) R(2) = 1. - R(1) PZ = R(1)*Z(N1) + R(2)*Z(N2) IER = 1 RETURN C C N .LT. 3 OR IST IS OUT OF RANGE C 5 PZ = 0. IER = -1 RETURN C C NODES ARE COLLINEAR C 6 PZ = 0. IER = -2 RETURN END SUBROUTINE INTRC1 (N,PX,PY,X,Y,Z,IADJ,IEND,IFLAG, . ZXZY, IST, PZ,IER) INTEGER N, IADJ(1), IEND(N), IFLAG, IST, IER REAL PX, PY, X(N), Y(N), Z(N), ZXZY(2,N), PZ C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C GIVEN A TRIANGULATION OF A SET OF POINTS IN THE PLANE, C THIS ROUTINE DETERMINES A PIECEWISE CUBIC FUNCTION F(X,Y) C WHICH INTERPOLATES A SET OF DATA VALUES AND PARTIAL C DERIVATIVES AT THE VERTICES. F HAS CONTINUOUS FIRST C DERIVATIVES OVER THE MESH AND EXTENDS BEYOND THE MESH C BOUNDARY ALLOWING EXTRAPOLATION. INTERPOLATION IS EXACT C FOR QUADRATIC DATA. THE VALUE OF F AT (PX,PY) IS C RETURNED. INTRC1 IS PART OF AN INTERPOLATION PACKAGE C WHICH PROVIDES ROUTINES TO GENERATE, UPDATE AND PLOT THE C MESH. C C INPUT PARAMETERS - N - NUMBER OF NODES IN THE MESH. C N .GE. 3. C C PX,PY - COORDINATES OF A POINT AT WHICH C F IS TO BE EVALUATED. C C X,Y - VECTORS OF COORDINATES OF THE C NODES IN THE MESH. C C Z - VECTOR OF DATA VALUES AT THE C NODES. C C IADJ - SET OF ADJACENCY LISTS OF NODES C IN THE MESH. C C IEND - POINTERS TO THE ENDS OF C ADJACENCY LISTS IN IADJ FOR C EACH NODE IN THE MESH. C C IFLAG - OPTION INDICATOR C IFLAG = 0 IF INTRC1 IS TO C PROVIDE DERIVATIVE C ESTIMATES (FROM C GRADL). C IFLAG = 1 IF DERIVATIVES ARE C USER PROVIDED. C C ZXZY - 2 BY N ARRAY WHOSE COLUMNS C CONTAIN ESTIMATED PARTIAL DER- C IVATIVES AT THE NODES (X PAR- C TIALS IN THE FIRST ROW) IF C IFLAG = 1, NOT USED IF IFLAG C = 0. C C IST - INDEX OF THE STARTING NODE IN C THE SEARCH FOR A TRIANGLE CON- C TAINING (PX,PY). 1 .LE. IST C .LE. N. THE OUTPUT VALUE OF C IST FROM A PREVIOUS CALL MAY C BE A GOOD CHOICE. C C IADJ AND IEND MAY BE CREATED BY TRMESH AND DERIVATIVE C ESTIMATES MAY BE COMPUTED BY GRADL OR GRADG. C C INPUT PARAMETERS OTHER THAN IST ARE NOT ALTERED BY THIS C ROUTINE. C C OUTPUT PARAMETERS - IST - INDEX OF ONE OF THE VERTICES OF C THE TRIANGLE CONTAINING (PX,PY) C UNLESS IER .LT. 0. C C PZ - VALUE OF F AT (PX,PY), OR 0 IF C IER .LT. 0. C C IER - ERROR INDICATOR C IER = 0 IF NO ERRORS WERE C ENCOUNTERED. C IER = 1 IF NO ERRORS WERE EN- C COUNTERED AND EXTRAPOLA- C TION WAS PERFORMED. C IER = -1 IF N, IFLAG, OR IST IS C OUT OF RANGE. C IER = -2 IF THE NODES ARE COL- C LINEAR. C C MODULES REFERENCED BY INTRC1 - TRFIND, TVAL, C (AND OPTIONALLY) GRADL, GETNP, SETUP, C GIVENS, ROTATE C C*********************************************************** C INTEGER NN, I1, I2, I3, IERR, N1, N2, INDX REAL XP, YP, ZX1, ZY1, ZX2, ZY2, ZX3, ZY3, X1, Y1, . X2, Y2, X3, Y3, Z1, Z2, Z3, DUM, DP, U, V, XQ, . YQ, R1, R2, A1, A2, B1, B2, C1, C2, F1, F2 C C LOCAL PARAMETERS - C C NN = LOCAL COPY OF N C I1,I2,I3 = VERTICES DETERMINED BY TRFIND C IERR = ERROR FLAG FOR CALLS TO GRADL C AND TVAL C N1,N2 = ENDPOINTS OF THE CLOSEST BOUND- C ARY EDGE TO P WHEN P IS OUT- C SIDE OF THE MESH BOUNDARY C INDX = IADJ INDEX OF N1 AS A NEIGHBOR C OF N2 C XP,YP = LOCAL COPIES OF THE COORDINATES C OF P=(PX,PY) C ZX1,ZY1,ZX2,ZY2,ZX3,ZY3 = X AND Y DERIVATIVES AT THE C VERTICES OF A TRIANGLE T WHICH C CONTAINS P OR AT N1 AND N2 C X1,Y1,X2,Y2,X3,Y3 = X,Y COORDINATES OF THE VERTICES C OF T OR OF N1 AND N2 C Z1,Z2,Z3 = DATA VALUES AT THE VERTICES OF T C DUM = DUMMY VARIABLE FOR CALL TO TVAL C DP = INNER PRODUCT OF N1-N2 AND P-N2 C U,V = X,Y COORDINATES OF THE VECTOR C N2-N1 C XQ,YQ = X,Y COORDINATES OF THE CLOSEST C BOUNDARY POINT TO P WHEN P IS C OUTSIDE OF THE MESH BOUNDARY C R1,R2 = BARYCENTRIC COORDINATES OF Q C WITH RESPECT TO THE LINE SEG- C MENT N2-N1 CONTAINING Q C A1,A2,B1,B2,C1,C2 = CARDINAL FUNCTIONS FOR EVALUAT- C ING THE INTERPOLATORY SURFACE C AT Q C F1,F2 = CUBIC FACTORS USED TO COMPUTE C THE CARDINAL FUNCTIONS C NN = N PZ = 0. IF (NN .LT. 3 .OR. IFLAG .LT. 0 .OR. IFLAG .GT. 1 . .OR. IST .LT. 1 .OR. IST .GT. NN) GO TO 11 XP = PX YP = PY C C FIND A TRIANGLE CONTAINING P IF P IS WITHIN THE MESH C BOUNDARY C CALL TRFIND(IST,XP,YP,X,Y,IADJ,IEND, I1,I2,I3)