faceAreaIntersect.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) 2011-2021 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 "faceAreaIntersect.H"
27 #include "polygonTriangulate.H"
28 
29 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
30 
31 namespace Foam
32 {
33  template<>
35  {
36  "fan",
37  "mesh"
38  };
39 }
40 
43 
44 Foam::scalar Foam::faceAreaIntersect::tol = 1e-6;
45 
46 // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
47 
48 void Foam::faceAreaIntersect::triSliceWithPlane
49 (
50  const triPoints& tri,
51  const plane& p,
53  label& nTris,
54  const scalar len
55 )
56 {
57  // distance to cutting plane
59 
60  // determine how many of the points are above the cutting plane
61  label nCoPlanar = 0;
62  label nPos = 0;
63  label posI = -1;
64  label negI = -1;
65  label copI = -1;
66  forAll(tri, i)
67  {
68  d[i] = ((tri[i] - p.refPoint()) & p.normal());
69 
70  if (mag(d[i]) < tol*len)
71  {
72  nCoPlanar++;
73  copI = i;
74  d[i] = 0.0;
75  }
76  else
77  {
78  if (d[i] > 0)
79  {
80  nPos++;
81  posI = i;
82  }
83  else
84  {
85  negI = i;
86  }
87  }
88  }
89 
90 
91  // Determine triangle area contribution
92 
93  if
94  (
95  (nPos == 3)
96  || ((nPos == 2) && (nCoPlanar == 1))
97  || ((nPos == 1) && (nCoPlanar == 2))
98  )
99  {
100  /*
101  /\ _____
102  / \ \ / /\
103  /____\ \ / / \
104  __________ ____v____ __/____\__
105 
106  all points above cutting plane
107  - add complete triangle to list
108  */
109 
110  tris[nTris++] = tri;
111  }
112  else if ((nPos == 2) && (nCoPlanar == 0))
113  {
114  /*
115  i1________i2
116  \ /
117  --\----/--
118  \ /
119  \/
120  i0
121 
122  2 points above plane, 1 below
123  - resulting quad above plane split into 2 triangles
124  - forget triangle below plane
125  */
126 
127  // point under the plane
128  label i0 = negI;
129 
130  // indices of remaining points
131  label i1 = d.fcIndex(i0);
132  label i2 = d.fcIndex(i1);
133 
134  // determine the two intersection points
135  point p01 = planeIntersection(d, tri, i0, i1);
136  point p02 = planeIntersection(d, tri, i0, i2);
137 
138  // forget triangle below plane
139  // - decompose quad above plane into 2 triangles and add to list
140  setTriPoints(tri[i1], tri[i2], p02, nTris, tris);
141  setTriPoints(tri[i1], p02, p01, nTris, tris);
142  }
143  else if (nPos == 1)
144  {
145  // point above the plane
146  label i0 = posI;
147 
148  if (nCoPlanar == 0)
149  {
150  /*
151  i0
152  /\
153  / \
154  --/----\--
155  /______\
156  i2 i1
157 
158  1 point above plane, 2 below
159  - keep triangle above intersection plane
160  - forget quad below plane
161  */
162 
163  // indices of remaining points
164  label i1 = d.fcIndex(i0);
165  label i2 = d.fcIndex(i1);
166 
167  // determine the two intersection points
168  point p01 = planeIntersection(d, tri, i1, i0);
169  point p02 = planeIntersection(d, tri, i2, i0);
170 
171  // add triangle above plane to list
172  setTriPoints(tri[i0], p01, p02, nTris, tris);
173  }
174  else
175  {
176  /*
177  i0
178  |\
179  | \
180  __|__\_i2_
181  | /
182  | /
183  |/
184  i1
185 
186  1 point above plane, 1 on plane, 1 below
187  - keep triangle above intersection plane
188  */
189 
190  // point indices
191  label i1 = negI;
192  label i2 = copI;
193 
194  // determine the intersection point
195  point p01 = planeIntersection(d, tri, i1, i0);
196 
197  // add triangle above plane to list - clockwise points
198  if (d.fcIndex(i0) == i1)
199  {
200  setTriPoints(tri[i0], p01, tri[i2], nTris, tris);
201  }
202  else
203  {
204  setTriPoints(tri[i0], tri[i2], p01, nTris, tris);
205  }
206  }
207  }
208  else
209  {
210  /*
211  _________ __________ ___________
212  /\ \ /
213  /\ / \ \ /
214  / \ /____\ \/
215  /____\
216 
217  all points below cutting plane - forget
218  */
219  }
220 }
221 
222 
223 Foam::scalar Foam::faceAreaIntersect::triangleIntersect
224 (
225  const triPoints& src,
226  const triPoints& tgt,
227  const vector& n
228 )
229 {
230  // Work storage
231  FixedList<triPoints, 10> workTris1;
232  label nWorkTris1 = 0;
233 
234  FixedList<triPoints, 10> workTris2;
235  label nWorkTris2 = 0;
236 
237  // cut source triangle with all inwards pointing faces of target triangle
238  // - triangles in workTris1 are inside target triangle
239 
240  scalar t = sqrt(triArea(src));
241 
242  // edge 0
243  {
244  // cut triangle src with plane and put resulting sub-triangles in
245  // workTris1 list
246 
247  scalar s = mag(tgt[1] - tgt[0]);
248  plane pl0(tgt[0], tgt[1], tgt[1] + s*n);
249  triSliceWithPlane(src, pl0, workTris1, nWorkTris1, t);
250  }
251 
252  if (nWorkTris1 == 0)
253  {
254  return 0.0;
255  }
256 
257  // edge1
258  {
259  // cut workTris1 with plane and put resulting sub-triangles in
260  // workTris2 list (re-use tris storage)
261 
262  scalar s = mag(tgt[2] - tgt[1]);
263  plane pl1(tgt[1], tgt[2], tgt[2] + s*n);
264 
265  nWorkTris2 = 0;
266 
267  for (label i = 0; i < nWorkTris1; i++)
268  {
269  triSliceWithPlane(workTris1[i], pl1, workTris2, nWorkTris2, t);
270  }
271 
272  if (nWorkTris2 == 0)
273  {
274  return 0.0;
275  }
276  }
277 
278  // edge2
279  {
280  // cut workTris2 with plane and put resulting sub-triangles in
281  // workTris1 list (re-use workTris1 storage)
282 
283  scalar s = mag(tgt[2] - tgt[0]);
284  plane pl2(tgt[2], tgt[0], tgt[0] + s*n);
285 
286  nWorkTris1 = 0;
287 
288  for (label i = 0; i < nWorkTris2; i++)
289  {
290  triSliceWithPlane(workTris2[i], pl2, workTris1, nWorkTris1, t);
291  }
292 
293  if (nWorkTris1 == 0)
294  {
295  return 0.0;
296  }
297  else
298  {
299  // calculate area of sub-triangles
300  scalar area = 0.0;
301  for (label i = 0; i < nWorkTris1; i++)
302  {
303  area += triArea(workTris1[i]);
304  }
305 
306  return area;
307  }
308  }
309 }
310 
311 
312 // * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * * //
313 
315 (
316  const pointField& pointsA,
317  const pointField& pointsB,
318  const bool reverseB
319 )
320 :
321  pointsA_(pointsA),
322  pointsB_(pointsB),
323  reverseB_(reverseB)
324 {}
325 
326 
327 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
328 
330 (
331  const face& f,
332  const pointField& points,
333  const triangulationMode& triMode,
334  triFaceList& faceTris
335 )
336 {
337  faceTris.resize(f.nTriangles());
338 
339  switch (triMode)
340  {
341  case tmFan:
342  {
343  for (label i = 0; i < f.nTriangles(); ++ i)
344  {
345  faceTris[i][0] = f[0];
346  faceTris[i][1] = f[i + 1];
347  faceTris[i][2] = f[i + 2];
348  }
349 
350  break;
351  }
352  case tmMesh:
353  {
354  polygonTriangulate triEngine;
355  triEngine.triangulate(UIndirectList<point>(points, f));
356  faceTris = triEngine.triPoints(f);
357  }
358  }
359 }
360 
361 
363 (
364  const face& faceA,
365  const face& faceB,
366  const vector& n,
367  const triangulationMode& triMode
368 )
369 {
370  // split faces into triangles
371  triFaceList trisA, trisB;
372  triangulate(faceA, pointsA_, triMode, trisA);
373  triangulate(faceB, pointsB_, triMode, trisB);
374 
375  // intersect triangles
376  scalar totalArea = 0.0;
377  forAll(trisA, tA)
378  {
379  triPoints tpA = getTriPoints(pointsA_, trisA[tA], false);
380 
381 // if (triArea(tpA) > rootVSmall)
382  {
383  forAll(trisB, tB)
384  {
385  triPoints tpB = getTriPoints(pointsB_, trisB[tB], !reverseB_);
386 
387 // if (triArea(tpB) > rootVSmall)
388  {
389  totalArea += triangleIntersect(tpA, tpB, n);
390  }
391  }
392  }
393  }
394 
395  return totalArea;
396 }
397 
398 
399 // ************************************************************************* //
#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
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
A 1D vector of objects of type <T> with a fixed size <Size>.
Definition: FixedList.H:54
dimensionedScalar sqrt(const dimensionedScalar &ds)
Geometric class that creates a 2D plane and can return the intersection point between a line and the ...
Definition: plane.H:60
static const NamedEnum< triangulationMode, 2 > triangulationModeNames_
void resize(const label)
Alias for setSize(const label)
Definition: ListI.H:138
static void triangulate(const face &f, const pointField &points, const triangulationMode &triMode, triFaceList &faceTris)
Triangulate a face using the given triangulation mode.
Initialise the NamedEnum HashTable from the static list of names.
Definition: NamedEnum.H:51
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.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
faceAreaIntersect(const pointField &pointsA, const pointField &pointsB, const bool reverseB=false)
Construct from components.
const point & refPoint() const
Return or return plane base point.
Definition: plane.C:241
scalar calc(const face &faceA, const face &faceB, const vector &n, const triangulationMode &triMode)
Return area of intersection of faceA with faceB.
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
Definition: FixedListI.H:125
const vector & normal() const
Return plane normal.
Definition: plane.C:235
Triangulation of three-dimensional polygons.
A List with indirect addressing.
Definition: fvMatrix.H:106
dimensioned< scalar > mag(const dimensioned< Type > &)
const doubleScalar e
Elementary charge.
Definition: doubleScalar.H:105
label nTriangles() const
Size of the face&#39;s triangulation.
Definition: faceI.H:112
const UList< triFace > & triPoints() const
Get the triangles&#39; points.
Namespace for OpenFOAM.