fluxSummary.C 30.7 KB
Newer Older
1
2
3
4
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
OpenFOAM bot's avatar
OpenFOAM bot committed
5
    \\  /    A nd           | www.openfoam.com
OpenFOAM bot's avatar
OpenFOAM bot committed
6
7
     \\/     M anipulation  |
-------------------------------------------------------------------------------
OpenFOAM bot's avatar
OpenFOAM bot committed
8
    Copyright (C) 2015 OpenFOAM Foundation
9
    Copyright (C) 2015-2020 OpenCFD Ltd.
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    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.

    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
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

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

#include "fluxSummary.H"
#include "surfaceFields.H"
31
#include "polySurfaceFields.H"
32
33
34
35
36
#include "dictionary.H"
#include "Time.H"
#include "syncTools.H"
#include "meshTools.H"
#include "PatchEdgeFaceWave.H"
37
#include "edgeTopoDistanceData.H"
38
39
#include "globalIndex.H"
#include "OBJstream.H"
40
#include "addToRunTimeSelectionTable.H"
41
42
43
44

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

namespace Foam
45
46
{
namespace functionObjects
47
48
{
    defineTypeNameAndDebug(fluxSummary, 0);
49
    addToRunTimeSelectionTable(functionObject, fluxSummary, dictionary);
50
}
51
52
53
}

const Foam::Enum
54
<
55
56
57
    Foam::functionObjects::fluxSummary::modeType
>
Foam::functionObjects::fluxSummary::modeTypeNames_
Mark Olesen's avatar
Mark Olesen committed
58
({
59
60
61
    { modeType::mdFaceZone , "faceZone" },
    { modeType::mdFaceZoneAndDirection, "faceZoneAndDirection" },
    { modeType::mdCellZoneAndDirection, "cellZoneAndDirection" },
62
    { modeType::mdSurface, "functionObjectSurface" },
63
64
    { modeType::mdSurface, "surface" },
    { modeType::mdSurfaceAndDirection, "surfaceAndDirection" },
Mark Olesen's avatar
Mark Olesen committed
65
});
66
67


68
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * * //
69

70
71
bool Foam::functionObjects::fluxSummary::isSurfaceMode() const
{
72
    return (mdSurface == mode_ || mdSurfaceAndDirection == mode_);
73
74
75
}


76
77
Foam::word Foam::functionObjects::fluxSummary::checkFlowType
(
78
    const dimensionSet& fieldDims,
79
80
81
    const word& fieldName
) const
{
82
    // Surfaces are multiplied by their area, so account for that
83
    // in the dimension checking
84
85
    const dimensionSet dims =
        (fieldDims * (isSurfaceMode() ? dimTime*dimArea : dimTime));
86
87

    if (dims == dimVolume)
88
89
90
    {
        return "volumetric";
    }
91
    else if (dims == dimMass)
92
93
94
95
    {
        return "mass";
    }

96
97
98
99
100
101
102
    FatalErrorInFunction
        << "Unsupported flux field " << fieldName << " with dimensions "
        << fieldDims
        << ".  Expected either mass flow or volumetric flow rate."
        << abort(FatalError);

    return word::null;
103
104
105
}


106
107
108
109
110
111
112
113
void Foam::functionObjects::fluxSummary::initialiseSurface
(
    const word& surfName,
    DynamicList<word>& names,
    DynamicList<vector>& directions,
    DynamicList<boolList>& faceFlip
) const
{
114
115
    const polySurface* surfptr =
        storedObjects().cfindObject<polySurface>(surfName);
116
117

    if (!surfptr)
118
119
120
    {
        FatalErrorInFunction
            << "Unable to find surface " << surfName
121
            << ".  Valid surfaces: "
122
            << storedObjects().sortedNames<polySurface>() << nl
123
124
125
126
            << exit(FatalError);
    }

    names.append(surfName);
127
128
    directions.append(Zero);        // Dummy value
    faceFlip.append(boolList());    // No flip-map
129
130
131
132
133
134
135
136
137
138
139
140
}


