/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:59 */
/*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */
#include
#include "fcrt.h"
#include "stgfi.h"
#include
/* PARAMETER translations */
#define ZERO 0.0e0
/* end of PARAMETER translations */
void /*FUNCTION*/ stgfi(
float x[],
float y[],
float z[],
float dz[][2],
long triang[],
long nt,
long bdry[][4-(1)+1],
long mb,
long ncont,
float q[],
float *zout,
LOGICAL32 wantdz,
float dzout[],
long *modefi,
float savwrk[])
{
LOGICAL32 newtri;
long int indsgn, indtri, modefnd, tri[8];
/* OFFSET Vectors w/subscript range: 1 to dimension */
float *const Dzout = &dzout[0] - 1;
float *const Q = &q[0] - 1;
float *const Savwrk = &savwrk[0] - 1;
long *const Tri = &tri[0] - 1;
float *const X = &x[0] - 1;
float *const Y = &y[0] - 1;
float *const Z = &z[0] - 1;
/* end of OFFSET VECTORS */
/* Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
* ALL RIGHTS RESERVED.
* Based on Government Sponsored Research NAS7-03001.
*>> 1997-07-01 STGFI Krogh Reversed subscripts in B (CLL suggestion)
*>> 1997-06-24 STGFI CLL Added names MODEFI and MODEFND.
*>> 1997-06-23 STGFI Krogh Removed "implicit none"
*>> 1997-06-22 STGFI CLL
*>> 1997-06-04 STGFI CLL
*>> 1996-02-02 STGFI CLL
*>> 1996-01-11 STGFI CLL
*>> 1995-10-30 STGFI CLL
*>> 1995-09-26 STGFI CLL Editing for inclusion into MATH77.
* This routine uses STGFND to do lookup, and either STGC0 or STGC1
* to do interpolation in a triangular grid, and uses STGEXT for
* extrapolation.
* ------------------------------------------------------------------
* Subroutine Arguments
*
* X(1:NP), Y(1:NP) [in] (x,y) coordinates of vertices of the
* triangular grid.
*
* Z(1:NP) [in] Z(i) is the value at (X(i),Y(i)) of the data to
* be interpolated.
*
* DZ(1:2, 1:NP) [in] DZ(1:2, i) are the values of the partial
* derivatives of the interpolation function at (X(i),Y(i)) with
* respect to x and y, respectively.
*
* {NP is the number of vertices in the triangular grid, but is
* not explicitly used in this subroutine.}
*
* TRIANG(1:6*NT) [in] Array of integer pointers defining the
* connectivity of the triangular grid.
*
* NT [in] No. of triangles in the triangular grid.
*
* Bdry(1:4, MB) [integer, in] Array containing pointers defining the
* boundary of the (convex) triangular grid.
* Bdry(1, K) = FWD POINTER. Points to next vertex in
* counterclockwise order.
* Bdry(2, K) = BACKWARD POINTER. Points to next vertex in
* clockwise order.
* Bdry(3, K) = A BOUNDARY POINT
* Bdry(4, K) = A BOUNDARY TRIANGLE
* The triangle Bdry(4, K) has a boundary edge that
* connects the points Bdry(3, K) and Bdry(3, K+1).
* In general not all elements of the array Bdry() are members of
* the linked list defining the boundary.
* The entry with K = 1 is a member of the boundary list. To scan
* the boundary, start at K = 1 and follow the forward or backward
* pointers.
*
* MB [in]
*
* NCONT [in] = 0 or 1 to request either C0 or C1 continuity.
*
* Q(1:2) [in] The (x,y) coordinates of the point for which this
* subr will attempt to find an enclosing triangle and then
* do interpolation.
*
* ZOUT [out] Interpolated value computed by this subroutine.
*
* WANTDZ [in] =.TRUE. means compute DZOUT() as well as ZOUT.
* =.FALSE. means compute only ZOUT and not DZOUT().
*
* DZOUT(1:2) [out] First partial derivs w.r.t. x and y of the
* interpolated surface at the interpolation point.
*
* MODEFI [out]
* =0 Means OK. Q is interior or almost so.
* =1 Q is exterior by more than the built-in tolerance.
* An extrapolated value will be returned.
* =2 Bad. Subroutine is cycling in the search. This should
* never happen.
*
* SAVWRK(1:28) [inout] This array is used as work space and to save
* quantities that may be reusable on a subsequent call.
* Let us say the settings of X(), Y(), Z(), DZ(), TRIANG(),
* and NT define a "grid and data". On the first call to this
* subroutine for use of a particular "grid and data" the user must
* set SAVWRK(1) = 0. This subroutine will generally reset
* SAVWRK(1) to a nonzero value on return, indicating there are
* saved quantities in SAVWRK() relating to this "grid and data".
* If the user will be switching back and forth between different
* "grid and data" specifications, some saving of computation time
* can be achieved by using distinct SAVWRK() arrays associated with
* distinct "grid and data" specifications. It is the user's
* responsibility to set SAVWRK(1) = 0. on any call for which the
* current SAVWRK() array is not the one that was used on the most
* recent previous call involving the current "grid and data".
*
* ------------------------------------------------------------------
* How the array SAVWRK() is used internally.
*
* The 1st location of SAVWRK() is regarded as containing a value
* of INDSGN, and the remaining 27 locations contain S(1:3, 1:9).
*
* We start by setting INDSGN = SAVWRK(1)
* and INDTRI = abs(INDSGN)
* We start the current lookup at INDTRI. If INDTRI is positive this */
/* causes the current lookup to start in the triangle at which the
* previous lookup for this "grid and data" ended, which in some
* cases significantly reduces the search time.
*
* If INDSGN > 0, and NCONT = 2, and the value of INDTRI does not
* change during the lookup, we assume there is useful saved info in
* S(1:3,7:9). Otherwise we compute this info as needed.
*
* Contents of S(,)
*
* (s(1:3, 1:9) [inout] Columns 1-6 contain data specifying the
* interpolation problem. Col 1 contains unnormalized
* barycentric coordinates of the interpolation point.
* The other cols contain data depending only on the
* triangle and its vertex data and not on the
* interpolation point. Cols 7-9 are used to save computed
* quantities depending only on the triangle and the
* partial deriv values from one call to the next.
* s( ,1) = Unnormalized barycentric coords of interp point.
* s( ,2) = u = x coord of edge vector
* s( ,3) = v = y coord of edge vector
* s( ,4) = z = Fcn value at vertex
* s( ,5) = zx = Partial deriv of z w.r.t. x
* s( ,6) = zy = Partial deriv of z w.r.t. y
* s( ,7) = htilda = Scaled p.d. tangential to edge at left end.
* s( ,8) = ktilda = Scaled p.d. tangential to edge at right end.
* s( ,9) = lfac = Factor involving lengths of edges.
*
* Let the indices of the vertices of this triangle, in counter-
* clockwise order be denoted by P(1), P(2), and P(3).
* For descriptive convenience we regard the subscript of P() and the
* 1st subsubscript of S(,) as always being reduced modulo 3 to 1, 2,
* or 3. Also for convenience we shall write P(i) to mean the vertex
* indexed by P(i).
*
* The unnormalized barycentric coordinate that is zero along the edge
* from P(i) to P(i+1) and has a positive value at P(i+2) is stored in
* S(i,1). The (x,y) coordinates of the vector from P(i) to P(i+1) are
* stored in (S(i,2), S(i,3)).
*
* The function value Z at P(i+2) is stored in S(i,4) for i = 1, 2, & 3.
* ZX, the partial deriv of Z w.r.t. X at P(i+2) is stored in S(i,5)
* for i = 1, 2, & 3.
* ZY, the partial deriv of Z w.r.t. Y at P(i+2) is stored in S(i,6)
* for i = 1, 2, & 3.
* ------------------------------------------------------------------
* Internal Variables
*
* NEWTRI [in] = .false. means that saved quantities are
* present in s(1:3, 7:9).9).
* If newtri=.true. the subr will compute s(1:3, 7:9).
*
* INDTRI On entry to STGFND, indtri designates the triangle at which
* the search will begin. On return from STGFND, indtri is the
* index of the last triangle tested.
* This is the triangle containing q if MODEFND=0.
* If MODEFND = 1, 2, or 3, this is a boundary triangle and q is
* outside a boundary edge of this triangle.
*
* TRI(1:8) Integer array of pointers defining one triangle. It is
* never assumed to contain any saved values on entry to
* STGFND. On return from STGFND, TRI(1:7) defines the triangle
* whose index is the returned value of INDTRI.
*
* The indices of the vertices of this triangle, in counter-
* clockwise order, are P(1), P(2), and P(3), which may be obtained
* as P(i) = TRI(3+i) for i = 1, 2, and 3. TRI(7) contains the same
* value as TRI(4). TRI(8) will be set in this subr to equal TRI(5).
* The triangle adjacent to this triangle across the edge from P(i)
* to P(i+1) is indexed by TRI(i) for i = 1, 2, and 3. If there is
* no adjacent triangle across this edge then TRI(i) = 0.
*
*
* ------------------------------------------------------------------
*--S replaces "?": ?TGFI, ?TGFND, ?TGC0, ?TGC1, ?TGEXT
* ------------------------------------------------------------------ */
/* ------------------------------------------------------------------ */
indsgn = Savwrk[1];
indtri = labs( indsgn );
stgfnd( x, y, triang, nt, q, &indtri, tri, (float(*)[3])&Savwrk[2],
&modefnd );
if (modefnd == 0)
{
/* The point q(1:2) was found ok.
* tri(1:7) and s(1:3, 1:3) contain info set by DTRFND.
* Next move info into s(1:3, 4:6) and call DTRC0 or DTRC1 for
* interpolation.
* Note that savwrk(-2+i + 3*j) contains s(i,j).
*
* S(1,4) = Z(TRI(6))
* S(2,4) = Z(TRI(4))
* S(3,4) = Z(TRI(5))
* */
Savwrk[1 + 10] = Z[Tri[6]];
Savwrk[1 + 11] = Z[Tri[4]];
Savwrk[1 + 12] = Z[Tri[5]];
if (ncont == 0)
{
stgc0( (float(*)[3])&Savwrk[2], zout, wantdz, dzout );
indsgn = -indtri;
}
else
{
/* S(1,5) = DZ(1,TRI(6))
* S(2,5) = DZ(1,TRI(4))
* S(3,5) = DZ(1,TRI(5))
*
* S(1,6) = DZ(2,TRI(6))
* S(2,6) = DZ(2,TRI(4))
* S(3,6) = DZ(2,TRI(5)) */
Savwrk[1 + 13] = dz[Tri[6] - 1][0];
Savwrk[1 + 14] = dz[Tri[4] - 1][0];
Savwrk[1 + 15] = dz[Tri[5] - 1][0];
Savwrk[1 + 16] = dz[Tri[6] - 1][1];
Savwrk[1 + 17] = dz[Tri[4] - 1][1];
Savwrk[1 + 18] = dz[Tri[5] - 1][1];
newtri = indtri != indsgn;
stgc1( newtri, (float(*)[3])&Savwrk[2], zout, wantdz,
dzout );
indsgn = indtri;
}
*modefi = 0;
}
else if (modefnd > 0)
{
/* Here MODEFND = 1, 2, or 3.
* Point q is outside the convex hull of the data.
* Use extrapolation.
* */
stgext( x, y, z, dz, triang, bdry, mb, ncont, q, indtri, modefnd,
zout, wantdz, dzout );
indsgn = -indtri;
*modefi = 1;
}
else
{
/* Here MODEFND = -1.
* This means cycling has happened in STGFND. this
* should not happen. Error message will have been printed
* from STGFND.
* */
*zout = ZERO;
if (wantdz)
{
Dzout[1] = ZERO;
Dzout[2] = ZERO;
}
indsgn = -indtri;
*modefi = 2;
}
Savwrk[1] = indsgn;
return;
} /* end of function */