filament.cc

Go to the documentation of this file.

/*
 * Filament detection routines
 */
 
#include "FMatrix.h"
#include "filament.h"
 
#ifdef I
#undef I
#endif
 
// this table is the mass matrix for the 1st subtriangle of the quadrilateral
// the rest of the mass matrices are just permutations
static const int TriFaces[1][3] = {
  { 0, 1, 2 }
};
 
static const int TetFaces[4][3] = {
  { 0, 1, 2 },
  { 0, 1, 3 },
  { 1, 2, 3 },
  { 0, 2, 3 }
};
 
 
#define FIL_PTS_FILE  0
#define FIL_ELEM_FILE 1
#define FIL_DATA_FILE 2
 
const int MAX_ELEM_NODES = 8;
const int MAX_ELEM_FACES = 6;
 
// prototypes
void output_filaments(filament *f, FILE *fh[]);
 
bool fl_chkIso(float *d, float th, Elem *e);
bool fl_IsoTagFaces   (float *d,float th, Elem *e, bool *faceTags );
bool fl_IsoTagTetFaces(float *d,float th, bool *faceTags);
bool fl_IsoTagTriFaces(float *d,float th, bool *faceTags);
bool fl_IsoTagTriangle(float *d, float th, const int idx[]);
bool fl_findIsoIntersect(float *d0, float *d1, float th, Point *p, Elem *e );
int  fl_FindFilament(float *d0, float *d1, float th, Point *lp, Elem *e, singularity *sng);
bool fl_triFindPS(float *d0, float *d1, float th, Point *p, Point *PS,int meth);
bool fl_triFindPS_Lines(float *d0, float *d1, float th, Point *p, Point *PS);
bool fl_triFindPS_ShpFnc(float *d0,float *d1,float th,Point *p,Point *PS);
bool fl_findEdgeSection(float d_p0, float d_p1, float th, Point *p0, Point *p1, Point *pSec);
void interpoate_vm(int numNodes,float **d,int M,float **rd,float rDist);
void interpolateTBuffer(tbuffer *tbf);
float fl_triArea( Point *pts );
 
// test point in triangle
bool SameSideOfLine (Point pTest, Point p0, Point p1, Point p2);
bool PointInTriangle(Point pTest, Point p0, Point p1, Point p2);
bool PointInTriangleBarycentric(Point pTest, Point p0, Point p1, Point p2);
bool PointInTriangleBarycentric(Point pTest, Point p0, Point p1, Point p2, Point &bcc);
bool PtInTriangleBarycentricSP(Point pTest, Point p1, Point p2);
bool PtInTriangleBarycentricSP(Point pTest, Point p1, Point p2, Point &bcc);
 
 
Point
UnTransformPoint( Point PTs, Point *opts, Point *spts)
{
  Real whole = fabs(spts[1].x*spts[2].y/2.);
  Point a[3] = { PTs, spts[1], spts[2] };
  Real c0 = fl_triArea( a )/whole;
  Real c2 = fabs(spts[1].x*PTs.y/2.)/whole;
  if( c0+c2 > 1.02 ){
    fprintf( stderr, "transformation back failed\n" );
  }
  return opts[0]*c0 + opts[2]*c2 + opts[1]*(1-c0-c2);
}
 
 
inline float
fl_triArea( Point *pts )
{
  return 0.5*mag(cross(pts[1] - pts[0],pts[2] - pts[0]));
}
 
