probes.C 13.4 KB
Newer Older
1
2
3
4
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
5
    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
6
     \\/     M anipulation  | Copyright (C) 2015-2017 OpenCFD Ltd.
7
8
9
10
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

11
12
13
14
    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.
15
16
17
18
19
20
21

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    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
27
28
29
30

\*---------------------------------------------------------------------------*/

#include "probes.H"
#include "volFields.H"
#include "dictionary.H"
#include "Time.H"
#include "IOmanip.H"
31
#include "mapPolyMesh.H"
32
#include "addToRunTimeSelectionTable.H"
33
34
35

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

36
37
namespace Foam
{
38
    defineTypeNameAndDebug(probes, 0);
39
40
41
42
43
44
45

    addToRunTimeSelectionTable
    (
        functionObject,
        probes,
        dictionary
    );
46
}
47

48

Henry Weller's avatar
Henry Weller committed
49
// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
50

51
void Foam::probes::findElements(const fvMesh& mesh)
52
{
53
    DebugInfo<< "probes: resetting sample locations" << endl;
54

55
56
    elementList_.clear();
    elementList_.setSize(size());
57

58
59
60
    faceList_.clear();
    faceList_.setSize(size());

61
62
63
    processor_.setSize(size());
    processor_ = -1;

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
75
            const labelList& cellFaces = mesh.cells()[celli];
            const vector& cellCentre = mesh.cellCentres()[celli];
76
77
            scalar minDistance = GREAT;
            label minFaceID = -1;
78
            forAll(cellFaces, i)
79
            {
80
81
                label facei = cellFaces[i];
                vector dist = mesh.faceCentres()[facei] - cellCentre;
82
83
84
                if (mag(dist) < minDistance)
                {
                    minDistance = mag(dist);
85
                    minFaceID = facei;
86
87
                }
            }
88
            faceList_[probei] = minFaceID;
89
90
91
        }
        else
        {
92
            faceList_[probei] = -1;
93
94
        }

95
        if (debug && (elementList_[probei] != -1 || faceList_[probei] != -1))
96
97
        {
            Pout<< "probes : found point " << location
98
99
                << " in cell " << elementList_[probei]
                << " and face " << faceList_[probei] << endl;
100
        }
101
    }
102
103


104
    // Check if all probes have been found.
105
    forAll(elementList_, probei)
106
    {
107
108
109
        const vector& location = operator[](probei);
        label celli = elementList_[probei];
        label facei = faceList_[probei];
110

111
112
        processor_[probei] = (celli != -1 ? Pstream::myProcNo() : -1);

113
        // Check at least one processor with cell.
114
115
        reduce(celli, maxOp<label>());
        reduce(facei, maxOp<label>());
116
        reduce(processor_[probei], maxOp<label>());
117

118
        if (celli == -1)
119
120
        {
            if (Pstream::master())
121
            {
122
                WarningInFunction
123
124
                    << "Did not find location " << location
                    << " in any cell. Skipping location." << endl;
125
            }
126
        }
127
        else if (facei == -1)
128
129
130
        {
            if (Pstream::master())
            {
131
                WarningInFunction
132
133
134
135
                    << "Did not find location " << location
                    << " in any face. Skipping location." << endl;
            }
        }
136
137
138
        else
        {
            // Make sure location not on two domains.
139
            if (elementList_[probei] != -1 && elementList_[probei] != celli)
140
            {
141
                WarningInFunction
142
143
                    << "Location " << location
                    << " seems to be on multiple domains:"
144
                    << " cell " << elementList_[probei]
145
                    << " on my domain " << Pstream::myProcNo()
146
147
                    << " and cell " << celli << " on some other domain."
                    << nl
148
149
150
                    << "This might happen if the probe location is on"
                    << " a processor patch. Change the location slightly"
                    << " to prevent this." << endl;
151
            }
152

153
            if (faceList_[probei] != -1 && faceList_[probei] != facei)
154
            {
155
                WarningInFunction
156
157
                    << "Location " << location
                    << " seems to be on multiple domains:"
158
                    << " cell " << faceList_[probei]
159
                    << " on my domain " << Pstream::myProcNo()
160
161
                    << " and face " << facei << " on some other domain."
                    << nl
162
163
164
165
                    << "This might happen if the probe location is on"
                    << " a processor patch. Change the location slightly"
                    << " to prevent this." << endl;
            }
166
167
168
169
170
        }
    }
}


