?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
fixed edge display for volume cells
[engrid-github.git] / src / libengrid / boundarylayeroperation.cpp
blob746f3c6aa334d39ebe489b0e809ac30bfdb4cb4b
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 //
22 #include "boundarylayeroperation.h"
23 #include "optimisenormalvector.h"
24 #include "guimainwindow.h"
25 #include "pointfinder.h"
26 #include "cgaltricadinterface.h"
27 #include "geometrysmoother.h"
28 #include "vtkEgPolyDataToUnstructuredGridFilter.h"
29
30 #include <vtkDataSetSurfaceFilter.h>
31 #include <vtkLoopSubdivisionFilter.h>
32 #include <vtkLinearSubdivisionFilter.h>
33 #include <vtkButterflySubdivisionFilter.h>
34 #include <vtkSmoothPolyDataFilter.h>
35 #include <vtkWindowedSincPolyDataFilter.h>
36
37 void BoundaryLayerOperation::readSettings()
38 {
39 EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
40 QString buffer = GuiMainWindow::pointer()->getXmlSection("engrid/blayer/settings");
41 if(!buffer.isEmpty()) {
42 QTextStream in(&buffer, QIODevice::ReadOnly);
43 int num_bcs;
44 in >> num_bcs;
45 QVector<bool> use_bc(num_bcs);
46 int N = 0;
47 for (int i = 0; i < num_bcs; ++i) {
48 int state;
49 in >> state;
50 use_bc[i] = bool(state);
51 if (use_bc[i]) {
52 ++N;
53 }
54 }
55 m_BoundaryLayerCodes.resize(N);
56 QVector<int> bcs;
57 mainWindow()->getAllBoundaryCodes(bcs);
58 int j = 0;
59 for (int i = 0; i < bcs.size(); ++i) {
60 if (use_bc[i]) {
61 m_BoundaryLayerCodes[j++] = bcs[i];
62 }
63 }
64 m_LayerAdjacentBoundaryCodes.clear();
65 for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
66 if (isSurface(id_cell, m_Grid)) {
67 if (m_BoundaryLayerCodes.contains(cell_code->GetValue(id_cell))) {
68 for (int i = 0; i < m_Part.c2cGSize(id_cell); ++i) {
69 vtkIdType id_neigh = m_Part.c2cGG(id_cell, i);
70 if (!m_BoundaryLayerCodes.contains(cell_code->GetValue(id_neigh))) {
71 m_LayerAdjacentBoundaryCodes.insert(cell_code->GetValue(id_neigh));
72 }
73 }
74 }
75 }
76 }
77 in >> m_FeatureAngle;
78 m_FeatureAngle = deg2rad(m_FeatureAngle);
79 in >> m_StretchingRatio;
80 in >> m_FarfieldRatio;
81 in >> m_NumBoundaryLayerVectorRelaxations;
82 in >> m_NumBoundaryLayerHeightRelaxations;
83 in >> m_ShellPassBand;
84 in >> m_FaceSizeLowerLimit;
85 in >> m_FaceSizeUpperLimit;
86 in >> m_FaceAngleLimit;
87 m_FaceAngleLimit = deg2rad(m_FaceAngleLimit);
88 in >> m_MaxHeightInGaps;
89 in >> m_RadarAngle;
90 int use_grouping;
91 in >> use_grouping;
92 m_UseGrouping = use_grouping;
93 in >> m_GroupingAngle;
94 m_GroupingAngle = deg2rad(m_GroupingAngle);
95 }
96 m_ELSManagerBLayer.clear();
97 m_ELSManagerBLayer.readBoundaryLayerRules(m_Grid);
98 m_ELSManagerSurface.clear();
99 m_ELSManagerSurface.read();
100 m_ELSManagerSurface.readRules(m_Grid);
101 }
102
103 void BoundaryLayerOperation::computeBoundaryLayerVectors()
104 {
105 EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
106 m_BoundaryLayerVectors.fill(vec3_t(0,0,0), m_Grid->GetNumberOfPoints());
107 QVector<int> num_bcs(m_Grid->GetNumberOfPoints());
108 QVector<OptimiseNormalVector> n_opt(m_Grid->GetNumberOfPoints(), OptimiseNormalVector(m_UseGrouping, m_GroupingAngle));
109 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
110 QSet<int> bcs;
111 for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
112 vtkIdType id_cell = m_Part.n2cGG(id_node, i);
113 if (isSurface(id_cell, m_Grid)) {
114 int bc = cell_code->GetValue(id_cell);
115 if (m_BoundaryLayerCodes.contains(bc)) {
116 bcs.insert(bc);
117 }
118 }
119 }
120 num_bcs[id_node] = bcs.size();
121 QVector<vec3_t> normal(num_bcs[id_node], vec3_t(0,0,0));
122 QMap<int,int> bcmap;
123 int i_bc = 0;
124 foreach (int bc, bcs) {
125 bcmap[bc] = i_bc;
126 ++i_bc;
127 }
128 for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
129 vtkIdType id_cell = m_Part.n2cGG(id_node, i);
130 if (isSurface(id_cell, m_Grid)) {
131 int bc = cell_code->GetValue(id_cell);
132 vtkIdType N_pts, *pts;
133 m_Grid->GetCellPoints(id_cell, N_pts, pts);
134 vec3_t a, b, c;
135 for (int j = 0; j < N_pts; ++j) {
136 if (pts[j] == id_node) {
137 m_Grid->GetPoint(pts[j], a.data());
138 if (j > 0) {
139 m_Grid->GetPoint(pts[j-1], b.data());
140 } else {
141 m_Grid->GetPoint(pts[N_pts-1], b.data());
142 }
143 if (j < N_pts - 1) {
144 m_Grid->GetPoint(pts[j+1], c.data());
145 } else {
146 m_Grid->GetPoint(pts[0], c.data());
147 }
148 }
149 }
150 vec3_t u = b - a;
151 vec3_t v = c - a;
152 double alpha = GeometryTools::angle(u, v);
153 vec3_t n = u.cross(v);
154 n.normalise();
155 if (m_BoundaryLayerCodes.contains(bc)) {
156 normal[bcmap[bc]] += alpha*n;
157 n_opt[id_node].addFace(n);
158 } else {
159 n_opt[id_node].addConstraint(n);
160 }
161 }
162 }
163 for (int i = 0; i < num_bcs[id_node]; ++i) {
164 normal[i].normalise();
165 }
166 if (num_bcs[id_node] > 0) {
167 if (num_bcs[id_node] > 1) {
168 if (num_bcs[id_node] == 3) {
169 for (int i = 0; i < num_bcs[id_node]; ++i) {
170 for (int j = i + 1; j < num_bcs[id_node]; ++j) {
171 vec3_t n = normal[i] + normal[j];
172 n.normalise();
173 m_BoundaryLayerVectors[id_node] += n;
174 }
175 }
176 } else {
177 for (int i = 0; i < num_bcs[id_node]; ++i) {
178 m_BoundaryLayerVectors[id_node] += normal[i];
179 }
180 }
181 } else {
182 m_BoundaryLayerVectors[id_node] = normal[0];
183 }
184 m_BoundaryLayerVectors[id_node].normalise();
185 m_BoundaryLayerVectors[id_node] = n_opt[id_node](m_BoundaryLayerVectors[id_node]);
186 m_BoundaryLayerVectors[id_node].normalise();
187 }
188 if (!checkVector(m_BoundaryLayerVectors[id_node])) {
189 EG_ERR_RETURN("invalid layer vector computed");
190 }
191 }
192
193 m_BoundaryLayerNode.fill(false, m_Grid->GetNumberOfPoints());
194 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
195 for (int i = 0; i < m_Part.n2bcGSize(id_node); ++i) {
196 if (m_BoundaryLayerCodes.contains(m_Part.n2bcG(id_node, i))) {
197 m_BoundaryLayerNode[id_node] = true;
198 }
199 }
200 }
201
202 computeNodeTypes();
203 smoothBoundaryLayerVectors(m_NumBoundaryLayerVectorRelaxations);
204
205 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
206 if (m_BoundaryLayerVectors[id_node].abs() < 0.1) {
207 vec3_t n(0,0,0);
208 int N = 0;
209 for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
210 vtkIdType id_cell = m_Part.n2cGG(id_node, i);
211 if (isSurface(id_cell, m_Grid)) {
212 n += GeometryTools::cellNormal(m_Grid, id_cell);
213 ++N;
214 }
215 }
216 if (N) {
217 n.normalise();
218 m_BoundaryLayerVectors[id_node] = n;
219 }
220 if (num_bcs[id_node] > 1) {
221 m_BoundaryLayerVectors[id_node] = n_opt[id_node](m_BoundaryLayerVectors[id_node]);
222 }
223 if (!checkVector(m_BoundaryLayerVectors[id_node])) {
224 EG_ERR_RETURN("invalid layer vector computed");
225 }
226 }
227 }
228
229 computeHeights();
230 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
231 if (m_BoundaryLayerNode[id_node]) {
232 m_BoundaryLayerVectors[id_node] *= m_Height[id_node];
233 if (!checkVector(m_BoundaryLayerVectors[id_node])) {
234 EG_ERR_RETURN("invalid layer vector computed");
235 }
236 }
237 }
238
239 }
240
241 void BoundaryLayerOperation::addToSnapPoints(vtkIdType id_node, vtkIdType id_snap)
242 {
243 if (m_NodeTypes[id_node] == NormalNode) {
244 m_SnapPoints[id_node].insert(id_snap);
245 } else if (m_NodeTypes[id_node] == EdgeNode && m_NodeTypes[id_snap] != NormalNode) {
246 m_SnapPoints[id_node].insert(id_snap);
247 }
248 }
249
250 void BoundaryLayerOperation::computeNodeTypes()
251 {
252 m_NodeTypes.fill(NormalNode, m_Grid->GetNumberOfPoints());
253 m_SnapPoints.resize(m_Grid->GetNumberOfPoints());
254 QVector<int> num_feature_edges(m_Grid->GetNumberOfPoints(), 0);
255
256 for (vtkIdType id_cell1 = 0; id_cell1 < m_Grid->GetNumberOfCells(); ++id_cell1) {
257 if (isSurface(id_cell1, m_Grid)) {
258 vec3_t n1 = cellNormal(m_Grid, id_cell1);
259 for (int i = 0; i < m_Part.c2cGSize(id_cell1); ++i) {
260 vtkIdType id_cell2 = m_Part.c2cGG(id_cell1, i);
261 if (id_cell2 > id_cell1) {
262 if (!isSurface(id_cell2, m_Grid)) {
263 EG_BUG;
264 }
265 vec3_t n2 = cellNormal(m_Grid, id_cell2);
266 if (angle(n1, n2) >= m_FeatureAngle) {
267 QVector<vtkIdType> nodes;
268 m_Part.getCommonNodes(id_cell1, id_cell2, nodes);
269 foreach (vtkIdType id_node, nodes) {
270 ++num_feature_edges[id_node];
271 }
272 }
273 }
274 }
275 }
276 }
277 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
278 if (num_feature_edges[id_node] > 2) m_NodeTypes[id_node] = CornerNode;
279 else if (num_feature_edges[id_node] > 1) m_NodeTypes[id_node] = EdgeNode;
280 else m_NodeTypes[id_node] = NormalNode;
281 }
282 for (vtkIdType id_cell1 = 0; id_cell1 < m_Grid->GetNumberOfCells(); ++id_cell1) {
283 if (isSurface(id_cell1, m_Grid)) {
284 vec3_t n1 = cellNormal(m_Grid, id_cell1);
285 for (int i = 0; i < m_Part.c2cGSize(id_cell1); ++i) {
286 vtkIdType id_cell2 = m_Part.c2cGG(id_cell1, i);
287 if (id_cell2 > id_cell1) {
288 if (!isSurface(id_cell2, m_Grid)) {
289 EG_BUG;
290 }
291 vec3_t n2 = cellNormal(m_Grid, id_cell2);
292 QVector<vtkIdType> nodes;
293 m_Part.getCommonNodes(id_cell1, id_cell2, nodes);
294 if (nodes.size() != 2) {
295 EG_BUG;
296 }
297 if (angle(n1, n2) >= m_FeatureAngle) {
298 addToSnapPoints(nodes[0], nodes[1]);
299 addToSnapPoints(nodes[1], nodes[0]);
300 } else {
301 if (m_NodeTypes[nodes[0]] == NormalNode) {
302 m_SnapPoints[nodes[0]].insert(nodes[1]);
303 }
304 if (m_NodeTypes[nodes[1]] == NormalNode) {
305 m_SnapPoints[nodes[1]].insert(nodes[0]);
306 }
307 }
308 }
309 }
310 }
311 }
312 }
313
314 void BoundaryLayerOperation::correctBoundaryLayerVectors()
315 {
316 EG_VTKDCC(vtkIntArray, bc, m_Grid, "cell_code");
317 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
318 if (m_BoundaryLayerVectors[id_node].abs() > 0.1) {
319 for (int iter = 0; iter < 20; ++iter) {
320 for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
321 vtkIdType id_cell = m_Part.n2cGG(id_node,i);
322 if (isSurface(id_cell, m_Grid)) {
323 if (!m_BoundaryLayerCodes.contains(bc->GetValue(id_cell))) {
324 vec3_t v = m_BoundaryLayerVectors[id_node];
325 vec3_t n = GeometryTools::cellNormal(m_Grid, id_cell);
326 n.normalise();
327 v -= (n*m_BoundaryLayerVectors[id_node])*n;
328 v.normalise();
329 m_BoundaryLayerVectors[id_node] = v;
330 }
331 }
332 }
333 }
334 }
335 }
336 }
337
338 void BoundaryLayerOperation::smoothBoundaryLayerVectors(int n_iter, double w_iso, double w_dir, QVector<bool> *node_fixed)
339 {
340 EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
341
342 for (int iter = 0; iter < n_iter; ++iter) {
343 QVector<vec3_t> v_new = m_BoundaryLayerVectors;
344 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
345 bool fixed = false;
346 if (node_fixed) {
347 fixed = (*node_fixed)[id_node];
348 }
349 if (m_BoundaryLayerNode[id_node] && !fixed) {
350 v_new[id_node] = vec3_t(0,0,0);
351 if (m_SnapPoints[id_node].size() > 0) {
352
353 // check for edge between corners
354 bool edge_between_corners = false;
355 if (m_NodeTypes[id_node] == EdgeNode) {
356 edge_between_corners = true;
357 foreach (vtkIdType id_snap, m_SnapPoints[id_node]) {
358 if (m_NodeTypes[id_snap] != CornerNode) {
359 edge_between_corners = false;
360 break;
361 }
362 }
363 }
364
365 bool smooth_node = !edge_between_corners;
366 if (!smooth_node && m_SnapPoints[id_node].size() > 0) {
367
368 // compute smoothed normal
369 vec3_t n_smooth(0,0,0);
370 foreach (vtkIdType id_snap, m_SnapPoints[id_node]) {
371 n_smooth += m_BoundaryLayerVectors[id_snap];
372 }
373 n_smooth.normalise();
374
375 // check if it has not been projected onto the plane of an adjacent face
376 double h_min = 1.0;
377 for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
378 vtkIdType id_face = m_Part.n2cGG(id_node, i);
379 if (m_BoundaryLayerCodes.contains(cell_code->GetValue(id_face))) {
380 vec3_t n_face = cellNormal(m_Grid, id_face);
381 n_face.normalise();
382 n_face *= -1;
383 h_min = min(h_min, n_face*n_smooth);
384 }
385 }
386 if (h_min > 0.2) {
387 smooth_node = true;
388 }
389
390 }
391
392 if (smooth_node) {
393 v_new[id_node] = vec3_t(0,0,0);
394 vec3_t x_node;
395 m_Grid->GetPoint(id_node, x_node.data());
396 double w_total = 0.0;
397 foreach (vtkIdType id_snap, m_SnapPoints[id_node]) {
398 vec3_t x_snap;
399 m_Grid->GetPoint(id_snap, x_snap.data());
400 QSet<vtkIdType> faces;
401 getEdgeCells(id_node, id_snap, faces);
402 vec3_t n_edge(0,0,0);
403 foreach (vtkIdType id_face, faces) {
404 n_edge += cellNormal(m_Grid, id_face);
405 }
406 n_edge.normalise();
407 vec3_t u = m_BoundaryLayerVectors[id_snap] - (m_BoundaryLayerVectors[id_snap]*n_edge)*m_BoundaryLayerVectors[id_snap];
408 vec3_t v = x_node - x_snap;
409 v.normalise();
410 double w = w_iso + w_dir*(u*v);
411 w_total += w;
412 v_new[id_node] += w*m_BoundaryLayerVectors[id_snap];
413 }
414 //v_new[id_node].normalise();
415 v_new[id_node] *= 1.0/w_total;
416 } else {
417 v_new[id_node] = m_BoundaryLayerVectors[id_node];
418 }
419 } else {
420 v_new[id_node] = m_BoundaryLayerVectors[id_node];
421 }
422 if (v_new[id_node].abs() < 0.1) {
423 EG_BUG;
424 }
425 }
426 }
427 m_BoundaryLayerVectors = v_new;
428 correctBoundaryLayerVectors();
429 }
430 }
431
432 void BoundaryLayerOperation::writeBoundaryLayerVectors(QString file_name, int counter)
433 {
434 if (counter >= 0) {
435 QString counter_txt;
436 counter_txt.setNum(counter);
437 counter_txt = counter_txt.rightJustified(3, '0');
438 file_name += "_" + counter_txt;
439 }
440 MeshPartition wall_part(m_Grid);
441 wall_part.setBCs(m_BoundaryLayerCodes);
442 file_name = GuiMainWindow::pointer()->getCwd() + "/" + file_name + ".vtk";
443 QFile file(file_name);
444 file.open(QIODevice::WriteOnly);
445 QTextStream f(&file);
446 f << "# vtk DataFile Version 2.0\n";
447 f << "m_BoundaryLayerVectors\n";
448 f << "ASCII\n";
449 f << "DATASET UNSTRUCTURED_GRID\n";
450 f << "POINTS " << wall_part.getNumberOfNodes() << " float\n";
451 for (int i = 0; i < wall_part.getNumberOfNodes(); ++i) {
452 vec3_t x;
453 vtkIdType id_node = wall_part.globalNode(i);
454 m_Grid->GetPoint(id_node, x.data());
455 f << x[0] << " " << x[1] << " " << x[2] << "\n";
456 }
457 f << "CELLS " << wall_part.getNumberOfCells() << " ";
458 int N = 0;
459 for (int i = 0; i < wall_part.getNumberOfCells(); ++ i) {
460 vtkIdType id_cell = wall_part.globalCell(i);
461 vtkIdType N_pts, *pts;
462 m_Grid->GetCellPoints(id_cell, N_pts, pts);
463 N += 1 + N_pts;
464 }
465 f << N << "\n";
466 for (int i = 0; i < wall_part.getNumberOfCells(); ++ i) {
467 vtkIdType id_cell = wall_part.globalCell(i);
468 vtkIdType N_pts, *pts;
469 m_Grid->GetCellPoints(id_cell, N_pts, pts);
470 f << N_pts;
471 for (int j = 0; j < N_pts; ++j) {
472 f << " " << wall_part.localNode(pts[j]);
473 }
474 f << "\n";
475 }
476
477 f << "CELL_TYPES " << wall_part.getNumberOfCells() << "\n";
478 for (int i = 0; i < wall_part.getNumberOfCells(); ++ i) {
479 vtkIdType id_cell = wall_part.globalCell(i);
480 vtkIdType type_cell = m_Grid->GetCellType(id_cell);
481 f << type_cell << "\n";
482 }
483 f << "POINT_DATA " << wall_part.getNumberOfNodes() << "\n";
484
485 f << "VECTORS N float\n";
486 for (int i = 0; i < wall_part.getNumberOfNodes(); ++i) {
487 vtkIdType id_node = wall_part.globalNode(i);
488 f << m_BoundaryLayerVectors[id_node][0] << " ";
489 f << m_BoundaryLayerVectors[id_node][1] << " ";
490 f << m_BoundaryLayerVectors[id_node][2] << "\n";
491 }
492
493 f << "SCALARS node_type int\n";
494 f << "LOOKUP_TABLE default\n";
495 for (int i = 0; i < wall_part.getNumberOfNodes(); ++i) {
496 vtkIdType id_node = wall_part.globalNode(i);
497 f << m_NodeTypes[id_node] << "\n";
498 }
499 }
500
501 void BoundaryLayerOperation::computeDesiredHeights()
502 {
503 EG_VTKDCN(vtkDoubleArray, cl, m_Grid, "node_meshdensity_desired");
504
505 // first pass (intial height)
506 m_Height.fill(0, m_Grid->GetNumberOfPoints());
507 int k = 1;
508 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
509 if (m_BoundaryLayerNode[id_node]) {
510 vec3_t x;
511 m_Grid->GetPoint(id_node, x.data());
512 double h0 = m_ELSManagerBLayer.minEdgeLength(x);
513 double h1 = cl->GetValue(id_node)*m_FarfieldRatio;
514 k = max(k, int(logarithm(m_StretchingRatio, h1/h0)));
515 }
516 }
517 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
518 if (m_BoundaryLayerNode[id_node]) {
519 vec3_t x;
520 m_Grid->GetPoint(id_node, x.data());
521 double h0 = m_ELSManagerBLayer.minEdgeLength(x);
522 double h1 = cl->GetValue(id_node)*m_FarfieldRatio;
523 double s = pow(h1/h0, 1.0/k);
524 double H = h0*(1 - pow(s, k))/(1 - s);
525 if (!checkReal(H)) {
526 EG_ERR_RETURN("floating point error while computing heights");
527 }
528 if (H < 0) {
529 EG_ERR_RETURN("negative height computed");
530 }
531 if (H > 1000*h1) {
532 cout << H << ", " << h1 << endl;
533 EG_ERR_RETURN("unrealistically large height computed");
534 }
535 if (H < 1e-3*h0) {
536 EG_ERR_RETURN("unrealistically small height computed");
537 }
538 if (h1 < h0) {
539 QString h0_txt, h1_txt, id_txt;
540 h0_txt.setNum(h0);
541 h1_txt.setNum(h1);
542 id_txt.setNum(id_node);
543 EG_ERR_RETURN("h1 < h0 (" + h1_txt + " < " + h0_txt + ", for node " + id_txt + ")");
544 }
545 m_Height[id_node] = H;
546 }
547 }
548
549 m_NumLayers = k;
550
551 // correct with angle between face normal and propagation direction (node normals)
552 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
553 if (m_BoundaryLayerNode[id_node]) {
554 int N = 0;
555 double scale = 0;
556 for (int j = 0; j < m_Part.n2cGSize(id_node); ++j) {
557 vtkIdType id_cell = m_Part.n2cGG(id_node, j);
558 if (isSurface(id_cell, m_Grid)) {
559 scale += m_BoundaryLayerVectors[id_node]*cellNormal(m_Grid, id_cell).normalise();
560 }
561 }
562 if (N > 0) {
563 scale /= N;
564 m_Height[id_node] /= scale;
565 }
566 }
567 }
568 }
569
570 bool BoundaryLayerOperation::faceFine(vtkIdType id_face, double scale)
571 {
572 EG_GET_CELL(id_face, m_Grid);
573 if (type_cell != VTK_TRIANGLE) {
574 EG_BUG;
575 }
576 QVector<vec3_t> x1(num_pts);
577 for (vtkIdType i = 0; i < num_pts; ++i) {
578 m_Grid->GetPoint(pts[i], x1[i].data());
579 }
580 vec3_t n1 = triNormal(x1[0], x1[1], x1[2]);
581 QVector<vec3_t> x2(num_pts);
582 for (vtkIdType i = 0; i < num_pts; ++i) {
583 x2[i] = x1[i] + scale*m_Height[pts[i]]*m_BoundaryLayerVectors[pts[i]];
584 }
585 vec3_t n2 = triNormal(x2[0], x2[1], x2[2]);
586 double A1 = n1.abs();
587 double A2 = n2.abs();
588 if (A2/A1 >= m_FaceSizeLowerLimit && A2/A1 <= m_FaceSizeUpperLimit){
589 if (angle(n1, n2) < m_FaceAngleLimit) {
590 return true;
591 }
592 }
593 return false;
594 }
595
596 void BoundaryLayerOperation::computeHeights()
597 {
598 EG_VTKDCC(vtkIntArray, bc, m_Grid, "cell_code");
599
600 computeDesiredHeights();
601 cout << "initial boundary layer heights computed" << endl;
602
603 // avoid collisions
604 limitHeights(1.0);
605
606 // limit face size and angle difference
607 //limitSizeAndAngleErrors();
608
609 // smoothing
610 {
611 QVector<double> h_safe = m_Height;
612 for (int iter = 0; iter < m_NumBoundaryLayerHeightRelaxations; ++iter) {
613 QVector<double> h_new = m_Height;
614 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
615 if (m_BoundaryLayerNode[id_node]) {
616 int count = 0;
617 h_new[id_node] = 0;
618 for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
619 vtkIdType id_neigh = m_Part.n2nGG(id_node, i);
620 if (m_BoundaryLayerNode[id_neigh]) {
621 ++count;
622 h_new[id_node] += min(h_safe[id_node], m_Height[id_neigh]);
623 }
624 }
625 if (count == 0) {
626 EG_BUG;
627 }
628 h_new[id_node] /= count;
629 }
630 }
631 m_Height = h_new;
632 }
633 }
634
635 // The mesh smoothing methods start here
636 // General variables kind of useful for push out, smoothing, etc
637 // Move to somewhere else later?
638 QVector<bool> on_boundary(m_Grid->GetNumberOfPoints(), false);
639 QVector<bool> is_convex(m_Grid->GetNumberOfPoints(), false);
640 QVector<vec3_t> grid_pnts(m_Grid->GetNumberOfPoints(), vec3_t(0,0,0));
641 {
642 EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
643 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
644 m_Grid->GetPoint(id_node, grid_pnts[id_node].data());
645 if (m_BoundaryLayerNode[id_node]) {
646 int n_faces = m_Part.n2cGSize(id_node);
647 for (int i = 0; i < n_faces; ++i) {
648 vtkIdType id_cell = m_Part.n2cGG(id_node, i);
649 if (!m_BoundaryLayerCodes.contains(cell_code->GetValue(id_cell))) {
650 on_boundary[id_node] = true;
651 break;
652 }
653 }
654 if ( (m_NodeTypes[id_node] == EdgeNode || m_NodeTypes[id_node] == CornerNode)
655 && m_Part.isConvexNode(id_node, m_BoundaryLayerCodes) )
656 {
657 is_convex[id_node] = true;
658 }
659 }
660 }
661 }
662
663 limitHeights(1.0);
664
665 //laplacianIntersectSmoother(on_boundary);
666 //angleSmoother(on_boundary, is_convex, grid_pnts);
667 smoothUsingBLVectors();
668
669 //limitHeights(1.0);
670
671 cout << "heights computed" << endl;
672 }
673
674 void BoundaryLayerOperation::createSmoothShell(vtkUnstructuredGrid *shell_grid, int num_iter)
675 {
676 // new version based on VTK technology
677
678 // Set grid to normal*height
679 QVector<vec3_t> x_org(m_Grid->GetNumberOfPoints(), vec3_t(0,0,0));
680 QVector<vec3_t> x_new(m_Grid->GetNumberOfPoints(), vec3_t(0,0,0));
681 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
682 vec3_t x(0,0,0);
683 m_Grid->GetPoint(id_node, x_org[id_node].data());
684 x_new[id_node] = x_org[id_node] + m_Height[id_node]*m_BoundaryLayerVectors[id_node];
685 m_Grid->GetPoints()->SetPoint(id_node, x_new[id_node].data());
686 }
687
688 // extract wall boundaries to a separate grid
689 EG_VTKSP(vtkEgBoundaryCodesFilter, extract_wall);
690 EG_VTKSP(vtkUnstructuredGrid, wall_grid);
691 extract_wall->SetInputData(m_Grid);
692 extract_wall->SetBoundaryCodes(m_BoundaryLayerCodes);
693 extract_wall->Update();
694 //writeGrid(extract_wall->GetOutput(), "walls_from_vtk");
695
696 EG_VTKSP(vtkDataSetSurfaceFilter, grid_to_pdata);
697 grid_to_pdata->SetInputConnection(extract_wall->GetOutputPort());
698 EG_VTKSP(vtkLinearSubdivisionFilter, subdiv);
699 //EG_VTKSP(vtkButterflySubdivisionFilter, subdiv);
700 subdiv->SetInputConnection(grid_to_pdata->GetOutputPort());
701 subdiv->SetNumberOfSubdivisions(1);
702
703 //EG_VTKSP(vtkSmoothPolyDataFilter, smooth);
704 EG_VTKSP(vtkWindowedSincPolyDataFilter, smooth);
705 smooth->SetInputConnection(subdiv->GetOutputPort());
706 smooth->BoundarySmoothingOff();
707 smooth->FeatureEdgeSmoothingOn();
708 smooth->SetFeatureAngle(180);
709 smooth->SetEdgeAngle(180);
710 //smooth->SetNumberOfIterations(num_iter/100);
711 smooth->SetNumberOfIterations(100);
712 smooth->NormalizeCoordinatesOn();
713 double pb = m_ShellPassBand;
714 cout << "pass-band = " << pb << endl;
715 smooth->SetPassBand(pb);
716 //smooth->SetRelaxationFactor(0.05);
717 smooth->Update();
718 //cout << "rf = " << smooth->GetRelaxationFactor() << endl;
719
720 EG_VTKSP(vtkEgPolyDataToUnstructuredGridFilter, pdata_to_grid);
721 //pdata_to_grid->SetInputConnection(subdiv->GetOutputPort());
722 pdata_to_grid->SetInputConnection(smooth->GetOutputPort());
723 pdata_to_grid->Update();
724 makeCopy(pdata_to_grid->GetOutput(), shell_grid);
725
726 /*
727 GeometrySmoother smooth;
728 smooth.setGrid(shell_grid);
729 smooth.setNumberOfIterations(num_iter);
730 smooth.setConvexRelaxation(0.0);
731 smooth();
732 */
733
734 // reset grid to original points
735 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
736 if (m_BoundaryLayerNode[id_node]) {
737 m_Grid->GetPoints()->SetPoint(id_node, x_org[id_node].data());
738 }
739 }
740
741
742
743 // Set grid to normal*height
744 // And set weight factor of edges and corners to 1.
745 /*
746 QVector<vec3_t> x_org(m_Grid->GetNumberOfPoints(), vec3_t(0,0,0));
747 QVector<vec3_t> x_new(m_Grid->GetNumberOfPoints(), vec3_t(0,0,0));
748 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
749 vec3_t x(0,0,0);
750 m_Grid->GetPoint(id_node, x_org[id_node].data());
751 x_new[id_node] = x_org[id_node] + m_Height[id_node]*m_BoundaryLayerVectors[id_node];
752 m_Grid->GetPoints()->SetPoint(id_node, x_new[id_node].data());
753 }
754
755 // Laplacian on points
756 for (int iter = 0; iter < num_iter; iter++) {
757 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
758 if (m_BoundaryLayerNode[id_node]) {
759
760 // !!!!!!! Warning, this might not be robust enough for complex geometries !!!!!!!
761 bool shared_boundary = false;
762 {
763 QVector<int> bc_ids;
764 for (int i = 0; i < m_Part.n2bcGSize(id_node); ++i) {
765 int bc = m_Part.n2bcG(id_node, i);
766 if (!m_BoundaryLayerCodes.contains(bc)
767 && !bc_ids.contains(bc) )
768 {
769 bc_ids.append(bc);
770 }
771 }
772
773 if ( bc_ids.size() > 1 ) continue;
774
775 if (bc_ids.size() == 1) {
776 shared_boundary = true;
777 }
778 }
779
780 vec3_t avg_pnt(0,0,0);
781 int count = 0;
782 if (!shared_boundary) {
783 double area = 0;
784 for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
785 vtkIdType id_cell = m_Part.n2cGG(id_node, i);
786 vec3_t cell_ctr = cellCentre(m_Grid, id_cell);
787 double cell_area = cellVA(m_Grid, id_cell);
788 avg_pnt += cell_area*cell_ctr;
789 area += cell_area;
790 ++count;
791 }
792 if (count == 0)
793 continue;
794 avg_pnt *= 1.0/area;
795 }
796 if (shared_boundary) {
797 for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
798 vtkIdType id_neigh = m_Part.n2nGG(id_node, i);
799 if (on_boundary[id_neigh]) {
800 vec3_t x(0,0,0);
801 m_Grid->GetPoint(id_neigh, x.data());
802 avg_pnt += x;
803 ++count;
804 }
805 }
806 if (count == 0)
807 continue;
808 avg_pnt *= 1.0/count;
809 }
810
811 if (on_boundary[id_node]) {
812 vec3_t node_pnt(0,0,0);
813 m_Grid->GetPoint(id_node, node_pnt.data());
814 double n = 1./3.;
815 avg_pnt = (1-n)*node_pnt + n*(avg_pnt);
816 }
817 //vec3_t new_vec = avg_pnt - org_grid_pts[id_node];
818 //m_Height[id_node] = new_vec.abs();
819 //new_vec.normalise();
820 //m_BoundaryLayerVectors[id_node] = new_vec;
821 x_new[id_node] = avg_pnt;
822 }
823 }
824 // Set grid to new points
825 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
826 if (m_BoundaryLayerNode[id_node]) {
827 m_Grid->GetPoints()->SetPoint(id_node, x_new[id_node].data());
828 }
829 }
830 }
831
832 //- Create subgrid from shell
833 EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
834 {
835 MeshPartition shell_part(m_Grid);
836 QList<vtkIdType> shell_faces;
837 for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
838 if (isSurface(id_cell, m_Grid)) {
839 if (m_BoundaryLayerCodes.contains(cell_code->GetValue(id_cell))) {
840 shell_faces.append(id_cell);
841 }
842 }
843 }
844 shell_part.setCells(shell_faces);
845 shell_part.extractToVtkGrid(shell_grid);
846 }
847
848 // reset grid to original points
849 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
850 if (m_BoundaryLayerNode[id_node]) {
851 m_Grid->GetPoints()->SetPoint(id_node, x_org[id_node].data());
852 }
853 }
854 */
855 }
856
857 double BoundaryLayerOperation::largestAngle(vtkIdType id_node1, vtkIdType id_node2)
858 {
859 double alpha = 0;
860 QSet<vtkIdType> faces;
861 if (id_node1 == id_node2) {
862 for (int i = 0; i < m_Part.n2cGSize(id_node1); ++i) {
863 faces.insert(m_Part.n2cGG(id_node1, i));
864 }
865 } else {
866 getEdgeCells(id_node1, id_node2, faces);
867 }
868 foreach (vtkIdType id_face, faces) {
869 vec3_t n = (-1)*cellNormal(m_Grid, id_face);
870 alpha = max(alpha, GeometryTools::angle(n, m_BoundaryLayerVectors[id_node1]));
871 }
872 return alpha;
873 }
874
875 void BoundaryLayerOperation::fixBoundaryLayerVectors(const QList<vtkIdType> &bad_cells, int num_smooth_iter)
876 {
877 QVector<vtkIdType> bad_nodes;
878 getNodesFromCells(bad_cells, bad_nodes, m_Grid);
879 int num_changes;
880 do {
881 num_changes = 0;
882 foreach (vtkIdType id_node, bad_nodes) {
883 foreach (vtkIdType id_snap, m_SnapPoints[id_node]) {
884 double node_alpha = largestAngle(id_node, id_snap);
885 double snap_alpha = largestAngle(id_snap, id_node);
886 if (snap_alpha > node_alpha) {
887 vec3_t dv = m_BoundaryLayerVectors[id_snap] - m_BoundaryLayerVectors[id_node];
888 if (dv.abs() > 1e-4) {
889 m_BoundaryLayerVectors[id_node] = m_BoundaryLayerVectors[id_snap];
890 ++num_changes;
891 }
892 }
893 }
894 }
895 } while (num_changes);
896 QVector<bool> node_fixed(m_Grid->GetNumberOfPoints(), false);
897 foreach (vtkIdType id_node, bad_nodes) {
898 node_fixed[id_node] = true;
899 }
900 correctBoundaryLayerVectors();
901 smoothBoundaryLayerVectors(num_smooth_iter, 1.0, 0.0, &node_fixed);
902 }
903
904 void BoundaryLayerOperation::writeWallGrid(QString file_name, int counter)
905 {
906 if (counter >= 0) {
907 QString counter_txt;
908 counter_txt.setNum(counter);
909 counter_txt = counter_txt.rightJustified(3, '0');
910 file_name += "_" + counter_txt;
911 }
912 MeshPartition wall_part(m_Grid);
913 wall_part.setBCs(m_BoundaryLayerCodes);
914 EG_VTKSP(vtkUnstructuredGrid, wall_grid);
915 wall_part.extractToVtkGrid(wall_grid);
916 writeGrid(wall_grid, file_name);
917 }
918
919 void BoundaryLayerOperation::smoothUsingBLVectors()
920 {
921 // create shell
922 EG_VTKSP(vtkUnstructuredGrid, shell_grid);
923 createSmoothShell(shell_grid, m_ShellPassBand);
924
925 newHeightFromShellIntersect(shell_grid, 1.0);
926 //writeGrid(shell_grid, "shell");
927
928 EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
929 //writeWallGrid("walls", 0);
930
931 QList<vtkIdType> bad_cells;
932
933 double w_iso = 1.0;
934 double w_dir = 0.0;
935 double dw = 0.01;
936
937 int iter = 0;
938
939 int last_num_bad = m_Grid->GetNumberOfCells();
940
941 //snapAllVectorsToShell(shell_grid);
942 //return;
943 smoothBoundaryLayerVectors(10, w_iso, w_dir);
944
945 do {
946 bad_cells.clear();
947 for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
948 if (m_BoundaryLayerCodes.contains(cell_code->GetValue(id_cell))) {
949 if (!faceFine(id_cell, 1.0)) {
950 bad_cells.append(id_cell);
951 }
952 }
953 }
954 cout << "found " << bad_cells.size() << " distorted faces" << endl;
955
956 if (bad_cells.size() > 0) {
957 if (bad_cells.size() <= last_num_bad) {
958 last_num_bad = bad_cells.size();
959 smoothBoundaryLayerVectors(10, w_iso, w_dir);
960 newHeightFromShellIntersect(shell_grid, 1.0);
961 } else {
962 cout << "cannot fix completely -- terminating the loop!" << endl;
963 //cout << "moving to global under relaxation now ..." << endl;
964 break;
965 }
966 }
967 w_iso -= dw;
968 w_dir += dw;
969 ++iter;
970 } while (bad_cells.size() > 0 && w_iso >= 0.0);
971
972 /*
973 double relax = 0.95;
974
975 do {
976
977 cout << "relaxation factor: " << relax << endl;
978 newHeightFromShellIntersect(shell_grid, relax);
979
980 bad_cells.clear();
981 for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
982 if (m_BoundaryLayerCodes.contains(cell_code->GetValue(id_cell))) {
983 if (!faceFine(id_cell, 1.0)) {
984 bad_cells.append(id_cell);
985 }
986 }
987 }
988 cout << "found " << bad_cells.size() << " distorted faces" << endl;
989 relax -= 0.05;
990 } while (bad_cells.size() > 0 && relax >= 0.25);
991 */
992
993 //swapEdgesToMatchShell(shell_grid, deg2rad(5.0));
994 }
995
996 bool BoundaryLayerOperation::checkVectorForNode(vec3_t v, vtkIdType id_node)
997 {
998 for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
999 vtkIdType id_face = m_Part.n2cGG(id_node, i);
1000 EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
1001 if (m_BoundaryLayerCodes.contains(cell_code->GetValue(id_face))) {
1002 vec3_t n = cellNormal(m_Grid, id_face);
1003 if (n*v > 0) {
1004 return false;
1005 }
1006 }
1007 }
1008 return true;
1009 }
1010
1011 vec3_t BoundaryLayerOperation::snapToShell(CadInterface* cad, vtkIdType id_node)
1012 {
1013 bool dbg = false;
1014
1015 if (id_node == 0) {
1016 cout << "break" << endl;
1017 dbg = true;
1018 }
1019
1020 vec3_t x_node;
1021 m_Grid->GetPoint(id_node, x_node.data());
1022 vec3_t x_snap = cad->snap(x_node);
1023 if (dbg) cout << x_snap << endl;
1024 if (checkVectorForNode(x_snap - x_node, id_node)) {
1025 return x_snap;
1026 }
1027
1028 // initial guess
1029 {
1030 x_snap = x_node + m_Height[id_node]*m_BoundaryLayerVectors[id_node];
1031 QVector<QPair<vec3_t, vtkIdType> > intersections;
1032 cad->computeIntersections(x_node, m_BoundaryLayerVectors[id_node], intersections);
1033 double h_min = EG_LARGE_REAL;
1034 bool found = false;
1035 vec3_t shell_vector = m_BoundaryLayerVectors[id_node];
1036 for (int i = 0; i < intersections.size(); ++i) {
1037 vec3_t xi = intersections[i].first;
1038 if (dbg) cout << "xi=" << xi << endl;
1039 vec3_t vi = xi - x_node;
1040 if (vi*m_BoundaryLayerVectors[id_node] > 0) {
1041 double h = vi.abs();
1042 if (h < h_min) {
1043 if (dbg) cout << "h=" << h << endl;
1044 found = true;
1045 shell_vector = vi;
1046 h_min = h;
1047 }
1048 }
1049 }
1050 if (found) {
1051 x_snap = x_node + shell_vector;
1052 }
1053 }
1054 if (dbg) cout << x_snap << endl;
1055
1056
1057 double dist_old = 0;
1058 double dist_new = m_Height[id_node];
1059 while (fabs(dist_new - dist_old) > 1e-3*m_Height[id_node]) {
1060 double w = 1.0;
1061 vec3_t x_new = x_snap;
1062 do {
1063 w *= 0.5;
1064 if (w < 1e-3) {
1065 QString msg;
1066 msg.setNum(id_node);
1067 msg = "unable to snap node " + msg + " to shell";
1068 EG_ERR_RETURN(msg);
1069 }
1070 x_new = cad->snap(w*x_node + (1-w)*x_snap);
1071 if (dbg) cout << "w=" << w << ", x_new=" << x_new << endl;
1072 } while (!checkVectorForNode(x_new - x_node, id_node));
1073 dist_old = dist_new;
1074 dist_new = (x_new - x_snap).abs();
1075 x_snap = x_new;
1076 }
1077
1078 return x_snap;
1079 }
1080
1081 void BoundaryLayerOperation::snapAllVectorsToShell(vtkUnstructuredGrid *shell_grid)
1082 {
1083 writeBoundaryLayerVectors("normals");
1084 CgalTriCadInterface cad(shell_grid);
1085 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
1086 if (m_BoundaryLayerNode[id_node]) {
1087 //if (m_NodeTypes[id_node] == CornerNode || m_NodeTypes[id_node] == EdgeNode) {
1088 vec3_t x1;
1089 m_Grid->GetPoint(id_node, x1.data());
1090 vec3_t x2 = snapToShell(&cad, id_node);
1091 m_BoundaryLayerVectors[id_node] = x2 - x1;
1092 m_Height[id_node] = m_BoundaryLayerVectors[id_node].abs();
1093 m_BoundaryLayerVectors[id_node].normalise();
1094 }
1095 }
1096 correctBoundaryLayerVectors();
1097 }
1098
1099 // Compute intersection points with a shell following m_BoundaryLayerVector
1100 // Updates m_Height
1101 void BoundaryLayerOperation::newHeightFromShellIntersect(vtkUnstructuredGrid* shell_grid, double relax)
1102 {
1103 CgalTriCadInterface cad(shell_grid);
1104 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
1105 if (m_BoundaryLayerNode[id_node]) {
1106 vec3_t x;
1107 m_Grid->GetPoint(id_node, x.data());
1108 QVector<QPair<vec3_t, vtkIdType> > intersections;
1109 cad.computeIntersections(x, m_BoundaryLayerVectors[id_node], intersections);
1110 double h = 2*m_Height[id_node];
1111 bool found = false;
1112 vec3_t layer_vector = m_BoundaryLayerVectors[id_node];
1113 for (int i = 0; i < intersections.size(); ++i) {
1114 vec3_t xi = intersections[i].first;
1115 vec3_t vi = xi - x;
1116 if (vi*m_BoundaryLayerVectors[id_node] > 0) {
1117 double hi = vi.abs();
1118 if (hi < h) {
1119 h = hi;
1120 layer_vector = vi.normalise();
1121 found = true;
1122 }
1123 }
1124 }
1125 if (found) {
1126 m_Height[id_node] = relax*h;
1127 m_BoundaryLayerVectors[id_node] = layer_vector;
1128 }
1129 }
1130 }
1131 }
1132
1133 void BoundaryLayerOperation::limitSizeAndAngleErrors()
1134 {
1135 bool done;
1136 int iter = 0;
1137 do {
1138 done = true;
1139 QVector<double> new_scale(m_Grid->GetNumberOfPoints(), 1.0);
1140 QVector<int> count(m_Grid->GetNumberOfPoints(), 1);
1141 for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
1142 if (isSurface(id_cell, m_Grid)) {
1143 vtkIdType num_pts, *pts;
1144 m_Grid->GetCellPoints(id_cell, num_pts, pts);
1145 bool check_face = true;
1146 for (vtkIdType i = 0; i < num_pts; ++i) {
1147 if (!m_BoundaryLayerNode[pts[i]]) {
1148 check_face = false;
1149 break;
1150 }
1151 }
1152 if (check_face) {
1153 if (!faceFine(id_cell, 1)) {
1154 done = false;
1155 double scale1 = 0.1;
1156 double scale2 = 1;
1157
1158 bool found;
1159 do {
1160 found = false;
1161 int num_steps = 10;
1162 for (int i = 1; i <= num_steps; ++i) {
1163 double s = scale2 - i*(scale2 - scale1)/num_steps;
1164 if (faceFine(id_cell, s)) {
1165 found = true;
1166 scale1 = s;
1167 scale2 -= (i-1)*(scale2 - scale1)/num_steps;
1168 break;
1169 }
1170 }
1171 if (!found) {
1172 scale1 = scale2;
1173 }
1174
1175 } while ((scale2 - scale1) > 1e-4 && found);
1176
1177 double scale = 0.5*(scale1 + scale2);
1178 for (vtkIdType i = 0; i < num_pts; ++i) {
1179 new_scale[pts[i]] += scale;
1180 ++count[pts[i]];
1181 }
1182 }
1183 }
1184 }
1185 }
1186 double relax = 0.2;
1187 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
1188 double h_new = m_Height[id_node]*new_scale[id_node]/count[id_node];
1189 m_Height[id_node] = relax*h_new + (1 - relax)*m_Height[id_node];
1190 }
1191 //done = true;
1192 ++iter;
1193 if (iter >= 10) {
1194 done = true;
1195 }
1196 } while (!done);
1197 }
1198
1199
1200 void BoundaryLayerOperation::angleSmoother(const QVector<bool>& on_boundary, const QVector<bool>& is_convex, QVector<vec3_t>& grid_pnts)
1201 {
1202 int n_iter = 20;
1203 double weight_const = 1.;
1204 const double PI = 3.14159265359;
1205
1206 EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
1207 for (int iter = 0; iter < n_iter; ++iter) {
1208 // Set points to bl_normal*height
1209 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
1210 if (m_BoundaryLayerNode[id_node]) {
1211 vec3_t x(0,0,0);
1212 m_Grid->GetPoint(id_node, x.data());
1213 x += m_Height[id_node]*m_BoundaryLayerVectors[id_node];
1214 m_Grid->GetPoints()->SetPoint(id_node, x.data());
1215 }
1216 }
1217
1218 QVector<int> move_count(grid_pnts.size());
1219 QVector<vec3_t> grid_smoothed(grid_pnts.size(), vec3_t(0,0,0));
1220 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
1221 if (m_BoundaryLayerNode[id_node]) {
1222 for (vtkIdType i = 0; i < m_Part.n2nGSize(id_node); ++i) {
1223 vtkIdType id_neigh = m_Part.n2nGG(id_node, i);
1224 if (!m_BoundaryLayerNode[id_neigh]) continue;
1225 if (id_neigh < id_node) continue;
1226 if (on_boundary[id_node] && on_boundary[id_neigh]) continue;
1227
1228 QList<vtkIdType> edge_faces;
1229 m_Part.getEdgeFaces(id_node, id_neigh, edge_faces);
1230 // Test for 2 faces
1231 int count = 0;
1232 for (int j = 0; j < edge_faces.size(); ++j) {
1233 if (m_BoundaryLayerCodes.contains(cell_code->GetValue(edge_faces[j]))) {
1234 ++count;
1235 }
1236 }
1237 if (count < 2 ) continue;
1238 if (edge_faces.size() > 2) EG_BUG;
1239
1240 // Prepare cell info
1241 vtkIdType id_cell_1 = edge_faces[0];
1242 vtkIdType id_cell_2 = edge_faces[1];
1243 vtkIdType npts_c1, *pts_c1;
1244 vtkIdType npts_c2, *pts_c2;
1245 m_Grid->GetCellPoints(id_cell_1, npts_c1, pts_c1);
1246 m_Grid->GetCellPoints(id_cell_2, npts_c2, pts_c2);
1247 if (npts_c1 != 3 || npts_c2 !=3) EG_BUG;
1248
1249 vtkIdType id_n3_c1;
1250 vtkIdType id_n3_c2;
1251 for (int j = 0; j < npts_c1; j++) {
1252 if ( pts_c1[j] != id_node && pts_c1[j] != id_neigh) {
1253 id_n3_c1 = pts_c1[j];
1254 }
1255 }
1256 for (int j = 0; j < npts_c2; j++) {
1257 if ( pts_c2[j] != id_node && pts_c2[j] != id_neigh) {
1258 id_n3_c2 = pts_c2[j];
1259 }
1260 }
1261
1262 vec3_t normal_c1 = cellNormal(m_Grid, id_cell_1);
1263 vec3_t normal_c2 = cellNormal(m_Grid, id_cell_2);
1264
1265 double angle = GeometryTools::angle(normal_c1, normal_c2);
1266 if (rad2deg(angle) < 1) continue;
1267
1268 double spring_angle = weight_const*angle*angle/(PI*PI);
1269
1270 vec3_t x_node(0,0,0);
1271 vec3_t x_neigh(0,0,0);
1272 vec3_t x_n3_c1(0,0,0);
1273 vec3_t x_n3_c2(0,0,0);
1274 m_Grid->GetPoint(id_node, x_node.data());
1275 m_Grid->GetPoint(id_neigh, x_neigh.data());
1276 m_Grid->GetPoint(id_n3_c1, x_n3_c1.data());
1277 m_Grid->GetPoint(id_n3_c2, x_n3_c2.data());
1278
1279 vec3_t axis = x_node.cross(x_neigh);
1280 vec3_t v1 = x_n3_c1 - x_node;
1281 vec3_t v2 = x_n3_c2 - x_node;
1282
1283 vec3_t cross_vector = axis.cross(v1);
1284 vec3_t v1_rot(0,0,0);
1285 vec3_t v2_rot(0,0,0);
1286 if (cross_vector*normal_c1 > 0) {
1287 v1_rot = GeometryTools::rotate(v1, axis, spring_angle);
1288 v2_rot = GeometryTools::rotate(v2, axis, -spring_angle);
1289 }
1290 else {
1291 v1_rot = GeometryTools::rotate(v1, axis, -spring_angle);
1292 v2_rot = GeometryTools::rotate(v2, axis, spring_angle);
1293 }
1294
1295 grid_smoothed[id_n3_c1] += x_node + v1_rot;
1296 move_count[id_n3_c1] += 1;
1297 grid_smoothed[id_n3_c2] += x_node + v2_rot;
1298 move_count[id_n3_c2] += 1;
1299 }
1300 }
1301 }
1302
1303 QVector<double> h_new = m_Height;
1304 QVector<vec3_t> new_BoundaryLayerVectors = m_BoundaryLayerVectors;
1305 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
1306 m_Grid->GetPoints()->SetPoint(id_node, grid_pnts[id_node].data());
1307 if (m_BoundaryLayerNode[id_node]) {
1308 if (move_count[id_node] > 0) {
1309 grid_smoothed[id_node] *= 1.0/move_count[id_node];
1310
1311 vec3_t new_norm = grid_smoothed[id_node] - grid_pnts[id_node];
1312 h_new[id_node] = new_norm.abs();
1313 new_norm.normalise();
1314 new_BoundaryLayerVectors[id_node] = new_norm;
1315 }
1316 }
1317 }
1318 double max_diff = 0;
1319 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
1320 double diff = m_Height[id_node] - h_new[id_node];
1321 diff = std::sqrt(diff*diff);
1322 max_diff = std::max(max_diff, diff);
1323 }
1324 cout << "========================== max diff->" << max_diff << endl;
1325 m_Height = h_new;
1326 m_BoundaryLayerVectors = new_BoundaryLayerVectors;
1327 }
1328 }
1329
1330 int BoundaryLayerOperation::limitHeights(double safety_factor)
1331 {
1332 EG_VTKDCC(vtkIntArray, bc, m_Grid, "cell_code");
1333
1334 QVector<bool> limited(m_Grid->GetNumberOfPoints(), false);
1335 // save original node positions to x_old
1336 QVector<vec3_t> x_old(m_Grid->GetNumberOfPoints());
1337 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
1338 m_Grid->GetPoint(id_node, x_old[id_node].data());
1339 }
1340
1341 for (int pass = 1; pass <= 2; ++pass) {
1342
1343 // move nodes for second pass
1344 if (pass == 2) {
1345 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
1346 if (m_BoundaryLayerNode[id_node]) {
1347 vec3_t x = x_old[id_node] + m_Height[id_node]*m_BoundaryLayerVectors[id_node];
1348 m_Grid->GetPoints()->SetPoint(id_node, x.data());
1349 }
1350 }
1351 }
1352
1353 CgalTriCadInterface cad(m_Grid);
1354
1355 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
1356 if (m_BoundaryLayerNode[id_node]) {
1357 QList<vtkIdType> cells_of_node;
1358 for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
1359 cells_of_node.append(m_Part.n2cGG(id_node, i));
1360 }
1361 QVector<QPair<vec3_t, vtkIdType> > intersections;
1362 cad.computeIntersections(x_old[id_node], m_BoundaryLayerVectors[id_node], intersections);
1363 for (int i = 0; i < intersections.size(); ++i) {
1364 QPair<vec3_t,vtkIdType> inters = intersections[i];
1365 vec3_t xi = inters.first;
1366 vtkIdType id_tri = inters.second;
1367 if (!cells_of_node.contains(id_tri)) {
1368
1369 double crit_angle = deg2rad(200.0); // consider all intersections
1370
1371 if (m_LayerAdjacentBoundaryCodes.contains(bc->GetValue(id_tri))) {
1372 crit_angle = deg2rad(85.0); // different angle for adjacent boundaries
1373 }
1374
1375 vec3_t dx = xi - x_old[id_node];
1376 double alpha = angle(dx, cellNormal(m_Grid, id_tri));
1377 if (dx*m_BoundaryLayerVectors[id_node] > 0 && alpha < crit_angle) {
1378 //double h_max = (1.0 - 0.5*safety_factor*(1.0 - m_MaxHeightInGaps))*dx.abs();
1379 double h_max = safety_factor*m_MaxHeightInGaps*dx.abs();
1380 if (h_max < m_Height[id_node]) {
1381 m_Height[id_node] = h_max;
1382 limited[id_node] = true;
1383 }
1384 }
1385 }
1386 }
1387 }
1388 }
1389 }
1390
1391 // reset node positions and count nodes where the height has been limited
1392 int num_limited = 0;
1393 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
1394 m_Grid->GetPoints()->SetPoint(id_node, x_old[id_node].data());
1395 if (limited[id_node]) {
1396 ++num_limited;
1397 }
1398 }
1399
1400 return num_limited;
1401 }
1402
1403 void BoundaryLayerOperation::laplacianIntersectSmoother(const QVector<bool>& on_boundary)
1404 {
1405 int n_loops = 1;
1406 int n_iter = 4;
1407
1408 // Set grid to normal*height
1409 // And set weight factor of edges and corners to 1.
1410 QVector<vec3_t> org_grid(m_Grid->GetNumberOfPoints(), vec3_t(0,0,0));
1411 QVector<vec3_t> new_grid(m_Grid->GetNumberOfPoints(), vec3_t(0,0,0));
1412 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
1413 vec3_t x(0,0,0);
1414 m_Grid->GetPoint(id_node, x.data());
1415 org_grid[id_node] = x;
1416 }
1417
1418 // Laplacian on points
1419 for (int loops = 0; loops < n_loops; ++loops) {
1420 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
1421 if (m_BoundaryLayerNode[id_node]) {
1422 vec3_t x(0,0,0);
1423 m_Grid->GetPoint(id_node, x.data());
1424 x += m_Height[id_node]*m_BoundaryLayerVectors[id_node];
1425 m_Grid->GetPoints()->SetPoint(id_node, x.data());
1426 new_grid[id_node] = x;
1427 }
1428 }
1429 for (int iter = 0; iter < n_iter; iter++) {
1430 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
1431 if (m_BoundaryLayerNode[id_node]) {
1432
1433 // !!!!!!! Warning, this might not be robust enough for complex geometries !!!!!!!
1434 bool shared_boundary = false;
1435 {
1436 QVector<int> bc_ids;
1437 for (int i = 0; i < m_Part.n2bcGSize(id_node); ++i) {
1438 int bc = m_Part.n2bcG(id_node, i);
1439 if (!m_BoundaryLayerCodes.contains(bc)
1440 && !bc_ids.contains(bc) )
1441 {
1442 bc_ids.append(bc);
1443 }
1444 }
1445
1446 if ( bc_ids.size() > 1 ) continue;
1447
1448 if (bc_ids.size() == 1) {
1449 shared_boundary = true;
1450 }
1451 }
1452
1453 vec3_t avg_pnt(0,0,0);
1454 int count = 0;
1455 if (!shared_boundary) {
1456 double area = 0;
1457 for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
1458 vtkIdType id_cell = m_Part.n2cGG(id_node, i);
1459 vec3_t cell_ctr = cellCentre(m_Grid, id_cell);
1460 double cell_area = cellVA(m_Grid, id_cell);
1461 avg_pnt += cell_area*cell_ctr;
1462 area += cell_area;
1463 ++count;
1464 }
1465 if (count == 0)
1466 continue;
1467 avg_pnt *= 1.0/area;
1468 }
1469 if (shared_boundary) {
1470 for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
1471 vtkIdType id_neigh = m_Part.n2nGG(id_node, i);
1472 if (on_boundary[id_neigh]) {
1473 vec3_t x(0,0,0);
1474 m_Grid->GetPoint(id_neigh, x.data());
1475 avg_pnt += x;
1476 ++count;
1477 }
1478 }
1479 if (count == 0)
1480 continue;
1481 avg_pnt *= 1.0/count;
1482 }
1483
1484 if (on_boundary[id_node]) {
1485 vec3_t node_pnt(0,0,0);
1486 m_Grid->GetPoint(id_node, node_pnt.data());
1487 double n = 1./3.;
1488 avg_pnt = (1-n)*node_pnt + n*(avg_pnt);
1489 }
1490 //vec3_t new_vec = avg_pnt - org_grid[id_node];
1491 //m_Height[id_node] = new_vec.abs();
1492 //new_vec.normalise();
1493 //m_BoundaryLayerVectors[id_node] = new_vec;
1494 new_grid[id_node] = avg_pnt;
1495 }
1496 }
1497 // Set grid to new points
1498 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
1499 if (m_BoundaryLayerNode[id_node]) {
1500 m_Grid->GetPoints()->SetPoint(id_node, new_grid[id_node].data());
1501 }
1502 }
1503 }
1504 return;
1505
1506 EG_VTKSP(vtkUnstructuredGrid, bl_grid);
1507 EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
1508 MeshPartition bl_part(m_Grid);
1509 QList<vtkIdType> bl_faces;
1510 for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
1511 if (isSurface(id_cell, m_Grid)) {
1512 if (m_BoundaryLayerCodes.contains(cell_code->GetValue(id_cell))) {
1513 bl_faces.append(id_cell);
1514 }
1515 }
1516 }
1517 bl_part.setCells(bl_faces);
1518 bl_part.extractToVtkGrid(bl_grid);
1519 CgalTriCadInterface cad(bl_grid);
1520 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
1521 if (m_BoundaryLayerNode[id_node]) {
1522 vec3_t org_pnt = org_grid[id_node];
1523 QVector<QPair<vec3_t, vtkIdType> > intersections;
1524 cad.computeIntersections(org_pnt, m_BoundaryLayerVectors[id_node], intersections);
1525
1526 QVector<QPair<double, vec3_t> > intersect_vec;
1527 for (int i = 0; i < intersections.size(); ++i) {
1528 vec3_t xi = intersections[i].first;
1529 vec3_t new_vec = xi - org_pnt;
1530 if (new_vec*m_BoundaryLayerVectors[id_node] < 0) {
1531 continue;
1532 }
1533 double length = new_vec.abs();
1534 intersect_vec.append(QPair<double, vec3_t>(length, xi));
1535 }
1536 vec3_t new_pnt(0,0,0);
1537 double length = 1e99;
1538 for (int i = 0; i < intersect_vec.size(); ++i) {
1539 if (intersect_vec[i].first < length) {
1540 length = intersect_vec[i].first;
1541 new_pnt = intersect_vec[i].second;
1542 }
1543 }
1544 //vec3_t new_vec = new_pnt - org_pnt;
1545 //m_Height[id_node] = new_vec.abs();
1546 //new_vec.normalise();
1547 //m_BoundaryLayerVectors[id_node] = new_vec;
1548 m_Height[id_node] = length;
1549 }
1550 }
1551
1552 // Set grid to original points before return
1553 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
1554 if (m_BoundaryLayerNode[id_node]) {
1555 vec3_t x = org_grid[id_node];
1556 m_Grid->GetPoints()->SetPoint(id_node, x.data());
1557 }
1558 }
1559
1560
1561 QVector<double> projected_height = m_Height;
1562 limitSizeAndAngleErrors();
1563 return;
1564
1565 QVector<double> height_diff(m_Height.size(), 0.0);
1566 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
1567 if (m_BoundaryLayerNode[id_node]) {
1568 height_diff[id_node] = projected_height[id_node] - m_Height[id_node];
1569 }
1570 }
1571
1572 int n_iter_smooth = 0;
1573 if (loops == 1) {
1574 n_iter_smooth = 1000;
1575 }
1576 else {
1577 n_iter_smooth = 1000;
1578 }
1579 for (int iter = 0; iter < n_iter_smooth; iter++) {
1580 QVector<double> new_height = height_diff;
1581 //QVector<double> new_height(m_Height.size(), 0.0);
1582 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
1583 if (m_BoundaryLayerNode[id_node]) {
1584 int count = 1;
1585 double org_height = height_diff[id_node];
1586 double avg_height = 0;
1587 for(int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
1588 vtkIdType id_neigh = m_Part.n2nGG(id_node, i);
1589 avg_height += height_diff[id_neigh];
1590 ++count;
1591 }
1592 avg_height /= count;
1593 if (org_height < avg_height) {
1594 new_height[id_node] = avg_height;
1595 }
1596 }
1597 }
1598 height_diff = new_height;
1599 }
1600 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
1601 if (m_BoundaryLayerNode[id_node]) {
1602 m_Height[id_node] = projected_height[id_node] - height_diff[id_node];
1603 }
1604 }
1605 // End of global pass loop:
1606 }
1607
1608 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
1609 if (m_BoundaryLayerNode[id_node]) {
1610 vec3_t x = org_grid[id_node];
1611 m_Grid->GetPoints()->SetPoint(id_node, x.data());
1612 }
1613 }
1614 }
1615
1616 bool BoundaryLayerOperation::swapRequired(stencil_t stencil, CadInterface *cad, double threshold_angle)
1617 {
1618
1619 QVector<vec3_t> x(4);
1620 m_Grid->GetPoint(stencil.id_node[0], x[0].data());
1621 m_Grid->GetPoint(stencil.p1, x[1].data());
1622 m_Grid->GetPoint(stencil.id_node[1], x[2].data());
1623 m_Grid->GetPoint(stencil.p2, x[3].data());
1624
1625 // centres of triangles
1626 vec3_t xc_tri_11 = 0.333333*(x[0] + x[1] + x[3]);
1627 vec3_t xc_tri_12 = 0.333333*(x[1] + x[2] + x[3]);
1628 vec3_t xc_tri_21 = 0.333333*(x[0] + x[1] + x[2]);
1629 vec3_t xc_tri_22 = 0.333333*(x[2] + x[3] + x[0]);
1630
1631 cad->snap(xc_tri_11);
1632 vec3_t n_snap_11 = cad->getLastNormal();
1633 cad->snap(xc_tri_12);
1634 vec3_t n_snap_12 = cad->getLastNormal();
1635 cad->snap(xc_tri_21);
1636 vec3_t n_snap_21 = cad->getLastNormal();
1637 cad->snap(xc_tri_22);
1638 vec3_t n_snap_22 = cad->getLastNormal();
1639
1640 vec3_t n_tri_11 = GeometryTools::triNormal(x[0], x[1], x[3]).normalise();
1641 vec3_t n_tri_12 = GeometryTools::triNormal(x[1], x[2], x[3]).normalise();
1642 vec3_t n_tri_21 = GeometryTools::triNormal(x[0], x[1], x[2]).normalise();
1643 vec3_t n_tri_22 = GeometryTools::triNormal(x[2], x[3], x[0]).normalise();
1644
1645 double alpha_11 = GeometryTools::angle(n_tri_11, n_snap_11);
1646 double alpha_12 = GeometryTools::angle(n_tri_12, n_snap_12);
1647 double alpha_21 = GeometryTools::angle(n_tri_21, n_snap_21);
1648 double alpha_22 = GeometryTools::angle(n_tri_22, n_snap_22);
1649
1650 if (alpha_11 > threshold_angle || alpha_12 > threshold_angle) {
1651 if (alpha_11 > alpha_21 && alpha_12 > alpha_22) {
1652 return true;
1653 }
1654 }
1655
1656 return false;
1657 }
1658
1659 void BoundaryLayerOperation::swapEdgesToMatchShell(vtkUnstructuredGrid *shell_grid, double threshold_angle)
1660 {
1661 // Set grid to normal*height
1662 QVector<vec3_t> x_org(m_Grid->GetNumberOfPoints());
1663 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
1664 if (m_BoundaryLayerNode[id_node]) {
1665 m_Grid->GetPoint(id_node, x_org[id_node].data());
1666 vec3_t x = x_org[id_node] + m_Height[id_node]*m_BoundaryLayerVectors[id_node];
1667 m_Grid->GetPoints()->SetPoint(id_node, x.data());
1668 }
1669 }
1670
1671 CgalTriCadInterface cad(shell_grid);
1672 int num_swaps = 0;
1673 int count = 0;
1674
1675 do {
1676 num_swaps = 0;
1677 m_Part.setAllCells();
1678 EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
1679 EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type");
1680 QVector<bool> swapped(m_Grid->GetNumberOfCells(), false);
1681 QVector<bool> marked(m_Grid->GetNumberOfCells(), false);
1682 for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
1683 if (m_BoundaryLayerCodes.contains(cell_code->GetValue(id_cell)) && m_Grid->GetCellType(id_cell) == VTK_TRIANGLE) { //if it is a selected triangle
1684 if (!marked[id_cell] && !swapped[id_cell]) {
1685 for (int j = 0; j < 3; ++j) {
1686 stencil_t stencil = getStencil(id_cell, j);
1687 if (stencil.id_cell.size() == 2 && stencil.sameBC) {
1688 if (swapRequired(stencil, &cad, threshold_angle)) {
1689 marked[stencil.id_cell[0]] = true;
1690 marked[stencil.id_cell[1]] = true;
1691 for (int k = 0; k < m_Part.n2cGSize(stencil.id_node[0]); ++k) {
1692 vtkIdType id_neigh = m_Part.n2cGG(stencil.id_node[0], k);
1693 marked[id_neigh] = true;
1694 }
1695 for (int k = 0; k < m_Part.n2cGSize(stencil.id_node[1]); ++k) {
1696 vtkIdType id_neigh = m_Part.n2cGG(stencil.id_node[1], k);
1697 marked[id_neigh] = true;
1698 }
1699 for (int k = 0; k < m_Part.n2cGSize(stencil.p1); ++k) {
1700 vtkIdType id_neigh = m_Part.n2cGG(stencil.p1, k);
1701 marked[id_neigh] = true;
1702 }
1703 for (int k = 0; k < m_Part.n2cGSize(stencil.p2); ++k) {
1704 vtkIdType id_neigh = m_Part.n2cGG(stencil.p2, k);
1705 marked[id_neigh] = true;
1706 }
1707
1708 vtkIdType new_pts1[3], new_pts2[3];
1709 new_pts1[0] = stencil.p1;
1710 new_pts1[1] = stencil.id_node[1];
1711 new_pts1[2] = stencil.id_node[0];
1712 new_pts2[0] = stencil.id_node[1];
1713 new_pts2[1] = stencil.p2;
1714 new_pts2[2] = stencil.id_node[0];
1715 /*
1716 vtkIdType old_pts1[3], old_pts2[3];
1717 old_pts1[0] = stencil.id_node[0];
1718 old_pts1[1] = stencil.p1;
1719 old_pts1[2] = stencil.p2;
1720 old_pts2[0] = stencil.id_node[1];
1721 old_pts2[1] = stencil.p2;
1722 old_pts2[2] = stencil.p1;
1723 */
1724 m_Grid->ReplaceCell(stencil.id_cell[0], 3, new_pts1);
1725 m_Grid->ReplaceCell(stencil.id_cell[1], 3, new_pts2);
1726
1727 swapped[stencil.id_cell[0]] = true;
1728 swapped[stencil.id_cell[1]] = true;
1729 ++num_swaps;
1730 break;
1731 }
1732 }
1733 }
1734 }
1735 }
1736 }
1737 ++count;
1738 cout << count << ": " << num_swaps << endl;
1739 } while (num_swaps > 0 && count < 3);
1740
1741
1742 // reset grid to original points
1743 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
1744 if (m_BoundaryLayerNode[id_node]) {
1745 m_Grid->GetPoints()->SetPoint(id_node, x_org[id_node].data());
1746 }
1747 }
1748 }
1749
1750
This is a mirror repository of the new location for engrid@github