BFLInterpolation Class Reference

#include <BFLInterpolation.h>

List of all members.

Public Member Functions

 BFLInterpolation (TObjArray *bfield, BFLWingedEdge *voronoi)
 ~BFLInterpolation ()
TVector3 GetB (const TVector3 &PosVector, const Bool_t updateSeedPolygon=kFALSE)
virtual BFLCacheGetCache (void)
virtual void SetInterpolant (BfldInterpMethod::InterpMethod_t intp)
virtual void SetGrid (BfldGrid::Grid_t grid)
virtual void SetVoronoi (BFLWingedEdge *vor)
virtual void SetBField (TObjArray *bfield)
virtual void SetCache (BFLCache *cache)
virtual void SetSearchSeed (Int_t seed)
virtual Int_t GetLastPolygon (void)

Private Member Functions

Int_t FormVoronoiCell (BFLNode *point)
Double_t CellArea (TObjArray *vertices)
TVector3 NNInterpolation (BFLNode *point)
TVector3 CNInterpolation (BFLNode *point)
TVector3 BilinearInterpolation (BFLNode *point)
TVector3 PlanarInterpolation (BFLNode *point)
void ResetForNextQuery (void)
Bool_t IsInsideANSYSCell (BFLNode2ACell *TableEntry, BFLNode *point)
void IntBI1stOrder (void)

Private Attributes

BfldInterpMethod::InterpMethod_t fHowToInterpolate
BfldGrid::Grid_t fGrid
Int_t fCurrentPolygID
Int_t fSearchSeed
Int_t fNPolygs
TObjArray * fBField
TObjArray * fGenerators
BFLWingedEdgefVor
BFLVorOperatorfOp
BFLCachefCache
TIntListfPolygons
TObjArray * fVertices


Detailed Description

Definition at line 34 of file BFLInterpolation.h.


Constructor & Destructor Documentation

BFLInterpolation::BFLInterpolation ( TObjArray *  bfield,
BFLWingedEdge voronoi 
)

Definition at line 90 of file BFLInterpolation.cxx.

References fBField, fCache, fNPolygs, fOp, fPolygons, fSearchSeed, fVertices, fVor, gBfldtestFile, BFLWingedEdge::GetNofPolygons(), BfldInterpMethod::kClosest, and SetInterpolant().

00092 { 
00093   fBField = bfield;                  /* BF Vectors */ 
00094   fVor = voronoi;                    /* Voronoi diagram */
00095   fOp = new BFLVorOperator(voronoi); /* Object that operates on fVor */  
00096   fNPolygs = fVor->GetNofPolygons() - 3;
00097 
00098   fCache = new BFLCache();
00099   SetInterpolant(BfldInterpMethod::kClosest);
00100   
00101   fPolygons = new TIntList();
00102   fVertices = new TObjArray(20); 
00103 
00104   fSearchSeed = 1; 
00105 
00106   gBfldtestFile = new ofstream("BfldPlanarInterpErrs.err",ios::app);
00107 } 

BFLInterpolation::~BFLInterpolation (  ) 

Definition at line 110 of file BFLInterpolation.cxx.

References TIntList::Delete(), fBField, fCache, fOp, fPolygons, fVertices, fVor, and gBfldtestFile.

00111 {
00112   fPolygons->Delete();
00113   delete fPolygons;    fPolygons=0;
00114   fVertices->Delete();
00115   delete fVertices;    fVertices=0;
00116   delete fOp;          fOp=0;
00117   delete fVor;         fVor=0;
00118   fBField->Delete();
00119   delete fBField;      fBField=0;
00120   delete fCache;       fCache=0;
00121 
00122   delete gBfldtestFile;
00123 }


Member Function Documentation

TVector3 BFLInterpolation::BilinearInterpolation ( BFLNode point  )  [private]

Definition at line 414 of file BFLInterpolation.cxx.

References fCurrentPolygID, fGrid, BFLVorOperator::FindCurrentPolygon(), BFLAnsysLookup::FindNode(), fOp, fSearchSeed, fVor, BFLWingedEdge::GetAnsysLookupTable(), BfldLoanPool::GetMesh(), BfldLoanPool::Instance(), IsInsideANSYSCell(), kNMaxNodesPerACell, and BFLVorOperator::SetCurrentGen().

Referenced by GetB().

00415 {
00416   // Find the current polygon ID
00417   fOp->SetCurrentGen(point); 
00418   fCurrentPolygID = fOp->FindCurrentPolygon(fSearchSeed); 
00419 
00420   // Get the ANSYS Node that generates this Voronoi Cell
00421   BfldLoanPool * loanpool = BfldLoanPool::Instance();
00422 
00423   Int_t NodeID = ((BfldMeshVoronoi*)(loanpool->GetMesh(fGrid,0)))
00424                                 ->FindANSYSGenerator(fCurrentPolygID);
00425 
00426   // Get the ANSYS Lookup Table
00427   BFLAnsysLookup lookup = fVor->GetAnsysLookupTable();
00428 
00429   // Get a list of all ANSYS Cells having the ANSYS Node
00430   // - we just found - as one of its corners
00431   TObjArray cells;
00432   lookup.FindNode(NodeID,&cells);
00433 
00434   // Select the ANSYS Cell that contains the user provided point
00435   TIter nextentry(&cells);
00436   BFLNode2ACell * TableEntry = 0, * SelectedTableEntry = 0;
00437   while( (TableEntry = (BFLNode2ACell*) nextentry() ) ) {
00438 
00439     Bool_t IsInside = IsInsideANSYSCell(TableEntry,point);
00440     if(IsInside) SelectedTableEntry = (BFLNode2ACell *) (TableEntry->Clone());
00441     if(IsInside) break;
00442   }
00443 
00444   // Get a list of all ANSYS Nodes of this cell
00445   // These are the points Jeff wants to use for 2-linear interpolation
00446   Int_t nodes[kNMaxNodesPerACell];
00447   SelectedTableEntry->GetCellNodes(nodes);
00448 
00449   // Elias should somehow calculate the correct weights
00450   // and return BField
00451 
00452 
00453   delete SelectedTableEntry;
00454 
00455   return TVector3(0,0,0); 
00456 }

