surfaceFieldValue.C 27.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) 2011-2017 OpenFOAM Foundation
9
    Copyright (C) 2017-2020 OpenCFD Ltd.
10
11
12
13
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

14
15
16
17
    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.
18
19
20
21
22
23
24

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

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

29
#include "surfaceFieldValue.H"
30
31
#include "fvMesh.H"
#include "emptyPolyPatch.H"
32
#include "coupledPolyPatch.H"
33
#include "sampledSurface.H"
34
35
#include "mergePoints.H"
#include "indirectPrimitivePatch.H"
36
#include "PatchTools.H"
37
#include "addToRunTimeSelectionTable.H"
38
39
40

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

41
namespace Foam
42
{
43
44
45
46
namespace functionObjects
{
namespace fieldValues
{
47
48
49
    defineTypeNameAndDebug(surfaceFieldValue, 0);
    addToRunTimeSelectionTable(fieldValue, surfaceFieldValue, dictionary);
    addToRunTimeSelectionTable(functionObject, surfaceFieldValue, dictionary);
50
51
}
}
52
53
54
}


55
const Foam::Enum
56
<
57
58
59
    Foam::functionObjects::fieldValues::surfaceFieldValue::regionTypes
>
Foam::functionObjects::fieldValues::surfaceFieldValue::regionTypeNames_
Mark OLESEN's avatar
Mark OLESEN committed
60
({
61
62
    { regionTypes::stFaceZone, "faceZone" },
    { regionTypes::stPatch, "patch" },
63
    { regionTypes::stObject, "functionObjectSurface" },
64
    { regionTypes::stSampled, "sampledSurface" },
Mark OLESEN's avatar
Mark OLESEN committed
65
});
66

67
68

const Foam::Enum
69
<
70
71
72
    Foam::functionObjects::fieldValues::surfaceFieldValue::operationType
>
Foam::functionObjects::fieldValues::surfaceFieldValue::operationTypeNames_
Mark OLESEN's avatar
Mark OLESEN committed
73
({
74
    // Normal operations
75
    { operationType::opNone, "none" },
76
77
    { operationType::opMin, "min" },
    { operationType::opMax, "max" },
78
79
80
81
82
83
84
85
86
87
    { operationType::opSum, "sum" },
    { operationType::opSumMag, "sumMag" },
    { operationType::opSumDirection, "sumDirection" },
    { operationType::opSumDirectionBalance, "sumDirectionBalance" },
    { operationType::opAverage, "average" },
    { operationType::opAreaAverage, "areaAverage" },
    { operationType::opAreaIntegrate, "areaIntegrate" },
    { operationType::opCoV, "CoV" },
    { operationType::opAreaNormalAverage, "areaNormalAverage" },
    { operationType::opAreaNormalIntegrate, "areaNormalIntegrate" },
88
    { operationType::opUniformity, "uniformity" },
89
90
91
92
93
94

    // Using weighting
    { operationType::opWeightedSum, "weightedSum" },
    { operationType::opWeightedAverage, "weightedAverage" },
    { operationType::opWeightedAreaAverage, "weightedAreaAverage" },
    { operationType::opWeightedAreaIntegrate, "weightedAreaIntegrate" },
95
    { operationType::opWeightedUniformity, "weightedUniformity" },
96
97
98
99
100
101

    // Using absolute weighting
    { operationType::opAbsWeightedSum, "absWeightedSum" },
    { operationType::opAbsWeightedAverage, "absWeightedAverage" },
    { operationType::opAbsWeightedAreaAverage, "absWeightedAreaAverage" },
    { operationType::opAbsWeightedAreaIntegrate, "absWeightedAreaIntegrate" },
102
    { operationType::opAbsWeightedUniformity, "absWeightedUniformity" },
Mark OLESEN's avatar
Mark OLESEN committed
103
});
104

105
const Foam::Enum
106
<
107
    Foam::functionObjects::fieldValues::surfaceFieldValue::postOperationType
