subroutine aptetrw (wa, wb, wc, wd, ax, ay, az, bx, by, bz,
& cx, cy, cz, dx, dy, dz, np, tol,
& px, py, pz, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTETRW
c
c call aptetrw (wa, wb, wc, wd, ax, ay, az, bx, by, bz, cx, cy, cz,
c & dx, dy, dz, np, tol, px, py, pz, nerr)
c
c Version: aptetrw Updated 1990 November 27 14:00.
c aptetrw Originated 1990 May 17 9:00.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To find, for each of np sets of vertex weights or local
c coordinates wa, wb, wc and wd, the point p = (px, py, pz),
c interpolated between the vertices a = (ax, ay, az),
c b = (bx, by, bz), c = (cx, cy, cz) and d = (dx, dy, dz)
c of a tetrahedron:
c
c p = (wa * a + wb * b + wc * c + wd * d) / (wa + wb + wc + wd).
c
c If points "a", "b", "c" and "d" are the vertices of a
c quadrangle, rather than a tetrahedron, the vertex weights
c wa, wb, wc and wd are related to the edge weights fdk and fdl
c (see aptqloc) as follows:
c
c wa = (1 - fdk) * (1 - fdl) fdk = wb + wc = 1 - wa - wd
c wb = fdk * (1 - fdl) fdl = wc + wd = 1 - wa - wb
c wc = fdk * fdl
c wd = (1 - fdk) * fdl wa + wb + wc + wd = 1
c
c where fdk is the fractional distance from edge "da" to edge "bc"
c and fdl is the fractional distance from edge "ab" to edge "cd".
c
c Flag nerr indicates any input error.
c
c Input: wa, wb, wc, wd, ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz,
c np, tol.
c
c Output: px, py, pz, nerr.
c
c Calls: None
c
c Glossary:
c
c ax,ay,az Input The x, y, z coordinates of vertex "a" of tetrahedron.
c Size np.
c
c bx,by,bz Input The x, y, z coordinates of vertex "b" of tetrahedron.
c Size np.
c
c cx,cy,cz Input The x, y, z coordinates of vertex "c" of tetrahedron.
c Size np.
c
c dx,dy,dz Input The x, y, z coordinates of vertex "d" of tetrahedron.
c Size np.
c
c nerr Output Indicates an input error, if not 0.
c 1 if np is not positive.
c
c np Input Size of arrays wa, wb, wc, wd, ax, ay, az, bx, by, bz,
c cx, cy, cz, dx, dy, dz, px, py, pz.
c
c px,py,pz Output Interpolated point "p". Size np. Coordinates may be
c truncated to zero, if less than the estimated error
c in their calculation, based on tol.
c
c tol Input Numerical tolerance limit.
c On Cray computers, recommend 1.e-5 to 1.e-11.
c
c wa Input Fractional distance of point "p" to vertex "a" from
c triangular face "bcd", when normalized. Size np.
c
c wb Input Fractional distance of point "p" to vertex "b" from
c triangular face "cda", when normalized. Size np.
c
c wc Input Fractional distance of point "p" to vertex "c" from
c triangular face "dab", when normalized. Size np.
c
c wd Input Fractional distance of point "p" to vertex "d" from
c triangular face "abc", when normalized. Size np.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments.
c---- Cooodinate 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---- Cooodinate 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---- Cooodinate 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---- Cooodinate 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---- Cooodinate 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---- Fractional distance to "a" from "bcd".
dimension wa (1)
c---- Fractional distance to "b" from "cda".
dimension wb (1)
c---- Fractional distance to "c" from "dab".
dimension wc (1)
c---- Fractional distance to "d" from "abc".
dimension wd (1)
c.... Local variables.
c---- A very small number.
common /laptetrw/ fuz
c---- Index in local array.
common /laptetrw/ n
c---- First index of subset of data.
common /laptetrw/ n1
c---- Last index of subset of data.
common /laptetrw/ n2
c---- Index in external array.
common /laptetrw/ nn
c---- Size of current subset of data.
common /laptetrw/ ns
c---- Estimated error in px.
common /laptetrw/ pxerr
c---- Estimated error in py.
common /laptetrw/ pyerr
c---- Estimated error in pz.
common /laptetrw/ pzerr
c---- Sum of wa + wb + wc + wd.
common /laptetrw/ sum
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptetrw interpolating in a tetrahedron.' /
cbug & (i3,' wa,wb= ',1p2e22.14 /
cbug & ' wc,wd= ',1p2e22.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) (n, wa(n), wb(n), wc(n), wd(n),
cbug & ax(n), ay(n), az(n), bx(n), by(n), bz(n),
cbug & cx(n), cy(n), cz(n), dx(n), dy(n), dz(n), 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
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 coordinates of the points.
c---- Loop over subset of data.
do 120 n = 1, ns
nn = n + n1 - 1
sum = wa(nn) + wb(nn) + wc(nn) + wd(nn)
px(nn) = (wa(nn) * ax(nn) + wb(nn) * bx(nn) +
& wc(nn) * cx(nn) + wd(nn) * dx(nn)) / (sum + fuz)
py(nn) = (wa(nn) * ay(nn) + wb(nn) * by(nn) +
& wc(nn) * cy(nn) + wd(nn) * dy(nn)) / (sum + fuz)
pz(nn) = (wa(nn) * az(nn) + wb(nn) * bz(nn) +
& wc(nn) * cz(nn) + wd(nn) * dz(nn)) / (sum + fuz)
c---- End of loop over subset of data.
120 continue
c.... See if truncation allowed.
c---- Truncation to zero allowed.
if (tol .gt. 0.0) then
c---- Loop over subset of data.
do 130 n = 1, ns
nn = n + n1 - 1
sum = wa(nn) + wb(nn) + wc(nn) + wd(nn)
pxerr = tol * (abs (wa(nn) * ax(nn)) +
& abs (wb(nn) * bx(nn)) +
& abs (wc(nn) * cx(nn)) +
& abs (wd(nn) * dx(nn))) / (sum + fuz)
pyerr = tol * (abs (wa(nn) * ay(nn)) +
& abs (wb(nn) * by(nn)) +
& abs (wc(nn) * cy(nn)) +
& abs (wd(nn) * dy(nn))) / (sum + fuz)
pzerr = tol * (abs (wa(nn) * az(nn)) +
& abs (wb(nn) * bz(nn)) +
& abs (wc(nn) * cz(nn)) +
& abs (wd(nn) * dz(nn))) / (sum + fuz)
if (abs(px(nn)) .lt. pxerr) then
px(nn) = 0.0
endif
if (abs(py(nn)) .lt. pyerr) then
py(nn) = 0.0
endif
if (abs(pz(nn)) .lt. pzerr) then
pz(nn) = 0.0
endif
c---- End of loop over subset of data.
130 continue
c---- Tested tol.
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 (/ 'aptetrw results:' /
cbug & (i3,' px,py,pz=',1p3e22.14))
cbug write ( 3, 9902) (n, px(n), py(n), pz(n), n = 1, np)
cbugc***DEBUG ends.
210 return
c.... End of subroutine aptetrw. (+1 line.)
end
UCRL-WEB-209832