FreeFOAM The Cross-Platform CFD Toolkit
sampledTriSurfaceMesh.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) 2010-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 "sampledTriSurfaceMesh.H"
28 #include <meshTools/meshSearch.H>
29 #include <OpenFOAM/Tuple2.H>
30 #include <OpenFOAM/globalIndex.H>
31 
33 
34 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
35 
36 namespace Foam
37 {
38  defineTypeNameAndDebug(sampledTriSurfaceMesh, 0);
40  (
41  sampledSurface,
42  sampledTriSurfaceMesh,
43  word
44  );
45 
46  //- Private class for finding nearest
47  // - global index
48  // - sqr(distance)
50 
52  {
53 
54  public:
55 
56  void operator()(nearInfo& x, const nearInfo& y) const
57  {
58  if (y.first() < x.first())
59  {
60  x = y;
61  }
62  }
63  };
64 }
65 
66 
67 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
68 
70 (
71  const word& name,
72  const polyMesh& mesh,
73  const word& surfaceName
74 )
75 :
76  sampledSurface(name, mesh),
77  surface_
78  (
79  IOobject
80  (
81  name,
82  mesh.time().constant(), // instance
83  "triSurface", // local
84  mesh, // registry
85  IOobject::MUST_READ,
86  IOobject::NO_WRITE,
87  false
88  )
89  ),
90  needsUpdate_(true),
91  cellLabels_(0),
92  pointToFace_(0)
93 {}
94 
95 
97 (
98  const word& name,
99  const polyMesh& mesh,
100  const dictionary& dict
101 )
102 :
103  sampledSurface(name, mesh, dict),
104  surface_
105  (
106  IOobject
107  (
108  dict.lookup("surface"),
109  mesh.time().constant(), // instance
110  "triSurface", // local
111  mesh, // registry
112  IOobject::MUST_READ,
113  IOobject::NO_WRITE,
114  false
115  )
116  ),
117  needsUpdate_(true),
118  cellLabels_(0),
119  pointToFace_(0)
120 {}
121 
122 
123 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
124 
126 {}
127 
128 
129 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
130 
132 {
133  return needsUpdate_;
134 }
135 
136 
138 {
139  // already marked as expired
140  if (needsUpdate_)
141  {
142  return false;
143  }
144 
147  cellLabels_.clear();
148  pointToFace_.clear();
149 
150  needsUpdate_ = true;
151  return true;
152 }
153 
154 
156 {
157  if (!needsUpdate_)
158  {
159  return false;
160  }
161 
162 
163  // Find the cells the triangles of the surface are in.
164  // Does approximation by looking at the face centres only
165  const pointField& fc = surface_.faceCentres();
166 
167  meshSearch meshSearcher(mesh(), false);
168 
169  const indexedOctree<treeDataPoint>& cellCentreTree =
170  meshSearcher.cellCentreTree();
171 
172 
173  // Global numbering for cells - only used to uniquely identify local cells.
174  globalIndex globalCells(mesh().nCells());
175  List<nearInfo> nearest(fc.size());
176  forAll(nearest, i)
177  {
178  nearest[i].first() = GREAT;
179  nearest[i].second() = labelMax;
180  }
181 
182  // Search triangles using nearest on local mesh
183  forAll(fc, triI)
184  {
185  pointIndexHit nearInfo = cellCentreTree.findNearest
186  (
187  fc[triI],
188  sqr(GREAT)
189  );
190  if (nearInfo.hit())
191  {
192  nearest[triI].first() = magSqr(nearInfo.hitPoint()-fc[triI]);
193  nearest[triI].second() = globalCells.toGlobal(nearInfo.index());
194  }
195  }
196 
197  // See which processor has the nearest.
200 
201  boolList include(surface_.size(), false);
202 
203  cellLabels_.setSize(fc.size());
204  cellLabels_ = -1;
205 
206  label nFound = 0;
207  forAll(nearest, triI)
208  {
209  if (nearest[triI].second() == labelMax)
210  {
211  // Not found on any processor. How to map?
212  }
213  else if (globalCells.isLocal(nearest[triI].second()))
214  {
215  cellLabels_[triI] = globalCells.toLocal(nearest[triI].second());
216 
217  include[triI] = true;
218  nFound++;
219  }
220  }
221 
222 
223  if (debug)
224  {
225  Pout<< "Local out of faces:" << cellLabels_.size()
226  << " keeping:" << nFound << endl;
227  }
228 
229  // Now subset the surface. Do not use triSurface::subsetMesh since requires
230  // original surface to be in compact numbering.
231 
232  const triSurface& s = surface_;
233 
234  // Compact to original triangle
235  labelList faceMap(s.size());
236  // Compact to original points
237  labelList pointMap(s.points().size());
238  // From original point to compact points
239  labelList reversePointMap(s.points().size(), -1);
240 
241  {
242  label newPointI = 0;
243  label newTriI = 0;
244 
245  forAll(s, triI)
246  {
247  if (include[triI])
248  {
249  faceMap[newTriI++] = triI;
250 
251  const labelledTri& f = s[triI];
252  forAll(f, fp)
253  {
254  if (reversePointMap[f[fp]] == -1)
255  {
256  pointMap[newPointI] = f[fp];
257  reversePointMap[f[fp]] = newPointI++;
258  }
259  }
260  }
261  }
262  faceMap.setSize(newTriI);
263  pointMap.setSize(newPointI);
264  }
265 
266  // Subset cellLabels
267  cellLabels_ = UIndirectList<label>(cellLabels_, faceMap)();
268 
269  // Store any face per point
270  pointToFace_.setSize(pointMap.size());
271 
272  // Create faces and points for subsetted surface
273  faceList& faces = this->storedFaces();
274  faces.setSize(faceMap.size());
275  forAll(faceMap, i)
276  {
277  const triFace& f = s[faceMap[i]];
278  triFace newF
279  (
280  reversePointMap[f[0]],
281  reversePointMap[f[1]],
282  reversePointMap[f[2]]
283  );
284  faces[i] = newF.triFaceFace();
285 
286  forAll(newF, fp)
287  {
288  pointToFace_[newF[fp]] = i;
289  }
290  }
291 
292  this->storedPoints() = pointField(s.points(), pointMap);
293 
294  if (debug)
295  {
296  print(Pout);
297  Pout<< endl;
298  }
299 
300  needsUpdate_ = false;
301  return true;
302 }
303 
304 
307 (
308  const volScalarField& vField
309 ) const
310 {
311  return sampleField(vField);
312 }
313 
314 
317 (
318  const volVectorField& vField
319 ) const
320 {
321  return sampleField(vField);
322 }
323 
326 (
327  const volSphericalTensorField& vField
328 ) const
329 {
330  return sampleField(vField);
331 }
332 
333 
336 (
337  const volSymmTensorField& vField
338 ) const
339 {
340  return sampleField(vField);
341 }
342 
343 
346 (
347  const volTensorField& vField
348 ) const
349 {
350  return sampleField(vField);
351 }
352 
353 
356 (
357  const interpolation<scalar>& interpolator
358 ) const
359 {
360  return interpolateField(interpolator);
361 }
362 
363 
366 (
367  const interpolation<vector>& interpolator
368 ) const
369 {
370  return interpolateField(interpolator);
371 }
372 
375 (
376  const interpolation<sphericalTensor>& interpolator
377 ) const
378 {
379  return interpolateField(interpolator);
380 }
381 
382 
385 (
386  const interpolation<symmTensor>& interpolator
387 ) const
388 {
389  return interpolateField(interpolator);
390 }
391 
392 
395 (
396  const interpolation<tensor>& interpolator
397 ) const
398 {
399  return interpolateField(interpolator);
400 }
401 
402 
404 {
405  os << "sampledTriSurfaceMesh: " << name() << " :"
406  << " surface:" << surface_.objectRegistry::name()
407  << " faces:" << faces().size()
408  << " points:" << points().size();
409 }
410 
411 
412 // ************************ vim: set sw=4 sts=4 et: ************************ //