?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
on the way to a generic geometry interface
[engrid.git] / src / libengrid / swaptriangles.cpp
blob3f41df51ada8b93c752f7799f6a54df22ae16dd0
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->snapNode(x, -1, false);
121 return (x - xp).abs();
122 }
123
124 double SwapTriangles::edgeAngle(vec3_t x1, vec3_t x2, vec3_t x3, vec3_t x4)
125 {
126 static const double f13 = 1.0/3.0;
127 vec3_t xf1 = f13*(x1 + x2 + x4);
128 vec3_t xf2 = f13*(x2 + x3 + x4);
129 vec3_t xe = 0.5*(x2 + x4);
130 vec3_t xfe = 0.5*(xf1 + xf2);
131 vec3_t n1 = triNormal(x1, x2, x4);
132 vec3_t n2 = triNormal(x2, x3, x4);
133 vec3_t ve = xe - xfe;
134 double alpha = angle(n1, n2);
135 if ((n1 + n2)*ve < 0) {
136 return -alpha;
137 }
138 return alpha;
139 }
140
141 vtkIdType SwapTriangles::neighbourNode(vtkIdType id_node0, vtkIdType id_node1, vtkIdType id_node2)
142 {
143 for (int i = 0; i < m_Part.n2nGSize(id_node1); ++i) {
144 vtkIdType id_neigh = m_Part.n2nGG(id_node1, i);
145 if (id_neigh != id_node0) {
146 if (m_Part.hasNeighNode(id_node2, id_neigh)) {
147 return id_neigh;
148 }
149 }
150 }
151 EG_BUG;
152 return -1;
153 }
154
155 bool SwapTriangles::swapDueToSurfaceNoise(stencil_t S)
156 {
157 vtkIdType p1 = S.id_node[0];
158 vtkIdType p2 = S.p1;
159 vtkIdType p3 = S.id_node[1];
160 vtkIdType p4 = S.p2;
161 vtkIdType p12 = neighbourNode(p4, p1, p2);
162 vtkIdType p23 = neighbourNode(p4, p2, p3);
163 vtkIdType p34 = neighbourNode(p2, p3, p4);
164 vtkIdType p41 = neighbourNode(p2, p4, p1);
165 vec3_t x1, x2, x3, x4, x12, x23, x34, x41;
166 m_Grid->GetPoint(p1, x1.data());
167 m_Grid->GetPoint(p2, x2.data());
168 m_Grid->GetPoint(p3, x3.data());
169 m_Grid->GetPoint(p4, x4.data());
170 m_Grid->GetPoint(p12, x12.data());
171 m_Grid->GetPoint(p23, x23.data());
172 m_Grid->GetPoint(p34, x34.data());
173 m_Grid->GetPoint(p41, x41.data());
174
175 double noise1, noise2;
176 {
177 double alpha_min = M_PI;
178 double alpha_max = -M_PI;
179 double alpha;
180
181 alpha = edgeAngle(x4, x1, x12, x2);
182 alpha_max = max(alpha, alpha_max);
183 alpha_min = min(alpha, alpha_min);
184
185 alpha = edgeAngle(x4, x2, x23, x3);
186 alpha_max = max(alpha, alpha_max);
187 alpha_min = min(alpha, alpha_min);
188
189 alpha = edgeAngle(x2, x3, x34, x4);
190 alpha_max = max(alpha, alpha_max);
191 alpha_min = min(alpha, alpha_min);
192
193 alpha = edgeAngle(x2, x4, x41, x1);
194 alpha_max = max(alpha, alpha_max);
195 alpha_min = min(alpha, alpha_min);
196
197 alpha = edgeAngle(x1, x2, x3, x4);
198 alpha_max = max(alpha, alpha_max);
199 alpha_min = min(alpha, alpha_min);
200
201 noise1 = alpha_max - alpha_min;
202 }
203 {
204 double alpha_min = M_PI;
205 double alpha_max = -M_PI;
206 double alpha;
207
208 alpha = edgeAngle(x3, x1, x12, x2);
209 alpha_max = max(alpha, alpha_max);
210 alpha_min = min(alpha, alpha_min);
211
212 alpha = edgeAngle(x1, x2, x23, x3);
213 alpha_max = max(alpha, alpha_max);
214 alpha_min = min(alpha, alpha_min);
215
216 alpha = edgeAngle(x1, x3, x34, x4);
217 alpha_max = max(alpha, alpha_max);
218 alpha_min = min(alpha, alpha_min);
219
220 alpha = edgeAngle(x3, x4, x41, x1);
221 alpha_max = max(alpha, alpha_max);
222 alpha_min = min(alpha, alpha_min);
223
224 alpha = edgeAngle(x2, x3, x4, x1);
225 alpha_max = max(alpha, alpha_max);
226 alpha_min = min(alpha, alpha_min);
227
228 noise2 = alpha_max - alpha_min;
229 }
230
231 if (noise1 > m_SurfErrorThreshold && noise1 > noise2) {
232 return true;
233 }
234
235 return false;
236 }
237
238 void SwapTriangles::computeSurfaceErrors(const QVector<vec3_t> &x, int bc, double &err1, double &err2)
239 {
240 using namespace GeometryTools;
241 err1 = 0;
242 err2 = 1e10;
243 SurfaceProjection* proj = GuiMainWindow::pointer()->getSurfProj(bc);
244 if (!proj) {
245 return;
246 }
247
248 double d013 = computeSurfaceDistance(x[0], x[1], x[3], proj);
249 double d123 = computeSurfaceDistance(x[1], x[2], x[3], proj);
250 double d012 = computeSurfaceDistance(x[0], x[1], x[2], proj);
251 double d023 = computeSurfaceDistance(x[0], x[2], x[3], proj);
252
253 double L = max((x[1] - x[3]).abs(), (x[0] - x[2]).abs());
254 err1 = max(d013, d123)/L;
255 err2 = max(d012, d023)/L;
256 if (err1 < m_SurfErrorThreshold) {
257 err1 = 0;
258 }
259 if (err2 < m_SurfErrorThreshold) {
260 err2 = 0;
261 }
262
263 /*
264 // new version based on surface normals
265 vec3_t n013 = triNormal(x[0], x[1], x[3]);
266 vec3_t n123 = triNormal(x[1], x[2], x[3]);
267 vec3_t n012 = triNormal(x[0], x[1], x[2]);
268 vec3_t n023 = triNormal(x[0], x[2], x[3]);
269 n013.normalise();
270 n123.normalise();
271 n012.normalise();
272 n023.normalise();
273 if (!checkVector(n013)) {
274 return;
275 }
276 if (!checkVector(n123)) {
277 return;
278 }
279 if (!checkVector(n012)) {
280 return;
281 }
282 if (!checkVector(n023)) {
283 return;
284 }
285
286 vec3_t xe13 = 0.5*(x[1] + x[3]);
287 vec3_t xe02 = 0.5*(x[0] + x[2]);
288
289 vec3_t n;
290
291 vec3_t xe13_p013 = proj->project(xe13, -1, true, n013);
292 n = proj->lastProjNormal();
293 double d13_013 = fabs((xe13_p013 - xe13)*n);
294
295 vec3_t xe13_p123 = proj->project(xe13, -1, true, n123);
296 n = proj->lastProjNormal();
297 double d13_123 = fabs((xe13_p123 - xe13)*n);
298
299 vec3_t xe02_p012 = proj->project(xe02, -1, true, n012);
300 n = proj->lastProjNormal();
301 double d02_012 = fabs((xe02_p012 - xe02)*n);
302
303 vec3_t xe02_p023 = proj->project(xe02, -1, true, n023);
304 n = proj->lastProjNormal();
305 double d02_023 = fabs((xe02_p023 - xe02)*n);
306
307 double L = max((x[1] - x[3]).abs(), (x[0] - x[2]).abs());
308
309 err1 = min(d13_013, d13_123)/L;
310 err2 = min(d02_012, d02_023)/L;
311
312 if (err1 < m_SurfErrorThreshold) {
313 err1 = 0;
314 }
315 if (err2 < m_SurfErrorThreshold) {
316 err2 = 0;
317 }
318 */
319 }
320
321 int SwapTriangles::swap()
322 {
323 int N_swaps = 0;
324 setAllSurfaceCells();
325 //computeAverageSurfaceError();
326 EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
327 EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type");
328 QVector<bool> marked(m_Grid->GetNumberOfCells(), false);
329 for (int i = 0; i < m_Part.getNumberOfCells(); ++i) {
330 vtkIdType id_cell = m_Part.globalCell(i);
331 if (!m_BoundaryCodes.contains(cell_code->GetValue(id_cell)) && m_Grid->GetCellType(id_cell) == VTK_TRIANGLE) { //if it is a selected triangle
332 if (!marked[id_cell] && !m_Swapped[id_cell]) {
333 for (int j = 0; j < 3; ++j) {
334 bool swap = false;
335 bool feature_improvement_swap = false;
336 stencil_t S = getStencil(id_cell, j);
337 if(S.id_cell.size() == 2 && S.sameBC) {
338 if (S.type_cell[1] == VTK_TRIANGLE) {
339 if(!isEdge(S.id_node[0], S.id_node[1]) ) {
340 if (!marked[S.id_cell[1]] && !m_Swapped[S.id_cell[1]]) {
341 QVector<vec3_t> x3(4);
342 vec2_t x[4];
343
344 m_Grid->GetPoint(S.id_node[0], x3[0].data());
345 m_Grid->GetPoint(S.p1, x3[1].data());
346 m_Grid->GetPoint(S.id_node[1], x3[2].data());
347 m_Grid->GetPoint(S.p2, x3[3].data());
348
349 vec3_t n1 = triNormal(x3[0], x3[1], x3[3]);
350 vec3_t n2 = triNormal(x3[1], x3[2], x3[3]);
351
352 bool force_swap = false;
353 if (m_SmallAreaSwap) {
354 double A1 = n1.abs();
355 double A2 = n2.abs();
356 if (isnan(A1) || isnan(A2)) {
357 force_swap = true;
358 } else {
359 force_swap = A1 < m_SmallAreaRatio*A2 || A2 < m_SmallAreaRatio*A1;
360 }
361 }
362 if (!isFeatureNode(S.p1) && !isFeatureNode(S.p2)) {
363
364 if (isFeatureNode(S.id_node[0]) && isFeatureNode(S.id_node[1])) {
365 //swap = true;
366 }
367
368 // Delaunay
369 if (!swap) {
370 vec3_t n = n1 + n2;
371 n.normalise();
372 vec3_t ex = orthogonalVector(n);
373 vec3_t ey = ex.cross(n);
374 for (int k = 0; k < 4; ++k) {
375 x[k] = vec2_t(x3[k]*ex, x3[k]*ey);
376 }
377 vec2_t r1, r2, r3, u1, u2, u3;
378 r1 = 0.5*(x[0] + x[1]); u1 = turnLeft(x[1] - x[0]);
379 r2 = 0.5*(x[1] + x[2]); u2 = turnLeft(x[2] - x[1]);
380 r3 = 0.5*(x[1] + x[3]); u3 = turnLeft(x[3] - x[1]);
381 double k, l;
382 vec2_t xm1, xm2;
383 bool ok = true;
384 if (intersection(k, l, r1, u1, r3, u3)) {
385 xm1 = r1 + k*u1;
386 if (intersection(k, l, r2, u2, r3, u3)) {
387 xm2 = r2 + k*u2;
388 } else {
389 ok = false;
390 }
391 } else {
392 ok = false;
393 swap = true;
394 }
395 if (ok) {
396 if ((xm1 - x[2]).abs() < (xm1 - x[0]).abs()) {
397 swap = true;
398 }
399 if ((xm2 - x[0]).abs() < (xm2 - x[2]).abs()) {
400 swap = true;
401 }
402 }
403 }
404 } //end of if feature angle
405 } //end of if l_marked
406 } //end of if TestSwap
407 }
408 } //end of S valid
409
410 if (swap) {
411 if (testOrientation(S)) {
412 marked[S.id_cell[0]] = true;
413 marked[S.id_cell[1]] = true;
414 for (int k = 0; k < m_Part.n2cGSize(S.id_node[0]); ++k) {
415 vtkIdType id_neigh = m_Part.n2cGG(S.id_node[0], k);
416 marked[id_neigh] = true;
417 }
418 for (int k = 0; k < m_Part.n2cGSize(S.id_node[1]); ++k) {
419 vtkIdType id_neigh = m_Part.n2cGG(S.id_node[1], k);
420 marked[id_neigh] = true;
421 }
422 for (int k = 0; k < m_Part.n2cGSize(S.p1); ++k) {
423 vtkIdType id_neigh = m_Part.n2cGG(S.p1, k);
424 marked[id_neigh] = true;
425 }
426 for (int k = 0; k < m_Part.n2cGSize(S.p2); ++k) {
427 vtkIdType id_neigh = m_Part.n2cGG(S.p2, k);
428 marked[id_neigh] = true;
429 }
430 vtkIdType new_pts1[3], new_pts2[3];
431 new_pts1[0] = S.p1;
432 new_pts1[1] = S.id_node[1];
433 new_pts1[2] = S.id_node[0];
434 new_pts2[0] = S.id_node[1];
435 new_pts2[1] = S.p2;
436 new_pts2[2] = S.id_node[0];
437 vtkIdType old_pts1[3], old_pts2[3];
438 old_pts1[0] = S.id_node[0];
439 old_pts1[1] = S.p1;
440 old_pts1[2] = S.p2;
441 old_pts2[0] = S.id_node[1];
442 old_pts2[1] = S.p2;
443 old_pts2[2] = S.p1;
444 m_Grid->ReplaceCell(S.id_cell[0], 3, new_pts1);
445 m_Grid->ReplaceCell(S.id_cell[1], 3, new_pts2);
446 m_Swapped[S.id_cell[0]] = true;
447 m_Swapped[S.id_cell[1]] = true;
448 ++N_swaps;
449 break;
450 }
451 } //end of if swap
452 } //end of loop through sides
453 } //end of if marked
454 } //end of if selected triangle
455 } //end of loop through cells
456 return N_swaps;
457 }
458
459 void SwapTriangles::computeAverageSurfaceError()
460 {
461 EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
462 QList<double> errors;
463 EG_FORALL_NODES(id_node1, m_Grid) {
464 vec3_t n1 = m_Part.globalNormal(id_node1);
465 vec3_t x1;
466 m_Grid->GetPoint(id_node1, x1.data());
467 for (int i = 0; i < m_Part.n2nGSize(id_node1); ++i) {
468 vtkIdType id_node2 = m_Part.n2nGG(id_node1, i);
469 if (id_node1 < id_node2) {
470 QVector<vtkIdType> faces;
471 getEdgeCells(id_node1, id_node2, faces);
472 if (faces.size() == 2) {
473 int bc1 = cell_code->GetValue(faces[0]);
474 int bc2 = cell_code->GetValue(faces[1]);
475 if (bc1 == bc2) {
476 SurfaceProjection* proj = GuiMainWindow::pointer()->getSurfProj(bc1);
477 if (proj) {
478 vec3_t n2 = m_Part.globalNormal(id_node2);
479 vec3_t n = 0.5*(n1 + n2);
480 vec3_t x2;
481 m_Grid->GetPoint(id_node2, x2.data());
482 vec3_t x = 0.5*(x1 + x2);
483 vec3_t xp = proj->snapNode(x, -1, false);
484 double err = (x - xp).abs()/(x1 - x2).abs();
485 errors.append(err);
486 }
487 }
488 }
489 }
490 }
491 }
492
493 m_AverageSurfaceError = Statistics::meanValue(errors);
494 m_SurfaceErrorDeviation = Statistics::standardDeviation(errors, m_AverageSurfaceError);
495 cout << "mean: " << m_AverageSurfaceError << " sigma: " << m_SurfaceErrorDeviation << endl;
496 }
497
498 void SwapTriangles::operate()
499 {
500 if (m_Verbose) cout << "swapping edges for surface triangles ..." << endl;
501 long int N_swaps = 100000000;
502 long int N_last_swaps = 100000001;
503 int loop = 1;
504 while ((N_swaps > 0) && (loop <= m_MaxNumLoops) && (N_swaps < N_last_swaps)) {
505 N_last_swaps = N_swaps;
506 if (m_Verbose) cout << " loop " << loop << "/" << m_MaxNumLoops << endl;
507 m_Swapped.fill(false, m_Grid->GetNumberOfCells());
508 N_swaps = 0;
509 int N;
510 int sub_loop = 0;
511 do {
512 ++sub_loop;
513 N = swap();
514 if (m_Verbose) cout << " sub-loop " << sub_loop << ": " << N << " swaps" << endl;
515 N_swaps += N;
516 } while (N > 0);
517 ++loop;
518 }
519 }
520
Clone of the official engrid repository