Point
fl_triNormalLocal( Point *pts, bool u )
{
  Point n=cross(pts[1] - pts[0],pts[2] - pts[0]);
  return u?normalize(n):n;
}
 
 
void
interpolateTBuffer(tbuffer *tbf)
{
  int     numNodes = tbf->nodes;
  float **d        = (float **)(tbf->d );
  float **rd       = (float **)(tbf->dr);
  int     M        = tbf->wdth;
  float   rdist    = tbf->rdist;
 
  interpolate_vm(numNodes,d,M,rd,rdist);
}
 
 
void
interpolate_vm(int numNodes,float **d,int M,float **rd,float rdist)
{
  for(int j=0;j<numNodes;j++)  {
    rd[0][j] = d[0  ][j]+rdist*(d[1  ][j]-d[0  ][j]);
    rd[1][j] = d[M-2][j]+rdist*(d[M-1][j]-d[M-2][j]);
  }
}
 
 
bool
compute_filament(float *d0, float *d1, float vth, ElemList *elst,
                 NodeList *nlst, filament *f, SingularityDetector_t meth)
{
  int   PScnt     = 0;
 
#pragma omp parallel for
  for(int i=0;i<elst->NElems; i++)  {
 
    float d0_elem[MAX_ELEM_NODES];
    float d1_elem[MAX_ELEM_NODES];
 
    Elem &e = elst->Ele[i];
 
    if((e.type!=Tri) && (e.type!=Tetra))
      log_msg(NULL,2,0, "%s: Type in Element %d not supported yet.\n",
                             __func__, i );
 
    // copy data at element nodes to local array
    for(int j=0;j<e.Nnode;j++)  {
      d0_elem[j] =  d0[e.N[j]];
      d1_elem[j] =  d1[e.N[j]];
    }
 
    // quick check, are there any isolines/surfaces on the element
    bool iso0 = fl_chkIso(d0_elem,vth,&e);
    bool iso1 = fl_chkIso(d1_elem,vth,&e);
 
    // if so, analyse in detail
    if( iso0 && iso1 )  {
 
      // possibly a filament, at least we have two isosurfaces/isolines
      // is this actually a filament?
      bool ftags0[MAX_ELEM_FACES] = {false};
      bool ftags1[MAX_ELEM_FACES] = {false};
 
      // all faces with with isolines tagged
      bool tags0 = fl_IsoTagFaces(d0_elem, vth, &e, ftags0 );
      bool tags1 = fl_IsoTagFaces(d1_elem, vth, &e, ftags1 );
 
      Point  lp[MAX_ELEM_NODES];
      for( int j=0; j<e.Nnode; j++ ) 
        lp[j] = nlst->Node[e.N[j]];
 
      if(tags0 && tags1) {
 
        if((e.type==Tri) && (meth==vol_based))  {
          lp[3] = lp[0] + fl_triNormalLocal(lp,false);
          d0_elem[3] = d0_elem[0];
          d1_elem[3] = d1_elem[0];
        }
 
        singularity sng;
        int FilamentType = fl_FindFilament(d0_elem,d1_elem,vth,lp,&e,&sng);
 
        if(FilamentType)  {
#pragma omp critical
          f->add_singularity(&sng,i);
#pragma omp atomic
          PScnt++;
        }
 
        //bool isFilament = fl_findIsoIntersect(d0_elem,d1_elem,vth,lp,e);
        //if(isFilament)
        //  add_singularity2filament(f, elst->ElemGlobInds[i] );
 
      }
    }
  }
  f->sort();
  return PScnt>0;
}
 