108
>
109
Foam::functionObjects::fieldValues::surfaceFieldValue::postOperationTypeNames_
Mark OLESEN's avatar
Mark OLESEN committed
110
({
111
    { postOperationType::postOpNone, "none" },
112
    { postOperationType::postOpMag, "mag" },
113
    { postOperationType::postOpSqrt, "sqrt" },
Mark OLESEN's avatar
Mark OLESEN committed
114
});
115

116

Henry Weller's avatar
Henry Weller committed
117
// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
118

119
120
121
const Foam::objectRegistry&
Foam::functionObjects::fieldValues::surfaceFieldValue::obr() const
{
122
    if (stObject == regionType_)
123
    {
124
        return storedObjects().lookupObject<polySurface>(regionName_);
125
    }
126
127

    return mesh_;
128
129
130
}


131
void Foam::functionObjects::fieldValues::surfaceFieldValue::setFaceZoneFaces()
132
{
133
    const label zoneId = mesh_.faceZones().findZoneID(regionName_);
134
135
136

    if (zoneId < 0)
    {
137
        FatalErrorInFunction
138
            << type() << " " << name() << ": "
139
            << regionTypeNames_[regionType_] << '(' << regionName_ << "):" << nl
140
            << "    Unknown face zone name: " << regionName_
141
            << ". Valid face zones are: " << mesh_.faceZones().names()
142
143
144
            << nl << exit(FatalError);
    }

145
    const faceZone& fZone = mesh_.faceZones()[zoneId];
146

147
148
    DynamicList<label> faceIds(fZone.size());
    DynamicList<label> facePatchIds(fZone.size());
149
    DynamicList<bool> faceFlip(fZone.size());
150
151
152

    forAll(fZone, i)
    {
153
        const label facei = fZone[i];
154
155
156

        label faceId = -1;
        label facePatchId = -1;
157
        if (mesh_.isInternalFace(facei))
158
        {
159
            faceId = facei;
160
161
162
163
            facePatchId = -1;
        }
        else
        {
164
165
            facePatchId = mesh_.boundaryMesh().whichPatch(facei);
            const polyPatch& pp = mesh_.boundaryMesh()[facePatchId];
166
            if (isA<coupledPolyPatch>(pp))
167
            {
168
                if (refCast<const coupledPolyPatch>(pp).owner())
169
                {
170
                    faceId = pp.whichFace(facei);
171
172
173
174
175
176
177
178
                }
                else
                {
                    faceId = -1;
                }
            }
            else if (!isA<emptyPolyPatch>(pp))
            {
179
                faceId = facei - pp.start();
180
181
182
183
184
185
186
187
188
189
            }
            else
            {
                faceId = -1;
                facePatchId = -1;
            }
        }

        if (faceId >= 0)
        {
190
191
            faceIds.append(faceId);
            facePatchIds.append(facePatchId);
192
            faceFlip.append(fZone.flipMap()[i] ? true : false);
193
194
195
        }
    }

196
197
    faceId_.transfer(faceIds);
    facePatchId_.transfer(facePatchIds);
198
    faceFlip_.transfer(faceFlip);
199
    nFaces_ = returnReduce(faceId_.size(), sumOp<label>());
200
201
202

    if (debug)
    {
203
        Pout<< "Original face zone size = " << fZone.size()
204
            << ", new size = " << faceId_.size() << endl;
205
206
207
208
    }
}


