/Users/craigcornelius/Projects/SPRING Mac Release 0.2/edge.cpp

Go to the documentation of this file.

/* $Id: edge.cpp,v 1.59 2006/01/06 19:28:11 craig Exp $ */
#include "edge.h"
#include "node.h"
#include "face.h"
#include "object.h"
#include "point3d.h"
#include "nodearray.h"
#include "edgearray.h"
#include "rasterfont.h"
#include <stdio.h>
#include <math.h>

#ifdef _WIN32
#include <windows.h>
#endif

#ifdef __APPLE__
#include <OpenGL/gl.h>
#else
#include <GL/gl.h>
#endif // __APPLE__

const char* Edge::rcsid = "@(#) $Id: edge.cpp,v 1.59 2006/01/06 19:28:11 craig Exp $ $Copyright: (c)2001 National Biocomputation Center, Stanford University $";


int Edge::debug = 0;

Edge::Edge()
{
    node[0] = node[1] = NULL;
    splitNode[0] = splitNode[1] = NULL;
    marker = -1.0;
    boundingsphere = NULL;
}

Edge::~Edge()
{
    node[0] = node[1] = NULL;
}

// this fills in values for an Edge from the edge data file
int Edge::init(NodeArray* nodearray_p, EdgeArray* edgearray_p_in, 
               char *initString)
{
    int n1, n2, index;
    double restLength_given, springConstant_given, dampConstant_given;
    
    // parse string
    if (sscanf(initString, "%d %d %d %lf %lf %lf", &index, &n1, &n2,
        &restLength_given, &springConstant_given, &dampConstant_given) !=6){ 
        printf("Wrong number of arguments in edge data line\n");
        return 0;
    }
    
    // call the main init routine
    int s = init(nodearray_p, edgearray_p_in, n1, n2, restLength_given, 
               springConstant_given, dampConstant_given);
    return(s);
}

// this fills in values for an Edge from the given data
int Edge::init(NodeArray* nodearray_p, EdgeArray* edgearray_p_in,
               int i0, int i1, double restLength_in, 
               double springConstant_in, double dampConstant_in)
{
    // check if node indices are in our node array
    int numnodes = nodearray_p->getNumNodes();
    if ((i0 < 0) || (i1 < 0) || (i0 > numnodes) || (i1 > numnodes)) {
        printf("Node numbers in edge data line are out of bounds\n");
        return(0);
    }
    else if (i0 == i1) {
        printf("Edge can't connect a node to itself\n");
        return(0);
    }
    
    // get node ptrs from the given indices
    Node* n0 = nodearray_p->getNode(i0);
    Node* n1 = nodearray_p->getNode(i1);
    
    // call the real init
    int s = init(nodearray_p, edgearray_p, n0, n1, restLength_in, 
               springConstant_in, dampConstant_in);
    
    // and return its value
    return(s);
}

// this fills in values for an Edge from the given data
int Edge::init(NodeArray* nodearray_p, EdgeArray* edgearray_p_in,
               Node* n0, Node* n1, double restLength_in, 
               double springConstant_in, double dampConstant_in)
{
    // set private vars
    edgearray_p = edgearray_p_in;
    node[0] = n0;
    node[1] = n1;
    
    // add our edge to the nodes
    node[0]->addEdgeLink(this);  // then update the edge arrays for these Nodes
    node[1]->addEdgeLink(this);
    
    // calculate the restlength of the spring.
    switch( (int) restLength_in ) {
    case -1: restLength = this->getLength(); break;
    case -2: restLength = this->getLength() * 0.95; break;
    case -3: restLength = this->getLength() * 1.07; break;
    default: restLength = restLength_in; break;
    }
    
    springConstant = springConstant_in;
    dampConstant = dampConstant_in;
    
    return (1);
}


