surfaceFieldValueTemplates.C 14.2 KB
Newer Older
1
2
3
4
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
5
    \\  /    A nd           | Copyright (C) 2015-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
#include "surfaceFields.H"
30
#include "polySurfaceFields.H"
31
#include "volFields.H"
32
#include "sampledSurface.H"
33
#include "surfaceWriter.H"
34
#include "interpolationCell.H"
35
#include "interpolationCellPoint.H"
36
37
38

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

39
40
41
42
43
44
45
46
47
48
49
50
51
52
template<class WeightType>
inline bool Foam::functionObjects::fieldValues::surfaceFieldValue::canWeight
(
    const Field<WeightType>& weightField
) const
{
    return
    (
        usesWeight()
     && returnReduce(!weightField.empty(), orOp<bool>()) // On some processor
    );
}


53
template<class Type>
54
bool Foam::functionObjects::fieldValues::surfaceFieldValue::validField
55
56
57
(
    const word& fieldName
) const
58
59
60
{
    typedef GeometricField<Type, fvsPatchField, surfaceMesh> sf;
    typedef GeometricField<Type, fvPatchField, volMesh> vf;
61
    typedef DimensionedField<Type, polySurfaceGeoMesh> smt;
62
63
64
65
66

    return
    (
        foundObject<smt>(fieldName)
     || foundObject<vf>(fieldName)
67
     || (withSurfaceFields() && foundObject<sf>(fieldName))
68
    );
69
}
70
71


72
template<class Type>
73
Foam::tmp<Foam::Field<Type>>
74
Foam::functionObjects::fieldValues::surfaceFieldValue::getFieldValues
75
(
76
    const word& fieldName,
77
    const bool mandatory
78
79
80
81
) const
{
    typedef GeometricField<Type, fvsPatchField, surfaceMesh> sf;
    typedef GeometricField<Type, fvPatchField, volMesh> vf;
82
    typedef DimensionedField<Type, polySurfaceGeoMesh> smt;
83

84
    if (foundObject<smt>(fieldName))
85
    {
86
        return lookupObject<smt>(fieldName);
87
    }
88
    else if (withSurfaceFields() && foundObject<sf>(fieldName))
89
    {
90
        return filterField(lookupObject<sf>(fieldName));
91
92
93
94
    }
    else if (foundObject<vf>(fieldName))
    {
        const vf& fld = lookupObject<vf>(fieldName);
95

96
        if (sampledPtr_.valid())
97
        {
98
            if (sampledPtr_().interpolate())
99
100
            {
                const interpolationCellPoint<Type> interp(fld);
101
102

                return sampledPtr_().sample(interp);
103
104
105
            }
            else
            {
106
107
                const interpolationCell<Type> interp(fld);

108
                return sampledPtr_().sample(interp);
109
            }
110
111
112
        }
        else
        {
113
            return filterField(fld);
114
        }
115
116
    }

117
    if (mandatory)
118
    {
119
120
        FatalErrorInFunction
            << "Field " << fieldName << " not found in database"
121
122
123
            << abort(FatalError);
    }

124
    return tmp<Field<Type>>::New();
125
126
127
}


