surfaceFieldValueTemplates.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2019 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "surfaceFieldValue.H"
27 #include "surfaceFields.H"
28 #include "volFields.H"
29 #include "sampledSurface.H"
30 #include "surfaceWriter.H"
31 #include "interpolationCellPoint.H"
32 
33 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
34 
35 template<class Type>
37 (
38  const word& fieldName
39 ) const
40 {
43 
44  if
45  (
46  regionType_ != regionTypes::sampledSurface
47  && obr_.foundObject<sf>(fieldName)
48  )
49  {
50  return true;
51  }
52  else if (obr_.foundObject<vf>(fieldName))
53  {
54  return true;
55  }
56 
57  return false;
58 }
59 
60 
61 template<class Type>
64 (
65  const word& fieldName,
66  const bool mustGet,
67  const bool applyOrientation
68 ) const
69 {
72 
73  if
74  (
75  regionType_ != regionTypes::sampledSurface
76  && obr_.foundObject<sf>(fieldName)
77  )
78  {
79  return filterField(obr_.lookupObject<sf>(fieldName), applyOrientation);
80  }
81  else if (obr_.foundObject<vf>(fieldName))
82  {
83  const vf& fld = obr_.lookupObject<vf>(fieldName);
84 
85  if (surfacePtr_.valid())
86  {
87  if (surfacePtr_().interpolate())
88  {
89  const interpolationCellPoint<Type> interp(fld);
90  tmp<Field<Type>> tintFld(surfacePtr_().interpolate(interp));
91  const Field<Type>& intFld = tintFld();
92 
93  // Average
94  const faceList& faces = surfacePtr_().faces();
95  tmp<Field<Type>> tavg
96  (
97  new Field<Type>(faces.size(), Zero)
98  );
99  Field<Type>& avg = tavg.ref();
100 
101  forAll(faces, facei)
102  {
103  const face& f = faces[facei];
104  forAll(f, fp)
105  {
106  avg[facei] += intFld[f[fp]];
107  }
108  avg[facei] /= f.size();
109  }
110 
111  return tavg;
112  }
113  else
114  {
115  return surfacePtr_().sample(fld);
116  }
117  }
118  else
119  {
120  return filterField(fld, applyOrientation);
121  }
122  }
123 
124  if (mustGet)
125  {
127  << "Field " << fieldName << " not found in database"
128  << abort(FatalError);
129  }
130 
131  return tmp<Field<Type>>(new Field<Type>(0));
132 }
133 
134 
135 template<class Type>
138 (
139  const Field<Type>& values,
140  const vectorField& Sf,
141  const scalarField& weightField
142 ) const
143 {
144  Type result = Zero;
145  switch (operation_)
146  {
147  case operationType::sum:
148  {
149  result = sum(values);
150  break;
151  }
152  case operationType::weightedSum:
153  {
154  if (weightField.size())
155  {
156  result = sum(weightField*values);
157  }
158  else
159  {
160  result = sum(values);
161  }
162  break;
163  }
165  {
166  result = sum(cmptMag(values));
167  break;
168  }
169  case operationType::sumDirection:
170  {
172  << "Operation " << operationTypeNames_[operation_]
173  << " not available for values of type "
175  << exit(FatalError);
176 
177  result = Zero;
178  break;
179  }
180  case operationType::sumDirectionBalance:
181  {
183  << "Operation " << operationTypeNames_[operation_]
184  << " not available for values of type "
186  << exit(FatalError);
187 
188  result = Zero;
189  break;
190  }
192  {
193  result = sum(values)/values.size();
194  break;
195  }
196  case operationType::weightedAverage:
197  {
198  if (weightField.size())
199  {
200  result =
201  sum(weightField*values)
202  /stabilise(sum(weightField), vSmall);
203  }
204  else
205  {
206  result = sum(values)/values.size();
207  }
208  break;
209  }
210  case operationType::areaAverage:
211  {
212  const scalarField magSf(mag(Sf));
213 
214  result = sum(magSf*values)/sum(magSf);
215  break;
216  }
217  case operationType::weightedAreaAverage:
218  {
219  const scalarField magSf(mag(Sf));
220 
221  if (weightField.size())
222  {
223  result =
224  sum(weightField*magSf*values)
225  /stabilise(sum(magSf*weightField), vSmall);
226  }
227  else
228  {
229  result = sum(magSf*values)/sum(magSf);
230  }
231  break;
232  }
233  case operationType::areaIntegrate:
234  {
235  const scalarField magSf(mag(Sf));
236 
237  result = sum(magSf*values);
238  break;
239  }
240  case operationType::weightedAreaIntegrate:
241  {
242  const scalarField magSf(mag(Sf));
243 
244  if (weightField.size())
245  {
246  result = sum(weightField*magSf*values);
247  }
248  else
249  {
250  result = sum(magSf*values);
251  }
252  break;
253  }
254  case operationType::min:
255  {
256  result = min(values);
257  break;
258  }
259  case operationType::max:
260  {
261  result = max(values);
262  break;
263  }
264  case operationType::CoV:
265  {
266  const scalarField magSf(mag(Sf));
267 
268  Type meanValue = sum(values*magSf)/sum(magSf);
269 
270  const label nComp = pTraits<Type>::nComponents;
271 
272  for (direction d=0; d<nComp; ++d)
273  {
274  scalarField vals(values.component(d));
275  scalar mean = component(meanValue, d);
276  scalar& res = setComponent(result, d);
277 
278  res = sqrt(sum(magSf*sqr(vals - mean))/sum(magSf))/mean;
279  }
280 
281  break;
282  }
283  case operationType::areaNormalAverage:
284  {}
285  case operationType::areaNormalIntegrate:
286  {}
287  case operationType::none:
288  {}
289  }
290 
291  return result;
292 }
293 
294 
295 template<class Type>
297 (
298  const Field<Type>& values,
299  const vectorField& Sf,
300  const scalarField& weightField
301 ) const
302 {
303  return processSameTypeValues(values, Sf, weightField);
304 }
305 
306 
307 
308 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
309 
310 template<class Type>
312 (
313  const word& fieldName,
314  const scalarField& weightField,
315  const bool orient
316 )
317 {
318  const bool ok = validField<Type>(fieldName);
319 
320  if (ok)
321  {
322  Field<Type> values(getFieldValues<Type>(fieldName, true, orient));
323 
324  vectorField Sf;
325  if (surfacePtr_.valid())
326  {
327  // Get oriented Sf
328  Sf = surfacePtr_().Sf();
329  }
330  else
331  {
332  // Get oriented Sf
333  Sf = filterField(mesh_.Sf(), true);
334  }
335 
336  // Combine onto master
337  combineFields(values);
338  combineFields(Sf);
339 
340  // Write raw values on surface if specified
341  if (surfaceWriterPtr_.valid())
342  {
343  faceList faces;
345 
346  if (surfacePtr_.valid())
347  {
348  combineSurfaceGeometry(faces, points);
349  }
350  else
351  {
352  combineMeshGeometry(faces, points);
353  }
354 
355  if (Pstream::master())
356  {
357  surfaceWriterPtr_->write
358  (
359  outputDir(),
360  regionTypeNames_[regionType_] + ("_" + regionName_),
361  points,
362  faces,
363  fieldName,
364  values,
365  false
366  );
367  }
368  }
369 
370  if (operation_ != operationType::none)
371  {
372  // Apply scale factor
373  values *= scaleFactor_;
374 
375  if (Pstream::master())
376  {
377  Type result = processValues(values, Sf, weightField);
378 
379  // Add to result dictionary, over-writing any previous entry
380  resultDict_.add(fieldName, result, true);
381 
382  file() << tab << result;
383 
384  Log << " " << operationTypeNames_[operation_]
385  << "(" << regionName_ << ") of " << fieldName
386  << " = " << result << endl;
387  }
388  }
389  }
390 
391  return ok;
392 }
393 
394 
395 template<class Type>
398 (
400  const bool applyOrientation
401 ) const
402 {
403  tmp<Field<Type>> tvalues(new Field<Type>(faceId_.size()));
404  Field<Type>& values = tvalues.ref();
405 
406  forAll(values, i)
407  {
408  label facei = faceId_[i];
409  label patchi = facePatchId_[i];
410  if (patchi >= 0)
411  {
412  values[i] = field.boundaryField()[patchi][facei];
413  }
414  else
415  {
417  << type() << " " << name() << ": "
418  << regionTypeNames_[regionType_] << "(" << regionName_ << "):"
419  << nl
420  << " Unable to process internal faces for volume field "
421  << field.name() << nl << abort(FatalError);
422  }
423  }
424 
425  if (applyOrientation)
426  {
427  forAll(values, i)
428  {
429  values[i] *= faceSign_[i];
430  }
431  }
432 
433  return tvalues;
434 }
435 
436 
437 template<class Type>
440 (
442  const bool applyOrientation
443 ) const
444 {
445  tmp<Field<Type>> tvalues(new Field<Type>(faceId_.size()));
446  Field<Type>& values = tvalues.ref();
447 
448  forAll(values, i)
449  {
450  label facei = faceId_[i];
451  label patchi = facePatchId_[i];
452  if (patchi >= 0)
453  {
454  values[i] = field.boundaryField()[patchi][facei];
455  }
456  else
457  {
458  values[i] = field[facei];
459  }
460  }
461 
462  if (applyOrientation)
463  {
464  forAll(values, i)
465  {
466  values[i] *= faceSign_[i];
467  }
468  }
469 
470  return tvalues;
471 }
472 
473 
474 // ************************************************************************* //
Foam::surfaceFields.
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
const word & name() const
Return name.
Definition: IOobject.H:303
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
static const char tab
Definition: Ostream.H:259
error FatalError
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
const Boundary & boundaryField() const
Return const-reference to the boundary field.
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:174
dimensionedSymmTensor sqr(const dimensionedVector &dv)
uint8_t direction
Definition: direction.H:45
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
Traits class for primitives.
Definition: pTraits.H:50
Generic GeometricField class.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().x()<< ' ';}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().y()<< ' ';}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().z()<< ' ';}gmvFile<< nl;forAll(lagrangianScalarNames, i){ const word &name=lagrangianScalarNames[i];IOField< scalar > fld(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
const pointField & points
Pre-declare SubField and related Field type.
Definition: Field.H:56
A class for handling words, derived from string.
Definition: word.H:59
bool validField(const word &fieldName) const
Return true if the field name is valid.
Type processValues(const Field< Type > &values, const vectorField &Sf, const scalarField &weightField) const
Apply the &#39;operation&#39; to the values. Wrapper around.
static const zero Zero
Definition: zero.H:97
tmp< Field< Type > > getFieldValues(const word &fieldName, const bool mustGet=false, const bool applyOrientation=false) const
Return field values by looking up field name.
errorManip< error > abort(error &err)
Definition: errorManip.H:131
tmp< Field< cmptType > > component(const direction) const
Return a component field of the field.
Definition: Field.C:460
void cmptMag(FieldField< Field, Type > &cf, const FieldField< Field, Type > &f)
static const char nl
Definition: Ostream.H:260
tmp< Field< Type > > filterField(const GeometricField< Type, fvsPatchField, surfaceMesh > &field, const bool applyOrientation) const
Filter a surface field according to faceIds.
Type processSameTypeValues(const Field< Type > &values, const vectorField &Sf, const scalarField &weightField) const
Apply the &#39;operation&#39; to the values. Operation has to.
labelList f(nPoints)
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
volScalarField sf(fieldObject, mesh)
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
static tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > interpolate(const GeometricField< Type, fvPatchField, volMesh > &tvf, const surfaceScalarField &faceFlux, Istream &schemeData)
Interpolate field onto faces using scheme given by Istream.
tmp< DimensionedField< scalar, GeoMesh > > stabilise(const DimensionedField< scalar, GeoMesh > &dsf, const dimensioned< scalar > &ds)
label patchi
Given cell centre values and point (vertex) values decompose into tetrahedra and linear interpolate w...
bool writeValues(const word &fieldName, const scalarField &weightField, const bool orient)
Templated helper function to output field values.
fileType type(const fileName &, const bool checkVariants=true, const bool followLink=true)
Return the file type: directory or file.
Definition: POSIX.C:488
#define Log
Report write to Foam::Info if the local log switch is true.
tmp< GeometricField< Type, fvPatchField, volMesh > > average(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Area-weighted average a surfaceField creating a volField.
Definition: fvcAverage.C:46
dimensioned< scalar > mag(const dimensioned< Type > &)
A class for managing temporary objects.
Definition: PtrList.H:53
label & setComponent(label &l, const direction)
Definition: label.H:86
void component(FieldField< Field, typename FieldField< Field, Type >::cmptType > &sf, const FieldField< Field, Type > &f, const direction d)
dimensioned< scalar > sumMag(const DimensionedField< Type, GeoMesh > &df)