probes.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 "probes.H"
27 #include "volFields.H"
28 #include "dictionary.H"
29 #include "Time.H"
30 #include "IOmanip.H"
31 #include "mapPolyMesh.H"
33 
34 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
35 
36 namespace Foam
37 {
38  defineTypeNameAndDebug(probes, 0);
39 
41  (
42  functionObject,
43  probes,
44  dictionary
45  );
46 }
47 
48 
49 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
50 
52 {
53  if (debug)
54  {
55  Info<< "probes: resetting sample locations" << endl;
56  }
57 
58  elementList_.clear();
59  elementList_.setSize(size());
60 
61  faceList_.clear();
62  faceList_.setSize(size());
63 
64  forAll(*this, probei)
65  {
66  const vector& location = operator[](probei);
67 
68  const label celli = mesh.findCell(location);
69 
70  elementList_[probei] = celli;
71 
72  if (celli != -1)
73  {
74  const labelList& cellFaces = mesh.cells()[celli];
75  const vector& cellCentre = mesh.cellCentres()[celli];
76  scalar minDistance = great;
77  label minFaceID = -1;
78  forAll(cellFaces, i)
79  {
80  label facei = cellFaces[i];
81  vector dist = mesh.faceCentres()[facei] - cellCentre;
82  if (mag(dist) < minDistance)
83  {
84  minDistance = mag(dist);
85  minFaceID = facei;
86  }
87  }
88  faceList_[probei] = minFaceID;
89  }
90  else
91  {
92  faceList_[probei] = -1;
93  }
94 
95  if (debug && (elementList_[probei] != -1 || faceList_[probei] != -1))
96  {
97  Pout<< "probes : found point " << location
98  << " in cell " << elementList_[probei]
99  << " and face " << faceList_[probei] << endl;
100  }
101  }
102 
103 
104  // Check if all probes have been found.
105  forAll(elementList_, probei)
106  {
107  const vector& location = operator[](probei);
108  label celli = elementList_[probei];
109  label facei = faceList_[probei];
110 
111  // Check at least one processor with cell.
112  reduce(celli, maxOp<label>());
113  reduce(facei, maxOp<label>());
114 
115  if (celli == -1)
116  {
117  if (Pstream::master())
118  {
120  << "Did not find location " << location
121  << " in any cell. Skipping location." << endl;
122  }
123  }
124  else if (facei == -1)
125  {
126  if (Pstream::master())
127  {
129  << "Did not find location " << location
130  << " in any face. Skipping location." << endl;
131  }
132  }
133  else
134  {
135  // Make sure location not on two domains.
136  if (elementList_[probei] != -1 && elementList_[probei] != celli)
137  {
139  << "Location " << location
140  << " seems to be on multiple domains:"
141  << " cell " << elementList_[probei]
142  << " on my domain " << Pstream::myProcNo()
143  << " and cell " << celli << " on some other domain."
144  << endl
145  << "This might happen if the probe location is on"
146  << " a processor patch. Change the location slightly"
147  << " to prevent this." << endl;
148  }
149 
150  if (faceList_[probei] != -1 && faceList_[probei] != facei)
151  {
153  << "Location " << location
154  << " seems to be on multiple domains:"
155  << " cell " << faceList_[probei]
156  << " on my domain " << Pstream::myProcNo()
157  << " and face " << facei << " on some other domain."
158  << endl
159  << "This might happen if the probe location is on"
160  << " a processor patch. Change the location slightly"
161  << " to prevent this." << endl;
162  }
163  }
164  }
165 }
166 
167 
169 {
170  const label nFields = classifyFields();
171 
172  // adjust file streams
173  if (Pstream::master())
174  {
175  wordHashSet currentFields;
176 
177  currentFields.insert(scalarFields_);
178  currentFields.insert(vectorFields_);
179  currentFields.insert(sphericalTensorFields_);
180  currentFields.insert(symmTensorFields_);
181  currentFields.insert(tensorFields_);
182 
183  currentFields.insert(surfaceScalarFields_);
184  currentFields.insert(surfaceVectorFields_);
185  currentFields.insert(surfaceSphericalTensorFields_);
186  currentFields.insert(surfaceSymmTensorFields_);
187  currentFields.insert(surfaceTensorFields_);
188 
189  if (debug)
190  {
191  Info<< "Probing fields: " << currentFields << nl
192  << "Probing locations: " << *this << nl
193  << endl;
194  }
195 
196 
197  fileName probeDir;
198  fileName probeSubDir = name();
199 
200  if (mesh_.name() != polyMesh::defaultRegion)
201  {
202  probeSubDir = probeSubDir/mesh_.name();
203  }
204  probeSubDir = "postProcessing"/probeSubDir/mesh_.time().timeName();
205 
206  if (Pstream::parRun())
207  {
208  // Put in undecomposed case
209  // (Note: gives problems for distributed data running)
210  probeDir = mesh_.time().path()/".."/probeSubDir;
211  }
212  else
213  {
214  probeDir = mesh_.time().path()/probeSubDir;
215  }
216  // Remove ".."
217  probeDir.clean();
218 
219  // ignore known fields, close streams for fields that no longer exist
220  forAllIter(HashPtrTable<OFstream>, probeFilePtrs_, iter)
221  {
222  if (!currentFields.erase(iter.key()))
223  {
224  if (debug)
225  {
226  Info<< "close probe stream: " << iter()->name() << endl;
227  }
228 
229  delete probeFilePtrs_.remove(iter);
230  }
231  }
232 
233  // currentFields now just has the new fields - open streams for them
234  forAllConstIter(wordHashSet, currentFields, iter)
235  {
236  const word& fieldName = iter.key();
237 
238  // Create directory if does not exist.
239  mkDir(probeDir);
240 
241  OFstream* fPtr = new OFstream(probeDir/fieldName);
242  OFstream& os = *fPtr;
243 
244  if (debug)
245  {
246  Info<< "open probe stream: " << os.name() << endl;
247  }
248 
249  probeFilePtrs_.insert(fieldName, fPtr);
250 
251  const unsigned int w = IOstream::defaultPrecision() + 7;
252  os << setf(ios_base::left);
253 
254  forAll(*this, probei)
255  {
256  os<< "# Probe " << probei << ' ' << operator[](probei)
257  << endl;
258  }
259 
260  os << setw(w) << "# Time";
261 
262  forAll(*this, probei)
263  {
264  os<< ' ' << setw(w) << probei;
265  }
266  os<< endl;
267  }
268  }
269 
270  return nFields;
271 }
272 
273 
274 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
275 
277 (
278  const word& name,
279  const Time& t,
280  const dictionary& dict
281 )
282 :
283  functionObject(name),
284  pointField(0),
285  mesh_
286  (
287  refCast<const fvMesh>
288  (
290  (
292  )
293  )
294  ),
295  loadFromFiles_(false),
296  fieldSelection_(),
297  fixedLocations_(true),
298  interpolationScheme_("cell")
299 {
300  read(dict);
301 }
302 
303 
305 (
306  const word& name,
307  const objectRegistry& obr,
308  const dictionary& dict,
309  const bool loadFromFiles
310 )
311 :
312  functionObject(name),
313  pointField(0),
314  mesh_(refCast<const fvMesh>(obr)),
315  loadFromFiles_(loadFromFiles),
316  fieldSelection_(),
317  fixedLocations_(true),
318  interpolationScheme_("cell")
319 {
320  read(dict);
321 }
322 
323 
324 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
325 
327 {}
328 
329 
330 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
331 
333 {
334  dict.lookup("probeLocations") >> *this;
335  dict.lookup("fields") >> fieldSelection_;
336 
337  dict.readIfPresent("fixedLocations", fixedLocations_);
338  if (dict.readIfPresent("interpolationScheme", interpolationScheme_))
339  {
340  if (!fixedLocations_ && interpolationScheme_ != "cell")
341  {
343  << "Only cell interpolation can be applied when "
344  << "not using fixedLocations. InterpolationScheme "
345  << "entry will be ignored";
346  }
347  }
348 
349  // Initialise cells to sample from supplied locations
350  findElements(mesh_);
351 
352  prepare();
353 
354  return true;
355 }
356 
357 
359 {
360  return true;
361 }
362 
363 
365 {
366  if (size() && prepare())
367  {
368  sampleAndWrite(scalarFields_);
369  sampleAndWrite(vectorFields_);
370  sampleAndWrite(sphericalTensorFields_);
371  sampleAndWrite(symmTensorFields_);
372  sampleAndWrite(tensorFields_);
373 
374  sampleAndWriteSurfaceFields(surfaceScalarFields_);
375  sampleAndWriteSurfaceFields(surfaceVectorFields_);
376  sampleAndWriteSurfaceFields(surfaceSphericalTensorFields_);
377  sampleAndWriteSurfaceFields(surfaceSymmTensorFields_);
378  sampleAndWriteSurfaceFields(surfaceTensorFields_);
379  }
380 
381  return true;
382 }
383 
384 
386 {
387  DebugInfo<< "probes: updateMesh" << endl;
388 
389  if (&mpm.mesh() != &mesh_)
390  {
391  return;
392  }
393 
394  if (fixedLocations_)
395  {
396  findElements(mesh_);
397  }
398  else
399  {
400  if (debug)
401  {
402  Info<< "probes: remapping sample locations" << endl;
403  }
404 
405  // 1. Update cells
406  {
407  DynamicList<label> elems(elementList_.size());
408 
409  const labelList& reverseMap = mpm.reverseCellMap();
410  forAll(elementList_, i)
411  {
412  label celli = elementList_[i];
413  label newCelli = reverseMap[celli];
414  if (newCelli == -1)
415  {
416  // cell removed
417  }
418  else if (newCelli < -1)
419  {
420  // cell merged
421  elems.append(-newCelli - 2);
422  }
423  else
424  {
425  // valid new cell
426  elems.append(newCelli);
427  }
428  }
429 
430  elementList_.transfer(elems);
431  }
432 
433  // 2. Update faces
434  {
435  DynamicList<label> elems(faceList_.size());
436 
437  const labelList& reverseMap = mpm.reverseFaceMap();
438  forAll(faceList_, i)
439  {
440  label facei = faceList_[i];
441  label newFacei = reverseMap[facei];
442  if (newFacei == -1)
443  {
444  // face removed
445  }
446  else if (newFacei < -1)
447  {
448  // face merged
449  elems.append(-newFacei - 2);
450  }
451  else
452  {
453  // valid new face
454  elems.append(newFacei);
455  }
456  }
457 
458  faceList_.transfer(elems);
459  }
460  }
461 }
462 
463 
465 {
466  DebugInfo<< "probes: movePoints" << endl;
467 
468  if (fixedLocations_ && &mesh == &mesh_)
469  {
470  findElements(mesh_);
471  }
472 }
473 
474 
475 // ************************************************************************* //
A HashTable with keys but without contents.
Definition: HashSet.H:59
#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
A class for handling file names.
Definition: fileName.H:79
virtual bool write()
Sample and write.
Definition: probes.C:364
label prepare()
Classify field type and Open/close file streams,.
Definition: probes.C:168
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:158
virtual void updateMesh(const mapPolyMesh &)
Update for changes of mesh.
Definition: probes.C:385
virtual bool execute()
Execute, currently does nothing.
Definition: probes.C:358
bool clean()
Cleanup file name.
Definition: fileName.C:81
#define forAllIter(Container, container, iter)
Iterate across all elements in the container object of type.
Definition: UList.H:459
Output to file stream.
Definition: OFstream.H:82
static unsigned int defaultPrecision()
Return the default precision.
Definition: IOstream.H:461
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:429
static word defaultRegion
Return the default region name.
Definition: polyMesh.H:309
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
static bool master(const label communicator=0)
Am I the master process.
Definition: UPstream.H:423
const cellList & cells() const
A HashTable specialization for hashing pointers.
Definition: HashPtrTable.H:50
Abstract base-class for Time/database functionObjects.
const fileName & name() const
Return the name of the stream.
Definition: OFstream.H:120
bool insert(const Key &key)
Insert a new entry.
Definition: HashSet.H:111
const Type & lookupObject(const word &name) const
Lookup and return the object of the given Type.
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
Definition: Time.H:68
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
Definition: mapPolyMesh.H:158
Macros for easy insertion into run-time selection tables.
bool erase(const iterator &)
Erase a hashedEntry specified by given iterator.
Definition: HashTable.C:371
bool read(const char *, int32_t &)
Definition: int32IO.C:85
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
dynamicFvMesh & mesh
A class for handling words, derived from string.
Definition: word.H:59
virtual void findElements(const fvMesh &)
Find cells and faces containing probes.
Definition: probes.C:51
word name() const
Return file name (part beyond last /)
Definition: fileName.C:183
bool readIfPresent(const word &, T &, bool recursive=false, bool patternMatch=true) const
Find an entry if present, and assign to T.
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:296
Smanip< ios_base::fmtflags > setf(const ios_base::fmtflags flags)
Definition: IOmanip.H:164
const vectorField & cellCentres() const
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:29
#define DebugInfo
Report an information message using Foam::Info.
Istream and Ostream manipulators taking arguments.
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
static const char nl
Definition: Ostream.H:260
const labelList & reverseCellMap() const
Reverse cell map.
Definition: mapPolyMesh.H:521
defineTypeNameAndDebug(combustionModel, 0)
bool mkDir(const fileName &, mode_t=0777)
Make a directory and return an error if it could not be created.
Definition: POSIX.C:290
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
virtual void movePoints(const polyMesh &)
Update for changes of mesh.
Definition: probes.C:464
const vectorField & faceCentres() const
static bool & parRun()
Is this a parallel run?
Definition: UPstream.H:399
T lookupOrDefault(const word &, const T &, bool recursive=false, bool patternMatch=true) const
Find and return a T,.
#define WarningInFunction
Report a warning using Foam::Warning.
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
const polyMesh & mesh() const
Return polyMesh.
Definition: mapPolyMesh.H:352
virtual bool read(const dictionary &)
Read the probes.
Definition: probes.C:332
label findCell(const point &p, const cellDecomposition=CELL_TETS) const
Find cell enclosing this location and return index.
Definition: polyMesh.C:1606
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
fileName path() const
Return directory path name (part before last /)
Definition: fileName.C:253
probes(const word &name, const Time &time, const dictionary &dict)
Construct from Time and dictionary.
Definition: probes.C:277
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
Omanip< int > setw(const int i)
Definition: IOmanip.H:199
Registry of regIOobjects.
virtual ~probes()
Destructor.
Definition: probes.C:326
const labelList & reverseFaceMap() const
Reverse face map.
Definition: mapPolyMesh.H:490
Namespace for OpenFOAM.
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:812