triSurfaceRegionSearch.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) 2011-2016 OpenFOAM Foundation
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 "triSurfaceRegionSearch.H"
27 #include "indexedOctree.H"
28 #include "triSurface.H"
29 #include "PatchTools.H"
30 
31 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
32 
33 Foam::triSurfaceRegionSearch::triSurfaceRegionSearch(const triSurface& surface)
34 :
35  triSurfaceSearch(surface),
36  indirectRegionPatches_(),
37  treeByRegion_()
38 {}
39 
40 
41 Foam::triSurfaceRegionSearch::triSurfaceRegionSearch
42 (
43  const triSurface& surface,
44  const dictionary& dict
45 )
46 :
47  triSurfaceSearch(surface, dict),
48  indirectRegionPatches_(),
49  treeByRegion_()
50 {}
51 
52 
53 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
54 
56 {
57  clearOut();
58 }
59 
60 
62 {
64  treeByRegion_.clear();
65 }
66 
67 
68 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
69 
72 {
73  if (treeByRegion_.empty())
74  {
75  Map<label> regionSizes;
76  forAll(surface(), fI)
77  {
78  const label regionI = surface()[fI].region();
79 
80  regionSizes(regionI)++;
81  }
82 
83  label nRegions = regionSizes.size();
84 
85  indirectRegionPatches_.setSize(nRegions);
86  treeByRegion_.setSize(nRegions);
87 
88  labelListList regionsAddressing(nRegions);
89 
90  forAll(regionsAddressing, regionI)
91  {
92  regionsAddressing[regionI] = labelList(regionSizes[regionI], -1);
93  }
94 
95  labelList nFacesInRegions(nRegions, 0);
96 
97  forAll(surface(), fI)
98  {
99  const label regionI = surface()[fI].region();
100 
101  regionsAddressing[regionI][nFacesInRegions[regionI]++] = fI;
102  }
103 
104  forAll(regionsAddressing, regionI)
105  {
106  scalar oldTol = treeType::perturbTol();
108 
109  indirectRegionPatches_.set
110  (
111  regionI,
113  (
115  (
116  surface(),
117  regionsAddressing[regionI]
118  ),
119  surface().points()
120  )
121  );
122 
123  // Calculate bb without constructing local point numbering.
124  treeBoundBox bb(Zero, Zero);
125 
126  if (indirectRegionPatches_[regionI].size())
127  {
128  label nPoints;
130  (
131  indirectRegionPatches_[regionI],
132  bb,
133  nPoints
134  );
135 
136  // if (nPoints != surface().points().size())
137  // {
138  // WarningInFunction
139  // << "Surface does not have compact point numbering. "
140  // << "Of " << surface().points().size()
141  // << " only " << nPoints
142  // << " are used."
143  // << " This might give problems in some routines."
144  // << endl;
145  // }
146 
147  // Random number generator. Bit dodgy since not exactly
148  // random ;-)
149  Random rndGen(65431);
150 
151  // Slightly extended bb. Slightly off-centred just so
152  // on symmetric geometry there are fewer face/edge
153  // aligned items.
154  bb = bb.extend(rndGen, 1e-4);
155  bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
156  bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
157  }
158 
159  treeByRegion_.set
160  (
161  regionI,
162  new treeType
163  (
165  (
166  true,
167  indirectRegionPatches_[regionI],
168  tolerance()
169  ),
170  bb,
171  maxTreeDepth(), // maxLevel
172  10, // leafsize
173  3.0 // duplicity
174  )
175  );
176 
177  treeType::perturbTol() = oldTol;
178  }
179  }
180 
181  return treeByRegion_;
182 }
183 
184 
186 (
187  const pointField& samples,
188  const scalarField& nearestDistSqr,
189  const labelList& regionIndices,
190  List<pointIndexHit>& info
191 ) const
192 {
193  if (regionIndices.empty())
194  {
195  triSurfaceSearch::findNearest(samples, nearestDistSqr, info);
196  }
197  else
198  {
199  scalar oldTol = treeType::perturbTol();
201 
202  const PtrList<treeType>& octrees = treeByRegion();
203 
204  info.setSize(samples.size());
205 
206  forAll(octrees, treeI)
207  {
208  if (findIndex(regionIndices, treeI) == -1)
209  {
210  continue;
211  }
212 
213  const treeType& octree = octrees[treeI];
214 
215  forAll(samples, i)
216  {
217 // if (!octree.bb().contains(samples[i]))
218 // {
219 // continue;
220 // }
221 
222  pointIndexHit currentRegionHit = octree.findNearest
223  (
224  samples[i],
225  nearestDistSqr[i],
227  );
228 
229  if
230  (
231  currentRegionHit.hit()
232  &&
233  (
234  !info[i].hit()
235  ||
236  (
237  magSqr(currentRegionHit.hitPoint() - samples[i])
238  < magSqr(info[i].hitPoint() - samples[i])
239  )
240  )
241  )
242  {
243  info[i] = currentRegionHit;
244  }
245  }
246  }
247 
248  treeType::perturbTol() = oldTol;
249  }
250 }
251 
252 
253 // ************************************************************************* //
label maxTreeDepth() const
Return max tree depth of octree.
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
bool empty() const
Return true if the UList is empty (ie, size() is zero)
Definition: UListI.H:313
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
const double e
Elementary charge.
Definition: doubleFloat.H:78
void findNearest(const pointField &samples, const scalarField &nearestDistSqr, const labelList &regionIndices, List< pointIndexHit > &info) const
Find the nearest point on the surface out of the regions.
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:163
static scalar & perturbTol()
Get the perturbation tolerance.
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
Definition: PointIndexHit.H:53
label size() const
Return number of elements in table.
Definition: HashTableI.H:65
void findNearest(const pointField &samples, const scalarField &nearestDistSqr, List< pointIndexHit > &info) const
Helper class to search on triSurface.
A list of faces which address into the list of points.
const Point & hitPoint() const
Return hit point.
const pointField & points
label nPoints
void clearOut()
Clear storage.
Encapsulation of data needed to search on PrimitivePatches.
cachedRandom rndGen(label(0), -1)
bool hit() const
Is there a hit.
static void calcBounds(const PrimitivePatch< Face, FaceList, PointField, PointType > &p, boundBox &bb, label &nPoints)
List< label > labelList
A List of labels.
Definition: labelList.H:56
static const zero Zero
Definition: zero.H:91
Simple random number generator.
Definition: Random.H:49
dimensioned< scalar > magSqr(const dimensioned< Type > &)
treeBoundBox extend(Random &, const scalar s) const
Return slightly wider bounding box.
const PtrList< treeType > & treeByRegion() const
Demand driven construction of octree for each region.
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurence of given element and return index,.
const triSurface & surface() const
Return reference to the surface.
void setSize(const label)
Reset size of List.
Definition: List.C:281
const point & max() const
Maximum describing the bounding box.
Definition: boundBoxI.H:60
vector point
Point is a vector.
Definition: point.H:41
Non-pointer based hierarchical recursive searching.
Definition: treeDataEdge.H:47
A templated 1D list of pointers to objects of type <T>, where the size of the array is known and used...
Definition: List.H:63
scalar tolerance() const
Return tolerance to use in searches.
Standard boundBox + extra functionality for use in octree.
Definition: treeBoundBox.H:87
const point & min() const
Minimum describing the bounding box.
Definition: boundBoxI.H:54
Triangulated surface description with patch information.
Definition: triSurface.H:65
A List with indirect addressing.
Definition: IndirectList.H:102