subroutine trifnd(j,tau,nedge,nadj,madj,x,y,ntot,eps,ntri,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)
dimension itmp(1)
integer tau(3)
logical adjace, anticl
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.
ktri = 0
1 continue
# Check that the vertices of the triangle listed in tau are
# in anticlockwise order. (If they aren't then reverse the order;
# if they are *still* not in anticlockwise order, theh alles
# upgefucken ist; throw an error.)
call acchk(tau(1),tau(2),tau(3),anticl,x,y,ntot,eps)
if(!anticl) {
call acchk(tau(3),tau(2),tau(1),anticl,x,y,ntot,eps)
if(!anticl) {
itmp(1) = j
call intpr("Point number =",-1,itmp,1)
call intpr("Previous triangle:",-1,tau,3)
call rexit("Both vertex orderings are clockwise. See help for deldir.")
} else {
ivtmp = tau(3)
tau(3) = tau(1)
tau(1) = ivtmp
}
}
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:
# case 0: All tests >= 0.; the point is inside; return.
if(ntau==0) 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"
if(ntau==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"
if(ntau==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"
if(ntau==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.
ktri = ktri + 1
if(ktri > ntri) {
itmp(1) = j
call intpr("Point being added:",-1,itmp,1)
call rexit("Cannot find an enclosing triangle. See help for deldir.")
}
go to 1
end
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.