BFLInterpolation.cxx

Go to the documentation of this file.

//                                                                   //
//authors: Costas Andreopoulos, Elias Athanasopoulos, George Tzanakos//
//                Physics Department, University of Athens           //

#include "TObjArray.h"
#include "TMath.h"

//#include "TObjectTable.h"

//#include <TVector3.h>
//#include <TROOT.h>

#include "BField/TIntList.h"
//#include "BField/BFLNode.h"
//#include "BField/BFLVtx.h"
//#include "BField/BFLEdge.h"
//#include "BField/BFLPolyg.h"
#include "BField/BFLWingedEdge.h"
#include "BField/BFLVorOperator.h"
#include "BField/BFLInterpolation.h"   
#include "BField/BFLCache.h"
#include "BField/BFLAnsysLookup.h"
#include "BField/BfldLoanPool.h"
#include "BField/BfldMeshVoronoi.h"
#include "Conventions/Munits.h"
#include "MessageService/MsgService.h"

#include <fstream>

CVSID("$Id: BFLInterpolation.cxx,v 1.12 2005/02/03 22:35:49 rhatcher Exp $");

ClassImp(BFLInterpolation)

class BfldHandyMath {
public:
  BfldHandyMath(void) {
  }

  static double Det(double a11,double a12,double a13,
                   double a21,double a22,double a23,
                   double a31,double a32,double a33) {
    return a11 * a22 * a33 - a11 * a23 * a32 - a12 * a21 * a33
      + a12 * a23 * a31 + a13 * a21 * a32 - a13 * a22 * a31;
  }

  static double Det(const TVector3 &v1,const TVector3 &v2,
                   const TVector3 &v3) {
    return Det(v1.X(),v1.Y(),v1.Z(),
               v2.X(),v2.Y(),v2.Z(),
               v3.X(),v3.Y(),v3.Z());
  }

  // Assuming v1,v2,v3 are in the xy plane, determine whether the
  // given order is counter-clockwise or clockwise by looking at
  // the sign of the z coordinate of (v2 - v1) x (v3 - v1).
  static bool IsCCW(const TVector3 &v1,const TVector3 &v2,
                    const TVector3 &v3) { 
    return (v2.X() - v1.X()) * (v3.Y() - v1.Y())
      - (v2.Y() - v1.Y()) * (v3.X() - v1.X())  > 0;
  }

  // Test to see if v is in the "wedge" formed by v1,v2,v3, with v1 as the
  // "pivot" of the wedge. Does the degenerate case in which any of these
  // vector triplets is collinear occur?

  static bool IsInWedge(const TVector3 &v,const TVector3 &v1,
                        const TVector3 &v2,const TVector3 &v3) {
    bool flag1 = IsCCW(v2,v1,v3),
      flag2 = IsCCW(v2,v1,v),
      flag3 = IsCCW(v1,v3,v);
    return ((flag1 && flag2 && flag3) || (!flag1 && !flag2 && !flag3));
  }

  static bool IsInTri(const TVector3 &v,const TVector3 &v1,
                      const TVector3 &v2,const TVector3 &v3) {
    bool flag1 = IsCCW(v2,v1,v3),
      flag2 = IsCCW(v2,v1,v),
      flag3 = IsCCW(v1,v3,v),
      flag4 = IsCCW(v2,v,v3);
    return ((flag1 && flag2 && flag3 && flag4)
            || (!flag1 && !flag2 && !flag3 && !flag4));

  }
};

ofstream *gBfldtestFile;

BFLInterpolation::BFLInterpolation(TObjArray * bfield, 
                                   BFLWingedEdge * voronoi)
{ 
  fBField = bfield;                  /* BF Vectors */ 
  fVor = voronoi;                    /* Voronoi diagram */
  fOp = new BFLVorOperator(voronoi); /* Object that operates on fVor */  
  fNPolygs = fVor->GetNofPolygons() - 3;

  fCache = new BFLCache();
  SetInterpolant(BfldInterpMethod::kClosest);
  
  fPolygons = new TIntList();
  fVertices = new TObjArray(20); 

  fSearchSeed = 1; 

  gBfldtestFile = new ofstream("BfldPlanarInterpErrs.err",ios::app);
} 

