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