FreeFOAM The Cross-Platform CFD Toolkit
distanceSurface.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 "distanceSurface.H"
27 #include <OpenFOAM/dictionary.H>
28 #include <finiteVolume/volFields.H>
31 #include <finiteVolume/fvMesh.H>
32 #include <sampling/isoSurface.H>
33 // #include "isoSurfaceCell.H"
34 
35 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
36 
37 namespace Foam
38 {
39  defineTypeNameAndDebug(distanceSurface, 0);
40  addToRunTimeSelectionTable(sampledSurface, distanceSurface, word);
41 }
42 
43 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
44 
45 void Foam::distanceSurface::createGeometry()
46 {
47  if (debug)
48  {
49  Pout<< "distanceSurface::createGeometry :updating geometry." << endl;
50  }
51 
52  // Clear any stored topologies
53  facesPtr_.clear();
54 
55  // Clear derived data
56  clearGeom();
57 
58  const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
59 
60  // Distance to cell centres
61  // ~~~~~~~~~~~~~~~~~~~~~~~~
62 
63  cellDistancePtr_.reset
64  (
65  new volScalarField
66  (
67  IOobject
68  (
69  "cellDistance",
70  fvm.time().timeName(),
71  fvm.time(),
74  false
75  ),
76  fvm,
77  dimensionedScalar("zero", dimLength, 0)
78  )
79  );
80  volScalarField& cellDistance = cellDistancePtr_();
81 
82  // Internal field
83  {
84  const pointField& cc = fvm.C();
85  scalarField& fld = cellDistance.internalField();
86 
87  List<pointIndexHit> nearest;
88  surfPtr_().findNearest
89  (
90  cc,
91  scalarField(cc.size(), GREAT),
92  nearest
93  );
94 
95  if (signed_)
96  {
97  vectorField normal;
98  surfPtr_().getNormal(nearest, normal);
99 
100  forAll(nearest, i)
101  {
102  vector d(cc[i]-nearest[i].hitPoint());
103 
104  if ((d&normal[i]) > 0)
105  {
106  fld[i] = Foam::mag(d);
107  }
108  else
109  {
110  fld[i] = -Foam::mag(d);
111  }
112  }
113  }
114  else
115  {
116  forAll(nearest, i)
117  {
118  fld[i] = Foam::mag(cc[i] - nearest[i].hitPoint());
119  }
120  }
121  }
122 
123  // Patch fields
124  {
125  forAll(fvm.C().boundaryField(), patchI)
126  {
127  const pointField& cc = fvm.C().boundaryField()[patchI];
128  fvPatchScalarField& fld = cellDistance.boundaryField()[patchI];
129 
130  List<pointIndexHit> nearest;
131  surfPtr_().findNearest
132  (
133  cc,
134  scalarField(cc.size(), GREAT),
135  nearest
136  );
137 
138  if (signed_)
139  {
140  vectorField normal;
141  surfPtr_().getNormal(nearest, normal);
142 
143  forAll(nearest, i)
144  {
145  vector d(cc[i]-nearest[i].hitPoint());
146 
147  if ((d&normal[i]) > 0)
148  {
149  fld[i] = Foam::mag(d);
150  }
151  else
152  {
153  fld[i] = -Foam::mag(d);
154  }
155  }
156  }
157  else
158  {
159  forAll(nearest, i)
160  {
161  fld[i] = Foam::mag(cc[i] - nearest[i].hitPoint());
162  }
163  }
164  }
165  }
166 
167 
168  // On processor patches the mesh.C() will already be the cell centre
169  // on the opposite side so no need to swap cellDistance.
170 
171 
172  // Distance to points
173  pointDistance_.setSize(fvm.nPoints());
174  {
175  const pointField& pts = fvm.points();
176 
177  List<pointIndexHit> nearest;
178  surfPtr_().findNearest
179  (
180  pts,
181  scalarField(pts.size(), GREAT),
182  nearest
183  );
184 
185  if (signed_)
186  {
187  vectorField normal;
188  surfPtr_().getNormal(nearest, normal);
189 
190  forAll(nearest, i)
191  {
192  vector d(pts[i]-nearest[i].hitPoint());
193 
194  if ((d&normal[i]) > 0)
195  {
196  pointDistance_[i] = Foam::mag(d);
197  }
198  else
199  {
200  pointDistance_[i] = -Foam::mag(d);
201  }
202  }
203  }
204  else
205  {
206  forAll(nearest, i)
207  {
208  pointDistance_[i] = Foam::mag(pts[i]-nearest[i].hitPoint());
209  }
210  }
211  }
212 
213 
214  if (debug)
215  {
216  Pout<< "Writing cell distance:" << cellDistance.objectPath() << endl;
217  cellDistance.write();
218  pointScalarField pDist
219  (
220  IOobject
221  (
222  "pointDistance",
223  fvm.time().timeName(),
224  fvm.time(),
226  IOobject::NO_WRITE,
227  false
228  ),
229  pointMesh::New(fvm),
230  dimensionedScalar("zero", dimLength, 0)
231  );
232  pDist.internalField() = pointDistance_;
233 
234  Pout<< "Writing point distance:" << pDist.objectPath() << endl;
235  pDist.write();
236  }
237 
238 
239  //- Direct from cell field and point field.
240  isoSurfPtr_.reset
241  (
242  new isoSurface
243  (
244  cellDistance,
245  pointDistance_,
246  distance_,
247  regularise_
248  )
249  );
250 
251  if (debug)
252  {
253  print(Pout);
254  Pout<< endl;
255  }
256 }
257 
258 
259 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
260 
262 (
263  const word& name,
264  const polyMesh& mesh,
265  const dictionary& dict
266 )
267 :
268  sampledSurface(name, mesh, dict),
269  surfPtr_
270  (
272  (
273  dict.lookup("surfaceType"),
274  IOobject
275  (
276  dict.lookupOrDefault("surfaceName", name), // name
277  mesh.time().constant(), // directory
278  "triSurface", // instance
279  mesh.time(), // registry
282  ),
283  dict
284  )
285  ),
286  distance_(readScalar(dict.lookup("distance"))),
287  signed_(readBool(dict.lookup("signed"))),
288  regularise_(dict.lookupOrDefault("regularise", true)),
289  zoneName_(word::null),
290  needsUpdate_(true),
291  isoSurfPtr_(NULL),
292  facesPtr_(NULL)
293 {
294 // dict.readIfPresent("zone", zoneName_);
295 //
296 // if (debug && zoneName_.size())
297 // {
298 // if (mesh.cellZones().findZoneID(zoneName_) < 0)
299 // {
300 // Info<< "cellZone \"" << zoneName_
301 // << "\" not found - using entire mesh" << endl;
302 // }
303 // }
304 }
305 
306 
307 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
308 
310 {}
311 
312 
313 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
314 
316 {
317  return needsUpdate_;
318 }
319 
320 
322 {
323  if (debug)
324  {
325  Pout<< "distanceSurface::expire :"
326  << " have-facesPtr_:" << facesPtr_.valid()
327  << " needsUpdate_:" << needsUpdate_ << endl;
328  }
329 
330  // Clear any stored topologies
331  facesPtr_.clear();
332 
333  // Clear derived data
334  clearGeom();
335 
336  // already marked as expired
337  if (needsUpdate_)
338  {
339  return false;
340  }
341 
342  needsUpdate_ = true;
343  return true;
344 }
345 
346 
348 {
349  if (debug)
350  {
351  Pout<< "distanceSurface::update :"
352  << " have-facesPtr_:" << facesPtr_.valid()
353  << " needsUpdate_:" << needsUpdate_ << endl;
354  }
355 
356  if (!needsUpdate_)
357  {
358  return false;
359  }
360 
361  createGeometry();
362 
363  needsUpdate_ = false;
364  return true;
365 }
366 
367 
370 (
371  const volScalarField& vField
372 ) const
373 {
374  return sampleField(vField);
375 }
376 
377 
380 (
381  const volVectorField& vField
382 ) const
383 {
384  return sampleField(vField);
385 }
386 
387 
390 (
391  const volSphericalTensorField& vField
392 ) const
393 {
394  return sampleField(vField);
395 }
396 
397 
400 (
401  const volSymmTensorField& vField
402 ) const
403 {
404  return sampleField(vField);
405 }
406 
407 
410 (
411  const volTensorField& vField
412 ) const
413 {
414  return sampleField(vField);
415 }
416 
417 
420 (
421  const interpolation<scalar>& interpolator
422 ) const
423 {
424  return interpolateField(interpolator);
425 }
426 
427 
430 (
431  const interpolation<vector>& interpolator
432 ) const
433 {
434  return interpolateField(interpolator);
435 }
436 
439 (
440  const interpolation<sphericalTensor>& interpolator
441 ) const
442 {
443  return interpolateField(interpolator);
444 }
445 
446 
449 (
450  const interpolation<symmTensor>& interpolator
451 ) const
452 {
453  return interpolateField(interpolator);
454 }
455 
456 
459 (
460  const interpolation<tensor>& interpolator
461 ) const
462 {
463  return interpolateField(interpolator);
464 }
465 
466 
468 {
469  os << "distanceSurface: " << name() << " :"
470  << " surface:" << surfPtr_().name()
471  << " distance:" << distance_
472  << " faces:" << faces().size()
473  << " points:" << points().size();
474 }
475 
476 
477 // ************************ vim: set sw=4 sts=4 et: ************************ //