/* this is called three times for each Face in the meshfile, once for each of
* the edges it contains */
// Note: Really only need the faces to get the normals to get the angle for
//                      torsion processing (which we don't use anymore)
void Edge::addFaceLink(Face *f)
{
    if (debug)
    cout << "\nEdge " << getIndex() << " ::addFaceLink(" << f->getIndex() << ")\n";
  face.Append(f);
    
#ifdef TORSION
    // handle torsion/initial angles
    if (face.NumElements()) {
        /* for now, calculate the initial angle between the faces here */
        // but we need normals to do that:
        face[0]->computenormal();
        face[1]->computenormal();
        init_angle = getAngle();
    }
#endif
}

void Edge::addTetraLink(Tetra *t)
{
  if (debug)
    cerr << "\nEdge " << getIndex() << " ::addTetraLink(" << t->getIndex() << ")\n";
  tetra.Append(t);
}

void Edge::deleteFaceLink(Face* f)
{
  if (debug)
    cout << "\nEdge " << getIndex() << " ::deleteFaceLink(" << f->getIndex() << ")\n";

  face.DeleteItem(f);
        if (debug) 
          cout << "ended up with " << face.NumElements() << " face left\n";
}

void Edge::deleteTetraLink(Tetra *t)
{
  if (debug)
    cout << "\nEdge " << getIndex() << " ::deleteTetraLink(" << t->getIndex() << ")\n";

  tetra.DeleteItem(t);
        if (debug) 
          cout << "ended up with " << tetra.NumElements() << " tetra left\n";
}

void Edge::replaceNodeLink(Node* old_node, Node* new_node)
{
    for (int i=0; i<2; i++) {
        if (node[i] == old_node) node[i] = new_node;
    }
}

int Edge::getIndex()
{ return(edgearray_p->getIndex(this)); }

/* this calculates the angle between this Edge's 2 adjacent faces.
* Normals must be consistent (hopefully outward facing), the angle is then
* starting at one face, in the direction of the normal, to the other face*/
double Edge::getAngle(void)
{
    Point3D norm1, norm2, vec;
    double dot, theta, dot2, retangle;
    if (face.NumElements()!= 2) { // to do: caller should check this?
        printf("called getAngle with face.NumElements() != 2, bad \n");
        return (-1);
    }
    // get the normals of the adjacent faces (they are unit length)
    norm1 = face[0]->getNormal();
    norm2 = face[1]->getNormal();
    dot = norm1.Dot(norm2);
    theta = (acos(dot))*(180.0/PI);
    /* theta is now the angle between the normals.  actual angle between the
    * faces is 180 +- theta.  decide by taking a 3rd vector, starting at this
    * edge, and going into face1, and dotting with face0's normal, norm1 */
    Node *n = face[1]->getThirdNode(node[0], node[1]);
    // find this vector in face1 by going from node to n
    // don't care about magnitude, just need sign of dot product
    vec = n->p - node[0]->p;
    dot2 = norm1.Dot(vec);
    if (dot2 < 0.0)
        retangle = 180.0 + theta;
    else
        retangle = 180.0 - theta;
    
    if (debug)
        printf("edge %d normal angle %f face angle %f \n", 
        getIndex(), theta, retangle);
    return (retangle);
}

// Spring force with FULL damping
void Edge::ComputeSpringForceFullDamp(void)
{       
        // Uses endpoints of the edge and spring characteristics to find
        // the force between the nodes.
        // Result computed is force from first to second node of the edge.
        
        // This is directional.  Computes the force from the first to the 2nd node.
        
        /* calculate vector forces as we should:
        * magnitude f = -k*delta_position
        * direction f = spring_unit = diff/||diff|| */
        Point3D diff = node[1]->p - node[0]->p;
        
        // Any way to avoid the square root computation?
        double normSquared = diff.Squared();
        
        // ??? Can we remove this test - will the two points ever be at the same location?
        if (normSquared > 0.0) {
                Point3D vdiff = node[1]->v - node[0]->v;
                double velocityDotSpring = vdiff.Dot(diff);
                
                double springLength = sqrt(normSquared);
                double deltaPosition = springLength - restLength;
                
                double totalDampMagnitude = 
                        (dampConstant * velocityDotSpring / normSquared);
                
                force = -((springConstant * deltaPosition / springLength) + totalDampMagnitude) * diff;
        }
        else {
                force = Point3D(0.0, 0.0, 0.0);
        }
}

