subroutine aptripc (pu, pv, au, av, bu, bv, cu, cv, np, tol,
& wa, wb, wc, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTRIPC
c
c call aptripc (pu, pv, au, av, bu, bv, cu, cv, np, tol,
c & wa, wb, wc, nerr)
c
c Version: aptripc Updated 1990 May 14 10:10.
c aptripc Originated 1990 May 14 10:10.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To find the local coordinates wa, wb and wc of the np points
c p = (pu, pv) in the triangles with vertices a = (au, av),
c b = (bu, bv) and c = (cu, cv) in the uv plane.
c
c p = wa * a + wb * b + wc * c
c wa + wb + wc = 1
c
c Flag nerr indicates any input error.
c
c Input: pu, pv, au, av, bu, bv, cu, cv, np, tol.
c
c Output: wa, wb, wc, nerr.
c
c Calls: aptptlc
c
c Glossary:
c
c au, av Input The u and v coordinates of vertex "a". Size np.
c
c bu, bv Input The u and v coordinates of vertex "b". Size np.
c
c cu, cv Input The u and v coordinates of vertex "c". Size np.
c
c nerr Output Indicates an input error, if not 0.
c
c np Input Size of arrays pu, pv, au, av bu, bv, cu, cv,
c wa, wb, wc.
c
c pu, pv Input The u and v coordinates of point "p". Size np.
c
c tol Input Numerical tolerance limit for wa, wb, wc.
c On Cray computers, recommend 1.e-5 to 1.e-11.
c
c wa Output Fractional distance of point "p" from edge "bc" to
c vertex "a". Size np. Meaningless if triangle "abc"
c has a zero area.
c
c wb Output Fractional distance of point "p" from edge "ca" to
c vertex "b". Size np. Meaningless if triangle "abc"
c has a zero area.
c
c wc Output Fractional distance of point "p" from edge "ab" to
c vertex "c". Size np. Meaningless if triangle "abc"
c has a zero area.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments.
c---- Coordinate u of point "a".
dimension au (1)
c---- Coordinate v of point "a".
dimension av (1)
c---- Coordinate u of point "b".
dimension bu (1)
c---- Coordinate v of point "b".
dimension bv (1)
c---- Coordinate u of point "c".
dimension cu (1)
c---- Coordinate v of point "c".
dimension cv (1)
c---- Coordinate u of point "p".
dimension pu (1)
c---- Coordinate v of point "p".
dimension pv (1)
c---- Local weight of vertex "a".
dimension wa (1)
c---- Local weight of vertex "b".
dimension wb (1)
c---- Local weight of vertex "c".
dimension wc (1)
c.... Local variables.
c---- Distance from edge "ab" to point "c".
common /laptripc/ dabc (64)
c---- Distance from edge "ab" to point "p".
common /laptripc/ dabp (64)
c---- Distance from edge "bc" to point "a".
common /laptripc/ dbca (64)
c---- Distance from edge "bc" to point "p".
common /laptripc/ dbcp (64)
c---- Distance from edge "ca" to point "b".
common /laptripc/ dcab (64)
c---- Distance from edge "ca" to point "p".
common /laptripc/ dcap (64)
c---- Dummy argument from aptptlc.
common /laptripc/ fdmin
c---- A very small number.
common /laptripc/ fuz
c---- Dummy argument from aptptlc.
common /laptripc/ itrun (64)
c---- Index in local array.
common /laptripc/ n
c---- First index of subset of data.
common /laptripc/ n1
c---- Last index of subset of data.
common /laptripc/ n2
c---- Dummy argument from aptptlc.
common /laptripc/ nlim
c---- Index in external array.
common /laptripc/ nn
c---- Size of current subset of data.
common /laptripc/ ns
cbugc***DEBUG begins.
cbug common /laptripc/ avgwa, avgwb, avgwc
cbug common /laptripc/ devwa, devwb, devwc
cbug common /laptripc/ avgwa2, avgwb2, avgwc2
cbug 9901 format (/ 'aptripc finding local coordinates in triangle.' /
cbug & (i3,' pu,pv=',1p2e22.14 /
cbug & ' au,av=',1p2e22.14 /
cbug & ' bu,bv=',1p2e22.14 /
cbug & ' cu,cv=',1p2e22.14))
cbug write ( 3, 9901) (n, pu(n), pv(n), au(n), av(n), bu(n), bv(n),
cbug & cu(n), cv(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 distances of point "p" from the triangle edges.
call aptptlc (pu(n1), pv(n1), au(n1), av(n1), bu(n1), bv(n1),
& ns, tol, -1, dabp, fdmin, nlim, itrun, nerr)
call aptptlc (pu(n1), pv(n1), bu(n1), bv(n1), cu(n1), cv(n1),
& ns, tol, -1, dbcp, fdmin, nlim, itrun, nerr)
call aptptlc (pu(n1), pv(n1), cu(n1), cv(n1), au(n1), av(n1),
& ns, tol, -1, dcap, fdmin, nlim, itrun, nerr)
c.... Find the distances of each vertex from the opposite triangle edges.
call aptptlc (cu(n1), cv(n1), au(n1), av(n1), bu(n1), bv(n1),
& ns, tol, -1, dabc, fdmin, nlim, itrun, nerr)
call aptptlc (au(n1), av(n1), bu(n1), bv(n1), cu(n1), cv(n1),
& ns, tol, -1, dbca, fdmin, nlim, itrun, nerr)
call aptptlc (bu(n1), bv(n1), cu(n1), cv(n1), au(n1), av(n1),
& ns, tol, -1, dcab, fdmin, nlim, itrun, nerr)
c.... Find the local triangular coordinates (vertex weights) of point "p".
c---- Loop over subset of data.
do 120 n = 1, ns
nn = n + n1 - 1
wa(nn) = dbcp(n) / (dbca(n) + fuz)
wb(nn) = dcap(n) / (dcab(n) + fuz)
wc(nn) = dabp(n) / (dabc(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 (/ 'aptripc results:')
cbug 9903 format (i3,' wa,wb,wc=',1p3e22.14)
cbug
cbug write ( 3, 9902)
cbug write ( 3, 9903) (n, wa(n), wb(n), wc(n), n = 1, np)
cbug
cbug avgwa = 0.0
cbug avgwb = 0.0
cbug avgwc = 0.0
cbug avgwa2 = 0.0
cbug avgwb2 = 0.0
cbug avgwc2 = 0.0
cbug
cbug do 130 n = 1, np
cbug avgwa = avgwa + wa(n)
cbug avgwb = avgwb + wb(n)
cbug avgwc = avgwc + wc(n)
cbug avgwa2 = avgwa2 + wa(n)**2
cbug avgwb2 = avgwb2 + wb(n)**2
cbug avgwc2 = avgwc2 + wc(n)**2
cbug 130 continue
cbug
cbug avgwa = avgwa / np
cbug avgwb = avgwb / np
cbug avgwc = avgwc / np
cbug avgwa2 = avgwa2 / np
cbug avgwb2 = avgwb2 / np
cbug avgwc2 = avgwc2 / np
cbug devwa = sqrt (avgwa2 - avgwa**2)
cbug devwb = sqrt (avgwb2 - avgwb**2)
cbug devwc = sqrt (avgwc2 - avgwc**2)
cbug
cbug 9904 format (/ 'aptripc mean and deviation:' /
cbug & ' avgw(a,b,c)=',1p3e22.14 /
cbug & ' devw(a,b,c)=',1p3e22.14)
cbug
cbug write ( 3, 9904) avgwa, avgwb, avgwc,
cbug & devwa, devwb, devwc
cbug
cbugc***DEBUG ends.
210 return
c.... End of subroutine aptripc. (+1 line.)
end
UCRL-WEB-209832