void Foam::functionObjects::fluxSummary::initialiseSurfaceAndDirection
(
    const word& surfName,
    const vector& dir,
    DynamicList<word>& names,
    DynamicList<vector>& directions,
    DynamicList<boolList>& faceFlip
) const
{
141
142
    const polySurface* surfptr =
        storedObjects().cfindObject<polySurface>(surfName);
143
144

    if (!surfptr)
145
146
147
    {
        FatalErrorInFunction
            << "Unable to find surface " << surfName
148
            << ".  Valid surfaces: "
149
            << storedObjects().sortedNames<polySurface>() << nl
150
151
152
            << exit(FatalError);
    }

153
    const auto& s = *surfptr;
154
155
156
157
    const vector refDir = dir/(mag(dir) + ROOTVSMALL);

    names.append(surfName);
    directions.append(refDir);
158
    faceFlip.append(boolList());    // No flip-map
159
160
161
162
163
164

    boolList& flips = faceFlip[faceFlip.size()-1];
    flips.setSize(s.size(), false);

    forAll(s, i)
    {
165
        // Orientation set by comparison with reference direction
166
167
168
169
170
171
172
173
174
175
176
177
178
179
        const vector& n = s.faceNormals()[i];

        if ((n & refDir) > tolerance_)
        {
            flips[i] = false;
        }
        else
        {
            flips[i] = true;
        }
    }
}


180
void Foam::functionObjects::fluxSummary::initialiseFaceZone
181
182
(
    const word& faceZoneName,
183
184
    DynamicList<word>& names,
    DynamicList<vector>& directions,
185
186
    DynamicList<labelList>& faceID,
    DynamicList<labelList>& facePatchID,
187
    DynamicList<boolList>& faceFlip
188
189
) const
{
190
    label zonei = mesh_.faceZones().findZoneID(faceZoneName);
Andrew Heather's avatar
Andrew Heather committed
191
    if (zonei == -1)
192
    {
193
        FatalErrorInFunction
194
            << "Unable to find faceZone " << faceZoneName
195
196
            << ".  Valid zones: "
            << mesh_.faceZones().sortedNames() << nl
197
198
            << exit(FatalError);
    }
199
    const faceZone& fZone = mesh_.faceZones()[zonei];
200

201
202
203
    names.append(faceZoneName);
    directions.append(Zero); // dummy value

204
205
    DynamicList<label> faceIDs(fZone.size());
    DynamicList<label> facePatchIDs(fZone.size());
206
    DynamicList<bool>  flips(fZone.size());
207
208
209

    forAll(fZone, i)
    {
Andrew Heather's avatar
Andrew Heather committed
210
        label facei = fZone[i];
211
212
213

        label faceID = -1;
        label facePatchID = -1;
214
        if (mesh_.isInternalFace(facei))
215
        {
Andrew Heather's avatar
Andrew Heather committed
216
            faceID = facei;
217
218
219
220
            facePatchID = -1;
        }
        else
        {
221
222
            facePatchID = mesh_.boundaryMesh().whichPatch(facei);
            const polyPatch& pp = mesh_.boundaryMesh()[facePatchID];
223
224
225
226
            if (isA<coupledPolyPatch>(pp))
            {
                if (refCast<const coupledPolyPatch>(pp).owner())
                {
Andrew Heather's avatar
Andrew Heather committed
227
                    faceID = pp.whichFace(facei);
228
229
230
231
232
233
234
235
                }
                else
                {
                    faceID = -1;
                }
            }
            else if (!isA<emptyPolyPatch>(pp))
            {
Andrew Heather's avatar
Andrew Heather committed
236
                faceID = facei - pp.start();
237
238
239
240
241
242
243
244
245
246
            }
            else
            {
                faceID = -1;
                facePatchID = -1;
            }
        }

        if (faceID >= 0)
        {
247
            // Orientation set by faceZone flip map
248
            if (fZone.flipMap()[i])
249
            {
250
                flips.append(true);
251
252
253
            }
            else
            {
254
                flips.append(false);
255
256
257
258
259
260
261
            }

            faceIDs.append(faceID);
            facePatchIDs.append(facePatchID);
        }
    }

262
    // could reduce some copying here
263
264
    faceID.append(faceIDs);
    facePatchID.append(facePatchIDs);
265
    faceFlip.append(flips);
266
267
268
}


