subroutine aptrcut (ax, ay, az, bx, by, bz, cx, cy, cz, tol,
& ncut, acut1, fcut1, dcut1, ex, ey, ez,
& acut2, fcut2, dcut2, fx, fy, fz, dcut3, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTRCUT
c
c call aptrcut (ax, ay, az, bx, by, bz, cx, cy, cz, tol,
c & ncut, acut1, fcut1, dcut1, ex, ey, ez,
c & acut2, fcut2, dcut2, fx, fy, fz, dcut3, nerr)
c
c Version: aptrcut Updated 1999 September 13 14:00.
c aptrcut Originated 1999 April 19 15:30.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: For the triangle with vertices A = (ax, ay, az),
c B = (bx, by, bz), and C = (cx, cy, cz), and opposite edges
c a, b and c, in 3-D space, to find the endpoints E = (ex, ey, ez)
c and F = (fx, fy, fz) of the ncut (1 to 3) straight line(s) that
c cut the triangle into sections with equal perimeters and equal
c areas, and acut1 and acut2, the names of the edges or vertices
c on which points E and F occur, to find the fractional
c distances fcut1 and fcut2 of the end points along the edges,
c and to find the lengths dcut1, dcut2 and dcut3 of the small
c triangular section.
c Flag nerr indicates any input error.
c
c Input: ax, ay, az, bx, by, bz, cx, cy, cz, tol.
c
c Output: ncut, acut1, fcut1, ex, ey, ez, acut2, fcut2, fx, fy, fz, nerr.
c
c Calls: aptqrts, aptrich, aptvdis, aptvsum
c
c
c Glossary:
c
c acut1 Output The names of the vertices or edges on which points
c E occur: A, B, C, a, b or c. Size 3.
c
c acut2 Output The names of the vertices or edges on which points
c F occur: A, B, C, a, b or c. Size 3.
c
c ax,ay,az Input The x, y, z coordinates of vertex A of the triangle.
c
c bx,by,bz Input The x, y, z coordinates of vertex B of the triangle.
c
c cx,cy,cz Input The x, y, z coordinates of vertex C of the triangle.
c
c ex,ey,ez Output The x, y, z coordinates of point E. Size 3.
c
c dcut1 Output The distance of the first endput, E, along the edge.
c Size 3.
c
c dcut2 Output The distance of the second endput, F, along the edge.
c Size 3.
c
c dcut3 Output The length of the cutting line, from E to F. Size 3.
c
c fcut1 Output The fractional distance of point E along the edge.
c Size 3.
c
c fcut2 Output The fractional distance of point F along the edge.
c Size 3.
C
c fx,fy,fz Output The x, y, z coordinates of point F. Size 3.
c
c ncut Output The number of cutting lines (1, 2 or 3).
c
c nerr Output Indicates an input error, if not 0.
c 1 if the triangle has too small an area (vertices are
c almost coincident or colinear).
c
c tol Input Numerical tolerance limit.
c On Cray computers, recommend 1.e-5 to 1.e-11.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
implicit none
c.... Arguments.
real ax ! Coordinate x of triangle vertex A.
real ay ! Coordinate y of triangle vertex A.
real az ! Coordinate z of triangle vertex A.
real bx ! Coordinate x of triangle vertex B.
real by ! Coordinate y of triangle vertex B.
real bz ! Coordinate z of triangle vertex B.
real cx ! Coordinate x of triangle vertex C.
real cy ! Coordinate y of triangle vertex C.
real cz ! Coordinate z of triangle vertex C.
real dcut1(3) ! Distance of point E along edge.
real dcut2(3) ! Distance of point F along edge.
real dcut3(3) ! Length of cutting line.
real ex(3) ! Coordinate x of point E.
real ey(3) ! Coordinate y of point E.
real ez(3) ! Coordinate z of point E.
real fcut1(3) ! Fractional distance of E along edge.
real fcut2(3) ! Fractional distance of F along edge.
real fx(3) ! Coordinate x of point F.
real fy(3) ! Coordinate y of point F.
real fz(3) ! Coordinate z of point F.
integer ncut ! Number of equipartition cutting lines.
integer nerr ! Error indicator.
real tol ! Numerical tolerence limit.
character*1 acut1(3) ! Names of edges or vertices for point E.
character*1 acut2(3) ! Names of edges or vertices for point F.
c.... Local variables.
real a ! Length of edge vector BC (edge a).
real abx ! Component x of vector AB (edge c).
real aby ! Component y of vector AB (edge c).
real abz ! Component z of vector AB (edge c).
real b ! Length of edge vector CA (edge b).
real bcx ! Component x of vector BC (edge a).
real bcy ! Component y of vector BC (edge a).
real bcz ! Component z of vector BC (edge a).
real c ! Length of edge vector AB (edge c).
real c0 ! Coefficient of 1 in quadratic eq.
real c1 ! Coefficient of t in quadratic eq.
real c2 ! Coefficient of t**2 in quadratic eq.
real cax ! Component x of vector CA (edge b).
real cay ! Component y of vector CA (edge b).
real caz ! Component z of vector CA (edge b).
real clen ! Distance of end of cut line from origin.
real fab ! Fractional distance from A to B.
real fac ! Fractional distance from A to C.
real fba ! Fractional distance from B to A.
real fbc ! Fractional distance from B to C.
real fca ! Fractional distance from C to A.
real fcb ! Fractional distance from C to B.
integer itrun ! Truncation flag from aptqrts.
integer kvera ! Cuts through vertex A.
integer kverb ! Cuts through vertex B.
integer kverc ! Cuts through vertex C.
integer n ! Index of cut line.
integer nerra ! Error flag from apt subroutine.
integer nroots ! # of real roots of quadratic equation.
real per ! Perimeter of triangle, a + b + c.
real qq ! c1**2 - 4*c0*c2.
real root1 ! First real root of quadatic equation.
real root2 ! Second real root of quadratic equation.
character*1 acut1x ! Names of edges or vertices for point E.
character*1 acut2x ! Names of edges or vertices for point F.
cbugcbugc***DEBUG begins.
cbugcbug real aex, aey, aez, afx, afy, afz, dae, daf
cbugcbug real bex, bey, bez, bfx, bfy, bfz, dbe, dbf
cbugcbug real cex, cey, cez, cfx, cfy, cfz, dce, dcf
cbugcbug real efx, efy, efz, def
cbugcbugc***DEBUG ends.
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptrcut finding triangle partition cuts.',
cbug & ' tol=',1pe22.14 /
cbug & ' ax,ay,az= ',1p3e22.14 /
cbug & ' bx,by,bz= ',1p3e22.14 /
cbug & ' cx,cy,cz= ',1p3e22.14 )
cbug write ( 3, 9901) tol, ax, ay, az, bx, by, bz, cx, cy, cz
cbugc***DEBUG ends.
c.... Initialize.
nerr = 0
ncut = 0
kvera = 0
kverb = 0
kverc = 0
do 110 n = 1, 3
acut1(n) = '?'
acut2(n) = '?'
ex(n) = -1.e99
ey(n) = -1.e99
ez(n) = -1.e99
fx(n) = -1.e99
fy(n) = -1.e99
fz(n) = -1.e99
110 continue
c=======================================================================********
c.... Find the edge vectors of the triangle, BC, CA and AB,
c.... and the edge lengths, a, b and c.
call aptvdis (bx, by, bz, cx, cy, cz, 1, tol,
& bcx, bcy, bcz, a, nerra)
call aptvdis (cx, cy, cz, ax, ay, az, 1, tol,
& cax, cay, caz, b, nerra)
call aptvdis (ax, ay, az, bx, by, bz, 1, tol,
& abx, aby, abz, c, nerra)
c.... Test for a good triangle.
call aptrich (a, b, c, tol, nerr)
if (nerr .ne. 0) go to 410
per = a + b + c ! Perimeter.
cbugcbugc***DEBUG begins.
cbugcbug 9902 format (/
cbugcbug & 'aptrcut results (edge vectors and lengths):' /
cbugcbug & ' abx,aby,abz= ',1p3e22.14 /
cbugcbug & ' bcx,bcy,bcz= ',1p3e22.14 /
cbugcbug & ' cax,cay,caz= ',1p3e22.14 /
cbugcbug & ' a, b, c = ',1p3e22.14 )
cbugcbug write ( 3, 9902) abx, aby, abz, bcx, bcy, bcz, cax, cay, caz,
cbugcbug & a, b, c
cbugcbugc***DEBUG ends.
c.... Find lines that cut triangle into equal areas, perimeters.
c=======================================================================********
c.... Try vertex C, find cuts from edge a to edge b.
c2 = 2.0 * a
c1 = -per
c0 = b
call aptqrts (0, c2, c1, c0, qq, tol,
& nroots, root1, root2, itrun)
if (nroots .ge. 1) then
if ((root1 .ge. (0.5 - tol)) .and.
& (root1 .le. (1.0 + tol))) then
acut1x = 'a'
acut2x = 'b'
fcb = root1
fca = 0.5 / fcb
fbc = 1.0 - fcb
if (abs (fbc) .le. tol) then
if (kverb .gt. 0) go to 150
acut1x = 'B'
kverb = 1
fbc = 0.0
fcb = 1.0
fac = 0.5
fca = 0.5
endif
fac = 1.0 - fca
if (abs (fac) .le. tol) then
if (kvera .gt. 0) go to 150
acut2x = 'A'
kvera = 1
fac = 0.0
fca = 1.0
fbc = 0.5
fcb = 0.5
endif
ncut = ncut + 1
acut1(ncut) = acut1x
acut2(ncut) = acut2x
fcut1(ncut) = fcb
fcut2(ncut) = fca
dcut1(ncut) = fcb * a
dcut2(ncut) = fca * b
dcut3(ncut) = sqrt (per * (c - 0.25 * per))
call aptvsum (0, fbc, cx, cy, cz, fcb, bx, by, bz, 1, tol,
& ex(ncut), ey(ncut), ez(ncut), clen, nerra)
call aptvsum (0, fac, cx, cy, cz, fca, ax, ay, az, 1, tol,
& fx(ncut), fy(ncut), fz(ncut), clen, nerra)
cbugcbugc***DEBUG begins.
cbugcbug call aptvdis (cx, cy, cz, ex(ncut), ey(ncut), ez(ncut), 1, tol,
cbugcbug $ cex, cey, cez, dce, nerra)
cbugcbug call aptvdis (cx, cy, cz, fx(ncut), fy(ncut), fz(ncut), 1, tol,
cbugcbug $ cfx, cfy, cfz, dcf, nerra)
cbugcbug call aptvdis (ex(ncut), ey(ncut), ez(ncut),
cbugcbug & fx(ncut), fy(ncut), fz(ncut), 1, tol,
cbugcbug $ efx, efy, efz, def, nerra)
cbugcbug 9301 format ('fac, fab =',1p2e22.14 /
cbugcbug & 'dae, daf, def=',1p3e22.14)
cbugcbug 9302 format ('fba, fbc =',1p2e22.14 /
cbugcbug & 'dbe, dbf, def=',1p3e22.14)
cbugcbug 9303 format ('fcb, fca =',1p2e22.14 /
cbugcbug & 'dce, dcf, def=',1p3e22.14)
cbugcbug write ( 3, 9303) fcb, fca, dce, dcf, def
cbugcbugc***DEBUG ends.
endif ! Cut endpoints are on edges.
endif ! Cut endpoints are real.
150 if (nroots .eq. 2) then
if ((root2 .ge. (0.5 - tol)) .and.
& (root2 .le. (1.0 + tol))) then
acut1x = 'a'
acut2x = 'b'
fcb = root2
fca = 0.5 / fcb
fbc = 1.0 - fcb
if (abs (fbc) .le. tol) then
if (kverb .gt. 0) go to 190
kverb = 1
acut1x = 'B'
fbc = 0.0
fcb = 1.0
fca = 0.5
fac = 0.5
endif
fac = 1.0 - fca
if (abs (fac) .le. tol) then
if (kvera .gt. 0) go to 190
kvera = 1
acut2x = 'A'
fac = 0.0
fca = 1.0
fbc = 0.5
fcb = 0.5
endif
ncut = ncut + 1
acut1(ncut) = acut1x
acut2(ncut) = acut2x
fcut1(ncut) = fcb
fcut2(ncut) = fca
dcut1(ncut) = fcb * a
dcut2(ncut) = fca * b
dcut3(ncut) = sqrt (per * (c - 0.25 * per))
call aptvsum (0, fbc, cx, cy, cz, fcb, bx, by, bz, 1, tol,
& ex(ncut), ey(ncut), ez(ncut), clen, nerra)
call aptvsum (0, fac, cx, cy, cz, fca, ax, ay, az, 1, tol,
& fx(ncut), fy(ncut), fz(ncut), clen, nerra)
cbugcbugc***DEBUG begins.
cbugcbug call aptvdis (cx, cy, cz, ex(ncut), ey(ncut), ez(ncut), 1, tol,
cbugcbug $ cex, cey, cez, dce, nerra)
cbugcbug call aptvdis (cx, cy, cz, fx(ncut), fy(ncut), fz(ncut), 1, tol,
cbugcbug $ cfx, cfy, cfz, dcf, nerra)
cbugcbug call aptvdis (ex(ncut), ey(ncut), ez(ncut),
cbugcbug & fx(ncut), fy(ncut), fz(ncut), 1, tol,
cbugcbug $ efx, efy, efz, def, nerra)
cbugcbug write ( 3, 9303) fcb, fca, dce, dcf, def
cbugcbugc***DEBUG ends.
endif ! Cut endpoints are on edges.
endif ! Cut endpoints are real.
190 continue
c=======================================================================********
c.... Try vertex A, find cuts from edge b to edge c.
c2 = 2.0 * b
c1 = -per
c0 = c
call aptqrts (0, c2, c1, c0, qq, tol,
& nroots, root1, root2, itrun)
if (nroots .ge. 1) then
if ((root1 .ge. (0.5 - tol)) .and.
& (root1 .le. (1.0 + tol))) then
acut1x = 'b'
acut2x = 'c'
fac = root1
fab = 0.5 / fac
fca = 1.0 - fac
if (abs (fca) .le. tol) then
if (kverc .gt. 0) go to 250
kverc = 1
acut1x = 'C'
fac = 1.0
fca = 0.0
fab = 0.5
fba = 0.5
endif
fba = 1.0 - fab
if (abs (fba) .le. tol) then
if (kverb .gt. 0) go to 250
kverb = 1
acut2x = 'B'
fab = 1.0
fba = 0.0
fac = 0.5
fca = 0.5
endif
ncut = ncut + 1
acut1(ncut) = acut1x
acut2(ncut) = acut2x
fcut1(ncut) = fac
fcut2(ncut) = fab
dcut1(ncut) = fac * b
dcut2(ncut) = fab * c
dcut3(ncut) = sqrt (per * (a - 0.25 * per))
call aptvsum (0, fca, ax, ay, az, fac, cx, cy, cz, 1, tol,
& ex(ncut), ey(ncut), ez(ncut), clen, nerra)
call aptvsum (0, fba, ax, ay, az, fab, bx, by, bz, 1, tol,
& fx(ncut), fy(ncut), fz(ncut), clen, nerra)
cbugcbugc***DEBUG begins.
cbugcbug call aptvdis (ax, ay, az, ex(ncut), ey(ncut), ez(ncut), 1, tol,
cbugcbug $ aex, aey, aez, dae, nerra)
cbugcbug call aptvdis (ax, ay, az, fx(ncut), fy(ncut), fz(ncut), 1, tol,
cbugcbug $ afx, afy, afz, daf, nerra)
cbugcbug call aptvdis (ex(ncut), ey(ncut), ez(ncut),
cbugcbug & fx(ncut), fy(ncut), fz(ncut), 1, tol,
cbugcbug $ efx, efy, efz, def, nerra)
cbugcbug write ( 3, 9301) fac, fab, dae, daf, def
cbugcbugc***DEBUG ends.
endif ! Cut endpoints are on edges.
endif ! Cut endpoints are real.
250 if (nroots .eq. 2) then
if ((root2 .ge. (0.5 - tol)) .and.
& (root2 .le. (1.0 + tol))) then
acut1x = 'b'
acut2x = 'c'
fac = root2
fab = 0.5 / fac
fca = 1.0 - fac
if (abs (fca) .le. tol) then
if (kverc .gt. 0) go to 290
kverc = 1
acut1x = 'C'
fac = 1.0
fca = 0.0
fab = 0.5
fba = 0.5
endif
fba = 1.0 - fab
if (abs (fba) .le. tol) then
if (kverb .gt. 0) go to 290
kverb = 1
acut2x = 'B'
fab = 1.0
fba = 0.0
fac = 0.5
fca = 0.5
endif
ncut = ncut + 1
acut1(ncut) = acut1x
acut2(ncut) = acut2x
fcut1(ncut) = fac
fcut2(ncut) = fab
dcut1(ncut) = fac * b
dcut2(ncut) = fab * c
dcut3(ncut) = sqrt (per * (a - 0.25 * per))
call aptvsum (0, fca, ax, ay, az, fac, cx, cy, cz, 1, tol,
& ex(ncut), ey(ncut), ez(ncut), clen, nerra)
call aptvsum (0, fba, ax, ay, az, fab, bx, by, bz, 1, tol,
& fx(ncut), fy(ncut), fz(ncut), clen, nerra)
cbugcbugc***DEBUG begins.
cbugcbug call aptvdis (ax, ay, az, ex(ncut), ey(ncut), ez(ncut), 1, tol,
cbugcbug $ aex, aey, aez, dae, nerra)
cbugcbug call aptvdis (ax, ay, az, fx(ncut), fy(ncut), fz(ncut), 1, tol,
cbugcbug $ afx, afy, afz, daf, nerra)
cbugcbug call aptvdis (ex(ncut), ey(ncut), ez(ncut),
cbugcbug & fx(ncut), fy(ncut), fz(ncut), 1, tol,
cbugcbug $ efx, efy, efz, def, nerra)
cbugcbug write ( 3, 9301) fac, fab, dae, daf, def
cbugcbugc***DEBUG ends.
endif ! Cut endpoints are on edges.
endif ! Cut endpoints are real.
290 continue
c=======================================================================********
c.... Try vertex B, find cuts from edge c to edge a.
c2 = 2.0 * c
c1 = -per
c0 = a
call aptqrts (0, c2, c1, c0, qq, tol,
& nroots, root1, root2, itrun)
if (nroots .ge. 1) then
if ((root1 .ge. (0.5 - tol)) .and.
& (root1 .le. (1.0 + tol))) then
acut1x = 'c'
acut2x = 'a'
fba = root1
fbc = 0.5 / fba
fab = 1.0 - fba
if (abs (fab) .le. tol) then
if (kvera .gt. 0) go to 350
kvera = 1
acut1x = 'A'
fab = 0.0
fba = 1.0
fbc = 0.5
fcb = 0.5
endif
fcb = 1.0 - fbc
if (abs (fcb) .le. tol) then
if (kverc .gt. 0) go to 350
acut2x = 'C'
kverc = 1
fbc = 1.0
fcb = 0.0
fab = 0.5
fba = 0.5
endif
ncut = ncut + 1
acut1(ncut) = acut1x
acut2(ncut) = acut2x
fcut1(ncut) = fba
fcut2(ncut) = fbc
dcut1(ncut) = fba * c
dcut2(ncut) = fbc * a
dcut3(ncut) = sqrt (per * (b - 0.25 * per))
call aptvsum (0, fab, bx, by, bz, fba, ax, ay, az, 1, tol,
& ex(ncut), ey(ncut), ez(ncut), clen, nerra)
call aptvsum (0, fcb, bx, by, bz, fbc, cx, cy, cz, 1, tol,
& fx(ncut), fy(ncut), fz(ncut), clen, nerra)
cbugcbugc***DEBUG begins.
cbugcbug call aptvdis (bx, by, bz, ex(ncut), ey(ncut), ez(ncut), 1, tol,
cbugcbug $ bex, bey, bez, dbe, nerra)
cbugcbug call aptvdis (bx, by, bz, fx(ncut), fy(ncut), fz(ncut), 1, tol,
cbugcbug $ bfx, bfy, bfz, dbf, nerra)
cbugcbug call aptvdis (ex(ncut), ey(ncut), ez(ncut),
cbugcbug & fx(ncut), fy(ncut), fz(ncut), 1, tol,
cbugcbug $ efx, efy, efz, def, nerra)
cbugcbug write ( 3, 9302) fba, fbc, dbe, dbf, def
cbugcbugc***DEBUG ends.
endif ! Cut endpoints are on edges.
endif ! Cut endpoints are real.
350 if (nroots .eq. 2) then
if ((root2 .ge. (0.5 - tol)) .and.
& (root2 .le. (1.0 + tol))) then
acut1x = 'c'
acut2x = 'a'
fba = root2
fbc = 0.5 / fba
fab = 1.0 - fba
if (abs (fab) .le. tol) then
if (kvera .gt. 0) go to 390
kvera = 1
acut1x = 'A'
fab = 0.0
fba = 1.0
fbc = 0.5
fcb = 0.5
endif
fcb = 1.0 - fbc
if (abs (fcb) .le. tol) fcb = 0.0
if (abs (fcb) .le. tol) then
if (kverc .gt. 0) go to 390
kverc = 1
acut2x = 'C'
fbc = 1.0
fcb = 0.0
fab = 0.5
fba = 0.5
endif
ncut = ncut + 1
acut1(ncut) = acut1x
acut2(ncut) = acut2x
fcut1(ncut) = fba
fcut2(ncut) = fbc
dcut1(ncut) = fba * c
dcut2(ncut) = fbc * a
dcut3(ncut) = sqrt (per * (b - 0.25 * per))
call aptvsum (0, fab, bx, by, bz, fba, ax, ay, az, 1, tol,
& ex(ncut), ey(ncut), ez(ncut), clen, nerra)
call aptvsum (0, fcb, bx, by, bz, fbc, cx, cy, cz, 1, tol,
& fx(ncut), fy(ncut), fz(ncut), clen, nerra)
cbugcbugc***DEBUG begins.
cbugcbug call aptvdis (bx, by, bz, ex(ncut), ey(ncut), ez(ncut), 1, tol,
cbugcbug $ bex, bey, bez, dbe, nerra)
cbugcbug call aptvdis (bx, by, bz, fx(ncut), fy(ncut), fz(ncut), 1, tol,
cbugcbug $ bfx, bfy, bfz, dbf, nerra)
cbugcbug call aptvdis (ex(ncut), ey(ncut), ez(ncut),
cbugcbug & fx(ncut), fy(ncut), fz(ncut), 1, tol,
cbugcbug $ efx, efy, efz, def, nerra)
cbugcbug write ( 3, 9302) fba, fbc, dbe, dbf, def
cbugcbugc***DEBUG ends.
endif ! Cut endpoints are on edges.
endif ! Cut endpoints are real.
390 continue
c=======================================================================********
410 continue
cbugc***DEBUG begins.
cbug 9916 format (/ 'aptrcut results. nerr = ',i2,' ncut = ',i1 )
cbug write ( 3, 9916) nerr, ncut
cbug
cbug 9918 format ('ncut = ',i2,' acut1 = ',a1,' acut2 = ',a1 /
cbug & 'fcut1, fcut2 =',1p2e22.14 /
cbug & 'ex, ey, ez =',1p3e22.14 /
cbug & 'fx, fy, fz =',1p3e22.14 )
cbug
cbug do 420 n = 1, ncut
cbug write ( 3, 9918) n, acut1(n), acut2(n), fcut1(n), fcut2(n),
cbug & ex(n), ey(n), ez(n), fx(n), fy(n), fz(n)
cbug 420 continue
cbugc***DEBUG ends.
return
c.... End of subroutine aptrcut. (+1 line.)
end
UCRL-WEB-209832