171
Foam::label Foam::probes::prepare()
172
{
173
    const label nFields = classifyFields();
174

175
176
    // adjust file streams
    if (Pstream::master())
177
    {
178
179
        wordHashSet currentFields;

180
181
182
183
184
        currentFields.insert(scalarFields_);
        currentFields.insert(vectorFields_);
        currentFields.insert(sphericalTensorFields_);
        currentFields.insert(symmTensorFields_);
        currentFields.insert(tensorFields_);
185

186
187
188
189
190
191
        currentFields.insert(surfaceScalarFields_);
        currentFields.insert(surfaceVectorFields_);
        currentFields.insert(surfaceSphericalTensorFields_);
        currentFields.insert(surfaceSymmTensorFields_);
        currentFields.insert(surfaceTensorFields_);

192
193
194
195
        DebugInfo
            << "Probing fields: " << currentFields << nl
            << "Probing locations: " << *this << nl
            << endl;
196
197
198


        fileName probeDir;
199
        fileName probeSubDir = name();
mattijs's avatar
mattijs committed
200

201
        if (mesh_.name() != polyMesh::defaultRegion)
mattijs's avatar
mattijs committed
202
        {
203
            probeSubDir = probeSubDir/mesh_.name();
mattijs's avatar
mattijs committed
204
        }
205
        probeSubDir = "postProcessing"/probeSubDir/mesh_.time().timeName();
mattijs's avatar
mattijs committed
206

207
208
209
210
        if (Pstream::parRun())
        {
            // Put in undecomposed case
            // (Note: gives problems for distributed data running)
211
            probeDir = mesh_.time().path()/".."/probeSubDir;
212
213
214
        }
        else
        {
215
            probeDir = mesh_.time().path()/probeSubDir;
216
217
        }

218
        // ignore known fields, close streams for fields that no longer exist
219
220
        forAllIter(HashPtrTable<OFstream>, probeFilePtrs_, iter)
        {
221
            if (!currentFields.erase(iter.key()))
222
            {
223
                DebugInfo<< "close probe stream: " << iter()->name() << endl;
224
225
226
227
228

                delete probeFilePtrs_.remove(iter);
            }
        }

229
230
        // currentFields now just has the new fields - open streams for them
        forAllConstIter(wordHashSet, currentFields, iter)
231
        {
232
            const word& fieldName = iter.key();
233

234
235
            // Create directory if does not exist.
            mkDir(probeDir);
236

237
238
239
            OFstream* fPtr = new OFstream(probeDir/fieldName);

            OFstream& fout = *fPtr;
240

241
            DebugInfo<< "open probe stream: " << fout.name() << endl;
242

243
            probeFilePtrs_.insert(fieldName, fPtr);
244

245
            unsigned int w = IOstream::defaultPrecision() + 7;
246

247
            forAll(*this, probei)
248
            {
249
250
251
252
253
254
255
                fout<< "# Probe " << probei << ' ' << operator[](probei);

                if (processor_[probei] == -1)
                {
                    fout<< "  # Not Found";
                }
                fout<< endl;
256
            }
257

258
259
260
            fout<< '#' << setw(IOstream::defaultPrecision() + 6)
                << "Probe";

261
            forAll(*this, probei)
262
            {
263
264
265
266
                if (includeOutOfBounds_ || processor_[probei] != -1)
                {
                    fout<< ' ' << setw(w) << probei;
                }
267
            }
268
            fout<< endl;
269

270
            fout<< '#' << setw(IOstream::defaultPrecision() + 6)
271
                << "Time" << endl;
272
273
274
        }
    }

275
    return nFields;
276
277
278
279
}


// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
280

281
282
283
Foam::probes::probes
(
    const word& name,
284
285
286
287
    const Time& runTime,
    const dictionary& dict,
    const bool loadFromFiles,
    const bool readFields
288
289
)
:
290
    stateFunctionObject(name, runTime),
291
    pointField(0),
292
293
294
295
296
297
298
299
300
301
302
    mesh_
    (
        refCast<const fvMesh>
        (
            runTime.lookupObject<objectRegistry>
            (
                dict.lookupOrDefault("region", polyMesh::defaultRegion)
            )
        )
    ),
    loadFromFiles_(loadFromFiles),
303
304
    fieldSelection_(),
    fixedLocations_(true),
305
306
    interpolationScheme_("cell"),
    includeOutOfBounds_(true)
307
{
308
309
310
311
    if (readFields)
    {
        read(dict);
    }
312
313
}

314
315
316
317
318
319
320
321
// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

Foam::probes::~probes()
{}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