void
filament::add_singularity(singularity *sng, int eIndx)
{
  s.push_back(*sng);
  numPts            += sng->ps.type;
  s.back().ps.eIndx = eIndx;
}
 
 
void
filament::sort()
{
  std::sort(s.begin(), s.end(), 
  [] (const singularity& lhs, const singularity& rhs)-> bool{return lhs.ps.eIndx<rhs.ps.eIndx;});
}
 
 
int
init_filament_IO(const char *meshname, int numT, FILE *fh[])
{
  int ok = -1;
 
  char pname[512], ename[512], dname[512];
 
  snprintf(pname, sizeof pname, "%s.pts_t",  meshname);
  snprintf(ename, sizeof ename, "%s.elem_t", meshname);
  snprintf(dname, sizeof dname, "%s.dat_t",  meshname);
 
  static int init = 0;
  if(!init)  {
 
    FILE *fhp = fopen(pname,"wt"); fh[0] = fhp;
    FILE *fhe = fopen(ename,"wt"); fh[1] = fhe;
    FILE *fhd = fopen(dname,"wt"); fh[2] = fhd;
 
    if( (fhp==NULL) || (fhe==NULL) || (fhd==NULL) )
      return ok;
    else
      ok = 0;
 
    // write comment lines
    fprintf(fhp,"## Format x y z # singularity type (0=PS || 1=FilSeg) Global element index \n");
    fprintf(fhp,"%d\n", numT);
    fprintf(fhe,"## Format Ln p0 p1 || PS p0 \n");
    fprintf(fhe,"%d\n", numT);
    fprintf(fhd,"## Format data (matching with pts_t file) \n");
    fprintf(fhd,"%d\n", numT);
 
    init = 1;
  }
 
  return ok;
}
 
 
// temporary debug code
int
write_connections(ElemList *elst, NodeList *nlst, filament *f, FILE *filHdls[],
                  float t)
{
  static int init = 0;
 
  if(!f->num_segs()) return 0;
 
  char fname[512];
  char pname[512], ename[512], dname[512];
 
  snprintf(fname, sizeof fname, "filament_%.2f.cnnx", t);
  FILE *fh = fopen(fname,"wt");
  fprintf( fh,"%d\n", f->num_segs());
  int i = 0;
 
  for( int i=0; i<f->num_segs(); i++ ) { 
    // pick the first edge of an element
    Elem &e = elst->Ele[f->s[i].ps.eIndx];
    //fl_getElem(elst,f->s[i].ps.eIndx,&e);
    int vtx0 = e.N[0];
    int vtx1 = e.N[1];
 
    fprintf( fh,"%d %d\n", vtx0, vtx1 );
  }
  fclose(fh);
 
  return f->num_segs();
}
 
 
void
write_filaments(filament *f,FILE *filHdls[],float t)
{
  int msg1    = 20;
  int msg2    = 30;
  int N       = f->num_segs();
  int numPts  = f->num_pts();
 
  for(int i=0; i<=FIL_DATA_FILE;i++) 
    fprintf(filHdls[i],"#%.6f\n",t);
 
  // preliminary, this needs to be fixed
  fprintf(filHdls[FIL_PTS_FILE ],"%d\n",numPts);
  //fprintf(filHdls[FIL_ELEM_FILE],"%d\n",N);
  fprintf(filHdls[FIL_DATA_FILE],"%d\n",numPts);
 
  fprintf(stderr, "%d filament segments found\n",N);
  output_filaments(f,filHdls);
}
 
 
void
output_filaments(filament *f, FILE *filHdls[])
{
  int nele=0;
  for( int i=0; i<f->num_segs(); i++ )
    if( f->s[i].ps.type != PhaseSingularity ) nele++;
  fprintf(filHdls[FIL_ELEM_FILE],"%d\n",nele);
 
  int ptCnt = 0;
  for( int i=0; i<f->num_segs(); i++ ){
    FILE *fhp = filHdls[0], *fhe = filHdls[1], *fhd = filHdls[2];
 
    FilSeg fs = f->s[i].fs;
    fprintf(fhp,"%.3f %.3f %.3f # %d %d\n", fs.p0.x,fs.p0.y,fs.p0.z,
                                            fs.type, fs.eIndx);
    // for now we use the enclosing element index as data
    // later we use a PS or filament ID
    fprintf(fhd,"%d\n",fs.eIndx);
 
    if(f->s[i].ps.type==PhaseSingularity)  {
      //fprintf(fhe,"PS %d \n", ptCnt++);
      ptCnt++;
    }
    else  {
      // filament segment
      fprintf(fhp,"%.3f %.3f %.3f # %d %d\n", fs.p1.x,fs.p1.y,fs.p1.z,
                                              fs.type, fs.eIndx);
      fprintf(fhd,"%d\n",fs.eIndx);
      fprintf(fhe,"Ln %d %d\n", ptCnt, ptCnt+1);
      ptCnt += 2;
    }
  }
}
 
 
void
close_filament_IO(FILE *fh[])
{
  for(int i=0;i<3;i++)
    fclose(fh[i]);
}
 
