subroutine apttetd (ax, ay, az, bx, by, bz, cx, cy, cz, & dx, dy, dz, np, tol, vtet, gx, gy, gz, nerr) ccbeg. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE APTTETD c c call apttetd (ax, ay, az, bx, by, bz, cx, cy, cz, c & dx, dy, dz, np, tol, vtet, gx, gy, gz, nerr) c c Version: apttetd Updated 1990 December 3 16:20. c apttetd Originated 1990 May 8 16:40. c c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123. c c c Purpose: To find the volumes vtet and the centers of gravity c g = (gx, gy, gz) of np tetrahedrons with vertices c a = (ax, ay, az), b = (bx, by, bz), c = (cx, cy, cz) and c d = (dx, dy, dz) in 3-D space. c Flag nerr indicates any input error. c c Input: ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, np, tol. c c Output: vtet, gx, gy, gz, nerr. c c Calls: aptvaxb, aptvdis, aptvdot 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 gx,gy,gz Output Center of gravity or centroid of tetrahedron. Size np. c Coordinates may be truncated to zero, if less than c the estimated error in their calculation, based on c tol. c c nerr Output Indicates an input error, if not 0. c 1 if np is not positive. c c np Input Size of external arrays. c c tol Input Numerical tolerance limit. c On Cray computers, recommend 1.e-5 to 1.e-11. c c vtet Output The volume of the tetrahedron. Size np. c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ccend. c.... Dimensioned arguments. c---- Coordinate x of vertex "a". dimension ax (1) c---- Coordinate y of vertex "a". dimension ay (1) c---- Coordinate z of vertex "a". dimension az (1) c---- Coordinate x of vertex "b". dimension bx (1) c---- Coordinate y of vertex "b". dimension by (1) c---- Coordinate z of vertex "b". dimension bz (1) c---- Coordinate x of vertex "c". dimension cx (1) c---- Coordinate y of vertex "c". dimension cy (1) c---- Coordinate z of vertex "c". dimension cz (1) c---- Coordinate x of vertex "d". dimension dx (1) c---- Coordinate y of vertex "d". dimension dy (1) c---- Coordinate z of vertex "d". dimension dz (1) c---- Cooodinate x of centroid "g". dimension gx (1) c---- Coordinate y of centroid "g". dimension gy (1) c---- Coordinate z of centroid "g". dimension gz (1) c---- Volume of tetrahedron. dimension vtet (1) c.... Local variables. c---- Component x of edge vector "ab". common /lapttetd/ abx (64) c---- Component y of edge vector "ab". common /lapttetd/ aby (64) c---- Component z of edge vector "ab". common /lapttetd/ abz (64) c---- Component x of edge vector "ac". common /lapttetd/ acx (64) c---- Component y of edge vector "ac". common /lapttetd/ acy (64) c---- Component z of edge vector "ac". common /lapttetd/ acz (64) c---- Component x of edge vector "ad". common /lapttetd/ adx (64) c---- Component y of edge vector "ad". common /lapttetd/ ady (64) c---- Component z of edge vector "ad". common /lapttetd/ adz (64) c---- Component x of vector "e". common /lapttetd/ ex (64) c---- Component y of vector "e". common /lapttetd/ ey (64) c---- Component z of vector "e". common /lapttetd/ ez (64) c---- Estimated error in gx. common /lapttetd/ gxerr c---- Estimated error in gy. common /lapttetd/ gyerr c---- Estimated error in gz. common /lapttetd/ gzerr c---- Index in local array. common /lapttetd/ n c---- First index of subset of data. common /lapttetd/ n1 c---- Last index of subset of data. common /lapttetd/ n2 c---- Index in external array. common /lapttetd/ nn c---- Size of current subset of data. common /lapttetd/ ns c---- Length of a vector. common /lapttetd/ vlen (64) cbugc***DEBUG begins. cbug 9901 format (/ 'apttetd finding tetrahedron data:' / cbug & (i3,' 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, 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. 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) 110 ns = n2 - n1 + 1 c.... Find the edge vectors of the tetrahedron, "ab", "ac" and "ad". call aptvdis (ax(n1), ay(n1), az(n1), bx(n1), by(n1), bz(n1), & ns, tol, abx, aby, abz, vlen, nerr) call aptvdis (ax(n1), ay(n1), az(n1), cx(n1), cy(n1), cz(n1), & ns, tol, acx, acy, acz, vlen, nerr) call aptvdis (ax(n1), ay(n1), az(n1), dx(n1), dy(n1), dz(n1), & ns, tol, adx, ady, adz, vlen, nerr) c.... Find six times the volume of the tetrahedron. call aptvaxb (abx, aby, abz, acx, acy, acz, ns, tol, & ex, ey, ez, vlen, nerr) call aptvdot (ex, ey, ez, adx, ady, adz, ns, tol, & vtet(n1), nerr) c.... Find the volume and the coordinates of the center of gravity. c---- Loop over subset of data. do 120 n = 1, ns nn = n + n1 - 1 vtet(nn) = abs (vtet(nn)) / 6.0 gx(nn) = 0.25 * (ax(nn) + bx(nn) + cx(nn) + dx(nn)) gy(nn) = 0.25 * (ay(nn) + by(nn) + cy(nn) + dy(nn)) gz(nn) = 0.25 * (az(nn) + bz(nn) + cz(nn) + dz(nn)) 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 gxerr = 0.25 * tol * (abs (ax(nn)) + abs (bx(nn)) + & abs (cx(nn)) + abs (dx(nn))) gyerr = 0.25 * tol * (abs (ay(nn)) + abs (by(nn)) + & abs (cy(nn)) + abs (dy(nn))) gzerr = 0.25 * tol * (abs (az(nn)) + abs (bz(nn)) + & abs (cz(nn)) + abs (dz(nn))) if (abs (gx(nn)) .lt. gxerr) then gx(nn) = 0.0 endif if (abs (gy(nn)) .lt. gyerr) then gy(nn) = 0.0 endif if (abs (gz(nn)) .lt. gzerr) then gz(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 (/ 'apttetd results:' / cbug & (i3,' vtet= ',1pe22.14 / cbug & ' gx,gy,gz=',1p3e22.14)) cbug cbug write ( 3, 9902) (n, vtet(n), cbug & gx(n), gy(n), gz(n), n = 1, np) cbugc***DEBUG ends. 210 return c.... End of subroutine apttetd. (+1 line.) end UCRL-WEB-209832