FreeFOAM The Cross-Platform CFD Toolkit
patchCloudSet.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-2011 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 "patchCloudSet.H"
27 #include <OpenFOAM/polyMesh.H>
29 #include <meshTools/treeBoundBox.H>
30 #include <meshTools/treeDataFace.H>
31 #include <OpenFOAM/Time.H>
32 #include <meshTools/meshTools.H>
33 #include <OpenFOAM/wordReList.H>
34 // For 'nearInfo' helper class only
36 
37 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
38 
39 namespace Foam
40 {
41  defineTypeNameAndDebug(patchCloudSet, 0);
42  addToRunTimeSelectionTable(sampledSet, patchCloudSet, word);
43 }
44 
45 
46 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
47 
48 void Foam::patchCloudSet::calcSamples
49 (
50  DynamicList<point>& samplingPts,
51  DynamicList<label>& samplingCells,
52  DynamicList<label>& samplingFaces,
53  DynamicList<label>& samplingSegments,
54  DynamicList<scalar>& samplingCurveDist
55 ) const
56 {
57  if (debug)
58  {
59  Info<< "patchCloudSet : sampling on patches :" << endl;
60  }
61 
62  // Construct search tree for all patch faces.
63  label sz = 0;
64  forAllConstIter(labelHashSet, patchSet_, iter)
65  {
66  const polyPatch& pp = mesh().boundaryMesh()[iter.key()];
67 
68  sz += pp.size();
69 
70  if (debug)
71  {
72  Info<< " " << pp.name() << " size " << pp.size() << endl;
73  }
74  }
75 
76  labelList patchFaces(sz);
77  sz = 0;
78  treeBoundBox bb(point::max, point::min);
79  forAllConstIter(labelHashSet, patchSet_, iter)
80  {
81  const polyPatch& pp = mesh().boundaryMesh()[iter.key()];
82 
83  forAll(pp, i)
84  {
85  patchFaces[sz++] = pp.start()+i;
86  }
87 
88  // Do not do reduction.
89  const boundBox patchBb(pp.localPoints(), false);
90 
91  bb.min() = min(bb.min(), patchBb.min());
92  bb.max() = max(bb.max(), patchBb.max());
93  }
94 
95  // Not very random
96  Random rndGen(123456);
97  // Make bb asymetric just to avoid problems on symmetric meshes
98  bb = bb.extend(rndGen, 1E-4);
99 
100  // Make sure bb is 3D.
101  bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
102  bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
103 
104 
105  indexedOctree<treeDataFace> patchTree
106  (
107  treeDataFace // all information needed to search faces
108  (
109  false, // do not cache bb
110  mesh(),
111  patchFaces // boundary faces only
112  ),
113  bb, // overall search domain
114  8, // maxLevel
115  10, // leafsize
116  3.0 // duplicity
117  );
118 
119 
120 
121  // All the info for nearest. Construct to miss
122  List<directMappedPatchBase::nearInfo> nearest(sampleCoords_.size());
123 
124  forAll(sampleCoords_, sampleI)
125  {
126  const point& sample = sampleCoords_[sampleI];
127 
128  pointIndexHit& nearInfo = nearest[sampleI].first();
129 
130  // Find the nearest locally
131  if (patchFaces.size())
132  {
133  nearInfo = patchTree.findNearest(sample, sqr(searchDist_));
134  }
135  else
136  {
137  nearInfo.setMiss();
138  }
139 
140 
141  // Fill in the distance field and the processor field
142  if (!nearInfo.hit())
143  {
144  nearest[sampleI].second().first() = Foam::sqr(GREAT);
145  nearest[sampleI].second().second() = Pstream::myProcNo();
146  }
147  else
148  {
149  // Set nearest to mesh face label
150  nearInfo.setIndex(patchFaces[nearInfo.index()]);
151 
152  nearest[sampleI].second().first() = magSqr
153  (
154  nearInfo.hitPoint()
155  - sample
156  );
157  nearest[sampleI].second().second() = Pstream::myProcNo();
158  }
159  }
160 
161 
162  // Find nearest.
163  Pstream::listCombineGather(nearest, directMappedPatchBase::nearestEqOp());
165 
166 
167  if (debug && Pstream::master())
168  {
169  OFstream str
170  (
171  mesh().time().path()
172  / name()
173  + "_nearest.obj"
174  );
175  Info<< "Dumping mapping as lines from supplied points to"
176  << " nearest patch face to file " << str.name() << endl;
177 
178  label vertI = 0;
179 
180  forAll(nearest, i)
181  {
182  if (nearest[i].first().hit())
183  {
184  meshTools::writeOBJ(str, sampleCoords_[i]);
185  vertI++;
186  meshTools::writeOBJ(str, nearest[i].first().hitPoint());
187  vertI++;
188  str << "l " << vertI-1 << ' ' << vertI << nl;
189  }
190  }
191  }
192 
193 
194  // Store the sampling locations on the nearest processor
195  forAll(nearest, sampleI)
196  {
197  const pointIndexHit& nearInfo = nearest[sampleI].first();
198 
199  if (nearInfo.hit())
200  {
201  if (nearest[sampleI].second().second() == Pstream::myProcNo())
202  {
203  label faceI = nearInfo.index();
204 
205  samplingPts.append(nearInfo.hitPoint());
206  samplingCells.append(mesh().faceOwner()[faceI]);
207  samplingFaces.append(faceI);
208  samplingSegments.append(0);
209  samplingCurveDist.append(1.0 * sampleI);
210  }
211  }
212  else
213  {
214  // No processor found point near enough. Mark with special value
215  // which is intercepted when interpolating
216  if (Pstream::master())
217  {
218  samplingPts.append(sampleCoords_[sampleI]);
219  samplingCells.append(-1);
220  samplingFaces.append(-1);
221  samplingSegments.append(0);
222  samplingCurveDist.append(1.0 * sampleI);
223  }
224  }
225  }
226 }
227 
228 
229 void Foam::patchCloudSet::genSamples()
230 {
231  // Storage for sample points
232  DynamicList<point> samplingPts;
233  DynamicList<label> samplingCells;
234  DynamicList<label> samplingFaces;
235  DynamicList<label> samplingSegments;
236  DynamicList<scalar> samplingCurveDist;
237 
238  calcSamples
239  (
240  samplingPts,
241  samplingCells,
242  samplingFaces,
243  samplingSegments,
244  samplingCurveDist
245  );
246 
247  samplingPts.shrink();
248  samplingCells.shrink();
249  samplingFaces.shrink();
250  samplingSegments.shrink();
251  samplingCurveDist.shrink();
252 
253  setSamples
254  (
255  samplingPts,
256  samplingCells,
257  samplingFaces,
258  samplingSegments,
259  samplingCurveDist
260  );
261 }
262 
263 
264 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
265 
267 (
268  const word& name,
269  const polyMesh& mesh,
270  meshSearch& searchEngine,
271  const word& axis,
272  const List<point>& sampleCoords,
273  const labelHashSet& patchSet,
274  const scalar searchDist
275 )
276 :
277  sampledSet(name, mesh, searchEngine, axis),
278  sampleCoords_(sampleCoords),
279  patchSet_(patchSet),
280  searchDist_(searchDist)
281 {
282  genSamples();
283 
284  if (debug)
285  {
286  write(Info);
287  }
288 }
289 
290 
292 (
293  const word& name,
294  const polyMesh& mesh,
295  meshSearch& searchEngine,
296  const dictionary& dict
297 )
298 :
299  sampledSet(name, mesh, searchEngine, dict),
300  sampleCoords_(dict.lookup("points")),
301  patchSet_
302  (
303  mesh.boundaryMesh().patchSet
304  (
305  wordList(dict.lookup("patches"))
306  )
307  ),
308  searchDist_(readScalar(dict.lookup("maxDistance")))
309 {
310  genSamples();
311 
312  if (debug)
313  {
314  write(Info);
315  }
316 }
317 
318 
319 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
320 
322 {}
323 
324 
325 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
326 
328 {
329  if (pts.size())
330  {
331  // Use first samplePt as starting point
332  return pts[0];
333  }
334  else
335  {
336  return vector::zero;
337  }
338 }
339 
340 
341 // ************************ vim: set sw=4 sts=4 et: ************************ //