#include <BFLInterpolation.h>
|
||||||||||||
|
Definition at line 90 of file BFLInterpolation.cxx. References fBField, fCache, fNPolygs, fOp, fPolygons, fSearchSeed, fVertices, fVor, gBfldtestFile, BFLWingedEdge::GetNofPolygons(), 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 }
|
|
|
Definition at line 110 of file BFLInterpolation.cxx. References TIntList::Delete(), fBField, fCache, fOp, fPolygons, fVertices, and fVor. 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 }
|
|
|
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(), 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 }
|
|
|
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 }
|
|
|
Definition at line 343 of file BFLInterpolation.cxx. References fBField, fCurrentPolygID, BFLVorOperator::FindCurrentPolygon(), fOp, fSearchSeed, fVor, BFLWingedEdge::GetGeneratorID(), ID, 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 }
|
|
|
Definition at line 162 of file BFLInterpolation.cxx. References TIntList::Add(), fCurrentPolygID, BFLVorOperator::FindCurrentPolygon(), BFLVorOperator::FindNewVtx(), fOp, fPolygons, fSearchSeed, fVertices, BFLVorOperator::GetFirstIntersectEdge(), BFLVorOperator::GetNextIntersectEdge(), 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 }
|
|
||||||||||||
|
Definition at line 126 of file BFLInterpolation.cxx. References BilinearInterpolation(), CNInterpolation(), FormVoronoiCell(), GetLastPolygon(), 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 }
|
|
|
Definition at line 45 of file BFLInterpolation.h. Referenced by BFLHandler::GetCache(). 00045 { return fCache; }
|
|
|
Definition at line 54 of file BFLInterpolation.h. Referenced by GetB(), and BFLHandler::SetLastPolygAsSeed(). 00054 { return fCurrentPolygID; }
|
|
|
Definition at line 754 of file BFLInterpolation.cxx. References 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 }
|
|
||||||||||||
|
Definition at line 821 of file BFLInterpolation.cxx. References BFLVorOperator::Clockwise(), fGrid, BFLNode2ACell::GetCellNodes(), BfldMeshVoronoi::GetGeneratorPosition(), BfldLoanPool::GetMesh(), and BfldLoanPool::Instance(). 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 }
|
|
|
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(), ID, 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 }
|
|
|
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, TIntList::NumberOfElements(), BFLWingedEdge::PolygonToGenerator(), BFLVorOperator::RetrieveEdgesSurrPolygon(), BFLVorOperator::SetCurrentGen(), BFLWingedEdge::SetX(), and BFLWingedEdge::SetY(). 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 }
|
|
|
Definition at line 372 of file BFLInterpolation.cxx. References TIntList::Delete(), fPolygons, and fVertices. Referenced by GetB().
|
|
|
Definition at line 51 of file BFLInterpolation.h. References fBField. 00051 { fBField = bfield; }
|
|
|
Definition at line 52 of file BFLInterpolation.h. References fCache. Referenced by BFLHandler::SetCache(). 00052 { fCache = cache; }
|
|
|
Definition at line 49 of file BFLInterpolation.h. References fGrid. 00049 { fGrid = grid; }
|
|
|
Definition at line 46 of file BFLInterpolation.h. References fHowToInterpolate. Referenced by BFLInterpolation(), and BFLHandler::SetInterpolant(). 00046 {
00047 fHowToInterpolate = intp;
00048 }
|
|
|
Definition at line 53 of file BFLInterpolation.h. References fSearchSeed. Referenced by GetB(), and BFLHandler::SetLastPolygAsSeed(). 00053 { fSearchSeed = seed; }
|
|
|
Definition at line 50 of file BFLInterpolation.h. References fVor. 00050 { fVor = vor; }
|
|
|
Definition at line 74 of file BFLInterpolation.h. Referenced by BFLInterpolation(), CNInterpolation(), NNInterpolation(), PlanarInterpolation(), SetBField(), and ~BFLInterpolation(). |
|
|
Definition at line 78 of file BFLInterpolation.h. Referenced by BFLInterpolation(), SetCache(), and ~BFLInterpolation(). |
|
|
Definition at line 71 of file BFLInterpolation.h. Referenced by BilinearInterpolation(), CNInterpolation(), FormVoronoiCell(), and PlanarInterpolation(). |
|
|
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 SetInterpolant(). |
|
|
Definition at line 73 of file BFLInterpolation.h. Referenced by BFLInterpolation(). |
|
|
Definition at line 77 of file BFLInterpolation.h. Referenced by BFLInterpolation(), BilinearInterpolation(), CNInterpolation(), FormVoronoiCell(), NNInterpolation(), PlanarInterpolation(), and ~BFLInterpolation(). |
|
|
Definition at line 79 of file BFLInterpolation.h. Referenced by BFLInterpolation(), FormVoronoiCell(), NNInterpolation(), ResetForNextQuery(), and ~BFLInterpolation(). |
|
|
Definition at line 72 of file BFLInterpolation.h. Referenced by BFLInterpolation(), BilinearInterpolation(), CNInterpolation(), FormVoronoiCell(), PlanarInterpolation(), and SetSearchSeed(). |
|
|
Definition at line 80 of file BFLInterpolation.h. Referenced by BFLInterpolation(), FormVoronoiCell(), NNInterpolation(), ResetForNextQuery(), and ~BFLInterpolation(). |
|
|
Definition at line 76 of file BFLInterpolation.h. Referenced by BFLInterpolation(), BilinearInterpolation(), CNInterpolation(), NNInterpolation(), PlanarInterpolation(), SetVoronoi(), and ~BFLInterpolation(). |
1.3.9.1