subroutine aptetrn (px, py, pz, ax, ay, az, bx, by, bz,
& cx, cy, cz, dx, dy, dz, np, tol, kloc,
& pabc, pbcd, pcda, pdab, pmin,
& wa, wb, wc, wd, nloc, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTETRN
c
c call aptetrn (px, py, pz, ax, ay, az, bx, by, bz,
c & cx, cy, cz, dx, dy, dz, np, tol, kloc,
c & pabc, pbcd, pcda, pdab, pmin,
c & wa, wb, wc, wd, nloc, nerr)
c
c Version: aptetrn Updated 1990 November 27 14:00.
c aptetrn Originated 1990 May 16 10:00.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To find, for each of np points p = (px, py, pz), the distances
c pabc, pbcd, pcda, pdab of the point "p" from the triangular
c faces "abc", "bcd", "cda" and "dab" of the tetrahedron with
c vertices a = (ax, ay, az), b = (bx, by, bz), c = (cx, cy, cz)
c and d = (dx, dy, dz) in 3-D space; the minimum of those
c distances, pmin; and the flag nloc indicating whether or not
c the point "p" is inside the tetrahedron. Optionally (kloc = 1),
c to return the local tetrahedral coordinates (vertex weights) of
c point "p", wa, wb, wc and wd:
c
c p = wa * a + wb * b + wc * c + wd * d
c
c Flag nerr indicates any input error.
c
c Input: px, py, pz, ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz,
c np, tol.
c
c Output: pabc, pbcd, pcda, pdab, pmin, nloc, nerr.
c
c Calls: aptfdav, aptript
c
c
c Glossary:
c
c ax,ay,az Input The x, y, z coordinates of tetrahedron vertex "a".
c Size np.
c
c bx,by,bz Input The x, y, z coordinates of tetrahedron vertex "b".
c Size np.
c
c cx,cy,cz Input The x, y, z coordinates of tetrahedron vertex "c".
c Size np.
c
c dx,dy,dz Input The x, y, z coordinates of tetrahedron vertex "d".
c Size np.
c
c kloc Input Indicates, if 1, that the local tetrahedral coordinates
c wa, wb, wc and wd are to be returned.
c
c nerr Output Indicates an input error, if not 0.
c 1 if np is not positive.
c 2 if kloc is not 0 or 1.
c
c nloc Output Indicates the location of point "p" relative to the
c tetrahedron "abcd":
c -1 if the tetrahedron has zero volume.
c 0 if point "p" is outside the tetrahedron "abcd".
c (one or more of pabc, pbcd, pcda or pdab are
c negative).
c 1 if point "p" is inside the tetrahedron "abcd"
c (all four distances pabc, pbcd, pcda and
c pdab are non-negative).
c Size np.
c
c np Input Size of arrays px, py, pz, ax, ay, az, bx, by, bz,
c cx, cy, cz, dx, dy, dz, pabc, pbcd, pcda, pdab,
c pmin, and nloc.
c
c pabc Output Distance of point "p" from triangular face "abc" of
c tetrahedron "abcd". Size np. See wd.
c
c pbcd Output Distance of point "p" from triangular face "bcd" of
c tetrahedron "abcd". Size np. See wa.
c
c pcda Output Distance of point "p" from triangular face "cda" of
c tetrahedron "abcd". Size np. See wb.
c
c pdab Output Distance of point "p" from triangular face "dab" of
c tetrahedron "abcd". Size np. See wc.
c
c pmin Output Minimum of distances pabc, pbcd, pcda, pdab.
c Size np.
c
c px,py,pz Input The x, y, z coordinates of point "p". Size np.
c
c tol Input Numerical tolerance limit.
c On Cray computers, recommend 1.e-5 to 1.e-11.
c Used to adjust values of distances for which local
c coordinates are within tol of 0.0 or 1.0.
c
c wa Output Fractional distance of point "p" from triangular face
c "bcd" to vertex "a". Size np, if kloc = 1.
c Values within tol of 0.0 or 1.0 will be adjusted
c to tol or 1.0 - tol, respectively.
c
c wb Output Fractional distance of point "p" from triangular face
c "cda" to vertex "b". Size np, if kloc = 1.
c Values within tol of 0.0 or 1.0 will be adjusted
c to tol or 1.0 - tol, respectively.
c
c wc Output Fractional distance of point "p" from triangular face
c "dab" to vertex "c". Size np, if kloc = 1.
c Values within tol of 0.0 or 1.0 will be adjusted
c to tol or 1.0 - tol, respectively.
c
c wd Output Fractional distance of point "p" from triangular face
c "abc" to vertex "d". Size np, if kloc = 1.
c Values within tol of 0.0 or 1.0 will be adjusted
c to tol or 1.0 - tol, respectively.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments.
c---- Coordinate x of point "a".
dimension ax (1)
c---- Coordinate y of point "a".
dimension ay (1)
c---- Coordinate z of point "a".
dimension az (1)
c---- Coordinate x of point "b".
dimension bx (1)
c---- Coordinate y of point "b".
dimension by (1)
c---- Coordinate z of point "b".
dimension bz (1)
c---- Coordinate x of point "c".
dimension cx (1)
c---- Coordinate y of point "c".
dimension cy (1)
c---- Coordinate z of point "c".
dimension cz (1)
c---- Coordinate x of point "d".
dimension dx (1)
c---- Coordinate y of point "d".
dimension dy (1)
c---- Coordinate z of point "d".
dimension dz (1)
c---- 1 if point "p" is inside tetrahedron.
dimension nloc (1)
c---- Distance from point "p" to face "abc".
dimension pabc (1)
c---- Distance from point "p" to face "bcd".
dimension pbcd (1)
c---- Distance from point "p" to face "cda".
dimension pcda (1)
c---- Distance from point "p" to face "dab".
dimension pdab (1)
c---- Min. dist. from point "p" to a face.
dimension pmin (1)
c---- Coordinate x of point "p".
dimension px (1)
c---- Coordinate y of point "p".
dimension py (1)
c---- Coordinate z of point "p".
dimension pz (1)
c---- Fract. dist. of "p" from "bcd" to "a".
dimension wa (1)
c---- Fract. dist. of "p" from "cda" to "b".
dimension wb (1)
c---- Fract. dist. of "p" from "dab" to "c".
dimension wc (1)
c---- Fract. dist. of "p" from "abc" to "d".
dimension wd (1)
c.... Local variables.
c---- Distance from "bcd" to vertex "a".
common /laptetrn/ dista (64)
c---- Distance from "cda" to vertex "b".
common /laptetrn/ distb (64)
c---- Distance from "dab" to vertex "c".
common /laptetrn/ distc (64)
c---- Distance from "abc" to vertex "d".
common /laptetrn/ distd (64)
c---- A very small number.
common /laptetrn/ fuz
c---- Index in local array.
common /laptetrn/ n
c---- First index of subset of data.
common /laptetrn/ n1
c---- Last index of subset of data.
common /laptetrn/ n2
c---- 1 if local coordinated adjusted.
common /laptetrn/ nlim (64)
c---- Index in external array.
common /laptetrn/ nn
c---- Size of current subset of data.
common /laptetrn/ ns
c---- Fract. dist. of "p" from "bcd" to "a".
common /laptetrn/ wal (64)
c---- Fract. dist. of "p" from "cda" to "b".
common /laptetrn/ wbl (64)
c---- Fract. dist. of "p" from "dab" to "c".
common /laptetrn/ wcl (64)
c---- Fract. dist. of "p" from "abc" to "d".
common /laptetrn/ wdl (64)
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptetrn finding distances in a tetrahedron.',
cbug & ' np=',i3 /
cbug & (i3,' px,py,pz=',1p3e22.14 /
cbug & ' ax,ay,az=',1p3e22.14 /
cbug & ' bx,by,bz=',1p3e22.14 /
cbug & ' cx,cy,cz=',1p3e22.14 /
cbug & ' dx,dy,dz=',1p3e22.14))
cbug write (3, 9901) np, (n, px(n), py(n), pz(n), ax(n), ay(n), az(n),
cbug & bx(n), by(n), bz(n), cx(n), cy(n), cz(n), dx(n), dy(n), dz(n),
cbug & n = 1, np)
cbugc***DEBUG ends.
c.... initialize.
c---- A very small number.
fuz = 1.e-99
nerr = 0
c.... Test for input errors.
if (np .le. 0) then
nerr = 1
go to 210
endif
if ((kloc .lt. 0) .or. (kloc .gt. 1)) then
nerr = 2
go to 210
endif
c.... Set up the indices of the first subset of data.
n1 = 1
n2 = min (np, 64)
c.... Loop over the data subsets.
110 ns = n2 - n1 + 1
c.... Find the distances of the tetrahedron vertices from the opposite faces.
call aptript (ax(n1), ay(n1), az(n1), bx(n1), by(n1), bz(n1),
& cx(n1), cy(n1), cz(n1), dx(n1), dy(n1), dz(n1),
& ns, tol, dista, nerr)
call aptript (bx(n1), by(n1), bz(n1), cx(n1), cy(n1), cz(n1),
& dx(n1), dy(n1), dz(n1), ax(n1), ay(n1), az(n1),
& ns, tol, distb, nerr)
call aptript (cx(n1), cy(n1), cz(n1), dx(n1), dy(n1), dz(n1),
& ax(n1), ay(n1), az(n1), bx(n1), by(n1), bz(n1),
& ns, tol, distc, nerr)
call aptript (dx(n1), dy(n1), dz(n1), ax(n1), ay(n1), az(n1),
& bx(n1), by(n1), bz(n1), cx(n1), cy(n1), cz(n1),
& ns, tol, distd, nerr)
c.... Find the distances of point "p" from the tetrahedron triangular faces.
call aptript (px(n1), py(n1), pz(n1), bx(n1), by(n1), bz(n1),
& cx(n1), cy(n1), cz(n1), dx(n1), dy(n1), dz(n1),
& ns, tol, pbcd(n1), nerr)
call aptript (px(n1), py(n1), pz(n1), cx(n1), cy(n1), cz(n1),
& dx(n1), dy(n1), dz(n1), ax(n1), ay(n1), az(n1),
& ns, tol, pcda(n1), nerr)
call aptript (px(n1), py(n1), pz(n1), dx(n1), dy(n1), dz(n1),
& ax(n1), ay(n1), az(n1), bx(n1), by(n1), bz(n1),
& ns, tol, pdab(n1), nerr)
call aptript (px(n1), py(n1), pz(n1), ax(n1), ay(n1), az(n1),
& bx(n1), by(n1), bz(n1), cx(n1), cy(n1), cz(n1),
& ns, tol, pabc(n1), nerr)
c.... Find the local tetrahedral coordinates (vertex weights) of point "p".
c---- Loop over subset of data.
do 120 n = 1, ns
nn = n + n1 - 1
wal(n) = pbcd(nn) / (dista(n) + fuz)
wbl(n) = pcda(nn) / (distb(n) + fuz)
wcl(n) = pdab(nn) / (distc(n) + fuz)
wdl(n) = pabc(nn) / (distd(n) + fuz)
c---- End of loop over subset of data.
120 continue
c.... See if the local tetrahedral coordinates should be adjusted.
c---- Adjust the local coordinates.
if (tol .gt. 0.0) then
call aptfdav (wal, ns, 1, tol, nlim, nerr)
call aptfdav (wbl, ns, 1, tol, nlim, nerr)
call aptfdav (wcl, ns, 1, tol, nlim, nerr)
call aptfdav (wdl, ns, 1, tol, nlim, nerr)
c---- Tested tol.
endif
c.... Find the (signed) distances of point "p" from the triangular faces,
c.... the minimum distance, and whether point "p" is in the tetrahedron.
c---- Loop over subset of data.
do 130 n = 1, ns
nn = n + n1 - 1
pabc(nn) = sign (pabc(nn), wdl(n))
pbcd(nn) = sign (pbcd(nn), wal(n))
pcda(nn) = sign (pcda(nn), wbl(n))
pdab(nn) = sign (pdab(nn), wcl(n))
pmin(nn) = amin1 (pabc(nn), pbcd(nn), pcda(nn), pdab(nn))
if (pmin(nn) .ge. 0.0) then
nloc(nn) = 1
else
nloc(nn) = 0
endif
if ((dista(n) .eq. 0.0) .or.
& (distb(n) .eq. 0.0) .or.
& (distc(n) .eq. 0.0) .or.
& (distd(n) .eq. 0.0) ) then
nloc(nn) = -1
endif
c---- End of loop over subset of data.
130 continue
c.... See if the local tetrahedral coordinates are to be returned.
c---- Return local coordinates.
if (kloc .eq. 1) then
c---- Loop over subset of data.
do 140 n = 1, ns
nn = n + n1 - 1
wa(nn) = wal(n)
wb(nn) = wbl(n)
wc(nn) = wcl(n)
wd(nn) = wdl(n)
c---- End of loop over subset of data.
140 continue
c---- Tested kloc.
endif
c.... See if all data subsets are done.
c---- Do another subset of data.
if (n2 .lt. np) then
n1 = n2 + 1
n2 = min (np, n1 + 63)
go to 110
endif
cbugc***DEBUG begins.
cbug 9902 format (/ 'aptetrn results:' /
cbug & (i3,' pmin= ',1pe22.14,' nloc=',i2 /
cbug & ' pabc,pbcd=',1p2e22.14 /
cbug & ' pcda,pdab=',1p2e22.14))
cbug 9903 format (/ 'aptetrn results: wa, wb, wc, wd' /
cbug & (i3,' wa,wb=',1p2e22.14 /
cbug & ' wc,wd=',1p2e22.14))
cbug
cbug write ( 3, 9902) (n, pmin(n), nloc(n), pabc(n), pbcd(n),
cbug & pcda(n), pdab(n), n = 1, np)
cbug
cbug if (kloc .eq. 1) then
cbug write ( 3, 9903) (n, wa(n), wb(n), wc(n), wd(n), n = 1, np)
cbug endif
cbugc***DEBUG ends.
210 return
c.... End of subroutine aptetrn. (+1 line.)
end
UCRL-WEB-209832