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-------------------------------------------------------------------------------
10License
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
34template<class Point, class PointRef>
35template<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
112template<class Point, class PointRef>
113template<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
125template<class Point, class PointRef>
126template<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// ************************************************************************* //
label n
A 1D vector of objects of type <T> with a fixed length <N>.
Definition: FixedList.H:81
Geometric class that creates a 3D plane and can return the intersection point between a line and the ...
Definition: plane.H:95
Triangle storage. Null constructable (unfortunately triangle<point, point> is not)
Definition: triPoints.H:55
Store resulting tris.
Definition: triangle.H:126
triIntersectionList & tris_
Definition: triangle.H:128
A triangle primitive used to calculate face normals and swept volumes.
Definition: triangle.H:80
void sliceWithPlane(const plane &pln, AboveOp &aboveOp, BelowOp &belowOp) const
Decompose triangle into triangles above and below plane.
Definition: triangle.C:115
const Point & a() const
Return first vertex.
Definition: triangleI.H:86
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
const Point & c() const
Return third vertex.
Definition: triangleI.H:98
const Point & b() const
Return second vertex.
Definition: triangleI.H:92
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))
vector point
Point is a vector.
Definition: point.H:43
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333