surfaceFieldValueTemplates.C
Go to the documentation of this file.
1/*---------------------------------------------------------------------------*\
2 ========= |
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4 \\ / O peration |
5 \\ / A nd | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2011-2017 OpenFOAM Foundation
9 Copyright (C) 2015-2021 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 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#include "surfaceFields.H"
31#include "polySurfaceFields.H"
32#include "volFields.H"
33#include "sampledSurface.H"
34#include "surfaceWriter.H"
35#include "interpolationCell.H"
37
38// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
39
40template<class WeightType>
43(
44 const Field<WeightType>& weightField,
45 const bool useMag /* ignore */
46)
47{
48 // The scalar form is specialized.
49 // Other types: use mag() to generate a scalar field.
50 return mag(weightField);
51}
52
53
54// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
55
56template<class WeightType>
58(
60) const
61{
62 // Non-empty on some processor
63 return returnReduce(!fld.empty(), orOp<bool>());
64}
65
66
67template<class Type>
69(
70 const word& fieldName
71) const
72{
76
77 return
78 (
79 foundObject<smt>(fieldName)
80 || foundObject<vf>(fieldName)
81 || (withSurfaceFields() && foundObject<sf>(fieldName))
82 );
83}
84
85
86template<class Type>
89(
90 const word& fieldName,
91 const bool mandatory
92) const
93{
97
98 if (foundObject<smt>(fieldName))
99 {
100 return lookupObject<smt>(fieldName);
101 }
102 else if (withSurfaceFields() && foundObject<sf>(fieldName))
103 {
104 return filterField(lookupObject<sf>(fieldName));
105 }
106 else if (foundObject<vf>(fieldName))
107 {
108 const vf& fld = lookupObject<vf>(fieldName);
109
110 if (sampledPtr_)
111 {
112 // Could be runtime selectable
113 // auto sampler = interpolation<Type>::New(sampleFaceScheme_, fld);
114
115 // const interpolationCellPoint<Type> interp(fld);
116 const interpolationCell<Type> interp(fld);
117
118 return sampledPtr_->sample(interp);
119 }
120 else
121 {
122 return filterField(fld);
123 }
124 }
125
126 if (mandatory)
127 {
129 << "Field " << fieldName << " not found in database" << nl
130 << abort(FatalError);
131 }
132
133 return tmp<Field<Type>>::New();
134}
135
136
137template<class Type, class WeightType>
140(
141 const Field<Type>& values,
142 const vectorField& Sf,
143 const Field<WeightType>& weightField
144) const
145{
146 Type result = Zero;
147 switch (operation_)
148 {
149 case opNone:
150 {
151 break;
152 }
153 case opMin:
154 {
155 result = gMin(values);
156 break;
157 }
158 case opMax:
159 {
160 result = gMax(values);
161 break;
162 }
163 case opSumMag:
164 {
165 result = gSum(cmptMag(values));
166 break;
167 }
168 case opSum:
169 case opWeightedSum:
170 case opAbsWeightedSum:
171 {
172 if (is_weightedOp() && canWeight(weightField))
173 {
174 tmp<scalarField> weight
175 (
176 weightingFactor(weightField, Sf, is_magOp())
177 );
178
179 result = gSum(weight*values);
180 }
181 else
182 {
183 // Unweighted form
184 result = gSum(values);
185 }
186 break;
187 }
188 case opSumDirection:
189 case opSumDirectionBalance:
190 {
192 << "Operation " << operationTypeNames_[operation_]
193 << " not available for values of type "
195 << exit(FatalError);
196
197 break;
198 }
199 case opAverage:
200 case opWeightedAverage:
201 case opAbsWeightedAverage:
202 {
203 if (is_weightedOp() && canWeight(weightField))
204 {
205 const scalarField factor
206 (
207 weightingFactor(weightField, Sf, is_magOp())
208 );
209
210 result = gSum(factor*values)/(gSum(factor) + ROOTVSMALL);
211 }
212 else
213 {
214 // Unweighted form
215 const label n = returnReduce(values.size(), sumOp<label>());
216
217 result = gSum(values)/(scalar(n) + ROOTVSMALL);
218 }
219 break;
220 }
221 case opAreaAverage:
222 case opWeightedAreaAverage:
223 case opAbsWeightedAreaAverage:
224 {
225 if (is_weightedOp() && canWeight(weightField))
226 {
227 const scalarField factor
228 (
229 areaWeightingFactor(weightField, Sf, is_magOp())
230 );
231
232 result = gSum(factor*values)/gSum(factor + ROOTVSMALL);
233 }
234 else
235 {
236 // Unweighted form
237 const scalarField factor(mag(Sf));
238
239 result = gSum(factor*values)/gSum(factor + ROOTVSMALL);
240 }
241 break;
242 }
243 case opAreaIntegrate:
244 case opWeightedAreaIntegrate:
245 case opAbsWeightedAreaIntegrate:
246 {
247 if (is_weightedOp() && canWeight(weightField))
248 {
249 tmp<scalarField> factor
250 (
251 areaWeightingFactor(weightField, Sf, is_magOp())
252 );
253
254 result = gSum(factor*values);
255 }
256 else
257 {
258 // Unweighted form
259 tmp<scalarField> factor(mag(Sf));
260
261 result = gSum(factor*values);
262 }
263 break;
264 }
265 case opCoV:
266 {
267 const scalarField magSf(mag(Sf));
268 const scalar gSumMagSf = gSum(magSf);
269
270 Type meanValue = gSum(values*magSf)/gSumMagSf;
271
272 for (direction d=0; d < pTraits<Type>::nComponents; ++d)
273 {
274 tmp<scalarField> vals(values.component(d));
275 const scalar mean = component(meanValue, d);
276 scalar& res = setComponent(result, d);
277
278 res =
279 sqrt(gSum(magSf*sqr(vals - mean))/gSumMagSf)
280 /(mean + ROOTVSMALL);
281 }
282
283 break;
284 }
285
286 case opAreaNormalAverage:
287 case opAreaNormalIntegrate:
288 case opUniformity:
289 {
290 // Handled in specializations only
291 break;
292 }
293
294 case opWeightedUniformity:
295 case opAbsWeightedUniformity:
296 {
297 if (is_weightedOp() && canWeight(weightField))
298 {
299 // Change weighting from vector -> scalar and dispatch again
300 return processValues<Type, scalar>
301 (
302 values,
303 Sf,
304 weightingFactor(weightField, is_magOp())
305 );
306 }
307
308 break;
309 }
310 }
311
312 return result;
313}
314
315
316template<class Type, class WeightType>
318(
319 const Field<Type>& values,
320 const vectorField& Sf,
321 const Field<WeightType>& weightField
322) const
323{
324 return processSameTypeValues(values, Sf, weightField);
325}
326
327
328template<class WeightType>
330(
331 const vectorField& Sf,
332 const Field<WeightType>& weightField,
333 const pointField& points,
334 const faceList& faces
335)
336{
337 label nProcessed = 0;
338
339 // If using the surface writer, the points/faces parameters have already
340 // been merged on the master and the writeValues routine will also gather
341 // all data onto the master before calling the writer.
342 // Thus only call the writer on master!
343
344 // Begin writer time step
345 if (Pstream::master() && surfaceWriterPtr_ && surfaceWriterPtr_->enabled())
346 {
347 auto& writer = *surfaceWriterPtr_;
348
349 writer.open
350 (
351 points,
352 faces,
353 (
354 outputDir()
355 / regionTypeNames_[regionType_] + ("_" + regionName_)
356 ),
357 false // serial - already merged
358 );
359
360 writer.beginTime(time_);
361 }
362
363 for (const word& fieldName : fields_)
364 {
365 if
366 (
367 writeValues<scalar>(fieldName, Sf, weightField, points, faces)
368 || writeValues<vector>(fieldName, Sf, weightField, points, faces)
369 || writeValues<sphericalTensor>
370 (
371 fieldName, Sf, weightField, points, faces
372 )
373 || writeValues<symmTensor>(fieldName, Sf, weightField, points, faces)
374 || writeValues<tensor>(fieldName, Sf, weightField, points, faces)
375 )
376 {
377 ++nProcessed;
378 }
379 else
380 {
382 << "Requested field " << fieldName
383 << " not found in database and not processed"
384 << endl;
385 }
386 }
387
388 // Finish writer time step
389 if (Pstream::master() && surfaceWriterPtr_ && surfaceWriterPtr_->enabled())
390 {
391 auto& writer = *surfaceWriterPtr_;
392
393 // Write geometry if no fields were written so that we still
394 // can have something to look at.
395
396 if (!writer.wroteData())
397 {
398 writer.write();
399 }
400
401 writer.endTime();
402 writer.clear();
403 }
404
405 return nProcessed;
406}
407
408
409template<class Type, class WeightType>
411(
412 const word& fieldName,
413 const vectorField& Sf,
414 const Field<WeightType>& weightField,
415 const pointField& points,
416 const faceList& faces
417)
418{
419 const bool ok = validField<Type>(fieldName);
420
421 if (ok)
422 {
423 Field<Type> values(getFieldValues<Type>(fieldName, true));
424
425 // Write raw values on surface if specified
426 if (surfaceWriterPtr_ && surfaceWriterPtr_->enabled())
427 {
428 Field<Type> allValues(values);
429 combineFields(allValues);
430
431 if (Pstream::master())
432 {
434 surfaceWriterPtr_->write(fieldName, allValues);
435
436 // Case-local file name with "<case>" to make relocatable
438 propsDict.add("file", time_.relativePath(outputName, true));
439 this->setProperty(fieldName, propsDict);
440 }
441 }
442
443 if (operation_ != opNone)
444 {
445 // Apply scale factor
446 values *= scaleFactor_;
447
448 Type result = processValues(values, Sf, weightField);
449
450 switch (postOperation_)
451 {
452 case postOpSqrt:
453 {
454 // sqrt: component-wise - does not change the type
455 for (direction d=0; d < pTraits<Type>::nComponents; ++d)
456 {
457 setComponent(result, d)
458 = sqrt(mag(component(result, d)));
459 }
460 break;
461 }
462 default:
463 {
464 break;
465 }
466 }
467
468 // Write state/results information
469 word prefix, suffix;
470 {
471 if (postOperation_ != postOpNone)
472 {
473 // Adjust result name to include post-operation
474 prefix += postOperationTypeNames_[postOperation_];
475 prefix += '(';
476 suffix += ')';
477 }
478
479 prefix += operationTypeNames_[operation_];
480 prefix += '(';
481 suffix += ')';
482 }
483
484 word resultName = prefix + regionName_ + ',' + fieldName + suffix;
485
486 // Write state/results information
487
488 Log << " " << prefix << regionName_ << suffix
489 << " of " << fieldName << " = ";
490
491
492 // Operation or post-operation returns scalar?
493
494 scalar sresult{0};
495
496 bool alwaysScalar(operation_ & typeScalar);
497
498 if (alwaysScalar)
499 {
500 sresult = component(result, 0);
501
502 if (postOperation_ == postOpMag)
503 {
504 sresult = mag(sresult);
505 }
506 }
507 else if (postOperation_ == postOpMag)
508 {
509 sresult = mag(result);
510 alwaysScalar = true;
511 }
512
513
514 if (alwaysScalar)
515 {
516 file()<< tab << sresult;
517
518 Log << sresult << endl;
519
520 this->setResult(resultName, sresult);
521 }
522 else
523 {
524 file()<< tab << result;
525
526 Log << result << endl;
527
528 this->setResult(resultName, result);
529 }
530 }
531 }
532
533 return ok;
534}
535
536
537template<class Type>
540(
542) const
543{
544 const labelList& own = field.mesh().faceOwner();
545 const labelList& nei = field.mesh().faceNeighbour();
546
547 auto tvalues = tmp<Field<Type>>::New(faceId_.size());
548 auto& values = tvalues.ref();
549
550 forAll(values, i)
551 {
552 const label facei = faceId_[i];
553 const label patchi = facePatchId_[i];
554
555 if (patchi >= 0)
556 {
557 // Boundary face - face id is the patch-local face id
558 values[i] = field.boundaryField()[patchi][facei];
559 }
560 else
561 {
562 // Internal face
563 values[i] = 0.5*(field[own[facei]] + field[nei[facei]]);
564 }
565 }
566
567 // No need to flip values - all boundary faces point outwards
568
569 return tvalues;
570}
571
572
573template<class Type>
576(
578) const
579{
580 auto tvalues = tmp<Field<Type>>::New(faceId_.size());
581 auto& values = tvalues.ref();
582
583 forAll(values, i)
584 {
585 const label facei = faceId_[i];
586 const label patchi = facePatchId_[i];
587
588 if (patchi >= 0)
589 {
590 values[i] = field.boundaryField()[patchi][facei];
591 }
592 else
593 {
594 values[i] = field[facei];
595 }
596 }
597
598 if (debug)
599 {
600 Pout<< "field " << field.name() << " oriented: "
601 << field.oriented()() << endl;
602 }
603
604 if (field.oriented()())
605 {
606 forAll(values, i)
607 {
608 if (faceFlip_[i])
609 {
610 values[i] *= -1;
611 }
612 }
613 }
614
615 return tvalues;
616}
617
618
619// ************************************************************************* //
#define Log
Definition: PDRblock.C:35
label n
vtk::internalMeshWriter writer(topoMesh, topoCells, vtk::formatType::INLINE_ASCII, runTime.path()/"blockTopology")
Info<< nl<< "Wrote faMesh in vtk format: "<< writer.output().name()<< nl;}{ vtk::lineWriter writer(aMesh.points(), aMesh.edges(), fileName(aMesh.mesh().time().globalPath()/"finiteArea-edges"));writer.writeGeometry();writer.beginCellData(4);writer.writeProcIDs();{ Field< scalar > fld(faMeshTools::flattenEdgeField(aMesh.magLe(), true))
Field with dimensions and associated with geometry type GeoMesh which is used to size the field and a...
Generic templated field type.
Definition: Field.H:82
Generic GeometricField class.
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
A class for handling file names.
Definition: fileName.H:76
Watches for presence of the named trigger file in the case directory and signals a simulation stop (o...
Definition: abort.H:128
static tmp< scalarField > weightingFactor(const Field< WeightType > &weightField, const bool useMag)
Weighting factor.
tmp< Field< Type > > getFieldValues(const word &fieldName, const bool mandatory=false) const
Return field values by looking up field name.
label writeAll(const vectorField &Sf, const Field< WeightType > &weightField, const pointField &points, const faceList &faces)
Templated helper function to output field values.
bool canWeight(const Field< WeightType > &fld) const
True if field is non-empty on any processor.
Type processSameTypeValues(const Field< Type > &values, const vectorField &Sf, const Field< WeightType > &weightField) const
Apply the 'operation' to the values. Operation must preserve Type.
bool validField(const word &fieldName) const
Return true if the field name is known and a valid type.
bool writeValues(const word &fieldName, const vectorField &Sf, const Field< WeightType > &weightField, const pointField &points, const faceList &faces)
Templated helper function to output field values.
tmp< Field< Type > > filterField(const GeometricField< Type, fvsPatchField, surfaceMesh > &field) const
Filter a surface field according to faceIds.
Type processValues(const Field< Type > &values, const vectorField &Sf, const Field< WeightType > &weightField) const
Apply the 'operation' to the values. Wrapper around.
Computes the magnitude of an input field.
Definition: mag.H:153
Uses the cell value for any location within the cell.
A traits class, which is primarily used for primitives.
Definition: pTraits.H:59
splitCell * master() const
Definition: splitCell.H:113
A class for managing temporary objects.
Definition: tmp.H:65
A class for handling words, derived from Foam::string.
Definition: word.H:68
rDeltaTY field()
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
word outputName("finiteArea-edges.obj")
const pointField & points
#define WarningInFunction
Report a warning using Foam::Warning.
Type gSum(const FieldField< Field, Type > &f)
dimensionedSymmTensor sqr(const dimensionedVector &dv)
void component(FieldField< Field, typename FieldField< Field, Type >::cmptType > &sf, const FieldField< Field, Type > &f, const direction d)
label & setComponent(label &val, const direction) noexcept
Non-const access to integer-type (has no components)
Definition: label.H:123
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
dimensionedScalar sqrt(const dimensionedScalar &ds)
uint8_t direction
Definition: direction.H:56
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh > > &tdf1, const word &name, const dimensionSet &dimensions)
Global function forwards to reuseTmpDimensionedField::New.
error FatalError
Type gMin(const FieldField< Field, Type > &f)
void cmptMag(FieldField< Field, Type > &cf, const FieldField< Field, Type > &f)
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce (copy) and return value.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
Type gMax(const FieldField< Field, Type > &f)
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
constexpr char tab
The tab '\t' character(0x09)
Definition: Ostream.H:52
Fields (face and point) for polySurface.
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
Foam::surfaceFields.
IOdictionary propsDict(dictIO)