subroutine aptetrp (ax, ay, az, bx, by, bz, cx, cy, cz, & dx, dy, dz, np, kloc, tol, & px, py, pz, wa, wb, wc, wd, nerr) ccbeg. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE APTETRP c c call aptetrp (ax, ay, az, bx, by, bz, cx, cy, cz, c & dx, dy, dz, np, kloc, tol, c & px, py, pz, wa, wb, wc, wd, nerr) c c Version: aptetrp Updated 1990 November 27 14:00. c aptetrp Originated 1990 May 9 9:30. c c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123. c c c Purpose: To randomly sample np points p = (px, py, pz) from a uniform c spatial distribution inside the tetrahedron 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, and optionally (kloc = 1), c to find the local coordinates or vertex weights wa, wb, c wc, wd, corresponding to the point "p" in the tetrahedron: c c p = wa * a + wb * b + wc * c + wd * d c c Coordinates of "p" less than the estimated error in their c calculation, based on tol, will be truncated to zero. c Flag nerr indicates any input error. c c Input: ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, np, c kloc, tol. c c Output: px, py, pz, wa, wb, wc, wd, nerr. c c Calls: aptspod, apttlod c c c Glossary: c c ax,ay,az Input The x, y, z coordinates of tetrahedron vertex "a". c c bx,by,bz Input The x, y, z coordinates of tetrahedron vertex "b". c c cx,cy,cz Input The x, y, z coordinates of tetrahedron vertex "c". c c dx,dy,dz Input The x, y, z coordinates of tetrahedron vertex "d". c c kloc Input Option flag to return local tetrahedron coordinates c for each point "p", if 1. 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 np Input Size of arrays px, py, pz. c Number of points "p" to sample. c c px,py,pz Output The coordinates of spatial points "p" randomly sampled c from a uniform distribution over the volume of the c tetrahedron. Size np. Values less than the c estimated error in their calculation, based on tol, c will be truncated to zero. c c tol Input Numerical tolerance limit. c On Cray computers, recommend 1.e-5 to 1.e-11. c c wa, ... Output Fractional distances of point "p" from the opposite c triangular faces to the tetrahedron vertices "a", ... c Size np only if kloc = 1. Otherwise, not returned. c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ccend. c.... Dimensioned arguments. 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---- Fractional distance of "p" to "a". dimension wa (1) c---- Fractional distance of "p" to "b". dimension wb (1) c---- Fractional distance of "p" to "c". dimension wc (1) c---- Fractional distance of "p" to "d". dimension wd (1) c.... Local variables. c---- Component x of vector "e". common /laptetrp/ ex (64) c---- Component y of vector "e". common /laptetrp/ ey (64) c---- Component z of vector "e". common /laptetrp/ ez (64) c---- Fractional distance from "d" to "abc". common /laptetrp/ fdis (64) c---- Index in local array. common /laptetrp/ n c---- First index of subset of data. common /laptetrp/ n1 c---- Last index of subset of data. common /laptetrp/ n2 c---- Index in external array. common /laptetrp/ nn c---- Size of current subset of data. common /laptetrp/ ns c---- Estimated error in px. common /laptetrp/ pxerr c---- Estimated error in py. common /laptetrp/ pyerr c---- Estimated error in pz. common /laptetrp/ pzerr c---- Weight of triangle vertex "a". common /laptetrp/ twa (64) c---- Weight of triangle vertex "b". common /laptetrp/ twb (64) c---- Weight of triangle vertex "c". common /laptetrp/ twc (64) cbugc***DEBUG begins. cbug common /laptetrp/ avgx, avgy, avgz, devx, devy, devz cbug common /laptetrp/ avgx2, avgy2, avgz2 cbug common /laptetrp/ avgwa, avgwb, avgwc, avgwd cbug common /laptetrp/ devwa, devwb, devwc, devwd cbug common /laptetrp/ avgwa2, avgwb2, avgwc2, avgwd2 cbug 9901 format (/ 'aptetrp sampling inside a tetrahedron. np=',i3 / 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, ax, ay, az, bx, by, bz, cbug & cx, cy, cz, dx, dy, dz cbugc***DEBUG ends. c.... initialize. 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.... Randomly sample points in the triangle "abc". call apttlod (ax, ay, az, bx, by, bz, cx, cy, cz, ns, 1, tol, & ex, ey, ez, twa, twb, twc, nerr) c.... Randomly sample fractional distances from a quadratic distribution. call aptspod (2., 0., 1., ns, fdis, nerr) c.... Find points "p" randomly sampled uniformly over the tetrahedron. c---- Loop over subset of data. do 120 n = 1, ns nn = n + n1 - 1 px(nn) = dx + fdis(n) * (ex(n) - dx) py(nn) = dy + fdis(n) * (ey(n) - dy) pz(nn) = dz + fdis(n) * (ez(n) - dz) c---- End of loop over subset of data. 120 continue c.... See if truncation to zero is 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 pxerr = tol * (abs (dx) + fdis(n) * abs (ex(n) - dx)) pyerr = tol * (abs (dy) + fdis(n) * abs (ey(n) - dy)) pzerr = tol * (abs (dz) + fdis(n) * abs (ez(n) - dz)) 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 local coordinates are to be returned. c---- Local coordinates needed. if (kloc .eq. 1) then c---- Loop over subset of data. do 140 n = 1, ns nn = n + n1 - 1 wa(nn) = fdis(n) * twa(n) wb(nn) = fdis(n) * twb(n) wc(nn) = fdis(n) * twc(n) wd(nn) = 1.0 - fdis(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 cbug 9902 format (/ 'aptetrp results:') cbug 9903 format (i3,' px,py,pz=',1p3e22.14) cbug 9904 format (/ 'aptetrp mean and deviation:' / cbug & ' avg(x,y,z)=',1p3e22.14 / cbug & ' dev(x,y,z)=',1p3e22.14) cbug 9905 format (i3,' wa,wb,wc,wd=',4f15.12) cbug 9906 format (/ 'aptetrp mean and deviation:' / cbug & ' avg(wa,wb,wc,wd)=',4f15.12 / cbug & ' dev(wa,wb,wc,wd)=',4f15.12) cbug cbug write ( 3, 9902) cbug cbug avgx = 0.0 cbug avgy = 0.0 cbug avgz = 0.0 cbug avgx2 = 0.0 cbug avgy2 = 0.0 cbug avgz2 = 0.0 cbug cbug do 150 n = 1, np cbug avgx = avgx + px(n) cbug avgy = avgy + py(n) cbug avgz = avgz + pz(n) cbug avgx2 = avgx2 + px(n)**2 cbug avgy2 = avgy2 + py(n)**2 cbug avgz2 = avgz2 + pz(n)**2 cbug write ( 3, 9903) n, px(n), py(n), pz(n) cbug 150 continue cbug cbug avgx = avgx / np cbug avgy = avgy / np cbug avgz = avgz / np cbug avgx2 = avgx2 / np cbug avgy2 = avgy2 / np cbug avgz2 = avgz2 / np cbug devx = sqrt (avgx2 - avgx**2) cbug devy = sqrt (avgy2 - avgy**2) cbug devz = sqrt (avgz2 - avgz**2) cbug cbug write ( 3, 9904) avgx, avgy, avgz, devx, devy, devz cbug cbugc---- Local coordinates returned. cbug if (kloc .eq. 1) then cbug cbug avgwa = 0.0 cbug avgwb = 0.0 cbug avgwc = 0.0 cbug avgwd = 0.0 cbug avgwa2 = 0.0 cbug avgwb2 = 0.0 cbug avgwc2 = 0.0 cbug avgwd2 = 0.0 cbug cbug do 160 n = 1, np cbug avgwa = avgwa + wa(n) cbug avgwb = avgwb + wb(n) cbug avgwc = avgwc + wc(n) cbug avgwd = avgwd + wd(n) cbug avgwa2 = avgwa2 + wa(n)**2 cbug avgwb2 = avgwb2 + wb(n)**2 cbug avgwc2 = avgwc2 + wc(n)**2 cbug avgwd2 = avgwd2 + wd(n)**2 cbug write ( 3, 9905) n, wa(n), wb(n), wc(n), wd(n) cbug 160 continue cbug cbug avgwa = avgwa / np cbug avgwb = avgwb / np cbug avgwc = avgwc / np cbug avgwd = avgwd / np cbug avgwa2 = avgwa2 / np cbug avgwb2 = avgwb2 / np cbug avgwc2 = avgwc2 / np cbug avgwd2 = avgwd2 / np cbug devwa = sqrt (avgwa2 - avgwa**2) cbug devwb = sqrt (avgwb2 - avgwb**2) cbug devwc = sqrt (avgwc2 - avgwc**2) cbug devwd = sqrt (avgwd2 - avgwd**2) cbug cbug write ( 3, 9906) avgwa, avgwb, avgwc, avgwd, cbug & devwa, devwb, devwc, devwd cbug cbugc---- Tested kloc. cbug endif cbugc***DEBUG ends. 210 return c.... End of subroutine aptetrp. (+1 line.) end UCRL-WEB-209832