BFLInterpolation::~BFLInterpolation() 
{
  fPolygons->Delete();
  delete fPolygons;    fPolygons=0;
  fVertices->Delete();
  delete fVertices;    fVertices=0;
  delete fOp;          fOp=0;
  delete fVor;         fVor=0;
  fBField->Delete();
  delete fBField;      fBField=0;
  delete fCache;       fCache=0;

  delete gBfldtestFile;
}

TVector3 BFLInterpolation::GetB(const TVector3& PosVector,
                                const Bool_t updateSeedPolygon)
{ 
 Int_t istatus;
 TVector3 BF(0.,0.,0.);
 BFLNode point(PosVector.X(),PosVector.Y()); 
 
 switch(fHowToInterpolate) {
   case(BfldInterpMethod::kNatural): 
     istatus = FormVoronoiCell(&point);
     if(istatus == kFail) { // for some reason
        BF = CNInterpolation(&point); // do the next best thing
     } else {
        BF = NNInterpolation(&point);
        ResetForNextQuery();
     }
     break;   
   case(BfldInterpMethod::kClosest): 
     BF = CNInterpolation(&point);
     break;
   case(BfldInterpMethod::kBilinear):
     BF = BilinearInterpolation(&point);
     break;
 case(BfldInterpMethod::kPlanar):
   BF = PlanarInterpolation(&point); //Added by Andrew Goldstone
   break;
   default:
        MSG("BFL",Msg::kWarning) 
           << "You did not specified a valid interpolation method" << endl;
 }

 if (updateSeedPolygon) SetSearchSeed(GetLastPolygon());
 return BF; 
}

Int_t BFLInterpolation::FormVoronoiCell(BFLNode * point)
{
  Int_t EdgeID, StartEdgeID, LastEdgeID, iedge, ivtx = 0; 
//rwh:  TIntList * IntersEdges = new TIntList();
 
  fOp->SetCurrentGen(point); 
  fCurrentPolygID = fOp->FindCurrentPolygon(fSearchSeed); 
  
  fOp->SetCurrentPolygID(fCurrentPolygID); 
  fPolygons->Add(fCurrentPolygID); 
  
  EdgeID = fOp->GetFirstIntersectEdge();
  if(EdgeID == -1) return kFail;
  
  StartEdgeID = EdgeID; //keep this id to terminate the next loop
  LastEdgeID = 0; // removes dependencies from the previous iteration.

  do {       
      // store intersect edges for this polygon
//rwh:      IntersEdges->Add(EdgeID); 

      // Get the intersection point -> new vtx.
      // this will get deleted in ResetForNextQuery
      fVertices->Add(new BFLNode(fOp->FindNewVtx(EdgeID)));
                                                                                                                
      // Get the next intersect edge
      iedge = fOp->GetNextIntersectEdge(EdgeID);
      if (iedge == LastEdgeID ) { 
          // It is time to move to the next polygon.

          // Get the next polygon ID.
          fCurrentPolygID = fOp->MoveToNextPolygon(fCurrentPolygID,EdgeID);
          fOp->SetCurrentPolygID(fCurrentPolygID);
          fPolygons->Add(fCurrentPolygID);
  
          LastEdgeID = EdgeID;
          EdgeID = fOp->GetNextIntersectEdge(EdgeID);                                         

       } else {       
          LastEdgeID = EdgeID;
          EdgeID = iedge;            
       }

       ivtx++;
       if(ivtx > kMaxVtx) return kFail;
       
  } while(EdgeID != StartEdgeID); 

  // this will get deleted in ResetForNextQuery
  fVertices->Add(new BFLNode(fOp->FindNewVtx(EdgeID))); 
  
  return 0;
} 

