/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:30:53 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "dtggrd.h" #include /* PARAMETER translations */ #define EIGHT 8.0e0 #define ZERO 0.0e0 /* end of PARAMETER translations */ void /*FUNCTION*/ dtggrd( double x[], double y[], long np, long ip[], double w[], long triang[], long mt, long b[][4], long mb, long *nt, long info[]) { LOGICAL32 ok, ok1; long int avail, ext, i, i1, i2, i3, i3prev, index, isave, jhi, jlo, jnew, k, kb1, kbhi, kblo, kbnew, kbz, mbused, mode, nb, ntlim, tri[6]; double aref, cross, dx12, dy12, pangle, xmean, ymean; /* OFFSET Vectors w/subscript range: 1 to dimension */ long *const Info = &info[0] - 1; long *const Ip = &ip[0] - 1; long *const Tri = &tri[0] - 1; long *const Triang = &triang[0] - 1; double *const W = &w[0] - 1; double *const X = &x[0] - 1; double *const Y = &y[0] - 1; /* end of OFFSET VECTORS */ /* Copyright (c) 1996 California Institute of Technology, Pasadena, CA. * ALL RIGHTS RESERVED. * Based on Government Sponsored Research NAS7-03001. * File: dtggrd.[f|for] * contains dtggrd, dtgang, dtgadj, dtgupd * *>> 1997-07-01 DTGGRD Krogh Reversed subscripts in B (CLL suggestion) *>> 1997-06-14 DTGGRD Krogh Removed implicit none, added external stmt. *>> 1997-06-14 DTGGRD CLL Argument INFO() replaces NBUSED, MODE, & OK1. *>> 1997-06-09 DTGGRD Krogh Converted SFTRAN code to Fortran. *>> 1995-09-26 DTGGRD CLL Editing for inclusion into MATH77. *>> 1991-11-20 DTGGRD CLL Using Math77 subroutine _SORTP for sorting. * * THIS SUBR CONSTRUCTS A LIST OF POINTERS IN THE ARRAY * TRIANG() DEFINING A TRIANGULAR GRID HAVING THE GIVEN (X,Y) * DATA AS VERTICES. * THE GRID WILL BE OPTIMAL IN THE SENSE OF THE MAX-MIN ANGLE * CRITERION. IT HAS BEEN SHOWN, INDEPENDENTLY BY LAWSON AND BY * R. SIBSON, THAT THIS IS EQUIVALENT TO SOME OTHER OPTIMALLITY * CRITERIA. THE GRID PRODUCED IS A DELAUNAY TRIANGULATION AND * IS THE GRAPH THEORETIC DUAL OF A TESSELLATION MADE UP OF * PROXIMITY REGIONS ASSOCIATED WITH THE NAMES DIRICHLET, * VORONOI, AND THIESSEN. * * C. L. LAWSON, JPL, 1976 NOV 28 * C.L.L., 1979 MAR 3. CONVERTED TO SFTRAN3. * C.L.L., 1979 JUL 22. ADDED CALL TO DTGSIZ. * * X(1:NP),Y(1:NP) [in] (X,Y) DATA POINTS. * NP [in] NO. OF DATA POINTS. REQUIRE NP .GE. 3. * * IP() [scratch] WORK SPACE OF LENGTH at least NP. * W() [scratch] WORK SPACE OF LENGTH at least NP. * TRIANG() [out] ARRAY IN WHICH THIS SUBR WILL BUILD A LIST DEFINING * A triangular grid. access to this list is via the * four subrs DTGGET, DTGPUT, DTGSIZ, and DTGSET. * MT [in] DIMENSION OF ARRAY TRIANG(). NO. OF TRIANGLES IT * CAN HOLD, NTLIM, IS DETERMINED BY USE OF SUBROUTINE * DTGSIZ. NTLIM WILL BE MT/6 IF POINTERS ARE NOT * PACKED. NTLIM = (3*MT)/6 IF POINTERS ARE PACKED * THREE PER WORD. * Bdry(1:4, MB) [integer, out] Array containing pointers defining the * boundary of the (convex) triangular grid. * Bdry(1, K) = FWD POINTER. Points to next vertex in * counterclockwise order. * Bdry(2, K) = BACKWARD POINTER. Points to next vertex in * clockwise order. * Bdry(3, K) = A BOUNDARY POINT * Bdry(4, K) = A BOUNDARY TRIANGLE * The triangle Bdry(4, K) has a boundary edge that * connects the points Bdry(3, K) and Bdry(3, K+1). * * On return from this subroutine, some elements of the array Bdry() * may not be members of the linked list defining the boundary. The * entry with K = 1 will always be a member of the boundary list. * Thus for another subroutine to scan the Bdry() list made by this * subroutine it should start with K = 1 and follow the forward or * backward pointers. * * MB [in] FIRST DIMENSION OF B(,). SETTING MB = NP+1 WILL * ALWAYS BE ADEQUATE, HOWEVER A MUCH SMALLER VALUE * WILL OFTEN SUFFICE. SUGGEST TRYING MB = 6 * (cube root * of NP) WHEN NP .GE. 64. * NT [out] OUTPUT.. NO. OF TRIANGLES * * INFO(1:3) [out, integer] Termination status information. * * INFO(1) INDICATES STATUS ON TERMINATION. * 0 NORMAL TERMINATION. The triangular grid is complete * and successfully optimized. * 1 ERROR. The triangular grid is complete but not * successfully optimized. Apparent looping in DTGADJ. * GRID MAY NOT BE OPTIMAL. * 2 ERROR. ALL GIVEN POINTS ARE COLINEAR. NO TRIANGLES * CONSTRUCTED. (OR ELSE THERE ARE SOME DUPLICATE PTS) * 3 ERROR. DUPLICATE POINTS. TRIANGULAR GRID NOT * COMPLETED. * 4 ERROR. NOT ENOUGH SPACE IN B(,). USER MUST * INCREASE MB. * 5 NOT ENOUGH SPACE IN TRIANG(). USER MUST MAKE * MT BIGGER. * INFO(2) NO. OF BOUNDARY POINTS. * INFO(3) SMALLEST SETTING OF MB THAT WOULD SUFFICE FOR THIS DATA. * * [Relations to internal variables or previously used variable * names: INFO(1) = MODE, INFO(2) = NB, INFO(3) = MBUSED, * The value of INFO(1) = MODE will be set to 0 when * OK1 = .true. and 1 when OK1 = .false.] * * USE SUBR DTGSIZ TO COMPUTE NTLIM, THE MAX NO. OF * TRIANGLES, GIVEN THE DIMENSION PARAMETER MT. * * ------------------------------------------------------------------ *--D replaces "?": ?TGGRD, ?TGANG, ?TGSET, ?TGGET, ?TGPUT, ?TGSIZ, *--& ?TGADJ, ?SORTP, ?TGUPD * ------------------------------------------------------------------ */ /* ------------------------------------------------------------------ */ dtgsiz( mt, &ntlim ); mode = 0; ok1 = TRUE; /* FIND AN EXTREME POINT (X(EXT),Y(EXT)) * */ ext = 1; for (i = 1; i <= np; i++) { if (X[i] < X[ext] || (X[i] == X[ext] && Y[i] < Y[ext])) ext = i; } /* COMPUTE SQUARED DISTANCES FROM (X(EXT),Y(EXT)). * ALSO INITIALIZE IP(). */ for (i = 1; i <= np; i++) { W[i] = powi(X[i] - X[ext],2) + powi(Y[i] - Y[ext],2); } /* Use _SORTP from the JPL MATH77 library to * determine a sorted order for the contents of W() and set * indices indicating the sorted order in IP(). * IP() will be set so that W(IP(1)) is the smallest number in * W(), W(IP(2)) is the next larger number, etc. The contents * of IP() will be used as indices into the X() and Y() arrays. * The contents of W() will not be reordered by _SORTP, however * this does not matter here since the next time we reference W() * will be to store other things into it. * */ dsortp( w, 1, np, ip ); /* BUILD THE FIRST TRIANGLE USING POINTS IP(1), IP(2), AND * THE NEXT POINT THAT IS NOT COLINEAR WITH THESE TWO. * */ i1 = Ip[1]; i2 = Ip[2]; dx12 = X[i2] - X[i1]; dy12 = Y[i2] - Y[i1]; i3prev = 0; for (i = 3; i <= np; i++) { i3 = Ip[i]; Ip[i] = i3prev; cross = -(X[i3] - X[i2])*dy12 + (Y[i3] - Y[i2])*dx12; if (cross != 0.) { if (cross < 0.) { /* SWAP I2 AND I3 */ isave = i2; i2 = i3; i3 = isave; } goto L_50; } i3prev = i3; } /* ERROR IF DROP THRU TO HERE. ALL NP POINTS ARE COLINEAR * OR ELSE POINTS I1 AND I2 ARE IDENTICAL. * */ mode = 2; goto L_999; L_50: ; /* HERE THE VERTICES OF THE FIRST TRIANGLE ARE I1,I2,I3 * * COMPUTE MEAN OF VERTICES OF FIRST TRIANGLE. * */ xmean = (X[i1] + X[i2] + X[i3])/3.; ymean = (Y[i1] + Y[i2] + Y[i3])/3.; /* COMPUTE THE PSEUDOANGLE OF VECTOR FROM (XMEAN,YMEAN) * TO I1. * */ aref = dtgang( X[i1] - xmean, Y[i1] - ymean, ZERO ); /* INITIALIZE BOUNDARY POINTER STRUCTURE B(,). * B(1, K) = FWD POINTER * B(2, K) = BACKWARD POINTER * B(3, K) = A BOUNDARY POINT * B(4, K) = A BOUNDARY TRIANGLE * W(K) = PSEUDOANGLE OF POINT B(3, K) * */ for (k = 1; k <= mb; k++) { b[k - 1][0] = k + 1; } b[mb - 1][0] = 0; b[3][0] = 0; avail = 5; for (k = 1; k <= 4; k++) { b[k - 1][1] = k - 1; b[k - 1][3] = 1; } b[3][3] = 0; b[0][2] = i1; b[1][2] = i2; b[2][2] = i3; b[3][2] = i1; W[1] = ZERO; W[2] = dtgang( X[i2] - xmean, Y[i2] - ymean, aref ); W[3] = dtgang( X[i3] - xmean, Y[i3] - ymean, aref ); W[4] = EIGHT; nb = 3; mbused = 3; /* BUILD FIRST TRIANGLE IN TRIANG(). * */ if (ntlim <= 0) { mode = 5; goto L_999; } dtgset( 1, 0, 0, 0, i1, i2, i3, triang, mt ); *nt = 1; /* ****************************************************************** * * BEGIN MAIN LOOP. ADD ONE NEW POINT AT A TIME INTO * THE TRIANGULAR GRID STRUCTURE. * * ****************************************************************** */ kbnew = 2; for (index = 4; index <= np; index++) { jnew = Ip[index]; pangle = dtgang( X[jnew] - xmean, Y[jnew] - ymean, aref ); /* USE PANGLE AS SEARCH KEY IN THE TABLE OF * PSEUDO ANGLES OF BOUNDARY POINTS. * START SEARCH AT PREVIOUS NEW BOUNDARY POINT. * */ kb1 = kbnew; if (pangle > W[kb1]) { kbhi = kb1; L_100: kblo = kbhi; kbhi = b[kblo - 1][0]; if (pangle > W[kbhi]) goto L_100; } else { /* HERE WE HAVE PANGLE .le. W(KB1) */ kblo = kb1; L_110: kbhi = kblo; kblo = b[kbhi - 1][1]; if (pangle < W[kblo]) goto L_110; } jlo = b[kblo - 1][2]; jhi = b[kbhi - 1][2]; /* TEST FOR DUPLICATE POINTS * */ if (pangle == W[kblo]) { if (X[jnew] == X[jlo] && Y[jnew] == Y[jlo]) { mode = 3; goto L_999; } } else if (pangle == W[kbhi]) { if (X[jnew] == X[jhi] && Y[jnew] == Y[jhi]) { mode = 3; goto L_999; } } /* Attach point JNEW to JHI and JLOW. */ if (avail == 0) { /* ERROR.. INSUFFICIENT STORAGE IN B() * TO REPRESENT THE BOUNDARY. */ mode = 4; goto L_999; } if (*nt >= ntlim) { /* ERROR.. INSUFFICIENT STORAGE IN * TRIANG() */ mode = 5; goto L_999; } /* UPDATE TRIANG(). */ *nt += 1; dtgupd( b[kblo - 1][2], b[kblo - 1][3], tri, triang, mt, *nt ); dtgset( *nt, b[kblo - 1][3], 0, 0, jhi, jlo, jnew, triang, mt ); /* UPDATE B() AND W(). */ nb += 1; if (nb > mbused) mbused = nb; kbnew = avail; avail = b[avail - 1][0]; b[kbnew - 1][0] = kbhi; b[kbnew - 1][1] = kblo; b[kbnew - 1][2] = jnew; b[kbnew - 1][3] = *nt; W[kbnew] = pangle; b[kblo - 1][0] = kbnew; b[kbhi - 1][1] = kbnew; b[kblo - 1][3] = *nt; dtgadj( *nt, 1, x, y, np, triang, mt, b, mb, kbnew, &ok ); ok1 = ok1 && ok; /* Done attaching point JNEW to JHI and JLOW. * We have just assigned JNEW the index KBNEW in the boundary list. * * LOOP BACKWARD ALONG BOUNDARY FROM KBNEW TO FIND OTHER * POINTS TO WHICH POINT JNEW CAN BE CONNECTED. * */ kbhi = kblo; kblo = b[kbhi - 1][1]; L_130: if (kblo != 0) { jhi = b[kbhi - 1][2]; jlo = b[kblo - 1][2]; cross = -(Y[jhi] - Y[jnew])*(X[jlo] - X[jhi]) + (X[jhi] - X[jnew])*(Y[jlo] - Y[jhi]); if (cross > 0.) { /* CONNECT POINT JNEW TO JLO */ if (*nt >= ntlim) { /* ERROR.. INSUFFICIENT SPACE IN TRIANG(). * */ mode = 5; goto L_999; } /* UPDATE TRIANG() */ *nt += 1; dtgupd( b[kblo - 1][2], b[kblo - 1][3], tri, triang, mt, *nt ); dtgupd( b[kbhi - 1][2], b[kbhi - 1][3], tri, triang, mt, *nt ); dtgset( *nt, b[kblo - 1][3], 0, b[kbhi - 1][3], jhi, jlo, jnew, triang, mt ); /* UPDATE B() */ b[kblo - 1][0] = kbnew; b[kbnew - 1][1] = kblo; b[kblo - 1][3] = *nt; nb -= 1; b[kbhi - 1][0] = avail; avail = kbhi; dtgadj( *nt, 1, x, y, np, triang, mt, b, mb, kbnew, &ok ); ok1 = ok1 && ok; /* Done connecting point JNEW to JLO * */ kbhi = kblo; kblo = b[kbhi - 1][1]; goto L_130; } } /* LOOP FORWARD ALONG BOUNDARY FROM KBNEW TO FIND OTHER * POINTS TO WHICH POINT JNEW CAN BE CONNECTED. * */ kblo = b[kbnew - 1][0]; kbhi = b[kblo - 1][0]; L_150: if (kbhi != 0) { jhi = b[kbhi - 1][2]; jlo = b[kblo - 1][2]; cross = -(Y[jlo] - Y[jhi])*(X[jnew] - X[jlo]) + (X[jlo] - X[jhi])*(Y[jnew] - Y[jlo]); if (cross > 0) { if (*nt >= ntlim) { /* ERROR.. INSUFFICIENT SPACE IN TRIANG(). */ mode = 5; goto L_999; } /* Connect point JNEW to JHI -- Update TRIANG() */ *nt += 1; dtgupd( b[kblo - 1][2], b[kblo - 1][3], tri, triang, mt, *nt ); dtgupd( b[kbnew - 1][2], b[kbnew - 1][3], tri, triang, mt, *nt ); dtgset( *nt, b[kblo - 1][3], b[kbnew - 1][3], 0, jhi, jlo, jnew, triang, mt ); /* UPDATE B() */ b[kbnew - 1][0] = kbhi; b[kbhi - 1][1] = kbnew; b[kbnew - 1][3] = *nt; nb -= 1; b[kblo - 1][0] = avail; avail = kblo; dtgadj( *nt, 1, x, y, np, triang, mt, b, mb, kbnew, &ok ); ok1 = ok1 && ok; /* Done connecting point JNEW to JHI * */ kblo = kbhi; kbhi = b[kblo - 1][0]; goto L_150; } } } mbused += 1; L_999: ; /* As built, the linked boundary list runs (counterclockwise) from * b(*, 1) to b(*, 4), with b(*, 1) and b(*, 4) representing the * same point. The back pointer in b(*, 1) and the forward * pointer in b(*, 4) are both zero, indicating ends of the list. * Change this so the linked list is fully circular, still * including b(*, 1), but detaching b(*, 4). Let kbz be the value * of the back pointer in b(*, 4). Set the back pointer in b(*, * 1) to be kbz, and the forward pointer in b(*, kbz) to be 1. * */ kbz = b[3][1]; b[0][1] = kbz; b[kbz - 1][0] = 1; if (mode == 0 && !ok1) mode = 1; Info[1] = mode; Info[2] = nb; Info[3] = mbused; return; } /* end of function */ /* ================================================================== */ void /*FUNCTION*/ dtgupd( long jfix, long kb4, long tri[], long triang[], long mt, long nt) { long int i; /* OFFSET Vectors w/subscript range: 1 to dimension */ long *const Tri = &tri[0] - 1; long *const Triang = &triang[0] - 1; /* end of OFFSET VECTORS */ /* Updates TRIANG, using JFIX and KB4 */ dtgget( kb4, tri, triang ); for (i = 1; i <= 3; i++) { if (Tri[i + 3] == jfix) { Tri[i] = nt; goto L_20; } } L_20: dtgput( kb4, tri, triang, mt ); return; } /* end of function */ /* ================================================================== */ /* PARAMETER translations */ #define FOUR 4.0e0 #define TWO 2.0e0 /* end of PARAMETER translations */ double /*FUNCTION*/ dtgang( double xx, double yy, double aref) { double a, dtgang_v; /*>> 1995-09-26 CLL Editing for inclusion into MATH77. * * Compute the pseudoangle DTGANG between the reference direction * whose pseudoangle is AREF and the vector (XX,YY). * Pseudoangle is measured counterclockwise with period 8. * Require 0. .le. AREF .lt. 8. on input. * On return DTGANG will be in the interval 0. .le. DTGANG .lt. 8. * * C.L.LAWSON, JPL, 1976 NOV 8 * * ------------------------------------------------------------------ */ /* ------------------------------------------------------------------ */ if (fabs( xx ) > fabs( yy )) { a = yy/xx; if (xx < ZERO) { a += FOUR; } else { if (yy < ZERO) { a += EIGHT; } } } else { if (yy == ZERO) { a = ZERO; } else { a = TWO - xx/yy; if (yy < ZERO) a += FOUR; } } /* . A IS NOW THE PSEUDOANGLE OF (XX,YY) * . RELATIVE TO (1.,0.) AND LYING IN * . THE INTERVAL 0. .le. A .lt. 8. */ dtgang_v = a - aref; if (dtgang_v < ZERO) dtgang_v += EIGHT; return( dtgang_v ); } /* end of function */ /* ================================================================== */ void /*FUNCTION*/ dtgadj( long tin, long nin, double x[], double y[], long np, long triang[], long mt, long b[][4], long mb, long kbnew, LOGICAL32 *ok) { LOGICAL32 swap; long int count, i, k, kbkwd, kfwd, more, nabor1, nabor2, qv1, qv2, qv3, qv4, t1, t11, t12, t2, t21, t22, tnew, told, tri1[6], tri2[6], _i, _r; static long int add1[3], qv1s[3], qv2s[3], qv3s[3], sub1[3]; double a123, a234, a341, a412, diag13, diag24, dx12, dx23, dx34, dx41, dy12, dy23, dy34, dy41, s12, s23, s34, s41, tau13, tau24; static int _aini = 1; /* OFFSET Vectors w/subscript range: 1 to dimension */ long *const Add1 = &add1[0] - 1; long *const Qv1s = &qv1s[0] - 1; long *const Qv2s = &qv2s[0] - 1; long *const Qv3s = &qv3s[0] - 1; long *const Sub1 = &sub1[0] - 1; long *const Tri1 = &tri1[0] - 1; long *const Tri2 = &tri2[0] - 1; long *const Triang = &triang[0] - 1; double *const X = &x[0] - 1; double *const Y = &y[0] - 1; /* end of OFFSET VECTORS */ if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */ Qv1s[1] = 5; Qv1s[2] = 6; Qv1s[3] = 4; Qv2s[1] = 6; Qv2s[2] = 4; Qv2s[3] = 5; Qv3s[1] = 4; Qv3s[2] = 5; Qv3s[3] = 6; Add1[1] = 2; Add1[2] = 3; Add1[3] = 1; Sub1[1] = 3; Sub1[2] = 1; Sub1[3] = 2; _aini = 0; } /*>> 1995-09-26 CLL Editing for inclusion into MATH77. * * C.L.LAWSON, JPL, 1976 NOV 30, Changed 1977 APR 7 * Changed 1977 June 16 to get rid of internal arrays * used to stack triangles to be tested. * * Given triangle TIN and the neighboring triangle that is in * position NIN relative to TIN, Test for possible exchange of * their common edge. * Use the criterion of maximizing the smallest angle, which is the * same as the criterion of the empty circumcircle. * If the swap is made then certain other edges must be tested. * This subr does the tests, swaps edges as needed, and does * tests and swaps on other edges as needed. * The counter COUNT is used to guard against infinite looping * which should never happen anyway. If count exceeds 100 the * subr sets OK = .false. indicating that the triangulation * produced may not be optimal. * ------------------------------------------------------------------ * Subroutine Arguments * * TIN [integer, in] * NIN [integer, in] * X(),Y() [floating, in] * TRIANG() [integer, inout] * MT [integer, in] * B(1:4, ) [integer, inout] * MB [integer, in] * KBNEW [integer, in] Index into the boundary array, B(,). * OK [logical, out] * ------------------------------------------------------------------ */ /* ------------------------------------------------------------------ */ t1 = tin; nabor1 = nin; count = 1; more = 1; *ok = TRUE; L_20: if ((more > 0) && *ok) { dtgget( t1, tri1, triang ); t2 = Tri1[nabor1]; if (t2 != 0) { /* T2=0 MEANS SIDE NABOR OF TRIANGLE T1 IS A BOUNDARY EDGE SO THERE IS NO * TRIANGLE T2. HERE WE HAVE T2 .ne. 0 */ dtgget( t2, tri2, triang ); /* SET QV1,...,QV4 AS INDICES OF THE 4 VERTICES OF THE QUADRILATERAL * FORMED BY THE 2 TRIANGLES T1 AND T2 */ i = Qv1s[nabor1]; qv1 = Tri1[i]; i = Qv2s[nabor1]; qv2 = Tri1[i]; i = Qv3s[nabor1]; qv3 = Tri1[i]; for (nabor2 = 1; nabor2 <= 3; nabor2++) { if (Tri2[nabor2] == t1) { i = Qv2s[nabor2]; qv4 = Tri2[i]; goto L_50; } } L_50: ; /* The max-min angle test below is made using the max-min angle criterion * as described in JPL internal TECH. MEMO. 914-299, FEB., 1972 by * C. L. LAWSON. * */ dx23 = X[qv3] - X[qv2]; dy23 = Y[qv3] - Y[qv2]; dx34 = X[qv4] - X[qv3]; dy34 = Y[qv4] - Y[qv3]; a234 = dx23*dy34 - dy23*dx34; if (a234 <= 0.) { swap = FALSE; } else { dx12 = X[qv2] - X[qv1]; dy12 = Y[qv2] - Y[qv1]; dx41 = X[qv1] - X[qv4]; dy41 = Y[qv1] - Y[qv4]; a412 = dx41*dy12 - dy41*dx12; if (a412 <= 0.) { swap = FALSE; } else { a341 = dx34*dy41 - dy34*dx41; a123 = a234 + a412 - a341; s12 = SQ(dx12) + SQ(dy12); s23 = SQ(dx23) + SQ(dy23); s34 = SQ(dx34) + SQ(dy34); s41 = SQ(dx41) + SQ(dy41); diag13 = powi(X[qv3] - X[qv1],2) + powi(Y[qv3] - Y[qv1],2); diag24 = powi(X[qv4] - X[qv2],2) + powi(Y[qv4] - Y[qv2],2); tau13 = fmin( SQ(a123)/fmax( s12, s23 ), SQ(a341)/ fmax( s34, s41 ) )/diag13; tau24 = fmin( SQ(a234)/fmax( s23, s34 ), SQ(a412)/ fmax( s41, s12 ) )/diag24; swap = tau24 > tau13; } } /* End of max-min angle test */ } else { swap = FALSE; } if (swap) { /* Remove the edge connecting QV1 and QV3. Insert new edge connecting * QV2 and QV4. Redefine triangle T1 to have vertices QV1, QV2, and * QV4. Redefine triangle T2 to have vertices QV2, QV3, AND QV4. */ i = Add1[nabor1]; t11 = Tri1[i]; i = Add1[i]; t12 = Tri1[i]; i = Add1[nabor2]; t21 = Tri2[i]; i = Add1[i]; t22 = Tri2[i]; dtgset( t1, t11, t2, t22, qv1, qv2, qv4, triang, mt ); dtgset( t2, t12, t21, t1, qv2, qv3, qv4, triang, mt ); tnew = t2; if (t12 != 0) { dtgget( t12, tri2, triang ); for (i = 1; i <= 3; i++) { if (Tri2[i] == t1) { Tri2[i] = tnew; goto L_80; } } L_80: dtgput( t12, tri2, triang, mt ); } else { /* The following test of (QV2 .eq. B(3, KBNEW)) will always be true if * the outer algorithm calling this subroutine is the one originally * planned, always adding new points from the exterior, not from the * interior. The false branch is provided to handle the more general * case. */ if (qv2 == b[kbnew - 1][2]) { b[kbnew - 1][3] = tnew; } else { /* WRITE(6,1001) QV2,B(3, KBNEW) * Search alternately fwd and bkwd from KBNEW in the boundary list B() * for the point QV2. When found set the associated triangle to TNEW. */ kfwd = kbnew; kbkwd = b[kfwd - 1][1]; L_90: ; if (kfwd != 0) { if (b[kfwd - 1][2] == qv2) { b[kfwd - 1][3] = tnew; goto L_100; } kfwd = b[kfwd - 1][0]; } if (kbkwd != 0) { if (b[kbkwd - 1][2] == qv2) { b[kbkwd - 1][3] = tnew; goto L_100; } kbkwd = b[kbkwd - 1][1]; } goto L_90; } } L_100: tnew = t1; if (t22 != 0) { dtgget( t22, tri2, triang ); for (i = 1; i <= 3; i++) { if (Tri2[i] == t2) { Tri2[i] = tnew; goto L_120; } } L_120: dtgput( t22, tri2, triang, mt ); } else { /* WRITE(6,1002) QV4 * Search alternately fwd and bkwd from KBNEW in the boundary list B() * for the point QV2. When found set the associated triangle to TNEW. */ kfwd = kbnew; kbkwd = b[kfwd - 1][1]; L_140: ; if (kfwd != 0) { if (b[kfwd - 1][2] == qv4) { b[kfwd - 1][3] = tnew; goto L_150; } kfwd = b[kfwd - 1][0]; } if (kbkwd != 0) { if (b[kbkwd - 1][2] == qv4) { b[kbkwd - 1][3] = tnew; goto L_150; } kbkwd = b[kbkwd - 1][1]; } goto L_140; } L_150: nabor1 = 3; more += 1; count += 1; *ok = count <= 100; } else { more -= 1; if (more > 0) { told = t1; i = Sub1[nabor1]; t1 = Tri1[i]; dtgget( t1, tri1, triang ); for (k = 1; k <= 3; k++) { if (Tri1[k] == told) { nabor1 = Sub1[k]; goto L_20; } } } } goto L_20; } return; /*1001 FORMAT(33H0ADJUST.. USING FIX-B WITH QV2 =,I5,5X, * * 12HB(3, KBNEW) =,I5) *1002 FORMAT(33H0ADJUST.. USING FIX-B WITH QV4 =,I5) */ } /* end of function */