sampledIsoSurface.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 "sampledIsoSurface.H"
27 #include "isoSurface.H"
29 
30 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
31 
32 namespace Foam
33 {
34 namespace sampledSurfaces
35 {
36  defineTypeNameAndDebug(isoSurface, 0);
37  addToRunTimeSelectionTable(sampledSurface, isoSurface, word);
38 }
39 }
40 
41 
42 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
43 
44 bool Foam::sampledSurfaces::isoSurface::updateGeometry() const
45 {
46  const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
47 
48  // no update needed
49  if (fvm.time().timeIndex() == prevTimeIndex_)
50  {
51  return false;
52  }
53 
54  prevTimeIndex_ = fvm.time().timeIndex();
55 
56  // Clear derived data
58 
59  // Optionally read volScalarField
60  autoPtr<volScalarField> readFieldPtr_;
61 
62  // 1. see if field in database
63  // 2. see if field can be read
64  const volScalarField* cellFldPtr = nullptr;
65  if (fvm.foundObject<volScalarField>(isoField_))
66  {
67  if (debug)
68  {
69  InfoInFunction << "Lookup " << isoField_ << endl;
70  }
71 
72  cellFldPtr = &fvm.lookupObject<volScalarField>(isoField_);
73  }
74  else
75  {
76  // Bit of a hack. Read field and store.
77 
78  if (debug)
79  {
81  << "Reading " << isoField_
82  << " from time " <<fvm.time().timeName()
83  << endl;
84  }
85 
86  readFieldPtr_.reset
87  (
88  new volScalarField
89  (
90  IOobject
91  (
92  isoField_,
93  fvm.time().timeName(),
94  fvm,
97  false
98  ),
99  fvm
100  )
101  );
102 
103  cellFldPtr = readFieldPtr_.operator->();
104  }
105  const volScalarField& cellFld = *cellFldPtr;
106 
107  tmp<pointScalarField> pointFld
108  (
110  );
111 
112  PtrList<Foam::isoSurface> isos(isoVals_.size());
113  forAll(isos, isoi)
114  {
115  isos.set
116  (
117  isoi,
118  new Foam::isoSurface
119  (
120  fvm,
121  cellFld.primitiveField(),
122  pointFld().primitiveField(),
123  isoVals_[isoi],
124  regularise_
127  )
128  );
129  }
130 
131  if (isos.size() == 1)
132  {
133  // Straight transfer
134  const_cast<isoSurface&>
135  (
136  *this
138  meshCells_ = isos[0].meshCells();
139  }
140  else
141  {
142  label nFaces = 0;
143  label nPoints = 0;
144  forAll(isos, isoi)
145  {
146  nFaces += isos[isoi].size();
147  nPoints += isos[isoi].points().size();
148  }
149 
150  faceList allFaces(nFaces);
151  labelList allMeshCells(nFaces);
152  pointField allPoints(nPoints);
153 
154  nFaces = 0;
155  nPoints = 0;
156  forAll(isos, isoi)
157  {
158  Foam::isoSurface& iso = isos[isoi];
159 
160  SubList<face> subAll(allFaces, iso.size(), nFaces);
161  subAll = iso;
162  // Offset point indices
163  if (nPoints > 0)
164  {
165  forAll(subAll, i)
166  {
167  face& f = subAll[i];
168  forAll(f, fp)
169  {
170  f[fp] += nPoints;
171  }
172  }
173  }
174  SubList<label>(allMeshCells, iso.size(), nFaces) = iso.meshCells();
175  nFaces += iso.size();
176 
177  const pointField& pts = iso.points();
178  SubList<point>(allPoints, pts.size(), nPoints) = pts;
179  nPoints += pts.size();
180 
181  // Clear out asap
182  isos.set(isoi, nullptr);
183  }
184 
185  if (nFaces != allFaces.size() || nPoints != allPoints.size())
186  {
187  FatalErrorInFunction << "nFaces:" << nFaces
188  << " nPoints:" << nPoints << exit(FatalError);
189  }
190 
191 
192  surfZoneList allZones(1);
193  allZones[0] = surfZone
194  (
195  "allFaces",
196  allFaces.size(), // size
197  0, // start
198  0 // index
199  );
200 
201  // Transfer
202  const_cast<isoSurface&>
203  (
204  *this
206  (
207  move(allPoints),
208  move(allFaces),
209  move(allZones)
210  );
211  meshCells_.transfer(allMeshCells);
212  }
213  if (debug)
214  {
215  Pout<< "sampledSurfaces::isoSurface::updateGeometry() : "
216  "constructed iso:"
217  << nl
218  << " regularise : " << regularise_ << nl
219  << " isoField : " << isoField_ << nl;
220  if (isoVals_.size() == 1)
221  {
222  Pout<< " isoValue : " << isoVals_[0] << nl;
223  }
224  else
225  {
226  Pout<< " isoValues : " << isoVals_ << nl;
227  }
228  Pout<< " points : " << points().size() << nl
229  << " faces : " << faces().size() << nl
230  << " cut cells : " << meshCells_.size() << endl;
231  }
232 
233  return true;
234 }
235 
236 
237 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
238 
240 (
241  const word& name,
242  const polyMesh& mesh,
243  const dictionary& dict
244 )
245 :
246  sampledSurface(name, mesh, dict),
247  isoField_(dict.lookup("isoField")),
248  isoVals_
249  (
250  dict.found("isoValues")
251  ? scalarField(dict.lookup("isoValues"))
252  : scalarField(1, readScalar(dict.lookup("isoValue")))
253  ),
254  regularise_(dict.lookupOrDefault("regularise", true)),
255  zoneKey_(keyType::null),
256  prevTimeIndex_(-1),
257  meshCells_(0)
258 {}
259 
260 
261 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
262 
264 {}
265 
266 
267 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
268 
270 {
271  const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
272 
273  return fvm.time().timeIndex() != prevTimeIndex_;
274 }
275 
276 
278 {
279  // Clear derived data
282 
283  // already marked as expired
284  if (prevTimeIndex_ == -1)
285  {
286  return false;
287  }
288 
289  // force update
290  prevTimeIndex_ = -1;
291  return true;
292 }
293 
294 
296 {
297  return updateGeometry();
298 }
299 
300 
303 (
304  const volScalarField& vField
305 ) const
306 {
307  return sampleField(vField);
308 }
309 
310 
313 (
314  const volVectorField& vField
315 ) const
316 {
317  return sampleField(vField);
318 }
319 
320 
323 (
324  const volSphericalTensorField& vField
325 ) const
326 {
327  return sampleField(vField);
328 }
329 
330 
333 (
334  const volSymmTensorField& vField
335 ) const
336 {
337  return sampleField(vField);
338 }
339 
340 
343 (
344  const volTensorField& vField
345 ) const
346 {
347  return sampleField(vField);
348 }
349 
350 
353 (
354  const interpolation<scalar>& interpolator
355 ) const
356 {
357  return interpolateField(interpolator);
358 }
359 
360 
363 (
364  const interpolation<vector>& interpolator
365 ) const
366 {
367  return interpolateField(interpolator);
368 }
369 
372 (
373  const interpolation<sphericalTensor>& interpolator
374 ) const
375 {
376  return interpolateField(interpolator);
377 }
378 
379 
382 (
383  const interpolation<symmTensor>& interpolator
384 ) const
385 {
386  return interpolateField(interpolator);
387 }
388 
389 
392 (
393  const interpolation<tensor>& interpolator
394 ) const
395 {
396  return interpolateField(interpolator);
397 }
398 
399 
401 {
402  os << "isoSurface: " << name() << " :"
403  << " field:" << isoField_;
404  if (isoVals_.size() == 1)
405  {
406  os << " value:" << isoVals_[0];
407  }
408  else
409  {
410  os << " values:" << isoVals_;
411  }
412  os << " faces:" << faces().size()
413  << " points:" << points().size();
414 }
415 
416 
417 // ************************************************************************* //
virtual tmp< scalarField > sample(const volScalarField &) const
Sample field on surface.
bool found(const word &, bool recursive=false, bool patternMatch=true) const
Search dictionary for given keyword.
Definition: dictionary.C:438
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
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
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
const word & name() const
Return the name of this functionObject.
virtual void print(Ostream &) const
Write.
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)
bool interpolate(const vector &p1, const vector &p2, const vector &o, vector &n, scalar l)
Definition: curveTools.C:75
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:163
bool interpolate() const
Interpolation requested for surface.
List< face > faceList
Definition: faceListFwd.H:43
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
virtual bool needsUpdate() const
Does the surface need an update?
void transfer(MeshedSurface< face > &)
Transfer the contents of the argument and annul the argument.
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:239
isoSurface(const word &name, const polyMesh &mesh, const dictionary &dict)
Construct from dictionary.
const labelList & meshCells() const
For every face original cell in mesh.
Definition: isoSurface.H:231
virtual bool expire()
Mark the surface as needing an update.
Macros for easy insertion into run-time selection tables.
static const volPointInterpolation & New(const fvMesh &mesh)
Definition: MeshObject.C:44
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:52
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
dynamicFvMesh & mesh
virtual void clearGeom() const
const pointField & points
A class for handling words, derived from string.
Definition: word.H:59
label nPoints
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
const Field< PointType > & points() const
Return reference to global points.
List< label > labelList
A List of labels.
Definition: labelList.H:56
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if successful.
Definition: doubleScalar.H:68
An Ostream is an abstract base class for all output systems (streams, files, token lists...
Definition: Ostream.H:53
virtual void reset(pointField &&points, List< face > &&faces, surfZoneList &&zones)
Reset primitive data (points, faces and zones)
static const char nl
Definition: Ostream.H:265
labelList f(nPoints)
addToRunTimeSelectionTable(sampledSurface, distanceSurface, word)
List< surfZone > surfZoneList
Definition: surfZoneList.H:45
label timeIndex() const
Return current time index.
Definition: TimeStateI.H:35
label size() const
The surface size is the number of faces.
T lookupOrDefault(const word &, const T &, bool recursive=false, bool patternMatch=true) const
Find and return a T,.
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
tmp< pointField > allPoints(const Triangulation &t)
Extract all points in vertex-index order.
virtual bool update()
Update the surface as required.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
Marching tet iso surface algorithm with optional filtering to keep only points originating from mesh ...
Definition: isoSurface.H:56
A class for managing temporary objects.
Definition: PtrList.H:53
static const keyType null
An empty keyType.
Definition: keyType.H:83
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
#define InfoInFunction
Report an information message using Foam::Info.