TVector3 BFLInterpolation::NNInterpolation(BFLNode * point)
{
  Int_t NatNeighbor = 0, ID = 0;
  Double_t AreaSum = 0, NeighborArea = 0, bx = 0, by = 0;

  // total area of new cell
  Double_t TotalArea = CellArea(fVertices); 
  
  // compute sub-areas
  for(Int_t inn=0; inn<fPolygons->NumberOfElements()-1; inn++) {
     NatNeighbor = fPolygons->At(inn);

     //rwh:     TObjArray * Area = new TObjArray(30); // vertices of sub-cell
     //rwh:  increase this size after seeing index out of bounds in AddAt()
     TObjArray * Area = new TObjArray(55); // vertices of sub-cell
     BFLNode * FirstVtx = (BFLNode *) fVertices->At(inn);
     BFLNode * SecondVtx = (BFLNode *) fVertices->At(inn+1);

     if (!FirstVtx || !SecondVtx) {
        //rwh: protect against SEGV
        MSG("BFL",Msg::kWarning) 
           << "failed fVertices->At(inn) " << inn 
           << " FirstVtx " << FirstVtx << " SecondVtx " << SecondVtx
           << endl
           << "   fPolygons->NumberOfElements " 
           << fPolygons->NumberOfElements() << endl;
        break;
     }

     Int_t Clockwiseness = fOp->Clockwise(FirstVtx,SecondVtx,point);

     // select other polygon's vertices on the same side with 'point' 
     TIntList * VtxAroundNeighbor = fOp->RetrieveVtxSurrPolygon(NatNeighbor);

     Int_t ivtxslot=0;
     for(Int_t ivtx=0; ivtx<VtxAroundNeighbor->NumberOfElements(); ivtx++) {
        Int_t VtxID = VtxAroundNeighbor->At(ivtx);
        BFLNode * Vtx = fVor->GetVtx(VtxID);       
        
        if(fOp->Clockwise(FirstVtx,SecondVtx,Vtx) == Clockwiseness) {
           Area->AddAt(Vtx,ivtxslot);
           ivtxslot++;
        } else {
           // fVor->GetVtx returns "new" BFLNode
           // if we're not saving it (to be deleted later) then toss it now
           delete Vtx;
        }
     }      

     VtxAroundNeighbor->Delete();
     delete VtxAroundNeighbor;  VtxAroundNeighbor=0;

     // ivtxslot = 0 means that the subcell was not formed
     if(ivtxslot == 0) {
        MSG("BFL",Msg::kWarning)
           << "BFLInterpolation::NNInterpolation --- subcell not formed" 
           << endl        
           << "    BFLNode 'point' " << point->GetNodeID() 
           << " x= " << point->GetX() << " y= " << point->GetY() << endl;
        return CNInterpolation(point);
     }

     // Make copies of the BFLNode's FirstVtx and SecondVtx
     // so that "Area" can _own_ them along with the entries
     // put in from a fVor->GetVtx call (which creates BFLNodes)
     BFLNode * FirstVtxCopy = new BFLNode(*FirstVtx);
     BFLNode * SecondVtxCopy = new BFLNode(*SecondVtx);

     if(ivtxslot == 1) {        
        Area->AddAt(FirstVtxCopy,1);
        Area->AddAt(SecondVtxCopy,2);
     } else {
        Int_t vcw = fOp->Clockwise((BFLNode *)Area->At(0),
                                   (BFLNode *)Area->At(1),FirstVtx);

        if(fOp->Clockwise((BFLNode *)Area->At(1),FirstVtx,SecondVtx) == vcw) {
           Area->AddAt(FirstVtxCopy,ivtxslot);
           Area->AddAt(SecondVtxCopy,ivtxslot+1);
        } else {
           Area->AddAt(SecondVtxCopy,ivtxslot);
           Area->AddAt(FirstVtxCopy,ivtxslot+1);
        }
     } 
     
     NeighborArea = CellArea(Area);
     AreaSum += NeighborArea;

     ID = fVor->GetGeneratorID(NatNeighbor);
     if (ID<0) {
        //rwh: protect against SEGV
        MSG("BFL",Msg::kWarning)
           << "BFLInterpolation::NNInterpolation bad ID: " << ID << endl;
     } else {
     bx += NeighborArea * ((TVector2 *) fBField->At(ID))->X();
     by += NeighborArea * ((TVector2 *) fBField->At(ID))->Y();
     }

     // Area owns some objects [created by fVor->GetVtx(VtxID)]
     Area->Delete();
     delete Area;  Area=0;
  }

  // For the last neighbor we use the fact that Sum(weights) = 1
  NatNeighbor = fPolygons->At(fPolygons->NumberOfElements() - 1);
  NeighborArea = TotalArea - AreaSum; 
  
  ID = fVor->GetGeneratorID(NatNeighbor);

  if (ID<0) {
     //rwh: protect against SEGV
     MSG("BFL",Msg::kWarning) 
        << "BFLInterpolation::NNInterpolation bad ID: " << ID << endl;
  } else {
  bx += NeighborArea * ((TVector2 *) fBField->At(ID))->X();
  by += NeighborArea * ((TVector2 *) fBField->At(ID))->Y();
  }

  bx /= TotalArea;
  by /= TotalArea;

  fOp->ResetVtxCache();

  return TVector3(bx,by,0); 
}