209
void Foam::functionObjects::fieldValues::surfaceFieldValue::setPatchFaces()
210
{
211
    const label patchid = mesh_.boundaryMesh().findPatchID(regionName_);
212

213
    if (patchid < 0)
214
    {
215
        FatalErrorInFunction
216
            << type() << " " << name() << ": "
217
            << regionTypeNames_[regionType_] << '(' << regionName_ << "):" << nl
218
            << "    Unknown patch name: " << regionName_
219
            << ". Valid patch names are: "
220
            << mesh_.boundaryMesh().names() << nl
221
222
223
            << exit(FatalError);
    }

224
    const polyPatch& pp = mesh_.boundaryMesh()[patchid];
225
226

    label nFaces = pp.size();
227
    if (isA<emptyPolyPatch>(pp))
228
229
230
231
232
233
    {
        nFaces = 0;
    }

    faceId_.setSize(nFaces);
    facePatchId_.setSize(nFaces);
234
    faceFlip_.setSize(nFaces);
235
    nFaces_ = returnReduce(faceId_.size(), sumOp<label>());
236

237
    forAll(faceId_, facei)
238
    {
239
240
        faceId_[facei] = facei;
        facePatchId_[facei] = patchid;
241
        faceFlip_[facei] = false;
242
243
244
245
    }
}


246
void Foam::functionObjects::fieldValues::surfaceFieldValue::combineMeshGeometry
247
248
249
250
251
252
253
254
255
256
257
258
259
(
    faceList& faces,
    pointField& points
) const
{
    List<faceList> allFaces(Pstream::nProcs());
    List<pointField> allPoints(Pstream::nProcs());

    labelList globalFacesIs(faceId_);
    forAll(globalFacesIs, i)
    {
        if (facePatchId_[i] != -1)
        {
260
            const label patchi = facePatchId_[i];
261
            globalFacesIs[i] += mesh_.boundaryMesh()[patchi].start();
262
263
264
265
266
267
        }
    }

    // Add local faces and points to the all* lists
    indirectPrimitivePatch pp
    (
268
269
        IndirectList<face>(mesh_.faces(), globalFacesIs),
        mesh_.points()
270
271
272
273
274
275
276
277
278
279
    );
    allFaces[Pstream::myProcNo()] = pp.localFaces();
    allPoints[Pstream::myProcNo()] = pp.localPoints();

    Pstream::gatherList(allFaces);
    Pstream::gatherList(allPoints);

    // Renumber and flatten
    label nFaces = 0;
    label nPoints = 0;
280
    forAll(allFaces, proci)
281
    {
282
283
        nFaces += allFaces[proci].size();
        nPoints += allPoints[proci].size();
284
285
286
287
288
289
290
291
292
293
294
    }

    faces.setSize(nFaces);
    points.setSize(nPoints);

    nFaces = 0;
    nPoints = 0;

    // My own data first
    {
        const faceList& fcs = allFaces[Pstream::myProcNo()];
295
        for (const face& f : fcs)
296
297
298
299
300
301
302
303
304
305
        {
            face& newF = faces[nFaces++];
            newF.setSize(f.size());
            forAll(f, fp)
            {
                newF[fp] = f[fp] + nPoints;
            }
        }

        const pointField& pts = allPoints[Pstream::myProcNo()];
306
        for (const point& pt : pts)
307
        {
308
            points[nPoints++] = pt;
309
310
311
312
        }
    }

    // Other proc data follows
313
    forAll(allFaces, proci)
314
    {
315
        if (proci != Pstream::myProcNo())
316
        {
317
            const faceList& fcs = allFaces[proci];
318
            for (const face& f : fcs)
319
320
321
322
323
324
325
326
327
            {
                face& newF = faces[nFaces++];
                newF.setSize(f.size());
                forAll(f, fp)
                {
                    newF[fp] = f[fp] + nPoints;
                }
            }

328
            const pointField& pts = allPoints[proci];
329
            for (const point& pt : pts)
330
            {
331
                points[nPoints++] = pt;
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
            }
        }
    }

    // Merge
    labelList oldToNew;
    pointField newPoints;
    bool hasMerged = mergePoints
    (
        points,
        SMALL,
        false,
        oldToNew,
        newPoints
    );

    if (hasMerged)
    {
        if (debug)
        {
            Pout<< "Merged from " << points.size()
                << " down to " << newPoints.size() << " points" << endl;
        }

        points.transfer(newPoints);
357
        for (face& f : faces)
358
        {
359
            inplaceRenumber(oldToNew, f);
360
361
362
363
364
        }
    }
}


