?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
updated to modern VTK
[engrid-github.git] / src / libengrid / swaptriangles.cpp
blobfcd1e00788758a20a22d01239dfa11075d6f0d76
1 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 // + +
3 // + This file is part of enGrid. +
4 // + +
5 // + Copyright 2008-2014 enGits GmbH +
6 // + +
7 // + enGrid is free software: you can redistribute it and/or modify +
8 // + it under the terms of the GNU General Public License as published by +
9 // + the Free Software Foundation, either version 3 of the License, or +
10 // + (at your option) any later version. +
11 // + +
12 // + enGrid is distributed in the hope that it will be useful, +
13 // + but WITHOUT ANY WARRANTY; without even the implied warranty of +
14 // + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +
15 // + GNU General Public License for more details. +
16 // + +
17 // + You should have received a copy of the GNU General Public License +
18 // + along with enGrid. If not, see <http://www.gnu.org/licenses/>. +
19 // + +
20 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
21 #include "swaptriangles.h"
22 #include "guimainwindow.h"
23
24 #include <vtkCellArray.h>
25
26 #include "checksurfaceintegrity.h"
27 #include "statistics.h"
28
29 using namespace GeometryTools;
30
31 SwapTriangles::SwapTriangles() : SurfaceOperation()
32 {
33 m_RespectBC = false;
34 m_FeatureSwap = false;
35 m_FeatureAngle = GeometryTools::deg2rad(30);
36 m_MaxNumLoops = 20;
37 m_SmallAreaSwap = false;
38 m_SmallAreaRatio = 1e-3;
39 m_Verbose = false;
40 m_DelaunayThreshold = 1.0;
41 getSet("surface meshing", "small area ratio for edge-swapping", 1e-3, m_SmallAreaRatio);
42 getSet("surface meshing", "surface error threshold", 0.05, m_SurfErrorThreshold);
43 m_SurfErrorThreshold = deg2rad(m_SurfErrorThreshold);
44 }
45
46 bool SwapTriangles::testOrientation(stencil_t S)
47 {
48 // old triangles
49 vec3_t n1_old = triNormal(m_Grid, S.id_node[0], S.p1, S.p2);
50 vec3_t n2_old = triNormal(m_Grid, S.id_node[1], S.p2, S.p1);
51
52 // new triangles
53 vec3_t n1_new = triNormal(m_Grid, S.p1, S.id_node[1], S.id_node[0]);
54 vec3_t n2_new = triNormal(m_Grid, S.p2, S.id_node[0], S.id_node[1]);
55
56 // top point
57 vec3_t x_summit(0,0,0);
58 vec3_t x[4];
59 double l_max = 0;
60 m_Grid->GetPoints()->GetPoint(S.id_node[0], x[0].data());
61 m_Grid->GetPoints()->GetPoint(S.p1, x[1].data());
62 m_Grid->GetPoints()->GetPoint(S.id_node[1], x[2].data());
63 m_Grid->GetPoints()->GetPoint(S.p2, x[3].data());
64 for (int k = 0; k < 4; ++k) {
65 x_summit += x[k];
66 }
67 for (int k = 0; k < 4; ++k) {
68 int i = k;
69 int j = k + 1;
70 if (j == 4) {
71 j = 0;
72 }
73 l_max = max(l_max, (x[i]-x[j]).abs());
74 }
75 x_summit *= 0.25;
76 vec3_t n = n1_old + n2_old;
77 n.normalise();
78 x_summit += 3*l_max*n;
79
80 // old volumes
81 double V1_old = tetraVol(x[0], x[1], x[3], x_summit, true);
82 double V2_old = tetraVol(x[2], x[3], x[1], x_summit, true);
83 // new volumes
84 double V1_new = tetraVol(x[1], x[2], x[0], x_summit, true);
85 double V2_new = tetraVol(x[3], x[0], x[2], x_summit, true);
86
87 bool swap_ok = false;
88 if (m_SmallAreaSwap) {
89 swap_ok = V1_new>0 && V2_new>0;
90 } else {
91 swap_ok = V1_old>0 && V2_old>0 && V1_new>0 && V2_new>0;
92 }
93 return swap_ok;
94 }
95
96 bool SwapTriangles::isEdge(vtkIdType id_node1, vtkIdType id_node2)
97 {
98 l2g_t nodes = getPartNodes();
99 g2l_t _nodes = getPartLocalNodes();
100 l2l_t n2n = getPartN2N();
101
102 bool ret = false;
103 foreach(int i_node, n2n[_nodes[id_node1]]) {
104 vtkIdType id_node = nodes[i_node];
105 if( id_node == id_node2 ) ret = true;
106 }
107 return(ret);
108 }
109
110 double SwapTriangles::computeSurfaceDistance(vec3_t x1, vec3_t x2, vec3_t x3, CadInterface *cad_interface)
111 {
112 static const double f13 = 1.0/3.0;
113 vec3_t x = f13*(x1 + x2 + x3);
114 vec3_t n = triNormal(x1, x2, x3);
115 n.normalise();
116 if (!checkVector(n)) {
117 return 0;
118 }
119 vec3_t xp = cad_interface->snap(x, false);
120 return (x - xp).abs();
121 }
122
123 double SwapTriangles::edgeAngle(vec3_t x1, vec3_t x2, vec3_t x3, vec3_t x4)
124 {
125 static const double f13 = 1.0/3.0;
126 vec3_t xf1 = f13*(x1 + x2 + x4);
127 vec3_t xf2 = f13*(x2 + x3 + x4);
128 vec3_t xe = 0.5*(x2 + x4);
129 vec3_t xfe = 0.5*(xf1 + xf2);
130 vec3_t n1 = triNormal(x1, x2, x4);
131 vec3_t n2 = triNormal(x2, x3, x4);
132 vec3_t ve = xe - xfe;
133 double alpha = angle(n1, n2);
134 if ((n1 + n2)*ve < 0) {
135 return -alpha;
136 }
137 return alpha;
138 }
139
140 vtkIdType SwapTriangles::neighbourNode(vtkIdType id_node0, vtkIdType id_node1, vtkIdType id_node2)
141 {
142 for (int i = 0; i < m_Part.n2nGSize(id_node1); ++i) {
143 vtkIdType id_neigh = m_Part.n2nGG(id_node1, i);
144 if (id_neigh != id_node0) {
145 if (m_Part.hasNeighNode(id_node2, id_neigh)) {
146 return id_neigh;
147 }
148 }
149 }
150 EG_BUG;
151 return -1;
152 }
153
154 bool SwapTriangles::swapDueToSurfaceNoise(stencil_t S)
155 {
156 vtkIdType p1 = S.id_node[0];
157 vtkIdType p2 = S.p1;
158 vtkIdType p3 = S.id_node[1];
159 vtkIdType p4 = S.p2;
160 vtkIdType p12 = neighbourNode(p4, p1, p2);
161 vtkIdType p23 = neighbourNode(p4, p2, p3);
162 vtkIdType p34 = neighbourNode(p2, p3, p4);
163 vtkIdType p41 = neighbourNode(p2, p4, p1);
164 vec3_t x1, x2, x3, x4, x12, x23, x34, x41;
165 m_Grid->GetPoint(p1, x1.data());
166 m_Grid->GetPoint(p2, x2.data());
167 m_Grid->GetPoint(p3, x3.data());
168 m_Grid->GetPoint(p4, x4.data());
169 m_Grid->GetPoint(p12, x12.data());
170 m_Grid->GetPoint(p23, x23.data());
171 m_Grid->GetPoint(p34, x34.data());
172 m_Grid->GetPoint(p41, x41.data());
173
174 double noise1, noise2;
175 {
176 double alpha_min = M_PI;
177 double alpha_max = -M_PI;
178 double alpha;
179
180 alpha = edgeAngle(x4, x1, x12, x2);
181 alpha_max = max(alpha, alpha_max);
182 alpha_min = min(alpha, alpha_min);
183
184 alpha = edgeAngle(x4, x2, x23, x3);
185 alpha_max = max(alpha, alpha_max);
186 alpha_min = min(alpha, alpha_min);
187
188 alpha = edgeAngle(x2, x3, x34, x4);
189 alpha_max = max(alpha, alpha_max);
190 alpha_min = min(alpha, alpha_min);
191
192 alpha = edgeAngle(x2, x4, x41, x1);
193 alpha_max = max(alpha, alpha_max);
194 alpha_min = min(alpha, alpha_min);
195
196 alpha = edgeAngle(x1, x2, x3, x4);
197 alpha_max = max(alpha, alpha_max);
198 alpha_min = min(alpha, alpha_min);
199
200 noise1 = alpha_max - alpha_min;
201 }
202 {
203 double alpha_min = M_PI;
204 double alpha_max = -M_PI;
205 double alpha;
206
207 alpha = edgeAngle(x3, x1, x12, x2);
208 alpha_max = max(alpha, alpha_max);
209 alpha_min = min(alpha, alpha_min);
210
211 alpha = edgeAngle(x1, x2, x23, x3);
212 alpha_max = max(alpha, alpha_max);
213 alpha_min = min(alpha, alpha_min);
214
215 alpha = edgeAngle(x1, x3, x34, x4);
216 alpha_max = max(alpha, alpha_max);
217 alpha_min = min(alpha, alpha_min);
218
219 alpha = edgeAngle(x3, x4, x41, x1);
220 alpha_max = max(alpha, alpha_max);
221 alpha_min = min(alpha, alpha_min);
222
223 alpha = edgeAngle(x2, x3, x4, x1);
224 alpha_max = max(alpha, alpha_max);
225 alpha_min = min(alpha, alpha_min);
226
227 noise2 = alpha_max - alpha_min;
228 }
229
230 if (noise1 > m_SurfErrorThreshold && noise1 > noise2) {
231 return true;
232 }
233
234 return false;
235 }
236
237 void SwapTriangles::computeSurfaceErrors(const QVector<vec3_t> &x, int bc, double &err1, double &err2)
238 {
239 using namespace GeometryTools;
240 err1 = 0;
241 err2 = 1e10;
242 CadInterface* cad_interface = GuiMainWindow::pointer()->getCadInterface(bc);
243 if (!cad_interface) {
244 return;
245 }
246
247 double d013 = computeSurfaceDistance(x[0], x[1], x[3], cad_interface);
248 double d123 = computeSurfaceDistance(x[1], x[2], x[3], cad_interface);
249 double d012 = computeSurfaceDistance(x[0], x[1], x[2], cad_interface);
250 double d023 = computeSurfaceDistance(x[0], x[2], x[3], cad_interface);
251
252 double L = max((x[1] - x[3]).abs(), (x[0] - x[2]).abs());
253 err1 = max(d013, d123)/L;
254 err2 = max(d012, d023)/L;
255 if (err1 < m_SurfErrorThreshold) {
256 err1 = 0;
257 }
258 if (err2 < m_SurfErrorThreshold) {
259 err2 = 0;
260 }
261
262 }
263
264 int SwapTriangles::swap()
265 {
266 int N_swaps = 0;
267 setAllSurfaceCells();
268 //computeAverageSurfaceError();
269 EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
270 EG_VTKDCN(vtkCharArray_t, node_type, m_Grid, "node_type");
271 QVector<bool> marked(m_Grid->GetNumberOfCells(), false);
272 for (int i = 0; i < m_Part.getNumberOfCells(); ++i) {
273 vtkIdType id_cell = m_Part.globalCell(i);
274 if (!m_BoundaryCodes.contains(cell_code->GetValue(id_cell)) && m_Grid->GetCellType(id_cell) == VTK_TRIANGLE) { //if it is a selected triangle
275 if (!marked[id_cell] && !m_Swapped[id_cell]) {
276 for (int j = 0; j < 3; ++j) {
277 bool swap = false;
278 stencil_t S = getStencil(id_cell, j);
279 if(S.id_cell.size() == 2 && S.sameBC) {
280 if (S.type_cell[1] == VTK_TRIANGLE) {
281 if(!isEdge(S.id_node[0], S.id_node[1]) ) {
282 //if (!marked[S.id_cell[1]] && !m_Swapped[S.id_cell[1]] && (!isFeatureNode(S.id_node[0]) || !isFeatureNode(S.id_node[1]))) {
283 if (!marked[S.id_cell[1]] && !m_Swapped[S.id_cell[1]] && (!isFeatureNode(S.p1) || !isFeatureNode(S.p2))) {
284 QVector<vec3_t> x3(4);
285 vec2_t x[4];
286
287 m_Grid->GetPoint(S.id_node[0], x3[0].data());
288 m_Grid->GetPoint(S.p1, x3[1].data());
289 m_Grid->GetPoint(S.id_node[1], x3[2].data());
290 m_Grid->GetPoint(S.p2, x3[3].data());
291
292 vec3_t n1 = triNormal(x3[0], x3[1], x3[3]);
293 vec3_t n2 = triNormal(x3[1], x3[2], x3[3]);
294
295 if (m_SmallAreaSwap) {
296 double A1 = n1.abs();
297 double A2 = n2.abs();
298 if (isnan(A1) || isnan(A2)) {
299 swap = true;
300 } else {
301 swap = A1 < m_SmallAreaRatio*A2 || A2 < m_SmallAreaRatio*A1;
302 }
303 }
304
305 if (GeometryTools::angle(n1, n2) < m_FeatureAngle) {
306
307 // Delaunay
308 if (!swap) {
309 vec3_t n = n1 + n2;
310 n.normalise();
311 vec3_t ex = orthogonalVector(n);
312 vec3_t ey = ex.cross(n);
313 for (int k = 0; k < 4; ++k) {
314 x[k] = vec2_t(x3[k]*ex, x3[k]*ey);
315 }
316 vec2_t r1, r2, r3, u1, u2, u3;
317 r1 = 0.5*(x[0] + x[1]); u1 = turnLeft(x[1] - x[0]);
318 r2 = 0.5*(x[1] + x[2]); u2 = turnLeft(x[2] - x[1]);
319 r3 = 0.5*(x[1] + x[3]); u3 = turnLeft(x[3] - x[1]);
320 double k, l;
321 vec2_t xm1, xm2;
322 bool ok = true;
323 if (intersection(k, l, r1, u1, r3, u3)) {
324 xm1 = r1 + k*u1;
325 if (intersection(k, l, r2, u2, r3, u3)) {
326 xm2 = r2 + k*u2;
327 } else {
328 ok = false;
329 }
330 } else {
331 ok = false;
332 swap = true;
333 }
334 if (ok) {
335 if ((xm1 - x[2]).abs() < (xm1 - x[0]).abs()/m_DelaunayThreshold) {
336 swap = true;
337 }
338 if ((xm2 - x[0]).abs() < (xm2 - x[2]).abs()/m_DelaunayThreshold) {
339 swap = true;
340 }
341 }
342 }
343 } //end of if feature angle
344 } //end of if l_marked
345 } //end of if TestSwap
346 }
347 } //end of S valid
348
349 if (swap) {
350 if (testOrientation(S)) {
351 marked[S.id_cell[0]] = true;
352 marked[S.id_cell[1]] = true;
353 for (int k = 0; k < m_Part.n2cGSize(S.id_node[0]); ++k) {
354 vtkIdType id_neigh = m_Part.n2cGG(S.id_node[0], k);
355 marked[id_neigh] = true;
356 }
357 for (int k = 0; k < m_Part.n2cGSize(S.id_node[1]); ++k) {
358 vtkIdType id_neigh = m_Part.n2cGG(S.id_node[1], k);
359 marked[id_neigh] = true;
360 }
361 for (int k = 0; k < m_Part.n2cGSize(S.p1); ++k) {
362 vtkIdType id_neigh = m_Part.n2cGG(S.p1, k);
363 marked[id_neigh] = true;
364 }
365 for (int k = 0; k < m_Part.n2cGSize(S.p2); ++k) {
366 vtkIdType id_neigh = m_Part.n2cGG(S.p2, k);
367 marked[id_neigh] = true;
368 }
369 vtkIdType new_pts1[3], new_pts2[3];
370 new_pts1[0] = S.p1;
371 new_pts1[1] = S.id_node[1];
372 new_pts1[2] = S.id_node[0];
373 new_pts2[0] = S.id_node[1];
374 new_pts2[1] = S.p2;
375 new_pts2[2] = S.id_node[0];
376 vtkIdType old_pts1[3], old_pts2[3];
377 old_pts1[0] = S.id_node[0];
378 old_pts1[1] = S.p1;
379 old_pts1[2] = S.p2;
380 old_pts2[0] = S.id_node[1];
381 old_pts2[1] = S.p2;
382 old_pts2[2] = S.p1;
383 m_Grid->ReplaceCell(S.id_cell[0], 3, new_pts1);
384 m_Grid->ReplaceCell(S.id_cell[1], 3, new_pts2);
385 m_Swapped[S.id_cell[0]] = true;
386 m_Swapped[S.id_cell[1]] = true;
387 ++N_swaps;
388 break;
389 }
390 } //end of if swap
391 } //end of loop through sides
392 } //end of if marked
393 } //end of if selected triangle
394 } //end of loop through cells
395 return N_swaps;
396 }
397
398 void SwapTriangles::computeAverageSurfaceError()
399 {
400 EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
401 QList<double> errors;
402 EG_FORALL_NODES(id_node1, m_Grid) {
403 vec3_t n1 = m_Part.globalNormal(id_node1);
404 vec3_t x1;
405 m_Grid->GetPoint(id_node1, x1.data());
406 for (int i = 0; i < m_Part.n2nGSize(id_node1); ++i) {
407 vtkIdType id_node2 = m_Part.n2nGG(id_node1, i);
408 if (id_node1 < id_node2) {
409 QVector<vtkIdType> faces;
410 getEdgeCells(id_node1, id_node2, faces);
411 if (faces.size() == 2) {
412 int bc1 = cell_code->GetValue(faces[0]);
413 int bc2 = cell_code->GetValue(faces[1]);
414 if (bc1 == bc2) {
415 CadInterface* cad_interface = GuiMainWindow::pointer()->getCadInterface(bc1);
416 if (cad_interface) {
417 vec3_t n2 = m_Part.globalNormal(id_node2);
418 vec3_t n = 0.5*(n1 + n2);
419 vec3_t x2;
420 m_Grid->GetPoint(id_node2, x2.data());
421 vec3_t x = 0.5*(x1 + x2);
422 vec3_t xp = cad_interface->snap(x);
423 double err = (x - xp).abs()/(x1 - x2).abs();
424 errors.append(err);
425 }
426 }
427 }
428 }
429 }
430 }
431
432 m_AverageSurfaceError = Statistics::meanValue(errors);
433 m_SurfaceErrorDeviation = Statistics::standardDeviation(errors, m_AverageSurfaceError);
434 cout << "mean: " << m_AverageSurfaceError << " sigma: " << m_SurfaceErrorDeviation << endl;
435 }
436
437 void SwapTriangles::operate()
438 {
439 if (m_Verbose) cout << "swapping edges for surface triangles ..." << endl;
440 long int N_swaps = 100000000;
441 long int N_last_swaps = 100000001;
442 int loop = 1;
443 m_NumSwaps = 0;
444 while ((N_swaps > 0) && (loop <= m_MaxNumLoops) && (N_swaps < N_last_swaps)) {
445 N_last_swaps = N_swaps;
446 if (m_Verbose) cout << " loop " << loop << "/" << m_MaxNumLoops << endl;
447 m_Swapped.fill(false, m_Grid->GetNumberOfCells());
448 N_swaps = 0;
449 int N;
450 int sub_loop = 0;
451 do {
452 ++sub_loop;
453 N = swap();
454 m_NumSwaps = max(N, m_NumSwaps);
455 if (m_Verbose) cout << " sub-loop " << sub_loop << ": " << N << " swaps" << endl;
456 N_swaps += N;
457 } while (N > 0);
458 ++loop;
459 }
460 }
461
This is a mirror repository of the new location for engrid@github