269
void Foam::functionObjects::fluxSummary::initialiseFaceZoneAndDirection
270
271
272
(
    const word& faceZoneName,
    const vector& dir,
273
274
    DynamicList<word>& names,
    DynamicList<vector>& directions,
275
276
    DynamicList<labelList>& faceID,
    DynamicList<labelList>& facePatchID,
277
    DynamicList<boolList>& faceFlip
278
279
) const
{
280
    const vector refDir = dir/(mag(dir) + ROOTVSMALL);
281

282
    label zonei = mesh_.faceZones().findZoneID(faceZoneName);
Andrew Heather's avatar
Andrew Heather committed
283
    if (zonei == -1)
284
    {
285
         FatalErrorInFunction
286
            << "Unable to find faceZone " << faceZoneName
287
288
            << ".  Valid zones: "
            << mesh_.faceZones().sortedNames() << nl
289
290
            << exit(FatalError);
    }
291
    const faceZone& fZone = mesh_.faceZones()[zonei];
292

293
294
295
    names.append(faceZoneName);
    directions.append(refDir);

296
297
    DynamicList<label> faceIDs(fZone.size());
    DynamicList<label> facePatchIDs(fZone.size());
298
    DynamicList<bool>  flips(fZone.size());
299

300
301
    const surfaceVectorField& Sf = mesh_.Sf();
    const surfaceScalarField& magSf = mesh_.magSf();
302

Andrew Heather's avatar
Andrew Heather committed
303
    vector n(Zero);
304
305
306

    forAll(fZone, i)
    {
Andrew Heather's avatar
Andrew Heather committed
307
        label facei = fZone[i];
308
309
310

        label faceID = -1;
        label facePatchID = -1;
311
        if (mesh_.isInternalFace(facei))
312
        {
Andrew Heather's avatar
Andrew Heather committed
313
            faceID = facei;
314
315
316
317
            facePatchID = -1;
        }
        else
        {
318
319
            facePatchID = mesh_.boundaryMesh().whichPatch(facei);
            const polyPatch& pp = mesh_.boundaryMesh()[facePatchID];
320
321
322
323
            if (isA<coupledPolyPatch>(pp))
            {
                if (refCast<const coupledPolyPatch>(pp).owner())
                {
Andrew Heather's avatar
Andrew Heather committed
324
                    faceID = pp.whichFace(facei);
325
326
327
328
329
330
331
332
                }
                else
                {
                    faceID = -1;
                }
            }
            else if (!isA<emptyPolyPatch>(pp))
            {
Andrew Heather's avatar
Andrew Heather committed
333
                faceID = facei - pp.start();
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
            }
            else
            {
                faceID = -1;
                facePatchID = -1;
            }
        }

        if (faceID >= 0)
        {
            // orientation set by comparison with reference direction
            if (facePatchID != -1)
            {
                n = Sf.boundaryField()[facePatchID][faceID]
                   /(magSf.boundaryField()[facePatchID][faceID] + ROOTVSMALL);
            }
            else
            {
                n = Sf[faceID]/(magSf[faceID] + ROOTVSMALL);
            }

            if ((n & refDir) > tolerance_)
            {
357
                flips.append(false);
358
359
360
            }
            else
            {
361
                flips.append(true);
362
363
364
365
366
367
368
            }

            faceIDs.append(faceID);
            facePatchIDs.append(facePatchID);
        }
    }

369
    // could reduce copying here
370
371
    faceID.append(faceIDs);
    facePatchID.append(facePatchIDs);
372
    faceFlip.append(flips);
373
374
375
}


