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