TVector3 BFLInterpolation::CNInterpolation(BFLNode * point)
{ 
  fOp->SetCurrentGen(point); 
  fCurrentPolygID = fOp->FindCurrentPolygon(fSearchSeed); 

  /*
  MSG("Bfldtest",Msg::kDebug)
    << "r is in polygon #" << fCurrentPolygID
    << " with surrounding polygons as follows:" << endl;

  TIntList *edgesList = fOp->RetrieveEdgesSurrPolygon(fCurrentPolygID);
  int n = edgesList->NumberOfElements();
  int i,edgeID, polyID,genID;
  for(i = 0;i < n;++i) {
    edgeID = edgesList->At(i);
    polyID = fVor->GetLeftPolyg(edgeID);
    if(polyID == fCurrentPolygID) {
      polyID = fVor->GetRightPolyg(edgeID);
    }
    MSG("Bfldtest",Msg::kDebug) << "Polygon #" << polyID << endl;
  }
  */
  if(fCurrentPolygID > 3) {
     Int_t ID = fVor->GetGeneratorID(fCurrentPolygID); 
     return *(TVector3 *) fBField->At(ID);
  } else return TVector3(0,0,0);
}

void BFLInterpolation::ResetForNextQuery(void)
{
  fVertices->Delete(); 
  fPolygons->Delete(); 
}  

Double_t BFLInterpolation::CellArea(TObjArray * vertices)
{
// Getting the area of a cell

  Double_t AreaOfCell = 0;
  BFLNode * Vtx;

  TMatrix area = TMatrix(3,3);
  for(Int_t i=0; i<3; i++) area(i,2)=1;

  area(0,0) = ((BFLNode *) vertices->At(0))->GetX();
  area(0,1) = ((BFLNode *) vertices->At(0))->GetY();
  
  Int_t ivtx=0;
  TIter nextv(vertices);
  while( (Vtx = (BFLNode *)nextv()) ) {    
    if(ivtx == 1) {
       area(1,0) = Vtx->GetX();
       area(1,1) = Vtx->GetY();
    }
    if(ivtx == 2) {
       area(2,0) = Vtx->GetX();
       area(2,1) = Vtx->GetY(); 
       AreaOfCell += 0.5 * area.Determinant();
       ivtx = 1;
       area(1,0) = area(2,0);
       area(1,1) = area(2,1);
    } 
    ivtx++;
  } 

  return TMath::Abs(AreaOfCell); 
}

TVector3 BFLInterpolation::BilinearInterpolation(BFLNode * point)       
{
  // Find the current polygon ID
  fOp->SetCurrentGen(point); 
  fCurrentPolygID = fOp->FindCurrentPolygon(fSearchSeed); 

  // Get the ANSYS Node that generates this Voronoi Cell
  BfldLoanPool * loanpool = BfldLoanPool::Instance();

  Int_t NodeID = ((BfldMeshVoronoi*)(loanpool->GetMesh(fGrid,0)))
                                ->FindANSYSGenerator(fCurrentPolygID);

  // Get the ANSYS Lookup Table
  BFLAnsysLookup lookup = fVor->GetAnsysLookupTable();

  // Get a list of all ANSYS Cells having the ANSYS Node
  // - we just found - as one of its corners
  TObjArray cells;
  lookup.FindNode(NodeID,&cells);

  // Select the ANSYS Cell that contains the user provided point
  TIter nextentry(&cells);
  BFLNode2ACell * TableEntry = 0, * SelectedTableEntry = 0;
  while( (TableEntry = (BFLNode2ACell*) nextentry() ) ) {

    Bool_t IsInside = IsInsideANSYSCell(TableEntry,point);
    if(IsInside) SelectedTableEntry = (BFLNode2ACell *) (TableEntry->Clone());
    if(IsInside) break;
  }

  // Get a list of all ANSYS Nodes of this cell
  // These are the points Jeff wants to use for 2-linear interpolation
  Int_t nodes[kNMaxNodesPerACell];
  SelectedTableEntry->GetCellNodes(nodes);

  // Elias should somehow calculate the correct weights
  // and return BField


  delete SelectedTableEntry;

  return TVector3(0,0,0); 
}


