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