Double_t BFLInterpolation::CellArea ( TObjArray *  vertices  )  [private]

Definition at line 379 of file BFLInterpolation.cxx.

References BFLNode::GetX(), and BFLNode::GetY().

Referenced by NNInterpolation().

00380 {
00381 // Getting the area of a cell
00382 
00383   Double_t AreaOfCell = 0;
00384   BFLNode * Vtx;
00385 
00386   TMatrix area = TMatrix(3,3);
00387   for(Int_t i=0; i<3; i++) area(i,2)=1;
00388 
00389   area(0,0) = ((BFLNode *) vertices->At(0))->GetX();
00390   area(0,1) = ((BFLNode *) vertices->At(0))->GetY();
00391   
00392   Int_t ivtx=0;
00393   TIter nextv(vertices);
00394   while( (Vtx = (BFLNode *)nextv()) ) {    
00395     if(ivtx == 1) {
00396        area(1,0) = Vtx->GetX();
00397        area(1,1) = Vtx->GetY();
00398     }
00399     if(ivtx == 2) {
00400        area(2,0) = Vtx->GetX();
00401        area(2,1) = Vtx->GetY(); 
00402        AreaOfCell += 0.5 * area.Determinant();
00403        ivtx = 1;
00404        area(1,0) = area(2,0);
00405        area(1,1) = area(2,1);
00406     } 
00407     ivtx++;
00408   } 
00409 
00410   return TMath::Abs(AreaOfCell); 
00411 }

TVector3 BFLInterpolation::CNInterpolation ( BFLNode point  )  [private]

Definition at line 343 of file BFLInterpolation.cxx.

References fBField, fCurrentPolygID, BFLVorOperator::FindCurrentPolygon(), fOp, fSearchSeed, fVor, BFLWingedEdge::GetGeneratorID(), and BFLVorOperator::SetCurrentGen().

Referenced by GetB(), NNInterpolation(), and PlanarInterpolation().

00344 { 
00345   fOp->SetCurrentGen(point); 
00346   fCurrentPolygID = fOp->FindCurrentPolygon(fSearchSeed); 
00347 
00348   /*
00349   MSG("Bfldtest",Msg::kDebug)
00350     << "r is in polygon #" << fCurrentPolygID
00351     << " with surrounding polygons as follows:" << endl;
00352 
00353   TIntList *edgesList = fOp->RetrieveEdgesSurrPolygon(fCurrentPolygID);
00354   int n = edgesList->NumberOfElements();
00355   int i,edgeID, polyID,genID;
00356   for(i = 0;i < n;++i) {
00357     edgeID = edgesList->At(i);
00358     polyID = fVor->GetLeftPolyg(edgeID);
00359     if(polyID == fCurrentPolygID) {
00360       polyID = fVor->GetRightPolyg(edgeID);
00361     }
00362     MSG("Bfldtest",Msg::kDebug) << "Polygon #" << polyID << endl;
00363   }
00364   */
00365   if(fCurrentPolygID > 3) {
00366      Int_t ID = fVor->GetGeneratorID(fCurrentPolygID); 
00367      return *(TVector3 *) fBField->At(ID);
00368   } else return TVector3(0,0,0);
00369 }

Int_t BFLInterpolation::FormVoronoiCell ( BFLNode point  )  [private]

Definition at line 162 of file BFLInterpolation.cxx.

References TIntList::Add(), fCurrentPolygID, BFLVorOperator::FindCurrentPolygon(), BFLVorOperator::FindNewVtx(), fOp, fPolygons, fSearchSeed, fVertices, BFLVorOperator::GetFirstIntersectEdge(), BFLVorOperator::GetNextIntersectEdge(), kFail, kMaxVtx, BFLVorOperator::MoveToNextPolygon(), BFLVorOperator::SetCurrentGen(), and BFLVorOperator::SetCurrentPolygID().

Referenced by GetB().

00163 {
00164   Int_t EdgeID, StartEdgeID, LastEdgeID, iedge, ivtx = 0; 
00165 //rwh:  TIntList * IntersEdges = new TIntList();
00166  
00167   fOp->SetCurrentGen(point); 
00168   fCurrentPolygID = fOp->FindCurrentPolygon(fSearchSeed); 
00169   
00170   fOp->SetCurrentPolygID(fCurrentPolygID); 
00171   fPolygons->Add(fCurrentPolygID); 
00172   
00173   EdgeID = fOp->GetFirstIntersectEdge();
00174   if(EdgeID == -1) return kFail;
00175   
00176   StartEdgeID = EdgeID; //keep this id to terminate the next loop
00177   LastEdgeID = 0; // removes dependencies from the previous iteration.
00178 
00179   do {       
00180       // store intersect edges for this polygon
00181 //rwh:      IntersEdges->Add(EdgeID); 
00182 
00183       // Get the intersection point -> new vtx.
00184       // this will get deleted in ResetForNextQuery
00185       fVertices->Add(new BFLNode(fOp->FindNewVtx(EdgeID)));
00186                                                                                                                 
00187       // Get the next intersect edge
00188       iedge = fOp->GetNextIntersectEdge(EdgeID);
00189       if (iedge == LastEdgeID ) { 
00190           // It is time to move to the next polygon.
00191 
00192           // Get the next polygon ID.
00193           fCurrentPolygID = fOp->MoveToNextPolygon(fCurrentPolygID,EdgeID);
00194           fOp->SetCurrentPolygID(fCurrentPolygID);
00195           fPolygons->Add(fCurrentPolygID);
00196   
00197           LastEdgeID = EdgeID;
00198           EdgeID = fOp->GetNextIntersectEdge(EdgeID);                                         
00199 
00200        } else {       
00201           LastEdgeID = EdgeID;
00202           EdgeID = iedge;            
00203        }
00204 
00205        ivtx++;
00206        if(ivtx > kMaxVtx) return kFail;
00207        
00208   } while(EdgeID != StartEdgeID); 
00209 
00210   // this will get deleted in ResetForNextQuery
00211   fVertices->Add(new BFLNode(fOp->FindNewVtx(EdgeID))); 
00212   
00213   return 0;
00214 } 

