geompack_resedg.cxx

hl2 source code. Do not use it illegal.

    }
    for (i__ = 1; i__ <= 2; ++i__) {
	ccw = 1;
	if (i__ == 1) {
	    la = lu;
	    b = faceverticeslist[u_in * 6 + 4];
	    if (lu > lv) {
		ccw = -1;
	    }
	} else {
	    la = lv;
	    b = v;
	    if (lu < lv) {
		ccw = -1;
	    }
	}
	first = 1;
	f = fl;
L30:
	if (! first) {
	    if (faceverticeslist[b * 6 + 1] == la) {
		b = faceverticeslist[b * 6 + 4];
	    } else {
		b = faceverticeslist[b * 6 + 3];
	    }
	}
	c__ = faceverticeslist[b * 6 + 3];
	if ((i__1 = facesdata[f * 3 + 2], abs(i__1)) == (i__2 = facesdata[f * 3 + 3], 
		abs(i__2))) {
	    if (faceverticeslist[b * 6 + 1] == la) {
		sp = -ccw * p;
	    } else {
		sp = ccw * p;
	    }
	} else if ((i__1 = facesdata[f * 3 + 2], abs(i__1)) == p) {
	    sp = facesdata[f * 3 + 2];
	} else {
	    sp = facesdata[f * 3 + 3];
	}
	if ((faceverticeslist[c__ * 6 + 1] - faceverticeslist[b * 6 + 1]) * sp > 0) {
	    b = faceverticeslist[b * 6 + 6];
	} else {
	    b = faceverticeslist[b * 6 + 5];
	}
	f = faceverticeslist[b * 6 + 2];
	if (f == fr) {
	    goto L50;
	}
	if (first) {
	    first = 0;
	} else {
	    if (faceverticeslist[b * 6 + 1] == la) {
		c__ = faceverticeslist[b * 6 + 3];
	    } else {
		c__ = b;
	    }
	    lc = faceverticeslist[c__ * 6 + 1];
	    ep[0] = vcl[lc * 3 + 1] - vcl[la * 3 + 1];
	    ep[1] = vcl[lc * 3 + 2] - vcl[la * 3 + 2];
	    ep[2] = vcl[lc * 3 + 3] - vcl[la * 3 + 3];
	    nrmlc[0] = ep[1] * e[2] - ep[2] * e[1];
	    nrmlc[1] = ep[2] * e[0] - ep[0] * e[2];
	    nrmlc[2] = ep[0] * e[1] - ep[1] * e[0];
/* Computing 2nd power */
	    d__1 = nrmlc[0];
/* Computing 2nd power */
	    d__2 = nrmlc[1];
/* Computing 2nd power */
	    d__3 = nrmlc[2];
	    leng = sqrt(d__1 * d__1 + d__2 * d__2 + d__3 * d__3);
/* Computing MAX */
	    d__1 = abs(ep[0]), d__2 = abs(ep[1]), d__1 = max(d__1,d__2), d__2 
		    = abs(ep[2]);
	    if (leng <= this->tolerance * max(d__1,d__2)) {
		goto L30;
	    }
	    nrmlc[0] /= leng;
	    nrmlc[1] /= leng;
	    nrmlc[2] /= leng;
	    enc[0] = e[1] * nrmlc[2] - e[2] * nrmlc[1];
	    enc[1] = e[2] * nrmlc[0] - e[0] * nrmlc[2];
	    enc[2] = e[0] * nrmlc[1] - e[1] * nrmlc[0];
	    nrmlc[0] = enr[1] * enc[2] - enr[2] * enc[1];
	    nrmlc[1] = enr[2] * enc[0] - enr[0] * enc[2];
	    nrmlc[2] = enr[0] * enc[1] - enr[1] * enc[0];
/* Computing 2nd power */
	    d__1 = nrmlc[0];
/* Computing 2nd power */
	    d__2 = nrmlc[1];
/* Computing 2nd power */
	    d__3 = nrmlc[2];
	    leng = sqrt(d__1 * d__1 + d__2 * d__2 + d__3 * d__3);
	    if (leng <= this->tolerance) {
		goto L30;
	    }
	    t = nrmlc[k - 1];
	    nrmlc[0] = enc[1] * enl[2] - enc[2] * enl[1];
	    nrmlc[1] = enc[2] * enl[0] - enc[0] * enl[2];
	    nrmlc[2] = enc[0] * enl[1] - enc[1] * enl[0];
/* Computing 2nd power */
	    d__1 = nrmlc[0];
/* Computing 2nd power */
	    d__2 = nrmlc[1];
/* Computing 2nd power */
	    d__3 = nrmlc[2];
	    leng = sqrt(d__1 * d__1 + d__2 * d__2 + d__3 * d__3);
	    if (leng <= this->tolerance) {
		goto L30;
	    }
	    if (nrmlc[k - 1] * t < 0.) {
		goto L30;
	    }
	    dotp = enr[0] * enc[0] + enr[1] * enc[1] + enr[2] * enc[2];
	    if (e[k - 1] * t < 0.) {
		dotp = -dotp;
	    }
	    if (abs(dotp) > 1. - this->tolerance) {
		dotp = d_sign(1.0, dotp);
	    }
	    t = acos(dotp);
	    if (t < anglo || t > anghi) {
		goto L30;
	    }
	    i__1 = nang;
	    for (j = 1; j <= i__1; ++j) {
		if ((d__1 = ang[j - 1] - t, abs(d__1)) <= this->tolerance) {
		    goto L30;
		}
/* L40: */
	    }
	    ++nang;
	    ang[nang - 1] = t;
	    if (nang >= 5) {
		goto L80;
	    }
	}
	goto L30;
L50:
	;
    }
    t = edge_angles[u_in] * .5;
    i__1 = nang;
    for (j = 1; j <= i__1; ++j) {
	if ((d__1 = ang[j - 1] - t, abs(d__1)) <= this->tolerance) {
	    goto L70;
	}
/* L60: */
    }
    ++nang;
    ang[nang - 1] = t;
