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