// Spring force with REL damping
void Edge::ComputeSpringForceRelDamp(void)
{
        // Uses endpoints of the edge and spring characteristics to find
        // the force between the nodes.
        // Result computed is force from first to second node of the edge.
        
        // This is directional.  Computes the force from the first to the 2nd node.
        
        /* calculate vector forces as we should:
        * magnitude f = -k*delta_position
        * direction f = spring_unit = diff/||diff|| */
        Point3D diff = node[1]->p - node[0]->p;
        
        // Any way to avoid the square root computation?
        double normSquared = diff.Squared();
        double springLength = sqrt(normSquared);
        double deltaPosition = springLength - restLength;

        // REL
        /* damping based on relative velocity for every Edge:
        * vxdiff:relative_velocity::fx:f
        * and f = -c*relative_velocity
        * so fx = -c*vxdiff */
        Point3D vdiff = node[1]->v - node[0]->v;

        force = -(springConstant * deltaPosition / springLength) * diff - dampConstant * vdiff;
}

// Computes the force on the vector from the first point to the second point,
// based on linear spring characteristics and damping.
// Note that the switch for dampingType slows this down.
// Perhaps make 3 specific versions to avoid the switch.
// Added 2-Nov-2005, cwc.

void Edge::ComputeSpringForce(int dampingType)
{
        // Uses endpoints of the edge and spring characteristics to find
        // the force between the nodes.
        // Result computed is force from first to second node of the edge.
        
        // Compute the force from the first to the 2nd node.
        
        /* calculate vector forces as we should:
        * magnitude f = -k*delta_position
        * direction f = spring_unit = diff/||diff|| */
        Point3D diff = node[1]->p - node[0]->p;
        
        // Any way to avoid the square root computation?
        double normSquared = diff.Squared();
        double springLength = sqrt(normSquared);
        double deltaPosition = springLength - restLength;

        // Note that ABS case is computed per node...
        // And NONE doesn't use any damping.
        switch (dampingType) {
        case ABS:
        case NONE:
                force = -(springConstant * deltaPosition / springLength) * diff;
                break;
        
        case FULL:
                // FULL
                /* damping based on relative velocity for every Edge, *BUT*
                * only in the direction of the spring:
                * f = -c*(relative velocity component along spring)
                *   = -c*(relative_velocity dot spring_unit)*spring_unit
                * relative_velocity = vdiff
                * spring_unit = diff/||diff||
                * f = -c*(vdiff dot diff)*diff/||diff||^2 */
                
                if (normSquared > 0.0) {
                        Point3D vdiff = node[1]->v - node[0]->v;
                        double velocityDotSpring = vdiff.Dot(diff);             // The fractional velocity along position vector.
                        
                        force = -(springConstant * deltaPosition / springLength + 
                                dampConstant * velocityDotSpring / normSquared) * diff;
                        
                }
                break;
                
        case REL:
                // REL
                /* damping based on relative velocity for every Edge:
                * vxdiff:relative_velocity::fx:f
                * and f = -c*relative_velocity
                * so fx = -c*vxdiff */
                Point3D vdiff = node[1]->v - node[0]->v;
                force = -(springConstant * deltaPosition / springLength) * diff;
                force -= dampConstant*vdiff;
                break;
                

        }
        
}