bool
fl_chkIso(float *d, float th, Elem *e)
{
  if(d[0]>=th)  {
    for(int i=1;i<e->Nnode;i++)
      if(d[i]<th)  {
        return true;
      }
  } else {
    for(int i=1;i<e->Nnode;i++)
      if(d[i]>=th)  {
        return true;
      }
  }
  return false;
}
 
 
bool
fl_IsoTagFaces(float *d, float th, Elem *e, bool *faceTags )
{
  bool flg = false;
  switch(e->type)  {
    case Tri:
      flg = fl_IsoTagTriFaces(d,th,faceTags);
      break;
    case Tetra:
      flg = fl_IsoTagTetFaces(d,th,faceTags);
      break;
    default:
      fprintf(stderr,"%s: Element type not implemented yet.\n", __func__);
  }
  return flg;
}
 
 
bool
fl_IsoTagTetFaces(float *d,float th, bool *faceTags)
{
  int tagged = 0;
  for(int i=0;i<4;i++)
    if((faceTags[i]=fl_IsoTagTriangle(d,th,TetFaces[i])))
      tagged++;
 
  return tagged;
}
 
 
bool
fl_IsoTagTriFaces(float *d,float th, bool *faceTags)
{
  int tagged = 0;
  for(int i=0;i<3;i++)
    if((faceTags[i]=fl_IsoTagTriangle(d,th,TriFaces[i])))
      tagged++;
 
  return tagged;
}
 
 
bool
fl_IsoTagTriangle(float *d, float th, const int idx[])
{
  bool tag = false;
 
  if(d[idx[0]]>=th)  {
    if((d[idx[1]]<th) || (d[idx[2]]<th)) tag = true;
  }
  else
    if((d[idx[1]]>=th) || (d[idx[2]]>=th)) tag = true;
  return tag;
}
 
 
int
fl_FindFilament(float *d0,float *d1,float th,Point *lp,Elem *e,singularity *sng)
{
  float d0_[3], d1_[3];
  Point lp_[3];
  int FilamentType = 0;
  Point PhaseS[4];
  int tagged[] = { 0, 0, 0, 0};
  int tagIdx[] = {-1,-1,-1,-1};
 
  int meth = 1;
  switch(e->type)  {
    case Tri:
      if(fl_triFindPS(d0,d1,th,lp,PhaseS,meth))  {
        int i = 0;
        tagged[i]++;
        tagIdx[FilamentType++] = i;
      }
      break;
    case Tetra:
      for(int i=0;i<4;i++)  {
        for(int j=0;j<3;j++)  {
          d0_[j] = d0[TetFaces[i][j]];
          d1_[j] = d1[TetFaces[i][j]];
          lp_[j] = lp[TetFaces[i][j]];
        }
        if(fl_triFindPS(d0_,d1_,th,lp_,PhaseS+i,meth))  {
          tagged[i]++;
          tagIdx[FilamentType++] = i;
        }
      }
      break;
    default:
      fprintf(stderr,"Element type not supported by %s.\n", __func__ );
  }
 
  if(FilamentType)  {
    sng->ps.type = FilamentType;
    sng->ps.eIndx = 0;
    sng->ps.p0 = PhaseS[tagIdx[0]];
    if(FilamentType==2)  // Filament segment
      sng->fs.p1 = PhaseS[tagIdx[1]];
  }
 
  return FilamentType;
}
 
 
bool
fl_triFindPS(float *d0, float *d1, float th, Point *lp, Point *PhaseS,
                   int meth)
{
  bool found = false;
  switch(meth)  {
    case 0:
      found = fl_triFindPS_ShpFnc(d0,d1,th,lp,PhaseS);
      break;
    case 1:
    default:
      found = fl_triFindPS_Lines(d0,d1,th,lp,PhaseS);
  }
  return found;
}
 
 
bool
fl_triFindPS_ShpFnc(float *d0, float *d1, float th, Point *lp, Point *PhaseS)
{
  fprintf(stderr, "Shape function method not implemented yet.\n");
  return false;
}
 