365
366
void Foam::functionObjects::fieldValues::surfaceFieldValue::
combineSurfaceGeometry
367
368
369
370
371
(
    faceList& faces,
    pointField& points
) const
{
372
    if (stObject == regionType_)
373
    {
374
        const polySurface& s = dynamicCast<const polySurface>(obr());
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401

        if (Pstream::parRun())
        {
            // Dimension as fraction of surface
            const scalar mergeDim = 1e-10*boundBox(s.points(), true).mag();

            labelList pointsMap;

            PatchTools::gatherAndMerge
            (
                mergeDim,
                primitivePatch
                (
                    SubList<face>(s.faces(), s.faces().size()),
                    s.points()
                ),
                points,
                faces,
                pointsMap
            );
        }
        else
        {
            faces = s.faces();
            points = s.points();
        }
    }
402
    else if (sampledPtr_.valid())
403
    {
404
        const sampledSurface& s = sampledPtr_();
405
406
407

        if (Pstream::parRun())
        {
408
            // Dimension as fraction of mesh bounding box
409
            const scalar mergeDim = 1e-10*mesh_.bounds().mag();
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434

            labelList pointsMap;

            PatchTools::gatherAndMerge
            (
                mergeDim,
                primitivePatch
                (
                    SubList<face>(s.faces(), s.faces().size()),
                    s.points()
                ),
                points,
                faces,
                pointsMap
            );
        }
        else
        {
            faces = s.faces();
            points = s.points();
        }
    }
}


435
Foam::scalar
436
Foam::functionObjects::fieldValues::surfaceFieldValue::totalArea() const
437
{
438
    scalar totalArea = 0;
439

440
    if (stObject == regionType_)
441
    {
442
        const polySurface& s = dynamicCast<const polySurface>(obr());
443
444
445

        totalArea = gSum(s.magSf());
    }
446
    else if (sampledPtr_.valid())
447
    {
448
        totalArea = gSum(sampledPtr_().magSf());
449
450
451
    }
    else
    {
452
        totalArea = gSum(filterField(mesh_.magSf()));
453
454
455
456
457
458
    }

    return totalArea;
}


459
460
// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //

461
bool Foam::functionObjects::fieldValues::surfaceFieldValue::usesSf() const
462
{
463
    // Only a few operations do not require the Sf field
464
465
466
    switch (operation_)
    {
        case opNone:
467
468
        case opMin:
        case opMax:
469
470
471
        case opSum:
        case opSumMag:
        case opAverage:
472
        {
473
            return false;
474
        }
475
        default:
476
        {
477
            return true;
478
        }
479
480
481
482
    }
}


483
bool Foam::functionObjects::fieldValues::surfaceFieldValue::update()
484
{
485
486
487
488
    if (sampledPtr_.valid())
    {
        sampledPtr_->update();
    }
489

490
491
492
493
    if (!needsUpdate_)
    {
        return false;
    }
494

495
    switch (regionType_)
496
497
498
499
500
501
502
503
504
    {
        case stFaceZone:
        {
            setFaceZoneFaces();
            break;
        }
        case stPatch:
        {
            setPatchFaces();
505
506
            break;
        }
507
        case stObject:
508
        {
509
            const polySurface& s = dynamicCast<const polySurface>(obr());
510
            nFaces_ = returnReduce(s.size(), sumOp<label>());
511
512
            break;
        }
513
        case stSampled:
514
        {
515
            nFaces_ = returnReduce(sampledPtr_->faces().size(), sumOp<label>());
516
517
            break;
        }
518
519

        // Compiler warning if we forgot an enumeration
520
521
    }

522
523
    if (nFaces_ == 0)
    {
524
        FatalErrorInFunction
525
            << type() << " " << name() << ": "
526
            << regionTypeNames_[regionType_] << '(' << regionName_ << "):" << nl
527
            << "    Region has no faces" << exit(FatalError);
528
529
    }

530
531
    totalArea_ = totalArea();

Andrew Heather's avatar
Andrew Heather committed
532
533
534
    Log << "    total faces   = " << nFaces_ << nl
        << "    total area    = " << totalArea_ << nl
        << endl;
535

536
    writeFileHeader(file());
537

538
539
    needsUpdate_ = false;
    return true;
540
541
542
}