void Edge::ComputeSpringForceNonlinear(int dampingType)
{
        // Uses endpoints of the edge and spring characteristics to find
        // the force between the nodes.
        // Result computed is force from first to second node of the edge.
        
        // Compute the force from the first to the 2nd node.
        
        /* calculate vector forces as we should:
        * magnitude f = -k*delta_position
        * direction f = spring_unit = diff/||diff|| */
        Point3D diff = node[1]->p - node[0]->p;
        
        // Any way to avoid the square root computation?
        double normSquared = diff.Squared();
        double springLength = sqrt(normSquared);
        double deltaPosition = springLength - restLength;

        if (normSquared <= 0.0) {
                // There's a chance that the nodes are very close.
                force.x = force.y = force.z = 0.0;
                return;
        }

        double nonlinear_factor = 1.0;
        {
                double ratio = springLength/restLength;
                if (ratio < 1.0 && ratio != 0.0)
                        // Force increases for compression too!
                        ratio = 1/ratio;
                if (ratio < 2)
                        ; // between 1/2 and 2 times normal length do nothing
                else if (ratio < 4)
                        nonlinear_factor = 1.2;
                else if (ratio < 8)
                        nonlinear_factor = 1.4;
                else if (ratio >= 8)
                        nonlinear_factor = 1.6;
        }

        // This is the force computation, based on the type of damping.

        // Note that ABS case is computed per node...
        switch (dampingType) {

        case FULL:
                {
                        // FULL
                        /* damping based on relative velocity for every Edge, *BUT*
                        * only in the direction of the spring:
                        * f = -c*(relative velocity component along spring)
                        *   = -c*(relative_velocity dot spring_unit)*spring_unit
                        * relative_velocity = vdiff
                        * spring_unit = diff/||diff||
                        * f = -c*(vdiff dot diff)*diff/||diff||^2 */
                        
                        Point3D vdiff = node[1]->v - node[0]->v;
                        
                        double velocityDotSpring = vdiff.Dot(diff);
                        
                        double totalDampMagnitude = (dampConstant * velocityDotSpring / normSquared);
                        
                        // Combine these to reduce multiplications by the diff vector.
                        force = -(nonlinear_factor * springConstant * deltaPosition / springLength + totalDampMagnitude) * diff;
                }
                break;
                
        case REL:
                {
                        // REL
                        /* damping based on relative velocity for every Edge:
                        * vxdiff:relative_velocity::fx:f
                        * and f = -c*relative_velocity
                        * so fx = -c*vxdiff */
                        Point3D vdiff = node[1]->v - node[0]->v;
                        force = -(nonlinear_factor * springConstant * deltaPosition / springLength) * diff;
                        force -= dampConstant*vdiff;
                }
                break;

        case ABS:
        case NONE:
        default:
                {
                        // Just the simple force computation
                        force = -(nonlinear_factor * springConstant * deltaPosition / springLength) * diff;
                }
                break;

        }
        
}

/* function to calculate the torsion force between face0 and face1, and apply
* this force to the exterior nodes of those faces */
void Edge::addTorsion(void)
{
    if (face.NumElements() == 2) {
        // for now we'll compute the normals here
        face[0]->computenormal();
        face[1]->computenormal();
        // if this is positive, apply force in direction of normal, else opposite
        double anglediff = getAngle() - init_angle;
        Point3D force1 = (anglediff*TORSION_FORCE)*(face[0]->getNormal());
        Point3D force2 = (anglediff*TORSION_FORCE)*(face[1]->getNormal());
        Node *n1 = face[0]->getThirdNode(node[0], node[1]);
        Node *n2 = face[1]->getThirdNode(node[0], node[1]);
        n1->addForce(force1);
        n2->addForce(force2);
    }
}

/* This function returns the potential energy of this Edge (spring) 
* E = 0.5*k*x^2 where k is the spring constant, x is the displacement from
* resting position, E calculated by integrating F = k*x from 0 to x
*/
double Edge::getEnergy(void)
{
    double length = getLength();
    double displacement = fabs(length - restLength);
    double energy = 0.5*springConstant*displacement*displacement;
    return (energy);
}


Point3D Edge::UnitDirectionVector(void)
{
    Point3D vect = getNodeLink(1)->p - getNodeLink(0)->p;
    vect.Normalize();
    return vect;
}

