#!/bin/sh
# This is a shell archive, meaning:
# 1. Remove everything above the #!/bin/sh line.
# 2. Save the resulting text in a file.
# 3. Execute the file with /bin/sh (not csh) to create the files:
# READ_ME
# Makefile
# acchk.r
# addpt.r
# adjchk.r
# binsrt.r
# circen.r
# cross.r
# delet.r
# delet1.r
# delout.r
# delseg.r
# dirout.r
# dirseg.r
# driver.r
# dumpts.r
# eldup.r
# initad.r
# inside.r
# insrt.r
# insrt1.r
# locn.r
# master.r
# pred.r
# qtest.r
# qtest1.r
# stoke.r
# succ.r
# swap.r
# testeq.r
# triar.r
# trifnd.r
# err.list
# ex.dat
# ex.out
# This archive created: Mon May 6 14:32:16 1996
export PATH; PATH=/bin:$PATH
if test -f 'READ_ME'
then
echo shar: over-writing existing file "'READ_ME'"
fi
cat << \SHAR_EOF > 'READ_ME'
*****************************
RATFOR CODE TO COMPUTE THE DIRICHLET TESSELLATION AND DELAUNAY
TRIANGULATION OF A SET OF DATA POINTS AND POSSIBLY A SUPERIMPOSED GRID
OF DUMMY POINTS.
THE TRIANGULATION IS CONSTRUCTED WITH REPECT TO THE WHOLE PLANE
BY SUSPENDING IT FROM IDEAL POINTS AT INFINITY.
Copyright (C) 1996 by T. Rolf Turner
Permission to use, copy, modify, and distribute this software and
its documentation for any purpose and without fee is hereby
granted, provided that the above copyright notice appear in all
copies and that both that copyright notice and this permission
notice appear in supporting documentation.
ORIGINALLY PROGRAMMED BY: Rolf Turner in 1987/88, while with the
Division of Mathematics and Statistics, CSIRO, Sydney, Australia.
Re-programmed by Rolf Turner while visiting the University of
Western Australia, May 1995. Revision completed at the University
of New Brunswick, May 1996.
Current address of the author:
Department of Mathematics and Statistics,
University of New Brunswick,
P.O. Box 4400, Fredericton, New Brunswick,
Canada E3B 5A3
Email:
rolf@math.unb.ca
The author gratefully acknowledges the contributions, assistance, and
guidance of Mark Berman, of D.M.S., CSIRO, in collaboration with whom
this project was originally undertaken. The author also acknowledges
much useful advice from Adrian Baddeley, formerly of D.M.S. CSIRO (now
Professor of Statistics at the University of Western Australia), and
a very useful conversation with Daryl Tingley of the Department of
Mathematics, University of New Brunswick. A data set provided by
Alan Johnson was extremely useful in tracking down some errors in the
original code.
*****************************
This archive contains:
(a) This READ_ME file.
(b) A Makefile.
(c) Many, many *.r files of ratfor subroutines, and a driver.
(d) A file err.list, containing a list of meanings of possible
error numbers which could be returned. NONE of these
errors should ever actually happen except for errors 4, 14,
and 15. These relate to insufficient dimensioning, and
if they occur, you should increase the relevant dimension
(either MADJ or NSGS) in the ``define'' statements in driver.r
and try again.
(e) A sample file of data, ex.dat. (Note that one of the
points in ex.dat lies outside the rectangular window
specified in the example driver.r file. Hence it gets
discarded by the program.)
(f) The output, ex.out, from running the program on this data
set, with parameters as set in the example driver.r
file. In particular the rectangular window was specified
to be (0,10,0,10) and the dummy points were specified to
coincide with the corners of this window.
The ratfor code is ponderous and quite possibly kludgy -- i.e. a good
programmer could make it ***much*** more efficient I'm sure. It
contains all sorts of checking for anomalous situations which probably
can/will never occur. These checks basically reflect my pessimism and
fervent belief in Murphy's Law.
The program was also designed with a particular application in mind,
in which we wished to superimpose a grid of dummy points onto the
actual data points which we were triangulating. This fact adds
slightly to the complication of the code.
Here follows a brief description of the package:
(1) The code computes the Delaunay Triangulation (and hence the
Dirichlet Tesselation) of a planar point set according to the second
(iterative) algorithm of Lee and Schacter, International Journal of
Computer and Information Sciences, Vol. 9, No. 3, 1980, pages 219 to
242.
The triangulation is made to be
**** with respect to the whole plane ****
by `suspending' it from `ideal' points (-Inf,-Inf), (Inf,-Inf)
(Inf,Inf), and (-Inf,Inf).
It is also enclosed in a finite rectangle
with corners
(xmin,ymax) * ------------------------ * (xmax,ymax)
| |
| |
| |
| |
| |
(xmin,ymin) * ------------------------ * (xmax,ymin)
The boundaries of this rectangle truncate some Dirichlet tiles, in
particular any infinite ones. This rectangle is referred to
elsewhere as `the' rectangular window.
===
(2) To use the program, edit the file ``driver.r'' and modify
the first 11 define statements to conform to your data set.
The quantities defined are:
N --- the number of (non-dummy) points in the data set
NDX --- The x-dimension of a rectangular grid of dummy points.
NDY --- The y-dimension of a rectangular grid of dummy points.
If either NDX or NDY is 0, no dummy points are added.
MADJ --- The maximum number of adjacencies in the triangulation
at any stage of its construction. Try the maximum of
20 and ceiling(3*sqrt(NTOT))) where NTOT = N+NDX*NDY+4
SORT --- Logical argument; if TRUE the data (including dummy
points) are sorted into a sequence of "bins" prior to
triangulation; this makes the algorithm slightly more
efficient. Normally one would set sort equal to FALSE only
if one wished to observe some of the fine detail of the way
in which adding a point to a data set affected the
triangulation, and therefore wished to make sure that the
point in question was added last. Essentially this argument
would get set to FALSE only in a de-bugging process.
EPS --- A value of epsilon used in testing whether a quantity is
zero, mainly in the context of whether points are collinear.
If anomalous errors arise, it is possible that these may
averted by adjusting the value of eps upward or downward.
FRAC --- A value specifying the tolerance used in eliminating
duplicate points. Points are considered duplicates if
abs(x1-x2) < frac*(xmax-xmin) AND abs(y1-y2) <
frac*(ymax-ymin).
RW1, RW2, RW3, and RW4 --- the coordinates of the corners of the
rectangular window enclosing the triangulation, in the order
(xmin, xmax, ymin, ymax). Any data points outside this
window are discarded.
Then type `make driver'. (The `Makefile' may need to be modified
for compatibility with your local system.) This will produce an
executable file (program) called `driver'.
(3) To execute the program type
driver < datafile > outputfile
(4) The file `datafile' should contain the coordinates of the point
set to be triangulated in (x,y) pairs, each pair on a separate
line.
(5) The output file contains:
(a) The total area of the Delaunay triangulation and
the total area of the Dirichlet tessellation. (These
will be the same if ndx >= 2 and ndy >= 2, in which
case the area is the area of `the' rectangular window.)
(b) The ``Delaunay segments''. In this section the first 4
entries of each row are the coordinates of the points joined
by an edge of a Delaunay triangle, in the order (x1,y1,x2,y2).
The last two entries are the indices of the two points which
are joined.
(c) The ``Dirichlet segments''. In this section the first 4
entries of each row are the coordinates of the endpoints of
one of the edges of a Dirichlet tile, in the order (x1,y1,x2,y2).
The fifth and sixth entries are the indices of the two
points, in the set being triangulated, which are separated by
that edge. The seventh and eighth entries are indicators of
``boundary points'' (points on the boundary of the
rectangular window). The seventh entry is 1 if the first
endpoint of the edge of the tile is a boundary point, and is
0 otherwise. Likewise for the eighth entry and the second
endpoint of the edge.
(d) The ``Delaunay summary''. In this section each row has
four entries: the x and y coordinates of the points in the
set being triangulated (including dummy points), the number
of Delaunay triangles emmanating from the given point, and
1/3 times the total area of these Delaunay triangles. Note
that the factor of 1/3 arises because each triangle occurs
three times --- once for each corner.
Note that the total area of the Delaunay triangulation is
equal to the sum of the last column of the ``Delaunay summary''.
(e) The ``Dirichlet summary''. In this section each row has
three entries: the number of sides --- within the rectangular
window --- of the Dirichlet tile surrounding a given point,
the number of points in which the Dirichlet tile intersects
the boundary of the rectangular window, and the area of the
Dirichlet tile surrounding the point. (The ``given point''
is the same as that in the corresponding row of the
``Delaunay summary''.)
Note that the total area of the Dirichlet tessellation is
equal to the sum of the last column of the ``Dirichlet summary''.
(f) An annotation saying how many of the rows (the first
n) correspond to ``real'' data points, and how many (the
last ndum) correspond to dummy points. Note that any
duplicated points, including dummy points which duplicate
real data points, will have been eliminated.
(6) Remove unwanted *.o files (if you're not going to be remaking
for a while) by typing ``make clean''.
SHAR_EOF
if test -f 'Makefile'
then
echo shar: over-writing existing file "'Makefile'"
fi
cat << \SHAR_EOF > 'Makefile'
F77 = f77
FFLAGS =
LFLAGS = -o
OBJECTS = acchk.o addpt.o adjchk.o binsrt.o circen.o cross.o \
delet.o delet1.o delout.o delseg.o dirout.o dirseg.o driver.o dumpts.o \
eldup.o initad.o inside.o insrt.o insrt1.o locn.o master.o pred.o \
qtest.o qtest1.o stoke.o succ.o swap.o testeq.o triar.o trifnd.o
FILES = READ_ME Makefile acchk.r addpt.r adjchk.r binsrt.r circen.r \
cross.r delet.r delet1.r delout.r delseg.r dirout.r dirseg.r \
driver.r dumpts.r eldup.r initad.r inside.r insrt.r insrt1.r \
locn.r master.r pred.r qtest.r qtest1.r stoke.r succ.r swap.r \
testeq.r triar.r trifnd.r err.list ex.dat ex.out
driver: $(OBJECTS)
$(F77) $(LFLAGS) driver $(OBJECTS)
clean:
/bin/rm -f $(OBJECTS) shark
shark: $(FILES)
shar $(FILES) > shark
SHAR_EOF
if test -f 'acchk.r'
then
echo shar: over-writing existing file "'acchk.r'"
fi
cat << \SHAR_EOF > 'acchk.r'
subroutine acchk(i,j,k,anticl,x,y,ntot,eps)
# Check whether vertices i, j, k, are in anti-clockwise order.
# Called by locn, qtest, qtest1.
implicit double precision(a-h,o-z)
dimension x(-3:ntot), y(-3:ntot), xt(3), yt(3)
logical anticl
# Create indicator telling which of i, j, and k are ideal points.
if(i<=0) i1 = 1
else i1 = 0
if(j<=0) j1 = 1
else j1 = 0
if(k<=0) k1 = 1
else k1 = 0
ijk = i1*4+j1*2+k1
# Get the coordinates of vertices i, j, and k. (Pseudo-coordinates for
# any ideal points.)
xt(1) = x(i)
yt(1) = y(i)
xt(2) = x(j)
yt(2) = y(j)
xt(3) = x(k)
yt(3) = y(k)
# Get the ``normalized'' cross product.
call cross(xt,yt,ijk,cprd)
# If cprd is positive then (ij-cross-ik) is directed ***upwards***
# and so i, j, k, are in anti-clockwise order; else not.
if(cprd > eps) anticl = .true.
else anticl = .false.
return
end
SHAR_EOF
if test -f 'addpt.r'
then
echo shar: over-writing existing file "'addpt.r'"
fi
cat << \SHAR_EOF > 'addpt.r'
subroutine addpt(j,nadj,madj,x,y,ntot,eps,nerror)
# Add point j to the triangulation.
# Called by master, dirseg.
implicit double precision(a-h,o-z)
dimension nadj(-3:ntot,0:madj), x(-3:ntot), y(-3:ntot)
logical didswp
# Put the new point in, joined to the vertices of its
# enclosing triangle.
call initad(j,nadj,madj,x,y,ntot,eps,nerror)
if(nerror > 0) return
# Look at each `gap', i.e. pair of adjacent segments
# emanating from the new point; they form two sides of a
# quadrilateral; see whether the extant diagonal of this
# quadrilateral should be swapped with its alternate
# (according to the LOP: local optimality principle).
now = nadj(j,1)
nxt = nadj(j,2)
ngap = 0
repeat {
call swap(j,now,nxt,didswp,nadj,madj,x,y,ntot,eps,nerror)
if(nerror > 0) return
n = nadj(j,0)
if(!didswp) { # If no swap of diagonals
now = nxt # move to the next gap.
ngap = ngap+1
}
call succ(nxt,j,now,nadj,madj,ntot,nerror)
if(nerror > 0) return
}
until(ngap==n)
return
end
SHAR_EOF
if test -f 'adjchk.r'
then
echo shar: over-writing existing file "'adjchk.r'"
fi
cat << \SHAR_EOF > 'adjchk.r'
subroutine adjchk(i,j,adj,nadj,madj,ntot,nerror)
# Check if vertices i and j are adjacent.
# Called by insrt, delet, trifnd, swap, delseg, dirseg.
dimension nadj(-3:ntot,0:madj)
logical adj
nerror = -1
# Check if j is in the adjacency list of i.
adj = .false.
ni = nadj(i,0)
if(ni>0) {
do k = 1,ni {
if(j==nadj(i,k)) {
adj = .true.
break
}
}
}
# Check if i is in the adjacency list of j.
nj = nadj(j,0)
if(nj>0) {
do k = 1,nj {
if(i==nadj(j,k)) {
if(adj) return # Have j in i's list and i in j's.
else {
nerror = 1
return
}
}
}
}
# If we get to here i is not in j's list.
if(adj) { # If adj is true, then j IS in i's list.
nerror = 1
return
}
return
end
SHAR_EOF
if test -f 'binsrt.r'
then
echo shar: over-writing existing file "'binsrt.r'"
fi
cat << \SHAR_EOF > 'binsrt.r'
subroutine binsrt(x,y,ntot,rw,npd,ind,tx,ty,ilst,nerror)
# Sort the data points into bins.
# Called by master.
implicit double precision(a-h,o-z)
dimension x(-3:ntot), y(-3:ntot), tx(npd), ty(npd)
dimension ind(npd), ilst(npd)
dimension rw(4)
nerror = -1
kdiv = 1+dble(npd)**0.25 # Round high.
xkdiv = dble(kdiv)
# Dig out the corners of the rectangular window.
xmin = rw(1)
xmax = rw(2)
ymin = rw(3)
ymax = rw(4)
w = xmax-xmin
h = ymax-ymin
# Number of bins is to be approx. sqrt(npd); thus number of subdivisions
# on each side of rectangle is approx. npd**(1/4).
dw = w/xkdiv
dh = h/xkdiv
# The width of each bin is dw; the height is dh. We shall move across
# the rectangle from left to right, then up, then back from right to
# left, then up, .... Note that kx counts the divisions from the left,
# ky counts the divisions from the bottom; kx is incremented by ink, which
# is +/- 1 and switches sign when ky is incremented; ky is always
# incremented by 1.
kx = 1
ky = 1
ink = 1
k = 0
do i = 1,npd { ilst(i) = 0 } # Keeps a list of those points already added
while(ky<=kdiv) { # to the new list.
do i = 1,npd {
if(ilst(i)==1) next # The i-th point has already been added
# to the new list.
# If the i-th point is in the current bin, add it to the list.
xt = x(i)
yt = y(i)
ix = 1+(xt-xmin)/dw
if(ix>kdiv) ix = kdiv
jy = 1+(yt-ymin)/dh
if(jy>kdiv) jy = kdiv
if(ix==kx&jy==ky) {
k = k+1
ind(i) = k # Index i is the pos'n. of (x,y) in the
tx(k) = xt # old list; k is its pos'n. in the new one.
ty(k) = yt
ilst(i) = 1 # Cross the i-th point off the old list.
}
}
# Move to the next bin.
kc = kx+ink
if((1<=kc)&(kc<=kdiv)) kx = kc
else {
ky = ky+1
ink = -ink
}
}
# Check that all points from old list have been added to the new,
# with no spurious additions.
if(k!=npd) {
nerror = 2
return
}
# Copy the new sorted vector back on top of the old ones.
do i = 1,npd {
x(i) = tx(i)
y(i) = ty(i)
}
return
end
SHAR_EOF
if test -f 'circen.r'
then
echo shar: over-writing existing file "'circen.r'"
fi
cat << \SHAR_EOF > 'circen.r'
subroutine circen(i,j,k,x0,y0,x,y,ntot,eps,collin,nerror)
# Find the circumcentre (x0,y0) of the triangle with
# vertices (x(i),y(i)), (x(j),y(j)), (x(k),y(k)).
# Called by qtest1, dirseg, dirout.
implicit double precision(a-h,o-z)
dimension x(-3:ntot), y(-3:ntot), xt(3), yt(3)
logical collin
nerror = -1
# Get the coordinates.
xt(1) = x(i)
yt(1) = y(i)
xt(2) = x(j)
yt(2) = y(j)
xt(3) = x(k)
yt(3) = y(k)
# Check for collinearity
ijk = 0
call cross(xt,yt,ijk,cprd)
if(abs(cprd) < eps) collin = .true.
else collin = .false.
# Form the vector u from i to j, and the vector v from i to k,
# and normalize them.
a = x(j) - x(i)
b = y(j) - y(i)
c = x(k) - x(i)
d = y(k) - y(i)
c1 = sqrt(a*a+b*b)
c2 = sqrt(c*c+d*d)
a = a/c1
b = b/c1
c = c/c2
d = d/c2
# If the points are collinear, make sure that they're in the right
# order --- i between j and k.
if(collin) {
alpha = a*c+b*d
# If they're not in the right order, bring things to
# a shuddering halt.
if(alpha>0) {
nerror = 3
return
}
# Collinear, but in the right order; think of this as meaning
# that the circumcircle in question has infinite radius.
return
}
# Not collinear; go ahead, make my circumcentre. (First, form
# the cross product of the ***unit*** vectors, instead of the
# ``normalized'' cross product produced by ``cross''.)
crss = a*d - b*c
x0 = x(i) + 0.5*(c1*d - c2*b)/crss
y0 = y(i) + 0.5*(c2*a - c1*c)/crss
return
end
SHAR_EOF
if test -f 'cross.r'
then
echo shar: over-writing existing file "'cross.r'"
fi
cat << \SHAR_EOF > 'cross.r'
subroutine cross(x,y,ijk,cprd)
implicit double precision(a-h,o-z)
dimension x(3), y(3)
# Calculates a ``normalized'' cross product of the vectors joining
# [x(1),y(1)] to [x(2),y(2)] and [x(3),y(3)] respectively.
# The normalizatiohn consists in dividing by the square of the
# shortest of the three sides of the triangle. This normalization is
# for the purposes of testing for collinearity; if the result is less
# than epsilon, the smallest of the sines of the angles is less than
# epsilon.
# Set constants
zero = 0.d0
one = 1.d0
two = 2.d0
four = 4.d0
# Adjust the coordinates depending upon which points are ideal,
# and calculate the squared length of the shortest side.
switch(ijk) {
case 0: # No ideal points; no adjustment necessary.
smin = -one
do i = 1,3 {
ip = i+1
if(ip==4) ip = 1
a = x(ip) - x(i)
b = y(ip) - y(i)
s = a*a+b*b
if(smin < zero | s < smin) smin = s
}
case 1: # Only k ideal.
x(2) = x(2) - x(1)
y(2) = y(2) - y(1)
x(1) = zero
y(1) = zero
cn = sqrt(x(2)**2+y(2)**2)
x(2) = x(2)/cn
y(2) = y(2)/cn
smin = one
case 2: # Only j ideal.
x(3) = x(3) - x(1)
y(3) = y(3) - y(1)
x(1) = zero
y(1) = zero
cn = sqrt(x(3)**2+y(3)**2)
x(3) = x(3)/cn
y(3) = y(3)/cn
smin = one
case 3: # Both j and k ideal (i not).
x(1) = zero
y(1) = zero
smin = 2
case 4: # Only i ideal.
x(3) = x(3) - x(2)
y(3) = y(3) - y(2)
x(2) = zero
y(2) = zero
cn = sqrt(x(3)**2+y(3)**2)
x(3) = x(3)/cn
y(3) = y(3)/cn
smin = one
case 5: # Both i and k ideal (j not).
x(2) = zero
y(2) = zero
smin = two
case 6: # Both i and j ideal (k not).
x(3) = zero
y(3) = zero
smin = two
case 7: # All three points ideal; no adjustment necessary.
smin = four
}
a = x(2)-x(1)
b = y(2)-y(1)
c = x(3)-x(1)
d = y(3)-y(1)
cprd = (a*d - b*c)/smin
return
end
SHAR_EOF
if test -f 'delet.r'
then
echo shar: over-writing existing file "'delet.r'"
fi
cat << \SHAR_EOF > 'delet.r'
subroutine delet(i,j,nadj,madj,ntot,nerror)
# Delete i and j from each other's adjacency lists.
# Called by initad, swap.
implicit double precision(a-h,o-z)
dimension nadj(-3:ntot,0:madj)
logical adj
# First check that they're IN each other's lists.
call adjchk(i,j,adj,nadj,madj,ntot,nerror)
if(nerror > 0) return
# Then do the actual deletion if they are.
if(adj) {
call delet1(i,j,nadj,madj,ntot)
call delet1(j,i,nadj,madj,ntot)
}
return
end
SHAR_EOF
if test -f 'delet1.r'
then
echo shar: over-writing existing file "'delet1.r'"
fi
cat << \SHAR_EOF > 'delet1.r'
subroutine delet1(i,j,nadj,madj,ntot)
# Delete j from the adjacency list of i.
# Called by delet.
implicit double precision(a-h,o-z)
dimension nadj(-3:ntot,0:madj)
n = nadj(i,0)
do k = 1,n {
if(nadj(i,k)==j) { # Find j in the list;
# then move everything back one notch.
do kk = k,n-1 { nadj(i,kk) = nadj(i,kk+1) }
nadj(i,n) = 0
nadj(i,0) = n-1
return
}
}
end
SHAR_EOF
if test -f 'delout.r'
then
echo shar: over-writing existing file "'delout.r'"
fi
cat << \SHAR_EOF > 'delout.r'
subroutine delout(delsum,nadj,madj,x,y,ntot,npd1,npd,ind,nerror)
# Put a summary of the Delaunay triangles with a vertex at point i,
# for i = 1, ..., npd, into the array delsum. Do this in the original
# order of the points, not the order into which they have been
# bin-sorted.
# Called by master.
implicit double precision(a-h,o-z)
dimension nadj(-3:ntot,0:madj), x(-3:ntot), y(-3:ntot)
dimension delsum(npd1,4), ind(npd1)
do i1 = 1,npd {
area = 0. # Initialize area of polygon consisting of triangles
# with a vertex at point i.
# Get the point number, its coordinates and the number of
# (real) triangles emanating from it.
i = ind(i1)
np = nadj(i,0)
xi = x(i)
yi = y(i)
npt = np
do k = 1,np {
kp = k+1
if(kp>np) kp = 1
if(nadj(i,k)<=0|nadj(i,kp)<=0) npt = npt-1
}
# For each point in the adjacency list of point i, find its
# successor, and the area of the triangle determined by these
# three points.
do j1 = 1,np {
j = nadj(i,j1)
if(j<=0) next
xj = x(j)
yj = y(j)
call succ(k,i,j,nadj,madj,ntot,nerror)
if(nerror > 0) return
if(k<=0) next
xk = x(k)
yk = y(k)
call triar(xi,yi,xj,yj,xk,yk,tmp)
# Downweight the area by 1/3, since each
# triangle eventually appears 3 times over.
area = area+tmp/3.
}
delsum(i1,1) = xi
delsum(i1,2) = yi
delsum(i1,3) = npt
delsum(i1,4) = area
}
return
end
SHAR_EOF
if test -f 'delseg.r'
then
echo shar: over-writing existing file "'delseg.r'"
fi
cat << \SHAR_EOF > 'delseg.r'
subroutine delseg(delsgs,ndel,nadj,madj,x,y,npd,ntot,ind,nerror)
# Output the endpoints of the line segments joining the
# vertices of the Delaunay triangles.
# Called by master.
implicit double precision(a-h,o-z)
dimension nadj(-3:ntot,0:madj), x(-3:ntot), y(-3:ntot)
dimension delsgs(6,1), ind(1)
logical value
# For each distinct pair of points i and j, if they are adjacent
# then put their endpoints into the output array.
kseg = 0
do i1 = 2,npd {
i = ind(i1)
do j1 = 1,i1-1 {
j = ind(j1)
call adjchk(i,j,value,nadj,madj,ntot,nerror)
if(nerror>0) return
if(value) {
kseg = kseg+1
if(kseg > ndel) {
nerror = 14
return
}
delsgs(1,kseg) = x(i)
delsgs(2,kseg) = y(i)
delsgs(3,kseg) = x(j)
delsgs(4,kseg) = y(j)
delsgs(5,kseg) = i1
delsgs(6,kseg) = j1
}
}
}
ndel = kseg
return
end
SHAR_EOF
if test -f 'dirout.r'
then
echo shar: over-writing existing file "'dirout.r'"
fi
cat << \SHAR_EOF > 'dirout.r'
subroutine dirout(dirsum,nadj,madj,x,y,ntot,npd1,npd,rw,ind,eps,nerror)
# Output the description of the Dirichlet tile centred at point
# i for i = 1, ..., npd. Do this in the original order of the
# points, not in the order into which they have been bin-sorted.
# Called by master.
implicit double precision(a-h,o-z)
dimension nadj(-3:ntot,0:madj), x(-3:ntot), y(-3:ntot)
dimension dirsum(npd1,3), ind(npd1), rw(4)
logical collin, intfnd, bptab, bptcd
# Note that at this point some Delaunay neighbors may be
# `spurious'; they are the corners of a `large' rectangle in which
# the rectangular window of interest has been suspended. This
# large window was brought in simply to facilitate output concerning
# the Dirichlet tesselation. They were added to the triangulation
# in the routine `dirseg' which ***must*** therefore be called before
# this routine (`dirout') is called. (Likewise `dirseg' must be called
# ***after*** `delseg' and `delout' have been called.)
# Dig out the corners of the rectangular window.
xmin = rw(1)
xmax = rw(2)
ymin = rw(3)
ymax = rw(4)
do i1 = 1,npd {
area = 0. # Initialize the area of the ith tile to zero.
nbpt = 0 # Initialize the number of boundary points of
# the ith tile to zero.
npt = 0 # Initialize the number of tile boundaries to zero.
i = ind(i1)
np = nadj(i,0)
xi = x(i)
yi = y(i)
# Output the point number, its coordinates, and the number of
# its Delaunay neighbors == the number of boundary segments in
# its Dirichlet tile.
# For each Delaunay neighbor, find the circumcentres of the
# triangles on each side of the segment joining point i to that
# neighbor.
do j1 = 1,np {
j = nadj(i,j1)
xj = x(j)
yj = y(j)
xij = 0.5*(xi+xj)
yij = 0.5*(yi+yj)
call pred(k,i,j,nadj,madj,ntot,nerror)
if(nerror > 0) return
call succ(l,i,j,nadj,madj,ntot,nerror)
if(nerror > 0) return
call circen(i,k,j,a,b,x,y,ntot,eps,collin,nerror)
if(nerror>0) return
if(collin) {
nerror = 13
return
}
call circen(i,j,l,c,d,x,y,ntot,eps,collin,nerror)
if(nerror>0) return
if(collin) {
nerror = 13
return
}
# Increment the area of the current Dirichlet
# tile (intersected with the rectangular window) by applying
# Stokes' Theorem to the segment of tile boundary joining
# (a,b) to (c,d). (Note that the direction is anti-clockwise.)
call stoke(a,b,c,d,rw,tmp,sn,eps,nerror)
if(nerror > 0) return
area = area+sn*tmp
# If a circumcentre is outside the rectangular window, replace
# it with the intersection of the rectangle boundary with the
# line joining the circumcentre to the midpoint of
# (xi,yi)->(xj,yj). Then output the number of the current
# Delaunay neighbor and the two circumcentres (or the points
# with which they have been replaced).
call inside(a,b,xij,yij,ai,bi,rw,intfnd,bptab)
if(intfnd) {
call inside(c,d,xij,yij,ci,di,rw,intfnd,bptcd)
if(!intfnd) {
nerror = 17
return
}
if(bptab & bptcd) {
xm = 0.5*(ai+ci)
ym = 0.5*(bi+di)
if(xmin 'dirseg.r'
subroutine dirseg(dirsgs,ndir,nadj,madj,x,y,npd,ntot,rw,eps,ind,nerror)
# Output the endpoints of the segments of boundary of Dirichlet
# tiles. (Do it economically; each such segment once and only once.)
# Called by master.
implicit double precision(a-h,o-z)
dimension nadj(-3:ntot,0:madj), x(-3:ntot), y(-3:ntot)
dimension dirsgs(8,ndir), rw(4), ind(1)
logical collin, adjace, intfnd, bptab, bptcd, goferit
nerror = -1
# Add in some dummy corner points, outside the actual window.
# Far enough out so that no resulting tile boundaries intersect the
# window.
# Note that these dummy corners are needed by the routine `dirout'
# but will screw things up for `delseg' and `delout'. Therefore
# this routine (`dirseg') must be called ***before*** dirout, and
# ***after*** delseg and delout.
# Dig out the corners of the rectangular window.
xmin = rw(1)
xmax = rw(2)
ymin = rw(3)
ymax = rw(4)
a = xmax-xmin
b = ymax-ymin
c = sqrt(a*a+b*b)
nstt = npd+1
i = nstt
x(i) = xmin-c
y(i) = ymin-c
i = i+1
x(i) = xmax+c
y(i) = ymin-c
i = i+1
x(i) = xmax+c
y(i) = ymax+c
i = i+1
x(i) = xmin-c
y(i) = ymax+c
nstp = npd+4
do j = nstt,nstp {
call addpt(j,nadj,madj,x,y,ntot,eps,nerror)
if(nerror > 0) return
}
# Put the segments into the array dirsgs.
# For each distinct pair of (genuine) data points, find out if they are
# adjacent. If so, find the circumcentres of the triangles lying on each
# side of the segment joining them.
kseg = 0
do i1 = 2,npd {
i = ind(i1)
do j1 = 1,i1-1 {
j = ind(j1)
call adjchk(i,j,adjace,nadj,madj,ntot,nerror)
if(nerror > 0) return
if(adjace) {
xi = x(i)
yi = y(i)
xj = x(j)
yj = y(j)
# Let (xij,yij) be the midpoint of the segment joining
# (xi,yi) to (xj,yj).
xij = 0.5*(xi+xj)
yij = 0.5*(yi+yj)
call pred(k,i,j,nadj,madj,ntot,nerror)
if(nerror > 0) return
call circen(i,k,j,a,b,x,y,ntot,eps,collin,nerror)
if(nerror > 0) return
if(collin) {
nerror = 12
return
}
# If a circumcentre is outside the rectangular window
# of interest, draw a line joining it to the mid-point
# of (xi,yi)->(xj,yj). Find the intersection of this
# line with the boundary of the window; for (a,b),
# call the point of intersection (ai,bi). For (c,d),
# call it (ci,di).
call inside(a,b,xij,yij,ai,bi,rw,intfnd,bptab)
if(!intfnd) {
nerror = 16
return
}
call succ(l,i,j,nadj,madj,ntot,nerror)
if(nerror > 0) return
call circen(i,j,l,c,d,x,y,ntot,eps,collin,nerror)
if(nerror > 0) return
if(collin) {
nerror = 12
return
}
call inside(c,d,xij,yij,ci,di,rw,intfnd,bptcd)
if(!intfnd) {
nerror = 16
return
}
goferit = .false.
if(bptab & bptcd) {
xm = 0.5*(ai+ci)
ym = 0.5*(bi+di)
if(xmin ndir) {
nerror = 15
return
}
dirsgs(1,kseg) = ai
dirsgs(2,kseg) = bi
dirsgs(3,kseg) = ci
dirsgs(4,kseg) = di
dirsgs(5,kseg) = i1
dirsgs(6,kseg) = j1
if(bptab) dirsgs(7,kseg) = 1.d0
else dirsgs(7,kseg) = 0.d0
if(bptcd) dirsgs(8,kseg) = 1.d0
else dirsgs(8,kseg) = 0.d0
}
}
}
}
ndir = kseg
return
end
SHAR_EOF
if test -f 'driver.r'
then
echo shar: over-writing existing file "'driver.r'"
fi
cat << \SHAR_EOF > 'driver.r'
# Adjust the following ***eleven*** (11) ``defines'' as necessary:
# If error number 4 is encountered, MADJ should be increased.
define N 7
define NDX 2
define NDY 2
define MADJ 20 # max(20,ceiling(3*sqrt(NTOT))), NTOT = N+NDX*NDY+4
define SORT .true.
define EPS 1.d-9
define FRAC 1.d-4
define RW1 0.00
define RW2 10.0
define RW3 0.00
define RW4 10.0
#
# This ``define'' may need to be adjusted in response to error 14 or 15.
#
define NSGS MADJ*(MADJ+1)/2
#
# These ``defines'' should ***NOT*** be changed.
#
define NDM NDX*NDY
define NPD N+NDM
define NTOT NPD+4
#
implicit double precision(a-h,o-z)
logical sort
dimension x(-3:NTOT), y(-3:NTOT), rw(4)
dimension nadj(-3:NTOT,0:MADJ)
dimension ind(NPD), tx(NPD), ty(NPD), ilst(NPD)
dimension delsgs(6,NSGS), dirsgs(8,NSGS)
dimension delsum(NPD,4), dirsum(NPD,3)
read(5,*) (x(i), y(i), i=1,N)
rw(1) = RW1
rw(2) = RW2
rw(3) = RW3
rw(4) = RW4
ndx = NDX
ndy = NDX
n = N
npd1 = NPD
ntot = NTOT
madj = MADJ
ndel = NSGS
ndir = NSGS
eps = EPS
frac = FRAC
sort = SORT
do i = -3,ntot {
do j = 0,madj{
nadj(i,j) = 0
}
}
call master(x,y,sort,rw,n,npd1,npd,ntot,ndx,ndy,nadj,madj,frac,ind,tx,ty,ilst,
eps,delsgs,ndel,delsum,dirsgs,ndir,dirsum,nerror)
if(nerror > 0) {
print *, 'nerror =', nerror
stop
}
adel = 0.d0
adir = 0.d0
do i = 1, npd {
adel = adel + delsum(i,4)
}
do i = 1, npd {
adir = adir + dirsum(i,3)
}
ndum = npd - n
print *, ''
print *, 'Delaunay area: ', adel
print *, 'Dirichlet area: ', adir
print *, ''
print *, 'Delaunay segments:'
write(6,50) ((delsgs(i,j),i=1,4),(int(delsgs(i,j)),i=5,6),j=1,ndel)
print *, ''
print *, 'Dirichlet segments:'
write(6,51) ((dirsgs(i,j),i=1,4),(int(dirsgs(i,j)),i=5,8),j=1,ndir)
print *, ''
print *, 'Delaunay summary:'
write(6,52) ((delsum(i,j),j=1,4),i=1,npd)
print *, ''
print *, 'Dirichlet summary:'
write(6,53) ((dirsum(i,j),j=1,3),i=1,npd)
print *, ''
print *, 'First', n, ' points are data points.'
print *, 'Last', ndum, ' points are dummy points.'
print *, ''
50 format(4f10.4,2i7)
51 format(4f10.4,2i7,' ',2i3)
52 format(4f10.4)
53 format(3f10.4)
stop
end
SHAR_EOF
if test -f 'dumpts.r'
then
echo shar: over-writing existing file "'dumpts.r'"
fi
cat << \SHAR_EOF > 'dumpts.r'
subroutine dumpts(x,y,ntot,rw,ndx,ndy,n)
# Add in an ndx-by-ndy grid of dummy points.
# Called by master.
implicit double precision(a-h,o-z)
dimension x(-3:ntot), y(-3:ntot), rw(4)
if(ndx<=0|ndy<=0) return
# Dig out the corners of the rectangular window.
xmin = rw(1)
xmax = rw(2)
ymin = rw(3)
ymax = rw(4)
# If only one grid-line perp. to the x-axis, made it the midline;
# else make the lines equi-spaced from side to side.
if(ndx==1) x0 = 0.5*(xmin+xmax)
else {
dx = (xmax-xmin)/(ndx-1)
x0 = xmin
}
# Likewise for the grid-lines perp. to the y-axis.
if(ndy==1) y0 = 0.5*(ymin+ymax)
else {
dy = (ymax-ymin)/(ndy-1)
y0 = ymin
}
xt = x0
yt = y0
nt = n
# Run along the grid-lines vertically; there are ndx such
# vertical lines; add ndy points on each.
do i = 1,ndx {
do j = 1,ndy {
nt = nt+1
x(nt) = xt
y(nt) = yt
yt = yt+dy
}
xt = xt+dx
yt = y0
}
return
end
SHAR_EOF
if test -f 'eldup.r'
then
echo shar: over-writing existing file "'eldup.r'"
fi
cat << \SHAR_EOF > 'eldup.r'
subroutine eldup(x,y,n,npd,rw,frac,tx,ty)
implicit double precision(a-h,o-z)
logical newpt
dimension x(npd), y(npd), tx(npd), ty(npd), rw(4)
xtol = frac*(rw(2)-rw(1))
ytol = frac*(rw(4)-rw(3))
tx(1) = x(1)
ty(1) = y(1)
k = 1
newpt = .true.
do i = 2,npd {
newpt = .true.
do j = 1,i-1 {
dx = abs(x(i)-x(j))
dy = abs(y(i)-y(j))
if(dx < xtol & dy < ytol) {
newpt = .false.
break
}
}
if(newpt) {
k = k + 1
tx(k) = x(i)
ty(k) = y(i)
}
if(i==n) n1 = k
}
n = n1
npd = k
do i = 1,npd {
x(i) = tx(i)
y(i) = ty(i)
}
return
end
SHAR_EOF
if test -f 'initad.r'
then
echo shar: over-writing existing file "'initad.r'"
fi
cat << \SHAR_EOF > 'initad.r'
subroutine initad(j,nadj,madj,x,y,ntot,eps,nerror)
# Intial adding-in of a new point j.
# Called by addpt.
implicit double precision(a-h,o-z)
dimension nadj(-3:ntot,0:madj), x(-3:ntot), y(-3:ntot)
integer tau(3)
# Find the triangle containing vertex j.
call trifnd(j,tau,nedge,nadj,madj,x,y,ntot,eps,nerror)
if(nerror > 0) return
# If the new point is on the edge of a triangle, detach the two
# vertices of that edge from each other. Also join j to the vertex
# of the triangle on the reverse side of that edge from the `found'
# triangle (defined by tau) -- given that there ***is*** such a triangle.
if(nedge!=0) {
ip = nedge
i = ip-1
if(i==0) i = 3 # Arithmetic modulo 3.
call pred(k,tau(i),tau(ip),nadj,madj,ntot,nerror)
if(nerror > 0) return
call succ(kk,tau(ip),tau(i),nadj,madj,ntot,nerror)
if(nerror > 0) return
call delet(tau(i),tau(ip),nadj,madj,ntot,nerror)
if(nerror > 0) return
if(k==kk) call insrt(j,k,nadj,madj,x,y,ntot,nerror,eps)
if(nerror > 0) return
}
# Join the new point to each of the three vertices.
do i = 1,3 {
call insrt(j,tau(i),nadj,madj,x,y,ntot,nerror,eps)
if(nerror > 0) return
}
return
end
SHAR_EOF
if test -f 'inside.r'
then
echo shar: over-writing existing file "'inside.r'"
fi
cat << \SHAR_EOF > 'inside.r'
subroutine inside(a,b,c,d,ai,bi,rw,intfnd,bpt)
# Get a point ***inside*** the rectangular window on the ray from
# one circumcentre to the next one. I.e. if the `next one' is
# inside, then that's it; else find the intersection of this ray with
# the boundary of the rectangle.
# Called by dirseg, dirout.
implicit double precision(a-h,o-z)
dimension rw(4)
logical intfnd, bpt
# Note that (a,b) is the circumcenter of a Delaunay triangle,
# and that (c,d) is the midpoint of one of its sides.
# When `inside' is called by `dirout' it is possible for (c,d) to
# lie ***outside*** the rectangular window, and for the ray not to
# intersect the window at all. (The point (c,d) might be the midpoint
# of a Delaunay edge connected to a `fake outer corner', added to facilitate
# constructing a tiling that completely covers the actual window.)
# The variable `intfnd' acts as an indicator as to whether such an
# intersection has been found.
# The variable `bpt' acts as an indicator as to whether the returned
# point (ai,bi) is a true circumcentre, inside the window (bpt == .false.),
# or is the intersection of a ray with the boundary of the window
# (bpt = .true.).
intfnd = .true.
bpt = .true.
# Dig out the corners of the rectangular window.
xmin = rw(1)
xmax = rw(2)
ymin = rw(3)
ymax = rw(4)
# Check if (a,b) is inside the rectangle.
if(xmin<=a&a<=xmax&ymin<=b&b<=ymax) {
ai = a
bi = b
bpt = .false.
return
}
# Look for appropriate intersections with the four lines forming
# the sides of the rectangular window.
# Line 1: x = xmin.
if(axmax) {
ai = xmax
s = (d-b)/(c-a)
t = b-s*a
bi = s*ai+t
if(ymin<=bi&bi<=ymax) return
}
# Line 4: y = ymax.
if(b>ymax) {
bi = ymax
s = (c-a)/(d-b)
t = a-s*b
ai = s*bi+t
if(xmin<=ai&ai<=xmax) return
}
intfnd = .false.
return
end
SHAR_EOF
if test -f 'insrt.r'
then
echo shar: over-writing existing file "'insrt.r'"
fi
cat << \SHAR_EOF > 'insrt.r'
subroutine insrt(i,j,nadj,madj,x,y,ntot,nerror,eps)
# Insert i and j into each other's adjacency list.
# Called by master, initad, swap.
implicit double precision(a-h,o-z)
dimension nadj(-3:ntot,0:madj), x(-3:ntot), y(-3:ntot)
logical adj
# Check whether i and are in each other's adjacency lists.
call adjchk(i,j,adj,nadj,madj,ntot,nerror)
if(nerror > 0) return
if(adj) return
# If not, find where in each list they should respectively be.
call locn(i,j,kj,nadj,madj,x,y,ntot,eps)
call locn(j,i,ki,nadj,madj,x,y,ntot,eps)
# Put them in each other's lists in the appropriate position.
call insrt1(i,j,kj,nadj,madj,ntot,nerror)
if(nerror >0) return
call insrt1(j,i,ki,nadj,madj,ntot,nerror)
if(nerror >0) return
return
end
SHAR_EOF
if test -f 'insrt1.r'
then
echo shar: over-writing existing file "'insrt1.r'"
fi
cat << \SHAR_EOF > 'insrt1.r'
subroutine insrt1(i,j,kj,nadj,madj,ntot,nerror)
# Insert j into the adjacency list of i.
# Called by insrt.
implicit double precision(a-h,o-z)
dimension nadj(-3:ntot,0:madj)
nerror = -1
# Variable kj is the index which j ***will***
# have when it is inserted into the adjacency list of i in
# the appropriate position.
# If the adjacency list of i had no points just stick j into the list.
n = nadj(i,0)
if(n==0) {
nadj(i,0) = 1
nadj(i,1) = j
return
}
# If the adjacency list had some points, move everything ahead of the
# kj-th place one place forward, and put j in position kj.
kk = n+1
if(kk>madj) { # Watch out for over-writing!!!
nerror = 4
return
}
while(kk>kj) {
nadj(i,kk) = nadj(i,kk-1)
kk = kk-1
}
nadj(i,kj) = j
nadj(i,0) = n+1
return
end
SHAR_EOF
if test -f 'locn.r'
then
echo shar: over-writing existing file "'locn.r'"
fi
cat << \SHAR_EOF > 'locn.r'
subroutine locn(i,j,kj,nadj,madj,x,y,ntot,eps)
# Find the appropriate location for j in the adjacency list
# of i. This is the index which j ***will*** have when
# it is inserted into the adjacency list of i in the
# appropriate place.
# Called by insrt.
implicit double precision(a-h,o-z)
dimension nadj(-3:ntot,0:madj), x(-3:ntot), y(-3:ntot)
logical before
n = nadj(i,0)
# If there is nothing already adjacent to i, then j will have place 1.
if(n==0) {
kj = 1
return
}
# Run through i's list, checking if j should come before each element
# of that list. (I.e. if i, j, and k are in anti-clockwise order.)
# If j comes before the kj-th item, but not before the (kj-1)-st, then
# j should have place kj.
do ks = 1,n {
kj = ks
k = nadj(i,kj)
call acchk(i,j,k,before,x,y,ntot,eps)
if(before) {
km = kj-1
if(km==0) km = n
k = nadj(i,km)
call acchk(i,j,k,before,x,y,ntot,eps)
if(before) next
# If j is before 1 and after n, then it should
# have place n+1.
if(kj==1) kj = n+1
return
}
}
# We've gone right through the list and haven't been before
# the kj-th item ***and*** after the (kj-1)-st on any occasion.
# Therefore j is before everything (==> place 1) or after
# everything (==> place n+1).
if(before) kj = 1
else kj = n+1
return
end
SHAR_EOF
if test -f 'master.r'
then
echo shar: over-writing existing file "'master.r'"
fi
cat << \SHAR_EOF > 'master.r'
subroutine master(x,y,sort,rw,n,npd1,npd,ntot,ndx,ndy,nadj,madj,frac,ind,tx,ty,
ilst,eps,delsgs,ndel,delsum,dirsgs,ndir,dirsum,nerror)
# Master subroutine, to organize all the others.
# Called by driver.
implicit double precision(a-h,o-z)
logical sort
dimension x(-3:ntot), y(-3:ntot)
dimension nadj(-3:ntot,0:madj)
dimension ind(npd1), tx(npd1), ty(npd1), ilst(npd1), rw(4)
dimension delsgs(6,ndel), dirsgs(8,ndir)
dimension delsum(npd1,4), dirsum(npd1,3)
# Define one.
one = 1.d0
# Eliminate data points outside the rectangular window.
k = 0
do i = 1,n {
if(x(i) >= rw(1) & x(i) <= rw(2) & y(i) >= rw(3) & y(i) <= rw(4)) {
k = k + 1
tx(k) = x(i)
ty(k) = y(i)
}
}
n = k
do i = 1,n {
x(i) = tx(i)
y(i) = ty(i)
}
# Add in dummy points if required.
call dumpts(x,y,ntot,rw,ndx,ndy,n)
# Eliminate duplicated points.
npd = n + ndx*ndy
call eldup(x(1),y(1),n,npd,rw,frac,tx,ty)
# Sort the points into bins, the number of such being approx. sqrt(n).
if(sort) {
call binsrt(x,y,ntot,rw,npd,ind,tx,ty,ilst,nerror)
if(nerror > 0) return
}
else {
do i = 1,npd {
ind(i) = i
}
}
# Initialize the adjacency list to 0.
do i = -3,ntot {
do j = 0,madj {
nadj(i,j) = 0
}
}
# Put the four ideal points into x and y and the adjacency list.
# The ideal points are given pseudo-coordinates
# (-1,-1), (1,-1), (1,1), and (-1,1). They are numbered as
# 0 -1 -2 -3
# i.e. the numbers decrease anticlockwise from the
# `bottom left corner'.
x(-3) = -one
y(-3) = one
x(-2) = one
y(-2) = one
x(-1) = one
y(-1) = -one
x(0) = -one
y(0) = -one
do i = 1,4 {
j = i-4
k = j+1
if(k>0) k = -3
call insrt(j,k,nadj,madj,x,y,ntot,nerror,eps)
if(nerror>0) return
}
# Put in the first of the point set into the adjacency list.
do i = 1,4 {
j = i-4
call insrt(1,j,nadj,madj,x,y,ntot,nerror,eps)
if(nerror>0) return
}
# Now add the rest of the point set
do j = 2,npd {
call addpt(j,nadj,madj,x,y,ntot,eps,nerror)
if(nerror>0) return
}
# Obtain the description of the triangulation.
call delseg(delsgs,ndel,nadj,madj,x,y,npd,ntot,ind,nerror)
if(nerror>0) return
call delout(delsum,nadj,madj,x,y,ntot,npd1,npd,ind,nerror)
if(nerror>0) return
call dirseg(dirsgs,ndir,nadj,madj,x,y,npd,ntot,rw,eps,ind,nerror)
if(nerror>0) return
call dirout(dirsum,nadj,madj,x,y,ntot,npd1,npd,rw,ind,eps,nerror)
return
end
SHAR_EOF
if test -f 'pred.r'
then
echo shar: over-writing existing file "'pred.r'"
fi
cat << \SHAR_EOF > 'pred.r'
subroutine pred(kpr,i,j,nadj,madj,ntot,nerror)
# Find the predecessor of j in the adjacency list of i.
# Called by initad, trifnd, swap, dirseg, dirout.
implicit double precision(a-h,o-z)
dimension nadj(-3:ntot,0:madj)
nerror = -1
n = nadj(i,0)
# If the adjacency list of i is empty, then clearly j has no predecessor
# in this adjacency list. Something's wrong; stop.
if(n==0) {
nerror = 5
return
}
# The adjacency list of i is non-empty; search through it until j is found;
# subtract 1 from the location of j, and find the contents of this new location
do k = 1,n {
if(j==nadj(i,k)) {
km = k-1
if(km<1) km = n # Take km modulo n. (The adjacency list
kpr = nadj(i,km) # is circular.)
return
}
}
# The adjacency list doesn't contain j. Something's wrong; stop.
nerror = 6
return
end
SHAR_EOF
if test -f 'qtest.r'
then
echo shar: over-writing existing file "'qtest.r'"
fi
cat << \SHAR_EOF > 'qtest.r'
subroutine qtest(h,i,j,k,shdswp,x,y,ntot,eps,nerror)
# Test whether the LOP is satisified; i.e. whether vertex j is
# outside the circumcircle of vertices h, i, and k of the
# quadrilateral. Vertex h is the vertex being added; i and k are the
# vertices of the quadrilateral which are currently joined; j is the
# vertex which is ``opposite'' the vertex being added. If the LOP is
# not satisfied, the shdswp ("should-swap") is true, i.e. h and j
# should be joined, rather than i and k.
# Called by swap.
implicit double precision(a-h,o-z)
dimension x(-3:ntot), y(-3:ntot)
integer h
logical shdswp
nerror = -1
# Look for ideal points.
if(i<=0) ii = 1
else ii = 0
if(j<=0) jj = 1
else jj = 0
if(k<=0) kk = 1
else kk = 0
ijk = ii*4+jj*2+kk
switch(ijk) {
# All three corners other than h (the point currently being
# added) are ideal --- so h, i, and k are co-linear; so don't swap.
case 7:
shdswp = .false.
return
# If i and j are ideal, find out which of h and k is closer to the
# intersection point of the two diagonals, and choose the diagonal
# emanating from that vertex. (I.e. if h is closer, swap.)
# Unless swapping yields a non-convex quadrilateral!!!
case 6:
ja = -j
xh = x(h)
yh = y(h)
xk = x(k)
yk = y(k)
test =(xh*yk+xk*yh-xh*yh-xk*yk)*(-1)**ja
if(test>0) shdswp = .true.
else shdswp = .false.
# Check for convexity:
if(shdswp) call acchk(j,k,h,shdswp,x,y,ntot,eps)
return
# Vertices i and k are ideal --- can't happen, but if it did, we'd
# increase the minimum angle ``from 0 to more than 2*0'' by swapping,
# and the quadrilateral would definitely be convex, so let's say swap.
case 5:
shdswp = .true.
return
# If i is ideal we'd increase the minimum angle ``from 0 to more than
# 2*0'' by swapping, so just check for convexity:
case 4:
call acchk(j,k,h,shdswp,x,y,ntot,eps)
return
# If j and k are ideal, this is like unto case 6.
case 3:
ja = -j
xi = x(i)
yi = y(i)
xh = x(h)
yh = y(h)
test = (xh*yi+xi*yh-xh*yh-xi*yi)*(-1)**ja
if(test>0.) shdswp = .true.
else shdswp = .false.
# Check for convexity:
if(shdswp) call acchk(j,i,h,shdswp,x,y,ntot,eps)
return
# If j is ideal we'd decrease the minimum angle ``from more than 2*0
# to 0'', by swapping; so don't swap.
case 2:
shdswp = .false.
return
# If k is ideal, this is like unto case 4.
case 1:
call acchk(i,j,h,shdswp,x,y,ntot,eps) # This checks
# for convexity.
return
# If none of the `other' three corners are ideal, do the Lee-Schacter
# test for the LOP.
case 0:
call qtest1(h,i,j,k,x,y,ntot,eps,shdswp,nerror)
return
default: # This CAN'T happen.
nerror = 7
return
}
end
SHAR_EOF
if test -f 'qtest1.r'
then
echo shar: over-writing existing file "'qtest1.r'"
fi
cat << \SHAR_EOF > 'qtest1.r'
subroutine qtest1(h,i,j,k,x,y,ntot,eps,shdswp,nerror)
# The Lee-Schacter test for the LOP (all points are
# real, i.e. non-ideal). If the LOP is ***not***
# satisfied then the diagonals should be swapped,
# i.e. shdswp ("should-swap") is true.
# Called by qtest.
implicit double precision(a-h,o-z)
dimension x(-3:ntot), y(-3:ntot)
integer h
logical shdswp
# The vertices of the quadrilateral are labelled
# h, i, j, k in the anticlockwise direction, h
# being the point of central interest.
# Make sure the quadrilateral is convex, so that
# it makes sense to swap the diagonal.
call acchk(i,j,k,shdswp,x,y,ntot,eps)
if(!shdswp) return
# Get the coordinates of vertices h and j.
xh = x(h)
yh = y(h)
xj = x(j)
yj = y(j)
# Find the centre of the circumcircle of vertices h, i, k.
call circen(h,i,k,x0,y0,x,y,ntot,eps,shdswp,nerror)
if(nerror>0) return
if(shdswp) return # The points h, i, and k are colinear, so
# the circumcircle has `infinite radius', so
# (xj,yj) is definitely inside.
# Check whether (xj,yj) is inside the circle of centre
# (x0,y0) and radius r = dist[(x0,y0),(xh,yh)]
a = x0-xh
b = y0-yh
r2 = a*a+b*b
a = x0-xj
b = y0-yj
ch = a*a+b*b
if(ch 'stoke.r'
subroutine stoke(x1,y1,x2,y2,rw,area,s1,eps,nerror)
# Apply Stokes' theorem to find the area of a polygon;
# we are looking at the boundary segment from (x1,y1)
# to (x2,y2), travelling anti-clockwise. We find the
# area between this segment and the horizontal base-line
# y = ymin, and attach a sign s1. (Positive if the
# segment is right-to-left, negative if left to right.)
# The area of the polygon is found by summing the result
# over all boundary segments.
# Just in case you thought this wasn't complicated enough,
# what we really want is the area of the intersection of
# the polygon with the rectangular window that we're using.
# Called by dirout.
implicit double precision(a-h,o-z)
dimension rw(4)
logical value
nerror = -1
# If the segment is vertical, the area is zero.
call testeq(x1,x2,eps,value)
if(value) {
area = 0.
s1 = 0.
return
}
# Find which is the right-hand end, and which is the left.
if(x1=xmax) {
area = 0.
return
}
# We're now looking at a trapezoidal region which may or may
# not protrude above or below the horizontal strip bounded by
# y = ymax and y = ymin.
ybot = min(yl,yr)
ytop = max(yl,yr)
# Case 1; ymax <= ybot:
# The `roof' of the trapezoid is entirely above the
# horizontal strip.
if(ymax<=ybot) {
area = (xr-xl)*(ymax-ymin)
return
}
# Case 2; ymin <= ybot <= ymax <= ytop:
# The `roof' of the trapezoid intersects the top of the
# horizontal strip (y = ymax) but not the bottom (y = ymin).
if(ymin<=ybot&ymax<=ytop) {
call testeq(slope,0.,eps,value)
if(value) {
w1 = 0.
w2 = xr-xl
}
else {
xit = xl+(ymax-yl)/slope
w1 = xit-xl
w2 = xr-xit
if(slope<0.) {
tmp = w1
w1 = w2
w2 = tmp
}
}
area = 0.5*w1*((ybot-ymin)+(ymax-ymin))+w2*(ymax-ymin)
return
}
# Case 3; ybot <= ymin <= ymax <= ytop:
# The `roof' intersects both the top (y = ymax) and
# the bottom (y = ymin) of the horizontal strip.
if(ybot<=ymin&ymax<=ytop) {
xit = xl+(ymax-yl)/slope
xib = xl+(ymin-yl)/slope
if(slope>0.) {
w1 = xit-xib
w2 = xr-xit
}
else {
w1 = xib-xit
w2 = xit-xl
}
area = 0.5*w1*(ymax-ymin)+w2*(ymax-ymin)
return
}
# Case 4; ymin <= ybot <= ytop <= ymax:
# The `roof' is ***between*** the bottom (y = ymin) and
# the top (y = ymax) of the horizontal strip.
if(ymin<=ybot&ytop<=ymax) {
area = 0.5*(xr-xl)*((ytop-ymin)+(ybot-ymin))
return
}
# Case 5; ybot <= ymin <= ytop <= ymax:
# The `roof' intersects the bottom (y = ymin) but not
# the top (y = ymax) of the horizontal strip.
if(ybot<=ymin&ymin<=ytop) {
call testeq(slope,0.,eps,value)
if(value) {
area = 0.
return
}
xib = xl+(ymin-yl)/slope
if(slope>0.) w = xr-xib
else w = xib-xl
area = 0.5*w*(ytop-ymin)
return
}
# Case 6; ytop <= ymin:
# The `roof' is entirely below the bottom (y = ymin), so
# there is no area contribution at all.
if(ytop<=ymin) {
area = 0.
return
}
# Default; all stuffed up:
nerror = 8
return
end
SHAR_EOF
if test -f 'succ.r'
then
echo shar: over-writing existing file "'succ.r'"
fi
cat << \SHAR_EOF > 'succ.r'
subroutine succ(ksc,i,j,nadj,madj,ntot,nerror)
# Find the successor of j in the adjacency list of i.
# Called by addpt, initad, trifnd, swap, delout, dirseg, dirout.
implicit double precision(a-h,o-z)
dimension nadj(-3:ntot,0:madj)
nerror = -1
n = nadj(i,0)
# If the adjacency list of i is empty, then clearly j has no successor
# in this adjacency list. Something's wrong; stop.
if(n==0) {
nerror = 9
return
}
# The adjacency list of i is non-empty; search through it until j is found;
# add 1 to the location of j, and find the contents of this new location.
do k = 1,n {
if(j==nadj(i,k)) {
kp = k+1
if(kp>n) kp = 1 # Take kp modulo n. (The adjacency list
ksc = nadj(i,kp) # is circular.)
return
}
}
# The adjacency list doesn't contain j. Something's wrong.
nerror = 10
return
end
SHAR_EOF
if test -f 'swap.r'
then
echo shar: over-writing existing file "'swap.r'"
fi
cat << \SHAR_EOF > 'swap.r'
subroutine swap(j,k1,k2,shdswp,nadj,madj,x,y,ntot,eps,nerror)
# The segment k1->k2 is a diagonal of a quadrilateral
# with a vertex at j (the point being added to the
# triangulation). If the LOP is not satisfied, swap
# it for the other diagonal.
# Called by addpt.
implicit double precision(a-h,o-z)
dimension nadj(-3:ntot,0:madj), x(-3:ntot), y(-3:ntot)
logical shdswp
# If vertices k1 and k2 are not connected there is no diagonal to swap.
# This could happen if vertices j, k1, and k2 were colinear, but shouldn't.
call adjchk(k1,k2,shdswp,nadj,madj,ntot,nerror)
if(nerror > 0) return
if(!shdswp) return
# Get the other vertex of the quadrilateral.
call pred(k,k1,k2,nadj,madj,ntot,nerror) # If these aren't the same, then
if(nerror > 0) return
call succ(kk,k2,k1,nadj,madj,ntot,nerror) # there is no other vertex.
if(nerror > 0) return
if(kk!=k) {
shdswp = .false.
return
}
# Check whether the LOP is satisified; i.e. whether
# vertex k is outside the circumcircle of vertices j, k1, and k2
call qtest(j,k1,k,k2,shdswp,x,y,ntot,eps,nerror)
if(nerror > 0) return
# Do the actual swapping.
if(shdswp) {
call delet(k1,k2,nadj,madj,ntot,nerror)
if(nerror > 0) return
call insrt(j,k,nadj,madj,x,y,ntot,nerror,eps)
if(nerror > 0) return
}
return
end
SHAR_EOF
if test -f 'testeq.r'
then
echo shar: over-writing existing file "'testeq.r'"
fi
cat << \SHAR_EOF > 'testeq.r'
subroutine testeq(a,b,eps,value)
# Test for the equality of a and b in a fairly
# robust way.
# Called by trifnd, circen, stoke.
implicit double precision(a-h,o-z)
logical value
# If b is essentially 0, check whether a is essentially zero also.
# The following is very sloppy! Must fix it!
if(abs(b)<=eps) {
if(abs(a)<=eps) value = .true.
else value = .false.
return
}
# Test if a is a `lot different' from b. (If it is
# they're obviously not equal.) This avoids under/overflow
# problems in dividing a by b.
if(abs(a)>10.*abs(b)|abs(a)<0.1*abs(b)) {
value = .false.
return
}
# They're non-zero and fairly close; compare their ratio with 1.
c = a/b
if(abs(c-1.)<=eps) value = .true.
else value = .false.
return
end
SHAR_EOF
if test -f 'triar.r'
then
echo shar: over-writing existing file "'triar.r'"
fi
cat << \SHAR_EOF > 'triar.r'
subroutine triar(x0,y0,x1,y1,x2,y2,area)
# Calculate the area of a triangle with given
# vertices, in anti clockwise direction.
# Called by delout.
implicit double precision(a-h,o-z)
area = 0.5*((x1-x0)*(y2-y0)-(x2-x0)*(y1-y0))
return
end
SHAR_EOF
if test -f 'trifnd.r'
then
echo shar: over-writing existing file "'trifnd.r'"
fi
cat << \SHAR_EOF > 'trifnd.r'
subroutine trifnd(j,tau,nedge,nadj,madj,x,y,ntot,eps,nerror)
# Find the triangle of the extant triangulation in which
# lies the point currently being added.
# Called by initad.
implicit double precision(a-h,o-z)
dimension nadj(-3:ntot,0:madj), x(-3:ntot), y(-3:ntot), xt(3), yt(3)
integer tau(3)
logical adjace
nerror = -1
# The first point must be added to the triangulation before
# calling trifnd.
if(j==1) {
nerror = 11
return
}
# Get the previous triangle:
j1 = j-1
tau(1) = j1
tau(3) = nadj(j1,1)
call pred(tau(2),j1,tau(3),nadj,madj,ntot,nerror)
if(nerror > 0) return
call adjchk(tau(2),tau(3),adjace,nadj,madj,ntot,nerror)
if(nerror>0) return
if(!adjace) {
tau(3) = tau(2)
call pred(tau(2),j1,tau(3),nadj,madj,ntot,nerror)
if(nerror > 0) return
}
# Move to the adjacent triangle in the direction of the new
# point, until the new point lies in this triangle.
1 continue
ntau = 0 # This number will identify the triangle to be moved to.
nedge = 0 # If the point lies on an edge, this number will identify that edge.
do i = 1,3 {
ip = i+1
if(ip==4) ip = 1 # Take addition modulo 3.
# Get the coordinates of the vertices of the current side,
# and of the point j which is being added:
xt(1) = x(tau(i))
yt(1) = y(tau(i))
xt(2) = x(tau(ip))
yt(2) = y(tau(ip))
xt(3) = x(j)
yt(3) = y(j)
# Create indicator telling which of tau(i), tau(ip), and j
# are ideal points. (The point being added, j, is ***never*** ideal.)
if(tau(i)<=0) i1 = 1
else i1 = 0
if(tau(ip)<=0) j1 = 1
else j1 = 0
k1 = 0
ijk = i1*4+j1*2+k1
# Calculate the ``normalized'' cross product; if this is positive
# then the point being added is to the left (as we move along the
# edge in an anti-clockwise direction). If the test value is positive
# for all three edges, then the point is inside the triangle. Note
# that if the test value is very close to zero, we might get negative
# values for it on both sides of an edge, and hence go into an
# infinite loop.
call cross(xt,yt,ijk,cprd)
if(cprd >= eps) continue
else if(cprd > -eps) nedge = ip
else {
ntau = ip
break
}
}
# We've played ring-around-the-triangle; now figure out the
# next move:
switch(ntau) {
case 0: # All tests >= 0.; the point is inside; return.
return
# The point is not inside; work out the vertices of the triangle to which
# to move. Notation: Number the vertices of the current triangle from 1 to 3,
# anti-clockwise. Then "triangle i+1" is adjacent to the side from vertex i to
# vertex i+1, where i+1 is taken modulo 3 (i.e. "3+1 = 1").
case 1:
# move to "triangle 1"
#tau(1) = tau(1)
tau(2) = tau(3)
call succ(tau(3),tau(1),tau(2),nadj,madj,ntot,nerror)
if(nerror > 0) return
case 2:
# move to "triangle 2"
#tau(1) = tau(1)
tau(3) = tau(2)
call pred(tau(2),tau(1),tau(3),nadj,madj,ntot,nerror)
if(nerror > 0) return
case 3:
# move to "triangle 3"
tau(1) = tau(3)
#tau(2) = tau(2)
call succ(tau(3),tau(1),tau(2),nadj,madj,ntot,nerror)
if(nerror > 0) return
}
# We've moved to a new triangle; check if the point being added lies
# inside this one.
go to 1
end
SHAR_EOF
if test -f 'err.list'
then
echo shar: over-writing existing file "'err.list'"
fi
cat << \SHAR_EOF > 'err.list'
Error list:
===========
nerror = 1: Contradictory adjacency lists. Error in adjchk.
nerror = 2: Number of points jumbled. Error in binsrt.
nerror = 3: Vertices of 'triangle' are collinear
and vertex 2 is not between 1 and 3. Error in circen.
nerror = 4: Number of adjacencies too large. Error in insrt.
(Increase MADJ in driver.r)
nerror = 5: Adjacency list of i is empty, and so cannot contain j.
Error in pred.
nerror = 6: Adjacency list of i does not contain j. Error in pred.
nerror = 7: Indicator ijk is out of range. (This CAN'T happen!)
Error in qtest.
nerror = 8: Fell through all six cases.
Something must be totally stuffed up. Error in stoke.
nerror = 9: Adjacency list of i is empty, and so cannot contain j.
Error in succ.
nerror = 10: Adjacency list of i does not contain j. Error in succ.
nerror = 11: No triangles to find. Error in trifnd.
nerror = 12: Vertices of triangle are collinear. Error in dirseg.
nerror = 13: Vertices of triangle are collinear. Error in dirout.
nerror = 14: Number of Delaunay segments exceeds alloted space.
Error in delseg. (Increase NSGS in driver.r)
nerror = 15: Number of Dirichlet segments exceeds alloted space.
Error in dirseg. (Increase NSGS in driver.r)
nerror = 16: Line from midpoint to circumcenter does not intersect
rectangle boundary; but it HAS to!!! Error in dirseg.
nerror = 17: Line from midpoint to circumcenter does not intersect
rectangle boundary; but it HAS to!!! Error in dirout.
SHAR_EOF
if test -f 'ex.dat'
then
echo shar: over-writing existing file "'ex.dat'"
fi
cat << \SHAR_EOF > 'ex.dat'
2.3 2.3
3 3
7 2
4 11
1 5
3 8
8 9
SHAR_EOF
if test -f 'ex.out'
then
echo shar: over-writing existing file "'ex.out'"
fi
cat << \SHAR_EOF > 'ex.out'
Delaunay area: 100.000000000000
Dirichlet area: 100.000000000000
Delaunay segments:
3.0000 3.0000 2.3000 2.3000 2 1
7.0000 2.0000 2.3000 2.3000 3 1
7.0000 2.0000 3.0000 3.0000 3 2
1.0000 5.0000 2.3000 2.3000 4 1
1.0000 5.0000 3.0000 3.0000 4 2
3.0000 8.0000 3.0000 3.0000 5 2
3.0000 8.0000 7.0000 2.0000 5 3
3.0000 8.0000 1.0000 5.0000 5 4
8.0000 9.0000 7.0000 2.0000 6 3
8.0000 9.0000 3.0000 8.0000 6 5
0.0000 0.0000 2.3000 2.3000 7 1
0.0000 0.0000 7.0000 2.0000 7 3
0.0000 0.0000 1.0000 5.0000 7 4
0.0000 10.0000 1.0000 5.0000 8 4
0.0000 10.0000 3.0000 8.0000 8 5
0.0000 10.0000 8.0000 9.0000 8 6
0.0000 10.0000 0.0000 0.0000 8 7
10.0000 0.0000 7.0000 2.0000 9 3
10.0000 0.0000 8.0000 9.0000 9 6
10.0000 0.0000 0.0000 0.0000 9 7
10.0000 10.0000 8.0000 9.0000 10 6
10.0000 10.0000 0.0000 10.0000 10 8
10.0000 10.0000 10.0000 0.0000 10 9
Dirichlet segments:
1.6500 3.6500 4.5600 0.7400 2 1 0 0
4.5600 0.7400 4.5128 0.0000 3 1 0 1
5.7500 5.5000 4.5600 0.7400 3 2 0 0
0.0000 2.8556 1.6500 3.6500 4 1 1 0
1.6500 3.6500 3.5000 5.5000 4 2 0 0
3.5000 5.5000 5.7500 5.5000 5 2 0 0
5.7500 5.5000 6.0588 5.7059 5 3 0 0
0.5000 7.5000 3.5000 5.5000 5 4 0 0
6.0588 5.7059 10.0000 5.1429 6 3 0 1
5.2000 10.0000 6.0588 5.7059 6 5 1 0
2.3000 0.0000 0.0000 2.3000 7 1 1 1
0.0000 7.4000 0.5000 7.5000 8 4 1 0
0.5000 7.5000 2.1667 10.0000 8 5 0 1
10.0000 3.2500 7.8333 0.0000 9 3 1 1
8.7500 10.0000 10.0000 7.5000 10 6 1 1
Delaunay summary:
2.3000 2.3000 4.0000 4.5000
3.0000 3.0000 4.0000 6.0500
7.0000 2.0000 6.0000 18.6667
1.0000 5.0000 5.0000 7.5000
3.0000 8.0000 5.0000 15.0000
8.0000 9.0000 5.0000 16.6667
0.0000 0.0000 4.0000 8.4500
0.0000 10.0000 4.0000 7.6667
10.0000 0.0000 3.0000 10.5000
10.0000 10.0000 2.0000 5.0000
Dirichlet summary:
4.0000 4.0000 9.0921
4.0000 0.0000 10.7385
5.0000 4.0000 23.3182
4.0000 2.0000 9.3942
5.0000 2.0000 18.0556
3.0000 4.0000 18.3148
1.0000 2.0000 2.6450
2.0000 2.0000 3.3583
1.0000 2.0000 3.5208
1.0000 2.0000 1.5625
First 6 points are data points.
Last 4 points are dummy points.
SHAR_EOF
# End of shell archive
exit 0