543
544
545
void Foam::functionObjects::fieldValues::surfaceFieldValue::writeFileHeader
(
    Ostream& os
546
)
547
{
548
    if (canWriteHeader() && (operation_ != opNone))
549
    {
550
551
        writeCommented(os, "Region type : ");
        os << regionTypeNames_[regionType_] << " " << regionName_ << endl;
552

553
554
555
556
        writeHeaderValue(os, "Faces", nFaces_);
        writeHeaderValue(os, "Area", totalArea_);
        writeHeaderValue(os, "Scale factor", scaleFactor_);

557
558
559
560
561
        if (weightFieldName_ != "none")
        {
            writeHeaderValue(os, "Weight field", weightFieldName_);
        }

562
        writeCommented(os, "Time");
563
564
        if (writeArea_)
        {
565
            os  << tab << "Area";
566
        }
567

568
569
        // TBD: add in postOperation information?

570
        for (const word& fieldName : fields_)
571
        {
572
            os  << tab << operationTypeNames_[operation_]
573
                << "(" << fieldName << ")";
574
        }
575

576
        os  << endl;
577
    }
578
579

    writtenHeader_ = true;
580
581
582
}


583
template<>
584
585
Foam::scalar
Foam::functionObjects::fieldValues::surfaceFieldValue::processValues
586
587
588
589
590
591
592
593
594
595
(
    const Field<scalar>& values,
    const vectorField& Sf,
    const scalarField& weightField
) const
{
    switch (operation_)
    {
        case opSumDirection:
        {
596
            const vector n(dict_.get<vector>("direction"));
Henry Weller's avatar
Henry Weller committed
597
            return gSum(pos0(values*(Sf & n))*mag(values));
598
599
600
        }
        case opSumDirectionBalance:
        {
601
            const vector n(dict_.get<vector>("direction"));
602
            const scalarField nv(values*(Sf & n));
603

Henry Weller's avatar
Henry Weller committed
604
            return gSum(pos0(nv)*mag(values) - neg(nv)*mag(values));
605
        }
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644

        case opUniformity:
        case opWeightedUniformity:
        case opAbsWeightedUniformity:
        {
            const scalar areaTotal = gSum(mag(Sf));
            tmp<scalarField> areaVal(values * mag(Sf));

            scalar mean, numer;

            if (canWeight(weightField))
            {
                // Weighted quantity = (Weight * phi * dA)

                tmp<scalarField> weight(weightingFactor(weightField));

                // Mean weighted value (area-averaged)
                mean = gSum(weight()*areaVal()) / areaTotal;

                // Abs. deviation from weighted mean value
                numer = gSum(mag(weight*areaVal - (mean * mag(Sf))));
            }
            else
            {
                // Unweighted quantity = (1 * phi * dA)

                // Mean value (area-averaged)
                mean = gSum(areaVal()) / areaTotal;

                // Abs. deviation from mean value
                numer = gSum(mag(areaVal - (mean * mag(Sf))));
            }

            // Uniformity index
            const scalar ui = 1 - numer/(2*mag(mean*areaTotal) + ROOTVSMALL);

            return min(max(ui, 0), 1);
        }

645
646
647
648
649
650
651
652
653
        default:
        {
            // Fall through to other operations
            return processSameTypeValues(values, Sf, weightField);
        }
    }
}