// distance between edge and point, as well as param for near point on edge
// if no near pointer is sent in, it will create a dummy value
double Edge::EdgeDistanceToPoint(Point3D pt, double *near1)
{
    Point3D dirVect;
    Node *n0, *n1;
    double t;
    double dummy;
    if (near1 == NULL)
      near1 = &dummy;
    *near1 = 0.0;
    
    n0 = getNodeLink(0);
    n1 = getNodeLink(1);
    dirVect = UnitDirectionVector();
    
    /* project n0 to pt vector onto this edge, it actually lies on this edge
    * if projected distance is between 0 and this edge's length */
    t = dirVect.Dot(pt - n0->p);
   
    if( t <= 0 ) {// closest to node 0
      // value of near1 is already 0.0
      return( (pt - n0->p).Length() );
    }
    double length = getLength();
    if( t >= length ) {// closest to node 1
      *near1 = 1.0;
      return( (pt - n1->p).Length() );
    }
    //closest to some point on edge, figure out its parameter in (0,1)
    if (length != 0.0)
      *near1 = t/length; // t can range from 0 to edge length in this case
    return((pt - n0->p - t*dirVect).Length());
}

// returns distance between 2 edges, as well as near point on each edge
// if no near pointer is set in, it will create dummy values
double Edge::EdgeDistanceToEdge(Edge* e2, double *near1, double *near2)
{
    Edge *e1 = this;
    double dummy1;
    double dummy2;
    if (near1 == NULL) {
      near1 = &dummy1;
      near2 = &dummy2;
    }
    /* get lines containing e1 and e2, then use Goldman's Gem (vol 1, p 304)
    * to get the points of closest approach.  if these points are on the
    * edges, we're done, otherwise I think the closest distance must involve
    * one of the 4 endpoints */
    // L1(t) = p1 + v1*t, L2(s) = p2 + v2*s, v1 and v2 unit direction vectors
    Point3D p1 = e1->getNodeLink(0)->p;
    Point3D p2 = e2->getNodeLink(0)->p;
    Point3D v1 = e1->UnitDirectionVector();
    Point3D v2 = e2->UnitDirectionVector();
    Point3D p2minusp1 = p2 - p1;
    Point3D v1crossv2 = v1.Cross(v2);
    double denom = v1crossv2.Squared();
    
    if (denom != 0.0) {  //not parallel, there is a point of closest approach
        double t = (p2minusp1.Determinant(v2, v1crossv2))/denom;
        double s = (p2minusp1.Determinant(v1, v1crossv2))/denom;
        if (t > 0 && s > 0 && t < e1->getLength() && s < e2->getLength()) {
            // point of closest approach is on the segments, return length
            *near1 = t/(e1->getLength());
            *near2 = s/(e2->getLength());
            return (((p1 + v1*t) - (p2 + v2*s)).Length());
        }
    }
  // find the distance from each endpoint to the other segment, return min
  double nearA, nearB, nearC, nearD;
  double d1 = e1->EdgeDistanceToPoint(p2, &nearA);
  double d2 = e1->EdgeDistanceToPoint(e2->getNodeLink(1)->p, &nearB);
  double d3 = e2->EdgeDistanceToPoint(p1, &nearC);
  double d4 = e2->EdgeDistanceToPoint(e1->getNodeLink(1)->p, &nearD);

  double min12 = MIN(d1, d2);
  double min34 = MIN(d3, d4);
  double min = MIN(min12, min34);
  
  if (d1 == min) {
    *near1 = nearA;
    *near2 = 0.0;
  }
  else if (d2 == min) {
    *near1 = nearB;
    *near2 = 1.0;
  }
  else if (d3 == min) {
    *near1 = 0.0;
    *near2 = nearC;
  }
  else if (d4 == min) {
    *near1 = 1.0;
    *near2 = nearD;
  }
  else // nothing is the minimum???
    printf("Um, problem in EdgeDistanceToEdge\n");

  return (min);
}