128
template<class Type, class WeightType>
129
130
Type Foam::functionObjects::fieldValues::surfaceFieldValue::
processSameTypeValues
131
(
132
    const Field<Type>& values,
133
    const vectorField& Sf,
134
    const Field<WeightType>& weightField
135
136
) const
{
137
    Type result = Zero;
138
139
    switch (operation_)
    {
140
        case opNone:
141
        {
142
            break;
143
        }
144
145
146
147
148
149
150
151
152
153
154
        case opMin:
        {
            result = gMin(values);
            break;
        }
        case opMax:
        {
            result = gMax(values);
            break;
        }
        case opSumMag:
155
        {
156
            result = gSum(cmptMag(values));
157
158
            break;
        }
159
        case opSum:
160
        case opWeightedSum:
161
        case opAbsWeightedSum:
162
        {
163
            if (canWeight(weightField))
164
            {
165
                tmp<scalarField> weight(weightingFactor(weightField));
166
167
                result = gSum(weight*values);
            }
168
169
170
171
172
            else
            {
                // Unweighted form
                result = gSum(values);
            }
173
174
            break;
        }
175
        case opSumDirection:
176
177
        case opSumDirectionBalance:
        {
178
            FatalErrorInFunction
179
180
181
182
183
184
185
                << "Operation " << operationTypeNames_[operation_]
                << " not available for values of type "
                << pTraits<Type>::typeName
                << exit(FatalError);

            break;
        }
186
        case opAverage:
187
        case opWeightedAverage:
188
        case opAbsWeightedAverage:
189
        {
190
            if (canWeight(weightField))
mattijs's avatar
mattijs committed
191
            {
192
                const scalarField factor(weightingFactor(weightField));
193
194

                result = gSum(factor*values)/(gSum(factor) + ROOTVSMALL);
mattijs's avatar
mattijs committed
195
            }
196
197
198
199
200
201
            else
            {
                // Unweighted form
                const label n = returnReduce(values.size(), sumOp<label>());
                result = gSum(values)/(scalar(n) + ROOTVSMALL);
            }
202
203
            break;
        }
204
        case opAreaAverage:
205
        case opWeightedAreaAverage:
206
        case opAbsWeightedAreaAverage:
207
        {
208
209
210
            if (canWeight(weightField))
            {
                const scalarField factor(weightingFactor(weightField, Sf));
211

212
213
214
215
216
217
218
219
                result = gSum(factor*values)/gSum(factor + ROOTVSMALL);
            }
            else
            {
                // Unweighted form
                const scalarField factor(mag(Sf));
                result = gSum(factor*values)/gSum(factor);
            }
220
221
222
            break;
        }
        case opAreaIntegrate:
223
        case opWeightedAreaIntegrate:
224
        case opAbsWeightedAreaIntegrate:
225
        {
226
227
228
229
230
231
232
233
234
235
236
            if (canWeight(weightField))
            {
                tmp<scalarField> factor(weightingFactor(weightField, Sf));
                result = gSum(factor*values);
            }
            else
            {
                // Unweighted form
                tmp<scalarField> factor(mag(Sf));
                result = gSum(factor*values);
            }
mattijs's avatar
mattijs committed
237
238
            break;
        }
239
240
        case opCoV:
        {
241
            const scalarField magSf(mag(Sf));
242
243
244
            const scalar gSumMagSf = gSum(magSf);

            Type meanValue = gSum(values*magSf)/gSumMagSf;
245

246
            for (direction d=0; d < pTraits<Type>::nComponents; ++d)
247
            {
248
249
                tmp<scalarField> vals(values.component(d));
                const scalar mean = component(meanValue, d);
250
251
                scalar& res = setComponent(result, d);

252
253
254
                res =
                    sqrt(gSum(magSf*sqr(vals - mean))/gSumMagSf)
                   /(mean + ROOTVSMALL);
255
256
257
258
            }

            break;
        }
259

260
261
        case opAreaNormalAverage:
        case opAreaNormalIntegrate:
262
        case opUniformity:
263
264
        {
            // Handled in specializations only
265
            break;
266
        }
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283

        case opWeightedUniformity:
        case opAbsWeightedUniformity:
        {
            if (canWeight(weightField))
            {
                // Change weighting from vector -> scalar and dispatch again
                return processValues<Type, scalar>
                (
                    values,
                    Sf,
                    weightingFactor(weightField)
                );
            }

            break;
        }
284
285
286
287
288
289
    }

    return result;
}


290
template<class Type, class WeightType>
291
Type Foam::functionObjects::fieldValues::surfaceFieldValue::processValues
292
293
294
(
    const Field<Type>& values,
    const vectorField& Sf,
295
    const Field<WeightType>& weightField
296
297
298
299
300
301
) const
{
    return processSameTypeValues(values, Sf, weightField);
}


302
303
304
305
306
template<class WeightType>
Foam::tmp<Foam::scalarField>
Foam::functionObjects::fieldValues::surfaceFieldValue::weightingFactor
(
    const Field<WeightType>& weightField
307
) const
308
{
309
310
    // The scalar form is specialized.
    // For other types always need mag() to generate a scalar field.
311
312
313
    return mag(weightField);
}

314

315
316
template<class WeightType>
Foam::label Foam::functionObjects::fieldValues::surfaceFieldValue::writeAll
317
(
318
    const vectorField& Sf,
319
    const Field<WeightType>& weightField,
320
321
    const pointField& points,
    const faceList& faces
322
)
323
{
324
    label nProcessed = 0;
325

326
    for (const word& fieldName : fields_)
327
    {
328
329
        if
        (
330
331
            writeValues<scalar>(fieldName, Sf, weightField, points, faces)
         || writeValues<vector>(fieldName, Sf, weightField, points, faces)
332
333
         || writeValues<sphericalTensor>
            (
334
                fieldName, Sf, weightField, points, faces
335
            )
336
337
         || writeValues<symmTensor>(fieldName, Sf, weightField, points, faces)
         || writeValues<tensor>(fieldName, Sf, weightField, points, faces)
338
        )
339
        {
340
            ++nProcessed;
341
342
343
        }
        else
        {
344
345
346
347
            WarningInFunction
                << "Requested field " << fieldName
                << " not found in database and not processed"
                << endl;
348
        }
349
350
351
352
353
354
355
356
357
358
359
360
    }

    return nProcessed;
}


