2 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
4 // + This file is part of enGrid. +
6 // + Copyright 2008,2009 Oliver Gloth +
8 // + enGrid is free software: you can redistribute it and/or modify +
9 // + it under the terms of the GNU General Public License as published by +
10 // + the Free Software Foundation, either version 3 of the License, or +
11 // + (at your option) any later version. +
13 // + enGrid is distributed in the hope that it will be useful, +
14 // + but WITHOUT ANY WARRANTY; without even the implied warranty of +
15 // + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +
16 // + GNU General Public License for more details. +
18 // + You should have received a copy of the GNU General Public License +
19 // + along with enGrid. If not, see <http://www.gnu.org/licenses/>. +
21 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
23 #include "geometrytools.h"
24 #include "containertricks.h"
26 #include <vtkCellType.h>
28 namespace GeometryTools
{
30 void rotate(vec3_t g1
, vec3_t g2
,vec3_t g3
, vec3_t
&b
, double theta
)
37 mat3_t g
= g_t
.transp();
41 mat3_t rot
= mat3_t::identity();
42 rot
[0][0] = cos(theta
);
43 rot
[0][1] = -sin(theta
);
44 rot
[1][0] = sin(theta
);
45 rot
[1][1] = cos(theta
);
46 //cout << rot << endl;
48 //cout << bbb << endl;
53 vec3_t
rotate(vec3_t v
, vec3_t axis
, double theta
)
57 // transposed base of rotate system
60 // compute projection of v in axis direction
61 vec3_t v_axis
= (axis
*v
)*axis
;
63 // compute first orthogonal vector (first base vector)
67 // second base vector is the normalised axis
70 // compute second orthogonal vector (third base vector)
71 g_t
[2] = g_t
[0].cross(g_t
[1]);
73 // base of rotate system
74 mat3_t g
= g_t
.transp();
76 // matrix for rotation around g_t[1];
77 mat3_t rot
= mat3_t::identity();
78 rot
[0][0] = cos(theta
);
79 rot
[0][2] = sin(theta
);
80 rot
[2][0] = -sin(theta
);
81 rot
[2][2] = cos(theta
);
83 // transfer v to rotate system
86 // rotate the vector and transfer it back
93 vec3_t
orthogonalVector(vec3_t v
)
99 for (int i
= 1; i
< 3; ++i
) {
100 if (v
[i
] > v
[i_max
]) i_max
= i
;
101 if (v
[i
] < v
[i_min
]) i_min
= i
;
112 double intersection(vec3_t x_straight
, vec3_t v_straight
, vec3_t x_plane
, vec3_t u_plane
, vec3_t v_plane
)
114 vec3_t n
= u_plane
.cross(v_plane
);
115 return intersection(x_straight
,v_straight
,x_plane
,n
);
118 double intersection(vec3_t x_straight
, vec3_t v_straight
, vec3_t x_plane
, vec3_t n_plane
)
120 double k
= (x_plane
*n_plane
- x_straight
*n_plane
)/(v_straight
*n_plane
);
124 bool intersection (double &k1
, double &k2
, vec2_t r1
, vec2_t u1
, vec2_t r2
, vec2_t u2
)
126 double ave_length
= .5*(sqrt(u1
[0]*u1
[0]+u1
[1]*u1
[1]) + sqrt(u2
[0]*u2
[0]+u2
[1]*u2
[1]));
127 double DET
= (u1
[0]*u2
[1]-u1
[1]*u2
[0]);
128 if (fabs(DET
) > 1e-6*ave_length
) {
129 k1
= -(u2
[0]*r2
[1]-u2
[0]*r1
[1]-r2
[0]*u2
[1]+r1
[0]*u2
[1])/DET
;
130 k2
= -(-u1
[1]*r2
[0]+u1
[0]*r2
[1]-u1
[0]*r1
[1]+u1
[1]*r1
[0])/DET
;
137 void sliceTriangle(const vector
<vec3_t
> &Tin
, vec3_t x
, vec3_t n
, vector
<vector
<vec3_t
> > &Tout
)
142 double kab
= intersection(a
,b
-a
,x
,n
);
143 double kbc
= intersection(b
,c
-b
,x
,n
);
144 double kca
= intersection(c
,a
-c
,x
,n
);
145 bool ab_cut
= ((kab
>= 0) && (kab
<= 1));
146 bool bc_cut
= ((kbc
>= 0) && (kbc
<= 1));
147 bool ca_cut
= ((kca
>= 0) && (kca
<= 1));
148 if (ab_cut
&& bc_cut
&& ca_cut
) {
149 //cerr << "invalid triangle (SliceTriangle) A" << endl;
150 //exit(EXIT_FAILURE);
151 if ((kab
<= kbc
) && (kab
<= kca
)) ab_cut
= false;
152 else if ((kbc
<= kab
) && (kbc
<= kca
)) bc_cut
= false;
155 if (ab_cut
&& bc_cut
) {
156 vec3_t ab
= a
+ kab
*(b
-a
);
157 vec3_t bc
= b
+ kbc
*(c
-b
);
158 Tout
.resize(3,vector
<vec3_t
>(3));
159 clinit(Tout
[0]) = a
,ab
,bc
;
160 clinit(Tout
[1]) = ab
,b
,bc
;
161 clinit(Tout
[2]) = bc
,c
,a
;
162 } else if (bc_cut
&& ca_cut
) {
163 vec3_t bc
= b
+ kbc
*(c
-b
);
164 vec3_t ca
= c
+ kca
*(a
-c
);
165 Tout
.resize(3,vector
<vec3_t
>(3));
166 clinit(Tout
[0]) = a
,bc
,ca
;
167 clinit(Tout
[1]) = a
,b
,bc
;
168 clinit(Tout
[2]) = bc
,c
,ca
;
169 } else if (ca_cut
&& ab_cut
) {
170 vec3_t ca
= c
+ kca
*(a
-c
);
171 vec3_t ab
= a
+ kab
*(b
-a
);
172 Tout
.resize(3,vector
<vec3_t
>(3));
173 clinit(Tout
[0]) = a
,ab
,ca
;
174 clinit(Tout
[1]) = ab
,b
,ca
;
175 clinit(Tout
[2]) = b
,c
,ca
;
177 Tout
.resize(1,vector
<vec3_t
>(3));
178 clinit(Tout
[0]) = a
,b
,c
;
183 double tetraVol(vec3_t x0
, vec3_t x1
, vec3_t x2
, vec3_t x3
, bool neg
)
188 double V
= v1
*(v2
.cross(v3
));
190 if (!neg
&& (V
< 0)) {
196 Lmin
= min(Lmin
,v1
.abs());
197 Lmin
= min(Lmin
,v2
.abs());
198 Lmin
= min(Lmin
,v3
.abs());
199 Lmin
= min(Lmin
,v4
.abs());
200 Lmin
= min(Lmin
,v5
.abs());
201 Lmin
= min(Lmin
,v6
.abs());
202 Lmax
= max(Lmax
,v1
.abs());
203 Lmax
= max(Lmax
,v2
.abs());
204 Lmax
= max(Lmax
,v3
.abs());
205 Lmax
= max(Lmax
,v4
.abs());
206 Lmax
= max(Lmax
,v5
.abs());
207 Lmax
= max(Lmax
,v6
.abs());
208 if (Lmin
/Lmax
< 1e-6) {
219 double pyraVol(vec3_t x0
, vec3_t x1
, vec3_t x2
, vec3_t x3
, vec3_t x4
, bool neg
)
221 double V1
= tetraVol(x0
, x1
, x3
, x4
, neg
) + tetraVol(x1
, x2
, x3
, x4
, neg
);
222 double V2
= tetraVol(x0
, x1
, x2
, x4
, neg
) + tetraVol(x2
, x3
, x0
, x4
, neg
);
226 vec3_t m0 = .25*(x0+x1+x2+x3);
227 V += tetraVol(x0, x1, m0, x4, neg);
228 V += tetraVol(x1, x2, m0, x4, neg);
229 V += tetraVol(x2, x3, m0, x4, neg);
230 V += tetraVol(x3, x0, m0, x4, neg);
235 double prismVol(vec3_t x0
, vec3_t x1
, vec3_t x2
, vec3_t x3
, vec3_t x4
, vec3_t x5
, bool neg
)
238 vec3_t p
= 1.0/6.0*(x0
+x1
+x2
+x3
+x4
+x5
);
239 V
+= tetraVol(x0
, x2
, x1
, p
, neg
);
240 V
+= tetraVol(x3
, x4
, x5
, p
, neg
);
241 V
+= pyraVol (x0
, x1
, x4
, x3
, p
, neg
);
242 V
+= pyraVol (x1
, x2
, x5
, x4
, p
, neg
);
243 V
+= pyraVol (x0
, x3
, x5
, x2
, p
, neg
);
247 double hexaVol(vec3_t x0
, vec3_t x1
, vec3_t x2
, vec3_t x3
, vec3_t x4
, vec3_t x5
, vec3_t x6
, vec3_t x7
, bool neg
)
250 vec3_t p
= 1.0/8.0*(x0
+x1
+x2
+x3
+x4
+x5
+x6
+x7
);
251 V
+= pyraVol(x0
, x1
, x3
, x2
, p
, neg
);
252 V
+= pyraVol(x0
, x4
, x5
, x1
, p
, neg
);
253 V
+= pyraVol(x4
, x6
, x7
, x5
, p
, neg
);
254 V
+= pyraVol(x2
, x3
, x7
, x6
, p
, neg
);
255 V
+= pyraVol(x1
, x5
, x7
, x3
, p
, neg
);
256 V
+= pyraVol(x0
, x2
, x6
, x4
, p
, neg
);
260 double triArea(vec3_t x0
, vec3_t x1
, vec3_t x2
)
264 double A
= 0.5*((a
.cross(b
)).abs());
268 double quadArea(vec3_t x0
, vec3_t x1
, vec3_t x2
, vec3_t x3
)
271 vec3_t p
= .25*(x0
+x1
+x2
+x3
);
272 A
+= triArea(x0
,x1
,p
);
273 A
+= triArea(x1
,x2
,p
);
274 A
+= triArea(x2
,x3
,p
);
275 A
+= triArea(x3
,x0
,p
);
279 vec3_t
triNormal(vec3_t x0
, vec3_t x1
, vec3_t x2
)
283 vec3_t n
= 0.5*(a
.cross(b
));
287 vec3_t
quadNormal(vec3_t x0
, vec3_t x1
, vec3_t x2
, vec3_t x3
)
291 vec3_t p
= .25*(x0
+x1
+x2
+x3
);
292 n
+= triNormal(x0
,x1
,p
);
293 n
+= triNormal(x1
,x2
,p
);
294 n
+= triNormal(x2
,x3
,p
);
295 n
+= triNormal(x3
,x0
,p
);
299 vec3_t
triNormal(vtkUnstructuredGrid
*grid
, vtkIdType p1
, vtkIdType p2
, vtkIdType p3
)
302 grid
->GetPoint(p1
,x1
.data());
303 grid
->GetPoint(p2
,x2
.data());
304 grid
->GetPoint(p3
,x3
.data());
305 return triNormal(x1
,x2
,x3
);
308 vec3_t
quadNormal(vtkUnstructuredGrid
*grid
, vtkIdType p1
, vtkIdType p2
, vtkIdType p3
, vtkIdType p4
)
310 vec3_t x1
, x2
, x3
, x4
;
311 grid
->GetPoint(p1
,x1
.data());
312 grid
->GetPoint(p2
,x2
.data());
313 grid
->GetPoint(p3
,x3
.data());
314 grid
->GetPoint(p4
,x4
.data());
315 return quadNormal(x1
,x2
,x3
,x4
);
318 vec3_t
cellNormal(vtkUnstructuredGrid
*grid
, vtkIdType i
)
324 grid
->GetCellPoints(i
,npts
,pts
);
325 if (npts
== 3) return triNormal(grid
,pts
[0],pts
[1],pts
[2]);
326 else if (npts
== 4) return quadNormal(grid
,pts
[0],pts
[1],pts
[2],pts
[3]);
330 double cellVA(vtkUnstructuredGrid
*grid
, vtkIdType cellId
, bool neg
)
335 grid
->GetCellPoints(cellId
, Npts
, pts
);
336 for (int i_pts
= 0; i_pts
< Npts
; ++i_pts
) {
337 grid
->GetPoints()->GetPoint(pts
[i_pts
], p
[i_pts
].data());
339 vtkIdType cellType
= grid
->GetCellType(cellId
);
340 if (cellType
== VTK_TRIANGLE
) return triArea (p
[0], p
[1], p
[2]);
341 else if (cellType
== VTK_QUAD
) return quadArea(p
[0], p
[1], p
[2], p
[3]);
342 else if (cellType
== VTK_TETRA
) return tetraVol(p
[0], p
[1], p
[2], p
[3], neg
);
343 else if (cellType
== VTK_PYRAMID
) return pyraVol (p
[0], p
[1], p
[2], p
[3], p
[4], neg
);
344 else if (cellType
== VTK_WEDGE
) return prismVol(p
[0], p
[1], p
[2], p
[3], p
[4], p
[5], neg
);
345 else if (cellType
== VTK_HEXAHEDRON
) return hexaVol (p
[0], p
[1], p
[2], p
[3], p
[4], p
[5], p
[6], p
[7], neg
);