// returns the distance between an edge and a face, as well as the
//      near points on each.
// ? does _not_ check for intersection, if they are intersecting, then the
//      motion limit has failed, and it's too late to deal with ?
double Edge::EdgeDistanceToFaceNoIntersect(Face* f, double *near1, 
                                           Point3D *nearpt)
{
  // 3 cases:
  // 1. This edge intersects face
  // 2. One of this edge's 2 endpoints closest to face
  // 3. This edge closest to one of face's 3 edges
  // We ignore 1., get minimum value from 2. and 3.
  // moreover, only care about 2. if endpoint projects inside face

  // early cutout?: if either endpoint is too far from face

  Point3D nearpt1, nearpt2;
  double nearA1, nearA2, nearB1, nearB2, nearC1, nearC2;
  Edge *fe0 = f->getEdgeLink(0);
  Edge *fe1 = f->getEdgeLink(1);
  Edge *fe2 = f->getEdgeLink(2);

  double d1 = f->FaceDistanceToPoint(getNodeLink(0)->p, 1, &nearpt1);
  double d2 = f->FaceDistanceToPoint(getNodeLink(1)->p, 1, &nearpt2);
  double d3 = EdgeDistanceToEdge(fe0, &nearA1, &nearA2);
  double d4 = EdgeDistanceToEdge(fe1, &nearB1, &nearB2);
  double d5 = EdgeDistanceToEdge(fe2, &nearC1, &nearC2);

  double min12 = MIN(d1, d2);
  double min34 = MIN(d3, d4);
  double min1234 = MIN(min12, min34);
  double min = MIN(min1234, d5);

  // near1 is param of edge near point, nearpt is face near point
  if (d1 == min) {
    *near1 = 0.0;
    *nearpt = nearpt1;
  }
  else if (d2 == min) {
    *near1 = 1.0;
    *nearpt = nearpt2;
  }
  else if (d3 == min) {
    *near1 = nearA1;
    *nearpt = fe0->getNodeLink(0)->p + nearA2*
      (fe0->getNodeLink(1)->p - fe0->getNodeLink(0)->p);
  }
  else if (d4 == min) {
    *near1 = nearB1;
    *nearpt = fe1->getNodeLink(0)->p + nearB2*
      (fe1->getNodeLink(1)->p - fe1->getNodeLink(0)->p);
  }
  else if (d5 == min) {
    *near1 = nearC1;
    *nearpt = fe2->getNodeLink(0)->p + nearC2*
      (fe2->getNodeLink(1)->p - fe2->getNodeLink(0)->p);
  }
  else // nothing is the minimum???
    printf("Um, problem in EdgeDistanceToFaceNoIntersect\n");

  return (min);
}

void Edge::draw()
{
    Node* n1 = getNodeLink(0);
    Node* n2 = getNodeLink(1);
    
    // needle edges are now handled in edgearray::DrawEdges
    
    if (marker >= 0.0) {
      // mark closest points of "colliding" edges in red
      glColor3d(1.0, 0.0, 0.0);
      glPointSize(5.0);
      Point3D pt = n1->p + marker*(n2->p - n1->p);
      glBegin(GL_POINTS);
      glVertex3d(pt.x, pt.y, pt.z);
      glEnd();
      // mark colliding edges as green
      glColor3d(0.0, 1.0, 0.0);
    }
    
    glLineWidth(2.0); // thread (if not textured) shows up better if wider
    glBegin(GL_LINES);
    n1->draw(0, 0, 0);
    n2->draw(0, 0, 0);
    glEnd();
    glLineWidth(1.0);
}

void Edge::drawover(double scale)
{
    // draw differently based on face.NumElements()
    switch (face.NumElements()) {
    case 0: glColor3d(1.0, 1.0, 0.0); break;    // yellow
    case 1: glColor3d(0.0, 1.0, 0.0); break;    // green
    case 2: glColor3d(0.0, 0.0, 0.0); break;    // black
    }
    
    Node* n1 = getNodeLink(0);
    Node* n2 = getNodeLink(1);
    glBegin(GL_LINES);
    n1->drawover(scale);
    n2->drawover(scale);
    glEnd();
}

void Edge::drawInitial()
{
    Node* n1 = getNodeLink(0);
    Node* n2 = getNodeLink(1);
    glColor3f(0.5, 0.5, 0.5);
    glBegin(GL_LINES);
    n1->drawinitial();
    n2->drawinitial();
    glEnd();
}

// make a label containing the face num
void Edge::drawlabel(float offset)
{
    extern RasterFont* font_p;
    char s[20]; sprintf(s, "%d", getIndex());
    Node* n1 = getNodeLink(0);
    Node* n2 = getNodeLink(1);
    Point3D center = (n1->p + n2->p)/2.0;
    
    Point3D avg_normal = (n1->getNormal() + n2->getNormal())/2.0;
    Point3D ll = center + offset*avg_normal;
    
    glBegin(GL_LINES);
    glVertex3d(center.x, center.y, center.z);
    glVertex3d(ll.x, ll.y, ll.z);
    glEnd();
    
    font_p->printString(s, ll.x,ll.y,ll.z);
}

