sampledCuttingPlane.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-2018 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 "sampledCuttingPlane.H"
27 #include "dictionary.H"
28 #include "volFields.H"
29 #include "volPointInterpolation.H"
31 #include "fvMesh.H"
32 
33 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
34 
35 namespace Foam
36 {
37  defineTypeNameAndDebug(sampledCuttingPlane, 0);
39  (
40  sampledSurface,
41  sampledCuttingPlane,
42  word,
43  cuttingPlane
44  );
45 }
46 
47 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
48 
49 void Foam::sampledCuttingPlane::createGeometry()
50 {
51  if (debug)
52  {
53  Pout<< "sampledCuttingPlane::createGeometry :updating geometry."
54  << endl;
55  }
56 
57  // Clear any stored topologies
58  facesPtr_.clear();
59  isoSurfPtr_.ptr();
60  pointDistance_.clear();
61  cellDistancePtr_.clear();
62 
63  // Clear derived data
64  clearGeom();
65 
66  // Get any subMesh
67  if (zoneID_.index() != -1 && !subMeshPtr_.valid())
68  {
69  const polyBoundaryMesh& patches = mesh().boundaryMesh();
70 
71  // Patch to put exposed internal faces into
72  const label exposedPatchi = patches.findPatchID(exposedPatchName_);
73 
74  DebugInfo
75  << "Allocating subset of size "
76  << mesh().cellZones()[zoneID_.index()].size()
77  << " with exposed faces into patch "
78  << patches[exposedPatchi].name() << endl;
79 
80  subMeshPtr_.reset
81  (
82  new fvMeshSubset(static_cast<const fvMesh&>(mesh()))
83  );
84  subMeshPtr_().setLargeCellSubset
85  (
86  labelHashSet(mesh().cellZones()[zoneID_.index()]),
87  exposedPatchi
88  );
89  }
90 
91 
92  // Select either the submesh or the underlying mesh
93  const fvMesh& mesh =
94  (
95  subMeshPtr_.valid()
96  ? subMeshPtr_().subMesh()
97  : static_cast<const fvMesh&>(this->mesh())
98  );
99 
100 
101  // Distance to cell centres
102  // ~~~~~~~~~~~~~~~~~~~~~~~~
103 
104  cellDistancePtr_.reset
105  (
106  new volScalarField
107  (
108  IOobject
109  (
110  "cellDistance",
111  mesh.time().timeName(),
112  mesh.time(),
115  false
116  ),
117  mesh,
118  dimensionedScalar("zero", dimLength, 0)
119  )
120  );
121  volScalarField& cellDistance = cellDistancePtr_();
122 
123  // Internal field
124  {
125  const pointField& cc = mesh.cellCentres();
126  scalarField& fld = cellDistance.primitiveFieldRef();
127 
128  forAll(cc, i)
129  {
130  // Signed distance
131  fld[i] = (cc[i] - plane_.refPoint()) & plane_.normal();
132  }
133  }
134 
135  volScalarField::Boundary& cellDistanceBf =
136  cellDistance.boundaryFieldRef();
137 
138  // Patch fields
139  {
140  forAll(cellDistanceBf, patchi)
141  {
142  if
143  (
144  isA<emptyFvPatchScalarField>
145  (
146  cellDistanceBf[patchi]
147  )
148  )
149  {
150  cellDistanceBf.set
151  (
152  patchi,
153  new calculatedFvPatchScalarField
154  (
155  mesh.boundary()[patchi],
156  cellDistance
157  )
158  );
159 
160  const polyPatch& pp = mesh.boundary()[patchi].patch();
161  pointField::subField cc = pp.patchSlice(mesh.faceCentres());
162 
163  fvPatchScalarField& fld = cellDistanceBf[patchi];
164  fld.setSize(pp.size());
165  forAll(fld, i)
166  {
167  fld[i] = (cc[i] - plane_.refPoint()) & plane_.normal();
168  }
169  }
170  else
171  {
172  const pointField& cc = mesh.C().boundaryField()[patchi];
173  fvPatchScalarField& fld = cellDistanceBf[patchi];
174 
175  forAll(fld, i)
176  {
177  fld[i] = (cc[i] - plane_.refPoint()) & plane_.normal();
178  }
179  }
180  }
181  }
182 
183 
184  // On processor patches the mesh.C() will already be the cell centre
185  // on the opposite side so no need to swap cellDistance.
186 
187 
188  // Distance to points
189  pointDistance_.setSize(mesh.nPoints());
190  {
191  const pointField& pts = mesh.points();
192 
193  forAll(pointDistance_, i)
194  {
195  pointDistance_[i] = (pts[i] - plane_.refPoint()) & plane_.normal();
196  }
197  }
198 
199 
200  if (debug)
201  {
202  Pout<< "Writing cell distance:" << cellDistance.objectPath() << endl;
203  cellDistance.write();
204  pointScalarField pDist
205  (
206  IOobject
207  (
208  "pointDistance",
209  mesh.time().timeName(),
210  mesh.time(),
213  false
214  ),
215  pointMesh::New(mesh),
216  dimensionedScalar("zero", dimLength, 0)
217  );
218  pDist.primitiveFieldRef() = pointDistance_;
219 
220  Pout<< "Writing point distance:" << pDist.objectPath() << endl;
221  pDist.write();
222  }
223 
224 
225  //- Direct from cell field and point field.
226  isoSurfPtr_.reset
227  (
228  new isoSurface
229  (
230  cellDistance,
231  pointDistance_,
232  0.0,
233  regularise_,
234  mergeTol_
235  )
236  // new isoSurfaceCell
237  //(
238  // mesh,
239  // cellDistance,
240  // pointDistance_,
241  // 0.0,
242  // regularise_,
243  // mergeTol_
244  //)
245  );
246 
247  if (debug)
248  {
249  print(Pout);
250  Pout<< endl;
251  }
252 }
253 
254 
255 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
256 
258 (
259  const word& name,
260  const polyMesh& mesh,
261  const dictionary& dict
262 )
263 :
264  sampledSurface(name, mesh, dict),
265  plane_(dict),
266  mergeTol_(dict.lookupOrDefault("mergeTol", 1e-6)),
267  regularise_(dict.lookupOrDefault("regularise", true)),
268  average_(dict.lookupOrDefault("average", false)),
269  zoneID_(dict.lookupOrDefault("zone", word::null), mesh.cellZones()),
270  exposedPatchName_(word::null),
271  needsUpdate_(true),
272  subMeshPtr_(nullptr),
273  cellDistancePtr_(nullptr),
274  isoSurfPtr_(nullptr),
275  facesPtr_(nullptr)
276 {
277  if (zoneID_.index() != -1)
278  {
279  dict.lookup("exposedPatchName") >> exposedPatchName_;
280 
281  if (mesh.boundaryMesh().findPatchID(exposedPatchName_) == -1)
282  {
284  << "Cannot find patch " << exposedPatchName_
285  << " in which to put exposed faces." << endl
286  << "Valid patches are " << mesh.boundaryMesh().names()
287  << exit(FatalError);
288  }
289 
290  if (debug && zoneID_.index() != -1)
291  {
292  Info<< "Restricting to cellZone " << zoneID_.name()
293  << " with exposed internal faces into patch "
294  << exposedPatchName_ << endl;
295  }
296  }
297 }
298 
299 
300 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
301 
303 {}
304 
305 
306 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
307 
309 {
310  return needsUpdate_;
311 }
312 
313 
315 {
316  if (debug)
317  {
318  Pout<< "sampledCuttingPlane::expire :"
319  << " have-facesPtr_:" << facesPtr_.valid()
320  << " needsUpdate_:" << needsUpdate_ << endl;
321  }
322 
323  // Clear any stored topologies
324  facesPtr_.clear();
325 
326  // Clear derived data
327  clearGeom();
328 
329  // Already marked as expired
330  if (needsUpdate_)
331  {
332  return false;
333  }
334 
335  needsUpdate_ = true;
336  return true;
337 }
338 
339 
341 {
342  if (debug)
343  {
344  Pout<< "sampledCuttingPlane::update :"
345  << " have-facesPtr_:" << facesPtr_.valid()
346  << " needsUpdate_:" << needsUpdate_ << endl;
347  }
348 
349  if (!needsUpdate_)
350  {
351  return false;
352  }
353 
354  createGeometry();
355 
356  needsUpdate_ = false;
357  return true;
358 }
359 
360 
363 (
364  const volScalarField& vField
365 ) const
366 {
367  return sampleField(vField);
368 }
369 
370 
373 (
374  const volVectorField& vField
375 ) const
376 {
377  return sampleField(vField);
378 }
379 
380 
383 (
384  const volSphericalTensorField& vField
385 ) const
386 {
387  return sampleField(vField);
388 }
389 
390 
393 (
394  const volSymmTensorField& vField
395 ) const
396 {
397  return sampleField(vField);
398 }
399 
400 
403 (
404  const volTensorField& vField
405 ) const
406 {
407  return sampleField(vField);
408 }
409 
410 
413 (
414  const interpolation<scalar>& interpolator
415 ) const
416 {
417  return interpolateField(interpolator);
418 }
419 
420 
423 (
424  const interpolation<vector>& interpolator
425 ) const
426 {
427  return interpolateField(interpolator);
428 }
429 
432 (
433  const interpolation<sphericalTensor>& interpolator
434 ) const
435 {
436  return interpolateField(interpolator);
437 }
438 
439 
442 (
443  const interpolation<symmTensor>& interpolator
444 ) const
445 {
446  return interpolateField(interpolator);
447 }
448 
449 
452 (
453  const interpolation<tensor>& interpolator
454 ) const
455 {
456  return interpolateField(interpolator);
457 }
458 
459 
461 {
462  os << "sampledCuttingPlane: " << name() << " :"
463  << " plane:" << plane_
464  << " faces:" << faces().size()
465  << " points:" << points().size();
466 }
467 
468 
469 // ************************************************************************* //
GeometricField< scalar, pointPatchField, pointMesh > pointScalarField
Definition: pointFields.H:49
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:424
virtual tmp< scalarField > sample(const volScalarField &) const
Sample field on surface.
virtual ~sampledCuttingPlane()
Destructor.
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
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
virtual bool expire()
Mark the surface as needing an update.
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
An abstract class for surfaces with sampling.
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
label findPatchID(const word &patchName) const
Find patch index given a name.
Macros for easy insertion into run-time selection tables.
static const pointMesh & New(const polyMesh &mesh)
Definition: MeshObject.C:44
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:52
HashSet< label, Hash< label > > labelHashSet
A HashSet with label keys.
Definition: HashSet.H:210
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
dynamicFvMesh & mesh
const pointField & points
fvPatchField< scalar > fvPatchScalarField
virtual bool needsUpdate() const
Does the surface need an update?
A class for handling words, derived from string.
Definition: word.H:59
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
const cellZoneMesh & cellZones() const
Return cell zone mesh.
Definition: polyMesh.H:472
wordList names() const
Return a list of patch names.
static const word null
An empty word.
Definition: word.H:77
#define DebugInfo
Report an information message using Foam::Info.
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.
defineTypeNameAndDebug(combustionModel, 0)
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
virtual void print(Ostream &) const
Write.
addNamedToRunTimeSelectionTable(GAMGProcAgglomeration, noneGAMGProcAgglomeration, GAMGAgglomeration, none)
label size() const
Return the number of elements in the UPtrList.
Definition: UPtrListI.H:29
label patchi
T lookupOrDefault(const word &, const T &, bool recursive=false, bool patternMatch=true) const
Find and return a T,.
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
sampledCuttingPlane(const word &name, const polyMesh &mesh, const dictionary &dict)
Construct from dictionary.
messageStream Info
SubField< vector > subField
Declare type of subField.
Definition: Field.H:89
virtual Ostream & write(const token &)=0
Write next token to stream.
const doubleScalar e
Elementary charge.
Definition: doubleScalar.H:98
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
A class for managing temporary objects.
Definition: PtrList.H:53
virtual bool update()
Update the surface as required.
Namespace for OpenFOAM.
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:576