#include <BFLInterpolation.h>
Definition at line 34 of file BFLInterpolation.h.
| 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 }
| 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().
| virtual void BFLInterpolation::SetBField | ( | TObjArray * | bfield | ) | [inline, virtual] |
| 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] |
| 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] |
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().
Int_t BFLInterpolation::fNPolygs [private] |
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().
1.4.7