FreeFOAM The Cross-Platform CFD Toolkit
meshSearch.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 \*---------------------------------------------------------------------------*/
25 
26 #include "meshSearch.H"
27 #include <OpenFOAM/polyMesh.H>
29 #include <OpenFOAM/DynamicList.H>
31 #include <meshTools/treeDataCell.H>
32 #include <meshTools/treeDataFace.H>
34 
35 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
36 
38 
39 Foam::scalar Foam::meshSearch::tol_ = 1E-3;
40 
41 
42 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
43 
44 bool Foam::meshSearch::findNearer
45 (
46  const point& sample,
47  const pointField& points,
48  label& nearestI,
49  scalar& nearestDistSqr
50 )
51 {
52  bool nearer = false;
53 
54  forAll(points, pointI)
55  {
56  scalar distSqr = magSqr(points[pointI] - sample);
57 
58  if (distSqr < nearestDistSqr)
59  {
60  nearestDistSqr = distSqr;
61  nearestI = pointI;
62  nearer = true;
63  }
64  }
65 
66  return nearer;
67 }
68 
69 
70 bool Foam::meshSearch::findNearer
71 (
72  const point& sample,
73  const pointField& points,
74  const labelList& indices,
75  label& nearestI,
76  scalar& nearestDistSqr
77 )
78 {
79  bool nearer = false;
80 
81  forAll(indices, i)
82  {
83  label pointI = indices[i];
84 
85  scalar distSqr = magSqr(points[pointI] - sample);
86 
87  if (distSqr < nearestDistSqr)
88  {
89  nearestDistSqr = distSqr;
90  nearestI = pointI;
91  nearer = true;
92  }
93  }
94 
95  return nearer;
96 }
97 
98 
99 // tree based searching
100 Foam::label Foam::meshSearch::findNearestCellTree(const point& location) const
101 {
102  const indexedOctree<treeDataPoint>& tree = cellCentreTree();
103 
104  scalar span = tree.bb().mag();
105 
106  pointIndexHit info = tree.findNearest(location, Foam::sqr(span));
107 
108  if (!info.hit())
109  {
110  info = tree.findNearest(location, Foam::sqr(GREAT));
111  }
112  return info.index();
113 }
114 
115 
116 // linear searching
117 Foam::label Foam::meshSearch::findNearestCellLinear(const point& location) const
118 {
119  const vectorField& centres = mesh_.cellCentres();
120 
121  label nearestIndex = 0;
122  scalar minProximity = magSqr(centres[nearestIndex] - location);
123 
124  findNearer
125  (
126  location,
127  centres,
128  nearestIndex,
129  minProximity
130  );
131 
132  return nearestIndex;
133 }
134 
135 
136 // walking from seed
137 Foam::label Foam::meshSearch::findNearestCellWalk
138 (
139  const point& location,
140  const label seedCellI
141 ) const
142 {
143  if (seedCellI < 0)
144  {
146  (
147  "meshSearch::findNearestCellWalk(const point&, const label)"
148  ) << "illegal seedCell:" << seedCellI << exit(FatalError);
149  }
150 
151  // Walk in direction of face that decreases distance
152 
153  label curCellI = seedCellI;
154  scalar distanceSqr = magSqr(mesh_.cellCentres()[curCellI] - location);
155 
156  bool closer;
157 
158  do
159  {
160  // Try neighbours of curCellI
161  closer = findNearer
162  (
163  location,
164  mesh_.cellCentres(),
165  mesh_.cellCells()[curCellI],
166  curCellI,
167  distanceSqr
168  );
169  } while (closer);
170 
171  return curCellI;
172 }
173 
174 
175 // tree based searching
176 Foam::label Foam::meshSearch::findNearestFaceTree(const point& location) const
177 {
178  // Search nearest cell centre.
179  const indexedOctree<treeDataPoint>& tree = cellCentreTree();
180 
181  scalar span = tree.bb().mag();
182 
183  // Search with decent span
184  pointIndexHit info = tree.findNearest(location, Foam::sqr(span));
185 
186  if (!info.hit())
187  {
188  // Search with desparate span
189  info = tree.findNearest(location, Foam::sqr(GREAT));
190  }
191 
192 
193  // Now check any of the faces of the nearest cell
194  const vectorField& centres = mesh_.faceCentres();
195  const cell& ownFaces = mesh_.cells()[info.index()];
196 
197  label nearestFaceI = ownFaces[0];
198  scalar minProximity = magSqr(centres[nearestFaceI] - location);
199 
200  findNearer
201  (
202  location,
203  centres,
204  ownFaces,
205  nearestFaceI,
206  minProximity
207  );
208 
209  return nearestFaceI;
210 }
211 
212 
213 // linear searching
214 Foam::label Foam::meshSearch::findNearestFaceLinear(const point& location) const
215 {
216  const vectorField& centres = mesh_.faceCentres();
217 
218  label nearestFaceI = 0;
219  scalar minProximity = magSqr(centres[nearestFaceI] - location);
220 
221  findNearer
222  (
223  location,
224  centres,
225  nearestFaceI,
226  minProximity
227  );
228 
229  return nearestFaceI;
230 }
231 
232 
233 // walking from seed
234 Foam::label Foam::meshSearch::findNearestFaceWalk
235 (
236  const point& location,
237  const label seedFaceI
238 ) const
239 {
240  if (seedFaceI < 0)
241  {
243  (
244  "meshSearch::findNearestFaceWalk(const point&, const label)"
245  ) << "illegal seedFace:" << seedFaceI << exit(FatalError);
246  }
247 
248  const vectorField& centres = mesh_.faceCentres();
249 
250 
251  // Walk in direction of face that decreases distance
252 
253  label curFaceI = seedFaceI;
254  scalar distanceSqr = magSqr(centres[curFaceI] - location);
255 
256  while (true)
257  {
258  label betterFaceI = curFaceI;
259 
260  findNearer
261  (
262  location,
263  centres,
264  mesh_.cells()[mesh_.faceOwner()[curFaceI]],
265  betterFaceI,
266  distanceSqr
267  );
268 
269  if (mesh_.isInternalFace(curFaceI))
270  {
271  findNearer
272  (
273  location,
274  centres,
275  mesh_.cells()[mesh_.faceNeighbour()[curFaceI]],
276  betterFaceI,
277  distanceSqr
278  );
279  }
280 
281  if (betterFaceI == curFaceI)
282  {
283  break;
284  }
285 
286  curFaceI = betterFaceI;
287  }
288 
289  return curFaceI;
290 }
291 
292 
293 Foam::label Foam::meshSearch::findCellLinear(const point& location) const
294 {
295  bool cellFound = false;
296  label n = 0;
297 
298  label cellI = -1;
299 
300  while ((!cellFound) && (n < mesh_.nCells()))
301  {
302  if (pointInCell(location, n))
303  {
304  cellFound = true;
305  cellI = n;
306  }
307  else
308  {
309  n++;
310  }
311  }
312  if (cellFound)
313  {
314  return cellI;
315  }
316  else
317  {
318  return -1;
319  }
320 }
321 
322 
323 Foam::label Foam::meshSearch::findNearestBoundaryFaceWalk
324 (
325  const point& location,
326  const label seedFaceI
327 ) const
328 {
329  if (seedFaceI < 0)
330  {
332  (
333  "meshSearch::findNearestBoundaryFaceWalk"
334  "(const point&, const label)"
335  ) << "illegal seedFace:" << seedFaceI << exit(FatalError);
336  }
337 
338  // Start off from seedFaceI
339 
340  label curFaceI = seedFaceI;
341 
342  const face& f = mesh_.faces()[curFaceI];
343 
344  scalar minDist = f.nearestPoint
345  (
346  location,
347  mesh_.points()
348  ).distance();
349 
350  bool closer;
351 
352  do
353  {
354  closer = false;
355 
356  // Search through all neighbouring boundary faces by going
357  // across edges
358 
359  label lastFaceI = curFaceI;
360 
361  const labelList& myEdges = mesh_.faceEdges()[curFaceI];
362 
363  forAll (myEdges, myEdgeI)
364  {
365  const labelList& neighbours = mesh_.edgeFaces()[myEdges[myEdgeI]];
366 
367  // Check any face which uses edge, is boundary face and
368  // is not curFaceI itself.
369 
370  forAll(neighbours, nI)
371  {
372  label faceI = neighbours[nI];
373 
374  if
375  (
376  (faceI >= mesh_.nInternalFaces())
377  && (faceI != lastFaceI)
378  )
379  {
380  const face& f = mesh_.faces()[faceI];
381 
382  pointHit curHit = f.nearestPoint
383  (
384  location,
385  mesh_.points()
386  );
387 
388  // If the face is closer, reset current face and distance
389  if (curHit.distance() < minDist)
390  {
391  minDist = curHit.distance();
392  curFaceI = faceI;
393  closer = true; // a closer neighbour has been found
394  }
395  }
396  }
397  }
398  } while (closer);
399 
400  return curFaceI;
401 }
402 
403 
404 Foam::vector Foam::meshSearch::offset
405 (
406  const point& bPoint,
407  const label bFaceI,
408  const vector& dir
409 ) const
410 {
411  // Get the neighbouring cell
412  label ownerCellI = mesh_.faceOwner()[bFaceI];
413 
414  const point& c = mesh_.cellCentres()[ownerCellI];
415 
416  // Typical dimension: distance from point on face to cell centre
417  scalar typDim = mag(c - bPoint);
418 
419  return tol_*typDim*dir;
420 }
421 
422 
423 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
424 
425 // Construct from components
426 Foam::meshSearch::meshSearch(const polyMesh& mesh, const bool faceDecomp)
427 :
428  mesh_(mesh),
429  faceDecomp_(faceDecomp),
430  cloud_(mesh_, IDLList<passiveParticle>()),
431  boundaryTreePtr_(NULL),
432  cellTreePtr_(NULL),
433  cellCentreTreePtr_(NULL)
434 {}
435 
436 
437 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
438 
440 {
441  clearOut();
442 }
443 
444 
445 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
446 
448  const
449 {
450  if (!boundaryTreePtr_)
451  {
452  //
453  // Construct tree
454  //
455 
456  // all boundary faces (not just walls)
457  labelList bndFaces(mesh_.nFaces()-mesh_.nInternalFaces());
458  forAll(bndFaces, i)
459  {
460  bndFaces[i] = mesh_.nInternalFaces() + i;
461  }
462 
463  treeBoundBox overallBb(mesh_.points());
464  Random rndGen(123456);
465  overallBb = overallBb.extend(rndGen, 1E-4);
466  overallBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
467  overallBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
468 
469  boundaryTreePtr_ = new indexedOctree<treeDataFace>
470  (
471  treeDataFace // all information needed to search faces
472  (
473  false, // do not cache bb
474  mesh_,
475  bndFaces // boundary faces only
476  ),
477  overallBb, // overall search domain
478  8, // maxLevel
479  10, // leafsize
480  3.0 // duplicity
481  );
482  }
483 
484  return *boundaryTreePtr_;
485 }
486 
487 
489  const
490 {
491  if (!cellTreePtr_)
492  {
493  //
494  // Construct tree
495  //
496 
497  treeBoundBox overallBb(mesh_.points());
498  Random rndGen(123456);
499  overallBb = overallBb.extend(rndGen, 1E-4);
500  overallBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
501  overallBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
502 
503  cellTreePtr_ = new indexedOctree<treeDataCell>
504  (
506  (
507  false, // not cache bb
508  mesh_
509  ),
510  overallBb, // overall search domain
511  8, // maxLevel
512  10, // leafsize
513  3.0 // duplicity
514  );
515  }
516 
517  return *cellTreePtr_;
518 
519 }
520 
521 
524 {
525  if (!cellCentreTreePtr_)
526  {
527  //
528  // Construct tree
529  //
530 
531  treeBoundBox overallBb(mesh_.cellCentres());
532  Random rndGen(123456);
533  overallBb = overallBb.extend(rndGen, 1E-4);
534  overallBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
535  overallBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
536 
537  cellCentreTreePtr_ = new indexedOctree<treeDataPoint>
538  (
539  treeDataPoint(mesh_.cellCentres()),
540  overallBb, // overall search domain
541  10, // max levels
542  10.0, // maximum ratio of cubes v.s. cells
543  100.0 // max. duplicity; n/a since no bounding boxes.
544  );
545  }
546 
547  return *cellCentreTreePtr_;
548 }
549 
550 
551 // Is the point in the cell
552 // Works by checking if there is a face inbetween the point and the cell
553 // centre.
554 // Check for internal uses proper face decomposition or just average normal.
555 bool Foam::meshSearch::pointInCell(const point& p, label cellI) const
556 {
557  if (faceDecomp_)
558  {
559  const point& ctr = mesh_.cellCentres()[cellI];
560 
561  vector dir(p - ctr);
562  scalar magDir = mag(dir);
563 
564  // Check if any faces are hit by ray from cell centre to p.
565  // If none -> p is in cell.
566  const labelList& cFaces = mesh_.cells()[cellI];
567 
568  // Make sure half_ray does not pick up any faces on the wrong
569  // side of the ray.
570  scalar oldTol = intersection::setPlanarTol(0.0);
571 
572  forAll(cFaces, i)
573  {
574  label faceI = cFaces[i];
575 
576  pointHit inter = mesh_.faces()[faceI].ray
577  (
578  ctr,
579  dir,
580  mesh_.points(),
583  );
584 
585  if (inter.hit())
586  {
587  scalar dist = inter.distance();
588 
589  if (dist < magDir)
590  {
591  // Valid hit. Hit face so point is not in cell.
593 
594  return false;
595  }
596  }
597  }
598 
600 
601  // No face inbetween point and cell centre so point is inside.
602  return true;
603  }
604  else
605  {
606  const labelList& f = mesh_.cells()[cellI];
607  const labelList& owner = mesh_.faceOwner();
608  const vectorField& cf = mesh_.faceCentres();
609  const vectorField& Sf = mesh_.faceAreas();
610 
611  forAll(f, facei)
612  {
613  label nFace = f[facei];
614  vector proj = p - cf[nFace];
615  vector normal = Sf[nFace];
616  if (owner[nFace] == cellI)
617  {
618  if ((normal & proj) > 0)
619  {
620  return false;
621  }
622  }
623  else
624  {
625  if ((normal & proj) < 0)
626  {
627  return false;
628  }
629  }
630  }
631 
632  return true;
633  }
634 }
635 
636 
638 (
639  const point& location,
640  const label seedCellI,
641  const bool useTreeSearch
642 ) const
643 {
644  if (seedCellI == -1)
645  {
646  if (useTreeSearch)
647  {
648  return findNearestCellTree(location);
649  }
650  else
651  {
652  return findNearestCellLinear(location);
653  }
654  }
655  else
656  {
657  return findNearestCellWalk(location, seedCellI);
658  }
659 }
660 
661 
663 (
664  const point& location,
665  const label seedFaceI,
666  const bool useTreeSearch
667 ) const
668 {
669  if (seedFaceI == -1)
670  {
671  if (useTreeSearch)
672  {
673  return findNearestFaceTree(location);
674  }
675  else
676  {
677  return findNearestFaceLinear(location);
678  }
679  }
680  else
681  {
682  return findNearestFaceWalk(location, seedFaceI);
683  }
684 }
685 
686 
687 Foam::label Foam::meshSearch::findCell
688 (
689  const point& location,
690  const label seedCellI,
691  const bool useTreeSearch
692 ) const
693 {
694  // Find the nearest cell centre to this location
695  label nearCellI = findNearestCell(location, seedCellI, useTreeSearch);
696 
697  if (debug)
698  {
699  Pout<< "findCell : nearCellI:" << nearCellI
700  << " ctr:" << mesh_.cellCentres()[nearCellI]
701  << endl;
702  }
703 
704  if (pointInCell(location, nearCellI))
705  {
706  return nearCellI;
707  }
708  else
709  {
710  if (useTreeSearch)
711  {
712  // Start searching from cell centre of nearCell
713  point curPoint = mesh_.cellCentres()[nearCellI];
714 
715  vector edgeVec = location - curPoint;
716  edgeVec /= mag(edgeVec);
717 
718  do
719  {
720  // Walk neighbours (using tracking) to get closer
721  passiveParticle tracker(cloud_, curPoint, nearCellI);
722 
723  if (debug)
724  {
725  Pout<< "findCell : tracked from curPoint:" << curPoint
726  << " nearCellI:" << nearCellI;
727  }
728 
729  tracker.track(location);
730 
731  if (debug)
732  {
733  Pout<< " to " << tracker.position()
734  << " need:" << location
735  << " onB:" << tracker.onBoundary()
736  << " cell:" << tracker.cell()
737  << " face:" << tracker.face() << endl;
738  }
739 
740  if (!tracker.onBoundary())
741  {
742  // stopped not on boundary -> reached location
743  return tracker.cell();
744  }
745 
746  // stopped on boundary face. Compare positions
747  scalar typDim = sqrt(mag(mesh_.faceAreas()[tracker.face()]));
748 
749  if ((mag(tracker.position() - location)/typDim) < SMALL)
750  {
751  return tracker.cell();
752  }
753 
754  // tracking stopped at boundary. Find next boundary
755  // intersection from current point (shifted out a little bit)
756  // in the direction of the destination
757 
758  curPoint =
759  tracker.position()
760  + offset(tracker.position(), tracker.face(), edgeVec);
761 
762  if (debug)
763  {
764  Pout<< "Searching for next boundary from curPoint:"
765  << curPoint
766  << " to location:" << location << endl;
767  }
768  pointIndexHit curHit = intersection(curPoint, location);
769  if (debug)
770  {
771  Pout<< "Returned from line search with ";
772  curHit.write(Pout);
773  Pout<< endl;
774  }
775 
776  if (!curHit.hit())
777  {
778  return -1;
779  }
780  else
781  {
782  nearCellI = mesh_.faceOwner()[curHit.index()];
783  curPoint =
784  curHit.hitPoint()
785  + offset(curHit.hitPoint(), curHit.index(), edgeVec);
786  }
787  }
788  while(true);
789  }
790  else
791  {
792  return findCellLinear(location);
793  }
794  }
795  return -1;
796 }
797 
798 
800 (
801  const point& location,
802  const label seedFaceI,
803  const bool useTreeSearch
804 ) const
805 {
806  if (seedFaceI == -1)
807  {
808  if (useTreeSearch)
809  {
810  const indexedOctree<treeDataFace>& tree = boundaryTree();
811 
812  scalar span = tree.bb().mag();
813 
814  pointIndexHit info = boundaryTree().findNearest
815  (
816  location,
817  Foam::sqr(span)
818  );
819 
820  if (!info.hit())
821  {
822  info = boundaryTree().findNearest
823  (
824  location,
825  Foam::sqr(GREAT)
826  );
827  }
828 
829  return tree.shapes().faceLabels()[info.index()];
830  }
831  else
832  {
833  scalar minDist = GREAT;
834 
835  label minFaceI = -1;
836 
837  for
838  (
839  label faceI = mesh_.nInternalFaces();
840  faceI < mesh_.nFaces();
841  faceI++
842  )
843  {
844  const face& f = mesh_.faces()[faceI];
845 
846  pointHit curHit =
847  f.nearestPoint
848  (
849  location,
850  mesh_.points()
851  );
852 
853  if (curHit.distance() < minDist)
854  {
855  minDist = curHit.distance();
856  minFaceI = faceI;
857  }
858  }
859  return minFaceI;
860  }
861  }
862  else
863  {
864  return findNearestBoundaryFaceWalk(location, seedFaceI);
865  }
866 }
867 
868 
870 (
871  const point& pStart,
872  const point& pEnd
873 ) const
874 {
875  pointIndexHit curHit = boundaryTree().findLine(pStart, pEnd);
876 
877  if (curHit.hit())
878  {
879  // Change index into octreeData into face label
880  curHit.setIndex(boundaryTree().shapes().faceLabels()[curHit.index()]);
881  }
882  return curHit;
883 }
884 
885 
887 (
888  const point& pStart,
889  const point& pEnd
890 ) const
891 {
893 
894  vector edgeVec = pEnd - pStart;
895  edgeVec /= mag(edgeVec);
896 
897  point pt = pStart;
898 
899  pointIndexHit bHit;
900  do
901  {
902  bHit = intersection(pt, pEnd);
903 
904  if (bHit.hit())
905  {
906  hits.append(bHit);
907 
908  const vector& area = mesh_.faceAreas()[bHit.index()];
909 
910  scalar typDim = Foam::sqrt(mag(area));
911 
912  if ((mag(bHit.hitPoint() - pEnd)/typDim) < SMALL)
913  {
914  break;
915  }
916 
917  // Restart from hitPoint shifted a little bit in the direction
918  // of the destination
919 
920  pt =
921  bHit.hitPoint()
922  + offset(bHit.hitPoint(), bHit.index(), edgeVec);
923  }
924 
925  } while(bHit.hit());
926 
927 
928  hits.shrink();
929 
930  return hits;
931 }
932 
933 
935 {
936  return
937  boundaryTree().getVolumeType(p)
939 }
940 
941 
942 // Delete all storage
944 {
945  deleteDemandDrivenData(boundaryTreePtr_);
946  deleteDemandDrivenData(cellTreePtr_);
947  deleteDemandDrivenData(cellCentreTreePtr_);
948 }
949 
950 
952 {
953  clearOut();
954 }
955 
956 
957 // ************************ vim: set sw=4 sts=4 et: ************************ //