// determines whether the plane ax+by+cz+d=0 crosses the edge [P1,P2]. If yes, 
// the function returns 1 (and the point of intersection is determined), else 0
int Edge::crossesPlane(double a, double b, double c, double d, 
                       Point3D *intersection_p)
{
    Point3D P1 = getNodeLink(0)->p;
    Point3D P2 = getNodeLink(1)->p;
    double alpha, t;
    
    alpha = a * (P1.x - P2.x) + b * (P1.y - P2.y) + c * (P1.z - P2.z);
    if(alpha == 0) return 0;
    
    t = -(a * P2.x + b * P2.y + c * P2.z + d) / alpha;
    if ((t < 0) || (t > 1)) return 0;
    
    if (intersection_p) {
        intersection_p->x = (P1.x - P2.x) * t + P2.x;
        intersection_p->y = (P1.y - P2.y) * t + P2.y;
        intersection_p->z = (P1.z - P2.z) * t + P2.z;
    }
    
    return 1;
}
int Edge::crossesPlane(Point3D p1, Point3D p2, Point3D p3, Point3D *intPt)
{
        int intersects = 0;
        double a,b,c,d;
                a = p1.y*(p2.z-p3.z) + p2.y*(p3.z - p1.z) + p3.y*(p1.z - p2.z);
                b = p1.z*(p2.x-p3.x) + p2.z*(p3.x - p1.x) + p3.z*(p1.x - p2.x);
                c = p1.x*(p2.y-p3.y) + p2.x*(p3.y - p1.y) + p3.x*(p1.y - p2.y);
                d = -(p1.x*(p2.y*p3.z - p2.z*p3.y) +
                                        p2.x*(p3.y*p1.z - p3.z*p1.y) + 
                                        p3.x*(p1.y*p2.z - p1.z*p2.y));
        intersects = crossesPlane(a,b,c,d,intPt);
        return intersects;
}

int Edge::SanityCheck(EdgeArray* real_ea_p)
{
    char s[1024]; s[0] = '\0';
    char tmps[1024]; tmps[0] = '\0';
    
    // look at my local stuff
    int my_index = -1;
    if (!edgearray_p) strcat(s,"edgearray_p is NULL!");
    if (edgearray_p != real_ea_p) strcat(s,"edgearray_p is incorrect!");
    else {
        my_index = getIndex();
        if (my_index < 0) { sprintf(tmps, "getIndex() = %d\n", my_index); strcat(s, tmps); }
    }
    
    // don't care about init_angle anymore (no more torsion)
    if (restLength <= 0) strcat(s, "restLength <= 0!");
    if (springConstant <= 0) strcat(s, "springConstant <= 0!");
    if (dampConstant <= 0) strcat(s, "dampConstant <= 0!");
    if (face.NumElements() < 0) strcat(s, "face.NumElements() < 0!");
    if (face.NumElements() == 0) strcat(s, "face.NumElements() = 0!");  // not an error in the future
    
    { // verify my node links
        for (int i=0; i<2; i++) {
            Node* n = getNodeLink(i);
            if (!n) {
                sprintf(tmps, "nodelink[%d] is NULL!", i); 
                strcat(s,tmps);
            }
            else if (!n->hasEdgeLink(this)) {
                sprintf(tmps, "nodelink[%d] (%d) does not have a link back to us!", i, n->getIndex() );
                strcat(s, tmps);
            }
        }
        
        // look for duplicates
        if (node[0] == node[1]) 
            strcat(s,"nodelinks point to the same node!\n");
    }
    
    { // verify my face links
        for (int i=0; i<face.NumElements(); i++) {
            Face* f = getFaceLink(i);
            if (!f) {
                sprintf(tmps, "facelink[%d] is NULL!", i);
                strcat(s, tmps);
            }
            else {
                if (!f->hasEdgeLink(this)) {
                    sprintf(tmps, "facelink[%d] (%d) does not have a link back to us!", i, f->getIndex() );
                    strcat(s, tmps);
                }
                
                // look for duplicate facelinks
                for (int j=0; j<i; j++) {
                    if (f == getFaceLink(j)) {
                        sprintf(tmps, "facelink[%d] is a duplicate of facelink[%d]!", i, j);
                        strcat(s, tmps);
                    }
                }
            }
        }
    }
    
    // if there's a problem, tell it
    if (strlen(s) > 1) {
        cerr << "Edge[" << my_index << "] is INSANE: " << s << endl;
        return(0);
    }
    
    //everything okay
    return(1);
}

