FreeFOAM The Cross-Platform CFD Toolkit
treeNode.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 Description
25 
26 \*---------------------------------------------------------------------------*/
27 
28 #include "treeNode.H"
29 #include "octree.H"
30 #include "treeLeaf.H"
31 #include "treeBoundBox.H"
32 #include <OpenFOAM/long.H>
33 #include <OpenFOAM/linePointRef.H>
34 
35 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
36 
37 template <class Type>
38 const Foam::label Foam::treeNode<Type>::leafOffset(100);
39 
40 
41 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
42 
43 template <class Type>
44 void Foam::treeNode<Type>::setAsNode(const label octant)
45 {
46  subNodeTypes_ |= (0x1 << octant);
47 }
48 
49 
50 template <class Type>
51 void Foam::treeNode<Type>::setAsLeaf(const label octant)
52 {
53  subNodeTypes_ &= ~(0x1 << octant);
54 }
55 
56 
57 // Set pointer to sub node
58 template <class Type>
60 (
61  const label octant,
62  treeElem<Type>* treeNodePtr
63 )
64 {
65  setAsNode(octant);
66  subNodes_[octant] = treeNodePtr;
67 }
68 
69 
70 // Set pointer to sub leaf
71 template <class Type>
73 (
74  const label octant,
75  treeElem<Type>* treeLeafPtr
76 )
77 {
78  setAsLeaf(octant);
79  subNodes_[octant] = treeLeafPtr;
80 }
81 
82 
83 template <class Type>
85 (
86  const label octant,
87  const label type
88 )
89 {
90  if ((type < 0) || (type > 3))
91  {
92  FatalErrorIn("treeNode<Type>::setVolType(const label, const label)")
93  << "Type " << type << " not within range 0..3" << endl;
94  }
95 
96  // Clear out two bits at position 2*octant
97  volType_ &= ~(0x3 << 2*octant);
98 
99  // Add the two bits of type
100  volType_ |= (type << 2*octant);
101 }
102 
103 
104 template <class Type>
105 void Foam::treeNode<Type>::space(Ostream& os, const label n)
106 {
107  for (label i=0; i<n; i++)
108  {
109  os<< ' ';
110  }
111 }
112 
113 
114 // look in single octant starting from <start>
115 template <class Type>
117 (
118  const int level,
119  const Type& shapes,
120  const label octant,
121  const vector& direction,
122  point& start,
123  const point& end
124 ) const
125 {
126  static const char* functionName =
127  "treeNode<Type>::findLeafLineOctant"
128  "(const int, const Type&, const label, const vector&,"
129  " point&, const point&)";
130 
131  if (debug & 2)
132  {
133  space(Pout, 2*level);
134  Pout<< "findLeafLineOctant : bb:" << this->bb()
135  << " start:" << start
136  << " end:" << end
137  << " mid:" << midpoint()
138  << " Searching octant:" << octant
139  << endl;
140  }
141 
142  if (subNodes()[octant])
143  {
144  if (isNode(octant))
145  {
146  // Node: recurse into subnodes
147  const treeNode<Type>* subNodePtr = getNodePtr(octant);
148 
149  if (subNodePtr->bb().contains(direction, start))
150  {
151  // Search on lower level
152  const treeLeaf<Type>* subLeafPtr = subNodePtr->findLeafLine
153  (
154  level + 1,
155  shapes,
156  start,
157  end
158  );
159 
160  if (debug & 2)
161  {
162  space(Pout, 2*level);
163  Pout<< "findLeafLineOctant : bb:" << this->bb()
164  << " returning from sub treeNode"
165  << " with start:" << start << " subLeaf:"
166  << long(subLeafPtr) << endl;
167  }
168 
169  return subLeafPtr;
170  }
171  else
172  {
173  FatalErrorIn(functionName)
174  << "Sub node " << subNodePtr->bb()
175  << " at octant " << octant
176  << " does not contain start " << start
177  << abort(FatalError);
178  }
179  }
180  else
181  {
182  // Leaf
183  const treeLeaf<Type>* subLeafPtr = getLeafPtr(octant);
184 
185  if (subLeafPtr->bb().contains(direction, start))
186  {
187  // Step to end of subleaf bb
188  point tmp;
189  if
190  (
191  !subLeafPtr->bb().intersects
192  (
193  end,
194  start,
195  tmp
196  )
197  )
198  {
199  FatalErrorIn(functionName)
200  << "Sub leaf contains start " << start
201  << " but line does not intersect its bb "
202  << subLeafPtr->bb()
203  << abort(FatalError);
204  }
205  start = tmp;
206 
207  if (debug & 2)
208  {
209  space(Pout, 2*level);
210  Pout<< "findLeafLineOctant : returning from intersecting"
211  << " treeLeaf " << subLeafPtr->bb()
212  << " with start:" << start << " subLeaf:"
213  << long(subLeafPtr) << endl;
214  }
215 
216  return subLeafPtr;
217  }
218  else
219  {
220  FatalErrorIn(functionName)
221  << "Sub leaf " << subLeafPtr->bb()
222  << " at octant " << octant
223  << " does not contain start " << start
224  << abort(FatalError);
225  }
226  }
227  }
228  else
229  {
230  // Empty subNode. Transfer across.
231  const treeBoundBox emptyBb = this->bb().subBbox(midpoint(), octant);
232 
233  if (emptyBb.contains(direction, start))
234  {
235  if (debug & 2)
236  {
237  space(Pout, 2*level);
238  Pout<< "findLeafLineOctant : Empty node. Octant:" << octant
239  << " start:" << start
240  << " bb:" << this->bb()
241  << " emptyBb:" << emptyBb << endl;
242  }
243 
244  // Update start by clipping to emptyBb
245  point tmp;
246  if
247  (
248  !emptyBb.intersects
249  (
250  end,
251  start,
252  tmp
253  )
254  )
255  {
256  FatalErrorIn(functionName)
257  << "Empty node contains start " << start
258  << " but line does not intersect its (calculated)"
259  << " bb " << emptyBb
260  << endl << "This might be due to truncation error"
261  << abort(FatalError);
262  }
263  start = tmp;
264 
265  if (debug & 2)
266  {
267  space(Pout, 2*level);
268  Pout<< "findLeafLineOctant : returning from intersecting with"
269  << " empty " << emptyBb
270  << " with start:" << start << " subLeaf:" << 0 << endl;
271  }
272 
273  return NULL;
274  }
275  else
276  {
277  FatalErrorIn(functionName)
278  << "Empty node " << emptyBb
279  << " at octant " << octant
280  << " does not contain start " << start
281  << abort(FatalError);
282  }
283  }
284 
285  FatalErrorIn(functionName)
286  << "Octant " << octant << " of cube " << this->bb()
287  << " does not contain start " << start
288  << abort(FatalError);
289 
290  return NULL;
291 }
292 
293 
294 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
295 
296 // Construct from components
297 template <class Type>
299 :
300  treeElem<Type>(bb),
301  treeNodeName(),
302  mid_(bb.midpoint()),
303  subNodeTypes_(0),
304  volType_(0)
305 {
306  for (label octantI=0; octantI<8; octantI++)
307  {
308  subNodes_[octantI] = NULL;
309  setVolType(octantI, octree<Type>::UNKNOWN);
310  }
311 }
312 
313 
314 // Construct from Istream
315 template <class Type>
317 {
318  is >> *this;
319 }
320 
321 
322 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
323 
324 template <class Type>
326 {
327  for (int octant=0; octant<8; octant++)
328  {
329  if (subNodes()[octant])
330  {
331  if (isNode(octant))
332  {
333  delete getNodePtr(octant);
334  }
335  else
336  {
337  delete getLeafPtr(octant);
338  }
339  }
340  }
341 }
342 
343 
344 // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
345 
346 // Distributes cells to subLeaves
347 template <class Type>
349 (
350  const label level,
351  octree<Type>& top,
352  const Type& shapes,
353  const labelList& indices
354 )
355 {
356  if (debug & 1)
357  {
358  space(Pout, level);
359  Pout<< "treeNode::distributing " << indices.size() << endl;
360  }
361 
362  // Create subLeaves if necessary
363  for (label octant=0; octant<8; octant++)
364  {
365  if (subNodes()[octant])
366  {
367  printNode(Pout, level);
369  (
370  "treeNode<Type>::distribute(const label, octree<Type>&, "
371  "const Type&, const labelList&)"
372  ) << "subNode already available at octant:" << octant
373  << abort(FatalError);
374  }
375  else
376  {
377  treeLeaf<Type>* subLeafPtr = new treeLeaf<Type>
378  (
379  this->bb().subBbox(midpoint(), octant),
380  indices.size()
381  );
382 
383  top.setLeaves(top.nLeaves() + 1);
384  setLeafPtr(octant, subLeafPtr);
385  }
386  }
387 
388 
389  // add cells to correct sub leaf
390  forAll(indices, i)
391  {
392  const label shapei = indices[i];
393 
394  for (label octant=0; octant<8; octant++)
395  {
396  treeLeaf<Type>* leafPtr = getLeafPtr(octant);
397 
398  if (shapes.overlaps(shapei, leafPtr->bb()))
399  {
400  if (debug == 1)
401  {
402  space(Pout, level);
403  Pout<< "inserting " << shapei;
404  shapes.write(Pout, shapei);
405  Pout<< " into " << leafPtr->bb() << endl;
406  }
407  leafPtr->insert(shapei);
408  top.setEntries(top.nEntries() + 1);
409  }
410  }
411  }
412 
413  // Trim size of subLeaves
414  for (label octant=0; octant<8; octant++)
415  {
416  treeLeaf<Type>* subLeafPtr = getLeafPtr(octant);
417 
418  if (subLeafPtr->size() == 0)
419  {
420  // Contains no data. Delete.
421  setLeafPtr(octant, NULL);
422  delete subLeafPtr;
423  top.setLeaves(top.nLeaves() - 1);
424  }
425  else
426  {
427  // Trim to actual size.
428  subLeafPtr->trim();
429  }
430  }
431 
432  if (debug & 1)
433  {
434  space(Pout, level);
435  Pout<< "end of treeNode::distribute" << endl;
436  }
437 }
438 
439 
440 // Descends to refineLevel and checks the subLeaves for redistribution
441 template <class Type>
443 (
444  const label level,
445  octree<Type>& top,
446  const Type& shapes,
447  const label refineLevel
448 )
449 {
450  if (debug & 1)
451  {
452  space(Pout, level);
453  Pout<< "treeNode::redistribute with level:" << level
454  << " refineLevel:" << refineLevel << endl;
455  }
456 
457  // Descend to correct level
458  if (level < refineLevel)
459  {
460  for (label octant=0; octant<8; octant++)
461  {
462  if (subNodes()[octant])
463  {
464  if (isNode(octant))
465  {
466  getNodePtr(octant)->redistribute
467  (
468  level+1,
469  top,
470  shapes,
471  refineLevel
472  );
473  }
474  }
475  }
476  }
477  else
478  {
479  // Reached correct (should also be deepest) level of treeNode
480  if (debug & 1)
481  {
482  space(Pout, level);
483  Pout<< "treeNode::redistribute : now at correct level" << endl;
484  }
485 
486  // handle redistribution of sub leaves
487  for (label octant=0; octant<8; octant++)
488  {
489  if (subNodes()[octant])
490  {
491  if (isNode(octant))
492  {
494  (
495  "treeNode<Type>::redistribute(const int, octree& top,"
496  "const int, const treeBoundBox&)"
497  ) << "found treeNode instead of treeLeaf" << endl
498  << abort(FatalError);
499  }
500  else
501  {
502  treeLeaf<Type>* leafPtr = getLeafPtr(octant);
503 
504  treeLeaf<Type>* newSubPtr = leafPtr->redistribute
505  (
506  level,
507  top,
508  shapes
509  );
510 
511  if (newSubPtr && (newSubPtr != leafPtr))
512  {
513  // redistribute has created nodePtr
514  // so delete no longer used subPtr and update info
515  if (debug & 1)
516  {
517  Pout<< "deleting "
518  << top.nEntries() - leafPtr->size()
519  << " entries" << endl;
520  }
521  top.setEntries(top.nEntries() - leafPtr->size());
522 
523  delete leafPtr;
524 
525  top.setLeaves(top.nLeaves() - 1);
526 
527  setNodePtr(octant, newSubPtr);
528  }
529  }
530  }
531  }
532  if (debug & 1)
533  {
534  space(Pout, level);
535  Pout<< "end of treeNode::redistribute for correct level" << endl;
536  }
537  }
538 
539  if (debug & 1)
540  {
541  space(Pout, level);
542  Pout<< "return from treeNode::redistribute with bb:" << this->bb()
543  << endl;
544  }
545 }
546 
547 
548 // Set type of node.
549 template <class Type>
551 (
552  const label level,
553  octree<Type>& top,
554  const Type& shapes
555 )
556 {
557  if (debug & 4)
558  {
559  space(Pout, level);
560  Pout<< "treeNode::setSubNodeType with level:" << level
561  << " bb:" << this->bb() << endl;
562  }
563 
564  label myType = -1;
565 
566  for (label octant=0; octant<8; octant++)
567  {
568  label subType = -1;
569 
570  if (subNodes()[octant])
571  {
572  if (isNode(octant))
573  {
574  subType = getNodePtr(octant)->setSubNodeType
575  (
576  level+1,
577  top,
578  shapes
579  );
580  }
581  else
582  {
583  subType = getLeafPtr(octant)->setSubNodeType
584  (
585  level+1,
586  top,
587  shapes
588  );
589  }
590  }
591  else
592  {
593  // No data in this one. Set type for octant acc. to its bounding
594  // box.
595  const treeBoundBox subBb = this->bb().subBbox(midpoint(), octant);
596 
597  subType = shapes.getSampleType(top, subBb.midpoint());
598  }
599 
600  if (debug & 4)
601  {
602  space(Pout, level);
603  Pout<< "treeNode::setSubNodeType : setting octant with bb:"
604  << this->bb().subBbox(midpoint(), octant)
605  << " to type:" << octree<Type>::volType(subType) << endl;
606  }
607  setVolType(octant, subType);
608 
609  // Combine sub node types into type for treeNode. Result is 'mixed' if
610  // types differ among subnodes.
611  if (myType == -1)
612  {
613  myType = subType;
614  }
615  else if (subType != myType)
616  {
617  myType = octree<Type>::MIXED;
618  }
619  }
620 
621  if (debug & 4)
622  {
623  space(Pout, level);
624  Pout<< "return from treeNode::setSubNodeType with type:"
625  << octree<Type>::volType(myType)
626  << " bb:" << this->bb() << endl;
627  }
628 
629  return myType;
630 }
631 
632 
633 // Get type of node.
634 template <class Type>
636 (
637  const label level,
638  const octree<Type>& top,
639  const Type& shapes,
640  const point& sample
641 ) const
642 {
643  if (debug & 4)
644  {
645  space(Pout, level);
646  Pout<< "treeNode::getSampleType with level:" << level
647  << " bb:" << this->bb() << " sample:" << sample << endl;
648  }
649 
650  // Determine octant of bb. If on edge just use whichever octant.
651  bool onEdge = false;
652 
653  label octant = this->bb().subOctant(midpoint(), sample, onEdge);
654 
655  label type = getVolType(octant);
656 
657  if (type == octree<Type>::MIXED)
658  {
659  // At this level multiple sub types. Recurse to resolve.
660  if (subNodes()[octant])
661  {
662  if (isNode(octant))
663  {
664  // Node: recurse into subnodes
665  type = getNodePtr(octant)->getSampleType
666  (
667  level + 1,
668  top,
669  shapes,
670  sample
671  );
672  }
673  else
674  {
675  // Leaf
676  type = getLeafPtr(octant)->getSampleType
677  (
678  level + 1,
679  top,
680  shapes,
681  sample
682  );
683  }
684  }
685  else
686  {
687  // Problem: empty subnode should have a type
689  (
690  "treeNode<Type>::getSampleType"
691  "(const label, octree<Type>&, const Type&, const point&)"
692  ) << "Empty node bb:" << this->bb().subBbox(midpoint(), octant)
693  << " has non-mixed type:"
694  << octree<Type>::volType(type)
695  << abort(FatalError);
696  }
697  }
698 
699  if (type == octree<Type>::MIXED)
700  {
702  (
703  "treeNode<Type>::getSampleType"
704  "(const label, octree<Type>&, const Type&, const point&)"
705  ) << "Type is MIXED when searching for " << sample
706  << " at level " << this->bb() << endl
707  << "This probably is because the octree has not been constructed"
708  << " with search facility." << exit(FatalError);
709  }
710 
711  if (debug & 4)
712  {
713  space(Pout, level);
714  Pout<< "return from treeNode::getSampleType with type:"
715  << octree<Type>::volType(type)
716  << " bb:" << this->bb()
717  << " sample:" << sample << endl;
718  }
719  return type;
720 }
721 
722 
723 template <class Type>
724 Foam::label Foam::treeNode<Type>::find
725 (
726  const Type& shapes,
727  const point& sample
728 ) const
729 {
730  // Find octant of sample. Don't care if on edge (since any item on edge
731  // will have been inserted in both subcubes)
732  bool onEdge = false;
733 
734  label octant = this->bb().subOctant(midpoint(), sample, onEdge);
735 
736  if (subNodes()[octant])
737  {
738  if (isNode(octant))
739  {
740  // Node: recurse into subnodes
741  return getNodePtr(octant)->find(shapes, sample);
742  }
743  else
744  {
745  // Leaf: let leaf::find handle this
746  return getLeafPtr(octant)->find(shapes, sample);
747  }
748  }
749  return -1;
750 }
751 
752 
753 template <class Type>
755 (
756  const Type& shapes,
757  const point& sample,
758  treeBoundBox& tightest
759 ) const
760 {
761  bool changed = false;
762  bool onEdge = false;
763  // Estimate for best place to start searching
764  label sampleOctant = this->bb().subOctant(midpoint(), sample, onEdge);
765 
766  // Go into all suboctants (one containing sample first) and update tightest.
767  // Order of visiting is if e.g. sampleOctant = 5:
768  // 5 1 2 3 4 0 6 7
769  for (label octantI=0; octantI<8; octantI++)
770  {
771  label octant;
772  if (octantI == 0)
773  {
774  // Use sampleOctant first
775  octant = sampleOctant;
776  }
777  else if (octantI == sampleOctant)
778  {
779  octant = 0;
780  }
781  else
782  {
783  octant = octantI;
784  }
785 
786  if (subNodes()[octant])
787  {
788  if (isNode(octant))
789  {
790  // Node: recurse into subnodes
791  const treeNode<Type>* subNodePtr = getNodePtr(octant);
792 
793  if (subNodePtr->bb().overlaps(tightest))
794  {
795  // there might be a better fit inside this subNode
796  changed |= subNodePtr->findTightest
797  (
798  shapes,
799  sample,
800  tightest
801  );
802  }
803  }
804  else
805  {
806  // Leaf: let leaf::find handle this
807  const treeLeaf<Type>* subLeafPtr = getLeafPtr(octant);
808 
809  if (subLeafPtr->bb().overlaps(tightest))
810  {
811  // there might be a better fit inside this subLeaf
812  changed |= subLeafPtr->findTightest
813  (
814  shapes,
815  sample,
816  tightest
817  );
818  }
819  }
820  }
821  }
822 
823  return changed;
824 }
825 
826 
827 template <class Type>
829 (
830  const Type& shapes,
831  const point& sample,
832  treeBoundBox& tightest,
833  label& tightestI,
834  scalar& tightestDist
835 ) const
836 {
837  if (debug & 8)
838  {
839  Pout<< "In findNearest with sample:" << sample << " cube:"
840  << this->bb() << " tightest:" << tightest << endl;
841  }
842 
843  bool changed = false;
844  bool onEdge = false;
845  // Estimate for best place to start searching
846  label sampleOctant = this->bb().subOctant(midpoint(), sample, onEdge);
847 
848  // Go into all suboctants (one containing sample first) and update tightest.
849  // Order of visiting is if e.g. sampleOctant = 5:
850  // 5 1 2 3 4 0 6 7
851  for (label octantI=0; octantI<8; octantI++)
852  {
853  label octant;
854  if (octantI == 0)
855  {
856  // Use sampleOctant first
857  octant = sampleOctant;
858  }
859  else if (octantI == sampleOctant)
860  {
861  octant = 0;
862  }
863  else
864  {
865  octant = octantI;
866  }
867 
868  if (subNodes()[octant])
869  {
870  if (isNode(octant))
871  {
872  // Node
873  const treeNode<Type>* subNodePtr = getNodePtr(octant);
874 
875  if (subNodePtr->bb().overlaps(tightest))
876  {
877  // there might be a better fit inside this subNode
878  changed |= subNodePtr->findNearest
879  (
880  shapes,
881  sample,
882  tightest,
883  tightestI,
884  tightestDist
885  );
886  }
887  }
888  else
889  {
890  // Leaf: let leaf::find handle this
891  const treeLeaf<Type>* subLeafPtr = getLeafPtr(octant);
892 
893  if (subLeafPtr->bb().overlaps(tightest))
894  {
895  // there might be a better fit inside this subNode
896  changed |= subLeafPtr->findNearest
897  (
898  shapes,
899  sample,
900  tightest,
901  tightestI,
902  tightestDist
903  );
904  }
905  }
906  }
907  }
908 
909  if (debug & 8)
910  {
911  Pout<< "Exiting findNearest for sample:" << sample << " cube:"
912  << this->bb() << " tightestI:" << tightestI << endl;
913  }
914 
915  return changed;
916 }
917 
918 
919 template <class Type>
921 (
922  const Type& shapes,
923  const linePointRef& ln,
924  treeBoundBox& tightest,
925  label& tightestI,
926  point& linePoint, // nearest point on line
927  point& shapePoint // nearest point on shape
928 ) const
929 {
930  bool changed = false;
931  bool onEdge = false;
932  // Estimate for best place to start searching
933  label sampleOctant = this->bb().subOctant(midpoint(), ln.centre(), onEdge);
934 
935  // Go into all suboctants (one containing sample first) and update tightest.
936  // Order of visiting is if e.g. sampleOctant = 5:
937  // 5 1 2 3 4 0 6 7
938  for (label octantI=0; octantI<8; octantI++)
939  {
940  label octant;
941  if (octantI == 0)
942  {
943  // Use sampleOctant first
944  octant = sampleOctant;
945  }
946  else if (octantI == sampleOctant)
947  {
948  octant = 0;
949  }
950  else
951  {
952  octant = octantI;
953  }
954 
955  if (subNodes()[octant])
956  {
957  if (isNode(octant))
958  {
959  // Node
960  const treeNode<Type>* subNodePtr = getNodePtr(octant);
961 
962  if (subNodePtr->bb().overlaps(tightest))
963  {
964  // there might be a better fit inside this subNode
965  changed |= subNodePtr->findNearest
966  (
967  shapes,
968  ln,
969  tightest,
970  tightestI,
971  linePoint,
972  shapePoint
973  );
974  }
975  }
976  else
977  {
978  // Leaf: let leaf::find handle this
979  const treeLeaf<Type>* subLeafPtr = getLeafPtr(octant);
980 
981  if (subLeafPtr->bb().overlaps(tightest))
982  {
983  // there might be a better fit inside this subNode
984  changed |= subLeafPtr->findNearest
985  (
986  shapes,
987  ln,
988  tightest,
989  tightestI,
990  linePoint,
991  shapePoint
992  );
993  }
994  }
995  }
996  }
997 
998  return changed;
999 }
1000 
1001 
1002 template <class Type>
1005  const Type& shapes,
1006  const boundBox& box,
1007  labelHashSet& elements
1008 ) const
1009 {
1010  bool changed = false;
1011  bool onEdge = false;
1012  // Estimate for best place to start searching
1013  label sampleOctant = this->bb().subOctant
1014  (
1015  midpoint(),
1016  box.midpoint(),
1017  onEdge
1018  );
1019 
1020  // Go into all suboctants (one containing sample first) and update tightest.
1021  // Order of visiting is if e.g. sampleOctant = 5:
1022  // 5 1 2 3 4 0 6 7
1023  for (label octantI=0; octantI<8; octantI++)
1024  {
1025  label octant;
1026  if (octantI == 0)
1027  {
1028  // Use sampleOctant first
1029  octant = sampleOctant;
1030  }
1031  else if (octantI == sampleOctant)
1032  {
1033  octant = 0;
1034  }
1035  else
1036  {
1037  octant = octantI;
1038  }
1039 
1040  if (subNodes()[octant])
1041  {
1042  if (isNode(octant))
1043  {
1044  // Node
1045  const treeNode<Type>* subNodePtr = getNodePtr(octant);
1046 
1047  if (subNodePtr->bb().overlaps(box))
1048  {
1049  // Visit sub node.
1050  changed |= subNodePtr->findBox(shapes, box, elements);
1051  }
1052  }
1053  else
1054  {
1055  // Leaf: let leaf::find handle this
1056  const treeLeaf<Type>* subLeafPtr = getLeafPtr(octant);
1057 
1058  if (subLeafPtr->bb().overlaps(box))
1059  {
1060  // Visit sub leaf.
1061  changed |= subLeafPtr->findBox(shapes, box, elements);
1062  }
1063  }
1064  }
1065  }
1066 
1067  return changed;
1068 }
1069 
1070 
1071 // look from <start> in current cube (given by this->bb()).
1072 template <class Type>
1075  const int level,
1076  const Type& shapes,
1077  point& start,
1078  const point& end
1079 ) const
1080 {
1081  if (debug & 2)
1082  {
1083  space(Pout, 2*level);
1084  Pout<< "findLeafLine : bb:" << this->bb() << " mid:" << midpoint()
1085  << " start:" << start << endl;
1086  }
1087 
1088  scalar typDim = this->bb().avgDim();
1089 
1090  const vector direction = end - start;
1091 
1092  // Loop on current level until start has been updated to be outside
1093  // of this->bb(). Note that max only four subcubes can be crossed so this is
1094  // check on whether there are any truncation error problems.
1095 
1096  label iter = 0;
1097 
1098  while (true)
1099  {
1100  if (!this->bb().contains(direction, start))
1101  {
1102  if (debug & 2)
1103  {
1104  space(Pout, 2*level);
1105  Pout<< "findLeafLine : Start not inside bb " << this->bb()
1106  << ". Returning with start:" << start << " subLeaf:"
1107  << 0 << endl;
1108  }
1109  return NULL;
1110  }
1111 
1112  // Check if start and <end> equal
1113  if ((mag(start - end)/typDim) < SMALL)
1114  {
1115  if (debug & 2)
1116  {
1117  space(Pout, 2*level);
1118  Pout<< "findLeafLine : start equals end"
1119  << ". Returning with start:" << start << " subLeaf:"
1120  << 0 << endl;
1121  }
1122  return NULL;
1123  }
1124 
1125  if (iter >= 4)
1126  {
1127  // Too many iterations. Is hanging. Handle outside of loop.
1128  break;
1129  }
1130 
1131  bool onEdge = false;
1132  label octant = this->bb().subOctant
1133  (
1134  midpoint(), direction, start, onEdge
1135  );
1136 
1137  // Try finding non-empty treeleaf in octant
1138  const treeLeaf<Type>* leafPtr = findLeafLineOctant
1139  (
1140  level,
1141  shapes,
1142  octant,
1143  direction,
1144  start,
1145  end
1146  );
1147 
1148  if (leafPtr)
1149  {
1150  // Found treeLeaf -> return
1151  if (debug & 2)
1152  {
1153  space(Pout, 2*level);
1154  Pout<< "findLeafLine : Found treeLeaf"
1155  << ". Returning with start:" << start << " subLeaf:"
1156  << long(leafPtr) << endl;
1157  }
1158 
1159  return leafPtr;
1160  }
1161 
1162  iter++;
1163  }
1164 
1165  // Check if is hanging. Max 4 octants can be crossed by a straight line
1166  FatalErrorIn
1167  (
1168  "treeNode<Type>::findLeafLine"
1169  "(const label, octree<Type>&, point&,"
1170  " const point&)"
1171  ) << "Did not leave bb " << this->bb()
1172  << " after " << iter
1173  << " iterations of updating starting point."
1174  << "start:" << start << " end:" << end
1175  << abort(FatalError);
1176 
1177  return NULL;
1178 }
1179 
1180 
1181 template <class Type>
1184  List<treeLeaf<Type>*>& leafArray,
1185  label& leafIndex
1186 ) const
1187 {
1188  // Go into all sub boxes
1189  for (label octant=0; octant<8; octant++)
1190  {
1191  if (subNodes()[octant])
1192  {
1193  if (isNode(octant))
1194  {
1195  // Node: recurse into subnodes
1196  const treeNode<Type>* subNodePtr = getNodePtr(octant);
1197  subNodePtr->findLeaves(leafArray, leafIndex);
1198  }
1199  else
1200  {
1201  // Leaf: store
1202  treeLeaf<Type>* subLeafPtr = getLeafPtr(octant);
1203  leafArray[leafIndex++] = subLeafPtr;
1204  }
1205  }
1206  }
1207 }
1208 
1209 
1210 template <class Type>
1213  List<const treeLeaf<Type>*>& leafArray,
1214  label& leafIndex
1215 ) const
1216 {
1217  // Go into all sub boxes
1218  for (label octant=0; octant<8; octant++)
1219  {
1220  if (subNodes()[octant])
1221  {
1222  if (isNode(octant))
1223  {
1224  // Node: recurse into subnodes
1225  const treeNode<Type>* subNodePtr = getNodePtr(octant);
1226  subNodePtr->findLeaves(leafArray, leafIndex);
1227  }
1228  else
1229  {
1230  // Leaf: store
1231  treeLeaf<Type>* subLeafPtr = getLeafPtr(octant);
1232  leafArray[leafIndex++] = subLeafPtr;
1233  }
1234  }
1235  }
1236 }
1237 
1238 
1239 template <class Type>
1242  Ostream& os,
1243  const label level
1244 ) const
1245 {
1246  space(os, 2*level);
1247 
1248  os << "node:" << this->bb() << endl;
1249 
1250  for (label octant=0; octant<8; octant++)
1251  {
1252  label type = getVolType(octant);
1253 
1254  string typeString = octree<Type>::volType(type);
1255 
1256  if (!subNodes_[octant])
1257  {
1258  space(os, level);
1259  os << octant << ":" << typeString << " : null" << endl;
1260  }
1261  else if (isNode(octant))
1262  {
1263  space(os, level);
1264  os << octant << ":" << typeString << " : node" << endl;
1265  getNodePtr(octant)->printNode(os, level+1);
1266  }
1267  else
1268  {
1269  space(os, level);
1270  os << octant << ":" << typeString << " : leaf" << endl;
1271 
1272  treeLeaf<Type>* leafPtr = getLeafPtr(octant);
1273  leafPtr->printLeaf(os, level+1);
1274  }
1275  }
1276 }
1277 
1278 
1279 template <class Type>
1282  Ostream& os,
1283  const label level,
1284  label& vertNo
1285 ) const
1286 {
1287  point midPoint(this->bb().midpoint());
1288 
1289  label midVertNo = vertNo;
1290  os << "v " << midPoint.x() << " " << midPoint.y() << " "
1291  << midPoint.z() << endl;
1292  vertNo++;
1293 
1294  for (label octant=0; octant<8; octant++)
1295  {
1296  if (subNodes_[octant])
1297  {
1298  if (isNode(octant))
1299  {
1300  treeNode<Type>* nodePtr = getNodePtr(octant);
1301 
1302  point subMidPoint(nodePtr->bb().midpoint());
1303  os << "v " << subMidPoint.x() << " " << subMidPoint.y() << " "
1304  << subMidPoint.z() << endl;
1305  os << "l " << midVertNo + 1<< " " << vertNo + 1 << endl;
1306  vertNo++;
1307 
1308  nodePtr->writeOBJ(os, level+1, vertNo);
1309  }
1310  else
1311  {
1312  treeLeaf<Type>* leafPtr = getLeafPtr(octant);
1313 
1314  point subMidPoint(leafPtr->bb().midpoint());
1315  os << "v " << subMidPoint.x() << " " << subMidPoint.y() << " "
1316  << subMidPoint.z() << endl;
1317  os << "l " << midVertNo + 1<< " " << vertNo + 1 << endl;
1318  vertNo++;
1319 
1320  //leafPtr->writeOBJ(os, level+1, vertNo);
1321  }
1322  }
1323  }
1324 }
1325 
1326 
1327 // * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
1328 
1329 template <class Type>
1331 {
1332  for (label octant = 0; octant < 8; octant++)
1333  {
1334  oc.subNodes_[octant] = NULL;
1335  }
1336 
1337  is >> oc.bb();
1338 
1339  label nPtrs;
1340 
1341  // Read number of entries folllowing
1342  is >> nPtrs;
1343 
1344  is.readBegin("treeNode");
1345  for (label octant = 0; octant < nPtrs; octant++)
1346  {
1347  label index;
1348  is >> index;
1349 
1350  if (index >= treeNode<Type>::leafOffset)
1351  {
1352  // Leaf recognized by out of range index
1353  treeLeaf<Type>* leafPtr = new treeLeaf<Type>(is);
1354  oc.setLeafPtr(index - treeNode<Type>::leafOffset, leafPtr);
1355  }
1356  else
1357  {
1358  oc.setNodePtr(index, new treeNode<Type>(is));
1359  }
1360  }
1361 
1362  // Read end of treeNode list
1363  is.readEnd("treeNode");
1364 
1365  // Check state of Istream
1366  is.check("Istream& operator>>(Istream&, treeNode&)");
1367 
1368  return is;
1369 }
1370 
1371 
1372 template <class Type>
1373 Foam::Ostream& Foam::operator<<(Ostream& os, const treeNode<Type>& tn)
1374 {
1375  // Count valid subnodes:
1376  // - treeNode
1377  // - treeLeafs with non-zero cell list.
1378  label nPtrs = 0;
1379  for (label octant = 0; octant < 8; octant++)
1380  {
1381  if (tn.subNodes_[octant])
1382  {
1383  if (tn.isNode(octant) || tn.getLeafPtr(octant)->indices().size())
1384  {
1385  nPtrs++;
1386  }
1387  }
1388  }
1389 
1390 
1391  // output subnodes as list of length nPtrs
1392  os << token::SPACE << tn.bb() << token::SPACE << nPtrs
1394 
1395  for (label octant = 0; octant < 8; octant++)
1396  {
1397  if (tn.subNodes_[octant])
1398  {
1399  if (tn.isNode(octant))
1400  {
1401  const treeNode<Type>* subNodePtr = tn.getNodePtr(octant);
1402 
1403  // Node: output index, value
1404  os << token::SPACE << octant << token::SPACE << *subNodePtr
1405  << token::NL;
1406  }
1407  else if (tn.getLeafPtr(octant)->indices().size())
1408  {
1409  // treeLeaf: mark by putting index invalid
1410  const treeLeaf<Type>* subLeafPtr = tn.getLeafPtr(octant);
1411 
1412  os << token::SPACE << octant + treeNode<Type>::leafOffset
1413  << token::SPACE << *subLeafPtr
1414  << token::NL;
1415  }
1416  }
1417  }
1418 
1419  os << token::SPACE << token::END_LIST;
1420 
1421  return os;
1422 }
1423 
1424 
1425 // ************************ vim: set sw=4 sts=4 et: ************************ //