654
template<>
655
656
Foam::vector
Foam::functionObjects::fieldValues::surfaceFieldValue::processValues
657
658
659
660
661
662
663
664
(
    const Field<vector>& values,
    const vectorField& Sf,
    const scalarField& weightField
) const
{
    switch (operation_)
    {
665
666
        case opSumDirection:
        {
667
            const vector n(dict_.get<vector>("direction").normalise());
668

669
            const scalarField nv(n & values);
Henry Weller's avatar
Henry Weller committed
670
            return gSum(pos0(nv)*n*(nv));
671
672
673
        }
        case opSumDirectionBalance:
        {
674
            const vector n(dict_.get<vector>("direction").normalise());
675

676
            const scalarField nv(n & values);
Henry Weller's avatar
Henry Weller committed
677
            return gSum(pos0(nv)*n*(nv));
678
        }
679
680
        case opAreaNormalAverage:
        {
681
682
            const scalar val = gSum(values & Sf)/gSum(mag(Sf));
            return vector(val, 0, 0);
683
684
685
        }
        case opAreaNormalIntegrate:
        {
686
687
            const scalar val = gSum(values & Sf);
            return vector(val, 0, 0);
688
        }
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727

        case opUniformity:
        case opWeightedUniformity:
        case opAbsWeightedUniformity:
        {
            const scalar areaTotal = gSum(mag(Sf));
            tmp<scalarField> areaVal(values & Sf);

            scalar mean, numer;

            if (canWeight(weightField))
            {
                // Weighted quantity = (Weight * phi . dA)

                tmp<scalarField> weight(weightingFactor(weightField));

                // Mean weighted value (area-averaged)
                mean = gSum(weight()*areaVal()) / areaTotal;

                // Abs. deviation from weighted mean value
                numer = gSum(mag(weight*areaVal - (mean * mag(Sf))));
            }
            else
            {
                // Unweighted quantity = (1 * phi . dA)

                // Mean value (area-averaged)
                mean = gSum(areaVal()) / areaTotal;

                // Abs. deviation from mean value
                numer = gSum(mag(areaVal - (mean * mag(Sf))));
            }

            // Uniformity index
            const scalar ui = 1 - numer/(2*mag(mean*areaTotal) + ROOTVSMALL);

            return vector(min(max(ui, 0), 1), 0, 0);
        }

728
729
730
731
732
733
734
735
736
        default:
        {
            // Fall through to other operations
            return processSameTypeValues(values, Sf, weightField);
        }
    }
}


737
738
739
740
741
template<>
Foam::tmp<Foam::scalarField>
Foam::functionObjects::fieldValues::surfaceFieldValue::weightingFactor
(
    const Field<scalar>& weightField
742
) const
743
{
744
745
746
747
    if (usesMag())
    {
        return mag(weightField);
    }
748
749
750

    // pass through
    return weightField;
751
752
753
754
755
756
757
758
759
}


template<>
Foam::tmp<Foam::scalarField>
Foam::functionObjects::fieldValues::surfaceFieldValue::weightingFactor
(
    const Field<scalar>& weightField,
    const vectorField& Sf
760
) const
761
762
{
    // scalar * Area
763

764
765
    if (returnReduce(weightField.empty(), andOp<bool>()))
    {
766
        // No weight field - revert to unweighted form
767
768
        return mag(Sf);
    }
769
770
771
772
    else if (usesMag())
    {
        return mag(weightField * mag(Sf));
    }
773
774

    return (weightField * mag(Sf));
775
776
777
778
779
780
781
782
783
}