ostream& operator<<(ostream& os, const Edge& e_in)
{
     Edge& e = (Edge&)e_in;

         if (e.edgearray_p)
        os << "Edge#" << ((Edge&)e).getIndex() << ": ";
    else os << "Edge#UNLINKED: ";
    
    // print out the edge links
    os << ", nodelinks:" 
        << e.node[0]->getIndex() << " "
        << e.node[1]->getIndex();
    
   os << ", facelinks:";
   int i;
  for ( i=0; i<e.face.NumElements(); i++) 
          os << e.face[i]->getIndex() << " ";
                
        // print out the tetra links
  os << ", tetralinks:";
  for ( i=0; i<e.tetra.NumElements(); i++) 
          os << e.tetra[i]->getIndex() << " ";
    
    return(os);
}
Face* Edge::getOtherFace(Face* f)
{
        if (face.NumElements() > 2) {
                cout << "can't use this function with this mesh" << endl;
                return  NULL;
        }
        
        if (face.NumElements() == 1) {
                if (face[0] == f) return (NULL);
        }
        
        if (face[0] == f) 
                return(face[1]);

  else if (face[1] == f)
                return(face[0]);

  else return(NULL);
}
int Edge::stepOffVerts(Point3D intPt) 
{
        for (int i=0; i<2; i++) {
                if (node[i]->p.Dist(intPt) < 0.00001) { 
                        //cout << "intersection on a vert" << endl;
                        intPt = .01*(node[(i+1)%3]->p - node[(i+0)%3]->p);                              
                        return i;
                }
        }
        return -1;
}
int Edge::intersectsPlaneSweptByCuttingEdge(Edge* e_in, Point3D* intersection_p)
{ return 0;}

int Edge::intersectsPlaneSweptFromEntryPt(Point3D entPt, Edge* cuttingEdge,Point3D *intPt)
{                               

        Point3D p2(cuttingEdge->getNodeLink(0)->p.x,
                                                 cuttingEdge->getNodeLink(0)->p.y,
                                                 cuttingEdge->getNodeLink(0)->p.z);
        Point3D p3(cuttingEdge->getNodeLink(1)->p.x,
                                                 cuttingEdge->getNodeLink(1)->p.y,
                                                 cuttingEdge->getNodeLink(1)->p.z);
                                
        int intersects = 0;
        intersects =crossesPlane(entPt, p2, p3, intPt);
//      cout << "crosses plane: " << intersects << endl;
        if (intersects) {
                intersects = intPt->insideTriangle(entPt, p2, p3);
        //      cout << "inside tri: " << intersects << endl;
        }
        return intersects;
}

// bump us back, based on them
void Edge::ResolveMyself(Edge* other_edge)
{
    // get locals
    double myradius = getRadius();
    double otherradius = other_edge->getRadius();
    double totalradius = otherradius + myradius;

    // loop for each of my nodes
    for (int i=0; i<2; i++) {
        Node* mynode = getNodeLink(i);
        // loop for each of their nodes
        for (int j=0; j<2; j++) {
            Node* othernode = other_edge->getNodeLink(j);

            // test distance
            double d = mynode->p.Dist(othernode->p);

            // if my half impinges on theirs, bump me back
            if (d < totalradius) 
                mynode->p += (mynode->p - othernode->p) * (totalradius-d);
        }
    }
}

---