?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
new vertex type for feature corners + geometric recogintion
[engrid.git] / src / libengrid / swaptriangles.cpp
blob5964237315d93db6bc9b265d69462079ed2c8d0d
1 //
2 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3 // + +
4 // + This file is part of enGrid. +
5 // + +
6 // + Copyright 2008-2012 enGits GmbH +
7 // + +
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. +
12 // + +
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. +
17 // + +
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/>. +
20 // + +
21 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
22 //
23 #include "swaptriangles.h"
24 #include "guimainwindow.h"
25
26 #include <vtkCellArray.h>
27
28 #include "checksurfaceintegrity.h"
29 #include "statistics.h"
30
31 using namespace GeometryTools;
32
33 SwapTriangles::SwapTriangles() : SurfaceOperation()
34 {
35 m_RespectBC = false;
36 m_FeatureSwap = false;
37 m_FeatureAngle = GeometryTools::deg2rad(30);
38 m_MaxNumLoops = 20;
39 m_SmallAreaSwap = false;
40 m_SmallAreaRatio = 1e-3;
41 m_Verbose = false;
42 getSet("surface meshing", "small area ratio for edge-swapping", 1e-3, m_SmallAreaRatio);
43 getSet("surface meshing", "surface error threshold", 0.05, m_SurfErrorThreshold);
44 m_SurfErrorThreshold = deg2rad(m_SurfErrorThreshold);
45 }
46
47 bool SwapTriangles::testOrientation(stencil_t S)
48 {
49 // old triangles
50 vec3_t n1_old = triNormal(m_Grid, S.id_node[0], S.p1, S.p2);
51 vec3_t n2_old = triNormal(m_Grid, S.id_node[1], S.p2, S.p1);
52
53 // new triangles
54 vec3_t n1_new = triNormal(m_Grid, S.p1, S.id_node[1], S.id_node[0]);
55 vec3_t n2_new = triNormal(m_Grid, S.p2, S.id_node[0], S.id_node[1]);
56
57 // top point
58 vec3_t x_summit(0,0,0);
59 vec3_t x[4];
60 double l_max = 0;
61 m_Grid->GetPoints()->GetPoint(S.id_node[0], x[0].data());
62 m_Grid->GetPoints()->GetPoint(S.p1, x[1].data());
63 m_Grid->GetPoints()->GetPoint(S.id_node[1], x[2].data());
64 m_Grid->GetPoints()->GetPoint(S.p2, x[3].data());
65 for (int k = 0; k < 4; ++k) {
66 x_summit += x[k];
67 }
68 for (int k = 0; k < 4; ++k) {
69 int i = k;
70 int j = k + 1;
71 if (j == 4) {
72 j = 0;
73 }
74 l_max = max(l_max, (x[i]-x[j]).abs());
75 }
76 x_summit *= 0.25;
77 vec3_t n = n1_old + n2_old;
78 n.normalise();
79 x_summit += 3*l_max*n;
80
81 // old volumes
82 double V1_old = tetraVol(x[0], x[1], x[3], x_summit, true);
83 double V2_old = tetraVol(x[2], x[3], x[1], x_summit, true);
84 // new volumes
85 double V1_new = tetraVol(x[1], x[2], x[0], x_summit, true);
86 double V2_new = tetraVol(x[3], x[0], x[2], x_summit, true);
87
88 bool swap_ok = false;
89 if (m_SmallAreaSwap) {
90 swap_ok = V1_new>0 && V2_new>0;
91 } else {
92 swap_ok = V1_old>0 && V2_old>0 && V1_new>0 && V2_new>0;
93 }
94 return swap_ok;
95 }
96
97 bool SwapTriangles::isEdge(vtkIdType id_node1, vtkIdType id_node2)
98 {
99 l2g_t nodes = getPartNodes();
100 g2l_t _nodes = getPartLocalNodes();
101 l2l_t n2n = getPartN2N();
102
103 bool ret = false;
104 foreach(int i_node, n2n[_nodes[id_node1]]) {
105 vtkIdType id_node = nodes[i_node];
106 if( id_node == id_node2 ) ret = true;
107 }
108 return(ret);
109 }
110
111 double SwapTriangles::computeSurfaceDistance(vec3_t x1, vec3_t x2, vec3_t x3, SurfaceProjection* proj)
112 {
113 static const double f13 = 1.0/3.0;
114 vec3_t x = f13*(x1 + x2 + x3);
115 vec3_t n = triNormal(x1, x2, x3);
116 n.normalise();
117 if (!checkVector(n)) {
118 return 0;
119 }
120 vec3_t xp = proj->project(x, -1, true, n);
121 n = proj->lastProjNormal();
122 xp = proj->project(x, -1, true, n);
123 return (x - xp).abs();
124 }
125
126 double SwapTriangles::edgeAngle(vec3_t x1, vec3_t x2, vec3_t x3, vec3_t x4)
127 {
128 static const double f13 = 1.0/3.0;
129 vec3_t xf1 = f13*(x1 + x2 + x4);
130 vec3_t xf2 = f13*(x2 + x3 + x4);
131 vec3_t xe = 0.5*(x2 + x4);
132 vec3_t xfe = 0.5*(xf1 + xf2);
133 vec3_t n1 = triNormal(x1, x2, x4);
134 vec3_t n2 = triNormal(x2, x3, x4);
135 vec3_t ve = xe - xfe;
136 double alpha = angle(n1, n2);
137 if ((n1 + n2)*ve < 0) {
138 return -alpha;
139 }
140 return alpha;
141 }
142
143 vtkIdType SwapTriangles::neighbourNode(vtkIdType id_node0, vtkIdType id_node1, vtkIdType id_node2)
144 {
145 for (int i = 0; i < m_Part.n2nGSize(id_node1); ++i) {
146 vtkIdType id_neigh = m_Part.n2nGG(id_node1, i);
147 if (id_neigh != id_node0) {
148 if (m_Part.hasNeighNode(id_node2, id_neigh)) {
149 return id_neigh;
150 }
151 }
152 }
153 EG_BUG;
154 return -1;
155 }
156
157 bool SwapTriangles::swapDueToSurfaceNoise(stencil_t S)
158 {
159 vtkIdType p1 = S.id_node[0];
160 vtkIdType p2 = S.p1;
161 vtkIdType p3 = S.id_node[1];
162 vtkIdType p4 = S.p2;
163 vtkIdType p12 = neighbourNode(p4, p1, p2);
164 vtkIdType p23 = neighbourNode(p4, p2, p3);
165 vtkIdType p34 = neighbourNode(p2, p3, p4);
166 vtkIdType p41 = neighbourNode(p2, p4, p1);
167 vec3_t x1, x2, x3, x4, x12, x23, x34, x41;
168 m_Grid->GetPoint(p1, x1.data());
169 m_Grid->GetPoint(p2, x2.data());
170 m_Grid->GetPoint(p3, x3.data());
171 m_Grid->GetPoint(p4, x4.data());
172 m_Grid->GetPoint(p12, x12.data());
173 m_Grid->GetPoint(p23, x23.data());
174 m_Grid->GetPoint(p34, x34.data());
175 m_Grid->GetPoint(p41, x41.data());
176
177 double noise1, noise2;
178 {
179 double alpha_min = M_PI;
180 double alpha_max = -M_PI;
181 double alpha;
182
183 alpha = edgeAngle(x4, x1, x12, x2);
184 alpha_max = max(alpha, alpha_max);
185 alpha_min = min(alpha, alpha_min);
186
187 alpha = edgeAngle(x4, x2, x23, x3);
188 alpha_max = max(alpha, alpha_max);
189 alpha_min = min(alpha, alpha_min);
190
191 alpha = edgeAngle(x2, x3, x34, x4);
192 alpha_max = max(alpha, alpha_max);
193 alpha_min = min(alpha, alpha_min);
194
195 alpha = edgeAngle(x2, x4, x41, x1);
196 alpha_max = max(alpha, alpha_max);
197 alpha_min = min(alpha, alpha_min);
198
199 alpha = edgeAngle(x1, x2, x3, x4);
200 alpha_max = max(alpha, alpha_max);
201 alpha_min = min(alpha, alpha_min);
202
203 noise1 = alpha_max - alpha_min;
204 }
205 {
206 double alpha_min = M_PI;
207 double alpha_max = -M_PI;
208 double alpha;
209
210 alpha = edgeAngle(x3, x1, x12, x2);
211 alpha_max = max(alpha, alpha_max);
212 alpha_min = min(alpha, alpha_min);
213
214 alpha = edgeAngle(x1, x2, x23, x3);
215 alpha_max = max(alpha, alpha_max);
216 alpha_min = min(alpha, alpha_min);
217
218 alpha = edgeAngle(x1, x3, x34, x4);
219 alpha_max = max(alpha, alpha_max);
220 alpha_min = min(alpha, alpha_min);
221
222 alpha = edgeAngle(x3, x4, x41, x1);
223 alpha_max = max(alpha, alpha_max);
224 alpha_min = min(alpha, alpha_min);
225
226 alpha = edgeAngle(x2, x3, x4, x1);
227 alpha_max = max(alpha, alpha_max);
228 alpha_min = min(alpha, alpha_min);
229
230 noise2 = alpha_max - alpha_min;
231 }
232
233 if (noise1 > m_SurfErrorThreshold && noise1 > noise2) {
234 return true;
235 }
236
237 return false;
238 }
239
240 void SwapTriangles::computeSurfaceErrors(const QVector<vec3_t> &x, int bc, double &err1, double &err2)
241 {
242 using namespace GeometryTools;
243 err1 = 0;
244 err2 = 1e10;
245 SurfaceProjection* proj = GuiMainWindow::pointer()->getSurfProj(bc);
246 if (!proj) {
247 return;
248 }
249
250 double d013 = computeSurfaceDistance(x[0], x[1], x[3], proj);
251 double d123 = computeSurfaceDistance(x[1], x[2], x[3], proj);
252 double d012 = computeSurfaceDistance(x[0], x[1], x[2], proj);
253 double d023 = computeSurfaceDistance(x[0], x[2], x[3], proj);
254
255 double L = max((x[1] - x[3]).abs(), (x[0] - x[2]).abs());
256 err1 = max(d013, d123)/L;
257 err2 = max(d012, d023)/L;
258 if (err1 < m_SurfErrorThreshold) {
259 err1 = 0;
260 }
261 if (err2 < m_SurfErrorThreshold) {
262 err2 = 0;
263 }
264
265 /*
266 // new version based on surface normals
267 vec3_t n013 = triNormal(x[0], x[1], x[3]);
268 vec3_t n123 = triNormal(x[1], x[2], x[3]);
269 vec3_t n012 = triNormal(x[0], x[1], x[2]);
270 vec3_t n023 = triNormal(x[0], x[2], x[3]);
271 n013.normalise();
272 n123.normalise();
273 n012.normalise();
274 n023.normalise();
275 if (!checkVector(n013)) {
276 return;
277 }
278 if (!checkVector(n123)) {
279 return;
280 }
281 if (!checkVector(n012)) {
282 return;
283 }
284 if (!checkVector(n023)) {
285 return;
286 }
287
288 vec3_t xe13 = 0.5*(x[1] + x[3]);
289 vec3_t xe02 = 0.5*(x[0] + x[2]);
290
291 vec3_t n;
292
293 vec3_t xe13_p013 = proj->project(xe13, -1, true, n013);
294 n = proj->lastProjNormal();
295 double d13_013 = fabs((xe13_p013 - xe13)*n);
296
297 vec3_t xe13_p123 = proj->project(xe13, -1, true, n123);
298 n = proj->lastProjNormal();
299 double d13_123 = fabs((xe13_p123 - xe13)*n);
300
301 vec3_t xe02_p012 = proj->project(xe02, -1, true, n012);
302 n = proj->lastProjNormal();
303 double d02_012 = fabs((xe02_p012 - xe02)*n);
304
305 vec3_t xe02_p023 = proj->project(xe02, -1, true, n023);
306 n = proj->lastProjNormal();
307 double d02_023 = fabs((xe02_p023 - xe02)*n);
308
309 double L = max((x[1] - x[3]).abs(), (x[0] - x[2]).abs());
310
311 err1 = min(d13_013, d13_123)/L;
312 err2 = min(d02_012, d02_023)/L;
313
314 if (err1 < m_SurfErrorThreshold) {
315 err1 = 0;
316 }
317 if (err2 < m_SurfErrorThreshold) {
318 err2 = 0;
319 }
320 */
321 }
322
323 int SwapTriangles::swap()
324 {
325 int N_swaps = 0;
326 setAllSurfaceCells();
327 //computeAverageSurfaceError();
328 EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
329 EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type");
330 QVector<bool> marked(m_Grid->GetNumberOfCells(), false);
331 for (int i = 0; i < m_Part.getNumberOfCells(); ++i) {
332 vtkIdType id_cell = m_Part.globalCell(i);
333 if (!m_BoundaryCodes.contains(cell_code->GetValue(id_cell)) && m_Grid->GetCellType(id_cell) == VTK_TRIANGLE) { //if it is a selected triangle
334 if (!marked[id_cell] && !m_Swapped[id_cell]) {
335 for (int j = 0; j < 3; ++j) {
336 bool swap = false;
337 bool feature_improvement_swap = false;
338 stencil_t S = getStencil(id_cell, j);
339 if(S.id_cell.size() == 2 && S.sameBC) {
340 if (S.type_cell[1] == VTK_TRIANGLE) {
341 if(!isEdge(S.id_node[0], S.id_node[1]) ) {
342 if (!marked[S.id_cell[1]] && !m_Swapped[S.id_cell[1]]) {
343 QVector<vec3_t> x3(4);
344 vec2_t x[4];
345
346 m_Grid->GetPoint(S.id_node[0], x3[0].data());
347 m_Grid->GetPoint(S.p1, x3[1].data());
348 m_Grid->GetPoint(S.id_node[1], x3[2].data());
349 m_Grid->GetPoint(S.p2, x3[3].data());
350
351 vec3_t n1 = triNormal(x3[0], x3[1], x3[3]);
352 vec3_t n2 = triNormal(x3[1], x3[2], x3[3]);
353
354 bool force_swap = false;
355 if (m_SmallAreaSwap) {
356 double A1 = n1.abs();
357 double A2 = n2.abs();
358 if (isnan(A1) || isnan(A2)) {
359 force_swap = true;
360 } else {
361 force_swap = A1 < m_SmallAreaRatio*A2 || A2 < m_SmallAreaRatio*A1;
362 }
363 }
364 if (node_type->GetValue(S.p1) != EG_FEATURE_EDGE_VERTEX && node_type->GetValue(S.p2) != EG_FEATURE_EDGE_VERTEX) {
365
366 if (node_type->GetValue(S.id_node[0]) == EG_FEATURE_EDGE_VERTEX && node_type->GetValue(S.id_node[1]) == EG_FEATURE_EDGE_VERTEX) {
367 //swap = true;
368 }
369
370 // Delaunay
371 if (!swap) {
372 vec3_t n = n1 + n2;
373 n.normalise();
374 vec3_t ex = orthogonalVector(n);
375 vec3_t ey = ex.cross(n);
376 for (int k = 0; k < 4; ++k) {
377 x[k] = vec2_t(x3[k]*ex, x3[k]*ey);
378 }
379 vec2_t r1, r2, r3, u1, u2, u3;
380 r1 = 0.5*(x[0] + x[1]); u1 = turnLeft(x[1] - x[0]);
381 r2 = 0.5*(x[1] + x[2]); u2 = turnLeft(x[2] - x[1]);
382 r3 = 0.5*(x[1] + x[3]); u3 = turnLeft(x[3] - x[1]);
383 double k, l;
384 vec2_t xm1, xm2;
385 bool ok = true;
386 if (intersection(k, l, r1, u1, r3, u3)) {
387 xm1 = r1 + k*u1;
388 if (intersection(k, l, r2, u2, r3, u3)) {
389 xm2 = r2 + k*u2;
390 } else {
391 ok = false;
392 }
393 } else {
394 ok = false;
395 swap = true;
396 }
397 if (ok) {
398 if ((xm1 - x[2]).abs() < (xm1 - x[0]).abs()) {
399 swap = true;
400 }
401 if ((xm2 - x[0]).abs() < (xm2 - x[2]).abs()) {
402 swap = true;
403 }
404 }
405 }
406 } //end of if feature angle
407 } //end of if l_marked
408 } //end of if TestSwap
409 }
410 } //end of S valid
411
412 if (swap) {
413 if (testOrientation(S)) {
414 marked[S.id_cell[0]] = true;
415 marked[S.id_cell[1]] = true;
416 for (int k = 0; k < m_Part.n2cGSize(S.id_node[0]); ++k) {
417 vtkIdType id_neigh = m_Part.n2cGG(S.id_node[0], k);
418 marked[id_neigh] = true;
419 }
420 for (int k = 0; k < m_Part.n2cGSize(S.id_node[1]); ++k) {
421 vtkIdType id_neigh = m_Part.n2cGG(S.id_node[1], k);
422 marked[id_neigh] = true;
423 }
424 for (int k = 0; k < m_Part.n2cGSize(S.p1); ++k) {
425 vtkIdType id_neigh = m_Part.n2cGG(S.p1, k);
426 marked[id_neigh] = true;
427 }
428 for (int k = 0; k < m_Part.n2cGSize(S.p2); ++k) {
429 vtkIdType id_neigh = m_Part.n2cGG(S.p2, k);
430 marked[id_neigh] = true;
431 }
432 vtkIdType new_pts1[3], new_pts2[3];
433 new_pts1[0] = S.p1;
434 new_pts1[1] = S.id_node[1];
435 new_pts1[2] = S.id_node[0];
436 new_pts2[0] = S.id_node[1];
437 new_pts2[1] = S.p2;
438 new_pts2[2] = S.id_node[0];
439 vtkIdType old_pts1[3], old_pts2[3];
440 old_pts1[0] = S.id_node[0];
441 old_pts1[1] = S.p1;
442 old_pts1[2] = S.p2;
443 old_pts2[0] = S.id_node[1];
444 old_pts2[1] = S.p2;
445 old_pts2[2] = S.p1;
446 m_Grid->ReplaceCell(S.id_cell[0], 3, new_pts1);
447 m_Grid->ReplaceCell(S.id_cell[1], 3, new_pts2);
448 m_Swapped[S.id_cell[0]] = true;
449 m_Swapped[S.id_cell[1]] = true;
450 ++N_swaps;
451 break;
452 }
453 } //end of if swap
454 } //end of loop through sides
455 } //end of if marked
456 } //end of if selected triangle
457 } //end of loop through cells
458 return N_swaps;
459 }
460
461 void SwapTriangles::computeAverageSurfaceError()
462 {
463 EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
464 QList<double> errors;
465 EG_FORALL_NODES(id_node1, m_Grid) {
466 vec3_t n1 = m_Part.globalNormal(id_node1);
467 vec3_t x1;
468 m_Grid->GetPoint(id_node1, x1.data());
469 for (int i = 0; i < m_Part.n2nGSize(id_node1); ++i) {
470 vtkIdType id_node2 = m_Part.n2nGG(id_node1, i);
471 if (id_node1 < id_node2) {
472 QVector<vtkIdType> faces;
473 getEdgeCells(id_node1, id_node2, faces);
474 if (faces.size() == 2) {
475 int bc1 = cell_code->GetValue(faces[0]);
476 int bc2 = cell_code->GetValue(faces[1]);
477 if (bc1 == bc2) {
478 SurfaceProjection* proj = GuiMainWindow::pointer()->getSurfProj(bc1);
479 if (proj) {
480 vec3_t n2 = m_Part.globalNormal(id_node2);
481 vec3_t n = 0.5*(n1 + n2);
482 vec3_t x2;
483 m_Grid->GetPoint(id_node2, x2.data());
484 vec3_t x = 0.5*(x1 + x2);
485 vec3_t xp = proj->project(x, -1, true, n);
486 double err = (x - xp).abs()/(x1 - x2).abs();
487 errors.append(err);
488 }
489 }
490 }
491 }
492 }
493 }
494
495 m_AverageSurfaceError = Statistics::meanValue(errors);
496 m_SurfaceErrorDeviation = Statistics::standardDeviation(errors, m_AverageSurfaceError);
497 cout << "mean: " << m_AverageSurfaceError << " sigma: " << m_SurfaceErrorDeviation << endl;
498 }
499
500 void SwapTriangles::operate()
501 {
502 if (m_Verbose) cout << "swapping edges for surface triangles ..." << endl;
503 long int N_swaps = 100000000;
504 long int N_last_swaps = 100000001;
505 int loop = 1;
506 while ((N_swaps > 0) && (loop <= m_MaxNumLoops) && (N_swaps < N_last_swaps)) {
507 N_last_swaps = N_swaps;
508 if (m_Verbose) cout << " loop " << loop << "/" << m_MaxNumLoops << endl;
509 m_Swapped.fill(false, m_Grid->GetNumberOfCells());
510 N_swaps = 0;
511 int N;
512 int sub_loop = 0;
513 do {
514 ++sub_loop;
515 N = swap();
516 if (m_Verbose) cout << " sub-loop " << sub_loop << ": " << N << " swaps" << endl;
517 N_swaps += N;
518 } while (N > 0);
519 ++loop;
520 }
521 }
522
Clone of the official engrid repository