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