TVector3 BFLInterpolation::GetB ( const TVector3 &  PosVector,
const Bool_t  updateSeedPolygon = kFALSE 
)

Definition at line 126 of file BFLInterpolation.cxx.

References BilinearInterpolation(), CNInterpolation(), fHowToInterpolate, FormVoronoiCell(), GetLastPolygon(), BfldInterpMethod::kBilinear, BfldInterpMethod::kClosest, kFail, BfldInterpMethod::kNatural, BfldInterpMethod::kPlanar, Msg::kWarning, MSG, NNInterpolation(), PlanarInterpolation(), ResetForNextQuery(), and SetSearchSeed().

Referenced by BFLHandler::GetBField().

00128 { 
00129  Int_t istatus;
00130  TVector3 BF(0.,0.,0.);
00131  BFLNode point(PosVector.X(),PosVector.Y()); 
00132  
00133  switch(fHowToInterpolate) {
00134    case(BfldInterpMethod::kNatural): 
00135      istatus = FormVoronoiCell(&point);
00136      if(istatus == kFail) { // for some reason
00137         BF = CNInterpolation(&point); // do the next best thing
00138      } else {
00139         BF = NNInterpolation(&point);
00140         ResetForNextQuery();
00141      }
00142      break;   
00143    case(BfldInterpMethod::kClosest): 
00144      BF = CNInterpolation(&point);
00145      break;
00146    case(BfldInterpMethod::kBilinear):
00147      BF = BilinearInterpolation(&point);
00148      break;
00149  case(BfldInterpMethod::kPlanar):
00150    BF = PlanarInterpolation(&point); //Added by Andrew Goldstone
00151    break;
00152    default:
00153         MSG("BFL",Msg::kWarning) 
00154            << "You did not specified a valid interpolation method" << endl;
00155  }
00156 
00157  if (updateSeedPolygon) SetSearchSeed(GetLastPolygon());
00158  return BF; 
00159 }

virtual BFLCache* BFLInterpolation::GetCache ( void   )  [inline, virtual]

Definition at line 45 of file BFLInterpolation.h.

References fCache.

Referenced by BFLHandler::GetCache().

00045 { return fCache; }

virtual Int_t BFLInterpolation::GetLastPolygon ( void   )  [inline, virtual]

Definition at line 54 of file BFLInterpolation.h.

References fCurrentPolygID.

Referenced by GetB(), and BFLHandler::SetLastPolygAsSeed().

00054 { return fCurrentPolygID; }

void BFLInterpolation::IntBI1stOrder ( void   )  [private]

Definition at line 754 of file BFLInterpolation.cxx.

References Msg::kFatal, and MSG.

00755 {
00756 // Implementation of a 1st order bilinear interpolation 
00757 // technique inside a triangle
00758 // For *notation/description of the technique* check my MINOS
00759 // BField web pages
00760 
00761 //_______ triangle vertices 
00762 //unused:  Float_t xA = 1;
00763 //unused:  Float_t yA = 1;
00764 
00765 //unused:  Float_t xB = 2;
00766 //unused:  Float_t yB = 2;
00767 
00768 //unused:  Float_t xC = 3;
00769 //unused:  Float_t yC = 3;
00770 
00771 //_______ BField @ triangle vertices
00772 //unused:  Float_t BXa = 4;
00773 //unused:  Float_t BYa = 4;
00774 
00775 //unused:  Float_t BXb = 5;
00776 //unused:  Float_t BYb = 5;
00777 
00778 //unused:  Float_t BXc = 6;
00779 //unused:  Float_t BYc = 6;
00780 
00781 //_______ user's input point 
00782 //unused:  Float_t xU = 10;
00783 //unused:  Float_t yU = 20;
00784 
00785 //_______ angle of rotation
00786 //unused:  Float_t theta = TMath::ATan((yB-yA)/(xB-xA));
00787 
00788 //_______ rotated coordinates
00789 //unused:  Float_t xa =   xA * TMath::Cos(theta) + yA * TMath::Sin(theta);
00790 //unused:  Float_t ya = - xA * TMath::Sin(theta) + yA * TMath::Cos(theta);
00791 
00792 //unused:  Float_t xb =   xB * TMath::Cos(theta) + yB * TMath::Sin(theta);
00793 //unused:  Float_t yb = - xB * TMath::Sin(theta) + yB * TMath::Cos(theta);
00794 
00795 //unused:  Float_t xc =   xC * TMath::Cos(theta) + yC * TMath::Sin(theta);
00796 //unused:  Float_t yc = - xC * TMath::Sin(theta) + yC * TMath::Cos(theta);
00797 
00798 //unused:  Float_t xu =   xU * TMath::Cos(theta) + yU * TMath::Sin(theta);
00799 //unused:  Float_t yu = - xU * TMath::Sin(theta) + yU * TMath::Cos(theta);
00800 
00801 //______ slopes and intermediate point D
00802 //unused:  Float_t dBX_dx = (BXb - BXa) / (xb - xa);
00803 //unused:  Float_t dBY_dx = (BYb - BYa) / (xb - xa);
00804 
00805 //unused:  Float_t BXd = BXa + dBX_dx * (xc - xa);  
00806 //unused:  Float_t BYd = BYa + dBY_dx * (xc - xa);  
00807 
00808 //unused:  Float_t dBX_dy = (BXc - BXd) / (yc - ya);
00809 //unused:  Float_t dBY_dy = (BYc - BYd) / (yc - ya);
00810 
00811 //______ BField @ user's point (2-D Taylor expansion, 1st order terms)
00812 //unused:  Float_t BX = BXa + dBX_dx * (xu - xa) + dBX_dy * (yu - ya);
00813 //unused:  Float_t BY = BYa + dBY_dx * (xu - xa) + dBY_dy * (yu - ya);
00814 
00815    MSG("BFL",Msg::kFatal)
00816       << "IntBI1stOrder is a completely void method" << endl;
00817 }

