subroutine addcst ( ncc, lcc, n, x, y, lwk, iwk, list, lptr, lend, ier ) !*****************************************************************************80 ! !! ADDCST adds constraint curves to a Delaunay triangulation. ! ! Discussion: ! ! This subroutine provides for creation of a constrained ! Delaunay triangulation which, in some sense, covers an ! arbitrary connected region R rather than the convex hull ! of the nodes. This is achieved simply by forcing the ! presence of certain adjacencies (triangulation arcs) ! corresponding to constraint curves. The union of triangles ! coincides with the convex hull of the nodes, but triangles ! in R can be distinguished from those outside of R. The ! only modification required to generalize the definition of ! the Delaunay triangulation is replacement of property 5 ! (refer to TRMESH) by the following: ! ! 5') If a node is contained in the interior of the ! circumcircle of a triangle, then every interior point ! of the triangle is separated from the node by a ! constraint arc. ! ! In order to be explicit, we make the following definitions. ! A constraint region is the open interior of a ! simple closed positively oriented polygonal curve defined ! by an ordered sequence of three or more distinct nodes ! (constraint nodes) P(1),P(2),...,P(K), such that P(I) is ! adjacent to P(I+1) for I = 1,...,K with P(K+1) = P(1). ! Thus, the constraint region is on the left (and may have ! nonfinite area) as the sequence of constraint nodes is ! traversed in the specified order. The constraint regions ! must not contain nodes and must not overlap. The region ! R is the convex hull of the nodes with constraint regions ! excluded. ! ! Note that the terms boundary node and boundary arc are ! reserved for nodes and arcs on the boundary of the convex ! hull of the nodes. ! ! The algorithm is as follows: given a triangulation ! which includes one or more sets of constraint nodes, the ! corresponding adjacencies (constraint arcs) are forced to ! be present (Subroutine EDGE). Any additional new arcs ! required are chosen to be locally optimal (satisfy the ! modified circumcircle property). ! ! Modified: ! ! 16 June 2007 ! ! Author: ! ! Robert Renka, ! Department of Computer Science, ! University of North Texas, ! renka@cs.unt.edu ! ! Reference: ! ! Robert Renka, ! Algorithm 751: TRIPACK, ! A Constrained Two-Dimensional Delaunay Triangulation Package, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 1, 1996. ! ! Parameters: ! ! Input, integer NCC, the number of constraint curves (constraint ! regions). 0 <= NCC. ! ! Input, integer LCC(NCC) (or dummy array of length 1 if NCC = 0) ! containing the index (for X, Y, and LEND) of the first node of ! constraint I in LCC(I) for I = 1 to NCC. Thus, constraint I ! contains K = LCC(I+1) - LCC(I) nodes, K >= 3, stored in (X,Y) ! locations LCC(I), ..., LCC(I+1)-1, where LCC(NCC+1) = N+1. ! ! Input, integer N, the number of nodes in the triangulation, including ! constraint nodes. 3 <= N. ! ! Input, real ( kind = 8 ) X(N), Y(N), the coordinates of the nodes with ! non-constraint nodes in the first LCC(1)-1 locations, followed by NCC ! sequences of constraint nodes. Only one of these sequences may be ! specified in clockwise order to represent an exterior constraint curve (a ! constraint region with nonfinite area). ! ! Input/output, integer LWK. On input, the length of IWK. This must be ! at least 2*NI where NI is the maximum number of arcs which intersect a ! constraint arc to be added. NI is bounded by N-3. On output, the ! required length of IWK unless IER = 1 or IER = 3. In the case of ! IER = 1, LWK is not altered from its input value. ! ! Output, integer IWK(LWK), the endpoint indexes of the new arcs which ! were swapped in by the last call to subroutine EDGE. ! ! Input/output, integer LIST(*), LPTR(*), LEND(N), the data structure ! defining the triangulation. Refer to subroutine TRMESH. On output, ! The structure has all constraint arcs present unless IER /= 0. ! These arrays are not altered if IER = 1. ! ! Output, integer IER = Error indicator: ! 0 if no errors were encountered. ! 1 if NCC, N, or an LCC entry is outside its valid range, or ! LWK < 0 on input. ! 2 if more space is required in IWK. ! 3 if the triangulation data structure is invalid, or failure (in ! EDGE or OPTIM) was caused by collinear nodes on the convex hull ! boundary. An error message is written to logical unit 6 in ! this case. ! 4 if intersecting constraint arcs were encountered. ! 5 if a constraint region contains a node. ! implicit none integer lwk integer n integer i integer ier integer ifrst integer ilast integer iwk(lwk) integer k integer kbak integer kfor integer kn integer lcc(*) integer lccip1 integer lend(n) integer list(*) integer lp integer lpb integer lpf integer lpl integer lptr(*) integer lw integer lwd2 integer n1 integer n2 integer ncc real ( kind = 8 ) x(n) real ( kind = 8 ) y(n) lwd2 = lwk / 2 ! ! Test for errors in input parameters. ! ier = 0 if ( ncc < 0 .or. lwk < 0 ) then ier = 1 return end if if ( ncc == 0 ) then if ( n < 3 ) then ier = 1 else ier = 0 lwk = 0 end if return else lccip1 = n + 1 do i = ncc, 1, -1 if ( lccip1 - lcc(i) < 3 ) then ier = 1 return end if lccip1 = lcc(i) end do if ( lccip1 < 1 ) then ier = 1 return end if end if ! ! Force the presence of constraint arcs. The outer loop is ! on constraints in reverse order. IFRST and ILAST are ! the first and last nodes of constraint I. ! lwk = 0 ifrst = n + 1 do i = ncc, 1, -1 ilast = ifrst - 1 ifrst = lcc(i) ! ! Inner loop on constraint arcs N1-N2 in constraint I. ! n1 = ilast do n2 = ifrst, ilast lw = lwd2 call edge ( n1, n2, x, y, lw, iwk, list, lptr, lend, ier ) lwk = max ( lwk, 2 * lw ) if ( ier == 4 ) then ier = 3 end if if ( ier /= 0 ) then return end if n1 = n2 end do end do ! ! Test for errors. The outer loop is on constraint I with ! first and last nodes IFRST and ILAST, and the inner loop ! is on constraint nodes K with (KBAK,K,KFOR) a subsequence ! of constraint I. ! ifrst = n + 1 do i = ncc, 1, -1 ilast = ifrst - 1 ifrst = lcc(i) kbak = ilast do k = ifrst,ilast kfor = k + 1 if ( k == ilast ) then kfor = ifrst end if ! ! Find the LIST pointers LPF and LPB of KFOR and KBAK as neighbors of K. ! lpf = 0 lpb = 0 lpl = lend(k) lp = lpl do lp = lptr(lp) kn = abs ( list(lp) ) if ( kn == kfor ) then lpf = lp end if if ( kn == kbak ) then lpb = lp end if if ( lp == lpl ) then exit end if end do ! ! A pair of intersecting constraint arcs was encountered ! if and only if a constraint arc is missing (introduction ! of the second caused the first to be swapped out). ! if ( lpf == 0 .or. lpb == 0 ) then ier = 4 return end if ! ! Loop on neighbors KN of node K which follow KFOR and ! precede KBAK. The constraint region contains no nodes ! if and only if all such nodes KN are in constraint I. ! lp = lpf do lp = lptr(lp) if ( lp == lpb ) then exit end if kn = abs ( list(lp) ) if ( kn < ifrst .or. ilast < kn ) then ier = 5 return end if end do kbak = k end do end do ier = 0 return end subroutine addnod ( k, xk, yk, ist, ncc, lcc, n, x, y, list, lptr, lend, lnew, & ier ) !*****************************************************************************80 ! !! ADDNOD adds a node to a triangulation. ! ! Discussion: ! ! Given a triangulation of N nodes in the plane created by ! subroutine TRMESH or TRMSHR, this subroutine updates the ! data structure with the addition of a new node in position ! K. If node K is inserted into X and Y (K <= N) rather ! than appended (K = N+1), then a corresponding insertion ! must be performed in any additional arrays associated ! with the nodes. For example, an array of data values Z ! must be shifted down to open up position K for the new ! value: set Z(I+1) to Z(I) for I = N,N-1,...,K. For ! optimal efficiency, new nodes should be appended whenever ! possible. Insertion is necessary, however, to add a non- ! constraint node when constraints are present (refer to ! subroutine ADDCST). ! ! Note that a constraint node cannot be added by this ! routine. In order to insert a constraint node, it is ! necessary to add the node with no constraints present ! (call this routine with NCC = 0), update LCC by increment- ! ing the appropriate entries, and then create (or restore) ! the constraints by a call to ADDCST. ! ! The algorithm consists of the following steps: node K ! is located relative to the triangulation (TRFIND), its ! index is added to the data structure (INTADD or BDYADD), ! and a sequence of swaps (SWPTST and SWAP) are applied to ! the arcs opposite K so that all arcs incident on node K ! and opposite node K (excluding constraint arcs) are local- ! ly optimal (satisfy the circumcircle test). Thus, if a ! (constrained) Delaunay triangulation is input, a (con- ! strained) Delaunay triangulation will result. All indexes ! are incremented as necessary for an insertion. ! ! Modified: ! ! 29 March 2002 ! ! Author: ! ! Robert Renka, ! Department of Computer Science, ! University of North Texas, ! renka@cs.unt.edu ! ! Reference: ! ! Robert Renka, ! Algorithm 751: TRIPACK, ! A Constrained Two-Dimensional Delaunay Triangulation Package, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 1, 1996. ! ! Parameters: ! ! Input, integer K, the nodal index (index for X, Y, and LEND) of the ! new node to be added. 1 <= K <= LCC(1). (K <= N+1 if NCC=0). ! ! Input, real ( kind = 8 ) XK, YK, the coordinates of the new node (to be ! stored in X(K) and Y(K)). The node must not lie in a constraint region. ! ! Input, integer IST, the index of a node at which TRFIND begins the ! search. Search time depends on the proximity ! of this node to node K. 1 <= IST <= N. ! ! Input, integer NCC, the number of constraint curves. 0 <= NCC. ! ! Input/output, integer LCC(*), list of constraint curve starting indexes ! (or dummy array of length 1 if NCC = 0). Refer to subroutine ADDCST. ! On output, starting indexes incremented by 1 to reflect the insertion ! of K unless NCC = 0 or (IER /= 0 and IER /= -4). ! ! Input/output, integer N, the number of nodes in the triangulation. ! 3 <= N. Note that N will be incremented following the addition of node K. ! ! Input, real ( kind = 8 ) X(N+1), real Y(N+1), containing the coordinates ! of the nodes in the first N positions with non-constraint nodes ! in the first LCC(1)-1 locations if 0 < NCC. On output, updated with ! the insertion of XK and YK in the K-th positions (node I+1 was node ! I before the insertion for I = K to N if K <= N) ! unless IER /= 0 and IER /= -4. ! ! Input/output, integer LIST(*), LPTR(*), LEND(N), LNEW, the data ! structure associated with the triangulation of nodes 1 to N. The ! arrays must have sufficient length for N+1 nodes. Refer to TRMESH. ! On output, updated with the addition of node K unless ! IER /= 0 and IER /= -4. ! ! Output, integer IER = Error indicator: ! 0 if no errors were encountered. ! -1 if K, IST, NCC, N, or an LCC entry is outside its valid range on input. ! -2 if all nodes (including K) are collinear. ! L if nodes L and K coincide for some L. ! -3 if K lies in a constraint region. ! -4 if an error flag is returned by SWAP implying that the triangulation ! (geometry) was bad on input. ! implicit none logical crtri integer i integer i1 integer i2 integer i3 integer ibk integer ier integer in1 integer indxcc integer io1 integer io2 integer ist integer k integer kk integer l integer lcc(*) integer lccip1 integer lend(*) integer list(*) integer lnew integer lp integer lpf integer lpo1 integer lptr(*) integer lstptr integer n integer ncc integer nm1 logical swptst real ( kind = 8 ) x(*) real ( kind = 8 ) xk real ( kind = 8 ) y(*) real ( kind = 8 ) yk kk = k ! ! Test for an invalid input parameter. ! if ( kk < 1 .or. ist < 1 .or. n < ist & .or. ncc < 0 .or. n < 3 ) then ier = -1 return end if lccip1 = n + 1 do i = ncc, 1, -1 if ( lccip1-lcc(i) < 3 ) then ier = -1 return end if lccip1 = lcc(i) end do if ( lccip1 < kk ) then ier = -1 return end if ! ! Find a triangle (I1,I2,I3) containing K or the rightmost ! (I1) and leftmost (I2) visible boundary nodes as viewed from node K. ! call trfind ( ist, xk, yk, n, x, y, list, lptr, lend, i1, i2, i3 ) ! ! Test for collinear nodes, duplicate nodes, and K lying in ! a constraint region. ! if ( i1 == 0 ) then ier = -2 return end if if ( i3 /= 0 ) then l = i1 if ( xk == x(l) .and. yk == y(l) ) then ier = l return end if l = i2 if ( xk == x(l) .and. yk == y(l) ) then ier = l return end if l = i3 if ( xk == x(l) .and. yk == y(l) ) then ier = l return end if if ( 0 < ncc .and. crtri(ncc,lcc,i1,i2,i3) ) then ier = -3 return end if else ! ! K is outside the convex hull of the nodes and lies in a ! constraint region iff an exterior constraint curve is present. ! if ( 0 < ncc .and. indxcc(ncc,lcc,n,list,lend) /= 0 ) then ier = -3 return end if end if ! ! No errors encountered. ! ier = 0 nm1 = n n = n + 1 if (kk < n) then ! ! Open a slot for K in X, Y, and LEND, and increment all ! nodal indexes which are greater than or equal to K. ! ! Note that LIST, LPTR, and LNEW are not yet updated with ! either the neighbors of K or the edges terminating on K. ! do ibk = nm1, kk, -1 x(ibk+1) = x(ibk) y(ibk+1) = y(ibk) lend(ibk+1) = lend(ibk) end do do i = 1, ncc lcc(i) = lcc(i) + 1 end do l = lnew - 1 do i = 1, l if ( kk <= list(i) ) then list(i) = list(i) + 1 end if if ( list(i) <= -kk ) then list(i) = list(i) - 1 end if end do if ( kk <= i1 ) then i1 = i1 + 1 end if if ( kk <= i2 ) then i2 = i2 + 1 end if if ( kk <= i3 ) then i3 = i3 + 1 end if end if ! ! Insert K into X and Y, and update LIST, LPTR, LEND, and ! LNEW with the arcs containing node K. ! x(kk) = xk y(kk) = yk if ( i3 == 0 ) then call bdyadd ( kk, i1, i2, list, lptr, lend, lnew ) else call intadd ( kk, i1, i2, i3, list, lptr, lend, lnew ) end if ! ! Initialize variables for optimization of the triangulation. ! lp = lend(kk) lpf = lptr(lp) io2 = list(lpf) lpo1 = lptr(lpf) io1 = abs ( list(lpo1) ) ! ! Begin loop: find the node opposite K. ! do lp = lstptr ( lend(io1), io2, list, lptr ) if ( 0 <= list(lp) ) then lp = lptr(lp) in1 = abs ( list(lp) ) ! ! Swap test: if a swap occurs, two new arcs are ! opposite K and must be tested. ! if ( .not. crtri ( ncc, lcc, io1, io2, in1 ) ) then if ( swptst(in1,kk,io1,io2,x,y) ) then call swap ( in1, kk, io1, io2, list, lptr, lend, lpo1 ) if ( lpo1 == 0 ) then ier = -4 exit end if io1 = in1 cycle end if end if end if ! ! No swap occurred. Test for termination and reset IO2 and IO1. ! if ( lpo1 == lpf .or. list(lpo1) < 0 ) then exit end if io2 = io1 lpo1 = lptr(lpo1) io1 = abs ( list(lpo1) ) end do return end function areap ( x, y, nb, nodes ) !*****************************************************************************80 ! !! AREAP computes the signed area of a polygonal curve. ! ! Discussion: ! ! Given a sequence of NB points in the plane, this function ! computes the signed area bounded by the closed polygonal ! curve which passes through the points in the ! specified order. Each simple closed curve is positively ! oriented (bounds positive area) if and only if the points ! are specified in counterclockwise order. The last point ! of the curve is taken to be the first point specified, and ! this point should therefore not be specified twice. ! ! The area of a triangulation may be computed by calling ! AREAP with values of NB and NODES determined by subroutine ! BNODES. ! ! Modified: ! ! 16 June 2007 ! ! Author: ! ! Robert Renka, ! Department of Computer Science, ! University of North Texas, ! renka@cs.unt.edu ! ! Reference: ! ! Robert Renka, ! Algorithm 751: TRIPACK, ! A Constrained Two-Dimensional Delaunay Triangulation Package, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 1, 1996. ! ! Parameters: ! ! Input, real ( kind = 8 ) X(*), Y(*), the Cartesian coordinates of a set ! of points. ! ! Input, integer NB, the number of points in the curve. ! ! Input, integer NODES(NB), the indices of the points that ! make up the closed curve. ! ! Output, real ( kind = 8 ) AREAP, the signed area bounded by the curve. ! implicit none integer nb real ( kind = 8 ) a real ( kind = 8 ) areap integer i integer nd1 integer nd2 integer nodes(nb) real ( kind = 8 ) x(*) real ( kind = 8 ) y(*) a = 0.0D+00 if ( nb < 3 ) then areap = 0.0D+00 return end if nd2 = nodes(nb) ! ! Loop on line segments NODES(I-1) -> NODES(I), where ! NODES(0) = NODES(NB), adding twice the signed trapezoid ! areas (integrals of the linear interpolants) to A. ! do i = 1, nb nd1 = nd2 nd2 = nodes(i) a = a + ( x(nd2) - x(nd1) ) * ( y(nd1) + y(nd2) ) end do ! ! A contains twice the negative signed area of the region. ! areap = -a / 2.0D+00 return end subroutine bdyadd ( kk, i1, i2, list, lptr, lend, lnew ) !*****************************************************************************80 ! !! BDYADD adds a boundary node to a triangulation. ! ! Discussion: ! ! This subroutine adds a boundary node to a triangulation ! of a set of points in the plane. The data structure is ! updated with the insertion of node KK, but no optimization ! is performed. ! ! Modified: ! ! 16 June 2007 ! ! Author: ! ! Robert Renka, ! Department of Computer Science, ! University of North Texas, ! renka@cs.unt.edu ! ! Reference: ! ! Robert Renka, ! Algorithm 751: TRIPACK, ! A Constrained Two-Dimensional Delaunay Triangulation Package, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 1, 1996. ! ! Parameters: ! ! Input, integer KK, the index of a node to be connected to the sequence ! of all visible boundary nodes. 1 <= KK and ! KK must not be equal to I1 or I2. ! ! Input, integer I1, the first (rightmost as viewed from KK) boundary ! node in the triangulation which is visible from ! node KK (the line segment KK-I1 intersects no arcs. ! ! Input, integer I2, the last (leftmost) boundary node which is visible ! from node KK. I1 and I2 may be determined by subroutine TRFIND. ! ! Input/output, integer LIST(*), LPTR(*), LEND(N), LNEW. The ! triangulation data structure created by TRMESH or TRMSHR. ! On input, nodes I1 and I2 must be included in the triangulation. ! On output, the data structure has been updated with the addition ! of node KK. Node KK is connected to I1, I2, and all boundary ! nodes in between. ! implicit none integer i1 integer i2 integer k integer kk integer lend(*) integer list(*) integer lnew integer lp integer lptr(*) integer lsav integer n1 integer n2 integer next integer nsav k = kk n1 = i1 n2 = i2 ! ! Add K as the last neighbor of N1. ! lp = lend(n1) lsav = lptr(lp) lptr(lp) = lnew list(lnew) = -k lptr(lnew) = lsav lend(n1) = lnew lnew = lnew + 1 next = -list(lp) list(lp) = next nsav = next ! ! Loop on the remaining boundary nodes between N1 and N2, ! adding K as the first neighbor. ! do lp = lend(next) call insert ( k, lp, list, lptr, lnew ) if ( next == n2 ) then exit end if next = -list(lp) list(lp) = next end do ! ! Add the boundary nodes between N1 and N2 as neighbors ! of node K. ! lsav = lnew list(lnew) = n1 lptr(lnew) = lnew + 1 lnew = lnew + 1 next = nsav do if ( next == n2 ) then exit end if list(lnew) = next lptr(lnew) = lnew + 1 lnew = lnew + 1 lp = lend(next) next = list(lp) end do list(lnew) = -n2 lptr(lnew) = lsav lend(k) = lnew lnew = lnew + 1 return end subroutine bnodes ( n, list, lptr, lend, nodes, nb, na, nt ) !*****************************************************************************80 ! !! BNODES returns a list of the boundary nodes. ! ! Discussion: ! ! Given a triangulation of N points in the plane, this ! subroutine returns an array containing the indexes, in ! counterclockwise order, of the nodes on the boundary of ! the convex hull of the set of points. ! ! Modified: ! ! 16 June 2007 ! ! Author: ! ! Robert Renka, ! Department of Computer Science, ! University of North Texas, ! renka@cs.unt.edu ! ! Reference: ! ! Robert Renka, ! Algorithm 751: TRIPACK, ! A Constrained Two-Dimensional Delaunay Triangulation Package, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 1, 1996. ! ! Parameters: ! ! Input, integer N, the number of nodes in the triangulation. 3 <= N. ! ! Input, integer LIST(*), LPTR(*), LEND(N), the data structure defining ! the triangulation. Refer to subroutine TRMESH. ! ! Output, integer NODES(NB), ordered sequence of boundary node indexes ! in the range 1 to N. ! ! Output, integer NB, the number of boundary nodes. ! ! Output, integer NA, NT, the number of arcs and triangles, respectively, ! in the triangulation. ! implicit none integer n integer k integer lend(n) integer list(*) integer lp integer lptr(*) integer n0 integer na integer nb integer nodes(*) integer nst integer nt ! ! Set NST to the first boundary node encountered. ! nst = 1 do lp = lend(nst) if ( list(lp) < 0 ) then exit end if nst = nst + 1 end do ! ! Initialization. ! nodes(1) = nst k = 1 n0 = nst ! ! Traverse the boundary in counterclockwise order. ! do lp = lend(n0) lp = lptr(lp) n0 = list(lp) if ( n0 == nst ) then exit end if k = k + 1 nodes(k) = n0 end do ! ! Termination. ! nb = k nt = 2 * n - nb - 2 na = nt + n - 1 return end subroutine circum ( x1, y1, x2, y2, x3, y3, ratio, xc, yc, cr, sa, ar ) !*****************************************************************************80 ! !! CIRCUM determines the circumcenter (and more) of a triangle. ! ! Discussion: ! ! Given three vertices defining a triangle, this routine ! returns the circumcenter, circumradius, signed ! triangle area, and, optionally, the aspect ratio of the ! triangle. ! ! Modified: ! ! 16 June 2007 ! ! Author: ! ! Robert Renka, ! Department of Computer Science, ! University of North Texas, ! renka@cs.unt.edu ! ! Reference: ! ! Robert Renka, ! Algorithm 751: TRIPACK, ! A Constrained Two-Dimensional Delaunay Triangulation Package, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 1, 1996. ! ! Parameters: ! ! Input, real ( kind = 8 ) X1, Y1, X2, Y2, X3, Y3, the coordinates of ! the vertices. ! ! Input, logical RATIO, is TRUE if and only if the aspect ratio is ! to be computed. ! ! Output, real ( kind = 8 ) XC, YC, coordinates of the circumcenter (center ! of the circle defined by the three points) unless SA = 0, in which XC ! and YC are not altered. ! ! Output, real ( kind = 8 ) CR, the circumradius (radius of the circle ! defined by the three points) unless SA = 0 (infinite radius), in which ! case CR is not altered. ! ! Output, real ( kind = 8 ) SA, the signed triangle area with positive value ! if and only if the vertices are specified in counterclockwise order: ! (X3,Y3) is strictly to the left of the directed line from (X1,Y1) ! toward (X2,Y2). ! ! Output, real ( kind = 8 ) AR, the aspect ratio r/CR, where r is the ! radius of the inscribed circle, unless RATIO = FALSE, in which case AR ! is not altered. AR is in the range 0 to 0.5, with value 0 iff SA = 0 and ! value 0.5 iff the vertices define an equilateral triangle. ! implicit none real ( kind = 8 ) ar real ( kind = 8 ) cr real ( kind = 8 ) ds(3) real ( kind = 8 ) fx real ( kind = 8 ) fy logical ratio real ( kind = 8 ) sa real ( kind = 8 ) u(3) real ( kind = 8 ) v(3) real ( kind = 8 ) x1 real ( kind = 8 ) x2 real ( kind = 8 ) x3 real ( kind = 8 ) xc real ( kind = 8 ) y1 real ( kind = 8 ) y2 real ( kind = 8 ) y3 real ( kind = 8 ) yc ! ! Set U(K) and V(K) to the x and y components, respectively, ! of the directed edge opposite vertex K. ! u(1) = x3 - x2 u(2) = x1 - x3 u(3) = x2 - x1 v(1) = y3 - y2 v(2) = y1 - y3 v(3) = y2 - y1 ! ! Set SA to the signed triangle area. ! sa = ( u(1) * v(2) - u(2) * v(1) ) / 2.0D+00 if ( sa == 0.0D+00 ) then if ( ratio ) then ar = 0.0D+00 end if return end if ! ! Set DS(K) to the squared distance from the origin to vertex K. ! ds(1) = x1 * x1 + y1 * y1 ds(2) = x2 * x2 + y2 * y2 ds(3) = x3 * x3 + y3 * y3 ! ! Compute factors of XC and YC. ! fx = - dot_product ( ds(1:3), v(1:3) ) fy = dot_product ( ds(1:3), u(1:3) ) xc = fx / ( 4.0D+00 * sa ) yc = fy / ( 4.0D+00 * sa ) cr = sqrt ( ( xc - x1 )**2 + ( yc - y1 )**2 ) if ( .not. ratio ) then return end if ! ! Compute the squared edge lengths and aspect ratio. ! ds(1:3) = u(1:3)**2 + v(1:3)**2 ar = 2.0D+00 * abs ( sa ) / & ( ( sqrt ( ds(1) ) + sqrt ( ds(2) ) + sqrt ( ds(3) ) ) * cr ) return end function crtri ( ncc, lcc, i1, i2, i3 ) !*****************************************************************************80 ! !! CRTRI determines if a triangle lies in a constraint region. ! ! Discussion: ! ! This function returns TRUE if and only if triangle (I1, ! I2,I3) lies in a constraint region. ! ! Modified: ! ! 16 June 2007 ! ! Author: ! ! Robert Renka, ! Department of Computer Science, ! University of North Texas, ! renka@cs.unt.edu ! ! Reference: ! ! Robert Renka, ! Algorithm 751: TRIPACK, ! A Constrained Two-Dimensional Delaunay Triangulation Package, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 1, 1996. ! ! Parameters: ! ! Input, integer NCC, LCC(*), onstraint data structure. Refer to ADDCST. ! ! Input, integer I1, I2, I3, Nodal indexes of the counterclockwise- ! ordered vertices of a triangle. ! ! Output, logical CRTRI, is TRUE if and only if (I1,I2,I3) is a ! constraint region triangle. ! implicit none logical crtri integer i integer i1 integer i2 integer i3 integer imax integer imin integer lcc(*) integer ncc imax = max ( i1, i2, i3 ) ! ! Find the index I of the constraint containing IMAX. ! i = ncc + 1 do i = i - 1 if ( i <= 0 ) then crtri = .false. return end if if ( lcc(i) <= imax ) then exit end if end do imin = min ( i1, i2, i3 ) ! ! P lies in a constraint region iff I1, I2, and I3 are nodes ! of the same constraint (LCC(I) <= IMIN), and (IMIN,IMAX) ! is (I1,I3), (I2,I1), or (I3,I2). ! crtri = lcc(i) <= imin .and. ( & ( imin == i1 .and. imax == i3 ) .or. & ( imin == i2 .and. imax == i1 ) .or. & ( imin == i3 .and. imax == i2 ) ) return end subroutine delarc ( n, io1, io2, list, lptr, lend, lnew, ier ) !*****************************************************************************80 ! !! DELARC deletes a boundary arc from a triangulation. ! ! Discussion: ! ! This subroutine deletes a boundary arc from a triangula- ! tion. It may be used to remove a null triangle from the ! convex hull boundary. Note, however, that if the union of ! triangles is rendered nonconvex, Subroutines DELNOD, EDGE, ! and TRFIND may fail. Thus, Subroutines ADDCST, ADDNOD, ! DELNOD, EDGE, and NEARND should not be called following ! an arc deletion. ! ! Modified: ! ! 16 June 2007 ! ! Author: ! ! Robert Renka, ! Department of Computer Science, ! University of North Texas, ! renka@cs.unt.edu ! ! Reference: ! ! Robert Renka, ! Algorithm 751: TRIPACK, ! A Constrained Two-Dimensional Delaunay Triangulation Package, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 1, 1996. ! ! Parameters: ! ! Input, integer N, the number of nodes in the triangulation. 4 <= N. ! ! Input, integer IO1, IO2, the indexes (in the range 1 to N) of a pair of ! adjacent boundary nodes defining the arc to be removed. ! ! Input/output, integer LIST(*), LPTR(*), LEND(N), LNEW, the ! triangulation data structure created by TRMESH or TRMSHR. ! On output, updated with the removal of arc IO1-IO2 unless 0 < IER. ! ! Output, integer IER = Error indicator: ! 0, if no errors were encountered. ! 1, if N, IO1, or IO2 is outside its valid range, or IO1 = IO2. ! 2, if IO1-IO2 is not a boundary arc. ! 3, if the node opposite IO1-IO2 is already a boundary node, and ! thus IO1 or IO2 has only two neighbors or a deletion would result ! in two triangulations sharing a single node. ! 4, if one of the nodes is a neighbor of the other, but not vice ! versa, implying an invalid triangulation data structure. ! implicit none integer n integer ier integer io1 integer io2 integer lend(n) integer list(*) integer lnew integer lp integer lph integer lpl integer lptr(*) integer lstptr integer n1 integer n2 integer n3 n1 = io1 n2 = io2 ! ! Test for errors, and set N1->N2 to the directed boundary ! edge associated with IO1-IO2: (N1,N2,N3) is a triangle ! for some N3. ! if ( n < 4 .or. n1 < 1 .or. n < n1 .or. & n2 < 1 .or. n < n2 .or. n1 == n2 ) then ier = 1 return end if lpl = lend(n2) if ( -list(lpl) /= n1 ) then n1 = n2 n2 = io1 lpl = lend(n2) if ( -list(lpl) /= n1 ) then ier = 2 return end if end if ! ! Set N3 to the node opposite N1->N2 (the second neighbor ! of N1), and test for error 3 (N3 already a boundary node). ! lpl = lend(n1) lp = lptr(lpl) lp = lptr(lp) n3 = abs ( list(lp) ) lpl = lend(n3) if ( list(lpl) <= 0 ) then ier = 3 return end if ! ! Delete N2 as a neighbor of N1, making N3 the first ! neighbor, and test for error 4 (N2 not a neighbor ! of N1). Note that previously computed pointers may ! no longer be valid following the call to DELNB. ! call delnb ( n1, n2, n, list, lptr, lend, lnew, lph ) if ( lph < 0 ) then ier = 4 return end if ! ! Delete N1 as a neighbor of N2, making N3 the new last neighbor. ! call delnb ( n2, n1, n, list, lptr, lend, lnew, lph ) ! ! Make N3 a boundary node with first neighbor N2 and last neighbor N1. ! lp = lstptr ( lend(n3), n1, list, lptr ) lend(n3) = lp list(lp) = -n1 ! ! No errors encountered. ! ier = 0 return end subroutine delnb ( n0, nb, n, list, lptr, lend, lnew, lph ) !*****************************************************************************80 ! !! DELNB deletes a neighbor from an adjacency list. ! ! Discussion: ! ! This subroutine deletes a neighbor NB from the adjacency ! list of node N0 (but N0 is not deleted from the adjacency ! list of NB) and, if NB is a boundary node, makes N0 a ! boundary node. For pointer (LIST index) LPH to NB as a ! neighbor of N0, the empty LIST,LPTR location LPH is filled ! in with the values at LNEW-1, pointer LNEW-1 (in LPTR and ! possibly in LEND) is changed to LPH, and LNEW is decremented ! This requires a search of LEND and LPTR entailing an ! expected operation count of O(N). ! ! Modified: ! ! 16 June 2007 ! ! Author: ! ! Robert Renka, ! Department of Computer Science, ! University of North Texas, ! renka@cs.unt.edu ! ! Reference: ! ! Robert Renka, ! Algorithm 751: TRIPACK, ! A Constrained Two-Dimensional Delaunay Triangulation Package, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 1, 1996. ! ! Parameters: ! ! Input, integer N0, NB, indexes, in the range 1 to N, of a pair of ! nodes such that NB is a neighbor of N0. ! (N0 need not be a neighbor of NB.) ! ! Input, integer N, the number of nodes in the triangulation. 3 <= N. ! ! Input/output, integer LIST(*), LPTR(*), LEND(N), LNEW, the data ! structure defining the triangulation. On output, updated with ! the removal of NB from the adjacency list of N0 unless LPH < 0. ! ! Output, integer LPH, list pointer to the hole (NB as a neighbor of ! N0) filled in by the values at LNEW-1 or error indicator: ! >0, if no errors were encountered. ! -1, if N0, NB, or N is outside its valid range. ! -2, if NB is not a neighbor of N0. ! ! Local parameters: ! ! I = DO-loop index ! LNW = LNEW-1 (output value of LNEW) ! LP = LIST pointer of the last neighbor of NB ! LPB = Pointer to NB as a neighbor of N0 ! LPL = Pointer to the last neighbor of N0 ! LPP = Pointer to the neighbor of N0 that precedes NB ! NN = Local copy of N ! implicit none integer n integer i integer lend(n) integer list(*) integer lnew integer lnw integer lp integer lpb integer lph integer lpl integer lpp integer lptr(*) integer n0 integer nb integer nn nn = n ! ! Test for error 1. ! if ( n0 < 1 .or. nn < n0 .or. nb < 1 .or. & nn < nb .or. nn < 3 ) then lph = -1 return end if ! ! Find pointers to neighbors of N0: ! ! LPL points to the last neighbor, ! LPP points to the neighbor NP preceding NB, and ! LPB points to NB. ! lpl = lend(n0) lpp = lpl lpb = lptr(lpp) do if ( list(lpb) == nb ) then go to 2 end if lpp = lpb lpb = lptr(lpp) if ( lpb == lpl ) then exit end if end do ! ! Test for error 2 (NB not found). ! if ( abs ( list(lpb) ) /= nb ) then lph = -2 return end if ! ! NB is the last neighbor of N0. Make NP the new last ! neighbor and, if NB is a boundary node, then make N0 ! a boundary node. ! lend(n0) = lpp lp = lend(nb) if ( list(lp) < 0 ) then list(lpp) = -list(lpp) end if go to 3 ! ! NB is not the last neighbor of N0. If NB is a boundary ! node and N0 is not, then make N0 a boundary node with ! last neighbor NP. ! 2 continue lp = lend(nb) if ( list(lp) < 0 .and. 0 < list(lpl) ) then lend(n0) = lpp list(lpp) = -list(lpp) end if ! ! Update LPTR so that the neighbor following NB now fol- ! lows NP, and fill in the hole at location LPB. ! 3 continue lptr(lpp) = lptr(lpb) lnw = lnew - 1 list(lpb) = list(lnw) lptr(lpb) = lptr(lnw) do i = nn, 1, -1 if (lend(i) == lnw) then lend(i) = lpb exit end if end do do i = 1, lnw-1 if (lptr(i) == lnw) then lptr(i) = lpb end if end do ! ! No errors encountered. ! lnew = lnw lph = lpb return end subroutine delnod ( k, ncc, lcc, n, x, y, list, lptr, lend, lnew, lwk, iwk, & ier ) !*****************************************************************************80 ! !! DELNOD deletes a node from a triangulation. ! ! Discussion: ! ! This subroutine deletes node K (along with all arcs ! incident on node K) from a triangulation of N nodes in the ! plane, and inserts arcs as necessary to produce a triangu- ! lation of the remaining N-1 nodes. If a Delaunay triangu- ! lation is input, a Delaunay triangulation will result, and ! thus, DELNOD reverses the effect of a call to subroutine ! ADDNOD. ! ! Note that a constraint node cannot be deleted by this ! routine. In order to delete a constraint node, it is ! necessary to call this routine with NCC = 0, decrement the ! appropriate LCC entries (LCC(I) such that K < LCC(I) ), and ! then create (or restore) the constraints by a call to sub- ! routine ADDCST. ! ! Note that the deletion may result in all remaining nodes ! being collinear. This situation is not flagged. ! ! Modified: ! ! 16 June 2007 ! ! Author: ! ! Robert Renka, ! Department of Computer Science, ! University of North Texas, ! renka@cs.unt.edu ! ! Reference: ! ! Robert Renka, ! Algorithm 751: TRIPACK, ! A Constrained Two-Dimensional Delaunay Triangulation Package, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 1, 1996. ! ! Parameters: ! ! Input, integer K, the index (for X and Y) of the node to be deleted. ! 1 <= K < LCC(1). (K <= N if NCC=0). ! ! Input, integer NCC, the number of constraint curves. 0 <= NCC. ! ! Input/output, integer LCC(*), list of constraint curve starting indexes ! (or dummy array of length 1 if NCC = 0). Refer to subroutine ADDCST. ! On output, decremented by 1 to reflect the deletion of K unless NCC = 0 ! or 1 <= IER <= 4. ! ! Input/output, integer N, the number of nodes in the triangulation. ! 4 <= N. Note that N will be decremented following the deletion. ! ! Input/output, real ( kind = 8 ) X(N), Y(N), the coordinates of the nodes ! with non-constraint nodes in the first LCC(1)-1 locations if 0 < NCC. ! On output, updated arrays of length N-1 containing nodal coordinates ! (with elements K+1,...,N shifted a position and thus overwriting ! element K) unless 1 <= IER <= 4. (N here denotes the input value.) ! ! Input/output, integer LIST(*), LPTR(*), LEND(N), LNEW, the data ! structure defining the triangulation. Refer to subroutine TRMESH. ! On output, updated to reflect the deletion unless IER /= 0. Note ! that the data structure may have been altered if 3 <= IER. ! ! Input/output, integer LWK. On input, the number of columns reserved ! for IWK. LWK must be at least NNB-3, where NNB is the number of ! neighbors of node K, including an extra pseudo-node if K is a ! boundary node. On output, the number of IWK columns required unless ! IER = 1 or IER = 3. ! ! Output, integer IWK(2,LWK), indexes of the endpoints of the new arcs added ! unless LWK = 0 or 1 <= IER <= 4. (Arcs are associated with columns, or ! pairs of adjacent elements if IWK is declared as a singly-subscripted ! array.) ! ! Output, integer IER = Error indicator: ! 0, if no errors were encountered. ! 1, if K, NCC, N, or an LCC entry is outside its valid range or LWK < 0 ! on input. ! 2, if more space is required in IWK. Refer to LWK. ! 3, if the triangulation data structure is invalid on input. ! 4, if K is an interior node with 4 or more neighbors, and the number of ! neighbors could not be reduced to 3 by swaps. This could be caused by ! floating point errors with collinear nodes or by an invalid data ! structure. ! 5, if an error flag was returned by OPTIM. An error message is written ! to the standard output unit in this event. ! implicit none logical bdry integer i integer ier integer ierr integer iwk(2,*) integer iwl integer j integer k integer lcc(*) integer lccip1 logical left integer lend(*) integer list(*) integer lnew integer lnw integer lp integer lp21 integer lpf integer lph integer lpl integer lpl2 integer lpn integer lptr(*) integer lstptr integer lwk integer lwkl integer n integer n1 integer n2 integer nbcnt integer ncc integer nfrst integer nit integer nl integer nn integer nnb integer nr real ( kind = 8 ) x(*) real ( kind = 8 ) x1 real ( kind = 8 ) x2 real ( kind = 8 ) xl real ( kind = 8 ) xr real ( kind = 8 ) y(*) real ( kind = 8 ) y1 real ( kind = 8 ) y2 real ( kind = 8 ) yl real ( kind = 8 ) yr ! ! Set N1 to K and NNB to the number of neighbors of N1 (plus ! one if N1 is a boundary node), and test for errors. LPF ! and LPL are LIST indexes of the first and last neighbors ! of N1, IWL is the number of IWK columns containing arcs, ! and BDRY is TRUE iff N1 is a boundary node. ! ier = 0 n1 = k nn = n if ( ncc < 0 .or. n1 < 1 .or. nn < 4 .or. lwk < 0 ) then ier = 1 return end if lccip1 = nn + 1 do i = ncc, 1, -1 if ( lccip1 - lcc(i) < 3 ) then ier = 1 return end if lccip1 = lcc(i) end do if ( lccip1 <= n1 ) then ier = 1 return end if lpl = lend(n1) lpf = lptr(lpl) nnb = nbcnt ( lpl, lptr ) bdry = list(lpl) < 0 if ( bdry ) then nnb = nnb + 1 end if if ( nnb < 3 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DELNOD - Fatal error!' write ( *, '(a)' ) ' NNB < 3.' write ( *, '(a,i6)' ) ' NNB = ', nnb ier = 3 return end if lwkl = lwk lwk = nnb - 3 if ( lwkl < lwk ) then ier = 2 return end if iwl = 0 if ( nnb == 3 ) then go to 5 end if ! ! Initialize for loop on arcs N1-N2 for neighbors N2 of N1, ! beginning with the second neighbor. NR and NL are the ! neighbors preceding and following N2, respectively, and ! LP indexes NL. The loop is exited when all possible ! swaps have been applied to arcs incident on N1. If N1 ! is interior, the number of neighbors will be reduced ! to 3. ! x1 = x(n1) y1 = y(n1) nfrst = list(lpf) nr = nfrst xr = x(nr) yr = y(nr) lp = lptr(lpf) n2 = list(lp) x2 = x(n2) y2 = y(n2) lp = lptr(lp) ! ! Top of loop: set NL to the neighbor following N2. ! 2 continue nl = abs ( list(lp) ) if ( nl == nfrst .and. bdry ) then go to 5 end if xl = x(nl) yl = y(nl) ! ! Test for a convex quadrilateral. To avoid an incorrect ! test caused by collinearity, use the fact that if N1 ! is a boundary node, then N1 LEFT NR->NL and if N2 is ! a boundary node, then N2 LEFT NL->NR. ! lpl2 = lend(n2) if ( (bdry .or. left(xr,yr,xl,yl,x1,y1)) .and. & (list(lpl2) < 0 .or. & left(xl,yl,xr,yr,x2,y2)) ) then go to 3 end if ! ! Nonconvex quadrilateral -- no swap is possible. ! nr = n2 xr = x2 yr = y2 go to 4 ! ! The quadrilateral defined by adjacent triangles ! (N1,N2,NL) and (N2,N1,NR) is convex. Swap in ! NL-NR and store it in IWK. Indexes larger than N1 ! must be decremented since N1 will be deleted from ! X and Y. ! 3 continue call swap ( nl, nr, n1, n2, list, lptr, lend, lp21 ) iwl = iwl + 1 if ( nl <= n1 ) then iwk(1,iwl) = nl else iwk(1,iwl) = nl - 1 end if if ( nr <= n1 ) then iwk(2,iwl) = nr else iwk(2,iwl) = nr - 1 end if ! ! Recompute the LIST indexes LPL,LP and decrement NNB. ! lpl = lend(n1) nnb = nnb - 1 if ( nnb == 3 ) then go to 5 end if lp = lstptr ( lpl, nl, list, lptr ) if ( nr == nfrst ) then go to 4 end if ! ! NR is not the first neighbor of N1. ! Back up and test N1-NR for a swap again: Set N2 to ! NR and NR to the previous neighbor of N1 -- the ! neighbor of NR which follows N1. LP21 points to NL ! as a neighbor of NR. ! n2 = nr x2 = xr y2 = yr lp21 = lptr(lp21) lp21 = lptr(lp21) nr = abs ( list(lp21) ) xr = x(nr) yr = y(nr) go to 2 ! ! Bottom of loop -- test for invalid termination. ! 4 continue if ( n2 == nfrst ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DELNOD - Fatal error!' write ( *, '(a)' ) ' N2 = NFRST.' ier = 4 return end if n2 = nl x2 = xl y2 = yl lp = lptr(lp) go to 2 ! ! Delete N1 from the adjacency list of N2 for all neighbors ! N2 of N1. LPL points to the last neighbor of N1. ! LNEW is stored in local variable LNW. ! 5 continue lp = lpl lnw = lnew ! ! Loop on neighbors N2 of N1, beginning with the first. ! 6 continue lp = lptr(lp) n2 = abs ( list(lp) ) call delnb ( n2, n1, n, list, lptr, lend, lnw, lph ) if ( lph < 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DELNOD - Fatal error!' write ( *, '(a)' ) ' LPH < 0.' ier = 3 return end if ! ! LP and LPL may require alteration. ! if ( lpl == lnw ) then lpl = lph end if if ( lp == lnw ) then lp = lph end if if ( lp /= lpl ) then go to 6 end if ! ! Delete N1 from X, Y, and LEND, and remove its adjacency ! list from LIST and LPTR. LIST entries (nodal indexes) ! which are larger than N1 must be decremented. ! nn = nn - 1 if ( nn < n1 ) then go to 9 end if do i = n1,nn x(i) = x(i+1) y(i) = y(i+1) lend(i) = lend(i+1) end do do i = 1,lnw-1 if (list(i) > n1) then list(i) = list(i) - 1 end if if (list(i) < -n1) then list(i) = list(i) + 1 end if end do ! ! For LPN = first to last neighbors of N1, delete the ! preceding neighbor (indexed by LP). ! ! Each empty LIST,LPTR location LP is filled in with the ! values at LNW-1, and LNW is decremented. All pointers ! (including those in LPTR and LEND) with value LNW-1 ! must be changed to LP. ! ! LPL points to the last neighbor of N1. ! 9 continue if ( bdry ) then nnb = nnb - 1 end if lpn = lpl do j = 1, nnb lnw = lnw - 1 lp = lpn lpn = lptr(lp) list(lp) = list(lnw) lptr(lp) = lptr(lnw) if ( lptr(lpn) == lnw ) then lptr(lpn) = lp end if if ( lpn == lnw ) then lpn = lp end if do i = nn,1,-1 if (lend(i) == lnw) then lend(i) = lp exit end if end do do i = lnw-1, 1, -1 if ( lptr(i) == lnw ) then lptr(i) = lp end if end do end do ! ! Decrement LCC entries. ! lcc(1:ncc) = lcc(1:ncc) - 1 ! ! Update N and LNEW, and optimize the patch of triangles ! containing K (on input) by applying swaps to the arcs in IWK. ! n = nn lnew = lnw if ( 0 < iwl ) then nit = 4 * iwl call optim ( x, y, iwl, list, lptr, lend, nit, iwk, ierr ) if ( ierr /= 0 ) then ier = 5 write (*,100) nit, ierr 100 format (//5x,'*** error in optim: nit = ',i4, & ', ier = ',i1,' ***'/) return end if end if ! ! Successful termination. ! ier = 0 return end subroutine edge ( in1, in2, x, y, lwk, iwk, list, lptr, lend, ier ) !*****************************************************************************80 ! !! EDGE swaps arcs to force two nodes to be adjacent. ! ! Discussion: ! ! Given a triangulation of N nodes and a pair of nodal ! indexes IN1 and IN2, this routine swaps arcs as necessary ! to force IN1 and IN2 to be adjacent. Only arcs which ! intersect IN1-IN2 are swapped out. If a Delaunay triangu- ! lation is input, the resulting triangulation is as close ! as possible to a Delaunay triangulation in the sense that ! all arcs other than IN1-IN2 are locally optimal. ! ! A sequence of calls to EDGE may be used to force the ! presence of a set of edges defining the boundary of a non- ! convex and/or multiply connected region (refer to Subrou- ! tine ADDCST), or to introduce barriers into the triangula- ! tion. Note that Subroutine GETNP will not necessarily ! return closest nodes if the triangulation has been con- ! strained by a call to EDGE. However, this is appropriate ! in some applications, such as triangle-based interpolation ! on a nonconvex domain. ! ! Modified: ! ! 16 June 2007 ! ! Author: ! ! Robert Renka, ! Department of Computer Science, ! University of North Texas, ! renka@cs.unt.edu ! ! Reference: ! ! Robert Renka, ! Algorithm 751: TRIPACK, ! A Constrained Two-Dimensional Delaunay Triangulation Package, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 1, 1996. ! ! Parameters: ! ! Input, integer IN1, IN2, indexes (of X and Y) in the range 1 to N ! defining a pair of nodes to be connected by an arc. ! ! Input, real ( kind = 8 ) X(N), Y(N), the coordinates of the nodes. ! ! Input/output, integer LWK. On input, the number of columns reserved ! for IWK. This must be at least NI, the number of arcs which intersect ! IN1-IN2. (NI is bounded by N-3.) On output, the number of arcs which ! intersect IN1-IN2 (but not more than the input value of LWK) unless ! IER = 1 or IER = 3. LWK = 0 if and only if IN1 and IN2 were adjacent ! (or LWK=0) on input. ! ! Output, integer IWK(2*LWK), the indexes of the endpoints of the new ! arcs other than IN1-IN2 unless IER > 0 or LWK = 0. New arcs to the ! left of IN2-IN1 are stored in the first K-1 columns (left portion of ! IWK), column K contains zeros, and new arcs to the right of IN2-IN1 ! occupy columns K+1,...,LWK. (K can be determined by searching IWK ! for the zeros.) ! ! Input, integer LIST(*), LPTR(*), LEND(N), the data structure defining ! the triangulation. Refer to subroutine TRMESH. On output, updated ! if necessary to reflect the presence of an arc connecting IN1 and ! IN2 unless IER /= 0. The data structure has been altered if IER = 4. ! ! Output, integer IER = Error indicator: ! 0, if no errors were encountered. ! 1, if IN1 < 1, IN2 .LT. 1, IN1 = IN2, or LWK < 0 on input. ! 2, if more space is required in IWK. ! 3, if IN1 and IN2 could not be connected due to either an invalid ! data structure or collinear nodes (and floating point error). ! 4, if an error flag was returned by OPTIM. ! ! Local parameters: ! ! DX,DY = Components of arc N1-N2. ! ! I = DO-loop index and column index for IWK ! IERR = Error flag returned by Subroutine OPTIM ! IWC = IWK index between IWF and IWL -- NL->NR is ! stored in IWK(1,IWC)->IWK(2,IWC) ! IWCP1 = IWC + 1 ! IWEND = Input or output value of LWK ! IWF = IWK (column) index of the first (leftmost) arc ! which intersects IN1->IN2 ! IWL = IWK (column) index of the last (rightmost) are ! which intersects IN1->IN2 ! LFT = Flag used to determine if a swap results in the ! new arc intersecting IN1-IN2 -- LFT = 0 iff ! N0 = IN1, LFT = -1 implies N0 LEFT IN1->IN2, ! and LFT = 1 implies N0 LEFT IN2->IN1 ! LP21 = Unused parameter returned by SWAP ! LP = List pointer (index) for LIST and LPTR ! LPL = Pointer to the last neighbor of IN1 or NL ! N0 = Neighbor of N1 or node opposite NR->NL ! N1,N2 = Local copies of IN1 and IN2 ! N1FRST = First neighbor of IN1 ! N1LST = (Signed) last neighbor of IN1 ! NEXT = Node opposite NL->NR ! NIT = Flag or number of iterations employed by OPTIM ! NL,NR = Endpoints of an arc which intersects IN1-IN2 ! with NL LEFT IN1->IN2 ! X0,Y0 = Coordinates of N0 ! X1,Y1 = Coordinates of IN1 ! X2,Y2 = Coordinates of IN2 ! implicit none real ( kind = 8 ) dx real ( kind = 8 ) dy integer i integer ier integer ierr integer in1 integer in2 integer iwc integer iwcp1 integer iwend integer iwf integer iwk(2,*) integer iwl logical left integer lend(*) integer lft integer list(*) integer lp integer lp21 integer lpl integer lptr(*) integer lwk integer n0 integer n1 integer n1frst integer n1lst integer n2 integer next integer nit integer nl integer nr real ( kind = 8 ) x(*) real ( kind = 8 ) x0 real ( kind = 8 ) x1 real ( kind = 8 ) x2 real ( kind = 8 ) y(*) real ( kind = 8 ) y0 real ( kind = 8 ) y1 real ( kind = 8 ) y2 ! ! Store IN1, IN2, and LWK in local variables and test for errors. ! n1 = in1 n2 = in2 iwend = lwk if ( n1 < 1 .or. n2 < 1 .or. n1 == n2 .or. iwend < 0 ) then ier = 1 return end if ! ! Test for N2 as a neighbor of N1. LPL points to the last neighbor of N1. ! lpl = lend(n1) n0 = abs ( list(lpl) ) lp = lpl do ! ! IN1 and IN2 were adjacent on input. ! if ( n0 == n2 ) then ier = 0 return end if lp = lptr(lp) n0 = list(lp) if ( lp == lpl ) then exit end if end do ! ! Initialize parameters. ! iwl = 0 nit = 0 ! ! Store the coordinates of N1 and N2. ! 2 continue x1 = x(n1) y1 = y(n1) x2 = x(n2) y2 = y(n2) ! ! Set NR and NL to adjacent neighbors of N1 such that ! NR LEFT N2->N1 and NL LEFT N1->N2, ! (NR Forward N1->N2 or NL Forward N1->N2), and ! (NR Forward N2->N1 or NL Forward N2->N1). ! ! Initialization: Set N1FRST and N1LST to the first and ! (signed) last neighbors of N1, respectively, and ! initialize NL to N1FRST. ! lpl = lend(n1) n1lst = list(lpl) lp = lptr(lpl) n1frst = list(lp) nl = n1frst if ( n1lst < 0 ) then go to 4 end if ! ! N1 is an interior node. Set NL to the first candidate ! for NR (NL LEFT N2->N1). ! 3 continue if ( left(x2,y2,x1,y1,x(nl),y(nl)) ) then go to 4 end if lp = lptr(lp) nl = list(lp) if ( nl /= n1frst ) then go to 3 end if ! ! All neighbors of N1 are strictly left of N1->N2. ! go to 5 ! ! NL = LIST(LP) LEFT N2->N1. Set NR to NL and NL to the ! following neighbor of N1. ! 4 continue nr = nl lp = lptr(lp) nl = abs ( list(lp) ) ! ! NL LEFT N1->N2 and NR LEFT N2->N1. The Forward tests ! are employed to avoid an error associated with ! collinear nodes. ! if ( left(x1,y1,x2,y2,x(nl),y(nl)) ) then dx = x2-x1 dy = y2-y1 if ((dx*(x(nl)-x1)+dy*(y(nl)-y1) >= 0. .or. & dx*(x(nr)-x1)+dy*(y(nr)-y1) >= 0.) .and. & (dx*(x(nl)-x2)+dy*(y(nl)-y2) <= 0. .or. & dx*(x(nr)-x2)+dy*(y(nr)-y2) <= 0.)) go to 6 ! ! NL-NR does not intersect N1-N2. However, there is ! another candidate for the first arc if NL lies on ! the line N1-N2. ! if ( .not. left(x2,y2,x1,y1,x(nl),y(nl)) ) go to 5 end if ! ! Bottom of loop. ! if ( nl /= n1frst ) then go to 4 end if ! ! Either the triangulation is invalid or N1-N2 lies on the ! convex hull boundary and an edge NR->NL (opposite N1 and ! intersecting N1-N2) was not found due to floating point ! error. Try interchanging N1 and N2 -- NIT > 0 iff this ! has already been done. ! 5 continue if ( 0 < nit ) then go to 33 end if nit = 1 n1 = n2 n2 = in1 go to 2 ! ! Store the ordered sequence of intersecting edges NL->NR in ! IWK(1,IWL)->IWK(2,IWL). ! 6 continue iwl = iwl + 1 if (iwl > iwend) then ier = 2 return end if iwk(1,iwl) = nl iwk(2,iwl) = nr ! ! Set NEXT to the neighbor of NL which follows NR. ! lpl = lend(nl) lp = lptr(lpl) ! ! Find NR as a neighbor of NL. The search begins with ! the first neighbor. ! 7 continue if (list(lp) == nr) go to 8 lp = lptr(lp) if (lp /= lpl) go to 7 ! ! NR must be the last neighbor, and NL->NR cannot be a ! boundary edge. ! if (list(lp) /= nr) then go to 33 end if ! ! Set NEXT to the neighbor following NR, and test for ! termination of the store loop. ! 8 continue lp = lptr(lp) next = abs ( list(lp) ) if (next == n2) then go to 9 end if ! ! Set NL or NR to NEXT. ! if ( left(x1,y1,x2,y2,x(next),y(next)) ) then nl = next else nr = next end if go to 6 ! ! IWL is the number of arcs which intersect N1-N2. ! Store LWK. ! 9 continue lwk = iwl iwend = iwl ! ! Initialize for edge swapping loop -- all possible swaps ! are applied (even if the new arc again intersects ! N1-N2), arcs to the left of N1->N2 are stored in the ! left portion of IWK, and arcs to the right are stored in ! the right portion. IWF and IWL index the first and last ! intersecting arcs. ! iwf = 1 ! ! Top of loop -- set N0 to N1 and NL->NR to the first edge. ! IWC points to the arc currently being processed. LFT ! <= 0 iff N0 LEFT N1->N2. ! 10 continue lft = 0 n0 = n1 x0 = x1 y0 = y1 nl = iwk(1,iwf) nr = iwk(2,iwf) iwc = iwf ! ! Set NEXT to the node opposite NL->NR unless IWC is the last arc. ! 11 continue if (iwc == iwl) go to 21 iwcp1 = iwc + 1 next = iwk(1,iwcp1) if (next /= nl) go to 16 next = iwk(2,iwcp1) ! ! NEXT RIGHT N1->N2 and IWC < IWL. Test for a possible swap. ! if ( .not. left(x0,y0,x(nr),y(nr),x(next),y(next)) ) then go to 14 end if if (lft >= 0) then go to 12 end if if ( .not. left(x(nl),y(nl),x0,y0,x(next),y(next)) ) then go to 14 end if ! ! Replace NL->NR with N0->NEXT. ! call swap (next,n0,nl,nr, list,lptr,lend, lp21) iwk(1,iwc) = n0 iwk(2,iwc) = next go to 15 ! ! Swap NL-NR for N0-NEXT, shift columns IWC+1,...,IWL to ! the left, and store N0-NEXT in the right portion of IWK. ! 12 continue call swap (next,n0,nl,nr, list,lptr,lend, lp21) do i = iwcp1,iwl iwk(1,i-1) = iwk(1,i) iwk(2,i-1) = iwk(2,i) end do iwk(1,iwl) = n0 iwk(2,iwl) = next iwl = iwl - 1 nr = next go to 11 ! ! A swap is not possible. Set N0 to NR. ! 14 continue n0 = nr x0 = x(n0) y0 = y(n0) lft = 1 ! ! Advance to the next arc. ! 15 continue nr = next iwc = iwc + 1 go to 11 ! ! NEXT LEFT N1->N2, NEXT /= N2, and IWC < IWL. ! Test for a possible swap. ! 16 continue if ( .not. left(x(nl),y(nl),x0,y0,x(next),y(next)) ) then go to 19 end if if (lft <= 0) then go to 17 end if if ( .not. left(x0,y0,x(nr),y(nr),x(next),y(next)) ) then go to 19 end if ! ! Replace NL->NR with NEXT->N0. ! call swap (next,n0,nl,nr, list,lptr,lend, lp21) iwk(1,iwc) = next iwk(2,iwc) = n0 go to 20 ! ! Swap NL-NR for N0-NEXT, shift columns IWF,...,IWC-1 to ! the right, and store N0-NEXT in the left portion of IWK. ! 17 continue call swap (next,n0,nl,nr, list,lptr,lend, lp21) do i = iwc-1,iwf,-1 iwk(1,i+1) = iwk(1,i) iwk(2,i+1) = iwk(2,i) end do iwk(1,iwf) = n0 iwk(2,iwf) = next iwf = iwf + 1 go to 20 ! ! A swap is not possible. Set N0 to NL. ! 19 continue n0 = nl x0 = x(n0) y0 = y(n0) lft = -1 ! ! Advance to the next arc. ! 20 continue nl = next iwc = iwc + 1 go to 11 ! ! N2 is opposite NL->NR (IWC = IWL). ! 21 continue if (n0 == n1) go to 24 if (lft < 0) go to 22 ! ! N0 RIGHT N1->N2. Test for a possible swap. ! if ( .not. left(x0,y0,x(nr),y(nr),x2,y2) ) go to 10 ! ! Swap NL-NR for N0-N2 and store N0-N2 in the right ! portion of IWK. ! call swap (n2,n0,nl,nr, list,lptr,lend, lp21) iwk(1,iwl) = n0 iwk(2,iwl) = n2 iwl = iwl - 1 go to 10 ! ! N0 LEFT N1->N2. Test for a possible swap. ! 22 continue if ( .not. left(x(nl),y(nl),x0,y0,x2,y2) ) then go to 10 end if ! ! Swap NL-NR for N0-N2, shift columns IWF,...,IWL-1 to the ! right, and store N0-N2 in the left portion of IWK. ! call swap ( n2, n0, nl, nr, list, lptr, lend, lp21 ) i = iwl 23 continue iwk(1,i) = iwk(1,i-1) iwk(2,i) = iwk(2,i-1) i = i - 1 if (i > iwf) then go to 23 end if iwk(1,iwf) = n0 iwk(2,iwf) = n2 iwf = iwf + 1 go to 10 ! ! IWF = IWC = IWL. Swap out the last arc for N1-N2 and ! store zeros in IWK. ! 24 continue call swap (n2,n1,nl,nr, list,lptr,lend, lp21) iwk(1,iwc) = 0 iwk(2,iwc) = 0 ! ! Optimization procedure. ! if ( 1 < iwc ) then ! ! Optimize the set of new arcs to the left of IN1->IN2. ! nit = 3*(iwc-1) call optim ( x, y, iwc-1, list, lptr, lend, nit, iwk, ierr ) if ( ierr /= 0 ) then go to 34 end if end if if ( iwc < iwend ) then ! ! Optimize the set of new arcs to the right of IN1->IN2. ! nit = 3*(iwend-iwc) call optim ( x, y, iwend-iwc, list, lptr, lend, nit, & iwk(1,iwc+1), ierr ) if (ierr /= 0) then go to 34 end if end if ! ! Successful termination. ! ier = 0 return ! ! Invalid triangulation data structure or collinear nodes ! on convex hull boundary. ! 33 ier = 3 write (*,130) in1, in2 130 format (//5x,'*** error in edge: invalid triangula', & 'tion or null triangles on boundary'/ & 9x,'in1 =',i4,', in2=',i4/) return ! ! Error flag returned by OPTIM. ! 34 ier = 4 write (*,140) nit, ierr 140 format (//5x,'*** error in optim: nit = ',i4, & ', ier = ',i1,' ***'/) return end subroutine getnp ( ncc, lcc, n, x, y, list, lptr, lend, l, npts, ds, ier ) !*****************************************************************************80 ! !! GETNP sets the next nearest node to a given node. ! ! Discussion: ! ! Given a triangulation of N nodes and an array NPTS con- ! taining the indexes of L-1 nodes ordered by distance from ! NPTS(1), this subroutine sets NPTS(L) to the index of the ! next node in the sequence -- the node, other than NPTS(1), ! ...,NPTS(L-1), which is closest to NPTS(1). Thus, the ! ordered sequence of K closest nodes to N1 (including N1) ! may be determined by K-1 calls to GETNP with NPTS(1) = N1 ! and L = 2,3,...,K for K >= 2. Note that NPTS must ! include constraint nodes as well as non-constraint nodes. ! Thus, a sequence of K1 closest non-constraint nodes to N1 ! must be obtained as a subset of the closest K2 nodes to N1 ! for some K2 >= K1. ! ! The terms closest and distance have special definitions ! when constraint nodes are present in the triangulation. ! Nodes N1 and N2 are said to be visible from each other if ! and only if the line segment N1-N2 intersects no constraint ! arc (except possibly itself) and is not an interi- ! or constraint arc (arc whose interior lies in a constraint ! region). A path from N1 to N2 is an ordered sequence of ! nodes, with N1 first and N2 last, such that adjacent path ! elements are visible from each other. The path length is ! the sum of the Euclidean distances between adjacent path ! nodes. Finally, the distance from N1 to N2 is defined to ! be the length of the shortest path from N1 to N2. ! ! The algorithm uses the property of a Delaunay triangulation ! that the K-th closest node to N1 is a neighbor of one ! of the K-1 closest nodes to N1. With the definition of ! distance used here, this property holds when constraints ! are present as long as non-constraint arcs are locally ! optimal. ! ! Modified: ! ! 16 June 2007 ! ! Author: ! ! Robert Renka, ! Department of Computer Science, ! University of North Texas, ! renka@cs.unt.edu ! ! Reference: ! ! Robert Renka, ! Algorithm 751: TRIPACK, ! A Constrained Two-Dimensional Delaunay Triangulation Package, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 1, 1996. ! ! Parameters: ! ! Input, integer NCC, the number of constraints. NCC >= 0. ! ! Input, integer LCC(*), a list of constraint curve starting indexes (or ! dummy array of length 1 if NCC = 0). Refer to subroutine ADDCST. ! ! Input, integer N, the number of nodes in the triangulation. N >= 3. ! ! Input, real ( kind = 8 ) X(N), Y(N), the coordinates of the nodes with ! non-constraint nodes in the first LCC(1)-1 locations if NCC > 0. ! ! Input, integer LIST(*), LPTR(*), LEND(N), the triangulation data ! structure. Refer to subroutine TRMESH. ! ! Input, integer L, the number of nodes in the sequence on output. ! 2 <= L <= N. ! ! Input/output, integer NPTS(L), on input, the indexes of the L-1 closest ! nodes to NPTS(1) in the first L-1 locations. On output, updated with ! the index of the L-th closest node to NPTS(1) in position L unless ! IER /= 0. ! ! Input/output, real ( kind = 8 ) DS(L), the distance (defined above) ! between NPTS(1) and NPTS(I) in the I-th position for I = 1,...,L-1. ! Thus, DS(1) = 0. On output, updated with the distance between NPTS(1) ! and NPTS(L) in position L unless IER /= 0. ! ! Output, integer IER = Error indicator: ! 0 if no errors were encountered. ! -1 if NCC, N, L, or an LCC entry is outside its valid range on input. ! K if NPTS(K) is not a valid index in the range 1 to N. ! implicit none integer l integer n real ( kind = 8 ) dc real ( kind = 8 ) dl real ( kind = 8 ) ds(l) integer i integer ier integer ifrst integer ilast logical intsec logical isw integer j integer k integer km1 integer lcc(*) integer lcc1 integer lend(n) logical lft1 logical lft2 logical lft12 integer list(*) integer lm1 integer lp integer lpcl integer lpk integer lpkl integer lptr(*) integer n1 integer nc integer ncc logical ncf integer nf1 integer nf2 integer nj logical njf integer nk integer nkbak integer nkfor integer nl integer nn integer npts(l) logical skip logical sksav logical vis real ( kind = 8 ) x(n) real ( kind = 8 ) x1 real ( kind = 8 ) xc real ( kind = 8 ) xj real ( kind = 8 ) xk real ( kind = 8 ) y(n) real ( kind = 8 ) y1 real ( kind = 8 ) yc real ( kind = 8 ) yj real ( kind = 8 ) yk ! ! Store parameters in local variables and test for errors. ! LCC1 indexes the first constraint node. ! ier = -1 nn = n lcc1 = nn + 1 lm1 = l - 1 if ( ncc < 0 ) then return end if if ( lm1 < 1 .or. lm1 >= nn) then return end if if (ncc == 0) then if (nn < 3) then return end if else do i = ncc,1,-1 if (lcc1 - lcc(i) < 3) return lcc1 = lcc(i) end do if (lcc1 < 1) then return end if end if ! ! Test for an invalid index in NPTS. ! do k = 1,lm1 nk = npts(k) if ( nk < 1 .or. nk > nn ) then ier = k return end if end do ! ! Store N1 = NPTS(1) and mark the elements of NPTS. ! n1 = npts(1) x1 = x(n1) y1 = y(n1) do k = 1,lm1 nk = npts(k) lend(nk) = -lend(nk) end do ! ! Candidates NC for NL = NPTS(L) are the unmarked visible ! neighbors of nodes NK in NPTS. ISW is an initialization ! switch set to .TRUE. when NL and its distance DL from N1 ! have been initialized with the first candidate encountered. ! isw = .false. dl = 0.0D+00 ! ! Loop on marked nodes NK = NPTS(K). LPKL indexes the last ! neighbor of NK in LIST. ! do k = 1,lm1 km1 = k - 1 nk = npts(k) xk = x(nk) yk = y(nk) lpkl = -lend(nk) nkfor = 0 nkbak = 0 vis = .true. ! ! NK is a constraint node. Set NKFOR and NKBAK to the ! constraint nodes which follow and precede NK. IFRST ! and ILAST are set to the first and last nodes in the ! constraint containing NK. ! if (nk >= lcc1) then ifrst = nn + 1 do i = ncc,1,-1 ilast = ifrst - 1 ifrst = lcc(i) if ( nk >= ifrst ) then exit end if end do if (nk < ilast) then nkfor = nk + 1 else nkfor = ifrst end if if (nk > ifrst) then nkbak = nk - 1 else nkbak = ilast end if ! ! Initialize VIS to TRUE iff NKFOR precedes NKBAK in the ! adjacency list for NK -- the first neighbor is visible and is not NKBAK. ! lpk = lpkl do lpk = lptr(lpk) nc = abs ( list(lpk) ) if ( nc == nkfor .or. nc == nkbak ) then exit end if end do vis = nc == nkfor end if ! ! Loop on neighbors NC of NK, bypassing marked and nonvisible neighbors. ! lpk = lpkl 7 continue lpk = lptr(lpk) nc = abs ( list(lpk) ) if ( nc == nkbak ) then vis = .true. end if ! ! VIS = .FALSE. iff NK-NC is an interior constraint arc ! (NK is a constraint node and NC lies strictly between ! NKFOR and NKBAK). ! if ( .not. vis ) go to 15 if ( nc == nkfor ) then vis = .false. end if if ( lend(nc) < 0 ) go to 15 ! ! Initialize distance DC between N1 and NC to Euclidean distance. ! xc = x(nc) yc = y(nc) dc = sqrt((xc-x1)*(xc-x1) + (yc-y1)*(yc-y1)) if (isw .and. dc >= dl) go to 15 if (k == 1) then go to 14 end if ! ! K >= 2. Store the pointer LPCL to the last neighbor of NC. ! lpcl = lend(nc) ! ! Set DC to the length of the shortest path from N1 to NC ! which has not previously been encountered and which is ! a viable candidate for the shortest path from N1 to NL. ! This is Euclidean distance iff NC is visible from N1. ! Since the shortest path from N1 to NL contains only ele- ! ments of NPTS which are constraint nodes (in addition to ! N1 and NL), only these need be considered for the path ! from N1 to NC. Thus, for distance function D(A,B) and ! J = 1,...,K, DC = min(D(N1,NJ) + D(NJ,NC)) over con- ! straint nodes NJ = NPTS(J) which are visible from NC. ! do j = 1,km1 nj = npts(j) if ( 1 < j .and. nj < lcc1 ) then go to 13 end if ! ! If NC is a visible neighbor of NJ, a path from N1 to NC ! containing NJ has already been considered. Thus, NJ may ! be bypassed if it is adjacent to NC. ! lp = lpcl 8 continue lp = lptr(lp) if ( nj == abs ( list(lp) ) ) then go to 12 end if if (lp /= lpcl) then go to 8 end if ! ! NJ is a constraint node (unless J=1) not adjacent to NC, ! and is visible from NC iff NJ-NC is not intersected by ! a constraint arc. Loop on constraints I in reverse ! order. ! xj = x(nj) yj = y(nj) ifrst = nn+1 do 11 i = ncc,1,-1 ilast = ifrst - 1 ifrst = lcc(i) nf1 = ilast ncf = nf1 == nc njf = nf1 == nj skip = ncf .or. njf ! ! Loop on boundary constraint arcs NF1-NF2 which contain ! neither NC nor NJ. NCF and NJF are TRUE iff NC (or NJ) ! has been encountered in the constraint, and SKIP = ! .TRUE. iff NF1 = NC or NF1 = NJ. ! do nf2 = ifrst,ilast if (nf2 == nc) ncf = .true. if (nf2 == nj) njf = .true. sksav = skip skip = nf2 == nc .or. nf2 == nj ! ! The last constraint arc in the constraint need not be ! tested if none of the arcs have been skipped. ! if ( sksav .or. skip .or. & (nf2 == ilast .and. & .not. ncf .and. .not. njf) ) then go to 9 end if if ( intsec(x(nf1),y(nf1),x(nf2),y(nf2), & xc,yc,xj,yj) ) then go to 12 end if 9 continue nf1 = nf2 end do ! ! NC and NJ are constraint nodes in the same constraint. ! NC-NJ is intersected by an interior constraint arc iff ! 1) NC LEFT NF2->NF1 and (NJ LEFT NF1->NC and NJ LEFT NC->NF2) or ! 2) NC .NOT. LEFT NF2->NF1 and (NJ LEFT NF1->NC or NJ LEFT NC->NF2), ! where NF1, NC, NF2 are consecutive constraint nodes. ! if (.not. ncf .or. .not. njf) go to 11 if (nc /= ifrst) then nf1 = nc - 1 else nf1 = ilast end if if (nc /= ilast) then nf2 = nc + 1 else nf2 = ifrst end if lft1 = (xc-x(nf1))*(yj-y(nf1)) >= (xj-x(nf1))*(yc-y(nf1)) lft2 = (x(nf2)-xc)*(yj-yc) >= (xj-xc)*(y(nf2)-yc) lft12 = (x(nf1)-x(nf2))*(yc-y(nf2)) >= (xc-x(nf2))*(y(nf1)-y(nf2)) if ( (lft1 .and. lft2) .or. (.not. lft12 & .and. (lft1 .or. lft2)) ) go to 12 11 continue ! ! NJ is visible from NC. Exit the loop with DC = Euclidean ! distance if J = 1. ! if (j == 1) then go to 14 end if dc = min(dc,ds(j) + sqrt((xc-xj)*(xc-xj) + & (yc-yj)*(yc-yj))) go to 13 ! ! NJ is not visible from NC or is adjacent to NC. Initialize DC ! with D(N1,NK) + D(NK,NC) if J = 1. ! 12 continue if (j == 1) then dc = ds(k) + sqrt((xc-xk)*(xc-xk) & + (yc-yk)*(yc-yk)) end if 13 continue end do ! ! Compare DC with DL. ! if ( isw .and. dc >= dl) then go to 15 end if ! ! The first (or a closer) candidate for NL has been encountered. ! 14 continue nl = nc dl = dc isw = .true. 15 continue if (lpk /= lpkl) then go to 7 end if end do ! ! Unmark the elements of NPTS and store NL and DL. ! do k = 1,lm1 nk = npts(k) lend(nk) = -lend(nk) end do npts(l) = nl ds(l) = dl ier = 0 return end function indxcc ( ncc, lcc, n, list, lend ) !*****************************************************************************80 ! !! INDXCC returns the index of an exterior constraint curve. ! ! Discussion: ! ! Given a constrained Delaunay triangulation, this ! function returns the index, if any, of an exterior constraint ! curve (an unbounded constraint region). An exterior ! constraint curve is assumed to be present if and only if the ! clockwise-ordered sequence of boundary nodes is a ! subsequence of a constraint node sequence. The triangulation ! adjacencies corresponding to constraint edges may or may ! not have been forced by a call to ADDCST, and the ! constraint region may or may not be valid (contain no nodes). ! ! Modified: ! ! 16 June 2007 ! ! Author: ! ! Robert Renka, ! Department of Computer Science, ! University of North Texas, ! renka@cs.unt.edu ! ! Reference: ! ! Robert Renka, ! Algorithm 751: TRIPACK, ! A Constrained Two-Dimensional Delaunay Triangulation Package, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 1, 1996. ! ! Parameters: ! ! Input, integer NCC, the number of constraints. NCC >= 0. ! ! Input, integer LCC(*), list of constraint curve starting indexes (or ! dummy array of length 1 if NCC = 0). Refer to subroutine ADDCST. ! ! Input, integer N, the number of nodes in the triangulation. N >= 3. ! ! Input, integer LIST(*), LEND(N), the data structure defining the ! triangulation. Refer to subroutine TRMESH. ! ! Output, integer INDXCC, index of the exterior constraint curve, if ! present, or 0 otherwise. ! implicit none integer n integer i integer ifrst integer ilast integer indxcc integer lcc(*) integer lend(n) integer list(*) integer lp integer n0 integer ncc integer nst integer nxt indxcc = 0 if ( ncc < 1 ) then return end if ! ! Set N0 to the boundary node with smallest index. ! n0 = 0 do n0 = n0 + 1 lp = lend(n0) if ( list(lp) <= 0 ) then exit end if end do ! ! Search in reverse order for the constraint I, if any, that ! contains N0. IFRST and ILAST index the first and last ! nodes in constraint I. ! i = ncc ilast = n do ifrst = lcc(i) if ( ifrst <= n0 ) then exit end if if ( i == 1 ) then return end if i = i - 1 ilast = ifrst - 1 end do ! ! N0 is in constraint I which indexes an exterior constraint ! curve iff the clockwise-ordered sequence of boundary ! node indexes beginning with N0 is increasing and bounded ! above by ILAST. ! nst = n0 do nxt = -list(lp) if ( nxt == nst ) then exit end if if ( nxt <= n0 .or. ilast < nxt ) then return end if n0 = nxt lp = lend(n0) end do ! ! Constraint I contains the boundary node sequence as a subset. ! indxcc = i return end subroutine insert ( k, lp, list, lptr, lnew ) !*****************************************************************************80 ! !! INSERT inserts K as a neighbor of N1. ! ! Discussion: ! ! This subroutine inserts K as a neighbor of N1 following ! N2, where LP is the LIST pointer of N2 as a neighbor of ! N1. Note that, if N2 is the last neighbor of N1, K will ! become the first neighbor (even if N1 is a boundary node). ! ! Modified: ! ! 16 June 2007 ! ! Author: ! ! Robert Renka, ! Department of Computer Science, ! University of North Texas, ! renka@cs.unt.edu ! ! Reference: ! ! Robert Renka, ! Algorithm 751: TRIPACK, ! A Constrained Two-Dimensional Delaunay Triangulation Package, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 1, 1996. ! ! Parameters: ! ! Input, integer K, the index of the node to be inserted. ! ! Input, integer LP, the LIST pointer of N2 as a neighbor of N1. ! ! Input/output, integer LIST(*), LPTR(*), LNEW, the data structure ! defining the triangulation. Refer to subroutine TRMESH. On output, ! the data structure has been updated to include node K. ! implicit none integer k integer list(*) integer lnew integer lp integer lptr(*) integer lsav lsav = lptr(lp) lptr(lp) = lnew list(lnew) = k lptr(lnew) = lsav lnew = lnew + 1 return end subroutine intadd ( kk, i1, i2, i3, list, lptr, lend, lnew ) !*****************************************************************************80 ! !! INTADD adds an interior point to a triangulation. ! ! Discussion: ! ! This subroutine adds an interior node to a triangulation ! of a set of points in the plane. The data structure is ! updated with the insertion of node KK into the triangle ! whose vertices are I1, I2, and I3. No optimization of the ! triangulation is performed. ! ! Modified: ! ! 16 June 2007 ! ! Author: ! ! Robert Renka, ! Department of Computer Science, ! University of North Texas, ! renka@cs.unt.edu ! ! Reference: ! ! Robert Renka, ! Algorithm 751: TRIPACK, ! A Constrained Two-Dimensional Delaunay Triangulation Package, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 1, 1996. ! ! Parameters: ! ! Input, integer KK, the index of the node to be inserted. 1 <= KK ! and KK must not be equal to I1, I2, or I3. ! ! Input, integer I1, I2, I3, indexes of the counterclockwise-ordered ! sequence of vertices of a triangle which contains node KK. ! ! Input/output, integer LIST(*), LPTR(*), LEND(N), LNEW, the data ! structure defining the triangulation. Refer to subroutine TRMESH. ! Triangle (I1,I2,I3) must be included in the triangulation. ! On output, updated with the addition of node KK. KK ! will be connected to nodes I1, I2, and I3. ! implicit none integer i1 integer i2 integer i3 integer k integer kk integer lend(*) integer list(*) integer lnew integer lp integer lptr(*) integer lstptr integer n1 integer n2 integer n3 k = kk ! ! Initialization. ! n1 = i1 n2 = i2 n3 = i3 ! ! Add K as a neighbor of I1, I2, and I3. ! lp = lstptr(lend(n1),n2,list,lptr) call insert (k,lp,list,lptr,lnew) lp = lstptr(lend(n2),n3,list,lptr) call insert (k,lp,list,lptr,lnew) lp = lstptr(lend(n3),n1,list,lptr) call insert (k,lp,list,lptr,lnew) ! ! Add I1, I2, and I3 as neighbors of K. ! list(lnew) = n1 list(lnew+1) = n2 list(lnew+2) = n3 lptr(lnew) = lnew + 1 lptr(lnew+1) = lnew + 2 lptr(lnew+2) = lnew lend(k) = lnew + 2 lnew = lnew + 3 return end function intsec ( x1, y1, x2, y2, x3, y3, x4, y4 ) !*****************************************************************************80 ! !! INTSEC determines if two line segments intersect. ! ! Discussion: ! ! Given a pair of line segments P1-P2 and P3-P4, this ! function returns the value .TRUE. if and only if P1-P2 ! shares one or more points with P3-P4. The line segments ! include their endpoints, and the four points need not be ! distinct. Thus, either line segment may consist of a ! single point, and the segments may meet in a V (which is ! treated as an intersection). Note that an incorrect ! decision may result from floating point error if the four ! endpoints are nearly collinear. ! ! Modified: ! ! 19 May 2005 ! ! Author: ! ! Robert Renka, ! Department of Computer Science, ! University of North Texas, ! renka@cs.unt.edu ! ! Reference: ! ! Robert Renka, ! Algorithm 751: TRIPACK, ! A Constrained Two-Dimensional Delaunay Triangulation Package, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 1, 1996. ! ! Parameters: ! ! Input, real ( kind = 8 ) X1, Y1 = Coordinates of P1. ! ! Input, real ( kind = 8 ) X2, Y2 = Coordinates of P2. ! ! Input, real ( kind = 8 ) X3, Y3 = Coordinates of P3. ! ! Input, real ( kind = 8 ) X4, Y4 = Coordinates of P4. ! ! Output, logical INTSEC, the logical value defined above. ! implicit none real ( kind = 8 ) a real ( kind = 8 ) b real ( kind = 8 ) d real ( kind = 8 ) dx12 real ( kind = 8 ) dx31 real ( kind = 8 ) dx34 real ( kind = 8 ) dy12 real ( kind = 8 ) dy31 real ( kind = 8 ) dy34 logical intsec real ( kind = 8 ) x1 real ( kind = 8 ) x2 real ( kind = 8 ) x3 real ( kind = 8 ) x4 real ( kind = 8 ) y1 real ( kind = 8 ) y2 real ( kind = 8 ) y3 real ( kind = 8 ) y4 ! ! Test for overlap between the smallest rectangles that ! contain the line segments and have sides parallel to ! the axes. ! if (( x1 < x3 .and. x1 < x4 .and. x2 < x3 .and. x2 < x4) .or. & ( x3 < x1 .and. x1 > x4 .and. x2 > x3 .and. x2 > x4) .or. & ( y1 < y3 .and. y1 < y4 .and. y2 < y3 .and. y2 < y4) .or. & ( y3 < y1 .and. y4 < y1 .and. y2 > y3 .and. y2 > y4)) then intsec = .false. return end if ! ! Compute A = P4-P3 X P1-P3, B = P2-P1 X P1-P3, and ! D = P2-P1 X P4-P3 (Z components). ! dx12 = x2 - x1 dy12 = y2 - y1 dx34 = x4 - x3 dy34 = y4 - y3 dx31 = x1 - x3 dy31 = y1 - y3 a = dx34 * dy31 - dx31 * dy34 b = dx12 * dy31 - dx31 * dy12 d = dx12 * dy34 - dx34 * dy12 ! ! D /= 0 and the point of intersection of the lines defined by the line ! segments is P = P1 + (A/D)*(P2-P1) = P3 + (B/D)*(P4-P3). ! if ( d /= 0.0D+00 ) then intsec = & 0.0D+00 <= a / d .and. a / d <= 1.0D+00 .and. & 0.0D+00 <= b / d .and. b / d <= 1.0D+00 ! ! D == 0 and thus either the line segments are parallel, ! or one (or both) of them is a single point. ! else intsec = ( a == 0.0D+00 .and. b == 0.0D+00 ) end if return end function jrand ( n, ix, iy, iz ) !*****************************************************************************80 ! !! JRAND returns a uniformly distributed random integer between 1 and N. ! ! Modified: ! ! 16 June 2007 ! ! Author: ! ! Robert Renka, ! Department of Computer Science, ! University of North Texas, ! renka@cs.unt.edu ! ! Reference: ! ! Brian Wichmann, David Hill, ! An Efficient and Portable Pseudo-random Number Generator, ! Applied Statistics, ! Volume 31, Number 2, 1982, pages 188-190. ! ! Parameters: ! ! Input, integer N, the maximum value to be returned. ! ! Input/output, integer IX, IY, IZ, seeds initialized to values in ! the range 1 to 30,000 before the first call to JRAND, and not altered ! by the user between subsequent calls (unless a sequence of random ! numbers is to be repeated by reinitializing the seeds). ! ! Output, integer JRAND, random integer in the range 1 to N. ! ! Local parameters: ! ! U = Pseudo-random number uniformly distributed in the interval (0,1). ! X = Pseudo-random number in the range 0 to 3 whose fractional part is U. ! implicit none integer ix integer iy integer iz integer jrand integer n real ( kind = 8 ) u real ( kind = 8 ) x ix = mod ( 171 * ix, 30269 ) iy = mod ( 172 * iy, 30307 ) iz = mod ( 170 * iz, 30323 ) x = ( real ( ix, kind = 8 ) / 30269.0D+00 ) & + ( real ( iy, kind = 8 ) / 30307.0D+00 ) & + ( real ( iz, kind = 8 ) / 30323.0D+00 ) u = x - int ( x ) jrand = real ( n, kind = 8 ) * u + 1.0D+00 return end function left ( x1, y1, x2, y2, x0, y0 ) !*****************************************************************************80 ! !! LEFT determines whether a node is to the left of a line. ! ! Discussion: ! ! This function determines whether node N0 is to the left ! or to the right of the line through N1-N2 as viewed by an ! observer at N1 facing N2. ! ! Modified: ! ! 16 June 2007 ! ! Author: ! ! Robert Renka, ! Department of Computer Science, ! University of North Texas, ! renka@cs.unt.edu ! ! Reference: ! ! Robert Renka, ! Algorithm 751: TRIPACK, ! A Constrained Two-Dimensional Delaunay Triangulation Package, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 1, 1996. ! ! Parameters: ! ! Input, real ( kind = 8 ) X1, Y1, coordinates of N1. ! ! Input, real ( kind = 8 ) X2, Y2, coordinates of N2. ! ! Input, real ( kind = 8 ) X0, Y0, coordinates of N0. ! ! Output, logical LEFT, is .TRUE. if and only if (X0,Y0) is on or ! to the left of the directed line N1->N2. ! ! Local parameters: ! ! DX1,DY1 = X,Y components of the vector N1->N2 ! DX2,DY2 = X,Y components of the vector N1->N0 ! implicit none real ( kind = 8 ) dx1 real ( kind = 8 ) dx2 real ( kind = 8 ) dy1 real ( kind = 8 ) dy2 logical left real ( kind = 8 ) x0 real ( kind = 8 ) x1 real ( kind = 8 ) x2 real ( kind = 8 ) y0 real ( kind = 8 ) y1 real ( kind = 8 ) y2 dx1 = x2 - x1 dy1 = y2 - y1 dx2 = x0 - x1 dy2 = y0 - y1 ! ! If the sign of the vector cross product of N1->N2 and ! N1->N0 is positive, then sin(A) > 0, where A is the ! angle between the vectors, and thus A is in the range ! (0,180) degrees. ! left = dx1 * dy2 >= dx2 * dy1 return end function lstptr ( lpl, nb, list, lptr ) !*****************************************************************************80 ! !! LSTPTR returns the index of NB in the adjacency list for N0. ! ! Discussion: ! ! This function returns the index (LIST pointer) of NB in ! the adjacency list for N0, where LPL = LEND(N0). ! ! Modified: ! ! 16 June 2007 ! ! Author: ! ! Robert Renka, ! Department of Computer Science, ! University of North Texas, ! renka@cs.unt.edu ! ! Reference: ! ! Robert Renka, ! Algorithm 751: TRIPACK, ! A Constrained Two-Dimensional Delaunay Triangulation Package, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 1, 1996. ! ! Parameters: ! ! Input, integer LPL = LEND(N0). ! ! Input, integer NB, the index of the node whose pointer is to be ! returned. NB must be connected to N0. ! ! Input, integer LIST(*), LPTR(*), the data structure defining the ! triangulation. Refer to subroutine TRMESH. ! ! Output, integer LSTPTR, pointer such that LIST(LSTPTR) = NB or ! LIST(LSTPTR) = -NB, unless NB is not a neighbor of N0, in which ! case LSTPTR = LPL. ! implicit none integer list(*) integer lp integer lpl integer lptr(*) integer lstptr integer nb integer nd lp = lptr(lpl) do nd = list(lp) if ( nd == nb ) then exit end if lp = lptr(lp) if ( lp == lpl ) then exit end if end do lstptr = lp return end function nbcnt ( lpl, lptr ) !*****************************************************************************80 ! !! NBCNT returns the number of neighbors of a node. ! ! Discussion: ! ! This function returns the number of neighbors of a node ! in a triangulation created by TRMESH or TRMSHR. ! ! Modified: ! ! 25 November 2002 ! ! Author: ! ! Robert Renka, ! Department of Computer Science, ! University of North Texas, ! renka@cs.unt.edu ! ! Reference: ! ! Robert Renka, ! Algorithm 751: TRIPACK, ! A Constrained Two-Dimensional Delaunay Triangulation Package, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 1, 1996. ! ! Parameters: ! ! Input, integer LPL, the LIST pointer to the last neighbor of N0. ! LPL = LEND(N0). ! ! Input, integer LPTR(*), pointers associated with LIST. ! ! Output, integer NBCNT, the number of neighbors of N0. ! implicit none integer k integer lp integer lpl integer lptr(*) integer nbcnt lp = lpl k = 1 do lp = lptr(lp) if ( lp == lpl ) then exit end if k = k + 1 end do nbcnt = k return end function nearnd ( xp, yp, ist, n, x, y, list, lptr, lend, dsq ) !*****************************************************************************80 ! !! NEARND finds the nearest triangulation node to a point. ! ! Discussion: ! ! Given a point P in the plane and a Delaunay triangulation created by ! subroutine TRMESH or TRMSHR, this function returns the index of the ! nearest triangulation node to P. ! ! The algorithm consists of implicitly adding P to the ! triangulation, finding the nearest neighbor to P, and ! implicitly deleting P from the triangulation. Thus, it ! is based on the fact that, if P is a node in a Delaunay ! triangulation, the nearest node to P is a neighbor of P. ! ! Note that the number of candidates for NEARND ! (neighbors of P) is limited to LMAX defined in ! the PARAMETER statement below. ! ! Modified: ! ! 16 June 2007 ! ! Author: ! ! Robert Renka, ! Department of Computer Science, ! University of North Texas, ! renka@cs.unt.edu ! ! Reference: ! ! Robert Renka, ! Algorithm 751: TRIPACK, ! A Constrained Two-Dimensional Delaunay Triangulation Package, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 1, 1996. ! ! Parameters: ! ! Input, real ( kind = 8 ) XP, YP, the coordinates of the point P. ! ! Input, integer IST, the index of a node at which TRFIND begins ! the search. Search time depends on the proximity ! of this node to P. ! ! Input, integer N, the number of nodes in the triangulation. ! ! Input, real ( kind = 8 ) X(N), Y(N), the coordinates of the nodes. ! ! Input, integer LIST(*), LPTR(*), LEND(N), a data structure ! defining the triangulation. Refer to TRMESH. ! ! Output, real ( kind = 8 ) DSQ, the square of the distance between P and ! node NEARND. ! ! Output, integer NEARND, the index of the nearest node to P, ! or 0 if N < 3 or the triangulation data structure is invalid. ! implicit none integer, parameter :: lmax = 25 integer n real ( kind = 8 ) cos1 real ( kind = 8 ) cos2 real ( kind = 8 ) ds1 real ( kind = 8 ) dsq real ( kind = 8 ) dsr real ( kind = 8 ) dx11 real ( kind = 8 ) dx12 real ( kind = 8 ) dx21 real ( kind = 8 ) dx22 real ( kind = 8 ) dy11 real ( kind = 8 ) dy12 real ( kind = 8 ) dy21 real ( kind = 8 ) dy22 integer i1 integer i2 integer i3 integer ist integer l integer lend(n) integer list(*) integer listp(lmax) integer lp integer lp1 integer lp2 integer lpl integer lptr(*) integer lptrp(lmax) integer lstptr integer n1 integer n2 integer n3 integer nearnd integer nr integer nst real ( kind = 8 ) sin1 real ( kind = 8 ) sin2 real ( kind = 8 ) x(n) real ( kind = 8 ) xp real ( kind = 8 ) y(n) real ( kind = 8 ) yp ! ! Store local parameters and test for N invalid. ! nearnd = 0 dsq = -1.0D+00 if ( n < 3 ) then return end if nst = ist if ( nst < 1 .or. n < nst ) then nst = 1 end if ! ! Find a triangle (I1,I2,I3) containing P, or the rightmost ! (I1) and leftmost (I2) visible boundary nodes as viewed from P. ! call trfind ( nst, xp, yp, n, x, y, list, lptr, lend, i1, i2, i3 ) ! ! Test for collinear nodes. ! if ( i1 == 0 ) then return end if ! ! Store the linked list of 'neighbors' of P in LISTP and ! LPTRP. I1 is the first neighbor, and 0 is stored as ! the last neighbor if P is not contained in a triangle. ! L is the length of LISTP and LPTRP, and is limited to LMAX. ! if ( i3 /= 0 ) then listp(1) = i1 lptrp(1) = 2 listp(2) = i2 lptrp(2) = 3 listp(3) = i3 lptrp(3) = 1 l = 3 else n1 = i1 l = 1 lp1 = 2 listp(l) = n1 lptrp(l) = lp1 ! ! Loop on the ordered sequence of visible boundary nodes ! N1 from I1 to I2. ! do lpl = lend(n1) n1 = -list(lpl) l = lp1 lp1 = l+1 listp(l) = n1 lptrp(l) = lp1 if ( n1 == i2 .or. lmax <= lp1 ) then exit end if end do l = lp1 listp(l) = 0 lptrp(l) = 1 end if ! ! Initialize variables for a loop on arcs N1-N2 opposite P ! in which new 'neighbors' are 'swapped' in. N1 follows ! N2 as a neighbor of P, and LP1 and LP2 are the LISTP ! indexes of N1 and N2. ! lp2 = 1 n2 = i1 lp1 = lptrp(1) n1 = listp(lp1) ! ! Begin loop: find the node N3 opposite N1->N2. ! do lp = lstptr ( lend(n1), n2, list, lptr ) if ( list(lp) < 0 ) then go to 4 end if lp = lptr(lp) n3 = abs ( list(lp) ) ! ! Swap test: Exit the loop if L = LMAX. ! if ( lmax <= l ) then exit end if dx11 = x(n1) - x(n3) dx12 = x(n2) - x(n3) dx22 = x(n2) - xp dx21 = x(n1) - xp dy11 = y(n1) - y(n3) dy12 = y(n2) - y(n3) dy22 = y(n2) - yp dy21 = y(n1) - yp cos1 = dx11 * dx12 + dy11 * dy12 cos2 = dx22 * dx21 + dy22 * dy21 if ( 0.0D+00 <= cos1 .and. cos2 >= 0.0D+00 ) then go to 4 end if if ( cos1 < 0.0D+00 .and. cos2 < 0.0D+00 ) then go to 3 end if sin1 = dx11 * dy12 - dx12 * dy11 sin2 = dx22 * dy21 - dx21 * dy22 if ( sin1 * cos2 + cos1 * sin2 >= 0.0D+00 ) then go to 4 end if ! ! Swap: Insert N3 following N2 in the adjacency list for P. ! The two new arcs opposite P must be tested. ! 3 continue l = l+1 lptrp(lp2) = l listp(l) = n3 lptrp(l) = lp1 lp1 = l n1 = n3 cycle ! ! No swap: Advance to the next arc and test for termination ! on N1 = I1 (LP1 = 1) or N1 followed by 0. ! 4 continue if ( lp1 == 1 ) then exit end if lp2 = lp1 n2 = n1 lp1 = lptrp(lp1) n1 = listp(lp1) if ( n1 == 0 ) then exit end if end do ! ! Set NR and DSR to the index of the nearest node to P and ! its squared distance from P, respectively. ! nr = i1 dsr = ( x(nr) - xp )**2 + ( y(nr) - yp )**2 do lp = 2, l n1 = listp(lp) if ( n1 == 0 ) then cycle end if ds1 = ( x(n1) - xp )**2 + ( y(n1) - yp )**2 if ( ds1 < dsr ) then nr = n1 dsr = ds1 end if end do dsq = dsr nearnd = nr return end subroutine optim ( x, y, na, list, lptr, lend, nit, iwk, ier ) !*****************************************************************************80 ! !! OPTIM optimizes the quadrilateral portion of a triangulation. ! ! Discussion: ! ! Given a set of NA triangulation arcs, this subroutine ! optimizes the portion of the triangulation consisting of ! the quadrilaterals (pairs of adjacent triangles) which ! have the arcs as diagonals by applying the circumcircle ! test and appropriate swaps to the arcs. ! ! An iteration consists of applying the swap test and ! swaps to all NA arcs in the order in which they are ! stored. The iteration is repeated until no swap occurs ! or NIT iterations have been performed. The bound on the ! number of iterations may be necessary to prevent an ! infinite loop caused by cycling (reversing the effect of a ! previous swap) due to floating point inaccuracy when four ! or more nodes are nearly cocircular. ! ! Modified: ! ! 20 October 2005 ! ! Author: ! ! Robert Renka, ! Department of Computer Science, ! University of North Texas, ! renka@cs.unt.edu ! ! Reference: ! ! Robert Renka, ! Algorithm 751: TRIPACK, ! A Constrained Two-Dimensional Delaunay Triangulation Package, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 1, 1996. ! ! Parameters: ! ! Input, real ( kind = 8 ) X(*), Y(*), the nodal coordinates. ! ! Input, integer NA, the number of arcs in the set. 0 <= NA. ! ! Input/output, integer LIST(*), LPTR(*), LEND(N), data structure ! defining the triangulation. Refer to subroutine TRMESH. ! On output, updated to reflect the swaps. ! ! Input/output, integer NIT. On input, the maximum number of iterations ! to be performed. A reasonable value is 3*NA. 1 <= NIT. On output, ! the number of iterations performed. ! ! Input/output, integer IWK(2,NA), containing the nodal indexes of the ! arc endpoints (pairs of endpoints are stored in columns). On output, ! the information has been updated to reflect the swaps. ! ! Output, integer IER = Error indicator: ! 0, if no errors were encountered. ! 1, if a swap occurred on the last of MAXIT iterations, where MAXIT is ! the value of NIT on input. The new set of arcs in not necessarily ! optimal in this case. ! 2, if NA < 0 or NIT < 1 on input. ! 3, if IWK(2,I) is not a neighbor of IWK(1,I) for some I in the range 1 ! to NA. A swap may have occurred in this case. ! 4, if a zero pointer was returned by subroutine SWAP. ! ! Local parameters: ! ! I = Column index for IWK ! IO1,IO2 = Nodal indexes of the endpoints of an arc in IWK ! ITER = Iteration count ! LP = LIST pointer ! LP21 = Parameter returned by SWAP (not used) ! LPL = Pointer to the last neighbor of IO1 ! LPP = Pointer to the node preceding IO2 as a neighbor of IO1 ! MAXIT = Input value of NIT ! N1,N2 = Nodes opposite IO1->IO2 and IO2->IO1, respectively ! NNA = Local copy of NA ! SWP = Flag set to TRUE iff a swap occurs in the optimization loop ! implicit none integer na integer i integer ier integer io1 integer io2 integer iter integer iwk(2,na) integer lend(*) integer list(*) integer lp integer lp21 integer lpl integer lpp integer lptr(*) integer maxit integer n1 integer n2 integer nit integer nna logical swp logical swptst real ( kind = 8 ) x(*) real ( kind = 8 ) y(*) nna = na maxit = nit if ( nna < 0 .or. maxit < 1 ) then nit = 0 ier = 2 return end if ! ! Initialize iteration count ITER and test for NA = 0. ! iter = 0 if ( nna == 0 ) then nit = 0 ier = 0 return end if ! ! Top of loop. ! SWP = TRUE iff a swap occurred in the current iteration. ! do if ( iter == maxit ) then nit = maxit ier = 1 return end if iter = iter + 1 swp = .false. ! ! Inner loop on arcs IO1-IO2. ! do i = 1, nna io1 = iwk(1,i) io2 = iwk(2,i) ! ! Set N1 and N2 to the nodes opposite IO1->IO2 and ! IO2->IO1, respectively. Determine the following: ! ! LPL = pointer to the last neighbor of IO1, ! LP = pointer to IO2 as a neighbor of IO1, and ! LPP = pointer to the node N2 preceding IO2. ! lpl = lend(io1) lpp = lpl lp = lptr(lpp) do if ( list(lp) == io2 ) then go to 3 end if lpp = lp lp = lptr(lpp) if ( lp == lpl ) then exit end if end do ! ! IO2 should be the last neighbor of IO1. Test for no ! arc and bypass the swap test if IO1 is a boundary node. ! if ( abs ( list(lp) ) /= io2 ) then nit = iter ier = 3 return end if if ( list(lp) < 0 ) then go to 4 end if ! ! Store N1 and N2, or bypass the swap test if IO1 is a ! boundary node and IO2 is its first neighbor. ! 3 continue n2 = list(lpp) if ( n2 < 0 ) then go to 4 end if lp = lptr(lp) n1 = abs ( list(lp) ) ! ! Test IO1-IO2 for a swap, and update IWK if necessary. ! if ( .not. swptst ( n1, n2, io1, io2, x, y ) ) then go to 4 end if call swap ( n1, n2, io1, io2, list, lptr, lend, lp21 ) if ( lp21 == 0 ) then nit = iter ier = 4 return end if swp = .true. iwk(1,i) = n1 iwk(2,i) = n2 4 continue end do if ( .not. swp ) then exit end if end do ! ! Successful termination. ! 5 continue nit = iter ier = 0 return end function store ( x ) !*****************************************************************************80 ! !! STORE forces its argument to be stored. ! ! Discussion: ! ! This function forces its argument X to be stored in a ! memory location, thus providing a means of determining ! floating point number characteristics (such as the machine ! precision) when it is necessary to avoid computation in ! high precision registers. ! ! Modified: ! ! 16 June 2007 ! ! Author: ! ! Robert Renka, ! Department of Computer Science, ! University of North Texas, ! renka@cs.unt.edu ! ! Reference: ! ! Robert Renka, ! Algorithm 751: TRIPACK, ! A Constrained Two-Dimensional Delaunay Triangulation Package, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 1, 1996. ! ! Parameters: ! ! Input, real ( kind = 8 ) X, the value to be stored. ! ! Output, real ( kind = 8 ) STORE, the value of X after it has been stored ! and possibly truncated or rounded to the single precision word length. ! implicit none real ( kind = 8 ) store real ( kind = 8 ) x real ( kind = 8 ) y common /stcom/ y y = x store = y return end subroutine swap ( in1, in2, io1, io2, list, lptr, lend, lp21 ) !*****************************************************************************80 ! !! SWAP adjusts a triangulation by swapping a diagonal arc. ! ! Discussion: ! ! Given a triangulation of a set of points on the unit ! sphere, this subroutine replaces a diagonal arc in a ! strictly convex quadrilateral (defined by a pair of adja- ! cent triangles) with the other diagonal. Equivalently, a ! pair of adjacent triangles is replaced by another pair ! having the same union. ! ! Modified: ! ! 16 June 2007 ! ! Author: ! ! Robert Renka, ! Department of Computer Science, ! University of North Texas, ! renka@cs.unt.edu ! ! Reference: ! ! Robert Renka, ! Algorithm 751: TRIPACK, ! A Constrained Two-Dimensional Delaunay Triangulation Package, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 1, 1996. ! ! Parameters: ! ! Input, integer IN1, IN2, IO1, IO2, the nodal indexes of the vertices of ! the quadrilateral. IO1-IO2 is replaced by IN1-IN2. (IO1,IO2,IN1) ! and (IO2,IO1,IN2) must be triangles on input. ! ! Input/output, integer LIST(*), LPTR(*), LEND(N), the data structure ! defining the triangulation. Refer to subroutine TRMESH. On output, ! updated with the swap; triangles (IO1,IO2,IN1) and (IO2,IO1,IN2) are ! replaced by (IN1,IN2,IO2) and (IN2,IN1,IO1) unless LP21 = 0. ! ! Output, integer LP21, the index of IN1 as a neighbor of IN2 after the ! swap is performed unless IN1 and IN2 are adjacent on input, in which ! case LP21 = 0. ! ! Local parameters: ! ! LP, LPH, LPSAV = LIST pointers ! implicit none integer in1 integer in2 integer io1 integer io2 integer lend(*) integer list(*) integer lp integer lp21 integer lph integer lpsav integer lptr(*) integer lstptr ! ! Test for IN1 and IN2 adjacent. ! lp = lstptr(lend(in1),in2,list,lptr) if ( abs ( list(lp) ) == in2 ) then lp21 = 0 return end if ! ! Delete IO2 as a neighbor of IO1. ! lp = lstptr(lend(io1),in2,list,lptr) lph = lptr(lp) lptr(lp) = lptr(lph) ! ! If IO2 is the last neighbor of IO1, make IN2 the last neighbor. ! if ( lend(io1) == lph ) then lend(io1) = lp end if ! ! Insert IN2 as a neighbor of IN1 following IO1 ! using the hole created above. ! lp = lstptr(lend(in1),io1,list,lptr) lpsav = lptr(lp) lptr(lp) = lph list(lph) = in2 lptr(lph) = lpsav ! ! Delete IO1 as a neighbor of IO2. ! lp = lstptr(lend(io2),in1,list,lptr) lph = lptr(lp) lptr(lp) = lptr(lph) ! ! If IO1 is the last neighbor of IO2, make IN1 the last neighbor. ! if ( lend(io2) == lph ) then lend(io2) = lp end if ! ! Insert IN1 as a neighbor of IN2 following IO2. ! lp = lstptr(lend(in2),io2,list,lptr) lpsav = lptr(lp) lptr(lp) = lph list(lph) = in1 lptr(lph) = lpsav lp21 = lph return end function swptst ( in1, in2, io1, io2, x, y ) !*****************************************************************************80 ! !! SWPTST applies the circumcircle test to a quadrilateral. ! ! Discussion: ! ! This function applies the circumcircle test to a quadri- ! lateral defined by a pair of adjacent triangles. The ! diagonal arc (shared triangle side) should be swapped for ! the other diagonl if and only if the fourth vertex is ! strictly interior to the circumcircle of one of the ! triangles (the decision is independent of the choice of ! triangle). Equivalently, the diagonal is chosen to maxi- ! mize the smallest of the six interior angles over the two ! pairs of possible triangles (the decision is for no swap ! if the quadrilateral is not strictly convex). ! ! When the four vertices are nearly cocircular (the ! neutral case), the preferred decision is no swap -- in ! order to avoid unnecessary swaps and, more important, to ! avoid cycling in subroutine OPTIM which is called by ! DELNOD and EDGE. Thus, a tolerance SWTOL (stored in ! SWPCOM by TRMESH or TRMSHR) is used to define 'nearness' ! to the neutral case. ! ! Modified: ! ! 16 June 2007 ! ! Author: ! ! Robert Renka, ! Department of Computer Science, ! University of North Texas, ! renka@cs.unt.edu ! ! Reference: ! ! Robert Renka, ! Algorithm 751: TRIPACK, ! A Constrained Two-Dimensional Delaunay Triangulation Package, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 1, 1996. ! ! Parameters: ! ! Input, integer IN1, IN2, IO1, IO2, the nodal indexes of the vertices of ! the quadrilateral. IO1-IO2 is the triangulation arc (shared triangle ! side) to be replaced by IN1-IN2 if the decision is to swap. The ! triples (IO1,IO2,IN1) and (IO2,IO1,IN2) must define triangles (be ! in counterclockwise order) on input. ! ! Input, real ( kind = 8 ) X(*), Y(*), the nodal coordinates. ! ! Output, logical SWPTST, .TRUE. if and only if the arc connecting ! IO1 and IO2 is to be replaced. ! ! Local parameters: ! ! DX11,DY11 = X,Y components of the vector IN1->IO1 ! DX12,DY12 = X,Y components of the vector IN1->IO2 ! DX22,DY22 = X,Y components of the vector IN2->IO2 ! DX21,DY21 = X,Y components of the vector IN2->IO1 ! SIN1 = Cross product of the vectors IN1->IO1 and ! IN1->IO2 -- proportional to sin(T1), where ! T1 is the angle at IN1 formed by the vectors ! COS1 = Inner product of the vectors IN1->IO1 and ! IN1->IO2 -- proportional to cos(T1) ! SIN2 = Cross product of the vectors IN2->IO2 and ! IN2->IO1 -- proportional to sin(T2), where ! T2 is the angle at IN2 formed by the vectors ! COS2 = Inner product of the vectors IN2->IO2 and ! IN2->IO1 -- proportional to cos(T2) ! SIN12 = SIN1*COS2 + COS1*SIN2 -- proportional to sin(T1+T2) ! implicit none real ( kind = 8 ) cos1 real ( kind = 8 ) cos2 real ( kind = 8 ) dx11 real ( kind = 8 ) dx12 real ( kind = 8 ) dx21 real ( kind = 8 ) dx22 real ( kind = 8 ) dy11 real ( kind = 8 ) dy12 real ( kind = 8 ) dy21 real ( kind = 8 ) dy22 integer in1 integer in2 integer io1 integer io2 real ( kind = 8 ) sin1 real ( kind = 8 ) sin12 real ( kind = 8 ) sin2 logical swptst real ( kind = 8 ) swtol real ( kind = 8 ) x(*) real ( kind = 8 ) y(*) ! ! Tolerance stored by TRMESH or TRMSHR. ! common /swpcom/ swtol ! ! Compute the vectors containing the angles T1 and T2. ! dx11 = x(io1) - x(in1) dx12 = x(io2) - x(in1) dx22 = x(io2) - x(in2) dx21 = x(io1) - x(in2) dy11 = y(io1) - y(in1) dy12 = y(io2) - y(in1) dy22 = y(io2) - y(in2) dy21 = y(io1) - y(in2) ! ! Compute inner products. ! cos1 = dx11 * dx12 + dy11 * dy12 cos2 = dx22 * dx21 + dy22 * dy21 ! ! The diagonals should be swapped iff 180 < (T1+T2) ! degrees. The following two tests ensure numerical ! stability: the decision must be FALSE when both ! angles are close to 0, and TRUE when both angles ! are close to 180 degrees. ! if ( 0.0D+00 <= cos1 .and. 0.0D+00 <= cos2 ) then swptst = .false. return end if if ( cos1 < 0.0D+00 .and. cos2 < 0.0D+00 ) then swptst = .true. return end if ! ! Compute vector cross products (Z-components). ! sin1 = dx11 * dy12 - dx12 * dy11 sin2 = dx22 * dy21 - dx21 * dy22 sin12 = sin1 * cos2 + cos1 * sin2 if ( -swtol <= sin12 ) then swptst = .false. else swptst = .true. end if return end subroutine trfind ( nst, px, py, n, x, y, list, lptr, lend, i1, i2, i3 ) !*****************************************************************************80 ! !! TRFIND locates a point relative to a triangulation. ! ! Discussion: ! ! This subroutine locates a point P relative to a triangu- ! lation created by subroutine TRMESH or TRMSHR. If P is ! contained in a triangle, the three vertex indexes are ! returned. Otherwise, the indexes of the rightmost and ! leftmost visible boundary nodes are returned. ! ! Modified: ! ! 16 June 2007 ! ! Author: ! ! Robert Renka, ! Department of Computer Science, ! University of North Texas, ! renka@cs.unt.edu ! ! Reference: ! ! Robert Renka, ! Algorithm 751: TRIPACK, ! A Constrained Two-Dimensional Delaunay Triangulation Package, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 1, 1996. ! ! Parameters: ! ! Input, integer NST, the index of a node at which TRFIND begins the ! search. Search time depends on the proximity of this node to P. ! ! Input, real ( kind = 8 ) PX, PY, the coordinates of the point P to be ! located. ! ! Input, integer N, the number of nodes in the triangulation. 3 <= N. ! ! Input, real ( kind = 8 ) X(N), Y(N), the coordinates of the nodes in ! the triangulation. ! ! Input, integer LIST(*), LPTR(*), LEND(N), the data structure defining ! the triangulation. Refer to subroutine TRMESH. ! ! Output, integer I1, I2, I3, nodal indexes, in counterclockwise order, ! of the vertices of a triangle containing P if P is contained in a ! triangle. If P is not in the convex hull of the nodes, I1 indexes ! the rightmost visible boundary node, I2 indexes the leftmost visible ! boundary node, and I3 = 0. Rightmost and leftmost are defined from ! the perspective of P, and a pair of points are visible from each ! other if and only if the line segment joining them intersects no ! triangulation arc. If P and all of the nodes lie on a common line, ! then I1 = I2 = I3 = 0 on output. ! ! Local parameters: ! ! B1,B2 = Unnormalized barycentric coordinates of P with respect ! to (N1,N2,N3) ! IX,IY,IZ = Integer seeds for JRAND ! LP = LIST pointer ! N0,N1,N2 = Nodes in counterclockwise order defining a ! cone (with vertex N0) containing P ! N1S,N2S = Saved values of N1 and N2 ! N3,N4 = Nodes opposite N1->N2 and N2->N1, respectively ! NB = Index of a boundary node -- first neighbor of ! NF or last neighbor of NL in the boundary traversal loops ! NF,NL = First and last neighbors of N0, or first ! (rightmost) and last (leftmost) nodes ! visible from P when P is exterior to the triangulation ! NP,NPP = Indexes of boundary nodes used in the boundary traversal loops ! XA,XB,XC = Dummy arguments for FRWRD ! YA,YB,YC = Dummy arguments for FRWRD ! XP,YP = Local variables containing the components of P ! implicit none integer n real ( kind = 8 ) b1 real ( kind = 8 ) b2 logical frwrd integer i1 integer i2 integer i3 integer, save :: ix = 1 integer, save :: iy = 2 integer, save :: iz = 3 integer jrand logical left integer lend(n) integer list(*) integer lp integer lptr(*) integer lstptr integer n0 integer n1 integer n1s integer n2 integer n2s integer n3 integer n4 integer nb integer nf integer nl integer np integer npp integer nst real ( kind = 8 ) px real ( kind = 8 ) py real ( kind = 8 ) store real ( kind = 8 ) x(n) real ( kind = 8 ) xa real ( kind = 8 ) xb real ( kind = 8 ) xc real ( kind = 8 ) xp real ( kind = 8 ) y(n) real ( kind = 8 ) ya real ( kind = 8 ) yb real ( kind = 8 ) yc real ( kind = 8 ) yp ! ! Statement function: ! ! FRWRD = TRUE iff C is forward of A->B iff B,A->C> >= 0. ! frwrd(xa,ya,xb,yb,xc,yc) = (xb-xa)*(xc-xa) + (yb-ya)*(yc-ya) >= 0.0D+00 ! ! Initialize variables. ! xp = px yp = py n0 = nst if ( n0 < 1 .or. n < n0 ) then n0 = jrand ( n, ix, iy, iz ) end if ! ! Set NF and NL to the first and last neighbors of N0, and ! initialize N1 = NF. ! 1 continue lp = lend(n0) nl = list(lp) lp = lptr(lp) nf = list(lp) n1 = nf ! ! Find a pair of adjacent neighbors N1,N2 of N0 that define ! a wedge containing P: P LEFT N0->N1 and P RIGHT N0->N2. ! if ( 0 < nl ) then go to 2 end if ! ! N0 is a boundary node. Test for P exterior. ! nl = -nl if ( .not. left ( x(n0), y(n0), x(nf), y(nf), xp, yp ) ) then nl = n0 go to 9 end if if ( .not. left(x(nl),y(nl),x(n0),y(n0),xp,yp) ) then nb = nf nf = n0 np = nl npp = n0 go to 11 end if go to 3 ! ! N0 is an interior node. Find N1. ! 2 continue do if ( left(x(n0),y(n0),x(n1),y(n1),xp,yp) ) then exit end if lp = lptr(lp) n1 = list(lp) if ( n1 == nl ) then go to 6 end if end do ! ! P is to the left of edge N0->N1. Initialize N2 to the ! next neighbor of N0. ! 3 continue lp = lptr(lp) n2 = abs ( list(lp) ) if ( .not. left(x(n0),y(n0),x(n2),y(n2),xp,yp) ) then go to 7 end if n1 = n2 if ( n1 /= nl ) then go to 3 end if if ( .not. left(x(n0),y(n0),x(nf),y(nf),xp,yp) ) then go to 6 end if if (xp == x(n0) .and. yp == y(n0)) then go to 5 end if ! ! P is left of or on edges N0->NB for all neighbors NB of N0. ! All points are collinear iff P is left of NB->N0 for ! all neighbors NB of N0. Search the neighbors of N0. ! NOTE: N1 = NL and LP points to NL. ! 4 continue if ( .not. left(x(n1),y(n1),x(n0),y(n0),xp,yp) ) then go to 5 end if lp = lptr(lp) n1 = abs ( list(lp) ) if ( n1 == nl ) then i1 = 0 i2 = 0 i3 = 0 return end if go to 4 ! ! P is to the right of N1->N0, or P=N0. Set N0 to N1 and start over. ! 5 continue n0 = n1 go to 1 ! ! P is between edges N0->N1 and N0->NF. ! 6 continue n2 = nf ! ! P is contained in the wedge defined by line segments ! N0->N1 and N0->N2, where N1 is adjacent to N2. Set ! N3 to the node opposite N1->N2, and save N1 and N2 to ! test for cycling. ! 7 continue n3 = n0 n1s = n1 n2s = n2 ! ! Top of edge hopping loop. Test for termination. ! 8 continue if ( left ( x(n1), y(n1), x(n2), y(n2), xp, yp ) ) then ! ! P LEFT N1->N2 and hence P is in (N1,N2,N3) unless an ! error resulted from floating point inaccuracy and ! collinearity. Compute the unnormalized barycentric ! coordinates of P with respect to (N1,N2,N3). ! b1 = (x(n3)-x(n2))*(yp-y(n2)) - (xp-x(n2))*(y(n3)-y(n2)) b2 = (x(n1)-x(n3))*(yp-y(n3)) - (xp-x(n3))*(y(n1)-y(n3)) if ( store ( b1 + 1.0D+00 ) >= 1.0D+00 .and. & store ( b2 + 1.0D+00 ) >= 1.0D+00 ) then go to 16 end if ! ! Restart with N0 randomly selected. ! n0 = jrand ( n, ix, iy, iz ) go to 1 end if ! ! Set N4 to the neighbor of N2 which follows N1 (node ! opposite N2->N1) unless N1->N2 is a boundary edge. ! lp = lstptr(lend(n2),n1,list,lptr) if ( list(lp) < 0 ) then nf = n2 nl = n1 go to 9 end if lp = lptr(lp) n4 = abs ( list(lp) ) ! ! Select the new edge N1->N2 which intersects the line ! segment N0-P, and set N3 to the node opposite N1->N2. ! if ( left(x(n0),y(n0),x(n4),y(n4),xp,yp) ) then n3 = n1 n1 = n4 n2s = n2 if (n1 /= n1s .and. n1 /= n0) go to 8 else n3 = n2 n2 = n4 n1s = n1 if ( n2 /= n2s .and. n2 /= n0 ) then go to 8 end if end if ! ! The starting node N0 or edge N1-N2 was encountered ! again, implying a cycle (infinite loop). Restart ! with N0 randomly selected. ! n0 = jrand ( n, ix, iy, iz ) go to 1 ! ! Boundary traversal loops. NL->NF is a boundary edge and ! P RIGHT NL->NF. Save NL and NF. 9 continue np = nl npp = nf ! ! Find the first (rightmost) visible boundary node NF. NB ! is set to the first neighbor of NF, and NP is the last neighbor. ! 10 continue lp = lend(nf) lp = lptr(lp) nb = list(lp) if ( .not. left(x(nf),y(nf),x(nb),y(nb),xp,yp) ) then go to 12 end if ! ! P LEFT NF->NB and thus NB is not visible unless an error ! resulted from floating point inaccuracy and collinear- ! ity of the 4 points NP, NF, NB, and P. ! 11 continue if ( frwrd(x(nf),y(nf),x(np),y(np),xp,yp) .or. & frwrd(x(nf),y(nf),x(np),y(np),x(nb),y(nb)) ) then i1 = nf go to 13 end if ! ! Bottom of loop. ! 12 continue np = nf nf = nb go to 10 ! ! Find the last (leftmost) visible boundary node NL. NB ! is set to the last neighbor of NL, and NPP is the first ! neighbor. ! 13 continue lp = lend(nl) nb = -list(lp) if ( .not. left(x(nb),y(nb),x(nl),y(nl),xp,yp) ) then go to 14 end if ! ! P LEFT NB->NL and thus NB is not visible unless an error ! resulted from floating point inaccuracy and collinear- ! ity of the 4 points P, NB, NL, and NPP. ! if ( frwrd(x(nl),y(nl),x(npp),y(npp),xp,yp) .or. & frwrd(x(nl),y(nl),x(npp),y(npp),x(nb),y(nb)) ) then go to 15 end if ! ! Bottom of loop. ! 14 continue npp = nl nl = nb go to 13 ! ! NL is the leftmost visible boundary node. ! 15 continue i2 = nl i3 = 0 return ! ! P is in the triangle (N1,N2,N3). ! 16 continue i1 = n1 i2 = n2 i3 = n3 return end subroutine trlist ( ncc, lcc, n, list, lptr, lend, nrow, nt, ltri, lct, ier ) !*****************************************************************************80 ! !! TRLIST converts a triangulation to triangle list form. ! ! Discussion: ! ! This subroutine converts a triangulation data structure ! from the linked list created by subroutine TRMESH or ! TRMSHR to a triangle list. ! ! Modified: ! ! 16 June 2007 ! ! Author: ! ! Robert Renka, ! Department of Computer Science, ! University of North Texas, ! renka@cs.unt.edu ! ! Reference: ! ! Robert Renka, ! Algorithm 751: TRIPACK, ! A Constrained Two-Dimensional Delaunay Triangulation Package, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 1, 1996. ! ! Parameters: ! ! Input, integer NCC, the number of constraints. NCC >= 0. ! ! Input, integer LCC(*), list of constraint curve starting indexes (or ! dummy array of length 1 if NCC = 0). Refer to subroutine ADDCST. ! ! Input, integer N, the number of nodes in the triangulation. N >= 3. ! ! Input, integer LIST(*), LPTR(*), LEND(N), linked list data structure ! defining the triangulation. Refer to subroutine TRMESH. ! ! Input, integer NROW, the number of rows (entries per triangle) ! reserved for the triangle list LTRI. The value must be 6 if only ! the vertex indexes and neighboring triangle indexes are to be ! stored, or 9 if arc indexes are also to be assigned and stored. ! Refer to LTRI. ! ! Input, integer LTRI(NROW*NT), where NT is at most 2N-5. (A sufficient ! length is 12 * N if NROW=6 or 18*N if NROW=9.) ! ! Output, integer NT, the number of triangles in the triangulation unless ! IER /= 0, in which case NT = 0. NT = 2N - NB- 2, where NB is the number ! of boundary nodes. ! ! Output, integer LTRI(NROW,NT), whose J-th column contains the vertex nodal ! indexes (first three rows), neighboring triangle indexes (second three ! rows), and, if NROW = 9, arc indexes (last three rows) associated with ! triangle J for J = 1,...,NT. The vertices are ordered counterclockwise ! with the first vertex taken to be the one with smallest index. Thus, ! LTRI(2,J) and LTRI(3,J) are larger than LTRI(1,J) and index adjacent ! neighbors of node LTRI(1,J). For I = 1,2,3, LTRI(I+3,J) and LTRI(I+6,J) ! index the triangle and arc, respectively, which are opposite (not shared ! by) node LTRI(I,J), with LTRI(I+3,J) = 0 if LTRI(I+6,J) indexes a boundary ! arc. Vertex indexes range from 1 to N, triangle indexes from 0 to NT, ! and, if included, arc indexes from 1 to NA = NT+N-1. The triangles are ! ordered on first (smallest) vertex indexes, except that the sets of ! constraint triangle (triangles contained in the closure of a constraint ! region) follow the non-constraint triangles. ! ! Output, integer LCT(NCC), containing the triangle index of the first ! triangle of constraint J in LCT(J). Thus, the number of non-constraint ! triangles is LCT(1)-1, and constraint J contains LCT(J+1)-LCT(J) ! triangles, where LCT(NCC+1) = NT+1. ! ! Output, integer IER = Error indicator. ! 0, if no errors were encountered. ! 1, if NCC, N, NROW, or an LCC entry is outside its valid range on input. ! 2, if the triangulation data structure (LIST,LPTR,LEND) is invalid. ! ! Local Parameters: ! ! ARCS = TRUE iff arc indexes are to be stored. ! KA,KT = Numbers of currently stored arcs and triangles. ! N1ST = Starting index for the loop on nodes (N1ST = 1 on ! pass 1, and N1ST = LCC1 on pass 2). ! NM2 = Upper bound on candidates for N1. ! PASS2 = TRUE iff constraint triangles are to be stored. ! implicit none integer n integer nrow logical arcs logical cstri integer i integer i1 integer i2 integer i3 integer ier integer isv integer j integer jlast integer ka integer kn integer kt integer l integer lcc(*) integer lcc1 integer lct(*) integer lend(n) integer list(*) integer lp integer lp2 integer lpl integer lpln1 integer lptr(*) integer ltri(nrow,*) integer n1 integer n1st integer n2 integer n3 integer ncc integer nm2 integer nn integer nt logical pass2 ! ! Test for invalid input parameters and store the index ! LCC1 of the first constraint node (if any). ! nn = n if ( ncc < 0 .or. ( nrow /= 6 .and. nrow /= 9 ) ) then nt = 0 ier = 1 return end if lcc1 = nn+1 if (ncc == 0) then if ( nn < 3 ) then nt = 0 ier = 1 return end if else do i = ncc, 1, -1 if ( lcc1 - lcc(i) < 3 ) then nt = 0 ier = 1 return end if lcc1 = lcc(i) end do if ( lcc1 < 1 ) then nt = 0 ier = 1 return end if end if ! ! Initialize parameters for loop on triangles KT = (N1,N2, ! N3), where N1 < N2 and N1 < N3. This requires two ! passes through the nodes with all non-constraint ! triangles stored on the first pass, and the constraint ! triangles stored on the second. ! arcs = nrow == 9 ka = 0 kt = 0 n1st = 1 nm2 = nn - 2 pass2 = .false. ! ! Loop on nodes N1: ! J = constraint containing N1, ! JLAST = last node in constraint J. ! 2 continue j = 0 jlast = lcc1 - 1 do n1 = n1st, nm2 if ( jlast < n1 ) then ! ! N1 is the first node in constraint J+1. Update J and ! JLAST, and store the first constraint triangle index ! if in pass 2. ! j = j + 1 if ( j < ncc ) then jlast = lcc(j+1) - 1 else jlast = nn end if if ( pass2 ) then lct(j) = kt + 1 end if end if ! ! Loop on pairs of adjacent neighbors (N2,N3). LPLN1 points ! to the last neighbor of N1, and LP2 points to N2. ! lpln1 = lend(n1) lp2 = lpln1 3 continue lp2 = lptr(lp2) n2 = list(lp2) lp = lptr(lp2) n3 = abs ( list(lp) ) if ( n2 < n1 .or. n3 < n1 ) then go to 10 end if ! ! (N1,N2,N3) is a constraint triangle iff the three nodes ! are in the same constraint and N2 < N3. Bypass con- ! straint triangles on pass 1 and non-constraint triangles ! on pass 2. ! cstri = n1 >= lcc1 .and. n2 < n3 .and. n3 <= jlast if ( ( cstri .and. .not. pass2 ) .or. & ( .not. cstri .and. pass2 ) ) then go to 10 end if ! ! Add a new triangle KT = (N1,N2,N3). ! kt = kt + 1 ltri(1,kt) = n1 ltri(2,kt) = n2 ltri(3,kt) = n3 ! ! Loop on triangle sides (I1,I2) with neighboring triangles ! KN = (I1,I2,I3). ! do i = 1,3 if ( i == 1 ) then i1 = n3 i2 = n2 else if ( i == 2 ) then i1 = n1 i2 = n3 else i1 = n2 i2 = n1 end if ! ! Set I3 to the neighbor of I1 which follows I2 unless ! I2->I1 is a boundary arc. ! lpl = lend(i1) lp = lptr(lpl) 4 continue if (list(lp) == i2) then go to 5 end if lp = lptr(lp) if ( lp /= lpl ) then go to 4 end if ! ! I2 is the last neighbor of I1 unless the data structure ! is invalid. Bypass the search for a neighboring ! triangle if I2->I1 is a boundary arc. ! if ( abs ( list(lp) ) /= i2 ) then go to 13 end if kn = 0 if (list(lp) < 0) then go to 8 end if ! ! I2->I1 is not a boundary arc, and LP points to I2 as ! a neighbor of I1. ! 5 continue lp = lptr(lp) i3 = abs ( list(lp) ) ! ! Find L such that LTRI(L,KN) = I3 (not used if KN > KT), ! and permute the vertex indexes of KN so that I1 is ! smallest. ! if ( i1 < i2 .and. i1 < i3 ) then l = 3 else if (i2 < i3) then l = 2 isv = i1 i1 = i2 i2 = i3 i3 = isv else l = 1 isv = i1 i1 = i3 i3 = i2 i2 = isv end if ! ! Test for KN > KT (triangle index not yet assigned). ! if ( i1 > n1 .and. .not. pass2 ) then go to 9 end if ! ! Find KN, if it exists, by searching the triangle list in ! reverse order. ! do kn = kt-1,1,-1 if ( ltri(1,kn) == i1 .and. ltri(2,kn) == & i2 .and. ltri(3,kn) == i3 ) then go to 7 end if end do go to 9 ! ! Store KT as a neighbor of KN. ! 7 continue ltri(l+3,kn) = kt ! ! Store KN as a neighbor of KT, and add a new arc KA. ! 8 continue ltri(i+3,kt) = kn if (arcs) then ka = ka + 1 ltri(i+6,kt) = ka if ( kn /= 0 ) then ltri(l+6,kn) = ka end if end if 9 continue end do ! ! Bottom of loop on triangles. ! 10 continue if ( lp2 /= lpln1 ) then go to 3 end if end do ! ! Bottom of loop on nodes. ! if ( .not. pass2 .and. 0 < ncc ) then pass2 = .true. n1st = lcc1 go to 2 end if ! ! No errors encountered. ! nt = kt ier = 0 return ! ! Invalid triangulation data structure: I1 is a neighbor of ! I2, but I2 is not a neighbor of I1. ! 13 continue nt = 0 ier = 2 return end subroutine trlprt ( ncc, lct, n, x, y, nrow, nt, ltri, prntx ) !*****************************************************************************80 ! !! TRLPRT prints the triangles in a triangulation. ! ! Discussion: ! ! Given a triangulation of a set of points in the plane, ! this subroutine prints the triangle list created by ! subroutine TRLIST and, optionally, the nodal coordinates ! on logical unit LOUT. The numbers of boundary nodes, ! triangles, and arcs, and the constraint region triangle ! indexes, if any, are also printed. ! ! All parameters other than PRNTX should be ! unaltered from their values on output from TRLIST. ! ! Modified: ! ! 16 June 2007 ! ! Author: ! ! Robert Renka, ! Department of Computer Science, ! University of North Texas, ! renka@cs.unt.edu ! ! Reference: ! ! Robert Renka, ! Algorithm 751: TRIPACK, ! A Constrained Two-Dimensional Delaunay Triangulation Package, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 1, 1996. ! ! Parameters: ! ! Input, integer NCC, the number of constraints. ! ! Input, integer LCT(NCC), the list of constraint triangle starting ! indexes (or dummy array of length 1 if NCC = 0). ! ! Input, integer N, the number of nodes in the triangulation. ! 3 <= N <= 9999. ! ! Input, real ( kind = 8 ) X(N), Y(N), the coordinates of the nodes in the ! triangulation; not used unless PRNTX = TRUE. ! ! Input, integer NROW, the number of rows (entries per triangle) ! reserved for the triangle list LTRI. The value must be 6 if only ! the vertex indexes and neighboring triangle indexes are stored, ! or 9 if arc indexes are also stored. ! ! Input, integer NT, the number of triangles in the triangulation. ! 1 <= NT <= 9999. ! ! Input, integer LTRI(NROW,NT), array whose J-th column contains ! the vertex nodal indexes (first three rows), neighboring triangle ! indexes (second three rows), and, if NROW = 9, arc indexes (last ! three rows) associated with triangle J for J = 1,...,NT. ! ! Input, logical PRNTX, is TRUE if and only if X and Y are to be printed. ! ! Local parameters: ! ! I = DO-loop, nodal index, and row index for LTRI ! K = DO-loop and triangle index ! LUN = Logical unit number for output ! NA = Number of triangulation arcs ! NB = Number of boundary nodes ! NL = Number of lines printed on the current page ! NLMAX = Maximum number of print lines per page ! NMAX = Maximum value of N and NT (4-digit format) !