template<>
Foam::tmp<Foam::scalarField>
Foam::functionObjects::fieldValues::surfaceFieldValue::weightingFactor
(
    const Field<vector>& weightField,
    const vectorField& Sf
784
) const
785
786
{
    // vector (dot) Area
787

788
789
    if (returnReduce(weightField.empty(), andOp<bool>()))
    {
790
        // No weight field - revert to unweighted form
791
792
        return mag(Sf);
    }
793
794
795
796
    else if (usesMag())
    {
        return mag(weightField & Sf);
    }
797
798

    return (weightField & Sf);
799
800
801
}


802
803
// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

804
Foam::functionObjects::fieldValues::surfaceFieldValue::surfaceFieldValue
805
806
807
808
809
810
811
(
    const word& name,
    const Time& runTime,
    const dictionary& dict
)
:
    fieldValue(name, runTime, dict, typeName),
Mark OLESEN's avatar
Mark OLESEN committed
812
813
    regionType_(regionTypeNames_.get("regionType", dict)),
    operation_(operationTypeNames_.get("operation", dict)),
814
815
    postOperation_
    (
816
        postOperationTypeNames_.getOrDefault
817
818
819
        (
            "postOperation",
            dict,
Mark OLESEN's avatar
Mark OLESEN committed
820
821
            postOperationType::postOpNone,
            true  // Failsafe behaviour
822
        )
823
    ),
824
    weightFieldName_("none"),
825
826
    needsUpdate_(true),
    writeArea_(false),
827
    totalArea_(0),
828
829
830
    nFaces_(0),
    faceId_(),
    facePatchId_(),
831
    faceFlip_()
832
833
834
835
{
    read(dict);
}

836

837
Foam::functionObjects::fieldValues::surfaceFieldValue::surfaceFieldValue
838
839
840
(
    const word& name,
    const objectRegistry& obr,
841
    const dictionary& dict
842
843
)
:
844
    fieldValue(name, obr, dict, typeName),
Mark OLESEN's avatar
Mark OLESEN committed
845
846
    regionType_(regionTypeNames_.get("regionType", dict)),
    operation_(operationTypeNames_.get("operation", dict)),
847
848
    postOperation_
    (
849
        postOperationTypeNames_.getOrDefault
850
851
852
        (
            "postOperation",
            dict,
Mark OLESEN's avatar
Mark OLESEN committed
853
854
            postOperationType::postOpNone,
            true  // Failsafe behaviour
855
        )
856
    ),
857
    weightFieldName_("none"),
858
859
    needsUpdate_(true),
    writeArea_(false),
860
    totalArea_(0),
861
    nFaces_(0),
862
863
    faceId_(),
    facePatchId_(),
864
    faceFlip_()
865
{
866
    read(dict);
867
868
869
870
871
}


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