Bool_t BFLInterpolation::IsInsideANSYSCell ( BFLNode2ACell TableEntry,
BFLNode point 
) [private]

Definition at line 821 of file BFLInterpolation.cxx.

References BFLVorOperator::Clockwise(), fGrid, BFLNode2ACell::GetCellNodes(), BfldMeshVoronoi::GetGeneratorPosition(), BfldLoanPool::GetMesh(), BfldLoanPool::Instance(), and kNMaxNodesPerACell.

Referenced by BilinearInterpolation().

00823 {
00824   // assuming that the ANSYS cell nodes are provided with a given
00825   // clockwiseness (-1 or 1), for any given sequential pair of cell
00826   // nodes (nA,nB) the clockwiseness of any triplet nA,nB,nUser (where
00827   // nUser is the user's input) will have the same sign IF AND ONLY IF
00828   // the user input is within the original cell
00829 
00830   Bool_t first = kTRUE;
00831   Bool_t IsInside = kTRUE;
00832   Int_t clockwiseness = 1;
00833   BfldMeshVoronoi * mesh;
00834   BFLVorOperator vorop;
00835   Int_t nodes[kNMaxNodesPerACell];
00836 
00837   BfldLoanPool * loanpool = BfldLoanPool::Instance();
00838 
00839   mesh = ((BfldMeshVoronoi*)(loanpool->GetMesh(fGrid,0)));
00840 
00841   TableEntry->GetCellNodes(nodes);
00842 
00843   for(Int_t inode = 0; inode <kNMaxNodesPerACell; inode++) {
00844 
00845      Int_t nodeIdA = nodes[inode];
00846      Int_t nodeIdB = (inode+1 == kNMaxNodesPerACell) ? nodes[0] : nodes[inode+1];
00847 
00848      BFLNode nA( (mesh->GetGeneratorPosition(nodeIdA)).X(),
00849                  (mesh->GetGeneratorPosition(nodeIdA)).Y()  );
00850 
00851      BFLNode nB( (mesh->GetGeneratorPosition(nodeIdB)).X(),
00852                  (mesh->GetGeneratorPosition(nodeIdB)).Y()  );
00853 
00854      Int_t icw = vorop.Clockwise(&nA,&nB,point);
00855 
00856      if(first) {
00857          clockwiseness = icw;
00858          first = kFALSE;
00859      } else {
00860          if (icw != clockwiseness) IsInside = kFALSE;
00861      }
00862 
00863   }
00864 
00865   return IsInside;
00866 }

TVector3 BFLInterpolation::NNInterpolation ( BFLNode point  )  [private]

Definition at line 217 of file BFLInterpolation.cxx.

References TIntList::At(), CellArea(), BFLVorOperator::Clockwise(), CNInterpolation(), TIntList::Delete(), fBField, fOp, fPolygons, fVertices, fVor, BFLWingedEdge::GetGeneratorID(), BFLNode::GetNodeID(), BFLWingedEdge::GetVtx(), BFLNode::GetX(), BFLNode::GetY(), Msg::kWarning, MSG, TIntList::NumberOfElements(), BFLVorOperator::ResetVtxCache(), and BFLVorOperator::RetrieveVtxSurrPolygon().

Referenced by GetB().