bool
fl_triFindPS_Lines(float *d0, float *d1, float th, Point *lp, Point *PhaseS)
{
  bool isPS = false;
  Point pt[3]; 
 
  // transform matrix into standard position first
  //static FMatrix fwd = ZERO_FMATRIX, bwd = ZERO_FMATRIX, A = ZERO_FMATRIX;
 
  FMatrix A(2,2);
  
  /*
  fl_getTransform3Dto2D(lp,&fwd,&bwd,&t);
  int numPts = 3;
  Point t;
  fl_AffineTransformPoints(lp,pt,numPts,&fwd,&t,true);
  */
 
  // transform to standard position
  pt[0].set(0,0,0);
  Point p01 = lp[1] - lp[0];
  pt[1].set(mag(p01), 0., 0.);
  Point p02 = lp[2] - lp[0];
  double ptx = dot(p02,p01)/pt[1].x;
  pt[2].set(ptx, sqrt(mag2(p02)-ptx*ptx),0.);
 
  int tEdg0[] = {0,0,0};
  int tEdg1[] = {0,0,0};
  Point tSec0[3], tSec1[3];
 
  // do our thingie
  tEdg0[0] = fl_findEdgeSection(d0[0],d0[1],th,pt+0,pt+1,tSec0+0);
  tEdg0[1] = fl_findEdgeSection(d0[1],d0[2],th,pt+1,pt+2,tSec0+1);
  tEdg0[2] = fl_findEdgeSection(d0[0],d0[2],th,pt+0,pt+2,tSec0+2);
 
  tEdg1[0] = fl_findEdgeSection(d1[0],d1[1],th,pt+0,pt+1,tSec1+0);
  tEdg1[1] = fl_findEdgeSection(d1[1],d1[2],th,pt+1,pt+2,tSec1+1);
  tEdg1[2] = fl_findEdgeSection(d1[0],d1[2],th,pt+0,pt+2,tSec1+2);
 
  // singularity found?
  int t0 = 0;
  int t1 = 0;
  int t0s[2], t1s[2];
  int c0 = 0;
  int c1 = 0;
  for(int i=0;i<3;i++)  {
    t0 += tEdg0[i];
    t1 += tEdg1[i];
    if(tEdg0[i]) t0s[c0++] = i;
    if(tEdg1[i]) t1s[c1++] = i;
  }
  if((t0==2) && (t1==2))  {
    // equation in parametrized form
    // a0, b0: intersection points of isoline 0 with tri
    // a1, b1: intersection points of isoline 1 with tri
    //
    // There are two lines:
    // l0 = a0 + lambda*(b0-a0)     (1)
    // l1 = a1 + gamma *(b1-a1)     (2)
    //
    // The intersection is found by solving for l0=l1.
    // this yields:
    //
    // (b0-a0)*lambda - (b1-a1)*gamma = (a1-a0)
    //
    // which gives a solution for lambda and gamma.
    // Either (1) or (2) can be used to find the intersection coordinates
 
    Point &a0 = tSec0[t0s[0]];
    Point &b0 = tSec0[t0s[1]];
    Point &a1 = tSec1[t1s[0]];
    Point &b1 = tSec1[t1s[1]];
 
    Point b0ma0 = b0 - a0;
    Point b1ma1 = a1 - b1;
    Point a1ma0 = a1 - a0;
 
    Real b[2];
    A.Ent[0][0] = b0ma0.x; A.Ent[0][1] = b1ma1.x;
    A.Ent[1][0] = b0ma0.y; A.Ent[1][1] = b1ma1.y;
    b[0]        = a1ma0.x; b[1]        = a1ma0.y;
 
    // find params lambda and gamma
    A.LUD();
    if(A.decomposed)  {
      // A is not singular
      A.LUsolve( b );
 
      // find intersection
      Point PSt = a0 + scal_X(b0ma0,b[0]);
 
      // check whether PS is within triangle
      Point bcc;
 
      if(PtInTriangleBarycentricSP(PSt,pt[1],pt[2],bcc))  {
         //fprintf(stderr, "PS transformed found at x = %.2f, y = %.2f.\n", PSt.x, PSt.y );
 
        // transform back
        // fl_AffineTransformPoints(&PSt,PhaseS,1,&bwd,&t,false);
        // PhaseS[0] = UnTransformPoint(PSt, lp, pt);
        
        PhaseS[0] = lp[0]*bcc.x + lp[1]*bcc.y + lp[2]*bcc.z;
 
        // fprintf(stderr, "PS real found at x = %.2f, y = %.2f.\n", PhaseS->x, PhaseS->y );
        isPS = true;
 
        bool debug = false;
        if(debug)  {
          // check back transform - debug only!
          Point ptb[3];
          //fl_AffineTransformPoints(pt,ptb,3,&bwd,&t,false);
          for( int j=0; j<3; j++ )
            ptb[j] = dot(bcc,pt[j]);
        }
      }
    }
  }
  return isPS;
}
 
 
bool
SameSideOfLine(Point pTest, Point p0, Point p1, Point p2)
{
  Point cp1 = cross(p1 - p0,pTest - p0);
  Point cp2 = cross(p1 - p0,p2 - p0);
  return  dot(cp1,cp2)>=0;
}
 
 
bool
PointInTriangle(Point pTest, Point p0, Point p1, Point p2)
{
  if(SameSideOfLine(pTest,p0,p1,p2) &&
     SameSideOfLine(pTest,p1,p2,p0) &&
     SameSideOfLine(pTest,p2,p0,p1)    )
    return true;
  else
    return false;
}
 
