MPLICfaceI.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 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 "MPLICface.H"
27 
28 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
29 
30 Foam::scalar Foam::MPLICface::alphaPhiU() const
31 {
32  if (flipped_)
33  {
34  return -alphaPhiU(subPointsU_, subPoints_);
35  }
36  else
37  {
38  return alphaPhiU(subPointsU_, subPoints_);
39  }
40 }
41 
42 
43 template<class VectorList, class PointList>
44 inline Foam::scalar Foam::MPLICface::alphaPhiU
45 (
46  const VectorList& pointsU,
47  const PointList& points
48 ) const
49 {
50  const point& baseP = points[0];
51  const vector& baseU = pointsU[0];
52 
53  scalar alphaPhiU = 0;
54  for(label i=1; i < points.size()-1; ++i)
55  {
56  // 1/2 Triangle area
57  const vector area = (points[i] - baseP) ^ (points[i + 1] - baseP);
58 
59  // Calculate face area weighted velocity field
60  // Average is missing 1/3 and area 1/2
61  alphaPhiU += (baseU + pointsU[i] + pointsU[i + 1]) & area;
62  }
63  alphaPhiU /= 6.0;
64 
65  return alphaPhiU;
66 }
67 
68 
69 template<class VectorList, class PointList>
70 inline Foam::scalar Foam::MPLICface::alphaPhiU
71 (
72  const VectorList& pointsU,
73  const PointList& points,
74  const labelList& f
75 ) const
76 {
77  const point& baseP = points[f[0]];
78  const vector& baseU = pointsU[f[0]];
79 
80  scalar alphaPhiU = 0;
81  for(label i=1; i < f.size()-1; ++i)
82  {
83  // 1/2 Triangle area
84  const vector area = (points[f[i]] - baseP) ^ (points[f[i + 1]] - baseP);
85 
86  // Calculate face area weighted velocity field
87  // Average is missing 1/3 and area 1/2
88  alphaPhiU += (baseU + pointsU[f[i]] + pointsU[f[i + 1]]) & area;
89  }
90  alphaPhiU /= 6.0;
91 
92  return alphaPhiU;
93 }
94 
95 
97 {
98  return cutPoints_;
99 }
100 
101 
102 inline const Foam::DynamicList<Foam::point>&
104 {
105  return subPoints_;
106 }
107 
108 
110 {
111  return cutEdges_;
112 }
113 
114 
115 inline const Foam::vector Foam::MPLICface::Sf() const
116 {
117  return face::area(subPoints_);
118 }
119 
120 
122 (
123  const vector& area
124 ) const
125 {
126  // If the face is a triangle, do a direct calculation
127  if (subPoints_.size() == 3)
128  {
129  return (1.0/3.0)*(subPoints_[0] + subPoints_[1] + subPoints_[2]);
130  }
131 
132  // For more complex faces, decompose into triangles ...
133 
134  // Compute an estimate of the centre as the average of the points
135  point pAvg = Zero;
136  forAll(subPoints_, pi)
137  {
138  pAvg += subPoints_[pi];
139  }
140  pAvg /= subPoints_.size();
141 
142  const vector sumAHat = normalised(area);
143 
144  // Compute the area-weighted sum of the triangle centres. Note use
145  // the triangle area projected in the direction of the face normal
146  // as the weight, *not* the triangle area magnitude. Only the
147  // former makes the calculation independent of the initial estimate.
148  scalar sumAn = 0;
149  vector sumAnc = Zero;
150  forAll(subPoints_, pi)
151  {
152  const point& p = subPoints_[pi];
153  const point& pNext = subPoints_[subPoints_.fcIndex(pi)];
154 
155  const vector a = (pNext - p)^(pAvg - p);
156  const vector c = p + pNext + pAvg;
157 
158  const scalar an = a & sumAHat;
159 
160  sumAn += an;
161  sumAnc += an*c;
162  }
163 
164  // Complete calculating centres and areas. If the face is too small
165  // for the sums to be reliably divided then just set the centre to
166  // the initial estimate.
167  if (sumAn > vSmall)
168  {
169  return (1.0/3.0)*sumAnc/sumAn;
170  }
171  else
172  {
173  return pAvg;
174  }
175 }
176 
177 
178 // ************************************************************************* //
const vector Sf() const
Return subface surface area vector.
Definition: MPLICfaceI.H:115
#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 vector Cf(const vector &area) const
Return subface centre.
Definition: MPLICfaceI.H:122
const DynamicList< point > & cutPoints() const
Access to cut points.
Definition: MPLICfaceI.H:96
const DynamicList< label > & cutEdges() const
Access to cut edges.
Definition: MPLICfaceI.H:109
static vector area(const PointField &ps)
Return vector area given face points.
const dimensionedScalar c
Speed of light in a vacuum.
Form normalised(const VectorSpace< Form, Cmpt, Ncmpts > &vs)
Definition: VectorSpaceI.H:413
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects...
Definition: DynamicList.H:56
static const zero Zero
Definition: zero.H:97
const DynamicList< point > & subPoints() const
Access to submerged face points.
Definition: MPLICfaceI.H:103
volScalarField & p
scalar alphaPhiU() const
Calculate and return alphaPhiU.
Definition: MPLICfaceI.H:30