00218 {
00219   Int_t NatNeighbor = 0, ID = 0;
00220   Double_t AreaSum = 0, NeighborArea = 0, bx = 0, by = 0;
00221 
00222   // total area of new cell
00223   Double_t TotalArea = CellArea(fVertices); 
00224   
00225   // compute sub-areas
00226   for(Int_t inn=0; inn<fPolygons->NumberOfElements()-1; inn++) {
00227      NatNeighbor = fPolygons->At(inn);
00228 
00229      //rwh:     TObjArray * Area = new TObjArray(30); // vertices of sub-cell
00230      //rwh:  increase this size after seeing index out of bounds in AddAt()
00231      TObjArray * Area = new TObjArray(55); // vertices of sub-cell
00232      BFLNode * FirstVtx = (BFLNode *) fVertices->At(inn);
00233      BFLNode * SecondVtx = (BFLNode *) fVertices->At(inn+1);
00234 
00235      if (!FirstVtx || !SecondVtx) {
00236         //rwh: protect against SEGV
00237         MSG("BFL",Msg::kWarning) 
00238            << "failed fVertices->At(inn) " << inn 
00239            << " FirstVtx " << FirstVtx << " SecondVtx " << SecondVtx
00240            << endl
00241            << "   fPolygons->NumberOfElements " 
00242            << fPolygons->NumberOfElements() << endl;
00243         break;
00244      }
00245 
00246      Int_t Clockwiseness = fOp->Clockwise(FirstVtx,SecondVtx,point);
00247 
00248      // select other polygon's vertices on the same side with 'point' 
00249      TIntList * VtxAroundNeighbor = fOp->RetrieveVtxSurrPolygon(NatNeighbor);
00250 
00251      Int_t ivtxslot=0;
00252      for(Int_t ivtx=0; ivtx<VtxAroundNeighbor->NumberOfElements(); ivtx++) {
00253         Int_t VtxID = VtxAroundNeighbor->At(ivtx);
00254         BFLNode * Vtx = fVor->GetVtx(VtxID);       
00255         
00256         if(fOp->Clockwise(FirstVtx,SecondVtx,Vtx) == Clockwiseness) {
00257            Area->AddAt(Vtx,ivtxslot);
00258            ivtxslot++;
00259         } else {
00260            // fVor->GetVtx returns "new" BFLNode
00261            // if we're not saving it (to be deleted later) then toss it now
00262            delete Vtx;
00263         }
00264      }      
00265 
00266      VtxAroundNeighbor->Delete();
00267      delete VtxAroundNeighbor;  VtxAroundNeighbor=0;
00268 
00269      // ivtxslot = 0 means that the subcell was not formed
00270      if(ivtxslot == 0) {
00271         MSG("BFL",Msg::kWarning)
00272            << "BFLInterpolation::NNInterpolation --- subcell not formed" 
00273            << endl        
00274            << "    BFLNode 'point' " << point->GetNodeID() 
00275            << " x= " << point->GetX() << " y= " << point->GetY() << endl;
00276         return CNInterpolation(point);
00277      }
00278 
00279      // Make copies of the BFLNode's FirstVtx and SecondVtx
00280      // so that "Area" can _own_ them along with the entries
00281      // put in from a fVor->GetVtx call (which creates BFLNodes)
00282      BFLNode * FirstVtxCopy = new BFLNode(*FirstVtx);
00283      BFLNode * SecondVtxCopy = new BFLNode(*SecondVtx);
00284 
00285      if(ivtxslot == 1) {        
00286         Area->AddAt(FirstVtxCopy,1);
00287         Area->AddAt(SecondVtxCopy,2);
00288      } else {
00289         Int_t vcw = fOp->Clockwise((BFLNode *)Area->At(0),
00290                                    (BFLNode *)Area->At(1),FirstVtx);
00291 
00292         if(fOp->Clockwise((BFLNode *)Area->At(1),FirstVtx,SecondVtx) == vcw) {
00293            Area->AddAt(FirstVtxCopy,ivtxslot);
00294            Area->AddAt(SecondVtxCopy,ivtxslot+1);
00295         } else {
00296            Area->AddAt(SecondVtxCopy,ivtxslot);
00297            Area->AddAt(FirstVtxCopy,ivtxslot+1);
00298         }
00299      } 
00300      
00301      NeighborArea = CellArea(Area);
00302      AreaSum += NeighborArea;
00303 
00304      ID = fVor->GetGeneratorID(NatNeighbor);
00305      if (ID<0) {
00306         //rwh: protect against SEGV
00307         MSG("BFL",Msg::kWarning)
00308            << "BFLInterpolation::NNInterpolation bad ID: " << ID << endl;
00309      } else {
00310      bx += NeighborArea * ((TVector2 *) fBField->At(ID))->X();
00311      by += NeighborArea * ((TVector2 *) fBField->At(ID))->Y();
00312      }
00313 
00314      // Area owns some objects [created by fVor->GetVtx(VtxID)]
00315      Area->Delete();
00316      delete Area;  Area=0;
00317   }
00318 
00319   // For the last neighbor we use the fact that Sum(weights) = 1
00320   NatNeighbor = fPolygons->At(fPolygons->NumberOfElements() - 1);
00321   NeighborArea = TotalArea - AreaSum; 
00322   
00323   ID = fVor->GetGeneratorID(NatNeighbor);
00324 
00325   if (ID<0) {
00326      //rwh: protect against SEGV
00327      MSG("BFL",Msg::kWarning) 
00328         << "BFLInterpolation::NNInterpolation bad ID: " << ID << endl;
00329   } else {
00330   bx += NeighborArea * ((TVector2 *) fBField->At(ID))->X();
00331   by += NeighborArea * ((TVector2 *) fBField->At(ID))->Y();
00332   }
00333 
00334   bx /= TotalArea;
00335   by /= TotalArea;
00336 
00337   fOp->ResetVtxCache();
00338 
00339   return TVector3(bx,by,0); 
00340 }

TVector3 BFLInterpolation::PlanarInterpolation ( BFLNode point  )  [private]

Definition at line 460 of file BFLInterpolation.cxx.

References TIntList::At(), CNInterpolation(), BfldHandyMath::Det(), fBField, fCurrentPolygID, BFLVorOperator::FindCurrentPolygon(), fOp, fSearchSeed, fVor, BFLWingedEdge::GeneratorToPolygon(), BFLWingedEdge::GetGeneratorID(), BFLWingedEdge::GetGeneratorX(), BFLWingedEdge::GetGeneratorY(), BFLWingedEdge::GetLeftPolyg(), BFLWingedEdge::GetRightPolyg(), BFLNode::GetX(), BFLNode::GetY(), BfldHandyMath::IsInTri(), BfldHandyMath::IsInWedge(), Msg::kDebug, Msg::kWarning, MSG, TIntList::NumberOfElements(), BFLWingedEdge::PolygonToGenerator(), BFLVorOperator::RetrieveEdgesSurrPolygon(), and BFLVorOperator::SetCurrentGen().

Referenced by GetB().

