distanceSurface.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2019 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 "distanceSurface.H"
28 
29 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
30 
31 namespace Foam
32 {
33 namespace sampledSurfaces
34 {
37 }
38 }
39 
40 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
41 
42 void Foam::sampledSurfaces::distanceSurface::createGeometry()
43 {
44  if (debug)
45  {
46  Pout<< "distanceSurface::createGeometry :updating geometry." << endl;
47  }
48 
49  // Clear any stored topologies
50  isoSurfPtr_.clear();
51 
52  // Clear derived data
53  clearGeom();
54 
55  const fvMesh& mesh = static_cast<const fvMesh&>(this->mesh());
56 
57  // Distance to cell centres
58  // ~~~~~~~~~~~~~~~~~~~~~~~~
59 
60  cellDistancePtr_.reset
61  (
62  new volScalarField
63  (
64  IOobject
65  (
66  "cellDistance",
67  mesh.time().timeName(),
68  mesh.time(),
71  false
72  ),
73  mesh,
75  )
76  );
77  volScalarField& cellDistance = cellDistancePtr_();
78 
79  // Internal field
80  {
81  const pointField& cc = mesh.C();
82  scalarField& fld = cellDistance.primitiveFieldRef();
83 
84  List<pointIndexHit> nearest;
85  surfPtr_().findNearest
86  (
87  cc,
88  scalarField(cc.size(), great),
89  nearest
90  );
91 
92  if (signed_)
93  {
94  List<volumeType> volType;
95 
96  surfPtr_().getVolumeType(cc, volType);
97 
98  forAll(volType, i)
99  {
100  volumeType vT = volType[i];
101 
102  if (vT == volumeType::outside)
103  {
104  fld[i] = Foam::mag(cc[i] - nearest[i].hitPoint());
105  }
106  else if (vT == volumeType::inside)
107  {
108  fld[i] = -Foam::mag(cc[i] - nearest[i].hitPoint());
109  }
110  else
111  {
113  << "getVolumeType failure, neither INSIDE or OUTSIDE"
114  << exit(FatalError);
115  }
116  }
117  }
118  else
119  {
120  forAll(nearest, i)
121  {
122  fld[i] = Foam::mag(cc[i] - nearest[i].hitPoint());
123  }
124  }
125  }
126 
127  volScalarField::Boundary& cellDistanceBf =
128  cellDistance.boundaryFieldRef();
129 
130  // Patch fields
131  {
132  forAll(mesh.C().boundaryField(), patchi)
133  {
134  const pointField& cc = mesh.C().boundaryField()[patchi];
135  fvPatchScalarField& fld = cellDistanceBf[patchi];
136 
137  List<pointIndexHit> nearest;
138  surfPtr_().findNearest
139  (
140  cc,
141  scalarField(cc.size(), great),
142  nearest
143  );
144 
145  if (signed_)
146  {
147  List<volumeType> volType;
148 
149  surfPtr_().getVolumeType(cc, volType);
150 
151  forAll(volType, i)
152  {
153  volumeType vT = volType[i];
154 
155  if (vT == volumeType::outside)
156  {
157  fld[i] = Foam::mag(cc[i] - nearest[i].hitPoint());
158  }
159  else if (vT == volumeType::inside)
160  {
161  fld[i] = -Foam::mag(cc[i] - nearest[i].hitPoint());
162  }
163  else
164  {
166  << "getVolumeType failure, "
167  << "neither INSIDE or OUTSIDE"
168  << exit(FatalError);
169  }
170  }
171  }
172  else
173  {
174  forAll(nearest, i)
175  {
176  fld[i] = Foam::mag(cc[i] - nearest[i].hitPoint());
177  }
178  }
179  }
180  }
181 
182 
183  // On processor patches the mesh.C() will already be the cell centre
184  // on the opposite side so no need to swap cellDistance.
185 
186 
187  // Distance to points
188  pointDistance_.setSize(mesh.nPoints());
189  {
190  const pointField& pts = mesh.points();
191 
192  List<pointIndexHit> nearest;
193  surfPtr_().findNearest
194  (
195  pts,
196  scalarField(pts.size(), great),
197  nearest
198  );
199 
200  if (signed_)
201  {
202  List<volumeType> volType;
203 
204  surfPtr_().getVolumeType(pts, volType);
205 
206  forAll(volType, i)
207  {
208  volumeType vT = volType[i];
209 
210  if (vT == volumeType::outside)
211  {
212  pointDistance_[i] =
213  Foam::mag(pts[i] - nearest[i].hitPoint());
214  }
215  else if (vT == volumeType::inside)
216  {
217  pointDistance_[i] =
218  -Foam::mag(pts[i] - nearest[i].hitPoint());
219  }
220  else
221  {
223  << "getVolumeType failure, neither INSIDE or OUTSIDE"
224  << exit(FatalError);
225  }
226  }
227  }
228  else
229  {
230  forAll(nearest, i)
231  {
232  pointDistance_[i] = Foam::mag(pts[i]-nearest[i].hitPoint());
233  }
234  }
235  }
236 
237 
238  if (debug)
239  {
240  Pout<< "Writing cell distance:" << cellDistance.objectPath() << endl;
241  cellDistance.write();
242  pointScalarField pDist
243  (
244  IOobject
245  (
246  "pointDistance",
247  mesh.time().timeName(),
248  mesh.time(),
251  false
252  ),
253  pointMesh::New(mesh),
255  );
256  pDist.primitiveFieldRef() = pointDistance_;
257 
258  Pout<< "Writing point distance:" << pDist.objectPath() << endl;
259  pDist.write();
260  }
261 
262 
263  //- Direct from cell field and point field.
264  isoSurfPtr_.reset
265  (
266  new isoSurface
267  (
268  mesh,
269  cellDistance,
270  pointDistance_,
271  distance_,
272  regularise_ ? isoSurface::DIAGCELL : isoSurface::NONE
273  )
274  );
275 
276  if (debug)
277  {
278  print(Pout);
279  Pout<< endl;
280  }
281 }
282 
283 
284 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
285 
287 (
288  const word& name,
289  const polyMesh& mesh,
290  const dictionary& dict
291 )
292 :
293  sampledSurface(name, mesh, dict),
294  surfPtr_
295  (
297  (
298  dict.lookup("surfaceType"),
299  IOobject
300  (
301  dict.lookupOrDefault("surfaceName", name), // name
302  mesh.time().constant(), // directory
303  "triSurface", // instance
304  mesh.time(), // registry
307  ),
308  dict
309  )
310  ),
311  distance_(readScalar(dict.lookup("distance"))),
312  signed_(readBool(dict.lookup("signed"))),
313  regularise_(dict.lookupOrDefault("regularise", true)),
314  average_(dict.lookupOrDefault("average", false)),
315  zoneKey_(keyType::null),
316  needsUpdate_(true),
317  isoSurfPtr_(nullptr)
318 {}
319 
320 
322 (
323  const word& name,
324  const polyMesh& mesh,
325  const bool interpolate,
326  const word& surfaceType,
327  const word& surfaceName,
328  const scalar distance,
329  const bool signedDistance,
330  const Switch regularise,
331  const Switch average
332 )
333 :
334  sampledSurface(name, mesh, interpolate),
335  surfPtr_
336  (
338  (
339  surfaceType,
340  IOobject
341  (
342  surfaceName, // name
343  mesh.time().constant(), // directory
344  "triSurface", // instance
345  mesh.time(), // registry
348  ),
349  dictionary()
350  )
351  ),
352  distance_(distance),
353  signed_(signedDistance),
354  regularise_(regularise),
355  average_(average),
356  zoneKey_(keyType::null),
357  needsUpdate_(true),
358  isoSurfPtr_(nullptr)
359 {}
360 
361 
362 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
363 
365 {}
366 
367 
368 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
369 
371 {
372  return needsUpdate_;
373 }
374 
375 
377 {
378  if (debug)
379  {
380  Pout<< "distanceSurface::expire :"
381  << " needsUpdate_:" << needsUpdate_ << endl;
382  }
383 
384  // Clear derived data
385  clearGeom();
386 
387  // already marked as expired
388  if (needsUpdate_)
389  {
390  return false;
391  }
392 
393  needsUpdate_ = true;
394  return true;
395 }
396 
397 
399 {
400  if (debug)
401  {
402  Pout<< "distanceSurface::update :"
403  << " needsUpdate_:" << needsUpdate_ << endl;
404  }
405 
406  if (!needsUpdate_)
407  {
408  return false;
409  }
410 
411  createGeometry();
412 
413  needsUpdate_ = false;
414  return true;
415 }
416 
417 
420 (
421  const volScalarField& vField
422 ) const
423 {
424  return sampleField(vField);
425 }
426 
427 
430 (
431  const volVectorField& vField
432 ) const
433 {
434  return sampleField(vField);
435 }
436 
437 
440 (
441  const volSphericalTensorField& vField
442 ) const
443 {
444  return sampleField(vField);
445 }
446 
447 
450 (
451  const volSymmTensorField& vField
452 ) const
453 {
454  return sampleField(vField);
455 }
456 
457 
460 (
461  const volTensorField& vField
462 ) const
463 {
464  return sampleField(vField);
465 }
466 
467 
470 (
471  const interpolation<scalar>& interpolator
472 ) const
473 {
474  return interpolateField(interpolator);
475 }
476 
477 
480 (
481  const interpolation<vector>& interpolator
482 ) const
483 {
484  return interpolateField(interpolator);
485 }
486 
489 (
490  const interpolation<sphericalTensor>& interpolator
491 ) const
492 {
493  return interpolateField(interpolator);
494 }
495 
496 
499 (
500  const interpolation<symmTensor>& interpolator
501 ) const
502 {
503  return interpolateField(interpolator);
504 }
505 
506 
509 (
510  const interpolation<tensor>& interpolator
511 ) const
512 {
513  return interpolateField(interpolator);
514 }
515 
516 
518 {
519  os << "distanceSurface: " << name() << " :"
520  << " surface:" << surfPtr_().name()
521  << " distance:" << distance_
522  << " faces:" << faces().size()
523  << " points:" << points().size();
524 }
525 
526 
527 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
const word & name() const
Return the name of this functionObject.
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:158
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
An abstract class for surfaces with sampling.
defineTypeNameAndDebug(distanceSurface, 0)
const Boundary & boundaryField() const
Return const-reference to the boundary field.
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: HashTable.H:59
virtual bool update()
Update the surface as required.
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:163
bool interpolate() const
Interpolation requested for surface.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
A simple wrapper around bool so that it can be read as a word: true/false, on/off, yes/no, y/n, t/f, or none/any.
Definition: Switch.H:60
Abstract base class with a fat-interface to all derived classes covering all possible ways in which t...
Definition: fvPatchField.H:66
static word timeName(const scalar, const int precision=precision_)
Return time name of given scalar time.
Definition: Time.C:626
bool readBool(Istream &)
Definition: boolIO.C:60
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:239
A sampledSurface defined by a distance to a surface.
Macros for easy insertion into run-time selection tables.
static const pointMesh & New(const polyMesh &mesh)
Definition: MeshObject.C:44
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1131
virtual void print(Ostream &) const
Write.
dynamicFvMesh & mesh
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().x()<< ' ';}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().y()<< ' ';}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().z()<< ' ';}gmvFile<< nl;forAll(lagrangianScalarNames, i){ const word &name=lagrangianScalarNames[i];IOField< scalar > fld(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
const pointField & points
A class for handling words, derived from string.
Definition: word.H:59
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
const word & constant() const
Return constant name.
Definition: TimePaths.H:124
virtual bool needsUpdate() const
Does the surface need an update?
virtual bool expire()
Mark the surface as needing an update.
A sampledSurface defined by a surface of iso value. To be used in sampleSurfaces / functionObjects...
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if successful.
Definition: doubleScalar.H:68
void reset(const label nPoints, const label nInternalFaces, const label nFaces, const label nCells)
Reset this primitiveMesh given the primitive array sizes.
An Ostream is an abstract base class for all output systems (streams, files, token lists...
Definition: Ostream.H:53
Internal::FieldType & primitiveFieldRef()
Return a reference to the internal field.
const Time & time() const
Return time.
addToRunTimeSelectionTable(sampledSurface, distanceSurface, word)
distanceSurface(const word &name, const polyMesh &mesh, const dictionary &dict)
Construct from dictionary.
label patchi
T lookupOrDefault(const word &, const T &, bool recursive=false, bool patternMatch=true) const
Find and return a T,.
Boundary & boundaryFieldRef()
Return a reference to the boundary field.
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:50
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
virtual tmp< scalarField > sample(const volScalarField &) const
Sample field on surface.
static autoPtr< searchableSurface > New(const word &surfaceType, const IOobject &io, const dictionary &dict)
Return a reference to the selected searchableSurface.
Set of surfaces to sample.
dimensioned< scalar > mag(const dimensioned< Type > &)
label nPoints() const
virtual Ostream & write(const token &)=0
Write next token to stream.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
const volVectorField & C() const
Return cell centres as volVectorField.
virtual bool write(const bool write=true) const
Write using setting from DB.
A class for managing temporary objects.
Definition: PtrList.H:53
static const keyType null
An empty keyType.
Definition: keyType.H:83
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
fileName objectPath() const
Return complete path + object name.
Definition: IOobject.H:404
Namespace for OpenFOAM.
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:583