MPLICcellI.H
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) 2020-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 "MPLICcell.H"
27 
28 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
29 
30 inline void Foam::MPLICcell::calcAddressing
31 (
32  const MPLICcellStorage& cellInfo
33 )
34 {
35  localEdgeFaces_.setSize(cellInfo.cellEdges().size());
36  localFaceEdges_.setSize(cellInfo.size());
37 
38  // Create edge map, setting edge faces to zero
39  Map<label> edgeMap;
40  forAll(cellInfo.cellEdges(), edgei)
41  {
42  edgeMap.set(cellInfo.cellEdges()[edgei], edgei);
43  localEdgeFaces_[edgei].clear();
44  }
45 
46  // Set size of each face edges to zero
47  forAll(localFaceEdges_, i)
48  {
49  localFaceEdges_[i].clear();
50  }
51 
52  forAll(cellInfo.cellFaces(), facei)
53  {
54  const label faceN = cellInfo.cellFaces()[facei];
55  const labelList& faceEdges = cellInfo.faceEdges()[faceN];
56  forAll(faceEdges, edgei)
57  {
58  const label edg = faceEdges[edgei];
59  const label localEdgeIndex = edgeMap[edg];
60  localEdgeFaces_[localEdgeIndex].append(facei);
61  localFaceEdges_[facei].append(localEdgeIndex);
62  }
63  }
64  addressingCalculated_ = true;
65 }
66 
67 
68 inline bool Foam::MPLICcell::cutStatusCalcSf()
69 {
70  bool cutOrientationDiffers = false;
71  const point pAvg = sum(cutPoints_)/cutPoints_.size();
72  for (label i=0; i<cutPoints_.size(); i+=2)
73  {
74  const vector area = (cutPoints_[i] - pAvg)^(cutPoints_[i+1] - pAvg);
75  cutSf_ += area;
76  if
77  (
78  sign(area.x()*cutSf_.x()) == -1
79  || sign(area.y()*cutSf_.y()) == -1
80  || sign(area.z()*cutSf_.z()) == -1
81  )
82  {
83  cutOrientationDiffers = true;
84  }
85  }
86  cutSf_ *= 0.5;
87 
88  return cutOrientationDiffers;
89 }
90 
91 
92 inline Foam::vector Foam::MPLICcell::calcCutSf() const
93 {
94  vector cutArea = Zero;
95  const point pAvg = sum(cutPoints_)/cutPoints_.size();
96  for (label i=0; i<cutPoints_.size(); i+=2)
97  {
98  cutArea += (cutPoints_[i] - pAvg)^(cutPoints_[i+1] - pAvg);
99  }
100  cutArea *= 0.5;
101 
102  return cutArea;
103 }
104 
105 
106 inline Foam::vector Foam::MPLICcell::calcCutCf(const vector& cutSf) const
107 {
108  const vector sumAHat = normalised(cutSf);
109  const point pAvg = sum(cutPoints_)/cutPoints_.size();
110 
111  scalar sumAn = 0;
112  vector sumAnc = Zero;
113 
114  for (label i=0; i < cutPoints_.size(); i+=2)
115  {
116  const vector a = (cutPoints_[i] - pAvg) ^ (cutPoints_[i+1] - pAvg);
117  const vector c = cutPoints_[i] + cutPoints_[i+1] + pAvg;
118  const scalar an = a & sumAHat;
119  sumAn += an;
120  sumAnc += an*c;
121  }
122 
123  if (sumAn > vSmall)
124  {
125  return (1.0/3.0)*sumAnc/sumAn;
126  }
127  else
128  {
129  return pAvg;
130  }
131 }
132 
133 
134 inline void Foam::MPLICcell::appendSfCf
135 (
136  const vector& Sf,
137  const vector& Cf,
138  const scalar magSf,
139  const bool own
140 )
141 {
142  if (magSf > 0)
143  {
144  if (own)
145  {
146  subFaceAreas_.append(Sf);
147  }
148  else
149  {
150  subFaceAreas_.append(-Sf);
151  }
152  subFaceCentres_.append(Cf);
153  }
154 }
155 
156 
157 inline void Foam::MPLICcell::phiU
158 (
159  const pointField& points,
160  const faceList& faces,
161  const labelList& cFaces,
162  const vectorField& pointsU
163 )
164 {
165  const label nFaces = cFaces.size();
166 
167  // Set size and initialise alphaPhiU
168  alphaPhiU_.setSize(nFaces);
169  alphaPhiU_ = 0;
170 
171  // Set size and initialise phiU
172  phiU_.setSize(nFaces);
173 
174  // Reconstruct fluxes
175  forAll(cFaces, facei)
176  {
177  phiU_[facei] =
178  MPLICface().alphaPhiU(pointsU, points, faces[cFaces[facei]]);
179  }
180 }
181 
182 
183 inline void Foam::MPLICcell::resetFaceFields(const label nFaces)
184 {
185  alphaf_.setSize(nFaces);
186  alphaf_ = 0;
187  if (unweighted_)
188  {
189  subFaceMagSf_.setSize(nFaces);
190  subFaceMagSf_ = 0;
191  }
192  else
193  {
194  alphaPhiU_.setSize(nFaces);
195  alphaPhiU_ = 0;
196  }
197 }
198 
199 
200 inline void Foam::MPLICcell::calcAlphaf
201 (
202  const UIndirectList<scalar>& magSfs
203 )
204 {
205  const label nFaces = magSfs.size();
206  alphaf_.setSize(nFaces);
207  forAll(alphaf_, facei)
208  {
209  if (magSfs[facei] > vSmall)
210  {
211  alphaf_[facei] = subFaceMagSf_[facei]/magSfs[facei];
212 
213  // Bound between 0 and 1 (it is always > 0)
214  alphaf_[facei] = (alphaf_[facei] > 1) ? 1 : alphaf_[facei];
215  }
216  else
217  {
218  alphaf_[facei] = 0;
219  }
220  }
221 }
222 
223 
224 inline void Foam::MPLICcell::calcAlphaUf()
225 {
226  const label nFaces = alphaPhiU_.size();
227  alphaf_.setSize(nFaces);
228  forAll(alphaf_, facei)
229  {
230  if (mag(phiU_[facei]) > vSmall)
231  {
232  alphaf_[facei] = alphaPhiU_[facei]/phiU_[facei];
233 
234  // Bound between 0 and 1
235  alphaf_[facei] = (alphaf_[facei] > 1) ? 1 : alphaf_[facei];
236  alphaf_[facei] = (alphaf_[facei] < 0) ? 0 : alphaf_[facei];
237  }
238  else
239  {
240  alphaf_[facei] = 0;
241  }
242  }
243 }
244 
245 
246 inline void Foam::MPLICcell::clearOneCut()
247 {
248  cutPoints_.clear();
249  cutEdges_.clear();
250  subFaceAreas_.clear();
251  subFaceCentres_.clear();
252 }
253 
254 
255 inline void Foam::MPLICcell::clear()
256 {
257  clearOneCut();
258  alphaPhiU_.clear();
259  subFaceMagSf_.clear();
260  cutAlpha_ = -1;
261  cutNormal_ = Zero;
262  cutSf_ = Zero;
263  subCellVolume_ = 0;
264 }
265 
266 
267 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
268 
270 {
271  return alphaf_;
272 }
273 
274 
276 {
277  return cutNormal_;
278 }
279 
280 
281 inline Foam::scalar Foam::MPLICcell::subCellVolume() const
282 {
283  return subCellVolume_;
284 }
285 
286 
287 inline Foam::scalar Foam::MPLICcell::cutAlpha() const
288 {
289  return cutAlpha_;
290 }
291 
292 
293 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:78
scalar cutAlpha() const
Return volume fraction corresponding to the cut.
Definition: MPLICcellI.H:287
const DynamicList< scalar > & alphaf() const
Return face volume fraction.
Definition: MPLICcellI.H:269
const vector & cutNormal() const
Return cut normal.
Definition: MPLICcellI.H:275
scalar subCellVolume() const
Return sub-cell volume.
Definition: MPLICcellI.H:281
const pointField & points
const dimensionedScalar c
Speed of light in a vacuum.
static const zero Zero
Definition: zero.H:97
List< label > labelList
A List of labels.
Definition: labelList.H:56
dimensionedScalar sign(const dimensionedScalar &ds)
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
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
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
dimensionSet normalised(const dimensionSet &)
Definition: dimensionSet.C:510
dimensioned< scalar > mag(const dimensioned< Type > &)
Field< vector > vectorField
Specialisation of Field<T> for vector.
List< face > faceList
Definition: faceListFwd.H:41