00460                                                              {
00461 
00462   fOp->SetCurrentGen(rNode);
00463   fCurrentPolygID = fOp->FindCurrentPolygon(fSearchSeed);
00464 
00465   if(fCurrentPolygID <= 3) {
00466     //This means the nearest generator is the point at infinity,
00467     //so we can safely return 0. See BFLInterpolation::CNInterpolation()
00468    return TVector3(0,0,0);
00469   }
00470 
00471   //Get the coordinates of the test point r
00472   TVector3 r(rNode->GetX(),rNode->GetY(),0);
00473 
00474   //Get generator coordinates
00475   TVector3 r1(fVor->GetGeneratorX(fVor->PolygonToGenerator(fCurrentPolygID)),
00476               fVor->GetGeneratorY(fVor->PolygonToGenerator(fCurrentPolygID)),
00477               0);
00478 
00479   //Get B field at closest generator. Note that we must map the ID of the
00480   //Voronoi cell into the ID of the B field for that cell's generator.
00481   //Note further that this is NOT the generator ID used by BFLWingedEdge.
00482   //In particular, BFLWingedEdge::GetGenerator[X,Y] takes the other ID,
00483   //which is, in fact, the same as the ID of the cell containing that
00484   //generator.
00485   int genID = fVor->GetGeneratorID(fCurrentPolygID);
00486   TVector3 b1 = *(TVector3 *)(fBField->At(genID));
00487 
00488   const double kDistanceThresholdSquared = 10.0; //In mm^2, alas
00489 
00490   if((r1-r).Mag2() < kDistanceThresholdSquared) {
00491     //r is essentially at a generator. In this case, it is not in the
00492     //interior of any Delaunay triangle, and we save ourselves the
00493     //trouble of a pointless interpolation.
00494     
00495     return b1;
00496   }
00497 
00498   //iterate through the generators of the neighboring cells until we
00499   //find the Delaunay triangle containing r
00500 
00501   TIntList *edgesList = fOp->RetrieveEdgesSurrPolygon(fCurrentPolygID);
00502 
00503   TVector3 r2,r3;
00504   int polyID,edgeID,gen2ID,gen3ID=-1;
00505   int i;
00506 
00507   //Apparently we get a list of surrounding edges in counterclockwise
00508   //order; see BFLVorOperator::RetrieveEdgesSurrPolygon(),
00509   //but apparently the edges do not have a common direction, so
00510   //we check which way they're directed to pick the right polygon.
00511 
00512   int n = edgesList->NumberOfElements();
00513   bool isFound = false;
00514 
00515   edgeID = edgesList->At(0);
00516   polyID = fVor->GetLeftPolyg(edgeID);
00517   if(polyID == fCurrentPolygID) {
00518     polyID = fVor->GetRightPolyg(edgeID);
00519   }
00520 
00521   //Take note: as mentioned in the discussion of generator ID's above,
00522   //the ID of the polygon and the ID of its generator are the same
00523   //for most calls to BFLWingedEdge. In fact,
00524   //BFLWingedEdge::PolygonToGenerator just returns its argument.
00525   gen2ID = fVor->PolygonToGenerator(polyID);
00526   r2.SetXYZ(fVor->GetGeneratorX(gen2ID),fVor->GetGeneratorY(gen2ID),0);
00527 
00528   /*
00529   MSG("Bfldtest",Msg::kDebug)
00530     << "Scanning generators around r = (" << r.X() << "," << r.Y()
00531     << ") in polygon #" << fCurrentPolygID << endl;
00532 
00533   MSG("Bfldtest",Msg::kDebug)
00534     << "ID's of cell edges: " << endl;
00535   for(i = 0;i < n;i++) {
00536     MSG("Bfldtest",Msg::kDebug) << "e" << i << ": #" << edgesList->At(i)
00537                                << endl;
00538   } 
00539   */
00540 
00541   for(i = 1;i <= n;++i) {
00542     edgeID = edgesList->At(i % n);
00543     polyID = fVor->GetLeftPolyg(edgeID);
00544 
00545     //As above, we don't know which way the edge is directed,
00546     //so we check to make sure we're looking at the right polygon.
00547     if(polyID == fCurrentPolygID)
00548       polyID = fVor->GetRightPolyg(edgeID);
00549 
00550     /*
00551     if(polyID == fCurrentPolygID) {
00552       polyID = fVor->GetRightPolyg(edgeID);
00553       MSG("Bfldtest",Msg::kDebug)
00554         << "Polygon " << i << ": #" << polyID << "(r)" << endl;
00555     }
00556     else {
00557       MSG("Bfldtest",Msg::kDebug)
00558         << "Polygon " << i << ": #" << polyID << "(l)" << endl;
00559     }
00560     */
00561 
00562     gen3ID = fVor->PolygonToGenerator(polyID);
00563 
00564     //get generator location
00565     r3.SetXYZ(fVor->GetGeneratorX(gen3ID),fVor->GetGeneratorY(gen3ID),0);
00566 
00567     /*
00568     MSG("Bfldtest",Msg::kDebug)
00569       << "g" << i % n << " = (" << r3.X() << "," << r3.Y() << ")" << endl;
00570     */
00571 
00572     //Test to see if r is in the "wedge" given--that is, the region spanned
00573     //by r2 - r1 and r3 - r1, with a vertex at r.
00574     //In most cases, if r is in this wedge, it is also in the triangle with
00575     //its vertices at r1, r2 and r3.
00576     if(BfldHandyMath::IsInWedge(r,r1,r2,r3)) {
00577       if(isFound) {
00578         MSG("Bfldtest",Msg::kDebug)
00579           << "r seems to be in several wedges!" << endl;
00580       }
00581       else {
00582         isFound = true;
00583         break;
00584       }
00585     }
00586 
00587     r2 = r3;
00588     gen2ID = gen3ID;
00589   }
00590 
00591   if(!isFound) {
00592     //r was in no wedge. This should be impossible.
00593     MSG("Bfldtest",Msg::kWarning)
00594       << "Couldn't pick a wedge for r = (" << r.X() << ","
00595       << r.Y() << ")!" << endl;
00596     return b1;
00597   }
00598 
00599   //Check to see if r is in fact in the Delaunay triangle we've chosen.
00600   //If it's not, we've come upon the rare (but not impossible) case
00601   //in which the Delaunay triangle containing r does not have the closest
00602   //generator, r1, among its vertices. We hope that we need only check
00603   //the Delaunay triangle sharing an edge with the triangle we have already
00604   //picked out.
00605   if(!BfldHandyMath::IsInTri(r,r1,r2,r3)) {
00606     /*
00607     MSG("Bfldtest",Msg::kDebug)
00608       << "Entering secondary search for r = (" << r.X() << ","
00609       << r.Y() << ")" << endl;
00610     */
00611 
00612     int homePolyID = fVor->GeneratorToPolygon(gen3ID);
00613     edgesList
00614       = fOp->RetrieveEdgesSurrPolygon(homePolyID);
00615     n = edgesList->NumberOfElements();
00616 
00617     //Look for the edge shared by the polygons of r1 and r3.
00618     isFound = false;
00619     for(i = 0;i < n;++i) {
00620       if(edgesList->At(i) == edgeID) {
00621         isFound = true;
00622         break;
00623       }
00624     }
00625     if(!isFound) {
00626       MSG("Bfldtest",Msg::kWarning)
00627         << "Couldn't find g1-g3 edge for secondary search." << endl;
00628       return b1;
00629     }
00630 
00631     int index = i - 1; //The poor man's modulus. % doesn't understand
00632     if(index < 0)      //negative numbers.
00633       index += n;
00634 
00635     int direction = -1;
00636 
00637     //Now, we want to find the edge shared by the polygons of r2 and r3
00638     //and then look at *that* edge's successor to pick out the last
00639     //generator we're going to look at for our new Delaunay triangle.
00640 
00641     edgeID = edgesList->At(index);
00642     if(fVor->GeneratorToPolygon(gen2ID) != fVor->GetLeftPolyg(edgeID)
00643        && fVor->GeneratorToPolygon(gen2ID) != fVor->GetRightPolyg(edgeID))
00644       direction = 1;
00645 
00646     index = i + 2 * direction;
00647     if(index < 0)
00648       index += n;
00649     if(index >= n)
00650       index -= n;
00651 
00652     edgeID = edgesList->At(index);
00653     polyID = fVor->GetLeftPolyg(edgeID);
00654     if(polyID == fVor->GeneratorToPolygon(gen3ID))
00655        polyID = fVor->GetRightPolyg(edgeID);
00656 
00657     r1.SetX(fVor->GetGeneratorX(fVor->PolygonToGenerator(polyID)));
00658     r1.SetY(fVor->GetGeneratorY(fVor->PolygonToGenerator(polyID)));
00659 
00660     genID = fVor->GetGeneratorID(polyID);
00661     if(genID == -1) //Point("s"?) at infinity. Zero B field is fine here.
00662       b1.SetXYZ(0,0,0);
00663     else
00664       b1 = *(TVector3 *)(fBField->At(genID));
00665   }
00666 
00667   //If we STILL haven't managed to pick the correct Delaunay triangle,
00668   //we're really in trouble. Bail out, and write the troublesome point
00669   //to a file for use in debugging.
00670   if(!BfldHandyMath::IsInTri(r,r1,r2,r3)) {
00671     MSG("Bfldtest",Msg::kWarning)
00672       << "Couldn't pick out Delaunay triangle for r = (" << r.X() << ","
00673       << r.Y() << ")!" << endl;
00674 
00675     (*gBfldtestFile) << r.X() << " " << r.Y() << endl;
00676 
00677     //Bail out and return the CNInterpolation
00678     return CNInterpolation(rNode);
00679   }
00680 
00681   //In principle, the vertices r1, r2, r3 should be in counterclockwise
00682   //order.
00683 
00684   /*
00685   MSG("Bfldtest",Msg::kDebug)
00686     << "Delaunay triangle about r:" << endl
00687     << "r1 = generator #" << fVor->PolygonToGenerator(fCurrentPolygID)
00688     << " = (" << r1.X() << "," << r1.Y() << ")" << endl
00689     << "r2 = generator #" << gen2ID
00690     << " = (" << r2.X() << "," << r2.Y() << ")" << endl
00691     << "r3 = generator #" << gen3ID
00692     << " = (" << r3.X() << "," << r3.Y() << ")" << endl;
00693   */
00694 
00695   TVector3 b2; //B field at other two generators
00696   TVector3 b3;
00697 
00698   //As before, if a generator is out at infinity, we can safely
00699   //give it zero B field without changing the interpolated B(r);
00700 
00701   if(fVor->GeneratorToPolygon(gen2ID) <= 3)
00702     b2.SetXYZ(0,0,0);
00703   else
00704     b2 = *(TVector3 *)(fBField->At(fVor->GetGeneratorID(gen2ID)));    
00705 
00706   if(fVor->GeneratorToPolygon(gen3ID) <= 3)
00707     b3.SetXYZ(0,0,0);
00708   else
00709     b3 = *(TVector3 *)(fBField->At(fVor->GetGeneratorID(gen3ID)));
00710 
00711   /*
00712   MSG("Bfldtest",Msg::kDebug)
00713     << "Field values at Delaunay triangle vertices: " << endl
00714     << "b1 = (" << b1.X() << "," << b1.Y() << ")" << endl
00715     << "b2 = (" << b2.X() << "," << b2.Y() << ")" << endl
00716     << "b3 = (" << b3.X() << "," << b3.Y() << ")" << endl;
00717   */
00718 
00719   //Interpolate b from b1,b2,b3...
00720 
00721   TVector3 b(0,0,0);
00722 
00723   // Invert the matrix
00724   //     (r1.X() r1.Y() 1)
00725   // M = (r2.X() r2.Y() 1)
00726   //     (r3.X() r3.Y() 1)
00727   double invDet = 1.0 / BfldHandyMath::Det(r1.X(),r1.Y(),1,
00728                                            r2.X(),r2.Y(),1,
00729                                            r3.X(),r3.Y(),1);
00730   TVector3 mInv1(r2.Y() - r3.Y(),r3.Y() - r1.Y(),r1.Y() - r2.Y());
00731   mInv1 *= invDet;
00732   TVector3 mInv2(r3.X() - r2.X(),r1.X() - r3.X(),r2.X() - r1.X());
00733   mInv2 *= invDet;
00734   TVector3 mInv3(r2.X() * r3.Y() - r3.X() * r2.Y(),
00735                  r3.X() * r1.Y() - r3.Y() * r1.X(),
00736                  r1.X() * r2.Y() - r1.Y() * r2.X());
00737   mInv3 *= invDet;
00738   
00739   //coeffs = M^-1 * bVals
00740 
00741   TVector3 bVals(b1.X(),b2.X(),b3.X());
00742   TVector3 coeffs(mInv1.Dot(bVals),mInv2.Dot(bVals),mInv3.Dot(bVals));
00743   b.SetX(r.X() * coeffs.X() + r.Y() * coeffs.Y() + coeffs.Z());
00744 
00745   bVals.SetXYZ(b1.Y(),b2.Y(),b3.Y());
00746   coeffs.SetXYZ(mInv1.Dot(bVals),mInv2.Dot(bVals),mInv3.Dot(bVals));
00747   b.SetY(r.X() * coeffs.X() + r.Y() * coeffs.Y() + coeffs.Z());
00748 
00749   return b;
00750 }