376
void Foam::functionObjects::fluxSummary::initialiseCellZoneAndDirection
377
378
379
(
    const word& cellZoneName,
    const vector& dir,
380
381
    DynamicList<word>& names,
    DynamicList<vector>& directions,
382
383
    DynamicList<labelList>& faceID,
    DynamicList<labelList>& facePatchID,
384
    DynamicList<boolList>& faceFlip
385
386
) const
{
387
    const vector refDir = dir/(mag(dir) + ROOTVSMALL);
388

389
    const label cellZonei = mesh_.cellZones().findZoneID(cellZoneName);
Andrew Heather's avatar
Andrew Heather committed
390
    if (cellZonei == -1)
391
    {
392
        FatalErrorInFunction
393
            << "Unable to find cellZone " << cellZoneName
394
395
            << ".  Valid zones: "
            << mesh_.cellZones().sortedNames() << nl
396
397
398
            << exit(FatalError);
    }

399
400
    const label nInternalFaces = mesh_.nInternalFaces();
    const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
401

402
403
    labelList cellAddr(mesh_.nCells(), -1);
    const labelList& cellIDs = mesh_.cellZones()[cellZonei];
404
    labelUIndList(cellAddr, cellIDs) = identity(cellIDs.size());
405
    labelList nbrFaceCellAddr(mesh_.nFaces() - nInternalFaces, -1);
406

Andrew Heather's avatar
Andrew Heather committed
407
    forAll(pbm, patchi)
408
    {
Andrew Heather's avatar
Andrew Heather committed
409
        const polyPatch& pp = pbm[patchi];
410
411
412
413
414

        if (pp.coupled())
        {
            forAll(pp, i)
            {
Andrew Heather's avatar
Andrew Heather committed
415
416
                label facei = pp.start() + i;
                label nbrFacei = facei - nInternalFaces;
417
                label own = mesh_.faceOwner()[facei];
Andrew Heather's avatar
Andrew Heather committed
418
                nbrFaceCellAddr[nbrFacei] = cellAddr[own];
419
420
421
422
            }
        }
    }

423
    // Correct boundary values for parallel running
424
    syncTools::swapBoundaryFaceList(mesh_, nbrFaceCellAddr);
425

426
    // Collect faces
427
    DynamicList<label> faceIDs(floor(0.1*mesh_.nFaces()));
428
429
    DynamicList<label> facePatchIDs(faceIDs.size());
    DynamicList<label> faceLocalPatchIDs(faceIDs.size());
430
    DynamicList<bool>  flips(faceIDs.size());
431

432
    // Internal faces
Andrew Heather's avatar
Andrew Heather committed
433
    for (label facei = 0; facei < nInternalFaces; facei++)
434
    {
435
436
        const label own = cellAddr[mesh_.faceOwner()[facei]];
        const label nbr = cellAddr[mesh_.faceNeighbour()[facei]];
437
438
439

        if (((own != -1) && (nbr == -1)) || ((own == -1) && (nbr != -1)))
        {
440
            vector n = mesh_.faces()[facei].unitNormal(mesh_.points());
441
442
443

            if ((n & refDir) > tolerance_)
            {
Andrew Heather's avatar
Andrew Heather committed
444
445
                faceIDs.append(facei);
                faceLocalPatchIDs.append(facei);
446
                facePatchIDs.append(-1);
447
                flips.append(false);
448
449
450
            }
            else if ((n & -refDir) > tolerance_)
            {
Andrew Heather's avatar
Andrew Heather committed
451
452
                faceIDs.append(facei);
                faceLocalPatchIDs.append(facei);
453
                facePatchIDs.append(-1);
454
                flips.append(true);
455
456
457
458
            }
        }
    }

459
    // Loop over boundary faces
Andrew Heather's avatar
Andrew Heather committed
460
    forAll(pbm, patchi)
461
    {
Andrew Heather's avatar
Andrew Heather committed
462
        const polyPatch& pp = pbm[patchi];
463

Andrew Heather's avatar
Andrew Heather committed
464
        forAll(pp, localFacei)
465
        {
Andrew Heather's avatar
Andrew Heather committed
466
            const label facei = pp.start() + localFacei;
467
            const label own = cellAddr[mesh_.faceOwner()[facei]];
Andrew Heather's avatar
Andrew Heather committed
468
            const label nbr = nbrFaceCellAddr[facei - nInternalFaces];
469
470
471

            if ((own != -1) && (nbr == -1))
            {
472
                vector n = mesh_.faces()[facei].unitNormal(mesh_.points());
473
474
475

                if ((n & refDir) > tolerance_)
                {
Andrew Heather's avatar
Andrew Heather committed
476
477
478
                    faceIDs.append(facei);
                    faceLocalPatchIDs.append(localFacei);
                    facePatchIDs.append(patchi);
479
                    flips.append(false);
480
481
482
                }
                else if ((n & -refDir) > tolerance_)
                {
Andrew Heather's avatar
Andrew Heather committed
483
484
485
                    faceIDs.append(facei);
                    faceLocalPatchIDs.append(localFacei);
                    facePatchIDs.append(patchi);
486
                    flips.append(true);
487
488
489
490
491
                }
            }
        }
    }

492
    // Convert into primitivePatch for convenience
493
494
    indirectPrimitivePatch patch
    (
495
496
        IndirectList<face>(mesh_.faces(), faceIDs),
        mesh_.points()
497
498
499
500
    );

    if (debug)
    {
501
        OBJstream os(mesh_.time().path()/"patch.obj");
502
        faceList faces(patch);
503
        os.write(faces, mesh_.points(), false);
504
505
506
    }


507
    // Data on all edges and faces
508
509
    List<edgeTopoDistanceData<label>> allEdgeInfo(patch.nEdges());
    List<edgeTopoDistanceData<label>> allFaceInfo(patch.size());
510
511
512

    bool search = true;

513
514
515
516
    DebugInfo
        << "initialiseCellZoneAndDirection: "
        << "Starting walk to split patch into faceZones"
        << endl;
517
518
519
520

    globalIndex globalFaces(patch.size());

    label oldFaceID = 0;
521
    label regioni = 0;
522
523
524
525

    // Dummy tracking data
    bool dummyData{false};

526
527
528
    while (search)
    {
        DynamicList<label> changedEdges;
529
        DynamicList<edgeTopoDistanceData<label>> changedInfo;
530

Andrew Heather's avatar
Andrew Heather committed
531
        label seedFacei = labelMax;
532
533
        for (; oldFaceID < patch.size(); oldFaceID++)
        {
534
            if (!allFaceInfo[oldFaceID].valid<bool>(dummyData))
535
            {
Andrew Heather's avatar
Andrew Heather committed
536
                seedFacei = globalFaces.toGlobal(oldFaceID);
537
538
539
                break;
            }
        }
Andrew Heather's avatar
Andrew Heather committed
540
        reduce(seedFacei, minOp<label>());
541

Andrew Heather's avatar
Andrew Heather committed
542
        if (seedFacei == labelMax)
543
544
545
546
        {
            break;
        }

Andrew Heather's avatar
Andrew Heather committed
547
        if (globalFaces.isLocal(seedFacei))
548
        {
549
            const label localFacei = globalFaces.toLocal(seedFacei);
Andrew Heather's avatar
Andrew Heather committed
550
            const labelList& fEdges = patch.faceEdges()[localFacei];
551

552
            for (const label edgei : fEdges)
553
            {
554
                if (allEdgeInfo[edgei].valid<bool>(dummyData))
555
                {
556
557
                    WarningInFunction
                        << "Problem in edge face wave: attempted to assign a "
558
                        << "value to an edge that has already been visited. "
559
                        << "Edge info: " << allEdgeInfo[edgei]
560
561
562
                        << endl;
                }

563
                changedEdges.append(edgei);
564
565
566
567
568
569
570
571
                changedInfo.append
                (
                    edgeTopoDistanceData<label>
                    (
                        0,          // distance
                        regioni
                    )
                );
572
573
574
575
576
577
578
            }
        }


        PatchEdgeFaceWave
        <
            indirectPrimitivePatch,
579
            edgeTopoDistanceData<label>
580
581
        > calc
        (
582
            mesh_,
583
584
585
586
587
588
589
590
            patch,
            changedEdges,
            changedInfo,
            allEdgeInfo,
            allFaceInfo,
            returnReduce(patch.nEdges(), sumOp<label>())
        );

591
        if (debug)
592
        {
593
            label nCells = 0;
Andrew Heather's avatar
Andrew Heather committed
594
            forAll(allFaceInfo, facei)
595
            {
596
                if (allFaceInfo[facei].data() == regioni)
597
598
599
                {
                    nCells++;
                }
600
601
            }

602
            Info<< "*** region:" << regioni
603
604
605
606
                << "  found:" << returnReduce(nCells, sumOp<label>())
                << " faces" << endl;
        }

607
        regioni++;
608
609
    }

610
    // Collect the data per region
611
612
613
614
615
616
617
618
619
620
621
    const label nRegion = regioni;

    #ifdef FULLDEBUG
    // May wish to have enabled always
    if (nRegion == 0)
    {
         FatalErrorInFunction
            << "Region split failed" << nl
            << exit(FatalError);
    }
    #endif
622

623
624
    List<DynamicList<label>> regionFaceIDs(nRegion);
    List<DynamicList<label>> regionFacePatchIDs(nRegion);
625
    List<DynamicList<bool>>  regionFaceFlips(nRegion);
626

627
    forAll(allFaceInfo, facei)
628
    {
629
        regioni = allFaceInfo[facei].data();
630

631
632
        regionFaceIDs[regioni].append(faceLocalPatchIDs[facei]);
        regionFacePatchIDs[regioni].append(facePatchIDs[facei]);
633
        regionFaceFlips[regioni].append(flips[facei]);
634
635
    }

636
637
    // Transfer to persistent storage
    forAll(regionFaceIDs, regioni)
638
    {
639
        const word zoneName = cellZoneName + ":faceZone" + Foam::name(regioni);
640
641
        names.append(zoneName);
        directions.append(refDir);
642
643
        faceID.append(regionFaceIDs[regioni]);
        facePatchID.append(regionFacePatchIDs[regioni]);
644
        faceFlip.append(regionFaceFlips[regioni]);
645

Andrew Heather's avatar
Andrew Heather committed
646
        // Write OBj of faces to file
647
648
        if (debug)
        {
649
650
651
            OBJstream os(mesh_.time().path()/zoneName + ".obj");
            faceList faces(mesh_.faces(), regionFaceIDs[regioni]);
            os.write(faces, mesh_.points(), false);
652
653
654
        }
    }

655
    if (log)
656
    {
657
        Info<< type() << " " << name() << " output:" << nl
658
659
660
            << "    Created " << faceID.size()
            << " separate face zones from cell zone " << cellZoneName << nl;

661
        forAll(names, i)
662
663
        {
            label nFaces = returnReduce(faceID[i].size(), sumOp<label>());
664
            Info<< "    " << names[i] << ": "
665
666
667
668
669
670
671
672
                << nFaces << " faces" << nl;
        }

        Info<< endl;
    }
}


