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 IF (ST0) print 1111,X0/1e1,X(N1)/1e1,Y0,Y(N1) IF (.NOT. ST0) print 1111,X(N1)/1e1,X0/1e1,Y(N1),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 1111 format(' oplot,[',f6.4,',',f6.4,'],[',f6.4,',',f6.4,']') 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) IF (I1 .EQ. 0) GO TO 12 IST = I1 IF (I3 .EQ. 0) GO TO 3 IF (IFLAG .NE. 1) GO TO 1 C C DERIVATIVES ARE USER PROVIDED C ZX1 = ZXZY(1,I1) ZX2 = ZXZY(1,I2) ZX3 = ZXZY(1,I3) ZY1 = ZXZY(2,I1) ZY2 = ZXZY(2,I2) ZY3 = ZXZY(2,I3) GO TO 2 C C COMPUTE DERIVATIVE ESTIMATES AT THE VERTICES C 1 CALL GRADL(NN,I1,X,Y,Z,IADJ,IEND, ZX1,ZY1,IERR) CALL GRADL(NN,I2,X,Y,Z,IADJ,IEND, ZX2,ZY2,IERR) CALL GRADL(NN,I3,X,Y,Z,IADJ,IEND, ZX3,ZY3,IERR) C C SET LOCAL PARAMETERS FOR CALL TO TVAL C 2 X1 = X(I1) Y1 = Y(I1) X2 = X(I2) Y2 = Y(I2) X3 = X(I3) Y3 = Y(I3) Z1 = Z(I1) Z2 = Z(I2) Z3 = Z(I3) CALL TVAL(XP,YP,X1,X2,X3,Y1,Y2,Y3,Z1,Z2,Z3,ZX1,ZX2, . ZX3,ZY1,ZY2,ZY3,0, PZ,DUM,DUM,IERR) IF (IERR .NE. 0) GO TO 12 IER = 0 RETURN C C P IS OUTSIDE OF THE MESH BOUNDARY. EXTRAPOLATE TO P BY C PASSING A LINEAR FUNCTION OF ONE VARIABLE THROUGH THE C VALUE AND DIRECTIONAL DERIVATIVE (IN THE DIRECTION C P-Q) OF THE INTERPOLATORY SURFACE (TVAL) AT Q WHERE C Q IS THE CLOSEST BOUNDARY POINT TO P. C C DETERMINE Q BY TRAVERSING THE BOUNDARY STARTING FROM C THE RIGHTMOST VISIBLE NODE I1. C 3 N2 = I1 C C SET N1 TO THE LAST NONZERO NEIGHBOR OF N2 AND COMPUTE DP C 4 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 5 IF ((XP-X1)*(X2-X1) + (YP-Y1)*(Y2-Y1) .GT. 0.) GO TO 8 N2 = N1 GO TO 4 C C N2 IS THE CLOSEST BOUNDARY POINT TO P. COMPUTE PARTIAL C DERIVATIVES AT N2. C 5 IF (IFLAG .NE. 1) GO TO 6 ZX2 = ZXZY(1,N2) ZY2 = ZXZY(2,N2) GO TO 7 6 CALL GRADL(NN,N2,X,Y,Z,IADJ,IEND, ZX2,ZY2,IERR) C C COMPUTE EXTRAPOLATED VALUE AT P C 7 PZ = Z(N2) + ZX2*(XP-X2) + ZY2*(YP-Y2) IER = 1 RETURN C C THE CLOSEST BOUNDARY POINT Q LIES ON N2-N1. COMPUTE C PARTIALS AT N1 AND N2. C 8 IF (IFLAG .NE. 1) GO TO 9 ZX1 = ZXZY(1,N1) ZY1 = ZXZY(2,N1) ZX2 = ZXZY(1,N2) ZY2 = ZXZY(2,N2) GO TO 10 9 CALL GRADL(NN,N1,X,Y,Z,IADJ,IEND, ZX1,ZY1,IERR) CALL GRADL(NN,N2,X,Y,Z,IADJ,IEND, ZX2,ZY2,IERR) C C COMPUTE Q, ITS BARYCENTRIC COORDINATES, AND THE CARDINAL C FUNCTIONS FOR EXTRAPOLATION C 10 U = X2-X1 V = Y2-Y1 R1 = DP/(U**2 + V**2) R2 = 1. - R1 XQ = R1*X1 + R2*X2 YQ = R1*Y1 + R2*Y2 F1 = R1*R1*R2 F2 = R1*R2*R2 A1 = R1 + (F1-F2) A2 = R2 - (F1-F2) B1 = U*F1 B2 = -U*F2 C1 = V*F1 C2 = -V*F2 C C COMPUTE THE VALUE OF THE INTERPOLATORY SURFACE (TVAL) C AT Q C PZ = A1*Z(N1) + A2*Z(N2) + B1*ZX1 + B2*ZX2 + . C1*ZY1 + C2*ZY2 C C COMPUTE THE EXTRAPOLATED VALUE AT P C PZ = PZ + (R1*ZX1 + R2*ZX2)*(XP-XQ) + . (R1*ZY1 + R2*ZY2)*(YP-YQ) IER = 1 RETURN C C N, IFLAG, OR IST OUT OF RANGE C 11 IER = -1 RETURN C C NODES ARE COLLINEAR C 12 IER = -2 RETURN END SUBROUTINE ROTATE (N,C,S, X,Y ) INTEGER N REAL C, S, X(N), Y(N) C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C ( C S) C THIS ROUTINE APPLIES THE GIVENS ROTATION ( ) TO THE C (-S C) C (X(1) ... X(N)) C 2 BY N MATRIX ( ). THIS ROUTINE WAS TAKEN C (Y(1) ... Y(N)) C LINPACK. C C INPUT PARAMETERS - N - NUMBER OF COLUMNS TO BE ROTATED. C C C,S - ELEMENTS OF THE GIVENS ROTATION. C THESE MAY BE DETERMINED BY C SUBROUTINE GIVENS. C C X,Y - VECTORS OF LENGTH .GE. N C CONTAINING THE 2-VECTORS TO BE C ROTATED. C C THE PARAMETERS N, C, AND S ARE NOT ALTERED BY THIS C ROUTINE. C C OUTPUT PARAMETERS - X,Y - ROTATED VECTORS C C MODULES REFERENCED BY ROTATE - NONE C C*********************************************************** C INTEGER I REAL XI, YI C C LOCAL PARAMETERS - C C I = DO-LOOP INDEX C XI,YI = X(I), Y(I) C IF (N .LE. 0 .OR. (C .EQ. 1. .AND. S .EQ. 0.)) RETURN DO 1 I = 1,N XI = X(I) YI = Y(I) X(I) = C*XI + S*YI Y(I) = -S*XI + C*YI 1 CONTINUE RETURN END SUBROUTINE SETUP (XK,YK,ZK,XI,YI,ZI,S1,S2,R, ROW) REAL XK, YK, ZK, XI, YI, ZI, S1, S2, R, ROW(6) C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C THIS ROUTINE SETS UP THE I-TH ROW OF AN AUGMENTED RE- C GRESSION MATRIX FOR A WEIGHTED LEAST-SQUARES FIT OF A C QUADRATIC FUNCTION Q(X,Y) TO A SET OF DATA VALUES Z WHERE C Q(XK,YK) = ZK. THE FIRST 3 COLUMNS (QUADRATIC TERMS) ARE C SCALED BY 1/S2 AND THE FOURTH AND FIFTH COLUMNS (LINEAR C TERMS) ARE SCALED BY 1/S1. THE WEIGHT IS (R-D)/(R*D) IF C R .GT. D AND 0 IF R .LE. D, WHERE D IS THE DISTANCE C BETWEEN NODES I AND K. C C INPUT PARAMETERS - XK,YK,ZK - COORDINATES AND DATA VALUE C AT NODE K -- INTERPOLATED C BY Q. C C XI,YI,ZI - COORDINATES AND DATA VALUE C AT NODE I. C C S1,S2 - INVERSE SCALE FACTORS. C C R - RADIUS OF INFLUENCE ABOUT C NODE K DEFINING THE WEIGHT. C C ROW - VECTOR OF LENGTH 6. C C INPUT PARAMETERS ARE NOT ALTERED BY THIS ROUTINE. C C OUTPUT PARAMETER - ROW - VECTOR CONTAINING A ROW OF THE C AUGMENTED REGRESSION MATRIX. C C MODULES REFERENCED BY SETUP - NONE C C INTRINSIC FUNCTION CALLED BY SETUP - SQRT C C*********************************************************** C INTEGER I REAL DX, DY, DXSQ, DYSQ, D, W, W1, W2 C C LOCAL PARAMETERS - C C I = DO-LOOP INDEX C DX = XI - XK C DY = YI - YK C DXSQ = DX*DX C DYSQ = DY*DY C D = DISTANCE BETWEEN NODES K AND I C W = WEIGHT ASSOCIATED WITH THE ROW C W1 = W/S1 C W2 = W/S2 C DX = XI - XK DY = YI - YK DXSQ = DX*DX DYSQ = DY*DY D = SQRT(DXSQ + DYSQ) IF (D .LE. 0. .OR. D .GE. R) GO TO 1 W = (R-D)/R/D W1 = W/S1 W2 = W/S2 ROW(1) = DXSQ*W2 ROW(2) = DX*DY*W2 ROW(3) = DYSQ*W2 ROW(4) = DX*W1 ROW(5) = DY*W1 ROW(6) = (ZI - ZK)*W RETURN C C NODES K AND I COINCIDE OR NODE I IS OUTSIDE OF THE RADIUS C OF INFLUENCE. SET ROW TO THE ZERO VECTOR. C 1 DO 2 I = 1,6 2 ROW(I) = 0. RETURN END FUNCTION TRVOL (X1,X2,X3,Y1,Y2,Y3,Z1,Z2,Z3) REAL X1, X2, X3, Y1, Y2, Y3, Z1, Z2, Z3 C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C THIS FUNCTION COMPUTES THE INTEGRAL OVER A TRIANGLE OF C THE PLANAR SURFACE WHICH INTERPOLATES DATA VALUES AT THE C VERTICES. C C INPUT PARAMETERS - X1,X2,X3 - X COORDINATES OF THE C VERTICES OF THE TRIANGLE IN C COUNTERCLOCKWISE ORDER. C C Y1,Y2,Y3 - Y COORDINATES OF THE C VERTICES OF THE TRIANGLE IN C THE SAME ORDER AS THE X C COORDINATES. C C Z1,Z2,Z3 - DATA VALUES AT THE VERTICES C IN THE SAME ORDER AS THE C COORDINATES. C C INPUT PARAMETERS ARE NOT ALTERED BY THIS FUNCTION. C C OUTPUT PARAMETER - TRVOL - VOLUME UNDER THE INTERPOLA- C TORY SURFACE ABOVE THE C TRIANGLE OR ZERO IF THE C COORDINATES ARE INCORRECTLY C ORDERED OR COLLINEAR. C C MODULES REFERENCED BY TRVOL - NONE C C*********************************************************** C T1 = X2*Y3 - X3*Y2 T2 = X3*Y1 - X1*Y3 T3 = X1*Y2 - X2*Y1 AREA = T1 + T2 + T3 IF (AREA .LT. 0.) AREA = 0. C C AREA IS TWICE THE AREA OF THE TRIANGLE. TRVOL IS THE MEAN C OF THE DATA VALUES TIMES THE AREA OF THE TRIANGLE. C TRVOL = (Z1 + Z2 + Z3)*AREA/6. RETURN END SUBROUTINE TVAL (X,Y,X1,X2,X3,Y1,Y2,Y3,Z1,Z2,Z3,ZX1, . ZX2,ZX3,ZY1,ZY2,ZY3,IFLAG, W,WX,WY, . IER) INTEGER IFLAG, IER REAL X, Y, X1, X2, X3, Y1, Y2, Y3, Z1, Z2, Z3, . ZX1, ZX2, ZX3, ZY1, ZY2, ZY3, W, WX, WY C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C GIVEN FUNCTION VALUES AND FIRST PARTIAL DERIVATIVES AT C THE THREE VERTICES OF A TRIANGLE, THIS ROUTINE DETERMINES C A FUNCTION W WHICH AGREES WITH THE GIVEN DATA, RETURNING C THE VALUE AND (OPTIONALLY) FIRST PARTIAL DERIVATIVES OF W C AT A POINT (X,Y) IN THE TRIANGLE. THE INTERPOLATION C METHOD IS EXACT FOR QUADRATIC POLYNOMIAL DATA. THE C TRIANGLE IS PARTITIONED INTO THREE SUBTRIANGLES WITH C EQUAL AREAS. W IS CUBIC IN EACH SUBTRIANGLE AND ALONG C THE EDGES, BUT HAS ONLY ONE CONTINUOUS DERIVATIVE ACROSS C EDGES. THE NORMAL DERIVATIVE OF W VARIES LINEARLY ALONG C EACH OUTER EDGE. THE VALUES AND PARTIAL DERIVATIVES OF W C ALONG A TRIANGLE EDGE DEPEND ONLY ON THE DATA VALUES AT C THE ENDPOINTS OF THE EDGE. THUS THE METHOD YIELDS C-1 C CONTINUITY WHEN USED TO INTERPOLATE OVER A TRIANGULAR C GRID. THIS ALGORITHM IS DUE TO C. L. LAWSON. C C INPUT PARAMETERS - X,Y - COORDINATES OF A POINT AT WHICH C W IS TO BE EVALUATED. C C X1,X2,X3,Y1,Y2,Y3 - COORDINATES OF THE VERTICES OF C A TRIANGLE CONTAINING (X,Y). C C Z1,Z2,Z3 - FUNCTION VALUES AT THE VERTICES C TO BE INTERPOLATED. C C ZX1,ZX2,ZX3 - X-DERIVATIVE VALUES AT THE C VERTICES. C C ZY1,ZY2,ZY3 - Y-DERIVATIVE VALUES AT THE C VERTICES. C C IFLAG - OPTION INDICATOR C IFLAG = 0 IF ONLY W IS TO BE C COMPUTED. C IFLAG = 1 IF W, WX, AND WY ARE C TO BE RETURNED. C C INPUT PARAMETERS ARE NOT ALTERED BY THIS ROUTINE. C C OUTPUT PARAMETERS - W - ESTIMATED VALUE OF THE INTERP- C OLATORY FUNCTION AT (X,Y) IF C IER = 0. OTHERWISE W = 0. C C WX,WY - PARTIAL DERIVATIVES OF W AT C (X,Y) IF IER = 0 AND IFLAG = 1, C UNCHANGED IF IFLAG .NE. 1, ZERO C IF IER .NE. 0 AND IFLAG = 1. 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 TVAL - NONE C C INTRINSIC FUNCTION CALLED BY TVAL - AMIN1 C C*********************************************************** C INTEGER I, IP1, IP2, IP3 REAL U(3), V(3), SL(3), AREA, XP, YP, R(3), RX(3), . RY(3), PHI(3), PHIX(3), PHIY(3), RMIN, C1, C2, . RO(3), ROX(3), ROY(3), F(3), G(3), GX(3), . GY(3), P(3), PX(3), PY(3), Q(3), QX(3), QY(3), . A(3), AX(3), AY(3), B(3), BX(3), BY(3), C(3), . CX(3), CY(3) C C LOCAL PARAMETERS - C C I = DO-LOOP INDEX C IP1,IP2,IP3 = PERMUTED INDICES FOR COMPUTING RO, ROX, C AND ROY C U(K) = X-COMPONENT OF THE VECTOR REPRESENTING C THE SIDE OPPOSITE VERTEX K C V(K) = Y-COMPONENT OF THE VECTOR REPRESENTING C THE SIDE OPPOSITE VERTEX K C SL(K) = SQUARE OF THE LENGTH OF THE SIDE C OPPOSITE VERTEX K C AREA = TWICE THE AREA OF THE TRIANGLE C XP,YP = X-X1, Y-Y1 C R(K) = K-TH BARYCENTRIC COORDINATE C RX(K),RY(K) = X,Y PARTIAL DERIVATIVES OF R(K) C PHI(K) R(K-1)*R(K+1) -- QUADRATIC C PHIX(K),PHIY(K) = X,Y PARTIALS OF PHI(K) C RMIN = MIN(R1,R2,R3) C C1,C2 = FACTORS FOR COMPUTING RO C RO(K) = FACTORS FOR COMPUTING G -- CUBIC C CORRECTION TERMS C ROX(K),ROY(K) = X,Y PARTIALS OF RO(K) C F(K) = FACTORS FOR COMPUTING G, GX, AND GY -- C CONSTANT C G(K) = FACTORS FOR COMPUTING THE CARDINAL C FUNCTIONS -- CUBIC C GX(K),GY(K) = X,Y PARTIALS OF G(K) C P(K) = G(K) + PHI(K) C PX(K),PY(K) = X,Y PARTIALS OF P(K) C Q(K) = G(K) - PHI(K) C QX(K),QY(K) = X,Y PARTIALS OF Q(K) C A(K) = CARDINAL FUNCTION WHOSE COEFFICIENT IS C Z(K) C AX(K),AY(K) = X,Y PARTIALS OF A(K) -- CARDINAL C FUNCTIONS FOR WX AND WY C B(K) = TWICE THE CARDINAL FUNCTION WHOSE C COEFFICIENT IS ZX(K) C BX(K),BY(K) = X,Y PARTIALS OF B(K) C C(K) = TWICE THE CARDINAL FUNCTION WHOSE C COEFFICIENT IS ZY(K) C CX(K),CY(K) = X,Y PARTIALS OF C(K) 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 DO 1 I = 1,3 SL(I) = U(I)*U(I) + V(I)*V(I) 1 CONTINUE C C AREA = 3-1 X 3-2 C AREA = U(1)*V(2) - U(2)*V(1) IF (AREA .EQ. 0.) GO TO 9 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 C PHI(1) = R(2)*R(3) PHI(2) = R(3)*R(1) PHI(3) = R(1)*R(2) C RMIN = AMIN1(R(1),R(2),R(3)) IF (RMIN .NE. R(1)) GO TO 3 IP1 = 1 IP2 = 2 IP3 = 3 GO TO 5 3 IF (RMIN .NE. R(2)) GO TO 4 IP1 = 2 IP2 = 3 IP3 = 1 GO TO 5 4 IP1 = 3 IP2 = 1 IP3 = 2 C 5 C1 = RMIN*RMIN/2. C2 = RMIN/3. RO(IP1) = (PHI(IP1) + 5.*C1/3.)*R(IP1) - C1 RO(IP2) = C1*(R(IP3) - C2) RO(IP3) = C1*(R(IP2) - C2) C F(1) = 3.*(SL(2)-SL(3))/SL(1) F(2) = 3.*(SL(3)-SL(1))/SL(2) F(3) = 3.*(SL(1)-SL(2))/SL(3) C G(1) = (R(2)-R(3))*PHI(1) + F(1)*RO(1) - RO(2) + RO(3) G(2) = (R(3)-R(1))*PHI(2) + F(2)*RO(2) - RO(3) + RO(1) G(3) = (R(1)-R(2))*PHI(3) + F(3)*RO(3) - RO(1) + RO(2) C DO 6 I = 1,3 P(I) = G(I) + PHI(I) Q(I) = G(I) - PHI(I) 6 CONTINUE C A(1) = R(1) + G(3) - G(2) A(2) = R(2) + G(1) - G(3) A(3) = R(3) + G(2) - G(1) C B(1) = U(3)*P(3) + U(2)*Q(2) B(2) = U(1)*P(1) + U(3)*Q(3) B(3) = U(2)*P(2) + U(1)*Q(1) C C(1) = V(3)*P(3) + V(2)*Q(2) C(2) = V(1)*P(1) + V(3)*Q(3) C(3) = V(2)*P(2) + V(1)*Q(1) C C W IS A LINEAR COMBINATION OF THE CARDINAL FUNCTIONS C W = A(1)*Z1 + A(2)*Z2 + A(3)*Z3 + (B(1)*ZX1 + B(2)*ZX2 . + B(3)*ZX3 + C(1)*ZY1 + C(2)*ZY2 + C(3)*ZY3)/2. IF (IFLAG .NE. 1) RETURN C C COMPUTE WX AND WY C DO 7 I = 1,3 RX(I) = -V(I)/AREA RY(I) = U(I)/AREA 7 CONTINUE PHIX(1) = R(2)*RX(3) + RX(2)*R(3) PHIY(1) = R(2)*RY(3) + RY(2)*R(3) PHIX(2) = R(3)*RX(1) + RX(3)*R(1) PHIY(2) = R(3)*RY(1) + RY(3)*R(1) PHIX(3) = R(1)*RX(2) + RX(1)*R(2) PHIY(3) = R(1)*RY(2) + RY(1)*R(2) C ROX(IP1) = RX(IP1)*(PHI(IP1) + 5.*C1) + . R(IP1)*(PHIX(IP1) - RX(IP1)) ROY(IP1) = RY(IP1)*(PHI(IP1) + 5.*C1) + . R(IP1)*(PHIY(IP1) - RY(IP1)) ROX(IP2) = RX(IP1)*(PHI(IP2) - C1) + C1*RX(IP3) ROY(IP2) = RY(IP1)*(PHI(IP2) - C1) + C1*RY(IP3) ROX(IP3) = RX(IP1)*(PHI(IP3) - C1) + C1*RX(IP2) ROY(IP3) = RY(IP1)*(PHI(IP3) - C1) + C1*RY(IP2) C GX(1) = (RX(2) - RX(3))*PHI(1) + (R(2) - R(3))*PHIX(1) . + F(1)*ROX(1) - ROX(2) + ROX(3) GY(1) = (RY(2) - RY(3))*PHI(1) + (R(2) - R(3))*PHIY(1) . + F(1)*ROY(1) - ROY(2) + ROY(3) GX(2) = (RX(3) - RX(1))*PHI(2) + (R(3) - R(1))*PHIX(2) . + F(2)*ROX(2) - ROX(3) + ROX(1) GY(2) = (RY(3) - RY(1))*PHI(2) + (R(3) - R(1))*PHIY(2) . + F(2)*ROY(2) - ROY(3) + ROY(1) GX(3) = (RX(1) - RX(2))*PHI(3) + (R(1) - R(2))*PHIX(3) . + F(3)*ROX(3) - ROX(1) + ROX(2) GY(3) = (RY(1) - RY(2))*PHI(3) + (R(1) - R(2))*PHIY(3) . + F(3)*ROY(3) - ROY(1) + ROY(2) C DO 8 I = 1,3 PX(I) = GX(I) + PHIX(I) PY(I) = GY(I) + PHIY(I) QX(I) = GX(I) - PHIX(I) QY(I) = GY(I) - PHIY(I) 8 CONTINUE C AX(1) = RX(1) + GX(3) - GX(2) AY(1) = RY(1) + GY(3) - GY(2) AX(2) = RX(2) + GX(1) - GX(3) AY(2) = RY(2) + GY(1) - GY(3) AX(3) = RX(3) + GX(2) - GX(1) AY(3) = RY(3) + GY(2) - GY(1) C BX(1) = U(3)*PX(3) + U(2)*QX(2) BY(1) = U(3)*PY(3) + U(2)*QY(2) BX(2) = U(1)*PX(1) + U(3)*QX(3) BY(2) = U(1)*PY(1) + U(3)*QY(3) BX(3) = U(2)*PX(2) + U(1)*QX(1) BY(3) = U(2)*PY(2) + U(1)*QY(1) C CX(1) = V(3)*PX(3) + V(2)*QX(2) CY(1) = V(3)*PY(3) + V(2)*QY(2) CX(2) = V(1)*PX(1) + V(3)*QX(3) CY(2) = V(1)*PY(1) + V(3)*QY(3) CX(3) = V(2)*PX(2) + V(1)*QX(1) CY(3) = V(2)*PY(2) + V(1)*QY(1) C C WX AND WY ARE LINEAR COMBINATIONS OF THE CARDINAL C FUNCTIONS C WX = AX(1)*Z1 + AX(2)*Z2 + AX(3)*Z3 + (BX(1)*ZX1 + . BX(2)*ZX2 + BX(3)*ZX3 + CX(1)*ZY1 + CX(2)*ZY2 + . CX(3)*ZY3)/2. WY = AY(1)*Z1 + AY(2)*Z2 + AY(3)*Z3 + (BY(1)*ZX1 + . BY(2)*ZX2 + BY(3)*ZX3 + CY(1)*ZY1 + CY(2)*ZY2 + . CY(3)*ZY3)/2. RETURN C C VERTICES ARE COLLINEAR C 9 IER = 1 W = 0. IF (IFLAG .NE. 1) RETURN WX = 0. WY = 0. RETURN END SUBROUTINE UNIF (N,X,Y,Z,IADJ,IEND,NROW,NX,NY,PX,PY, . IFLAG, ZXZY, ZZ,IER) INTEGER N, IADJ(1), IEND(N), NROW, NX, NY, IFLAG, IER REAL X(N), Y(N), Z(N), PX(NX), PY(NY), ZXZY(2,N), . ZZ(NROW,NY) C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C GIVEN A THIESSEN TRIANGULATION OF A SET OF POINTS IN THE C PLANE WITH CORRESPONDING DATA VALUES, THIS ROUTINE INTERP- C OLATES THE DATA VALUES TO A SET OF RECTANGULAR GRID POINTS C FOR SUCH APPLICATIONS AS CONTOURING. THE INTERPOLANT IS C ONCE CONTINUOUSLY DIFFERENTIABLE. EXTRAPOLATION IS PER- C FORMED AT GRID POINTS EXTERIOR TO THE TRIANGULATION. C C INPUT PARAMETERS - N - NUMBER OF NODES IN THE TRIANGULA- C TION. N .GE. 3 C C X,Y,Z - N-VECTORS OF NODAL COORDINATES AND C DATA VALUES. C C IADJ,IEND - TRIANGULATION DATA STRUCTURE -- MAY C BE CREATED BY TRMESH. C C NROW - NUMBER OF ROWS IN THE DIMENSION C STATEMENT OF ZZ. C C NX,NY - NUMBER OF ROWS AND COLUMNS, RESPEC- C TIVELY, IN THE RECTANGULAR GRID. C 1 .LE. NX .LE. NROW AND 1 .LE. NY. C C PX,PY - VECTORS OF LENGTH NX AND NY, RE- C SPECTIVELY, CONTAINING THE COORDI- C NATES OF THE GRID LINES. C C IFLAG - OPTION INDICATOR C IFLAG = 0 IF DERIVATIVE ESTIMATES C AT THE VERTICES OF A C TRIANGLE ARE TO BE RECOM- C PUTED FOR EACH GRID POINT C IN THE TRIANGLE AND NOT C SAVED. C IFLAG = 1 IF DERIVATIVE ESTIMATES C ARE INPUT IN ZXZY. C IFLAG = 2 IF DERIVATIVE ESTIMATES C ARE TO BE COMPUTED ONCE C FOR EACH NODE (BY GRADL) C AND SAVED IN ZXZY. C C ZXZY - NOT USED IF IFLAG = 0, 2 BY N ARRAY C WHOSE COLUMNS CONTAIN THE X AND Y C PARTIAL DERIVATIVE ESTIMATES (X C PARTIALS IN THE FIRST ROW) IF C IFLAG = 1, OR 2 BY N ARRAY (OR WORK C SPACE OF LENGTH .GE. 2*N) IF IFLAG C = 2. C C DERIVATIVE ESTIMATES MAY BE COMPUTED BY GRADL OR GRADG. C C ZZ - NROW BY NCOL ARRAY WITH NROW .GE. C NX AND NCOL .GE. NY. C C NONE OF THE INPUT PARAMETERS ARE ALTERED EXCEPT AS C NOTED BELOW. C C OUTPUT PARAMETERS - ZXZY - 2 BY N ARRAY WHOSE COLUMNS CON- C TAIN X AND Y PARTIAL DERIVATIVE C ESTIMATES AT THE NODES IF IFLAG C = 2 AND IER .GE. 0, NOT ALTERED C IF IFLAG .NE. 2. C C ZZ - INTERPOLATED VALUES AT THE GRID C POINTS IF IER .GE. 0. C ZZ(I,J) = F(PX(I),PY(J)) FOR C I = 1,...,NX AND J = 1,...,NY. C C IER - ERROR INDICATOR C IER .GE. 0 IF NO ERRORS WERE C ENCOUNTERED. IER C CONTAINS THE NUMBER C OF GRID POINTS EXT- C ERIOR TO THE TRIAN- C GULATION BOUNDARY. C IER = -1 IF N, NX, NY, OR C IFLAG IS OUT OF C RANGE. C IER = -2 IF THE NODES ARE C COLLINEAR. C C MODULES REFERENCED BY UNIF - INTRC1, TRFIND, TVAL, C (AND OPTIONALLY) GRADL, GETNP, SETUP, GIVENS, C AND ROTATE C C*********************************************************** C INTEGER NST, IST, NN, NI, NJ, IFL, I, J, NIT, IERR, . NEX DATA NST/1/ C C LOCAL PARAMETERS - C C IST = PARAMETER FOR INTRC1 C NST = INITIAL VALUE FOR IST C NN = LOCAL COPY OF N C NI,NJ = LOCAL COPIES OF NX AND NY C IFL = LOCAL COPY OF IFLAG FOR INTRC1 C I,J = DO-LOOP INDICES C IERR = ERROR FLAG FOR CALLS TO GRADL AND INTRC1 C NEX = NUMBER OF GRID POINTS EXTERIOR TO THE TRIANGULA- C TION BOUNDARY (NUMBER OF EXTRAPOLATED VALUES) C NN = N NI = NX NJ = NY IFL = IFLAG IF (NN .LT. 3 .OR. NI .LT. 1 .OR. NI .GT. NROW . .OR. NJ .LT. 1 .OR. IFL .LT. 0 .OR. . IFL .GT. 2) GO TO 4 IST = NST IF (IFL .NE. 2) GO TO 2 C C COMPUTE DERIVATIVE ESTIMATES AT THE NODES. C IFL = 1 DO 1 I = 1,NN CALL GRADL(NN,I,X,Y,Z,IADJ,IEND, ZXZY(1,I), . ZXZY(2,I),IERR) IF (IERR .LT. 0) GO TO 5 1 CONTINUE C C COMPUTE INTERPOLATED VALUES C 2 NEX = 0 DO 3 J = 1,NJ DO 3 I = 1,NI CALL INTRC1(NN,PX(I),PY(J),X,Y,Z,IADJ,IEND,IFL, . ZXZY, IST, ZZ(I,J),IERR) IF (IERR .LT. 0) GO TO 5 NEX = NEX + IERR 3 CONTINUE IER = NEX RETURN C C N, NX, NY, OR IFLAG IS OUT OF RANGE C 4 IER = -1 RETURN C C TRIANGULATION NODES ARE COLLINEAR C 5 IER = -2 RETURN END FUNCTION VOLUME (N,X,Y,Z,IADJ,IEND) INTEGER N, IADJ(1), IEND(N) REAL X(N), Y(N), Z(N) C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C GIVEN A SET OF N DATA POINTS (X(I),Y(I)) AND FUNCTION C VALUES Z(I)=F(X(I),Y(I)) AND A TRIANGULATION COVERING THE C CONVEX HULL H OF THE DATA POINTS, THIS FUNCTION APPROXI- C MATES THE INTEGRAL OF F OVER H BY INTEGRATING THE PIECE- C WISE LINEAR INTERPOLANT OF THE DATA VALUES. VOLUME IS C PART OF AN INTERPOLATION PACKAGE WHICH PROVIDES ROUTINES C TO CREATE, UPDATE, AND PLOT THE TRIANGULAR MESH. C C INPUT PARAMETERS - N - NUMBER OF NODES IN THE MESH. C N .GE. 3. C C X,Y - VECTORS OF COORDINATES OF C THE NODES IN THE MESH. C C Z - VECTOR OF DATA VALUES AT THE C NODES. 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 FUNCTION. C C OUTPUT PARAMETER - VOLUME - SUM OF THE VOLUMES OF THE C LINEAR INTERPOLANTS ON THE C TRIANGLES. C C MODULE REFERENCED BY VOLUME - TRVOL C C*********************************************************** C INTEGER NN, NM2, N1, N2, N3, INDF, INDL, INDX REAL SUM, XN1, YN1, ZN1 C C LOCAL PARAMETERS - C C NN = LOCAL COPY OF N C NM2 = N-2 C N1,N2,N3 = VERTICES OF A TRIANGLE IN COUNTERCLOCKWISE C ORDER C INDF = IADJ INDEX OF THE FIRST NEIGHBOR OF N1 C INDL = IADJ INDEX OF THE LAST NEIGHBOR OF N1 C INDX = IADJ INDEX VARYING FROM INDF TO INDL C SUM = TEMPORARY STORAGE FOR ACCUMULATED VOLUME C XN1,YN1,ZN1 = X(N1), Y(N1), Z(N1) -- STORED LOCALLY FOR C EFFICIENCY C NN = N IF (NN .LT. 3) GO TO 3 C C INITIALIZATION C NM2 = NN-2 INDF = 1 SUM = 0. C C LOOP ON TRIANGLES (N1,N2,N3) SUCH THAT N2 AND N3 ARE C ADJACENT NEIGHBORS OF N1 WHICH ARE BOTH LARGER THAN N1 C DO 2 N1 = 1,NM2 XN1 = X(N1) YN1 = Y(N1) ZN1 = Z(N1) INDL = IEND(N1) DO 1 INDX = INDF,INDL N2 = IADJ(INDX) N3 = IADJ(INDX+1) IF (INDX .EQ. INDL) N3 = IADJ(INDF) IF (N2 .LT. N1 .OR. N3 .LT. N1) GO TO 1 SUM = SUM + TRVOL(XN1,X(N2),X(N3),YN1,Y(N2),Y(N3), . ZN1,Z(N2),Z(N3)) 1 CONTINUE INDF = INDL + 1 2 CONTINUE C VOLUME = SUM RETURN C C N IS OUT OF RANGE C 3 VOLUME = 0. RETURN END C SUBROUTINE ARCTST (N,X,Y,IADJ,IEND, LL, LIST,IER) INTEGER N, IADJ(1), IEND(N), LL, LIST(2,1), IER REAL X(N), Y(N) LOGICAL SWPTST C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C THIS ROUTINE DETERMINES THE NUMBER OF ARCS IN A TRIANGU- C LATION WHICH ARE NOT LOCALLY OPTIMAL AS DEFINED BY LOGICAL C FUNCTION SWPTST. A THIESSEN TRIANGULATION (SEE SUBROUTINE C TRMESH) IS CHARACTERIZED BY THE PROPERTY THAT ALL ARCS ARE C LOCALLY OPTIMAL. THE ALGORITHM CONSISTS OF APPLYING THE C SWAP TEST TO ALL INTERIOR ARCS. C C INPUT PARAMETERS - C C N - NUMBER OF NODES IN THE TRIANGULATION. N .GE. 3. C C X,Y - COORDINATES OF THE NODES. C C IADJ,IEND - DATA STRUCTURE CONTAINING THE TRIANGULA- C TION. SEE SUBROUTINE TRMESH. C C THE ABOVE PARAMETERS ARE NOT ALTERED BY THIS ROUTINE. C C LL - NUMBER OF COLUMNS RESERVED FOR LIST (SEE THE C OUTPUT VALUE OF LIST). LL .GE. 0. C C LIST - INTEGER ARRAY DIMENSIONED 2 BY LL (OR VECTOR C OF LENGTH .GE. 2*LL) IF LL .GT. 0. C C OUTPUT PARAMETERS - C C LL - NUMBER OF ARCS WHICH ARE NOT LOCALLY OPTIMAL C UNLESS IER = 1. C C LIST - COLUMNS CONTAIN THE ENDPOINT NODAL INDICES C (SMALLER INDEX IN THE FIRST ROW) OF THE FIRST C K NONOPTIMAL ARCS ENCOUNTERED, WHERE K IS THE C MINIMUM OF THE INPUT AND OUTPUT VALUES OF LL, C IF IER = 0 AND K .GT. 0. THE NUMBER OF INTE- C RIOR ARCS IS 3N-2NB-3 .LE. 3(N-3) WHERE NB IS C THE NUMBER OF BOUNDARY NODES. BOUNDARY ARCS C ARE OPTIMAL BY DEFINITION. C C IER - ERROR INDICATOR C IER = 0 IF NO ERRORS WERE ENCOUNTERED. C IER = 1 IF N .LT. 3 OR LL .LT. 0. C C MODULES REQUIRED BY ARCTST - SWPTST C C*********************************************************** C NM1 = N - 1 LMAX = LL C C TEST FOR ERRORS AND INITIALIZE FOR LOOP ON INTERIOR ARCS. C IF (NM1 .LT. 2 .OR. LMAX .LT. 0) GO TO 3 IER = 0 L = 0 INDF = 1 C C OUTER LOOP ON NODES IO1 C DO 2 IO1 = 1,NM1 INDL = IEND(IO1) IN2 = IADJ(INDL) C C INNER LOOP ON NEIGHBORS IO2 OF IO1 SUCH THAT IO2 .GT. IO1 C AND IO1-IO2 IS AN INTERIOR ARC -- (IO1,IO2,IN1) AND C (IO2,IO1,IN2) ARE TRIANGLES. C DO 1 INDX = INDF,INDL IO2 = IADJ(INDX) IN1 = IADJ(INDX+1) IF (INDX .EQ. INDL) IN1 = IADJ(INDF) IF (IO2 .LT. IO1 .OR. IN1 .EQ. 0 .OR. . IN2 .EQ. 0) GO TO 1 C C TEST FOR A SWAP. C IF (.NOT. SWPTST(IN1,IN2,IO1,IO2,X,Y)) GO TO 1 L = L + 1 IF (L .GT. LMAX) GO TO 1 LIST(1,L) = IO1 LIST(2,L) = IO2 1 IN2 = IO2 2 INDF = INDL + 1 LL = L RETURN C C N OR LL OUT OF RANGE C 3 IER = 1 RETURN END SUBROUTINE CIRCUM (X1,X2,X3,Y1,Y2,Y3, CX,CY,IER) INTEGER IER REAL X1, X2, X3, Y1, Y2, Y3, CX, CY C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C THIS SUBROUTINE COMPUTES THE COORDINATES OF THE CENTER C OF A CIRCLE DEFINED BY THREE POINTS IN THE PLANE. C C INPUT PARAMETERS - C C X1,...,Y3 - X AND Y COORDINATES OF THREE POINTS IN C THE PLANE. C C INPUT PARAMETERS ARE NOT ALTERED BY THIS ROUTINE. C C OUTPUT PARAMETERS - C C CX,CY - COORDINATES OF THE CENTER OF THE CIRCLE C UNLESS IER = 1. C C IER - ERROR INDICATOR C IER = 0 IF NO ERRORS WERE ENCOUNTERED. C IER = 1 IF THE POINTS ARE COLLINEAR. C C MODULES REQUIRED BY CIRCUM - NONE C C*********************************************************** C REAL U(3), V(3), DS(3) C C SET U(K) AND V(K) TO THE X AND Y COMPONENTS OF THE EDGE C OPPOSITE VERTEX K, TREATING THE POINTS AS VERTICES OF C A TRIANGLE. C U(1) = X3 - X2 U(2) = X1 - X3 U(3) = X2 - X1 V(1) = Y3 - Y2 V(2) = Y1 - Y3 V(3) = Y2 - Y1 C C SET A TO TWICE THE SIGNED AREA OF THE TRIANGLE. A .GT. 0 C IFF (X3,Y3) IS STRICTLY TO THE LEFT OF THE EDGE FROM C (X1,Y1) TO (X2,Y2). C A = U(1)*V(2) - U(2)*V(1) IF (A .EQ. 0.) GO TO 2 C C SET DS(K) TO THE SQUARED DISTANCE FROM THE ORIGIN TO C VERTEX K. C DS(1) = X1**2 + Y1**2 DS(2) = X2**2 + Y2**2 DS(3) = X3**2 + Y3**2 C C COMPUTE FACTORS OF CX AND CY. C FX = 0. FY = 0. DO 1 I = 1,3 FX = FX - DS(I)*V(I) 1 FY = FY + DS(I)*U(I) CX = FX/2./A CY = FY/2./A IER = 0 RETURN C C COLLINEAR POINTS C 2 IER = 1 RETURN END INTEGER FUNCTION LOPTST (N1,N2,X,Y,IADJ,IEND) INTEGER N1, N2, IADJ(1), IEND(1) REAL X(1), Y(1) C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C GIVEN A PAIR OF INDICES DEFINING A TRIANGULATION ARC, C THIS FUNCTION DETERMINES WHETHER OR NOT THE ARC IS LOCALLY C OPTIMAL AS DEFINED BY LOGICAL FUNCTION SWPTST. C C INPUT PARAMETERS - C C N1,N2 - X,Y INDICES OF A PAIR OF ADJACENT NODES. C C X,Y - NODAL COORDINATES. C C IADJ,IEND - DATA STRUCTURE CONTAINING THE TRIANGULA- C TION. SEE SUBROUTINE TRMESH. C C INPUT PARAMETERS ARE NOT ALTERED BY THIS FUNCTION. C C OUTPUT PARAMETER - C C LOPTST = -2 IF N1 AND N2 ARE NOT ADJACENT, C = -1 IF N1-N2 IS AN INTERIOR ARC WHICH IS C NOT LOCALLY OPTIMAL, C = 0 IF N1-N2 SATISFIES THE NEUTRAL CASE (THE C VERTICES OF THE CORRESPONDING QUADRILAT- C ERAL LIE ON A COMMON CIRCLE), C = 1 IF N1-N2 IS A LOCALLY OPTIMAL INTERIOR C ARC, C = 2 IF N1-N2 IS A BOUNDARY ARC. C NOTE THAT N1-N2 IS LOCALLY OPTIMAL IFF LOPTST .GE. 0. C C MODULES REQUIRED BY LOPTST - NONE C C*********************************************************** C IO1 = N1 IO2 = N2 C C FIND THE INDEX OF IO2 AS A NEIGHBOR OF IO1 C INDF = 1 IF (IO1 .GT. 1) INDF = IEND(IO1-1) + 1 INDL = IEND(IO1) 1 INDX = INDX - 1 IF (IADJ(INDX) .EQ. IO2) GO TO 2 IF (INDX .NE. INDF) GO TO 1 C C N1 AND N2 ARE NOT ADJACENT C LOPTST = -2 RETURN C C DETERMINE IN1 AND IN2 SUCH THAT (IO1,IO2,IN1) AND C (IO2,IO1,IN2) ARE TRIANGLES. C 2 IF (INDX .NE. INDL) IN1 = IADJ(INDX+1) IF (INDX .EQ. INDL) IN1 = IADJ(INDF) IF (INDX .NE. INDF) IN2 = IADJ(INDX-1) IF (INDX .EQ. INDF) IN2 = IADJ(INDL) IF (IN1 .NE. 0 .AND. IN2 .NE. 0) GO TO 3 C C N1-N2 IS A BOUNDARY ARC C LOPTST = 2 RETURN C C COMPUTE COMPONENTS OF THE QUADRILATERAL SIDES. C 3 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 IO1-IO2 IS LOCALLY OPTIMAL IFF A1+A2 .LE. 180 DEGREES C WHERE A1 AND A2 DENOTE THE ANGLES AT IN1 AND IN2. C IF (COS1 .LT. 0. .AND. COS2 .LT. 0.) GO TO 4 IF (COS1 .GT. 0. .AND. COS2 .GT. 0.) GO TO 5 C C COMPUTE A QUANTITY WITH THE SIGN OF SIN(A1+A2). C SIN12 = (DX11*DY12 - DX12*DY11)*COS2 + . (DX22*DY21 - DX21*DY22)*COS1 IF (SIN12 .LT. 0.) GO TO 4 IF (SIN12 .GT. 0.) GO TO 5 C C NEUTRAL CASE C LOPTST = 0 RETURN C C N1-N2 NOT LOCALLY OPTIMAL C 4 LOPTST = -1 RETURN C C N1-N2 LOCALLY OPTIMAL C 5 LOPTST = 1 RETURN END * FUNCTION STORE (X) * REAL X C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C THIS FUNCTION FORCES ITS ARGMENT X TO BE STORED IN MAIN C MEMORY. THIS IS USEFUL FOR COMPUTING MACHINE DEPENDENT C PARAMETERS (SUCH AS THE MACHINE PRECISION) WHERE IT IS C NECESSARY TO AVOID COMPUTATION IN HIGH PRECISION REG- C ISTERS. C C INPUT PARAMETER - C C X - VALUE TO BE STORED. C C X IS NOT ALTERED BY THIS FUNCTION. C C OUTPUT PARAMETER - C C STORE - VALUE OF X AFTER IT HAS BEEN STORED AND C (POSSIBLY) TRUNCATED OR ROUNDED TO THE C MACHINE WORD LENGTH. C C MODULES REQUIRED BY STORE - NONE C C*********************************************************** C * COMMON/STCOM/Y * Y = X * STORE = Y * RETURN * END SUBROUTINE TRMTST (N,X,Y,IADJ,IEND,TOL,LUN, IER) INTEGER N, IADJ(1), IEND(N), LUN, IER REAL X(N), Y(N), TOL C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C THIS ROUTINE TESTS THE VALIDITY OF THE DATA STRUCTURE C REPRESENTING A THIESSEN TRIANGULATION CREATED BY SUBROU- C TINE TRMESH. THE FOLLOWING PROPERTIES ARE TESTED -- C 1) IEND(1) .GE. 3 AND IEND(K) .GE. IEND(K-1)+3 FOR K = C 2,...,N (EACH NODE HAS AT LEAST THREE NEIGHBORS). C 2) 0 .LE. IADJ(K) .LE. N FOR K = 1,...,IEND(N) (IADJ C ENTRIES ARE NODAL INDICES OR ZEROS REPRESENTING THE C BOUNDARY). C 3) NB .GE. 3, NT = 2N-NB-2, AND NA = 3N-NB-3 WHERE NB, C NT, AND NA ARE THE NUMBERS OF BOUNDARY NODES, TRI- C ANGLES, AND ARCS, RESPECTIVELY. C 4) EACH CIRCUMCIRCLE DEFINED BY THE VERTICES OF A TRI- C ANGLE CONTAINS NO NODES IN ITS INTERIOR. THIS PROP- C ERTY DISTINGUISHES A THIESSEN TRIANGULATION FROM AN C ARBITRARY TRIANGULATION OF THE NODES. C NOTE THAT NO TEST IS MADE FOR THE PROPERTY THAT A TRIANGU- C LATION COVERS THE CONVEX HULL OF THE NODES, AND THUS A C TEST ON A DATA STRUCTURE ALTERED BY SUBROUTINE DELETE C SHOULD NOT RESULT IN AN ERROR. C C INPUT PARAMETERS - C C N - NUMBER OF NODES. N .GE. 3. C C X,Y - NODAL COORDINATES. C C IADJ,IEND - TRIANGULATION DATA STRUCTURE. SEE SUB- C ROUTINE TRMESH. C C TOL - NONNEGATIVE TOLERANCE TO ALLOW FOR FLOATING- C POINT ERRORS IN THE CIRCUMCIRCLE TEST. AN C ERROR SITUATION IS DEFINED AS R**2 - D**2 .GT. C TOL WHERE R IS THE RADIUS OF A CIRCUMCIRCLE C AND D IS THE DISTANCE FROM THE CIRCUMCENTER C TO THE NEAREST NODE. A REASONABLE VALUE C FOR TOL IS 10*EPS WHERE EPS IS THE MACHINE C PRECISION. THE TEST IS EFFECTIVELY BYPASSED C BY MAKING TOL LARGER THAN THE DIAMETER OF THE C CONVEX HULL OF THE NODES. C C LUN - LOGICAL UNIT NUMBER FOR PRINTING ERROR MES- C SAGES. IF LUN .LT. 1 OR LUN .GT. 99, NO MES- C SAGES ARE PRINTED. C C INPUT PARAMETERS ARE NOT ALTERED BY THIS ROUTINE. C C OUTPUT PARAMETER - C C IER - ERROR INDICATOR C IER = -1 IF ONE OR MORE NULL TRIANGLES (AREA = C 0) ARE PRESENT BUT NO (OTHER) ERRORS C WERE ENCOUNTERED. A NULL TRIANGLE IS C AN ERROR ONLY IF IT OCCURS IN THE C THE INTERIOR. C IER = 0 IF NO ERRORS OR NULL TRIANGLES WERE C ENCOUNTERED. C IER = 1 IF N .LT. 3 OR TOL .LT. 0. C IER = 2 IF AN IEND OR IADJ ENTRY IS OUT OF C RANGE. C IER = 3 IF THE TRIANGULATION PARAMETERS (NB, C NT, AND NA) ARE INCONSISTENT. C IER = 4 IF A TRIANGLE CONTAINS A NODE INTERIOR C TO ITS CIRCUMCIRCLE. C THE ERROR SITUATIONS ARE TESTED IN THE ORDER C DEFINED BY THE (POSITIVE) IER VALUES. C C MODULE REQUIRED BY TRMTST - CIRCUM C C*********************************************************** C LOGICAL RITE C C SET LOCAL VARIABLES, TEST FOR ERRORS IN INPUT, AND C INITIALIZE COUNTS. C NN = N RTOL = TOL RITE = LUN .GE. 1 .AND. LUN .LE. 99 IF (NN .LT. 3 .OR. RTOL .LT. 0.) GO TO 5 NB = 0 NT = 0 NULL = 0 NFAIL = 0 C C LOOP ON TRIANGLES (N1,N2,N3) SUCH THAT N2 AND N3 INDEX C ADJACENT NEIGHBORS OF N1 AND ARE BOTH LARGER THAN N1 C (TRIANGLES ARE ASSOCIATED WITH THEIR SMALLEST INDEX). C INDF = 1 DO 4 N1 = 1,NN INDL = IEND(N1) IF (INDL .LT. INDF+2) GO TO 6 IF (IADJ(INDL) .EQ. 0) NB = NB + 1 C C LOOP ON NEIGHBORS OF N1 C DO 3 INDX = INDF,INDL N2 = IADJ(INDX) IF (N2 .LT. 0 .OR. N2 .GT. NN .OR. . (INDX .LT. INDL .AND. N2 .EQ. 0)) GO TO 7 IF (INDX .LT. INDL) N3 = IADJ(INDX+1) IF (INDX .EQ. INDL) N3 = IADJ(INDF) IF (N2 .LT. N1 .OR. N3 .LT. N1) GO TO 3 NT = NT + 1 C C COMPUTE THE COORDINATES OF THE CIRCUMCENTER OF C (N1,N2,N3). C CALL CIRCUM (X(N1),X(N2),X(N3),Y(N1),Y(N2),Y(N3), . CX,CY,IERR) IF (IERR .NE. 0) NULL = NULL + 1 IF (IERR .NE. 0) GO TO 3 C C TEST FOR NODES WITHIN THE CIRCUMCIRCLE. C R = (CX-X(N1))**2 + (CY-Y(N1))**2 - RTOL IF (R .LE. 0.) GO TO 3 DO 1 I = 1,NN IF (I .EQ. N1 .OR. I .EQ. N2 .OR. . I .EQ. N3) GO TO 1 IF ((CX-X(I))**2 + (CY-Y(I))**2 .LT. R) GO TO 2 1 CONTINUE GO TO 3 C C NODE I IS INTERIOR TO THE CIRCUMCIRCLE OF (N1,N2,N3). C 2 NFAIL = NFAIL + 1 3 CONTINUE 4 INDF = INDL + 1 C C CHECK PARAMETERS FOR CONSISTENCY AND TEST FOR NFAIL = 0. C NA = (IEND(NN)-NB)/2 IF (NB .LT. 3 .OR. NT .NE. 2*NN-NB-2 .OR. . NA .NE. 3*NN-NB-3) GO TO 8 IF (NFAIL .NE. 0) GO TO 9 C C NO ERRORS WERE ENCOUNTERED. C IER = 0 IF (NULL .EQ. 0) RETURN IER = -1 IF (RITE) WRITE (LUN,100) NULL 100 FORMAT (1H ,10HTRMTST -- ,I5,16H NULL TRIANGLES , . 11HARE PRESENT/1H ,10X,16H(NULL TRIANGLES , . 32HON THE BOUNDARY ARE UNAVOIDABLE)//) RETURN C C N OR TOL IS OUT OF RANGE. C 5 IER = 1 IF (RITE) WRITE (LUN,110) N, RTOL 110 FORMAT (1H ,33HTRMTST -- INVALID INPUT PARAMETER/ . 1H ,10X,4HN = ,I5,8H, TOL = ,E8.1) RETURN C C IEND ENTRY OUT OF RANGE C 6 IER = 2 IF (RITE) WRITE (LUN,120) N1 120 FORMAT (1H ,15HTRMTST -- NODE ,I5, . 26H HAS LESS THAN 3 NEIGHBORS) RETURN C C IADJ ENTRY OUT OF RANGE C 7 IER = 2 IF (RITE) WRITE (LUN,130) INDX 130 FORMAT (1H ,33HTRMTST -- IADJ(K) IS NOT A VALID , . 14HINDEX FOR K = ,I5) RETURN C C INCONSISTENT TRIANGULATION PARAMETERS C 8 IER = 3 IF (RITE) WRITE (LUN,140) NB, NT, NA 140 FORMAT (1H ,33HTRMTST -- INCONSISTENT PARAMETERS/ . 1H ,10X,I5,15H BOUNDARY NODES,3X,I5, . 10H TRIANGLES,3X,I5,5H ARCS) RETURN C C CIRCUMCIRCLE TEST FAILURE C 9 IER = 4 IF (RITE) WRITE (LUN,150) NFAIL 150 FORMAT (1H ,10HTRMTST -- ,I5,15H CIRCUMCIRCLES , . 32HCONTAIN NODES IN THEIR INTERIORS) RETURN END