void BFLInterpolation::ResetForNextQuery ( void   )  [private]

Definition at line 372 of file BFLInterpolation.cxx.

References TIntList::Delete(), fPolygons, and fVertices.

Referenced by GetB().

00373 {
00374   fVertices->Delete(); 
00375   fPolygons->Delete(); 
00376 }  

virtual void BFLInterpolation::SetBField ( TObjArray *  bfield  )  [inline, virtual]

Definition at line 51 of file BFLInterpolation.h.

References fBField.

00051 { fBField = bfield; }

virtual void BFLInterpolation::SetCache ( BFLCache cache  )  [inline, virtual]

Definition at line 52 of file BFLInterpolation.h.

References fCache.

Referenced by BFLHandler::SetCache().

00052 { fCache = cache; }

virtual void BFLInterpolation::SetGrid ( BfldGrid::Grid_t  grid  )  [inline, virtual]

Definition at line 49 of file BFLInterpolation.h.

References fGrid.

00049 { fGrid = grid;  }

virtual void BFLInterpolation::SetInterpolant ( BfldInterpMethod::InterpMethod_t  intp  )  [inline, virtual]

Definition at line 46 of file BFLInterpolation.h.

References fHowToInterpolate.

Referenced by BFLInterpolation(), and BFLHandler::SetInterpolant().

