surfaceFieldValue.C 27.4 KB
Newer Older
1
2
3
4
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
5
    \\  /    A nd           | Copyright (C) 2017-2019 OpenCFD Ltd.
OpenFOAM bot's avatar
OpenFOAM bot committed
6
7
8
     \\/     M anipulation  |
-------------------------------------------------------------------------------
                            | Copyright (C) 2011-2017 OpenFOAM Foundation
9
10
11
12
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

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

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

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

28
#include "surfaceFieldValue.H"
29
30
#include "fvMesh.H"
#include "emptyPolyPatch.H"
31
#include "coupledPolyPatch.H"
32
#include "sampledSurface.H"
33
34
#include "mergePoints.H"
#include "indirectPrimitivePatch.H"
35
#include "PatchTools.H"
36
#include "meshedSurfRef.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
112
    { postOperationType::postOpNone, "none" },
    { postOperationType::postOpSqrt, "sqrt" },
Mark OLESEN's avatar
Mark OLESEN committed
113
});
114

115

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

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

    return mesh_;
127
128
129
}


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

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

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

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

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

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

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

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

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


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

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

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

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

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

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


245
void Foam::functionObjects::fieldValues::surfaceFieldValue::combineMeshGeometry
246
247
248
249
250
251
252
253
254
255
256
257
258
(
    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)
        {
259
            const label patchi = facePatchId_[i];
260
            globalFacesIs[i] += mesh_.boundaryMesh()[patchi].start();
261
262
263
264
265
266
        }
    }

    // Add local faces and points to the all* lists
    indirectPrimitivePatch pp
    (
267
268
        IndirectList<face>(mesh_.faces(), globalFacesIs),
        mesh_.points()
269
270
271
272
273
274
275
276
277
278
    );
    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;
279
    forAll(allFaces, proci)
280
    {
281
282
        nFaces += allFaces[proci].size();
        nPoints += allPoints[proci].size();
283
284
285
286
287
288
289
290
291
292
293
    }

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

    nFaces = 0;
    nPoints = 0;

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

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

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

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

    // 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);
356
        for (face& f : faces)
357
        {
358
            inplaceRenumber(oldToNew, f);
359
360
361
362
363
        }
    }
}


364
365
void Foam::functionObjects::fieldValues::surfaceFieldValue::
combineSurfaceGeometry
366
367
368
369
370
(
    faceList& faces,
    pointField& points
) const
{
371
    if (stObject == regionType_)
372
    {
373
        const polySurface& s = dynamicCast<const polySurface>(obr());
374
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

        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();
        }
    }
401
    else if (sampledPtr_.valid())
402
    {
403
        const sampledSurface& s = sampledPtr_();
404
405
406

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

            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();
        }
    }
}


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

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

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

    return totalArea;
}


458
459
// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //

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


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

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

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

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

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

531
532
    totalArea_ = totalArea();

533
534
    Log
        << "    total faces   = " << nFaces_ << nl
535
536
        << "    total area    = " << totalArea_ << nl;

537
    writeFileHeader(file());
538

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


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

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

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

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

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

575
        os  << endl;
576
    }
577
578
579
}


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

Henry Weller's avatar
Henry Weller committed
601
            return gSum(pos0(nv)*mag(values) - neg(nv)*mag(values));
602
        }
603
604
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

        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);
        }

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


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

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

673
            const scalarField nv(n & values);
Henry Weller's avatar
Henry Weller committed
674
            return gSum(pos0(nv)*n*(nv));
675
        }
676
677
        case opAreaNormalAverage:
        {
678
679
            const scalar val = gSum(values & Sf)/gSum(mag(Sf));
            return vector(val, 0, 0);
680
681
682
        }
        case opAreaNormalIntegrate:
        {
683
684
            const scalar val = gSum(values & Sf);
            return vector(val, 0, 0);
685
        }
686
687
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

        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);
        }

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


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

    // pass through
    return weightField;
748
749
750
751
752
753
754
755
756
}


template<>
Foam::tmp<Foam::scalarField>
Foam::functionObjects::fieldValues::surfaceFieldValue::weightingFactor
(
    const Field<scalar>& weightField,
    const vectorField& Sf
757
) const
758
759
760
761
{
    // scalar * Area
    if (returnReduce(weightField.empty(), andOp<bool>()))
    {
762
        // No weight field - revert to unweighted form
763
764
        return mag(Sf);
    }
765
766
767
768
    else if (usesMag())
    {
        return mag(weightField * mag(Sf));
    }
769
770

    return (weightField * mag(Sf));
771
772
773
774
775
776
777
778
779
}


