/*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 */