TVector3 BFLInterpolation::PlanarInterpolation(BFLNode *rNode) {

  fOp->SetCurrentGen(rNode);
  fCurrentPolygID = fOp->FindCurrentPolygon(fSearchSeed);

  if(fCurrentPolygID <= 3) {
    //This means the nearest generator is the point at infinity,
    //so we can safely return 0. See BFLInterpolation::CNInterpolation()
   return TVector3(0,0,0);
  }

  //Get the coordinates of the test point r
  TVector3 r(rNode->GetX(),rNode->GetY(),0);

  //Get generator coordinates
  TVector3 r1(fVor->GetGeneratorX(fVor->PolygonToGenerator(fCurrentPolygID)),
              fVor->GetGeneratorY(fVor->PolygonToGenerator(fCurrentPolygID)),
              0);

  //Get B field at closest generator. Note that we must map the ID of the
  //Voronoi cell into the ID of the B field for that cell's generator.
  //Note further that this is NOT the generator ID used by BFLWingedEdge.
  //In particular, BFLWingedEdge::GetGenerator[X,Y] takes the other ID,
  //which is, in fact, the same as the ID of the cell containing that
  //generator.
  int genID = fVor->GetGeneratorID(fCurrentPolygID);
  TVector3 b1 = *(TVector3 *)(fBField->At(genID));

  const double kDistanceThresholdSquared = 10.0; //In mm^2, alas

  if((r1-r).Mag2() < kDistanceThresholdSquared) {
    //r is essentially at a generator. In this case, it is not in the
    //interior of any Delaunay triangle, and we save ourselves the
    //trouble of a pointless interpolation.
    
    return b1;
  }

  //iterate through the generators of the neighboring cells until we
  //find the Delaunay triangle containing r

  TIntList *edgesList = fOp->RetrieveEdgesSurrPolygon(fCurrentPolygID);

  TVector3 r2,r3;
  int polyID,edgeID,gen2ID,gen3ID=-1;
  int i;

  //Apparently we get a list of surrounding edges in counterclockwise
  //order; see BFLVorOperator::RetrieveEdgesSurrPolygon(),
  //but apparently the edges do not have a common direction, so
  //we check which way they're directed to pick the right polygon.

  int n = edgesList->NumberOfElements();
  bool isFound = false;

  edgeID = edgesList->At(0);
  polyID = fVor->GetLeftPolyg(edgeID);
  if(polyID == fCurrentPolygID) {
    polyID = fVor->GetRightPolyg(edgeID);
  }

  //Take note: as mentioned in the discussion of generator ID's above,
  //the ID of the polygon and the ID of its generator are the same
  //for most calls to BFLWingedEdge. In fact,
  //BFLWingedEdge::PolygonToGenerator just returns its argument.
  gen2ID = fVor->PolygonToGenerator(polyID);
  r2.SetXYZ(fVor->GetGeneratorX(gen2ID),fVor->GetGeneratorY(gen2ID),0);

  /*
  MSG("Bfldtest",Msg::kDebug)
    << "Scanning generators around r = (" << r.X() << "," << r.Y()
    << ") in polygon #" << fCurrentPolygID << endl;

  MSG("Bfldtest",Msg::kDebug)
    << "ID's of cell edges: " << endl;
  for(i = 0;i < n;i++) {
    MSG("Bfldtest",Msg::kDebug) << "e" << i << ": #" << edgesList->At(i)
                               << endl;
  } 
  */

  for(i = 1;i <= n;++i) {
    edgeID = edgesList->At(i % n);
    polyID = fVor->GetLeftPolyg(edgeID);

    //As above, we don't know which way the edge is directed,
    //so we check to make sure we're looking at the right polygon.
    if(polyID == fCurrentPolygID)
      polyID = fVor->GetRightPolyg(edgeID);

    /*
    if(polyID == fCurrentPolygID) {
      polyID = fVor->GetRightPolyg(edgeID);
      MSG("Bfldtest",Msg::kDebug)
        << "Polygon " << i << ": #" << polyID << "(r)" << endl;
    }
    else {
      MSG("Bfldtest",Msg::kDebug)
        << "Polygon " << i << ": #" << polyID << "(l)" << endl;
    }
    */

    gen3ID = fVor->PolygonToGenerator(polyID);

    //get generator location
    r3.SetXYZ(fVor->GetGeneratorX(gen3ID),fVor->GetGeneratorY(gen3ID),0);

    /*
    MSG("Bfldtest",Msg::kDebug)
      << "g" << i % n << " = (" << r3.X() << "," << r3.Y() << ")" << endl;
    */

    //Test to see if r is in the "wedge" given--that is, the region spanned
    //by r2 - r1 and r3 - r1, with a vertex at r.
    //In most cases, if r is in this wedge, it is also in the triangle with
    //its vertices at r1, r2 and r3.
    if(BfldHandyMath::IsInWedge(r,r1,r2,r3)) {
      if(isFound) {
        MSG("Bfldtest",Msg::kDebug)
          << "r seems to be in several wedges!" << endl;
      }
      else {
        isFound = true;
        break;
      }
    }

    r2 = r3;
    gen2ID = gen3ID;
  }

  if(!isFound) {
    //r was in no wedge. This should be impossible.
    MSG("Bfldtest",Msg::kWarning)
      << "Couldn't pick a wedge for r = (" << r.X() << ","
      << r.Y() << ")!" << endl;
    return b1;
  }

  //Check to see if r is in fact in the Delaunay triangle we've chosen.
  //If it's not, we've come upon the rare (but not impossible) case
  //in which the Delaunay triangle containing r does not have the closest
  //generator, r1, among its vertices. We hope that we need only check
  //the Delaunay triangle sharing an edge with the triangle we have already
  //picked out.
  if(!BfldHandyMath::IsInTri(r,r1,r2,r3)) {
    /*
    MSG("Bfldtest",Msg::kDebug)
      << "Entering secondary search for r = (" << r.X() << ","
      << r.Y() << ")" << endl;
    */

    int homePolyID = fVor->GeneratorToPolygon(gen3ID);
    edgesList
      = fOp->RetrieveEdgesSurrPolygon(homePolyID);
    n = edgesList->NumberOfElements();

    //Look for the edge shared by the polygons of r1 and r3.
    isFound = false;
    for(i = 0;i < n;++i) {
      if(edgesList->At(i) == edgeID) {
        isFound = true;
        break;
      }
    }
    if(!isFound) {
      MSG("Bfldtest",Msg::kWarning)
        << "Couldn't find g1-g3 edge for secondary search." << endl;
      return b1;
    }

    int index = i - 1; //The poor man's modulus. % doesn't understand
    if(index < 0)      //negative numbers.
      index += n;

    int direction = -1;

    //Now, we want to find the edge shared by the polygons of r2 and r3
    //and then look at *that* edge's successor to pick out the last
    //generator we're going to look at for our new Delaunay triangle.

    edgeID = edgesList->At(index);
    if(fVor->GeneratorToPolygon(gen2ID) != fVor->GetLeftPolyg(edgeID)
       && fVor->GeneratorToPolygon(gen2ID) != fVor->GetRightPolyg(edgeID))
      direction = 1;

    index = i + 2 * direction;
    if(index < 0)
      index += n;
    if(index >= n)
      index -= n;

    edgeID = edgesList->At(index);
    polyID = fVor->GetLeftPolyg(edgeID);
    if(polyID == fVor->GeneratorToPolygon(gen3ID))
       polyID = fVor->GetRightPolyg(edgeID);

    r1.SetX(fVor->GetGeneratorX(fVor->PolygonToGenerator(polyID)));
    r1.SetY(fVor->GetGeneratorY(fVor->PolygonToGenerator(polyID)));

    genID = fVor->GetGeneratorID(polyID);
    if(genID == -1) //Point("s"?) at infinity. Zero B field is fine here.
      b1.SetXYZ(0,0,0);
    else
      b1 = *(TVector3 *)(fBField->At(genID));
  }

  //If we STILL haven't managed to pick the correct Delaunay triangle,
  //we're really in trouble. Bail out, and write the troublesome point
  //to a file for use in debugging.
  if(!BfldHandyMath::IsInTri(r,r1,r2,r3)) {
    MSG("Bfldtest",Msg::kWarning)
      << "Couldn't pick out Delaunay triangle for r = (" << r.X() << ","
      << r.Y() << ")!" << endl;

    (*gBfldtestFile) << r.X() << " " << r.Y() << endl;

    //Bail out and return the CNInterpolation
    return CNInterpolation(rNode);
  }

  //In principle, the vertices r1, r2, r3 should be in counterclockwise
  //order.

  /*
  MSG("Bfldtest",Msg::kDebug)
    << "Delaunay triangle about r:" << endl
    << "r1 = generator #" << fVor->PolygonToGenerator(fCurrentPolygID)
    << " = (" << r1.X() << "," << r1.Y() << ")" << endl
    << "r2 = generator #" << gen2ID
    << " = (" << r2.X() << "," << r2.Y() << ")" << endl
    << "r3 = generator #" << gen3ID
    << " = (" << r3.X() << "," << r3.Y() << ")" << endl;
  */

  TVector3 b2; //B field at other two generators
  TVector3 b3;

  //As before, if a generator is out at infinity, we can safely
  //give it zero B field without changing the interpolated B(r);

  if(fVor->GeneratorToPolygon(gen2ID) <= 3)
    b2.SetXYZ(0,0,0);
  else
    b2 = *(TVector3 *)(fBField->At(fVor->GetGeneratorID(gen2ID)));    

  if(fVor->GeneratorToPolygon(gen3ID) <= 3)
    b3.SetXYZ(0,0,0);
  else
    b3 = *(TVector3 *)(fBField->At(fVor->GetGeneratorID(gen3ID)));

  /*
  MSG("Bfldtest",Msg::kDebug)
    << "Field values at Delaunay triangle vertices: " << endl
    << "b1 = (" << b1.X() << "," << b1.Y() << ")" << endl
    << "b2 = (" << b2.X() << "," << b2.Y() << ")" << endl
    << "b3 = (" << b3.X() << "," << b3.Y() << ")" << endl;
  */

  //Interpolate b from b1,b2,b3...

  TVector3 b(0,0,0);

  // Invert the matrix
  //     (r1.X() r1.Y() 1)
  // M = (r2.X() r2.Y() 1)
  //     (r3.X() r3.Y() 1)
  double invDet = 1.0 / BfldHandyMath::Det(r1.X(),r1.Y(),1,
                                           r2.X(),r2.Y(),1,
                                           r3.X(),r3.Y(),1);
  TVector3 mInv1(r2.Y() - r3.Y(),r3.Y() - r1.Y(),r1.Y() - r2.Y());
  mInv1 *= invDet;
  TVector3 mInv2(r3.X() - r2.X(),r1.X() - r3.X(),r2.X() - r1.X());
  mInv2 *= invDet;
  TVector3 mInv3(r2.X() * r3.Y() - r3.X() * r2.Y(),
                 r3.X() * r1.Y() - r3.Y() * r1.X(),
                 r1.X() * r2.Y() - r1.Y() * r2.X());
  mInv3 *= invDet;
  
  //coeffs = M^-1 * bVals

  TVector3 bVals(b1.X(),b2.X(),b3.X());
  TVector3 coeffs(mInv1.Dot(bVals),mInv2.Dot(bVals),mInv3.Dot(bVals));
  b.SetX(r.X() * coeffs.X() + r.Y() * coeffs.Y() + coeffs.Z());

  bVals.SetXYZ(b1.Y(),b2.Y(),b3.Y());
  coeffs.SetXYZ(mInv1.Dot(bVals),mInv2.Dot(bVals),mInv3.Dot(bVals));
  b.SetY(r.X() * coeffs.X() + r.Y() * coeffs.Y() + coeffs.Z());

  return b;
}