322
bool Foam::probes::read(const dictionary& dict)
Andrew Heather's avatar
Andrew Heather committed
323
{
324
325
    dict.lookup("probeLocations") >> *this;
    dict.lookup("fields") >> fieldSelection_;
Andrew Heather's avatar
Andrew Heather committed
326

327
328
329
330
331
332
333
334
    dict.readIfPresent("fixedLocations", fixedLocations_);
    if (dict.readIfPresent("interpolationScheme", interpolationScheme_))
    {
        if (!fixedLocations_ && interpolationScheme_ != "cell")
        {
            WarningInFunction
                << "Only cell interpolation can be applied when "
                << "not using fixedLocations.  InterpolationScheme "
335
336
                << "entry will be ignored"
                << endl;
337
338
        }
    }
339
    dict.readIfPresent("includeOutOfBounds", includeOutOfBounds_);
Andrew Heather's avatar
Andrew Heather committed
340

341
342
    // Initialise cells to sample from supplied locations
    findElements(mesh_);
Andrew Heather's avatar
Andrew Heather committed
343

344
    prepare();
345

346
    return true;
347
348
349
}


350
bool Foam::probes::execute()
351
{
352
    return true;
353
354
355
}


356
bool Foam::probes::write()
357
{
358
    if (size() && prepare())
359
360
361
362
363
364
    {
        sampleAndWrite(scalarFields_);
        sampleAndWrite(vectorFields_);
        sampleAndWrite(sphericalTensorFields_);
        sampleAndWrite(symmTensorFields_);
        sampleAndWrite(tensorFields_);
365
366
367
368
369
370

        sampleAndWriteSurfaceFields(surfaceScalarFields_);
        sampleAndWriteSurfaceFields(surfaceVectorFields_);
        sampleAndWriteSurfaceFields(surfaceSphericalTensorFields_);
        sampleAndWriteSurfaceFields(surfaceSymmTensorFields_);
        sampleAndWriteSurfaceFields(surfaceTensorFields_);
371
    }
372

373
    return true;
374
375
376
}


377
378
void Foam::probes::updateMesh(const mapPolyMesh& mpm)
{
379
380
381
    DebugInfo<< "probes: updateMesh" << endl;

    if (&mpm.mesh() != &mesh_)
382
    {
383
        return;
384
385
386
387
388
389
390
391
    }

    if (fixedLocations_)
    {
        findElements(mesh_);
    }
    else
    {
392
        DebugInfo<< "probes: remapping sample locations" << endl;
393
394
395
396
397
398
399
400

        // 1. Update cells
        {
            DynamicList<label> elems(elementList_.size());

            const labelList& reverseMap = mpm.reverseCellMap();
            forAll(elementList_, i)
            {
401
                label celli = elementList_[i];
402
                if (celli != -1)
403
                {
404
                    label newCelli = reverseMap[celli];
405
                    if (newCelli == -1)
406
407
408
                    {
                        // cell removed
                    }
409
                    else if (newCelli < -1)
410
411
                    {
                        // cell merged
412
                        elems.append(-newCelli - 2);
413
414
415
416
                    }
                    else
                    {
                        // valid new cell
417
                        elems.append(newCelli);
418
                    }
419
420
421
                }
                else
                {
422
423
                    // Keep -1 elements so the size stays the same
                    elems.append(-1);
424
425
426
427
428
429
430
431
432
433
434
435
436
                }
            }

            elementList_.transfer(elems);
        }

        // 2. Update faces
        {
            DynamicList<label> elems(faceList_.size());

            const labelList& reverseMap = mpm.reverseFaceMap();
            forAll(faceList_, i)
            {
437
                label facei = faceList_[i];
438
                if (facei != -1)
439
                {
440
441
                    label newFacei = reverseMap[facei];
                    if (newFacei == -1)
442
443
444
                    {
                        // face removed
                    }
445
                    else if (newFacei < -1)
446
447
                    {
                        // face merged
448
                        elems.append(-newFacei - 2);
449
450
451
452
                    }
                    else
                    {
                        // valid new face
453
                        elems.append(newFacei);
454
                    }
455
456
457
                }
                else
                {
458
459
                    // Keep -1 elements
                    elems.append(-1);
460
461
462
463
464
465
466
467
468
                }
            }

            faceList_.transfer(elems);
        }
    }
}


469
void Foam::probes::movePoints(const polyMesh& mesh)
470
{
471
    DebugInfo<< "probes: movePoints" << endl;
472

473
    if (fixedLocations_ && &mesh == &mesh_)
474
475
476
477
478
479
    {
        findElements(mesh_);
    }
}


480
// ************************************************************************* //