L70:
    if (nang <= 3 && edge_angles[u_in] <= IVP_PI * 1.33) {
	nang += 2;
	ang[nang - 2] = edge_angles[u_in] * .25;
	ang[nang - 1] = edge_angles[u_in] * .75;
	if (IVP_PI - ang[nang - 1] <= .1) {
	    ang[nang - 2] = edge_angles[u_in] * .375;
	    ang[nang - 1] = edge_angles[u_in] * .625;
	}
	if (ang[nang - 2] < anglo || ang[nang - 2] > anghi) {
	    ang[nang - 2] = ang[nang - 1];
	    --nang;
	}
	if (ang[nang - 1] < anglo || ang[nang - 1] > anghi) {
	    --nang;
	}
    } else if (nang < 5 && edge_angles[u_in] <= IVP_PI * 1.4) {
	t = edge_angles[u_in] * 0.375;
	if (t >= anglo && t <= anghi) {
	    ++nang;
	    ang[nang - 1] = t;
	}
    }
L80:
    IVP_IF(1) {
	IVP_IFDEBUG(IVP_DM_GEOMPACK_LEVEL3) {
	    ivp_debugmanager.dprint(IVP_DM_GEOMPACK_LEVEL3, "Candidate angles (dihedral angle = %f ", edge_angles[u_in] * 180.0 / IVP_PI);
	    for (i__ = 1; i__ <= nang; ++i__) {
		ivp_debugmanager.dprint(IVP_DM_GEOMPACK_LEVEL3, "%f ", ang[i__ - 1] * 180.0 / IVP_PI);
	    }
	}
    }