673
674
Foam::scalar Foam::functionObjects::fluxSummary::totalArea
(
675
    const label idx
676
) const
677
{
678
679
680
681
    scalar sumMagSf = 0;

    if (isSurfaceMode())
    {
682
683
684
        const polySurface& s =
            storedObjects().lookupObject<polySurface>(zoneNames_[idx]);

685
686
687
688
689
        sumMagSf = sum(s.magSf());
    }
    else
    {
        const surfaceScalarField& magSf = mesh_.magSf();
690

691
692
        const labelList& faceIDs = faceID_[idx];
        const labelList& facePatchIDs = facePatchID_[idx];
693

694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
        forAll(faceIDs, i)
        {
            label facei = faceIDs[i];

            if (facePatchIDs[i] == -1)
            {
                sumMagSf += magSf[facei];
            }
            else
            {
                label patchi = facePatchIDs[i];
                sumMagSf += magSf.boundaryField()[patchi][facei];
            }
        }
    }

    return returnReduce(sumMagSf, sumOp<scalar>());
}


bool Foam::functionObjects::fluxSummary::surfaceModeWrite()
{
716
    for (const word& surfName : zoneNames_)
717
    {
718
719
720
721
        const polySurface& s =
            storedObjects().lookupObject<polySurface>(surfName);

        const auto& phi = s.lookupObject<polySurfaceVectorField>(phiName_);
722
723
724
725
726
727
728

        Log << type() << ' ' << name() << ' '
            << checkFlowType(phi.dimensions(), phi.name()) << " write:" << nl;
    }


    forAll(zoneNames_, surfi)
729
    {
730
731
732
733
        const polySurface& s =
            storedObjects().lookupObject<polySurface>(zoneNames_[surfi]);

        const auto& phi = s.lookupObject<polySurfaceVectorField>(phiName_);
734
735

        checkFlowType(phi.dimensions(), phi.name());
736

737
738
        const boolList& flips = faceFlip_[surfi];

739
740
        scalar phiPos(0);
        scalar phiNeg(0);
741
742
743
744
745

        tmp<scalarField> tphis = phi & s.Sf();
        const scalarField& phis = tphis();

        forAll(s, i)
746
        {
747
748
749
750
751
752
753
754
755
756
757
758
759
760
            scalar phif = phis[i];
            if (flips[i])
            {
                phif *= -1;
            }

            if (phif > 0)
            {
                phiPos += phif;
            }
            else
            {
                phiNeg += phif;
            }
761
        }
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779

        reduce(phiPos, sumOp<scalar>());
        reduce(phiNeg, sumOp<scalar>());

        phiPos *= scaleFactor_;
        phiNeg *= scaleFactor_;

        scalar netFlux = phiPos + phiNeg;
        scalar absoluteFlux = phiPos - phiNeg;

        Log << "    surface " << zoneNames_[surfi] << ':' << nl
            << "        positive : " << phiPos << nl
            << "        negative : " << phiNeg << nl
            << "        net      : " << netFlux << nl
            << "        absolute : " << absoluteFlux
            << nl << endl;

        if (writeToFile())
780
        {
781
782
783
784
785
786
787
            filePtrs_[surfi]
                << time_.value() << token::TAB
                << phiPos << token::TAB
                << phiNeg << token::TAB
                << netFlux << token::TAB
                << absoluteFlux
                << endl;
788
789
        }
    }
790

791
792
793
    Log << endl;

    return true;
794
795
796
}


