cutFace.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 | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2016-2017 DHI
9 Copyright (C) 2018-2019 Johan Roenby
10 Copyright (C) 2019-2020 DLR
11-------------------------------------------------------------------------------
12License
13 This file is part of OpenFOAM.
14
15 OpenFOAM is free software: you can redistribute it and/or modify it
16 under the terms of the GNU General Public License as published by
17 the Free Software Foundation, either version 3 of the License, or
18 (at your option) any later version.
19
20 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
21 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
22 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
23 for more details.
24
25 You should have received a copy of the GNU General Public License
26 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
27
28\*---------------------------------------------------------------------------*/
29
30#include "cutFace.H"
31
32// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
33
35
36// * * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * *
37
39(
40 const label faceI,
41 const scalarList& pointStatus,
42 label firstFullySubmergedPoint,
43 DynamicList<point>& subFacePoints,
44 DynamicList<point>& surfacePoints,
45 label& faceStatus,
46 vector& subFaceCentre,
47 vector& subFaceArea
48)
49{
50 const pointField& points = mesh_.points();
51 const face& f = mesh_.faces()[faceI];
52
53 if (firstFullySubmergedPoint == -1) // is in gasPhase
54 {
55 faceStatus = 1;
56 subFaceCentre = Zero;
57 subFaceArea = Zero;
58 return;
59 }
60
61 // loop face and append the cuts
62 // loop starts at firstFullySubmergedPoint
63 for
64 (
65 label i = firstFullySubmergedPoint;
66 i < firstFullySubmergedPoint + f.size();
67 ++i
68 )
69 {
70 // max two points are appended during one cycle
71 label idx = i % f.size();
72 label nextIdx = (i + 1) % f.size();
73
74 if (pointStatus[idx] > 0) // append fluid point
75 {
76 subFacePoints.append(points[f[idx]]);
77 }
78 else if (pointStatus[idx] == 0) // append cut point
79 {
80 subFacePoints.append(points[f[idx]]);
81 surfacePoints.append(points[f[idx]]);
82 }
83
84 if
85 (
86 (pointStatus[idx] < 0 && pointStatus[nextIdx] > 0)
87 || (pointStatus[idx] > 0 && pointStatus[nextIdx] < 0)
88 ) // cut on edge cut Value is zero
89 {
90 label nextP = f.nextLabel(idx);
91 vector dir = points[nextP] - points[f[idx]];
92 scalar weight =
93 (0.0 - pointStatus[idx]) /
94 (pointStatus[nextIdx] - pointStatus[idx]); // cutValue is zero
95
96 point p = points[f[idx]] + weight * dir;
97
98 subFacePoints.append(p);
99 surfacePoints.append(p);
100 }
101 }
102
103 if (subFacePoints.size() >= 3)
104 {
105 faceStatus = 0;
106 calcSubFaceCentreAndArea(subFacePoints, subFaceCentre, subFaceArea);
107 }
108 else
109 {
110 faceStatus = -1;
111 }
112}
113
114
116(
117 const label faceI,
118 const scalarList& pointStatus,
119 const scalarList& weights,
120 label firstFullySubmergedPoint,
121 DynamicList<point>& subFacePoints,
122 DynamicList<point>& surfacePoints,
123 label& faceStatus,
124 vector& subFaceCentre,
125 vector& subFaceArea
126)
127{
128 const pointField& points = mesh_.points();
129 const face& f = mesh_.faces()[faceI];
130
131 if (firstFullySubmergedPoint == -1) // is in gasPhase
132 {
133 faceStatus = 1;
134 subFaceCentre = Zero;
135 subFaceArea = Zero;
136 return;
137 }
138
139 // loop face and append the cuts
140 // loop starts at firstFullySubmergedPoint
141 for
142 (
143 label i = firstFullySubmergedPoint;
144 i < firstFullySubmergedPoint + f.size();
145 ++i
146 )
147 {
148 // max two points are appended during one cycle
149 label idx = i % f.size();
150 label nextIdx = (i + 1) % f.size();
151
152 if (pointStatus[idx] > 0) // append fluid point
153 {
154 subFacePoints.append(points[f[idx]]);
155 }
156 else if (pointStatus[idx] == 0) // append cut point
157 {
158 subFacePoints.append(points[f[idx]]);
159 surfacePoints.append(points[f[idx]]);
160 }
161
162 if
163 (
164 (pointStatus[idx] < 0 && pointStatus[nextIdx] > 0)
165 || (pointStatus[idx] > 0 && pointStatus[nextIdx] < 0)
166 ) // cut on edge cut Value is zero
167 {
168 label nextP = f.nextLabel(idx);
169 vector dir = points[nextP] - points[f[idx]];
170
171 point p = points[f[idx]] + weights[idx] * dir;
172
173 subFacePoints.append(p);
174 surfacePoints.append(p);
175 }
176 }
177
178 if (subFacePoints.size() >= 3)
179 {
180 faceStatus = 0;
181 calcSubFaceCentreAndArea(subFacePoints, subFaceCentre, subFaceArea);
182 }
183 else
184 {
185 faceStatus = -1;
186 }
187}
188
189
191(
192 const face& f,
193 const pointField& points,
194 const scalarList& pointStatus,
195 label firstFullySubmergedPoint,
196 DynamicList<point>& subFacePoints,
197 DynamicList<point>& surfacePoints,
198 label& faceStatus,
199 vector& subFaceCentre,
200 vector& subFaceArea
201)
202{
203 if (firstFullySubmergedPoint == -1) // in Gas
204 {
205 faceStatus = 1;
206 subFaceCentre = Zero;
207 subFaceArea = Zero;
208 return;
209 }
210
211 // loop face and append the cuts
212 for
213 (
214 label i = firstFullySubmergedPoint;
215 i < firstFullySubmergedPoint + f.size();
216 ++i
217 )
218 {
219 // max two points are appended during one cycle
220 label idx = i % f.size();
221 label nextIdx = (i + 1) % f.size();
222
223 if (pointStatus[idx] > 0) // append fluid point
224 {
225 subFacePoints.append(points[f[idx]]);
226 }
227 else if (pointStatus[idx] == 0) // append cut point
228 {
229 subFacePoints.append(points[f[idx]]);
230 surfacePoints.append(points[f[idx]]);
231 }
232
233 if
234 (
235 (pointStatus[idx] < 0 && pointStatus[nextIdx] > 0)
236 || (pointStatus[idx] > 0 && pointStatus[nextIdx] < 0)
237 )
238 {
239 label nextP = f.nextLabel(idx);
240 vector dir = points[nextP] - points[f[idx]];
241 scalar weight =
242 (0.0 - pointStatus[idx]) /
243 (pointStatus[nextIdx] - pointStatus[idx]);
244
245 point p = points[f[idx]] + weight * dir;
246
247 subFacePoints.append(p);
248 surfacePoints.append(p);
249 }
250 }
251
252 if (subFacePoints.size() >= 3)
253 {
254 faceStatus = 0;
255 calcSubFaceCentreAndArea(subFacePoints, subFaceCentre, subFaceArea);
256 }
257 else
258 {
259 faceStatus = -1;
260 }
261}
262
263
265(
266 DynamicList<point>& subFacePoints,
267 vector& subFaceCentre,
268 vector& subFaceArea
269)
270{
271 const label nPoints = subFacePoints.size();
272
273 // If the face is a triangle, do a direct calculation for efficiency
274 // and to avoid round-off error-related problems
275 if (nPoints == 3)
276 {
277 subFaceCentre =
278 (1.0/3.0)*(subFacePoints[0] + subFacePoints[1] + subFacePoints[2]);
279
280 subFaceArea = 0.5*((subFacePoints[1] - subFacePoints[0]) ^
281 (subFacePoints[2] - subFacePoints[0]));
282 }
283 else if (nPoints > 0)
284 {
285 vector sumN{Zero};
286 scalar sumA{0};
287 vector sumAc{Zero};
288
289 point fCentre = subFacePoints[0];
290 // initial guess of centre as average of subFacePoints
291 for (label pi = 1; pi < nPoints; pi++)
292 {
293 fCentre += subFacePoints[pi];
294 }
295
296 fCentre /= nPoints;
297
298 // loop sub triangles
299 for (label pi = 0; pi < nPoints; pi++)
300 {
301 const point& nextPoint = subFacePoints[(pi + 1) % nPoints];
302
303 vector c = subFacePoints[pi] + nextPoint + fCentre;
304 vector n =
305 (nextPoint - subFacePoints[pi]) ^ (fCentre - subFacePoints[pi]);
306 scalar a = mag(n);
307
308 sumN += n;
309 sumA += a;
310 sumAc += a * c;
311 }
312
313 // This is to deal with zero-area faces. Mark very small faces
314 // to be detected in e.g., processorPolyPatch.
315 if (sumA < ROOTVSMALL)
316 {
317 subFaceCentre = fCentre;
318 subFaceArea = Zero;
319 }
320 else
321 {
322 subFaceCentre = (1.0 / 3.0) * sumAc / sumA;
323 subFaceArea = 0.5 * sumN;
324 }
325 }
326}
327
328
329// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
330
332(
333 const fvMesh& mesh
334)
335:
336 mesh_(mesh)
337{}
338
339
340// ************************************************************************* //
label n
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:72
void append(const T &val)
Copy append an element to the end of this list.
Definition: DynamicListI.H:503
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
Base class for cutting a face, faceI, of an fvMesh, mesh_, at its intersections.
Definition: cutFace.H:60
void calcSubFace(const label faceI, const scalarList &pointStatus, label firstFullySubmergedPoint, DynamicList< point > &subFacePoints, DynamicList< point > &surfacePoints, label &faceStatus, vector &subFaceCentre, vector &subFaceArea)
Definition: cutFace.C:39
void calcSubFaceCentreAndArea(DynamicList< point > &subFacePoints, vector &subFaceCentre, vector &subFaceArea)
Calculates centre and normal of the face.
Definition: cutFace.C:265
static int debug
Definition: cutFace.H:134
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:91
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1108
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1083
volScalarField & p
dynamicFvMesh & mesh
const pointField & points
label nPoints
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
labelList f(nPoints)