00046                                                                    {
00047       fHowToInterpolate = intp; 
00048   }

virtual void BFLInterpolation::SetSearchSeed ( Int_t  seed  )  [inline, virtual]

Definition at line 53 of file BFLInterpolation.h.

References fSearchSeed.

Referenced by GetB(), and BFLHandler::SetLastPolygAsSeed().

00053 { fSearchSeed = seed; }

virtual void BFLInterpolation::SetVoronoi ( BFLWingedEdge vor  )  [inline, virtual]

Definition at line 50 of file BFLInterpolation.h.

References fVor.

00050 { fVor = vor; }


Member Data Documentation

TObjArray* BFLInterpolation::fBField [private]

Definition at line 74 of file BFLInterpolation.h.

Referenced by BFLInterpolation(), CNInterpolation(), NNInterpolation(), PlanarInterpolation(), SetBField(), and ~BFLInterpolation().

BFLCache* BFLInterpolation::fCache [private]

Definition at line 78 of file BFLInterpolation.h.

Referenced by BFLInterpolation(), GetCache(), SetCache(), and ~BFLInterpolation().

Int_t BFLInterpolation::fCurrentPolygID [private]

Definition at line 71 of file BFLInterpolation.h.

Referenced by BilinearInterpolation(), CNInterpolation(), FormVoronoiCell(), GetLastPolygon(), and PlanarInterpolation().

TObjArray* BFLInterpolation::fGenerators [private]

Definition at line 75 of file BFLInterpolation.h.

BfldGrid::Grid_t BFLInterpolation::fGrid [private]

Definition at line 70 of file BFLInterpolation.h.

Referenced by BilinearInterpolation(), IsInsideANSYSCell(), and SetGrid().

BfldInterpMethod::InterpMethod_t BFLInterpolation::fHowToInterpolate [private]

Definition at line 69 of file BFLInterpolation.h.

Referenced by GetB(), and SetInterpolant().

Int_t BFLInterpolation::fNPolygs [private]

Definition at line 73 of file BFLInterpolation.h.

Referenced by BFLInterpolation().

BFLVorOperator* BFLInterpolation::fOp [private]

Definition at line 77 of file BFLInterpolation.h.

Referenced by BFLInterpolation(), BilinearInterpolation(), CNInterpolation(), FormVoronoiCell(), NNInterpolation(), PlanarInterpolation(), and ~BFLInterpolation().

TIntList* BFLInterpolation::fPolygons [private]

Definition at line 79 of file BFLInterpolation.h.

Referenced by BFLInterpolation(), FormVoronoiCell(), NNInterpolation(), ResetForNextQuery(), and ~BFLInterpolation().

Int_t BFLInterpolation::fSearchSeed [private]

Definition at line 72 of file BFLInterpolation.h.

Referenced by BFLInterpolation(), BilinearInterpolation(), CNInterpolation(), FormVoronoiCell(), PlanarInterpolation(), and SetSearchSeed().

TObjArray* BFLInterpolation::fVertices [private]

Definition at line 80 of file BFLInterpolation.h.

Referenced by BFLInterpolation(), FormVoronoiCell(), NNInterpolation(), ResetForNextQuery(), and ~BFLInterpolation().

BFLWingedEdge* BFLInterpolation::fVor [private]

Definition at line 76 of file BFLInterpolation.h.

Referenced by BFLInterpolation(), BilinearInterpolation(), CNInterpolation(), NNInterpolation(), PlanarInterpolation(), SetVoronoi(), and ~BFLInterpolation().


The documentation for this class was generated from the following files:
Generated on Fri Oct 10 22:45:21 2014 for loon by  doxygen 1.4.7