subroutine aptrlop (wa, wb, wc, ax, ay, az, bx, by, bz,
& cx, cy, cz, np, tol, px, py, pz, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTRLOP
c
c call aptrlop (wa, wb, wc, ax, ay, az, bx, by, bz,
c & cx, cy, cz, np, tol, px, py, pz, nerr)
c
c Version: aptrlop Updated 1990 December 3 14:20.
c aptrlop Originated 1990 May 16 16: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 wa, wb and wc,
c the point p = (px, py, pz), interpolated between the vertices
c a = (ax, ay, az), b = (bx, by, bz) and c = (cx, cy, cz) of a
c triangle in 3-D xyz space:
c
c p = (wa * a + wb * b + wc * c) / (wa + wb + wc).
c
c Flag nerr indicates any input error.
c
c Input: wa, wb, wc, ax, ay, az, bx, by, bz, cx, cy, cz, 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 triangle.
c Size np.
c
c bx,by,bz Input The x, y, z coordinates of vertex "b" of triangle.
c Size np.
c
c cx,cy,cz Input The x, y, z coordinates of vertex "c" of triangle.
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, ax, ay, az, bx, by, bz,
c cx, cy, cz, 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 edge "bc", when normalized. Size np.
c
c wb Input Fractional distance of point "p" to vertex "b" from
c edge "ca", when normalized. Size np.
c
c wc Input Fractional distance of point "p" to vertex "c" from
c edge "ab", 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 "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 "bc".
dimension wa (1)
c---- Fractional distance to "b" from "ca".
dimension wb (1)
c---- Fractional distance to "c" from "ab".
dimension wc (1)
c.... Local variables.
c---- A very small number.
common /laptrlop/ fuz
c---- Index in local array.
common /laptrlop/ n
c---- First index of subset of data.
common /laptrlop/ n1
c---- Last index of subset of data.
common /laptrlop/ n2
c---- Index in external array.
common /laptrlop/ nn
c---- Size of current subset of data.
common /laptrlop/ ns
c---- Estimated error in px.
common /laptrlop/ pxerr
c---- Estimated error in py.
common /laptrlop/ pyerr
c---- Estimated error in pz.
common /laptrlop/ pzerr
c---- Sum of wa + wb + wc.
common /laptrlop/ sum
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptrlop interpolating in a triangle.' /
cbug & (i3,' wa,wb,wc=',1p3e22.14 /
cbug & ' ax,ay,az=',1p3e22.14 /
cbug & ' bx,by,bz=',1p3e22.14 /
cbug & ' cx,cy,cz=',1p3e22.14))
cbug write (3, 9901) (n, wa(n), wb(n), wc(n), ax(n), ay(n), az(n),
cbug & bx(n), by(n), bz(n), cx(n), cy(n), cz(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)
px(nn) = (wa(nn) * ax(nn) + wb(nn) * bx(nn) + wc(nn) * cx(nn)) /
& (sum + fuz)
py(nn) = (wa(nn) * ay(nn) + wb(nn) * by(nn) + wc(nn) * cy(nn)) /
& (sum + fuz)
pz(nn) = (wa(nn) * az(nn) + wb(nn) * bz(nn) + wc(nn) * cz(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)
pxerr = tol * (abs (wa(nn) * ax(nn)) +
& abs (wb(nn) * bx(nn)) +
& abs (wc(nn) * cx(nn))) / (sum + fuz)
pyerr = tol * (abs (wa(nn) * ay(nn)) +
& abs (wb(nn) * by(nn)) +
& abs (wc(nn) * cy(nn))) / (sum + fuz)
pzerr = tol * (abs (wa(nn) * az(nn)) +
& abs (wb(nn) * bz(nn)) +
& abs (wc(nn) * cz(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 (/ 'aptrlop 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 aptrlop. (+1 line.)
end
UCRL-WEB-209832