geompack_cutfac.cxx

hl2 source code. Do not use it illegal.

/* cutfac.f -- translated by f2c (version 19990311).
*/


#include <ivp_physics.hxx>
#include <ivp_betterdebugmanager.hxx>
#include <geompack.hxx>


void IVP_Geompack::cutfac_(
			   int		*p,
			   double	*nrmlc,
			   double	*dtol,
			   int		*nedgc,
			   int		*pedge,
			   int		*nce,
			   int		*cedge,
			   double	*cdang,
			   long int	*rflag
			  ) {

    double *vcl              = this->g_vcl;
    int    *facesdata        = this->g_facesdata;
    double *normals          = this->g_normals;
    int    *faceverticeslist = this->g_faceverticeslist;
    double *edge_angles      = this->g_edge_angles;

    /* System generated locals */
    int i__1, i__2;
    double d__1, d__2, d__3, d__4;


    /* Local variables */
    double angr, cmax;
    int imin, kmax;
    double dist, nmax, dotp, tmin, ntol;
    int a, e, f, i__, j, k, l, n;
    long int eflag;
    double s, t;
    int w;
    double iamin;
    int ccwfl;
    double dsave[3];
    int isave;
    double dirsq;
    int estop, inout, estrt;
    double dir1sq;
    int ca, cb;
    double de[3];
    int ee, la, lb, fl;
    double cp[3];
    int fr, sf, lv, lw, sp;
    double intang;
    double pi2;
    int lw1;
    double dee[3], ang;
    long int dof;
    double dir[3];
    int nev, sgn;
    double rhs[3], dir1[3];


/*     Written and copyright by: */
/*        Barry Joe, Dept. of Computing Science, Univ. of Alberta */
/*        Edmonton, Alberta, Canada  T6G 2H1 */
/*        Phone: [phone]      Email: [e-mail] */

/*     Purpose: Trace out cut face in polyhedron P given a starting edge, */
/*        e.g. a reflex edge or an edge in the interior of a face. */
/*        Accept cut face if it creates no small angles, it is an outer */
/*        boundary and there are no holes (inner polygons) in it. It is */
/*        assumed starting edge does not lie on a double-occurring face. */

/*     Input parameters: */
/*        P - polyhedron index */
/*        NRMLC(1:4) - unit normal vector of cut plane plus right hand */
/*              side constant term of plane equation */
/*        ANGACC - min acceptable dihedral angle in radians produced by */
/*              a cut face */
/*        DTOL - absolute tolerance to determine if point is on cut plane */
/*        ADDR_OF_N_ORIGINAL_VERTICES - number of vertex coordinates */
/*        SIZE_VCL - maximum size available for VCL array */
/*        VCL(1:3,1:ADDR_OF_N_ORIGINAL_VERTICES) - vertex coordinate list */
/*        FACESDATA(1:3,1:*) - face pointer list */
/*        NRML(1:3,1:*) - unit normal vectors for faces */
/*        FACEVERTICESLIST(1:6,1:*) - face vertex list */
/*        EDGE_ANGLES(1:*) - edge angles */
/*        NEDGC - number of edges which intersect cut plane */
/*        PEDGE(1:3,1:NEDGC) - edges of polyh P that intersect cut plane */
/*              excluding starting edge if it is an edge of P or the edge */
/*              containing CEDGE(1,1) if CEDGE(1,1) > ADDR_OF_N_ORIGINAL_VERTICES; */
/*              PEDGE(1,I), PEDGE(2,I) are indices of FACEVERTICESLIST; if PEDGE(1,I) */
/*              = A and B = FACEVERTICESLIST(SUCC,A) then PEDGE(3,I) = 10*CA+CB where */
/*              CA = 1, 2, or 3 if vertex A in negative half-space, on */
/*              cut plane, or in positive half-space determined by cut */
/*              plane, and similarly for CB */
/*        CEDGE(1:2,0:1) - CEDGE(1,0) = LV, CEDGE(1,1) = LU where LU, LV */
/*              are indices of VCL and are vertices of starting edge; if */
/*              called by RESEDG, LU < LV <= ADDR_OF_N_ORIGINAL_VERTICES are on reflex edge and */
/*              CEDGE(2,1) = U is index of FACEVERTICESLIST indicating reflex edge; */
/*              LV may be ADDR_OF_N_ORIGINAL_VERTICES+1 and LU may be ADDR_OF_N_ORIGINAL_VERTICES+1 or ADDR_OF_N_ORIGINAL_VERTICES+2 to indicate */
/*              a point on interior of an edge; CEDGE(2,1) is specified */
/*              as described for output below; if LV > ADDR_OF_N_ORIGINAL_VERTICES, CEDGE(2,0) */
/*              takes on the value ABS(CEDGE(2,NCE)) described for output */
/*              below, else CEDGE(2,0) is not used */
/*        CDANG(1) - dihedral angle at starting edge determined by cut */
/*              plane in positive half-space */

/*     Updated parameters: */
/*        VCL - some temporary or permanent entries may be added at end */
/*              of array */
/*        NEDGC - may be decreased to indicate number of unchecked edges */
/*        PEDGE(1:3,1:NEDGC) - columns may be permuted (NEDGC=input val) */

/*     Output parameters: */
/*        NCE - number of edges in cut polygon; it is assumed there is */
/*              enough space in the following two arrays */
/*        CEDGE(1:2,0:NCE) - CEDGE(1,I) is an index of VCL, indices > ADDR_OF_N_ORIGINAL_VERTICES */
/*              are new points; CEDGE(2,I) = J indicates that edge of cut */
/*              face ending at CEDGE(1,I) is edge from J to FACEVERTICESLIST(SUCC,J) */
/*              if J > 0; else if J < 0 then edge of cut face ending at */
/*              CEDGE(1,I) is a new edge and CEDGE(1,I) lies on edge from */
/*              -J to FACEVERTICESLIST(SUC,-J) and new edge lies in face FACEVERTICESLIST(FACN,-J); */
/*              CEDGE(2,I) always refers to an edge in the subpolyhedron */
/*              in negative half-space; CEDGE(1,NCE) = CEDGE(1,0) */
/*        CDANG(1:NCE) - dihedral angles created by edges of cut polygon */
/*              in pos. half-space; neg. sign for angle I indicates that */
/*              face containing edge I is oriented CW in polyhedron P */
/*        RFLAG - .TRUE. iff reflex or starting edge is resolved */

/*     Routines called: */
/*        PTPOLG */

/*     Abnormal return: */
/*        IERR is set to 14, 325, 326, or 328 */




    /* Parameter adjustments */
    --nrmlc;
    vcl -= 4;
    facesdata -= 4;
    normals -= 4;
    faceverticeslist -= 7;
    --edge_angles;
    --cdang;

    /* Function Body */
    *rflag = 0;
    pi2 = IVP_PI + IVP_PI;
/* Computing MAX */
    i__1 = max(this->n_original_vertices, cedge[1]);
    n = max(i__1,cedge[3]);
    *nce = 1;
    w = cedge[4];
    lv = cedge[1];
    lw = lv;
    lw1 = cedge[3];
    if (lw1 <= this->n_original_vertices ) {
	nev = 1;
	this->g_ev[nev - 1] = lw1;
    } else {
	nev = 0;
    }
    if (w > 0) {
	fl = faceverticeslist[w * 6 + 2];
	if (lw > lw1) {
	    fr = faceverticeslist[faceverticeslist[w * 6 + 6] * 6 + 2];
	} else {
	    fr = faceverticeslist[faceverticeslist[w * 6 + 5] * 6 + 2];
	}
    } else {
	fl = faceverticeslist[-w * 6 + 2];
	fr = fl;
    }
    dir[0] = vcl[lw1 * 3 + 1] - vcl[lw * 3 + 1];
    dir[1] = vcl[lw1 * 3 + 2] - vcl[lw * 3 + 2];
    dir[2] = vcl[lw1 * 3 + 3] - vcl[lw * 3 + 3];
    kmax = 1;
    if (abs(nrmlc[2]) > abs(nrmlc[1])) {
	kmax = 2;
    }
    if (abs(nrmlc[3]) > (d__1 = nrmlc[kmax], abs(d__1))) {
	kmax = 3;
    }
    if ((i__1 = facesdata[fl * 3 + 2], abs(i__1)) == *p) {
	ccwfl = facesdata[fl * 3 + 2];
    } else {
	ccwfl = facesdata[fl * 3 + 3];
    }
    if (ccwfl < 0) {
	cdang[1] = -cdang[1];
    }

/*     LW, LW1, FL, FR, DIR are updated before each iteration of loop. */
/*     CCWFL = P (-P) if FL is CCW (CW) according to SUCC traversal. */

L10:
    if (lw1 > this->n_original_vertices ) {

/*           LW1 is new vertex on interior of edge E. FL is used for */
/*           previous and next faces containing edges of cut polygon. */

	e = -cedge[(*nce << 1) + 2];
	la = faceverticeslist[e * 6 + 1];
	lb = faceverticeslist[faceverticeslist[e * 6 + 3] * 6 + 1];
	if ((lb - la) * ccwfl > 0) {
	    ee = faceverticeslist[e * 6 + 6];
	} else {
	    ee = faceverticeslist[e * 6 + 5];
	}
	fl = faceverticeslist[ee * 6 + 2];
	dof = (i__1 = facesdata[fl * 3 + 2], abs(i__1)) == (i__2 = facesdata[fl * 3 + 
		3], abs(i__2));
	if (dof) {
	    l = faceverticeslist[ee * 6 + 1];
	    if (l == la) {
		ccwfl = -ccwfl;
	    }
	} else if ((i__1 = facesdata[fl * 3 + 2], abs(i__1)) == *p) {
	    ccwfl = facesdata[fl * 3 + 2];
	} else {
	    ccwfl = facesdata[fl * 3 + 3];
	}
	dir[0] = nrmlc[2] * normals[fl * 3 + 3] - nrmlc[3] * normals[fl * 3 + 2];
	dir[1] = nrmlc[3] * normals[fl * 3 + 1] - nrmlc[1] * normals[fl * 3 + 3];
	dir[2] = nrmlc[1] * normals[fl * 3 + 2] - nrmlc[2] * normals[fl * 3 + 1];
	if ((i__1 = facesdata[fl * 3 + 2], abs(i__1)) != *p || dof && ccwfl < 0) {
	    dir[0] = -dir[0];
	    dir[1] = -dir[1];
	    dir[2] = -dir[2];
	}
	goto L70;
    } else {

/*           LW1 is existing vertex of polyhedron P. FL (FL and FR) is */
/*           previous face if edge ending at LW1 is new (already exists). */
/*           In former case, -CEDGE(2,NCE) is the edge of FL incident on */
/*           LW1 which will lie only in subpolyhedron PL. In latter case, */
/*           CEDGE(2,NCE) is an edge of FL. Cycle thru edges, faces CCW */
/*           (from outside) between edges ESTRT, ESTOP. */
/*           If LW1 lies on a doubly-occurring face, there are 2 cycles */
/*           around LW1 and the correct one is chosen based on CCWFL. */

	iamin = pi2;
	imin = 0;
/* Computing 2nd power */
	d__1 = dir[0];
/* Computing 2nd power */
	d__2 = dir[1];
/* Computing 2nd power */
	d__3 = dir[2];
	dirsq = d__1 * d__1 + d__2 * d__2 + d__3 * d__3;
	eflag = cedge[(*nce << 1) + 2] > 0;
	if (! eflag) {
	    estrt = -cedge[(*nce << 1) + 2];
	    sp = ccwfl;
	    if (ccwfl > 0) {
		estop = faceverticeslist[estrt * 6 + 3];
	    } else {
		estop = faceverticeslist[estrt * 6 + 4];
	    }
	} else {
	    w = cedge[(*nce << 1) + 2];
	    la = faceverticeslist[w * 6 + 1];
	    if (ccwfl > 0) {
		estrt = faceverticeslist[w * 6 + 4];
	    } else {
		estrt = faceverticeslist[w * 6 + 3];
	    }
	    if (la == lw) {
		l = lw1 - lw;
	    } else {
		l = lw - lw1;
	    }
	    if (l * ccwfl > 0) {
		w = faceverticeslist[w * 6 + 6];
	    } else {
		w = faceverticeslist[w * 6 + 5];
	    }
	    if ((i__1 = facesdata[fr * 3 + 2], abs(i__1)) == (i__2 = facesdata[fr * 3 
		    + 3], abs(i__2))) {
		lb = faceverticeslist[w * 6 + 1];
		if (la == lb) {
		    sp = -ccwfl;
		} else {
		    sp = ccwfl;
		}
	    } else if ((i__1 = facesdata[fr * 3 + 2], abs(i__1)) == *p) {
		sp = facesdata[fr * 3 + 2];
	    } else {
		sp = facesdata[fr * 3 + 3];
	    }
	    if (sp > 0) {
		estop = faceverticeslist[w * 6 + 3];
	    } else {
		estop = faceverticeslist[w * 6 + 4];
	    }
	}
	la = faceverticeslist[estop * 6 + 1];
	lb = faceverticeslist[faceverticeslist[estop * 6 + 3] * 6 + 1];
	if ((lb - la) * sp > 0) {
	    estop = faceverticeslist[estop * 6 + 6];
	} else {
	    estop = faceverticeslist[estop * 6 + 5];
	}
	e = estrt;
	sf = ccwfl;
L20:
	if (eflag || e != estrt && e != estop) {
	    if (faceverticeslist[e * 6 + 1] == lw1) {
		l = faceverticeslist[faceverticeslist[e * 6 + 3] * 6 + 1];
	    } else {
		l = faceverticeslist[e * 6 + 1];
	    }
	    dist = nrmlc[1] * vcl[l * 3 + 1] + nrmlc[2] * vcl[l * 3 + 2] + 
		    nrmlc[3] * vcl[l * 3 + 3] - nrmlc[4];
	    if (abs(dist) <= *dtol) {
		dir1[0] = vcl[l * 3 + 1] - vcl[lw1 * 3 + 1];
		dir1[1] = vcl[l * 3 + 2] - vcl[lw1 * 3 + 2];
		dir1[2] = vcl[l * 3 + 3] - vcl[lw1 * 3 + 3];
/* Computing 2nd power */
		d__1 = dir1[0];
/* Computing 2nd power */
		d__2 = dir1[1];
/* Computing 2nd power */
		d__3 = dir1[2];
		dir1sq = d__1 * d__1 + d__2 * d__2 + d__3 * d__3;
		dotp = -(dir[0] * dir1[0] + dir[1] * dir1[1] + dir[2] * dir1[
			2]) / sqrt(dirsq * dir1sq);
		if (abs(dotp) > 1. - this->tolerance) {
		    dotp = d_sign(1.0, dotp);
		}
		if (kmax == 1) {
		    cp[0] = dir[1] * dir1[2] - dir[2] * dir1[1];
		} else if (kmax == 2) {
		    cp[1] = dir[2] * dir1[0] - dir[0] * dir1[2];
		} else {
		    cp[2] = dir[0] * dir1[1] - dir[1] * dir1[0];
		}
/* Computing MAX */
		d__2 = abs(dir[0]), d__3 = abs(dir[1]), d__2 = max(d__2,d__3),
			 d__3 = abs(dir[2]), d__2 = max(d__2,d__3), d__3 = 
			abs(dir1[0]), d__2 = max(d__2,d__3), d__3 = abs(dir1[
			1]), d__2 = max(d__2,d__3), d__3 = abs(dir1[2]);
		if ((d__1 = cp[kmax - 1], abs(d__1)) <= this->tolerance * max(
			d__2,d__3)) {
		    intang = IVP_PI;
		} else if (cp[kmax - 1] * nrmlc[kmax] > 0.) {
		    intang = acos(dotp);
		} else {
		    intang = pi2 - acos(dotp);
		}
		if (intang < iamin) {
		    iamin = intang;