void BFLInterpolation::IntBI1stOrder(void)
{
// Implementation of a 1st order bilinear interpolation 
// technique inside a triangle
// For *notation/description of the technique* check my MINOS
// BField web pages

//_______ triangle vertices 
//unused:  Float_t xA = 1;
//unused:  Float_t yA = 1;

//unused:  Float_t xB = 2;
//unused:  Float_t yB = 2;

//unused:  Float_t xC = 3;
//unused:  Float_t yC = 3;

//_______ BField @ triangle vertices
//unused:  Float_t BXa = 4;
//unused:  Float_t BYa = 4;

//unused:  Float_t BXb = 5;
//unused:  Float_t BYb = 5;

//unused:  Float_t BXc = 6;
//unused:  Float_t BYc = 6;

//_______ user's input point 
//unused:  Float_t xU = 10;
//unused:  Float_t yU = 20;

//_______ angle of rotation
//unused:  Float_t theta = TMath::ATan((yB-yA)/(xB-xA));

//_______ rotated coordinates
//unused:  Float_t xa =   xA * TMath::Cos(theta) + yA * TMath::Sin(theta);
//unused:  Float_t ya = - xA * TMath::Sin(theta) + yA * TMath::Cos(theta);

//unused:  Float_t xb =   xB * TMath::Cos(theta) + yB * TMath::Sin(theta);
//unused:  Float_t yb = - xB * TMath::Sin(theta) + yB * TMath::Cos(theta);

//unused:  Float_t xc =   xC * TMath::Cos(theta) + yC * TMath::Sin(theta);
//unused:  Float_t yc = - xC * TMath::Sin(theta) + yC * TMath::Cos(theta);

//unused:  Float_t xu =   xU * TMath::Cos(theta) + yU * TMath::Sin(theta);
//unused:  Float_t yu = - xU * TMath::Sin(theta) + yU * TMath::Cos(theta);

//______ slopes and intermediate point D
//unused:  Float_t dBX_dx = (BXb - BXa) / (xb - xa);
//unused:  Float_t dBY_dx = (BYb - BYa) / (xb - xa);

//unused:  Float_t BXd = BXa + dBX_dx * (xc - xa);  
//unused:  Float_t BYd = BYa + dBY_dx * (xc - xa);  

//unused:  Float_t dBX_dy = (BXc - BXd) / (yc - ya);
//unused:  Float_t dBY_dy = (BYc - BYd) / (yc - ya);

//______ BField @ user's point (2-D Taylor expansion, 1st order terms)
//unused:  Float_t BX = BXa + dBX_dx * (xu - xa) + dBX_dy * (yu - ya);
//unused:  Float_t BY = BYa + dBY_dx * (xu - xa) + dBY_dy * (yu - ya);

   MSG("BFL",Msg::kFatal)
      << "IntBI1stOrder is a completely void method" << endl;
}