template<class Type, class WeightType>
bool Foam::functionObjects::fieldValues::surfaceFieldValue::writeValues
(
    const word& fieldName,
    const vectorField& Sf,
    const Field<WeightType>& weightField,
361
362
    const pointField& points,
    const faceList& faces
363
364
365
366
367
368
)
{
    const bool ok = validField<Type>(fieldName);

    if (ok)
    {
369
        Field<Type> values(getFieldValues<Type>(fieldName, true));
370

371
        // Write raw values on surface if specified
372
        if (surfaceWriterPtr_.valid() && surfaceWriterPtr_->enabled())
373
        {
374
375
376
            Field<Type> allValues(values);
            combineFields(allValues);

377
378
            if (Pstream::master())
            {
379
                surfaceWriterPtr_->open
380
                (
381
382
                    points,
                    faces,
383
384
385
386
387
                    (
                        outputDir()
                      / regionTypeNames_[regionType_] + ("_" + regionName_)
                    ),
                    false  // serial - already merged
388
                );
389
390
391
392

                surfaceWriterPtr_->write(fieldName, allValues);

                surfaceWriterPtr_->clear();
393
            }
394
395
        }

396
        if (operation_ != opNone)
397
        {
398
399
            // Apply scale factor
            values *= scaleFactor_;
400

401
            Type result = processValues(values, Sf, weightField);
402

403
404
405
            switch (postOperation_)
            {
                case postOpNone:
406
                {
407
                    break;
408
                }
409
                case postOpSqrt:
410
411
412
                {
                    // sqrt: component-wise - doesn't change the type
                    for (direction d=0; d < pTraits<Type>::nComponents; ++d)
413
                    {
414
415
                        setComponent(result,  d)
                            = sqrt(mag(component(result,  d)));
416
417
                    }
                    break;
418
                }
419
420
            }

421
            file()<< tab << result;
422

423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
            // Write state/results information
            word prefix, suffix;
            {
                if (postOperation_ != postOpNone)
                {
                    // Adjust result name to include post-operation
                    prefix += postOperationTypeNames_[postOperation_];
                    prefix += '(';
                    suffix += ')';
                }

                prefix += operationTypeNames_[operation_];
                prefix += '(';
                suffix += ')';
            }

            Log << "    " << prefix << regionName_ << suffix
                << " of " << fieldName
441
                <<  " = " << result << endl;
442

443
            // Write state/results information
444
            word resultName = prefix + regionName_ + ',' + fieldName + suffix;
445
            this->setResult(resultName, result);
446
447
448
        }
    }

449
    return ok;
450
451
452
453
}


template<class Type>
454
Foam::tmp<Foam::Field<Type>>
455
Foam::functionObjects::fieldValues::surfaceFieldValue::filterField
456
(
457
    const GeometricField<Type, fvPatchField, volMesh>& field
458
459
) const
{
460
461
    auto tvalues = tmp<Field<Type>>::New(faceId_.size());
    auto& values = tvalues.ref();
462
463
464

    forAll(values, i)
    {
465
466
        const label facei = faceId_[i];
        const label patchi = facePatchId_[i];
467
        if (patchi >= 0)
468
        {
469
            values[i] = field.boundaryField()[patchi][facei];
470
471
472
        }
        else
        {
473
            FatalErrorInFunction
474
                << type() << " " << name() << ": "
475
                << regionTypeNames_[regionType_] << "(" << regionName_ << "):"
476
477
478
479
                << nl
                << "    Unable to process internal faces for volume field "
                << field.name() << nl << abort(FatalError);
        }
480
    }
481

482
    // No need to flip values - all boundary faces point outwards
483
484
485
486
487
488

    return tvalues;
}


template<class Type>
489
Foam::tmp<Foam::Field<Type>>
490
Foam::functionObjects::fieldValues::surfaceFieldValue::filterField
491
(
492
    const GeometricField<Type, fvsPatchField, surfaceMesh>& field
493
494
) const
{
495
496
    auto tvalues = tmp<Field<Type>>::New(faceId_.size());
    auto& values = tvalues.ref();
497
498
499

    forAll(values, i)
    {
500
501
        const label facei = faceId_[i];
        const label patchi = facePatchId_[i];
502
        if (patchi >= 0)
503
        {
504
            values[i] = field.boundaryField()[patchi][facei];
505
506
507
        }
        else
        {
508
            values[i] = field[facei];
509
        }
510
    }
511

512
513
514
515
516
517
518
    if (debug)
    {
        Pout<< "field " << field.name() << " oriented: "
            << field.oriented()() << endl;
    }

    if (field.oriented()())
519
520
521
    {
        forAll(values, i)
        {
522
523
524
525
            if (faceFlip_[i])
            {
                values[i] *= -1;
            }
526
        }
527
528
529
530
531
532
533
    }

    return tvalues;
}


// ************************************************************************* //