797
bool Foam::functionObjects::fluxSummary::update()
798
{
799
800
801
802
    if (!needsUpdate_)
    {
        return false;
    }
803

804
805
806
807
808
809
    // Initialise with capacity == number of input names
    DynamicList<word> faceZoneName(zoneNames_.size());
    DynamicList<vector>    refDir(faceZoneName.capacity());
    DynamicList<labelList> faceID(faceZoneName.capacity());
    DynamicList<labelList> facePatchID(faceZoneName.capacity());
    DynamicList<boolList>  faceFlips(faceZoneName.capacity());
810

811
812
813
    switch (mode_)
    {
        case mdFaceZone:
814
        {
815
            forAll(zoneNames_, zonei)
816
            {
817
818
                initialiseFaceZone
                (
819
                    zoneNames_[zonei],
820
                    faceZoneName,
821
                    refDir, // fill with dummy value
822
823
                    faceID,
                    facePatchID,
824
                    faceFlips
825
                );
826
            }
827
828
829
830
            break;
        }
        case mdFaceZoneAndDirection:
        {
831
            forAll(zoneNames_, zonei)
832
            {
833
                initialiseFaceZoneAndDirection
834
                (
835
836
                    zoneNames_[zonei],
                    zoneDirections_[zonei],
837
                    faceZoneName,
838
                    refDir,
839
840
                    faceID,
                    facePatchID,
841
                    faceFlips
842
                );
843
            }
844
            break;
845
        }
846
        case mdCellZoneAndDirection:
847
        {
848
            forAll(zoneNames_, zonei)
849
            {
850
                initialiseCellZoneAndDirection
851
                (
852
853
                    zoneNames_[zonei],
                    zoneDirections_[zonei],
854
                    faceZoneName,
855
                    refDir,
856
857
                    faceID,
                    facePatchID,
858
                    faceFlips
859
860
                );
            }
861
862
            break;
        }
863
864
        case mdSurface:
        {
865
            forAll(zoneNames_, zonei)
866
867
868
            {
                initialiseSurface
                (
869
                    zoneNames_[zonei],
870
871
872
873
874
875
876
877
878
                    faceZoneName,
                    refDir,
                    faceFlips
                );
            }
            break;
        }
        case mdSurfaceAndDirection:
        {
879
            forAll(zoneNames_, zonei)
880
881
882
            {
                initialiseSurfaceAndDirection
                (
883
884
                    zoneNames_[zonei],
                    zoneDirections_[zonei],
885
886
887
888
889
890
891
                    faceZoneName,
                    refDir,
                    faceFlips
                );
            }
            break;
        }
892
893

        // Compiler warning if we forgot an enumeration
894
    }
895

896
    zoneNames_.transfer(faceZoneName);
897
898
    faceID_.transfer(faceID);
    facePatchID_.transfer(facePatchID);
899
    faceFlip_.transfer(faceFlips);
900

901
902
903
904
905
906
907
908
909
    Info<< type() << " " << name() << " output:" << nl;

    // Calculate and report areas
    List<scalar> areas(zoneNames_.size());
    forAll(zoneNames_, zonei)
    {
        const word& zoneName = zoneNames_[zonei];
        areas[zonei] = totalArea(zonei);

910
911
912
913
914
915
916
917
918
919
        if (isSurfaceMode())
        {
            Info<< "    Surface: " << zoneName
                << ", area: " << areas[zonei] << nl;
        }
        else
        {
            Info<< "    Zone: " << zoneName
                << ", area: " << areas[zonei] << nl;
        }
920
921
    }
    Info<< endl;
922
923
924

    if (writeToFile())
    {
925
        filePtrs_.resize(zoneNames_.size());
926

927
        forAll(filePtrs_, zonei)
928
        {
929
930
            const word& zoneName = zoneNames_[zonei];
            filePtrs_.set(zonei, createFile(zoneName));
931
932
            writeFileHeader
            (
933
934
935
936
                zoneName,
                areas[zonei],
                refDir[zonei],
                filePtrs_[zonei]
937
938
939
940
            );
        }
    }

941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
    Info<< endl;

    needsUpdate_ = false;

    return true;
}


// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

Foam::functionObjects::fluxSummary::fluxSummary
(
    const word& name,
    const Time& runTime,
    const dictionary& dict
)
:
    fvMeshFunctionObject(name, runTime, dict),
    writeFile(obr_, name),
    needsUpdate_(true),
    mode_(mdFaceZone),
    scaleFactor_(1),
    phiName_("phi"),
    zoneNames_(),
    zoneDirections_(),
    faceID_(),
    facePatchID_(),
    faceFlip_(),
    filePtrs_(),
    tolerance_(0.8)
{
    read(dict);
}


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

bool Foam::functionObjects::fluxSummary::read(const dictionary& dict)
{
    fvMeshFunctionObject::read(dict);
    writeFile::read(dict);

    needsUpdate_ = true;
    mode_ = modeTypeNames_.get("mode", dict);
985
986
    phiName_ = dict.getOrDefault<word>("phi", "phi");
    scaleFactor_ = dict.getOrDefault<scalar>("scaleFactor", 1);
987
    tolerance_ = dict.getOrDefault<scalar>("tolerance", 0.8);
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052

    zoneNames_.clear();
    zoneDirections_.clear();

    List<Tuple2<word, vector>> nameAndDirection;

    switch (mode_)
    {
        case mdFaceZone:
        {
            dict.readEntry("faceZones", zoneNames_);
            break;
        }
        case mdFaceZoneAndDirection:
        {
            dict.readEntry("faceZoneAndDirection", nameAndDirection);
            break;
        }
        case mdCellZoneAndDirection:
        {
            dict.readEntry("cellZoneAndDirection", nameAndDirection);
            break;
        }
        case mdSurface:
        {
            dict.readEntry("surfaces", zoneNames_);
            break;
        }
        case mdSurfaceAndDirection:
        {
            dict.readEntry("surfaceAndDirection", nameAndDirection);
            break;
        }
        default:
        {
            FatalIOErrorInFunction(dict)
                << "unhandled enumeration " << modeTypeNames_[mode_]
                << abort(FatalIOError);
        }
    }


    // Split name/vector into separate lists
    if (nameAndDirection.size())
    {
        zoneNames_.resize(nameAndDirection.size());
        zoneDirections_.resize(nameAndDirection.size());

        label zonei = 0;

        for (const Tuple2<word, vector>& nameDirn : nameAndDirection)
        {
            zoneNames_[zonei] = nameDirn.first();
            zoneDirections_[zonei] = nameDirn.second();
            ++zonei;
        }

        nameAndDirection.clear();
    }


    Info<< type() << ' ' << name() << " ("
        << modeTypeNames_[mode_] << ") with selection:\n    "
        << flatOutput(zoneNames_) << endl;

1053
    return !zoneNames_.empty();
1054
1055
1056
}


1057
void Foam::functionObjects::fluxSummary::writeFileHeader
1058
(
1059
    const word& zoneName,
1060
1061
1062
1063
1064
1065
    const scalar area,
    const vector& refDir,
    Ostream& os
) const
{
    writeHeader(os, "Flux summary");
1066
1067
1068
1069
1070
1071
1072
1073
    if (isSurfaceMode())
    {
        writeHeaderValue(os, "Surface", zoneName);
    }
    else
    {
        writeHeaderValue(os, "Face zone", zoneName);
    }
1074
1075
1076
1077
1078
1079
    writeHeaderValue(os, "Total area", area);

    switch (mode_)
    {
        case mdFaceZoneAndDirection:
        case mdCellZoneAndDirection:
1080
        case mdSurfaceAndDirection:
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
        {
            writeHeaderValue(os, "Reference direction", refDir);
            break;
        }
        default:
        {}
    }

    writeHeaderValue(os, "Scale factor", scaleFactor_);

    writeCommented(os, "Time");
    os  << tab << "positive"
        << tab << "negative"
        << tab << "net"
        << tab << "absolute"
        << endl;
}