Bool_t BFLInterpolation::IsInsideANSYSCell(BFLNode2ACell * TableEntry,
                                           BFLNode * point)
{
  // assuming that the ANSYS cell nodes are provided with a given
  // clockwiseness (-1 or 1), for any given sequential pair of cell
  // nodes (nA,nB) the clockwiseness of any triplet nA,nB,nUser (where
  // nUser is the user's input) will have the same sign IF AND ONLY IF
  // the user input is within the original cell

  Bool_t first = kTRUE;
  Bool_t IsInside = kTRUE;
  Int_t clockwiseness = 1;
  BfldMeshVoronoi * mesh;
  BFLVorOperator vorop;
  Int_t nodes[kNMaxNodesPerACell];

  BfldLoanPool * loanpool = BfldLoanPool::Instance();

  mesh = ((BfldMeshVoronoi*)(loanpool->GetMesh(fGrid,0)));

  TableEntry->GetCellNodes(nodes);

  for(Int_t inode = 0; inode <kNMaxNodesPerACell; inode++) {

     Int_t nodeIdA = nodes[inode];
     Int_t nodeIdB = (inode+1 == kNMaxNodesPerACell) ? nodes[0] : nodes[inode+1];

     BFLNode nA( (mesh->GetGeneratorPosition(nodeIdA)).X(),
                 (mesh->GetGeneratorPosition(nodeIdA)).Y()  );

     BFLNode nB( (mesh->GetGeneratorPosition(nodeIdB)).X(),
                 (mesh->GetGeneratorPosition(nodeIdB)).Y()  );

     Int_t icw = vorop.Clockwise(&nA,&nB,point);

     if(first) {
         clockwiseness = icw;
         first = kFALSE;
     } else {
         if (icw != clockwiseness) IsInside = kFALSE;
     }

  }

  return IsInside;
}

---