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