subroutine aptetrl (px, py, pz, ax, ay, az, bx, by, bz, & cx, cy, cz, dx, dy, dz, np, tol, & wa, wb, wc, wd, nerr) ccbeg. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE APTETRL c c call aptetrl (px, py, pz, ax, ay, az, bx, by, bz, c & cx, cy, cz, dx, dy, dz, np, tol, c & wa, wb, wc, wd, nerr) c c Version: aptetrl Updated 1990 May 11 12:40. c aptetrl Originated 1990 May 11 12:40. c c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123. c c c Purpose: To find, for each of np points p = (px, py, pz), the local c coordinates or vertex weights wa, wb, wc, wd in the tetrahedron c with vertices a = (ax, ay, az), b = (bx, by, bz), c c = (cx, cy, cz) and d = (dx, dy, dz) in 3-D space, such that: c c p = wa * a + wb * b + wc * c + wd * d c c Flag nerr indicates any input error. c c Input: px, py, pz, ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, c np, tol. c c Output: wa, wb, wc, wd, nerr. c c Calls: aptript 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 nerr Output Indicates an input error, if not 0. c 1 if np is not positive. c c np Input Size of arrays px, py, pz, ax, ay, az, bx, by, bz, c cx, cy, cz, dx, dy, dz, wa, wb, wc, wd. c c px,py,pz Input The x, y, z coordinates of point "p". Size np. c c tol Input Numerical tolerance limit. c On Cray computers, recommend 1.e-5 to 1.e-11. c c wa Output Fractional distance of point "p" from triangular face c "bcd" to vertex "a". Weight of vertex "a". Size np. c c wb Output Fractional distance of point "p" from triangular face c "cda" to vertex "b". Weight of vertex "b". Size np. c c wc Output Fractional distance of point "p" from triangular face c "dab" to vertex "c". Weight of vertex "c". Size np. c c wd Output Fractional distance of point "p" from triangular face c "abc" to vertex "d". Weight of vertex "d". Size np. c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ccend. c.... Dimensioned arguments. c---- Coordinate 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---- Coordinate 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---- Coordinate 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---- Coordinate 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---- 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---- Fract. dist. of "p" from "bcd" to "a". dimension wa (1) c---- Fract. dist. of "p" from "cda" to "b". dimension wb (1) c---- Fract. dist. of "p" from "dab" to "c". dimension wc (1) c---- Fract. dist. of "p" from "abc" to "d". dimension wd (1) c.... Local variables. c---- Distance from "bcd" to vertex "a". common /laptetrl/ dista (64) c---- Distance from "cda" to vertex "b". common /laptetrl/ distb (64) c---- Distance from "dab" to vertex "c". common /laptetrl/ distc (64) c---- Distance from "abc" to vertex "d". common /laptetrl/ distd (64) c---- Distance from "bcd" to point "p". common /laptetrl/ distpa (64) c---- Distance from "cda" to point "p". common /laptetrl/ distpb (64) c---- Distance from "dab" to point "p". common /laptetrl/ distpc (64) c---- Distance from "abc" to point "p". common /laptetrl/ distpd (64) c---- A very small number. common /laptetrl/ fuz c---- Index in local array. common /laptetrl/ n c---- First index of subset of data. common /laptetrl/ n1 c---- Last index of subset of data. common /laptetrl/ n2 c---- Index in external array. common /laptetrl/ nn c---- Size of current subset of data. common /laptetrl/ ns cbugc***DEBUG begins. cbug common /laptetrl/ avgwa, avgwb, avgwc, avgwd cbug common /laptetrl/ devwa, devwb, devwc, devwd cbug common /laptetrl/ avgwa2, avgwb2, avgwc2, avgwd2 cbug 9901 format (/ 'aptetrl finding local coordinates in a tetrahedron.', cbug & ' np=',i3 / cbug & (i3,' px,py,pz=',1p3e22.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) np, (n, px(n), py(n), pz(n), ax(n), ay(n), az(n), cbug & bx(n), by(n), bz(n), cx(n), cy(n), cz(n), dx(n), dy(n), dz(n), cbug & 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 distances of the tetrahedron vertices from the opposite faces. call aptript (ax(n1), ay(n1), az(n1), bx(n1), by(n1), bz(n1), & cx(n1), cy(n1), cz(n1), dx(n1), dy(n1), dz(n1), & ns, tol, dista, nerr) call aptript (bx(n1), by(n1), bz(n1), cx(n1), cy(n1), cz(n1), & dx(n1), dy(n1), dz(n1), ax(n1), ay(n1), az(n1), & ns, tol, distb, nerr) call aptript (cx(n1), cy(n1), cz(n1), dx(n1), dy(n1), dz(n1), & ax(n1), ay(n1), az(n1), bx(n1), by(n1), bz(n1), & ns, tol, distc, nerr) call aptript (dx(n1), dy(n1), dz(n1), ax(n1), ay(n1), az(n1), & bx(n1), by(n1), bz(n1), cx(n1), cy(n1), cz(n1), & ns, tol, distd, nerr) c.... Find the distances of point "p" from the tetrahedron triangular faces. call aptript (px(n1), py(n1), pz(n1), bx(n1), by(n1), bz(n1), & cx(n1), cy(n1), cz(n1), dx(n1), dy(n1), dz(n1), & ns, tol, distpa, nerr) call aptript (px(n1), py(n1), pz(n1), cx(n1), cy(n1), cz(n1), & dx(n1), dy(n1), dz(n1), ax(n1), ay(n1), az(n1), & ns, tol, distpb, nerr) call aptript (px(n1), py(n1), pz(n1), dx(n1), dy(n1), dz(n1), & ax(n1), ay(n1), az(n1), bx(n1), by(n1), bz(n1), & ns, tol, distpc, nerr) call aptript (px(n1), py(n1), pz(n1), ax(n1), ay(n1), az(n1), & bx(n1), by(n1), bz(n1), cx(n1), cy(n1), cz(n1), & ns, tol, distpd, nerr) c.... Find the local tetrahedral coordinates (vertex weights) of point "p". c---- Loop over subset of data. do 120 n = 1, ns nn = n + n1 - 1 wa(nn) = distpa(n) / (dista(n) + fuz) wb(nn) = distpb(n) / (distb(n) + fuz) wc(nn) = distpc(n) / (distc(n) + fuz) wd(nn) = distpd(n) / (distd(n) + fuz) c---- End of loop over subset of data. 120 continue 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 (/ 'aptetrl results:') cbug 9903 format (i3,' wa,wb,wc,wd=',1p2e22.14 / 16x,1p2e22.14) cbug cbug write ( 3, 9902) cbug write ( 3, 9903) (n, wa(n), wb(n), wc(n), wd(n), n = 1, np) 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 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 9904 format (/ 'aptetrl mean and deviation:' / cbug & ' avgw(a,b,c,d)=',1p2e22.14 / 16x,1p2e22.14 / cbug & ' devw(a,b,c,d)=',1p2e22.14 / 16x,1p2e22.14) cbug cbug write ( 3, 9904) avgwa, avgwb, avgwc, avgwd, cbug & devwa, devwb, devwc, devwd cbug cbugc***DEBUG ends. 210 return c.... End of subroutine aptetrl. (+1 line.) end UCRL-WEB-209832