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, n, 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 78 of file BFLInterpolation.h.

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

TObjArray* BFLInterpolation::fGenerators [private]

Definition at line 75 of file BFLInterpolation.h.

Definition at line 70 of file BFLInterpolation.h.

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

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().

TObjArray* BFLInterpolation::fVertices [private]

The documentation for this class was generated from the following files:

Generated on 13 Sep 2017 for loon by  doxygen 1.6.1