bool
PtInTriangleBarycentricSP(Point pTest, Point p1, Point p2)
{
  Point null;
  return PtInTriangleBarycentricSP(pTest, p1, p2, null);
}
 
bool
PtInTriangleBarycentricSP(Point v2, Point v1, Point v0, Point &bcc)
{
  // Compute dot products
  double dot00 = v0.x*v0.x+v0.y*v0.y;
  double dot01 = v0.x*v1.x;
  double dot02 = v0.x*v2.x+v0.y*v2.y;
  double dot11 = v1.x*v1.x;
  double dot12 = v1.x*v2.x;
  
  // Compute barycentric coordinates
  double invDenom = 1/(dot00*dot11-dot01*dot01);
  bcc.z = (dot11*dot02-dot01*dot12)*invDenom;
  if( bcc.z < 0 ) return false;  
  bcc.y = (dot00*dot12-dot01*dot02)*invDenom;
  bcc.x = 1. - bcc.z - bcc.y;
  
  // Check if point is in triangle
  return (bcc.x>=0) && (bcc.y>=0);
}
 
 
bool
PointInTriangleBarycentric(Point pTest, Point p0, Point p1, Point p2)
{
  Point null;
  return PointInTriangleBarycentric(pTest, p0, p1, p2, null);
}
 
