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 | Copyright (C) 2011-2017 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 (regionType_ != stSampledSurface && obr_.foundObject<sf>(fieldName))
45  {
46  return true;
47  }
48  else if (obr_.foundObject<vf>(fieldName))
49  {
50  return true;
51  }
52 
53  return false;
54 }
55 
56 
57 template<class Type>
60 (
61  const word& fieldName,
62  const bool mustGet,
63  const bool applyOrientation
64 ) const
65 {
68 
69  if (regionType_ != stSampledSurface && obr_.foundObject<sf>(fieldName))
70  {
71  return filterField(obr_.lookupObject<sf>(fieldName), applyOrientation);
72  }
73  else if (obr_.foundObject<vf>(fieldName))
74  {
75  const vf& fld = obr_.lookupObject<vf>(fieldName);
76 
77  if (surfacePtr_.valid())
78  {
79  if (surfacePtr_().interpolate())
80  {
81  const interpolationCellPoint<Type> interp(fld);
82  tmp<Field<Type>> tintFld(surfacePtr_().interpolate(interp));
83  const Field<Type>& intFld = tintFld();
84 
85  // Average
86  const faceList& faces = surfacePtr_().faces();
87  tmp<Field<Type>> tavg
88  (
89  new Field<Type>(faces.size(), Zero)
90  );
91  Field<Type>& avg = tavg.ref();
92 
93  forAll(faces, facei)
94  {
95  const face& f = faces[facei];
96  forAll(f, fp)
97  {
98  avg[facei] += intFld[f[fp]];
99  }
100  avg[facei] /= f.size();
101  }
102 
103  return tavg;
104  }
105  else
106  {
107  return surfacePtr_().sample(fld);
108  }
109  }
110  else
111  {
112  return filterField(fld, applyOrientation);
113  }
114  }
115 
116  if (mustGet)
117  {
119  << "Field " << fieldName << " not found in database"
120  << abort(FatalError);
121  }
122 
123  return tmp<Field<Type>>(new Field<Type>(0));
124 }
125 
126 
127 template<class Type>
130 (
131  const Field<Type>& values,
132  const vectorField& Sf,
133  const scalarField& weightField
134 ) const
135 {
136  Type result = Zero;
137  switch (operation_)
138  {
139  case opSum:
140  {
141  result = sum(values);
142  break;
143  }
144  case opWeightedSum:
145  {
146  if (weightField.size())
147  {
148  result = sum(weightField*values);
149  }
150  else
151  {
152  result = sum(values);
153  }
154  break;
155  }
156  case opSumMag:
157  {
158  result = sum(cmptMag(values));
159  break;
160  }
161  case opSumDirection:
162  {
164  << "Operation " << operationTypeNames_[operation_]
165  << " not available for values of type "
167  << exit(FatalError);
168 
169  result = Zero;
170  break;
171  }
172  case opSumDirectionBalance:
173  {
175  << "Operation " << operationTypeNames_[operation_]
176  << " not available for values of type "
178  << exit(FatalError);
179 
180  result = Zero;
181  break;
182  }
183  case opAverage:
184  {
185  result = sum(values)/values.size();
186  break;
187  }
188  case opWeightedAverage:
189  {
190  if (weightField.size())
191  {
192  result = sum(weightField*values)/sum(weightField);
193  }
194  else
195  {
196  result = sum(values)/values.size();
197  }
198  break;
199  }
200  case opAreaAverage:
201  {
202  const scalarField magSf(mag(Sf));
203 
204  result = sum(magSf*values)/sum(magSf);
205  break;
206  }
207  case opWeightedAreaAverage:
208  {
209  const scalarField magSf(mag(Sf));
210 
211  if (weightField.size())
212  {
213  result = sum(weightField*magSf*values)/sum(magSf*weightField);
214  }
215  else
216  {
217  result = sum(magSf*values)/sum(magSf);
218  }
219  break;
220  }
221  case opAreaIntegrate:
222  {
223  const scalarField magSf(mag(Sf));
224 
225  result = sum(magSf*values);
226  break;
227  }
228  case opWeightedAreaIntegrate:
229  {
230  const scalarField magSf(mag(Sf));
231 
232  if (weightField.size())
233  {
234  result = sum(weightField*magSf*values);
235  }
236  else
237  {
238  result = sum(magSf*values);
239  }
240  break;
241  }
242  case opMin:
243  {
244  result = min(values);
245  break;
246  }
247  case opMax:
248  {
249  result = max(values);
250  break;
251  }
252  case opCoV:
253  {
254  const scalarField magSf(mag(Sf));
255 
256  Type meanValue = sum(values*magSf)/sum(magSf);
257 
258  const label nComp = pTraits<Type>::nComponents;
259 
260  for (direction d=0; d<nComp; ++d)
261  {
262  scalarField vals(values.component(d));
263  scalar mean = component(meanValue, d);
264  scalar& res = setComponent(result, d);
265 
266  res = sqrt(sum(magSf*sqr(vals - mean))/sum(magSf))/mean;
267  }
268 
269  break;
270  }
271  case opAreaNormalAverage:
272  {}
273  case opAreaNormalIntegrate:
274  {}
275  case opNone:
276  {}
277  }
278 
279  return result;
280 }
281 
282 
283 template<class Type>
285 (
286  const Field<Type>& values,
287  const vectorField& Sf,
288  const scalarField& weightField
289 ) const
290 {
291  return processSameTypeValues(values, Sf, weightField);
292 }
293 
294 
295 
296 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
297 
298 template<class Type>
300 (
301  const word& fieldName,
302  const scalarField& weightField,
303  const bool orient
304 )
305 {
306  const bool ok = validField<Type>(fieldName);
307 
308  if (ok)
309  {
310  Field<Type> values(getFieldValues<Type>(fieldName, true, orient));
311 
312  vectorField Sf;
313  if (surfacePtr_.valid())
314  {
315  // Get oriented Sf
316  Sf = surfacePtr_().Sf();
317  }
318  else
319  {
320  // Get oriented Sf
321  Sf = filterField(mesh_.Sf(), true);
322  }
323 
324  // Combine onto master
325  combineFields(values);
326  combineFields(Sf);
327 
328  // Write raw values on surface if specified
329  if (surfaceWriterPtr_.valid())
330  {
331  faceList faces;
333 
334  if (surfacePtr_.valid())
335  {
336  combineSurfaceGeometry(faces, points);
337  }
338  else
339  {
340  combineMeshGeometry(faces, points);
341  }
342 
343  if (Pstream::master())
344  {
345  surfaceWriterPtr_->write
346  (
347  outputDir(),
348  regionTypeNames_[regionType_] + ("_" + regionName_),
349  points,
350  faces,
351  fieldName,
352  values,
353  false
354  );
355  }
356  }
357 
358  if (operation_ != opNone)
359  {
360  // Apply scale factor
361  values *= scaleFactor_;
362 
363  if (Pstream::master())
364  {
365  Type result = processValues(values, Sf, weightField);
366 
367  // Add to result dictionary, over-writing any previous entry
368  resultDict_.add(fieldName, result, true);
369 
370  file() << tab << result;
371 
372  Log << " " << operationTypeNames_[operation_]
373  << "(" << regionName_ << ") of " << fieldName
374  << " = " << result << endl;
375  }
376  }
377  }
378 
379  return ok;
380 }
381 
382 
383 template<class Type>
386 (
388  const bool applyOrientation
389 ) const
390 {
391  tmp<Field<Type>> tvalues(new Field<Type>(faceId_.size()));
392  Field<Type>& values = tvalues.ref();
393 
394  forAll(values, i)
395  {
396  label facei = faceId_[i];
397  label patchi = facePatchId_[i];
398  if (patchi >= 0)
399  {
400  values[i] = field.boundaryField()[patchi][facei];
401  }
402  else
403  {
405  << type() << " " << name() << ": "
406  << regionTypeNames_[regionType_] << "(" << regionName_ << "):"
407  << nl
408  << " Unable to process internal faces for volume field "
409  << field.name() << nl << abort(FatalError);
410  }
411  }
412 
413  if (applyOrientation)
414  {
415  forAll(values, i)
416  {
417  values[i] *= faceSign_[i];
418  }
419  }
420 
421  return tvalues;
422 }
423 
424 
425 template<class Type>
428 (
430  const bool applyOrientation
431 ) const
432 {
433  tmp<Field<Type>> tvalues(new Field<Type>(faceId_.size()));
434  Field<Type>& values = tvalues.ref();
435 
436  forAll(values, i)
437  {
438  label facei = faceId_[i];
439  label patchi = facePatchId_[i];
440  if (patchi >= 0)
441  {
442  values[i] = field.boundaryField()[patchi][facei];
443  }
444  else
445  {
446  values[i] = field[facei];
447  }
448  }
449 
450  if (applyOrientation)
451  {
452  forAll(values, i)
453  {
454  values[i] *= faceSign_[i];
455  }
456  }
457 
458  return tvalues;
459 }
460 
461 
462 // ************************************************************************* //
Foam::surfaceFields.
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
type
Types of root.
Definition: Roots.H:52
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:291
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:261
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:163
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
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:57
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:91
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:657
void cmptMag(FieldField< Field, Type > &cf, const FieldField< Field, Type > &f)
static const char nl
Definition: Ostream.H:262
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.
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.
#define Log
Report write to Foam::Info if the local log switch is true.
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:79
void component(FieldField< Field, typename FieldField< Field, Type >::cmptType > &sf, const FieldField< Field, Type > &f, const direction d)