872
bool Foam::functionObjects::fieldValues::surfaceFieldValue::read
873
874
875
(
    const dictionary& dict
)
876
{
Andrew Heather's avatar
Andrew Heather committed
877
    fieldValue::read(dict);
878
879
880

    weightFieldName_ = "none";
    needsUpdate_ = true;
881
    writeArea_ = dict.getOrDefault("writeArea", false);
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
    totalArea_ = 0;
    nFaces_ = 0;
    faceId_.clear();
    facePatchId_.clear();
    faceFlip_.clear();
    sampledPtr_.clear();
    surfaceWriterPtr_.clear();

    dict.readEntry("name", regionName_);

    // Create sampled surface, but leave 'expired' (ie, no update) since it
    // may depend on fields or data that do not yet exist
    if (stSampled == regionType_)
    {
        sampledPtr_ = sampledSurface::New
        (
            name(),
            mesh_,
            dict.subDict("sampledSurfaceDict")
        );
    }

    Info<< type() << " " << name() << ":" << nl
        << "    operation     = ";

    if (postOperation_ != postOpNone)
    {
        Info<< postOperationTypeNames_[postOperation_] << '('
            << operationTypeNames_[operation_] << ')'  << nl;
    }
    else
    {
        Info<< operationTypeNames_[operation_] << nl;
    }

    if (usesWeight())
    {
        if (stSampled == regionType_)
        {
            FatalIOErrorInFunction(dict)
                << "Cannot use weighted operation '"
                << operationTypeNames_[operation_]
                << "' for sampledSurface"
                << exit(FatalIOError);
        }

        if (dict.readIfPresent("weightField", weightFieldName_))
        {
            Info<< "    weight field  = " << weightFieldName_ << nl;
        }
        else
        {
            // Suggest possible alternative unweighted operation?
            FatalIOErrorInFunction(dict)
                << "The '" << operationTypeNames_[operation_]
                << "' operation is missing a weightField." << nl
                << "Either provide the weightField, "
                << "use weightField 'none' to suppress weighting," << nl
                << "or use a different operation."
                << exit(FatalIOError);
        }
    }

    // Backwards compatibility for v1612 and older
    List<word> orientedFields;
    if (dict.readIfPresent("orientedFields", orientedFields))
    {
        fields_.append(orientedFields);

        WarningInFunction
            << "The 'orientedFields' option is deprecated.  These fields can "
            << "and have been added to the standard 'fields' list."
            << endl;
    }

    if (writeFields_)
    {
        const word formatName(dict.get<word>("surfaceFormat"));

961
962
963
        surfaceWriterPtr_.reset
        (
            surfaceWriter::New
964
            (
965
966
967
968
                formatName,
                dict.subOrEmptyDict("formatOptions").subOrEmptyDict(formatName)
            )
        );
969

970
971
972
973
974
        if (debug)
        {
            surfaceWriterPtr_->verbose() = true;
        }

975
976
        if (surfaceWriterPtr_->enabled())
        {
977
978
            Info<< "    surfaceFormat = " << formatName << nl;
        }
979
980
981
982
        else
        {
            surfaceWriterPtr_->clear();
        }
983
984
985
    }

    Info<< nl << endl;
986

987
    return true;
988
989
990
}


991
bool Foam::functionObjects::fieldValues::surfaceFieldValue::write()
992
{
Andrew Heather's avatar
Andrew Heather committed
993
994
995
996
997
    if (needsUpdate_ || operation_ != opNone)
    {
        fieldValue::write();
    }

998
    update();
999

1000
    if (operation_ != opNone)
1001
    {
1002
        writeCurrentTime(file());
1003
    }
1004

1005
1006
1007
    if (writeArea_)
    {
        totalArea_ = totalArea();
1008
1009
        Log << "    total area = " << totalArea_ << endl;

1010
        if (operation_ != opNone && Pstream::master())
1011
        {
1012
            file() << tab << totalArea_;
1013
        }
1014
    }
1015

1016
1017
    // Many operations use the Sf field
    vectorField Sf;
1018
    if (usesSf())
1019
    {
1020
        if (stObject == regionType_)
1021
        {
1022
            const polySurface& s = dynamicCast<const polySurface>(obr());
1023
1024
            Sf = s.Sf();
        }
1025
        else if (sampledPtr_.valid())
1026
        {
1027
            Sf = sampledPtr_().Sf();
1028
1029
1030
        }
        else
        {
1031
            Sf = filterField(mesh_.Sf());
1032
1033
1034
1035
1036
1037
1038
1039
1040
        }
    }

    // Faces and points for surface format (if specified)
    faceList faces;
    pointField points;

    if (surfaceWriterPtr_.valid())
    {
1041
        if (withTopologicalMerge())
1042
        {
1043
            combineMeshGeometry(faces, points);
1044
1045
1046
        }
        else
        {
1047
            combineSurfaceGeometry(faces, points);
1048
1049
1050
1051
        }
    }

    // Only a few weight types (scalar, vector)
1052
1053
    if (weightFieldName_ != "none")
    {
1054
1055
        if (validField<scalar>(weightFieldName_))
        {
1056
            scalarField weightField
1057
            (
1058
                getFieldValues<scalar>(weightFieldName_, true)
1059
            );
1060

1061
            /