Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/advanced/dna/moleculardna/src/OctreeNode.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /examples/advanced/dna/moleculardna/src/OctreeNode.cc (Version 11.3.0) and /examples/advanced/dna/moleculardna/src/OctreeNode.cc (Version 7.1)


  1 //                                                  1 
  2 // *******************************************    
  3 // * License and Disclaimer                       
  4 // *                                              
  5 // * The  Geant4 software  is  copyright of th    
  6 // * the Geant4 Collaboration.  It is provided    
  7 // * conditions of the Geant4 Software License    
  8 // * LICENSE and available at  http://cern.ch/    
  9 // * include a list of copyright holders.         
 10 // *                                              
 11 // * Neither the authors of this software syst    
 12 // * institutes,nor the agencies providing fin    
 13 // * work  make  any representation or  warran    
 14 // * regarding  this  software system or assum    
 15 // * use.  Please see the license in the file     
 16 // * for the full disclaimer and the limitatio    
 17 // *                                              
 18 // * This  code  implementation is the result     
 19 // * technical work of the GEANT4 collaboratio    
 20 // * By using,  copying,  modifying or  distri    
 21 // * any work based  on the software)  you  ag    
 22 // * use  in  resulting  scientific  publicati    
 23 // * acceptance of all terms of the Geant4 Sof    
 24 // *******************************************    
 25 //                                                
 26 //                                                
 27 #include "OctreeNode.hh"                          
 28                                                   
 29 #include "G4VPhysicalVolume.hh"                   
 30                                                   
 31 #include <cassert>                                
 32 #include <utility>                                
 33                                                   
 34 //....oooOO0OOooo........oooOO0OOooo........oo    
 35                                                   
 36 OctreeNode::OctreeNode(const G4ThreeVector& po    
 37                        G4int maxContents, Octr    
 38   : fPosition(position),                          
 39     fHalfLengths(halfLengths),                    
 40     fMaxContents(maxContents),                    
 41     fParent(parent),                              
 42     fChildren(std::array<OctreeNode*, 8>())       
 43 {                                                 
 44   fHalfLengthsMag = fHalfLengths.mag();           
 45 }                                                 
 46                                                   
 47 OctreeNode::~OctreeNode()                         
 48 {                                                 
 49   // Delete children                              
 50   for (auto& it : fChildren) {                    
 51     delete it;                                    
 52   }                                               
 53 }                                                 
 54                                                   
 55 //....oooOO0OOooo........oooOO0OOooo........oo    
 56                                                   
 57 const OctreeNode* OctreeNode::GetChildFromPosi    
 58 {                                                 
 59   G4ThreeVector localPosition = pos - fPositio    
 60                                                   
 61   // Get the right child                          
 62   // The scheme is defined by quadrant as foll    
 63   // Zero is treated as positive                  
 64   // X Y Z Index | X Y Z Index                    
 65   // + + +   0   | - + +   4                      
 66   // + + -   1   | - + -   5                      
 67   // + - +   2   | - - +   6                      
 68   // + - -   3   | - - -   7                      
 69   int bit2 = (localPosition.getX() < 0) ? 4 :     
 70   int bit1 = (localPosition.getY() < 0) ? 2 :     
 71   int bit0 = (localPosition.getZ() < 0) ? 1 :     
 72   return fChildren[bit0 + bit1 + bit2];           
 73 }                                                 
 74                                                   
 75 //....oooOO0OOooo........oooOO0OOooo........oo    
 76                                                   
 77 // Non-const version of above function. Ugly b    
 78 OctreeNode* OctreeNode::GetChildFromPosition(c    
 79 {                                                 
 80   G4ThreeVector localPosition = pos - fPositio    
 81                                                   
 82   // Get the right child                          
 83   // The scheme is defined by quadrant as foll    
 84   // Zero is treated as positive                  
 85   // X Y Z Index | X Y Z Index                    
 86   // + + +   0   | - + +   4                      
 87   // + + -   1   | - + -   5                      
 88   // + - +   2   | - - +   6                      
 89   // + - -   3   | - - -   7                      
 90   int bit2 = (localPosition.getX() < 0) ? 4 :     
 91   int bit1 = (localPosition.getY() < 0) ? 2 :     
 92   int bit0 = (localPosition.getZ() < 0) ? 1 :     
 93   int idx = bit0 + bit1 + bit2;                   
 94                                                   
 95   assert(idx < (int)fChildren.size());            
 96   return fChildren[idx];                          
 97 }                                                 
 98                                                   
 99 //....oooOO0OOooo........oooOO0OOooo........oo    
100                                                   
101 void OctreeNode::AddPhysicalVolume(G4VPhysical    
102 {                                                 
103   // Consider adding test for the validity of     
104   if (this->HasChildren()) {                      
105     OctreeNode* child = this->GetChildFromPosi    
106     child->AddPhysicalVolume(pv);                 
107   }                                               
108   else {                                          
109     // if there are fMaxContents elements in t    
110     // the bin before adding a new element.       
111     if ((G4int)fContents.size() == fMaxContent    
112       this->Split();                              
113       this->AddPhysicalVolume(pv);                
114     }                                             
115     else {                                        
116       fContents.push_back(pv);                    
117     }                                             
118   }                                               
119 }                                                 
120                                                   
121 //....oooOO0OOooo........oooOO0OOooo........oo    
122                                                   
123 void OctreeNode::Split()                          
124 {                                                 
125   G4ThreeVector hl = 0.5 * fHalfLengths;          
126   G4double xp = fPosition.getX();                 
127   G4double yp = fPosition.getY();                 
128   G4double zp = fPosition.getZ();                 
129   G4double xl = hl.getX();                        
130   G4double yl = hl.getY();                        
131   G4double zl = hl.getZ();                        
132   G4ThreeVector newpos;                           
133                                                   
134   // Construct children                           
135   newpos = G4ThreeVector(xp + xl, yp + yl, zp     
136   fChildren[0] = new OctreeNode(newpos, hl, fM    
137   newpos = G4ThreeVector(xp + xl, yp + yl, zp     
138   fChildren[1] = new OctreeNode(newpos, hl, fM    
139   newpos = G4ThreeVector(xp + xl, yp - yl, zp     
140   fChildren[2] = new OctreeNode(newpos, hl, fM    
141   newpos = G4ThreeVector(xp + xl, yp - yl, zp     
142   fChildren[3] = new OctreeNode(newpos, hl, fM    
143   newpos = G4ThreeVector(xp - xl, yp + yl, zp     
144   fChildren[4] = new OctreeNode(newpos, hl, fM    
145   newpos = G4ThreeVector(xp - xl, yp + yl, zp     
146   fChildren[5] = new OctreeNode(newpos, hl, fM    
147   newpos = G4ThreeVector(xp - xl, yp - yl, zp     
148   fChildren[6] = new OctreeNode(newpos, hl, fM    
149   newpos = G4ThreeVector(xp - xl, yp - yl, zp     
150   fChildren[7] = new OctreeNode(newpos, hl, fM    
151                                                   
152   // Distribute contents to children              
153                                                   
154   for (int i = fContents.size() - 1; i >= 0; -    
155     G4VPhysicalVolume* pv = fContents[i];         
156     OctreeNode* child = this->GetChildFromPosi    
157     assert(child != this);                        
158     child->AddPhysicalVolume(pv);                 
159   }                                               
160   fContents.clear();                              
161 }                                                 
162                                                   
163 //....oooOO0OOooo........oooOO0OOooo........oo    
164                                                   
165 std::vector<G4VPhysicalVolume*> OctreeNode::Ge    
166 {                                                 
167   if (this->HasChildren())  // if has sub-node    
168   {                                               
169     std::vector<G4VPhysicalVolume*> vec;          
170     std::vector<G4VPhysicalVolume*> childCont;    
171     for (auto it = fChildren.begin(); it != fC    
172       childCont = (*it)->GetContents();           
173       for (auto jt = childCont.begin(); jt !=     
174         vec.push_back(*jt);                       
175       }                                           
176     }                                             
177     return vec;                                   
178   }                                               
179   else  // if no sub-nodes                        
180   {                                               
181     return fContents;                             
182   }                                               
183 }                                                 
184                                                   
185 //....oooOO0OOooo........oooOO0OOooo........oo    
186                                                   
187 // Search for Octree nodes in a volume            
188 const std::vector<G4VPhysicalVolume*> OctreeNo    
189                                                   
190 {                                                 
191   // Need to search based on absolute position    
192   // relative  position, and maintain a list o    
193   // final nodes within the radius                
194                                                   
195   std::vector<const OctreeNode*> nodes;           
196   std::vector<const OctreeNode*> candidates;      
197   nodes.reserve(512);                             
198   candidates.reserve(512);                        
199                                                   
200   if (this->HasChildren()) {                      
201     for (auto it = this->GetChildren().begin()    
202       candidates.push_back(*it);                  
203     }                                             
204   }                                               
205   else {                                          
206     candidates.push_back(this);                   
207   }                                               
208                                                   
209   const OctreeNode* aNode;                        
210   G4double d;                                     
211   while (!candidates.empty()) {                   
212     aNode = candidates.back();                    
213     candidates.pop_back();                        
214     d = (aNode->GetPosition() - pos).mag() - a    
215     // if node within circle                      
216     if (d < _rad) {                               
217       if (aNode->HasChildren()) {                 
218         for (auto it = aNode->GetChildren().be    
219           candidates.push_back(*it);              
220         }                                         
221       }                                           
222       else {                                      
223         nodes.push_back(aNode);                   
224       }                                           
225     }                                             
226   }                                               
227                                                   
228   // G4cout << "Found " << nodes.size() << " n    
229                                                   
230   // Get the physical volumes                     
231   G4double r2 = _rad * _rad;                      
232   std::vector<G4VPhysicalVolume*> pvols;          
233   for (auto it = nodes.begin(); it != nodes.en    
234     std::vector<G4VPhysicalVolume*> conts = (*    
235     for (auto jt = conts.begin(); jt != conts.    
236       if ((pos - (*jt)->GetTranslation()).mag2    
237         pvols.push_back((*jt));                   
238       }                                           
239     }                                             
240   }                                               
241   // G4cout << "Found " << pvols.size() << " c    
242   return pvols;                                   
243 }                                                 
244                                                   
245 //....oooOO0OOooo........oooOO0OOooo........oo    
246                                                   
247 // Search for Octree nodes in a volume            
248 void OctreeNode::SearchOctree(const G4ThreeVec    
249                               G4double _rad) c    
250 {                                                 
251   // Need to search based on absolute position    
252   // relative  position, and maintain a list o    
253   // final nodes within the radius                
254                                                   
255   std::vector<const OctreeNode*> nodes;           
256   std::vector<const OctreeNode*> candidates;      
257   nodes.reserve(512);                             
258   candidates.reserve(512);                        
259                                                   
260   if (this->HasChildren()) {                      
261     for (auto it = this->GetChildren().begin()    
262       candidates.push_back(*it);                  
263     }                                             
264   }                                               
265   else {                                          
266     candidates.push_back(this);                   
267   }                                               
268                                                   
269   const OctreeNode* aNode;                        
270   G4double d;                                     
271   while (!candidates.empty()) {                   
272     aNode = candidates.back();                    
273     candidates.pop_back();                        
274     d = (aNode->GetPosition() - pos).mag() - a    
275     // if node within circle                      
276     if (d < _rad) {                               
277       if (aNode->HasChildren()) {                 
278         for (auto it = aNode->GetChildren().be    
279           candidates.push_back(*it);              
280         }                                         
281       }                                           
282       else {                                      
283         nodes.push_back(aNode);                   
284       }                                           
285     }                                             
286   }                                               
287                                                   
288   // Get the physical volumes                     
289   G4double r2 = _rad * _rad;                      
290   for (auto it = nodes.begin(); it != nodes.en    
291     std::vector<G4VPhysicalVolume*> conts = (*    
292     for (auto jt = conts.begin(); jt != conts.    
293       if ((pos - (*jt)->GetTranslation()).mag2    
294         out.push_back((*jt));                     
295       }                                           
296     }                                             
297   }                                               
298 }                                                 
299                                                   
300 //....oooOO0OOooo........oooOO0OOooo........oo    
301                                                   
302 const std::vector<G4VPhysicalVolume*> OctreeNo    
303 {                                                 
304   const OctreeNode* child = GetChildFromPositi    
305   return child->GetContents();                    
306 }                                                 
307                                                   
308 //....oooOO0OOooo........oooOO0OOooo........oo    
309                                                   
310 G4int OctreeNode::GetNumberOfTerminalNodes()      
311 {                                                 
312   if (!this->HasChildren()) {                     
313     return 1;                                     
314   }                                               
315   // sum the number of nodes of each child        
316   G4int n = 0;                                    
317   for (auto it = this->GetChildren().begin();     
318     n = n + (*it)->GetNumberOfTerminalNodes();    
319   }                                               
320                                                   
321   return n;                                       
322 }                                                 
323 //....oooOO0OOooo........oooOO0OOooo........oo    
324