bool
PointInTriangleBarycentric(Point pTest, Point p0, Point p1, Point p2, Point &bcc)
{
  // Compute vectors
  Point v0 = p2 - p0;
  Point v1 = p1 - p0;
  Point v2 = pTest - p0;
 
  // Compute dot products
  double dot00 = dot(v0,v0);
  double dot01 = dot(v0,v1);
  double dot02 = dot(v0,v2);
  double dot11 = dot(v1,v1);
  double dot12 = dot(v1,v2);
  
  // Compute barycentric coordinates
  double invDenom = 1/(dot00*dot11-dot01*dot01);
  bcc.z = (dot11*dot02-dot01*dot12)*invDenom;
  bcc.y = (dot00*dot12-dot01*dot02)*invDenom;
  bcc.x = 1. - bcc.z - bcc.y;
  
  // Check if point is in triangle
  return (bcc.x>=0) && (bcc.y>=0) && (bcc.z>=0);
}
 
 
bool
fl_findEdgeSection(float d_p0, float d_p1, float th, Point *p0, Point *p1, Point *pSec)
{
  bool hasIso = false;
 
  if(d_p0>=th)  {
    if (d_p1<th)
      hasIso = true;
  }
  else  {
    if (d_p1>=th)
      hasIso = true;
  }
 
  if(hasIso)  {
    Point l   = *p1 - *p0;
    float len = mag(l);
    float sec = (d_p0-th)/(d_p0-d_p1);
    *pSec = *p0 + scal_X(l,sec);
  }
 
  return hasIso;
}
 
 
bool
fl_findIsoIntersect(float *d0, float *d1, float th, Point *p, Elem *e )
{
  bool isFilament = false;
 
  FMatrix M(4,4);
  FMatrix Id(3,3);
 
  for (int i=0;i<4;i++) {
    M.Ent[i][0] = 1.;
    M.Ent[i][1] = p[i].x;
    M.Ent[i][2] = p[i].y;
    M.Ent[i][3] = p[i].z;
  }
 
  double v0[4], v1[4], vf[4];
  for(int i=0;i<4;i++) v0[i] = d0[i];
  for(int i=0;i<4;i++) v1[i] = d1[i];
//  for(int i=0;i<4;i++) vf[i] = phiIsoFace[i];
 
  M.LUD();
 
  // find base coefficients V(x,y,z) = a + bx +cy + dz
  M.LUsolve( v0 );
  M.LUsolve( v1 );
 
  // now we try to intesect the two hyperplanes with isosurface planes
  // spanned by the faces of the tet, i.e. we assume a face is an isosurface hyperplane
  int intersect = 0, NFaces = 4;
  for(int i=0;i<NFaces;i++)  {
    // assign threshold values to nodes which span isosurface
    for(int j=0;j<4;j++) vf[j] = th+1;
    for(int j=0;j<3;j++) vf[TetFaces[i][j]] = th;
 
    // get base coeffs for this hyperplane
    M.LUsolve( vf );
 
    // fill in rhs which is vth-a
    double r[3];
    r[0] = th - v0[0];
    r[1] = th - v1[0];
    r[2] = th - vf[0];
 
    for(int j=0;j<3;j++)  {
      Id.Ent[0][j] = v0[j+1];
      Id.Ent[1][j] = v1[j+1];
      Id.Ent[2][j] = vf[j+1];
    }
 
    Id.LUD( );
    if(Id.decomposed)  {
      Id.LUsolve( r );
 
      // this is primitive and not exactly correct,
      // a proper test of whether the intersection is within the tet or not
      // should go here instead
      float xmn,xmx,ymn,ymx,zmn,zmx;
      xmn = xmx = p[0].x;
      ymn = ymx = p[0].y;
      zmn = zmx = p[0].z;
 
      for(int j=1;j<4;j++)  {
        if(p[j].x<xmn) xmn = p[j].x;
        if(p[j].x>xmx) xmx = p[j].x;
        if(p[j].y<ymn) ymn = p[j].y;
        if(p[j].y>ymx) ymx = p[j].y;
        if(p[j].z<zmn) zmn = p[j].z;
        if(p[j].z>zmx) zmx = p[j].z;
      }
 
      if((r[0]>=xmn) && (r[0]<=xmx) &&
         (r[1]>=ymn) && (r[1]<=ymx) &&
         (r[2]>=zmn) && (r[2]<=zmx) )  {
        intersect++;
        log_msg(NULL,0,0, "Intersection found at [x,y,z]=[%.2f,%.2f,%.2f]",
                               r[0],r[1],r[2]);
      }
    }
 
  }
  if( ((e->type==Tri)   && (intersect==1)) ||
      ((e->type==Tetra) && (intersect==2))    )
      isFilament = true;
 
  return isFilament;
}