cutPolyIntegralTemplates.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) 2018-2022 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 "cutPolyIntegral.H"
27 
28 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
29 
30 namespace Foam
31 {
32 namespace cutPoly
33 {
34 
35 template<class Type>
37 {
38  Type result = pTraits<Type>::zero;
39  forAll(l, i)
40  {
41  result += l[i];
42  }
43  return result/l.size();
44 }
45 
46 template<class Type>
47 Type average(const List<Type>& xs, const labelUList& is)
48 {
49  return average(UIndirectList<Type>(xs, is));
50 }
51 
52 template<class Container>
53 typename Container::value_type iterableAverage(const Container& xs)
54 {
55  label n = 0;
56  typename Container::value_type nResult =
58 
59  forAllConstIter(typename Container, xs, iter)
60  {
61  ++ n;
62  nResult += *iter;
63  }
64 
65  return nResult/n;
66 }
67 
68 } // End namespace cutPoly
69 } // End namespace Foam
70 
71 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
72 
73 namespace Foam
74 {
75 namespace cutPoly
76 {
77 
78 template<class Op, class Tuple, class Int, Int ... Is>
79 auto tupleOp
80 (
81  const Tuple& tuple,
82  const std::integer_sequence<Int, Is ...>&,
83  const Op& op
84 )
85 {
86  return std::make_tuple(op(std::get<Is>(tuple)) ...);
87 }
88 
89 
90 template<class Op, class ... Types>
91 auto tupleOp(const std::tuple<Types...>& tuple, const Op& op)
92 {
93  return tupleOp(tuple, std::make_index_sequence<sizeof ... (Types)>(), op);
94 }
95 
96 
97 struct OpBegin
98 {
99  template<class Type>
100  auto operator()(const Type& x) const
101  {
102  return x.begin();
103  }
104 };
105 
106 
108 {
109  template<class Type>
110  auto operator()(const Type& x) const
111  {
112  return *x;
113  }
114 };
115 
116 
117 struct OpNext
118 {
119  template<class Type>
120  auto operator()(const Type& x) const
121  {
122  return x.next();
123  }
124 };
125 
126 
127 template<class ScaleType>
128 struct OpScaled
129 {
130  const ScaleType s_;
131 
132  OpScaled(const ScaleType& s)
133  :
134  s_(s)
135  {}
136 
137  template<class Type>
138  auto operator()(const Type& x) const
139  {
140  return s_*x;
141  }
142 };
143 
144 
145 template<class Op, class Tuple, class Int, Int ... Is>
147 (
148  Tuple& tuple,
149  const std::integer_sequence<Int, Is ...>&,
150  const Op& op
151 )
152 {
153  (void)std::initializer_list<nil>
154  {(
155  op(std::get<Is>(tuple)),
156  nil()
157  ) ... };
158 }
159 
160 
161 template<class Op, class ... Types>
162 void tupleInPlaceOp(std::tuple<Types...>& tuple, const Op& op)
163 {
164  tupleInPlaceOp(tuple, std::make_index_sequence<sizeof ... (Types)>(), op);
165 }
166 
167 
169 {
170  template<class Type>
171  void operator()(Type& x) const
172  {
173  ++ x;
174  }
175 };
176 
177 
178 template<class BinaryOp, class Tuple, class Int, Int ... Is>
180 (
181  const Tuple& tupleA,
182  const Tuple& tupleB,
183  const std::integer_sequence<Int, Is ...>&,
184  const BinaryOp& bop
185 )
186 {
187  return std::make_tuple(bop(std::get<Is>(tupleA), std::get<Is>(tupleB)) ...);
188 }
189 
190 
191 template<class BinaryOp, class ... TypesA, class ... TypesB>
193 (
194  const std::tuple<TypesA ...>& tupleA,
195  const std::tuple<TypesB ...>& tupleB,
196  const BinaryOp& bop
197 )
198 {
199  return tupleBinaryOp
200  (
201  tupleA,
202  tupleB,
203  std::make_index_sequence<sizeof ... (TypesA)>(),
204  bop
205  );
206 }
207 
208 
210 {
211  template<class Type>
212  auto operator()(const Type& a, const Type& b) const
213  {
214  return a + b;
215  }
216 };
217 
218 } // End namespace cutPoly
219 } // End namespace Foam
220 
221 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
222 
223 template<template<class> class FaceValues, class ... Types>
224 std::tuple
225 <
226  Foam::vector,
228 >
230 (
231  const FaceValues<point>& fPs,
232  const point& fPAvg,
233  const std::tuple<FaceValues<Types> ...>& fPsis,
234  const std::tuple<Types ...>& fPsiAvg
235 )
236 {
237  vector fCutsArea = Zero;
238  auto fCutsAreaPsi =
239  std::make_tuple(typename outerProduct<vector, Types>::type(Zero) ...);
240 
241  typename FaceValues<point>::const_iterator fPIter(fPs.begin());
242  auto fPsiIter = tupleOp(fPsis, OpBegin());
243 
244  for(; fPIter != fPs.end(); ++ fPIter)
245  {
246  const point p0 = *fPIter;
247  const point p1 = fPIter.next();
248  auto psi0 = tupleOp(fPsiIter, OpDereference());
249  auto psi1 = tupleOp(fPsiIter, OpNext());
250 
251  const vector a = ((p1 - p0)^(fPAvg - p0))/2;
252  auto v =
254  (
255  psi0,
257  (
258  psi1,
259  fPsiAvg,
260  BinaryOpAdd()
261  ),
262  BinaryOpAdd()
263  );
264 
265  auto av = tupleOp(v, OpScaled<vector>(a/3));
266 
267  fCutsArea += a;
268  fCutsAreaPsi = tupleBinaryOp(fCutsAreaPsi, av, BinaryOpAdd());
269 
270  tupleInPlaceOp(fPsiIter, InPlaceOpAdvance());
271  }
272 
273  return std::tuple_cat(std::tuple<vector>(fCutsArea), fCutsAreaPsi);
274 }
275 
276 
277 template<class Type>
279 <
280  Foam::vector,
282 >
284 (
285  const face& f,
286  const vector& fArea,
287  const Type& fPsi,
288  const List<labelPair>& fCuts,
289  const pointField& ps,
290  const Field<Type>& pPsis,
291  const scalarField& pAlphas,
292  const scalar isoAlpha,
293  const bool below
294 )
295 {
296  typedef typename outerProduct<vector, Type>::type IntegralType;
297 
298  // If there are no cuts return either the entire face or zero, depending on
299  // which side of the iso-surface the face is
300  if (fCuts.size() == 0)
301  {
302  if ((pAlphas[f[0]] < isoAlpha) == below)
303  {
304  return Tuple2<vector, IntegralType>(fArea, fArea*fPsi);
305  }
306  else
307  {
308  return Tuple2<vector, IntegralType>(Zero, Zero);
309  }
310  }
311 
312  auto result =
314  (
315  cutPoly::FaceCutValues<vector>
316  (
317  f,
318  fCuts,
319  ps,
320  pAlphas,
321  isoAlpha,
322  below
323  ),
324  average(ps, f),
325  std::make_tuple
326  (
327  cutPoly::FaceCutValues<Type>
328  (
329  f,
330  fCuts,
331  pPsis,
332  pAlphas,
333  isoAlpha,
334  below
335  )
336  ),
337  std::make_tuple(average(pPsis, f))
338  );
339 
340  return
341  Tuple2<vector, IntegralType>
342  (
343  std::get<0>(result),
344  std::get<1>(result)
345  );
346 }
347 
348 
349 template<class Type>
351 (
352  const cell& c,
353  const cellEdgeAddressing& cAddr,
354  const scalar cVolume,
355  const Type& cPsi,
356  const labelListList& cCuts,
357  const faceUList& fs,
358  const vectorField& fAreas,
359  const vectorField& fCentres,
360  const vectorField& fPsis,
361  const vectorField& fCutAreas,
362  const vectorField& fCutPsis,
363  const pointField& ps,
364  const Field<Type>& pPsis,
365  const scalarField& pAlphas,
366  const scalar isoAlpha,
367  const bool below
368 )
369 {
370  // If there are no cuts return either the entire cell or zero, depending on
371  // which side of the iso-surface the cell is
372  if (cCuts.size() == 0)
373  {
374  if ((pAlphas[fs[c[0]][0]] < isoAlpha) == below)
375  {
376  return Tuple2<scalar, Type>(cVolume, cVolume*cPsi);
377  }
378  else
379  {
380  return Tuple2<scalar, Type>(0, Zero);
381  }
382  }
383 
384  // Averages
385  const point cPAvg = average(fCentres, c);
386  const Type cPsiAvg = average(fPsis, c);
387 
388  // Initialise totals
389  scalar cCutsVolume = 0;
390  Type cCutsVolumePsi = Zero;
391 
392  // Face contributions
393  forAll(c, cfi)
394  {
395  // We use the un-cut face's centroid as the base of the pyramid formed
396  // by the cut face. This is potentially less exact, but it is
397  // consistent with the un-cut cell volume calculation. This means if
398  // you run this function with both values of "below" then the result
399  // will exactly sum to the volume of the overall cell. If we used the
400  // cut-face's centroid this would not be the case.
401  const vector& fBaseP = fCentres[c[cfi]];
402 
403  const scalar pyrVolume =
404  (cAddr.cOwns()[cfi] ? +1 : -1)
405  *(fCutAreas[c[cfi]] & (fBaseP - cPAvg))/3;
406 
407  cCutsVolume += pyrVolume;
408  cCutsVolumePsi += pyrVolume*(3*fCutPsis[c[cfi]] + cPsiAvg)/4;
409  }
410 
411  // Cut contributions
412  {
413  const cutPoly::CellCutValues<point> cPValues
414  (
415  c,
416  cAddr,
417  cCuts,
418  fs,
419  ps,
420  pAlphas,
421  isoAlpha
422  );
423 
424  const cutPoly::CellCutValues<Type> cPsiValues
425  (
426  c,
427  cAddr,
428  cCuts,
429  fs,
430  pPsis,
431  pAlphas,
432  isoAlpha
433  );
434 
435  // This method is more exact, as it uses the true centroid of the cell
436  // cut. However, to obtain that centroid we have to divide by the area
437  // magnitude, so this can't be generalised to types that do not support
438  // division.
439  /*
440  auto fSumPPsis =
441  faceAreaIntegral
442  (
443  cPValues,
444  iterableAverage(cPValues),
445  std::make_tuple
446  (
447  cPValues,
448  cPsiValues
449  ),
450  std::make_tuple
451  (
452  iterableAverage(cPValues),
453  iterableAverage(cPsiValues)
454  )
455  );
456 
457  const vector& fCutArea = std::get<0>(fSumPPsis);
458  const scalar fMagSqrCutArea = magSqr(fCutArea);
459 
460  const point fCutCentre =
461  fMagSqrCutArea > vSmall
462  ? (fCutArea & std::get<1>(fSumPPsis))/fMagSqrCutArea
463  : cPAvg;
464  const Type fCutPsi =
465  fMagSqrCutArea > vSmall
466  ? (fCutArea & std::get<2>(fSumPPsis))/fMagSqrCutArea
467  : cPsiAvg;
468  */
469 
470  // This method is more approximate, as it uses point averages rather
471  // than centroids. This does not involve division, though, so this
472  // could be used with types like polynomials that only support addition
473  // and multiplication.
474  auto fSumPPsis =
476  (
477  cPValues,
478  iterableAverage(cPValues),
479  std::make_tuple(),
480  std::make_tuple()
481  );
482 
483  const vector& fCutArea = std::get<0>(fSumPPsis);
484  const point fCutCentre = iterableAverage(cPValues);
485  const Type fCutPsi = iterableAverage(cPsiValues);
486 
487  const scalar pyrVolume =
488  (below ? +1 : -1)*(fCutArea & (fCutCentre - cPAvg))/3;
489 
490  cCutsVolume += pyrVolume;
491  cCutsVolumePsi += pyrVolume*(3*fCutPsi + cPsiAvg)/4;
492  }
493 
494  return Tuple2<scalar, Type>(cCutsVolume, cCutsVolumePsi);
495 }
496 
497 
498 // ************************************************************************* //
label n
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
#define forAllConstIter(Container, container, iter)
Iterate across all elements in the container object of type.
Definition: UList.H:477
A 2-tuple for storing two objects of different types.
Definition: Tuple2.H:66
A List with indirect addressing.
Definition: UIndirectList.H:60
label size() const
Return the number of elements in the list.
A zero-sized class without any storage. Used, for example, in HashSet.
Definition: nil.H:59
typeOfRank< typename pTraits< arg1 >::cmptType, direction(pTraits< arg1 >::rank)+direction(pTraits< arg2 >::rank) >::type type
Definition: products.H:90
Traits class for primitives.
Definition: pTraits.H:53
volScalarField scalarField(fieldObject, mesh)
volVectorField vectorField(fieldObject, mesh)
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){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.name(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
volScalarField & b
Definition: createFields.H:25
const dimensionedScalar c
Speed of light in a vacuum.
auto tupleOp(const std::tuple< Types... > &tuple, const Op &op)
Type average(const UIndirectList< Type > &l)
std::tuple< vector, typename outerProduct< vector, Types >::type ... > faceAreaIntegral(const FaceValues< point > &fPs, const point &fPAvg, const std::tuple< FaceValues< Types > ... > &fPsis, const std::tuple< Types ... > &fPsiAvg)
Compute the face-area and face-area-integrals of the given properties over.
Container::value_type iterableAverage(const Container &xs)
Tuple2< scalar, Type > cellCutVolumeIntegral(const cell &c, const cellEdgeAddressing &cAddr, const scalar cVolume, const Type &cPsi, const labelListList &cCuts, const faceUList &fs, const vectorField &fAreas, const vectorField &fCentres, const vectorField &fPsis, const vectorField &fCutAreas, const vectorField &fCutPsis, const pointField &ps, const Field< Type > &pPsis, const scalarField &pAlphas, const scalar isoAlpha, const bool below)
Compute the cell-volume and cell-volume-integral of the given property over.
auto tupleBinaryOp(const Tuple &tupleA, const Tuple &tupleB, const std::integer_sequence< Int, Is ... > &, const BinaryOp &bop)
void tupleInPlaceOp(Tuple &tuple, const std::integer_sequence< Int, Is ... > &, const Op &op)
void tupleInPlaceOp(std::tuple< Types... > &tuple, const Op &op)
Type average(const List< Type > &xs, const labelUList &is)
auto tupleOp(const Tuple &tuple, const std::integer_sequence< Int, Is ... > &, const Op &op)
auto tupleBinaryOp(const std::tuple< TypesA ... > &tupleA, const std::tuple< TypesB ... > &tupleB, const BinaryOp &bop)
Tuple2< vector, typename outerProduct< vector, Type >::type > faceCutAreaIntegral(const face &f, const vector &fArea, const Type &fPsi, const List< labelPair > &fCuts, const pointField &ps, const Field< Type > &pPsis, const scalarField &pAlphas, const scalar isoAlpha, const bool below)
Compute the face-area and face-area-integral of the given property over.
Namespace for OpenFOAM.
static const zero Zero
Definition: zero.H:97
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
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
vector point
Point is a vector.
Definition: point.H:41
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
fileType type(const fileName &, const bool checkVariants=true, const bool followLink=true)
Return the file type: directory or file.
Definition: POSIX.C:488
UList< face > faceUList
Definition: faceListFwd.H:39
#define Op(opName, op)
Definition: ops.H:100
labelList f(nPoints)
auto operator()(const Type &a, const Type &b) const
auto operator()(const Type &x) const
auto operator()(const Type &x) const
auto operator()(const Type &x) const
auto operator()(const Type &x) const