triangle.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) 2015 OpenFOAM Foundation
9 -------------------------------------------------------------------------------
10 License
11  This file is part of OpenFOAM.
12 
13  OpenFOAM is free software: you can redistribute it and/or modify it
14  under the terms of the GNU General Public License as published by
15  the Free Software Foundation, either version 3 of the License, or
16  (at your option) any later version.
17 
18  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21  for more details.
22 
23  You should have received a copy of the GNU General Public License
24  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25 
26 \*---------------------------------------------------------------------------*/
27 
28 #include "triangle.H"
29 #include "triPoints.H"
30 #include "plane.H"
31 
32 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
33 
34 template<class Point, class PointRef>
35 template<class AboveOp, class BelowOp>
37 (
38  const plane& pln,
39  const triPoints& tri,
40  AboveOp& aboveOp,
41  BelowOp& belowOp
42 )
43 {
44  // distance to cutting plane
45  FixedList<scalar, 3> d;
46 
47  // determine how many of the points are above the cutting plane
48  label nPos = 0;
49  label posI = -1;
50  label negI = -1;
51  forAll(tri, i)
52  {
53  d[i] = pln.signedDistance(tri[i]);
54 
55  if (d[i] > 0)
56  {
57  nPos++;
58  posI = i;
59  }
60  else
61  {
62  negI = i;
63  }
64  }
65 
66  if (nPos == 3)
67  {
68  aboveOp(tri);
69  }
70  else if (nPos == 2)
71  {
72  // point under the plane
73  label i0 = negI;
74 
75  // indices of remaining points
76  label i1 = d.fcIndex(i0);
77  label i2 = d.fcIndex(i1);
78 
79  // determine the two intersection points
80  point p01 = planeIntersection(d, tri, i0, i1);
81  point p02 = planeIntersection(d, tri, i0, i2);
82 
83  aboveOp(triPoints(tri[i1], tri[i2], p02));
84  aboveOp(triPoints(tri[i1], p02, p01));
85  belowOp(triPoints(tri[i0], p01, p02));
86  }
87  else if (nPos == 1)
88  {
89  // point above the plane
90  label i0 = posI;
91 
92  // indices of remaining points
93  label i1 = d.fcIndex(i0);
94  label i2 = d.fcIndex(i1);
95 
96  // determine the two intersection points
97  point p01 = planeIntersection(d, tri, i0, i1);
98  point p02 = planeIntersection(d, tri, i0, i2);
99 
100  belowOp(triPoints(tri[i1], tri[i2], p02));
101  belowOp(triPoints(tri[i1], p02, p01));
102  aboveOp(triPoints(tri[i0], p01, p02));
103  }
104  else
105  {
106  // All below
107  belowOp(tri);
108  }
109 }
110 
111 
112 template<class Point, class PointRef>
113 template<class AboveOp, class BelowOp>
115 (
116  const plane& pl,
117  AboveOp& aboveOp,
118  BelowOp& belowOp
119 ) const
120 {
121  triSliceWithPlane(pl, triPoints(a_, b_, c_), aboveOp, belowOp);
122 }
123 
124 
125 template<class Point, class PointRef>
126 template<class InsideOp, class OutsideOp>
128 (
129  const vector& n,
130  const triangle<Point, PointRef>& tgt,
131  InsideOp& insideOp,
132  OutsideOp& outsideOp
133 ) const
134 {
135  // There are two possibilities with this algorithm - we either cut
136  // the outside triangles with all the edges or not (and keep them
137  // as disconnected triangles). In the first case
138  // we cannot do any evaluation short cut so we've chosen not to re-cut
139  // the outside triangles.
140 
141 
142  triIntersectionList insideTrisA;
143  label nInsideA = 0;
144  storeOp insideOpA(insideTrisA, nInsideA);
145 
146  triIntersectionList outsideTrisA;
147  label nOutsideA = 0;
148  storeOp outsideOpA(outsideTrisA, nOutsideA);
149 
150 
151  const triPoints thisTri(a_, b_, c_);
152 
153 
154  // Cut original triangle with tgt edge 0.
155  // From *this to insideTrisA, outsideTrisA.
156  {
157  scalar s = Foam::mag(tgt.b() - tgt.a());
158  const plane pl0(tgt.a(), tgt.b(), tgt.b() + s*n);
159  triSliceWithPlane(pl0, thisTri, insideOpA, outsideOpA);
160  }
161 
162  // Shortcut if nothing cut
163 
164  if (insideOpA.nTris_ == 0)
165  {
166  outsideOp(thisTri);
167  return;
168  }
169 
170  if (outsideOpA.nTris_ == 0)
171  {
172  insideOp(thisTri);
173  return;
174  }
175 
176 
177  // Cut all triangles with edge 1.
178  // From insideTrisA to insideTrisB, outsideTrisA
179 
180  triIntersectionList insideTrisB;
181  label nInsideB = 0;
182  storeOp insideOpB(insideTrisB, nInsideB);
183 
184  //triIntersectionList outsideTrisB;
185  //label nOutsideB = 0;
186  //storeOp outsideOpB(outsideTrisB, nOutsideB);
187 
188  {
189  scalar s = Foam::mag(tgt.c() - tgt.b());
190  const plane pl0(tgt.b(), tgt.c(), tgt.c() + s*n);
191 
192  for (label i = 0; i < insideOpA.nTris_; i++)
193  {
194  const triPoints& tri = insideOpA.tris_[i];
195  triSliceWithPlane(pl0, tri, insideOpB, outsideOpA);
196  }
197 
200  //for (label i = 0; i < outsideOpA.nTris_; i++)
201  //{
202  // const triPoints& tri = outsideOpA.tris_[i];
203  // triSliceWithPlane(pl0, tri, outsideOpB, outsideOpB);
204  //}
205  }
206 
207 
208  // Cut all triangles with edge 2.
209  // From insideTrisB to insideTrisA, outsideTrisA
210  {
211  scalar s = Foam::mag(tgt.a() - tgt.c());
212  const plane pl0(tgt.c(), tgt.a(), tgt.a() + s*n);
213 
214  insideOpA.nTris_ = 0;
215  //outsideOpA.nTris_ = 0;
216  for (label i = 0; i < insideOpB.nTris_; i++)
217  {
218  const triPoints& tri = insideOpB.tris_[i];
219  triSliceWithPlane(pl0, tri, insideOpA, outsideOpA);
220  }
221 
224  //for (label i = 0; i < outsideOpB.nTris_; i++)
225  //{
226  // const triPoints& tri = outsideOpB.tris_[i];
227  // triSliceWithPlane(pl0, tri, outsideOpA, outsideOpA);
228  //}
229  }
230 
231  // Transfer from A to argument
232  for (label i = 0; i < insideOpA.nTris_; i++)
233  {
234  insideOp(insideOpA.tris_[i]);
235  }
236 
237  for (label i = 0; i < outsideOpA.nTris_; i++)
238  {
239  outsideOp(outsideOpA.tris_[i]);
240  }
241 }
242 
243 
244 // ************************************************************************* //
s
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.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))
Definition: gmvOutputSpray.H:25
Foam::triangle::storeOp
Store resulting tris.
Definition: triangle.H:125
triangle.H
triPoints.H
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
Foam::plane
Geometric class that creates a 3D plane and can return the intersection point between a line and the ...
Definition: plane.H:89
Foam::triangle
A triangle primitive used to calculate face normals and swept volumes.
Definition: triangle.H:59
n
label n
Definition: TABSMDCalcMethod2.H:31
Foam::triPoints
Triangle storage. Null constructable (unfortunately triangle<point, point> is not)
Definition: triPoints.H:52
Foam::triangle::c
const Point & c() const
Return third vertex.
Definition: triangleI.H:98
plane.H
Foam::triangle::a
const Point & a() const
Return first vertex.
Definition: triangleI.H:86
Foam::triangle::b
const Point & b() const
Return second vertex.
Definition: triangleI.H:92
Foam::Vector< scalar >
Foam::triangle::storeOp::nTris_
label & nTris_
Definition: triangle.H:129
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
Foam::FixedList
A 1D vector of objects of type <T> with a fixed length <N>.
Definition: HashTable.H:104
Foam::triangle::storeOp::tris_
triIntersectionList & tris_
Definition: triangle.H:128
Foam::point
vector point
Point is a vector.
Definition: point.H:43
Foam::triangle::sliceWithPlane
void sliceWithPlane(const plane &pln, AboveOp &aboveOp, BelowOp &belowOp) const
Decompose triangle into triangles above and below plane.
Definition: triangle.C:115
Foam::triangle::triangleOverlap
void triangleOverlap(const vector &n, const triangle< Point, PointRef > &tri, InsideOp &insideOp, OutsideOp &outsideOp) const
Decompose triangle into triangles inside and outside.
Definition: triangle.C:128