/*     For each angle in ANG array, try to resolve reflex edge as */
/*     follows. Compute unit normal NRMLC and equation of cut plane; this */
/*     normal is outward wrt subpolyhedron containing face FL; NRMLC(4) */
/*     is right hand side constant term of plane equation. */
/*     Then determine which edges of polyhedron intersect cut plane; */
/*     reject plane if there is a vertex within distance MINDIS of plane */
/*     which does not lie on plane. */

    i__1 = nang;
    for (i__ = 1; i__ <= i__1; ++i__) {
	cn = cos(ang[i__ - 1]);
	if ((i__2 = facesdata[fr * 3 + 2], abs(i__2)) != p) {
	    cn = -cn;
	}
	ce = sin(ang[i__ - 1]);
	nrmlc[0] = cn * normals[fr * 3 + 1] + ce * enr[0];
	nrmlc[1] = cn * normals[fr * 3 + 2] + ce * enr[1];
	nrmlc[2] = cn * normals[fr * 3 + 3] + ce * enr[2];
	nrmlc[3] = nrmlc[0] * vcl[lu * 3 + 1] + nrmlc[1] * vcl[lu * 3 + 2] + 
		nrmlc[2] * vcl[lu * 3 + 3];
	nedgc = n_edges;
	j = 1;
L90:
	la = faceverticeslist[intworkarray[j * 3 + 1] * 6 + 1];
	da = nrmlc[0] * vcl[la * 3 + 1] + nrmlc[1] * vcl[la * 3 + 2] + nrmlc[
		2] * vcl[la * 3 + 3] - nrmlc[3];
	if (abs(da) <= dtol) {
	    ca = 0;
	} else if (abs(da) <= mindis) {
	    goto L110;
	} else if (da < 0.) {
	    ca = -1;
	} else {
	    ca = 1;
	}
	b = faceverticeslist[intworkarray[j * 3 + 1] * 6 + 3];
	lb = faceverticeslist[b * 6 + 1];
	db = nrmlc[0] * vcl[lb * 3 + 1] + nrmlc[1] * vcl[lb * 3 + 2] + nrmlc[
		2] * vcl[lb * 3 + 3] - nrmlc[3];
	if (abs(db) <= dtol) {
	    cb = 0;
	} else if (abs(db) <= mindis) {
	    goto L110;
	} else if (db < 0.) {
	    cb = -1;
	} else {
	    cb = 1;
	}
	if (ca * cb > 0) {
	    for (k = 1; k <= 2; ++k) {
		c__ = intworkarray[k + j * 3];
		intworkarray[k + j * 3] = intworkarray[k + nedgc * 3];
		intworkarray[k + nedgc * 3] = c__;
/* L100: */
	    }
	    --nedgc;
	} else {
	    intworkarray[j * 3 + 3] = (ca + 2) * 10 + (cb + 2);
	    ++j;
	}
	if (j <= nedgc) {
	    goto L90;
	}
	intworkarray[(n_edges + 1) * 3 + 1] = lvo;
	intworkarray[(n_edges + 1) * 3 + 3] = luo;
	intworkarray[(n_edges + 2) * 3 + 1] = u_in;
	doubleworkarray[1] = ang[i__ - 1];
	IVP_IF(1) {
	    IVP_IFDEBUG(IVP_DM_GEOMPACK_LEVEL3) {
		ivp_debugmanager.dprint(IVP_DM_GEOMPACK_LEVEL3, "Trying angle %d   NEDGC = %d\n", i__, nedgc);
	    }
	}

	// *******************************************************************
	cutfac_(&p,
		nrmlc,
		&dtol,
		&nedgc,
		intworkarray,
		&nce,
		&intworkarray[((n_edges+1)*3)],
		&doubleworkarray[1],
		rflag
		);

	vcl = this->g_vcl;
	vcl -= 4;

	if (this->ierr != 0) {
	    return;
	}
	if (*rflag) {

	    // ***************************************************************
	    insfac_(&p,
		    nrmlc,
		    &nce,
		    &intworkarray[(n_edges+1)*3],
		    &doubleworkarray[1]
		    );

	    polyhedronfirstfaceoffset = this->g_polyhedronfirstfaceoffset;
	    --polyhedronfirstfaceoffset;

	    polyhedronfaceindices = this->g_polyhedronfaceindices;
	    polyhedronfaceindices -= 3;

	    facesdata = this->g_facesdata;
	    facesdata -= 4;

	    normals = this->g_normals;
	    normals -= 4;

	    faceverticeslist = this->g_faceverticeslist;
	    faceverticeslist -= 7;

	    edge_angles = this->g_edge_angles;
	    --edge_angles;

	    return;
	}
L110:
	;
    }

    return;
}