template<>
Foam::tmp<Foam::scalarField>
Foam::functionObjects::fieldValues::surfaceFieldValue::weightingFactor
(
    const Field<vector>& weightField,
    const vectorField& Sf
780
) const
781
782
783
784
{
    // vector (dot) Area
    if (returnReduce(weightField.empty(), andOp<bool>()))
    {
785
        // No weight field - revert to unweighted form
786
787
        return mag(Sf);
    }
788
789
790
791
    else if (usesMag())
    {
        return mag(weightField & Sf);
    }
792
793

    return (weightField & Sf);
794
795
796
}


797
798
// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

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

831

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


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

867
bool Foam::functionObjects::fieldValues::surfaceFieldValue::read
868
869
870
(
    const dictionary& dict
)
871
{
Andrew Heather's avatar
Andrew Heather committed
872
    fieldValue::read(dict);
873
874
875
876
877
878
879
880
881
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

    weightFieldName_ = "none";
    needsUpdate_ = true;
    writeArea_ = dict.lookupOrDefault("writeArea", false);
    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"));

956
957
958
        surfaceWriterPtr_.reset
        (
            surfaceWriter::New
959
            (
960
961
962
963
                formatName,
                dict.subOrEmptyDict("formatOptions").subOrEmptyDict(formatName)
            )
        );
964

965
966
        if (surfaceWriterPtr_->enabled())
        {
967
968
            Info<< "    surfaceFormat = " << formatName << nl;
        }
969
970
971
972
        else
        {
            surfaceWriterPtr_->clear();
        }
973
974
975
    }

    Info<< nl << endl;
976

977
    return true;
978
979
980
}


981
bool Foam::functionObjects::fieldValues::surfaceFieldValue::write()
982
{
983
    update();
984

985
    if (operation_ != opNone)
986
    {
987
988
989
990
991
992
        fieldValue::write();

        if (Pstream::master())
        {
            writeTime(file());
        }
993
    }
994

995
996
997
    if (writeArea_)
    {
        totalArea_ = totalArea();
998
999
        Log << "    total area = " << totalArea_ << endl;

1000
        if (operation_ != opNone && Pstream::master())
1001
        {
1002
            file() << tab << totalArea_;
1003
        }
1004
    }
1005

1006
1007
    // Many operations use the Sf field
    vectorField Sf;
1008
    if (usesSf())
1009
    {
1010
        if (stObject == regionType_)
1011
        {
1012
            const polySurface& s = dynamicCast<const polySurface>(obr());
1013
1014
            Sf = s.Sf();
        }
1015
        else if (sampledPtr_.valid())
1016
        {
1017
            Sf = sampledPtr_().Sf();
1018
1019
1020
        }
        else
        {
1021
            Sf = filterField(mesh_.Sf());
1022
1023
1024
1025
1026
1027
1028
1029
1030
        }
    }

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

    if (surfaceWriterPtr_.valid())
    {
1031
        if (withTopologicalMerge())
1032
        {
1033
            combineMeshGeometry(faces, points);
1034
1035
1036
        }
        else
        {
1037
            combineSurfaceGeometry(faces, points);
1038
1039
1040
1041
1042
1043
        }
    }

    meshedSurfRef surfToWrite(points, faces);

    // Only a few weight types (scalar, vector)
1044
1045
    if (weightFieldName_ != "none")
    {
1046
1047
        if (validField<scalar>(weightFieldName_))
        {
1048
            scalarField weightField
1049
            (
1050
                getFieldValues<scalar>(weightFieldName_, true)
1051
            );
1052

1053
1054
1055
1056
1057
            // Process the fields
            writeAll(Sf, weightField, surfToWrite);
        }
        else if (validField<vector>(weightFieldName_))
        {
1058
            vectorField weightField
1059
            (
1060
                getFieldValues<vector>(weightFieldName_, true)
1061
            );
1062

1063
1064
1065
1066
            // Process the fields
            writeAll(Sf, weightField, surfToWrite);
        }
        else
1067
        {
1068
1069
1070
1071
            FatalErrorInFunction
                << "weightField " << weightFieldName_
                << " not found or an unsupported type"
                << abort(FatalError);
1072
        }
1073
    }
1074
1075
1076
1077
1078
1079
1080
1081
    else
    {
        // Default is a zero-size scalar weight field (ie, weight = 1)
        scalarField weightField;

        // Process the fields
        writeAll(Sf, weightField, surfToWrite);
    }