24 subroutine addnod ( mpl, nst, k, x, y, z, list, lptr, lend, lnew, ier )
140 character(len=1024),
parameter :: subr =
'addnod'
146 write ( *,
'(a)' )
' '
147 write ( *,
'(a)' )
'ADDNOD - Fatal error!'
148 write ( *,
'(a)' )
' K < 4.'
149 call mpl%abort(subr,
'stop in stripack')
169 call trfind ( ist, p, km1, x, y, z, list, lptr, lend, b1, b2, b3, &
174 if ( i1 == 0 )
call mpl%abort(subr,
'colinear points')
179 if (
eq(p(1),x(l)) .and.
eq(p(2),y(l)) .and.
eq(p(3),z(l)) )
then
186 if (
eq(p(1),x(l)) .and.
eq(p(2),y(l)) .and.
eq(p(3),z(l)) )
then
192 if (
eq(p(1),x(l)) .and.
eq(p(2),y(l)) .and.
eq(p(3),z(l)) )
then
197 call intadd ( kk, i1, i2, i3, list, lptr, lend, lnew )
202 call bdyadd ( kk, i1,i2, list, lptr, lend, lnew )
204 call covsph ( kk, i1, list, lptr, lend, lnew )
217 io1 = abs( list(lpo1) )
223 call lstptr ( lend(io1), io2, list, lptr, lp )
225 if ( 0 <= list(lp) )
then
228 in1 = abs( list(lp) )
235 call swptst ( in1, kk, io1, io2, x, y, z, output )
236 if ( .not. output )
then
238 if ( lpo1 == lpf .or. list(lpo1) < 0 )
then
244 io1 = abs( list(lpo1) )
249 call swap ( in1, kk, io1, io2, list, lptr, lend, lpo1 )
255 if ( lpo1 /= 0 )
then
266 if ( lpo1 == lpf .or. list(lpo1) < 0 )
then
272 io1 = abs( list(lpo1) )
278 subroutine bdyadd ( kk, i1, i2, list, lptr, lend, lnew )
380 call insert ( k, lp, list, lptr, lnew )
382 if ( next == n2 )
then
395 lptr(lnew) = lnew + 1
401 if ( next == n2 )
then
406 lptr(lnew) = lnew + 1
420 subroutine bnodes ( n, list, lptr, lend, nodes, nb, na, nt )
486 integer list(6*(n-2))
488 integer lptr(6*(n-2))
507 if ( list(lp) < 0 )
then
539 if ( n0 == nst )
then
556 subroutine covsph ( kk, n0, list, lptr, lend, lnew )
636 call insert ( k, lp, list, lptr, lnew )
640 if ( next == nst )
then
654 lptr(lnew) = lnew + 1
658 if ( next == nst )
then
669 subroutine det ( x1, y1, z1, x2, y2, z2, x0, y0, z0, output )
682 real ( kind_real ) x1
683 real ( kind_real ) y1
684 real ( kind_real ) z1
685 real ( kind_real ) x2
686 real ( kind_real ) y2
687 real ( kind_real ) z2
688 real ( kind_real ) x0
689 real ( kind_real ) y0
690 real ( kind_real ) z0
691 real ( kind_real ) t1
692 real ( kind_real ) t2
693 real ( kind_real ) t3
694 real ( kind_real ) output
696 t1 = x0*(y1*z2-y2*z1)
697 t2 = y0*(x1*z2-x2*z1)
698 t3 = z0*(x1*y2-x2*y1)
699 output = t1 - t2 + t3
702 if (
small(output,t1).or.
small(output,t2).or.
small(output,t3)) output = 0.0
704 subroutine insert ( k, lp, list, lptr, lnew )
763 function inside ( p, lv, xv, yv, zv, nv, listv, ier )
913 real (
kind_real ),
parameter :: eps = 0.001_kind_real
954 if ( n < 3 .or. imx < n )
then
964 if ( i1 < 1 .or. imx < i1 )
then
990 if ( i2 < 1 .or. imx < i2 )
then
995 vn(1) = yv(i1) * zv(i2) - zv(i1) * yv(i2)
996 vn(2) = zv(i1) * xv(i2) - xv(i1) * zv(i2)
997 vn(3) = xv(i1) * yv(i2) - yv(i1) * xv(i2)
998 vnrm = sqrt( sum( vn(1:3)**2 ) )
1000 if (
eq(vnrm,0.0_kind_real) )
then
1004 q(1) = xv(i1) + xv(i2) + eps * vn(1) / vnrm
1005 q(2) = yv(i1) + yv(i2) + eps * vn(2) / vnrm
1006 q(3) = zv(i1) + zv(i2) + eps * vn(3) / vnrm
1008 qnrm = sqrt( sum( q(1:3)**2 ) )
1016 cn(1) = q(2) * p(3) - q(3) * p(2)
1017 cn(2) = q(3) * p(1) - q(1) * p(3)
1018 cn(3) = q(1) * p(2) - q(2) * p(1)
1020 if (
eq(cn(1),0.0_kind_real) .and.
eq(cn(2),0.0_kind_real) .and.
eq(cn(2),0.0_kind_real) )
then
1024 pn(1) = p(2) * cn(3) - p(3) * cn(2)
1025 pn(2) = p(3) * cn(1) - p(1) * cn(3)
1026 pn(3) = p(1) * cn(2) - p(2) * cn(1)
1027 qn(1) = cn(2) * q(3) - cn(3) * q(2)
1028 qn(2) = cn(3) * q(1) - cn(1) * q(3)
1029 qn(3) = cn(1) * q(2) - cn(2) * q(1)
1041 if ( i2 < 1 .or. imx < i2 )
then
1046 lft2 =
sup(cn(1) * xv(i2) + cn(2) * yv(i2) + cn(3) * zv(i2),0.0_kind_real)
1056 if ( i2 < 1 .or. imx < i2 )
then
1061 lft2 =
sup(cn(1) * xv(i2) + cn(2) * yv(i2) + cn(3) * zv(i2),0.0_kind_real)
1063 if ( lft1 .eqv. lft2 )
then
1078 call intrsc ( v1, v2, cn, b, ierr )
1084 if (
sup(dot_product( b(1:3), qn(1:3) ),0.0_kind_real) .and. &
1085 sup(dot_product( b(1:3), pn(1:3) ),0.0_kind_real) )
then
1090 d = dot_product( b(1:3), q(1:3) )
1092 if (
sup(d,bq) )
then
1097 d = dot_product( b(1:3), p(1:3) )
1099 if (
sup(d,bq) )
then
1110 if ( ni /= 2 * ( ni / 2 ) .or. .not. qinr )
then
1116 if ( pinr .neqv. even )
then
1126 subroutine intadd ( kk, i1, i2, i3, list, lptr, lend, lnew )
1206 call lstptr ( lend(n1), n2, list, lptr, lp )
1207 call insert ( k, lp, list, lptr, lnew )
1209 call lstptr ( lend(n2), n3, list, lptr, lp )
1210 call insert ( k, lp, list, lptr, lnew )
1212 call lstptr ( lend(n3), n1, list, lptr, lp )
1213 call insert ( k, lp, list, lptr, lnew )
1220 lptr(lnew) = lnew + 1
1221 lptr(lnew+1) = lnew + 2
1292 real ( kind_real ) cn(3)
1293 real ( kind_real ) d1
1294 real ( kind_real ) d2
1296 real ( kind_real ) p(3)
1297 real ( kind_real ) p1(3)
1298 real ( kind_real ) p2(3)
1299 real ( kind_real ) pp(3)
1300 real ( kind_real ) ppn
1301 real ( kind_real ) t
1303 d1 = dot_product( cn(1:3), p1(1:3) )
1304 d2 = dot_product( cn(1:3), p2(1:3) )
1306 if (
eq(d1,d2) )
then
1313 t = d1 / ( d1 - d2 )
1315 pp(1:3) = p1(1:3) + t * ( p2(1:3) - p1(1:3) )
1317 ppn = dot_product( pp(1:3), pp(1:3) )
1321 if (
eq(ppn,0.0_kind_real) )
then
1330 p(1:3) = pp(1:3) / ppn
1336 subroutine jrand ( n, ix, iy, iz, output )
1386 real ( kind_real ) u
1387 real ( kind_real ) x
1389 ix = mod( 171 * ix, 30269 )
1390 iy = mod( 172 * iy, 30307 )
1391 iz = mod( 170 * iz, 30323 )
1393 x = ( real( ix, kind_real ) / 30269.0_kind_real ) &
1394 + ( real( iy, kind_real ) / 30307.0_kind_real ) &
1395 + ( real( iz, kind_real ) / 30323.0_kind_real )
1398 output = int( real( n, kind_real ) * u ) + 1
1401 end subroutine jrand
1402 subroutine left ( x1, y1, z1, x2, y2, z2, x0, y0, z0, output )
1446 real ( kind_real ) x0
1447 real ( kind_real ) x1
1448 real ( kind_real ) x2
1449 real ( kind_real ) y0
1450 real ( kind_real ) y1
1451 real ( kind_real ) y2
1452 real ( kind_real ) z0
1453 real ( kind_real ) z1
1454 real ( kind_real ) z2
1455 real ( kind_real ) zz
1460 call det ( x1, y1, z1, x2, y2, z2, x0, y0, z0, zz )
1461 output =
sup(zz,0.0_kind_real)
1465 subroutine lstptr ( lpl, nb, list, lptr, output )
1530 if ( nd == nb )
then
1536 if ( lp == lpl )
then
1546 subroutine swap ( in1, in2, io1, io2, list, lptr, lend, lp21 )
1613 call lstptr ( lend(in1), in2, list, lptr, lp )
1615 if ( abs( list(lp) ) == in2 )
then
1622 call lstptr ( lend(io1), in2, list, lptr, lp )
1624 lptr(lp) = lptr(lph)
1628 if ( lend(io1) == lph )
then
1634 call lstptr ( lend(in1), io1, list, lptr, lp )
1642 call lstptr ( lend(io2), in1, list, lptr, lp )
1644 lptr(lp) = lptr(lph)
1648 if ( lend(io2) == lph )
then
1654 call lstptr ( lend(in2), io2, list, lptr, lp )
1663 subroutine swptst ( n1, n2, n3, n4, x, y, z, output )
1718 real ( kind_real ) dx1
1719 real ( kind_real ) dx2
1720 real ( kind_real ) dx3
1721 real ( kind_real ) dy1
1722 real ( kind_real ) dy2
1723 real ( kind_real ) dy3
1724 real ( kind_real ) dz1
1725 real ( kind_real ) dz2
1726 real ( kind_real ) dz3
1732 real ( kind_real ) x(*)
1733 real ( kind_real ) x4
1734 real ( kind_real ) y(*)
1735 real ( kind_real ) y4
1736 real ( kind_real ) z(*)
1737 real ( kind_real ) z4
1738 real ( kind_real ) zz
1758 call det ( dx2, dy2, dz2, dx1, dy1, dz1, dx3, dy3, dz3, zz )
1759 output =
sup(zz,0.0_kind_real)
1763 subroutine trfind ( nst, p, n, x, y, z, list, lptr, lend, b1, b2, b3, i1, &
1869 integer ( kind = 4 ),
save :: ix = 1
1870 integer ( kind = 4 ),
save :: iy = 2
1871 integer ( kind = 4 ),
save :: iz = 3
1873 integer list(6*(n-2))
1875 integer lptr(6*(n-2))
1907 if ( n0 < 1 .or. n < n0 )
then
1908 call jrand ( n, ix, iy, iz, n0 )
1913 eps = epsilon( eps )
1914 tol = 100.0_kind_real * eps
1935 call det(x(n0),y(n0),z(n0),x(n1),y(n1),z(n1),xp,yp,zp,output)
1936 if (
inf(output,0.0_kind_real) )
then
1939 if ( n1 == nl )
then
1952 call det(x(n0),y(n0),z(n0),x(nf),y(nf),z(nf), xp,yp,zp,output)
1953 if (
inf(output,0.0_kind_real) )
then
1961 call det(x(nl),y(nl),z(nl),x(n0),y(n0),z(n0),xp,yp,zp,output)
1962 if (
inf(output,0.0_kind_real) )
then
1976 n2 = abs( list(lp) )
1978 call det(x(n0),y(n0),z(n0),x(n2),y(n2),z(n2),xp,yp,zp,output)
1979 if (
inf(output,0.0_kind_real) )
then
1985 if ( n1 /= nl )
then
1989 call det ( x(n0), y(n0), z(n0), x(nf), y(nf), z(nf), xp, yp, zp, output )
1990 if (
inf(output,0.0_kind_real) )
then
1997 if (
inf(abs( x(n0 ) * xp + y(n0) * yp + z(n0) * zp),1.0-4.0*eps) )
then
2005 call det(x(n1),y(n1),z(n1),x(n0),y(n0),z(n0),xp,yp,zp,output)
2006 if (
inf(output,0.0_kind_real) )
then
2011 n1 = abs( list(lp) )
2013 if ( n1 == nl )
then
2049 call det ( x(n1),y(n1),z(n1),x(n2),y(n2),z(n2),xp,yp,zp,b3 )
2051 if (
inf(b3,0.0_kind_real) )
then
2056 call lstptr ( lend(n2), n1, list, lptr, lp )
2058 if ( list(lp) < 0 )
then
2063 n4 = abs( list(lp) )
2067 call det ( x(n0),y(n0),z(n0),x(n4),y(n4),z(n4),xp,yp,zp,output )
2068 if (
inf(output,0.0_kind_real) )
then
2072 if ( n2 /= n2s .and. n2 /= n0 )
then
2079 if ( n1 /= n1s .and. n1 /= n0 )
then
2088 call jrand ( n, ix, iy, iz, n0 )
2096 if (
supeq(b3,eps) )
then
2100 call det(x(n2),y(n2),z(n2),x(n3),y(n3),z(n3),xp,yp,zp,b1)
2101 call det(x(n3),y(n3),z(n3),x(n1),y(n1),z(n1),xp,yp,zp,b2)
2105 if (
inf(b1,-tol) .or.
inf(b2,-tol) )
then
2106 call jrand ( n, ix, iy, iz, n0 )
2116 s12 = x(n1) * x(n2) + y(n1) * y(n2) + z(n1) * z(n2)
2117 ptn1 = xp * x(n1) + yp * y(n1) + zp * z(n1)
2118 ptn2 = xp * x(n2) + yp * y(n2) + zp * z(n2)
2119 b1 = ptn1 - s12 * ptn2
2120 b2 = ptn2 - s12 * ptn1
2124 if (
inf(b1,-tol) .or.
inf(b2,-tol) )
then
2125 call jrand ( n, ix, iy, iz, n0 )
2136 b1 = max( b1, 0.0_kind_real )
2137 b2 = max( b2, 0.0_kind_real )
2158 call det(x(n2),y(n2),z(n2),x(next),y(next),z(next),xp,yp,zp,output)
2159 if (
supeq(output,0.0_kind_real) )
then
2164 s12 = x(n1) * x(n2) + y(n1) * y(n2) + z(n1) * z(n2)
2166 q(1) = x(n1) - s12 * x(n2)
2167 q(2) = y(n1) - s12 * y(n2)
2168 q(3) = z(n1) - s12 * z(n2)
2170 if (
sup(xp * q(1) + yp * q(2) + zp * q(3),0.0_kind_real) )
then
2174 if (
sup(x(next) * q(1) + y(next) * q(2) + z(next) * q(3),0.0_kind_real) )
then
2189 if ( n2 /= n1s )
then
2221 call det ( x(next), y(next), z(next), x(n1), y(n1), z(n1), xp, yp, zp, output )
2222 if (
supeq(output,0.0_kind_real) )
then
2227 s12 = x(n1) * x(n2) + y(n1) * y(n2) + z(n1) * z(n2)
2228 q(1) = x(n2) - s12 * x(n1)
2229 q(2) = y(n2) - s12 * y(n1)
2230 q(3) = z(n2) - s12 * z(n1)
2232 if (
sup(xp * q(1) + yp * q(2) + zp * q(3),0.0_kind_real) )
then
2236 if (
sup(x(next) * q(1) + y(next) * q(2) + z(next) * q(3),0.0_kind_real) )
then
2251 if ( n1 /= n1s )
then
2278 subroutine trlist ( n, list, lptr, lend, nrow, nt, ltri, ier )
2379 integer list(6*(n-2))
2384 integer lptr(6*(n-2))
2385 integer ltri(nrow,*)
2394 if ( n < 3 .or. ( nrow /= 6 .and. nrow /= 9 ) )
then
2427 n3 = abs( list(lp) )
2429 if ( n2 < n1 .or. n3 < n1 )
then
2448 else if ( i == 2 )
then
2465 if ( list(lp) == i2 )
then
2471 if ( lp == lpl )
then
2480 if ( abs( list(lp) ) /= i2 )
then
2491 if ( list(lp) < 0 )
then
2501 i3 = abs( list(lp) )
2506 if ( i1 < i2 .and. i1 < i3 )
then
2508 else if ( i2 < i3 )
then
2532 if ( ltri(1,kn) == i1 .and. &
2533 ltri(2,kn) == i2 .and. &
2534 ltri(3,kn) == i3 )
then
2567 if ( lp2 /= lpln1 )
then
2578 subroutine trmesh ( mpl, n, x, y, z, list, lptr, lend, lnew, near, next, dist, ier )
2738 type(
mpl_type),
intent(inout) :: mpl
2752 logical output_1, output_2
2754 integer list(6*(n-2))
2758 integer lptr(6*(n-2))
2766 character(len=1024),
parameter :: subr =
'trmesh'
2772 write ( *,
'(a)' )
' '
2773 write ( *,
'(a)' )
'TRMESH - Fatal error!'
2774 write ( *,
'(a)' )
' N < 3.'
2775 call mpl%abort(subr,
'stop in stripack')
2786 call left (x(1),y(1),z(1),x(2),y(2),z(2),x(3),y(3),z(3),output_1)
2787 call left (x(2),y(2),z(2),x(1),y(1),z(1),x(3),y(3),z(3),output_2)
2789 if ( .not. output_1 )
then
2811 else if ( .not. output_2 )
then
2884 d1 = -( x(k) * x(1) + y(k) * y(1) + z(k) * z(1) )
2885 d2 = -( x(k) * x(2) + y(k) * y(2) + z(k) * z(2) )
2886 d3 = -( x(k) * x(3) + y(k) * y(3) + z(k) * z(3) )
2888 if (
inf(d1,d2) .and.
inf(d1,d3) )
then
2893 else if (
inf(d2,d1) .and.
inf(d2,d3) )
then
2910 call addnod ( mpl, near(k), k, x, y, z, list, lptr, lend, lnew, ier )
2912 if ( ier /= 0 )
then
2913 write ( *,
'(a)' )
'TRMESH - ADDNOD - Fatal error!'
2914 write ( *,
'(a,i8,a,i8)' )
' Node ', k,
' is equal to already used node ', ier
2915 call mpl%abort(subr,
'stop in stripack')
2923 if ( near(i) == k )
then
2977 d = - ( x(i) * x(k) + y(i) * y(k) + z(i) * z(k) )
2978 if (
inf(d,dist(i)) )
then
2987 if ( i == near(j) )
then
3005 if ( lp /= lpl )
then