slider for smoothing
[engrid-github.git] / src / libengrid / createboundarylayershell.cpp
blob26ba7b9d898715eba458cfe02aa014bb84d62167
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 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
22 #include "createboundarylayershell.h"
24 #include "createvolumemesh.h"
25 #include "swaptriangles.h"
26 #include "deletetetras.h"
27 #include "deletecells.h"
28 #include "meshpartition.h"
29 #include "deletevolumegrid.h"
30 #include "laplacesmoother.h"
31 #include "surfacemeshsmoother.h"
32 #include "deletecells.h"
33 #include "stitchholes.h"
35 CreateBoundaryLayerShell::CreateBoundaryLayerShell()
37 m_RestGrid = vtkSmartPointer<vtkUnstructuredGrid>::New();
38 m_OriginalGrid = vtkSmartPointer<vtkUnstructuredGrid>::New();
39 m_PrismaticGrid = vtkSmartPointer<vtkUnstructuredGrid>::New();
40 m_Success = false;
43 void CreateBoundaryLayerShell::prepare()
45 m_Part.trackGrid(m_Grid);
47 DeleteVolumeGrid delete_volume;
48 delete_volume.setGrid(m_Grid);
49 delete_volume.setAllCells();
50 delete_volume();
52 readSettings();
53 setAllCells();
54 getSurfaceCells(m_BoundaryLayerCodes, layer_cells, m_Grid);
56 // fill m_LayerAdjacentBoundaryCodes
57 EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
58 foreach (vtkIdType id_cell, layer_cells) {
59 for (int i = 0; i < m_Part.c2cGSize(id_cell); ++i) {
60 vtkIdType id_neigh = m_Part.c2cGG(id_cell, i);
61 int bc = cell_code->GetValue(id_neigh);
62 if (!m_BoundaryLayerCodes.contains(bc)) {
63 m_LayerAdjacentBoundaryCodes.insert(bc);
68 // compute normals and origins of adjacent planes
69 m_LayerAdjacentNormals.clear();
70 m_LayerAdjacentOrigins.clear();
71 foreach (int bc, m_LayerAdjacentBoundaryCodes) {
72 double L = EG_LARGE_REAL;
73 vec3_t n0(0, 0, 0);
74 vec3_t x0(0, 0, 0);
75 double total_area = 0;
76 for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
77 if (isSurface(id_cell, m_Grid) && cell_code->GetValue(id_cell) == bc) {
78 vec3_t n = cellNormal(m_Grid, id_cell);
79 double A = n.abs();
80 total_area += A;
81 n0 += n;
82 x0 += A*cellCentre(m_Grid, id_cell);
83 L = min(L, sqrt(4*A/sqrt(3.0)));
86 n0.normalise();
87 x0 *= 1.0/total_area;
88 m_LayerAdjacentNormals[bc] = n0;
89 m_LayerAdjacentOrigins[bc] = x0;
92 computeBoundaryLayerVectors();
93 makeCopy(m_Grid, m_OriginalGrid);
96 QList<vtkIdType> CreateBoundaryLayerShell::correctAdjacentBC(int bc, int num_levels)
98 cout << "correcting boundary \"" << qPrintable(GuiMainWindow::pointer()->getBC(bc).getName()) << "\" with " << num_levels << " levels" << endl;
100 EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
101 EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type");
102 double scal_min = -1;
103 int count = 0;
105 SurfaceMeshSmoother smooth;
106 smooth.setGrid(m_Grid);
107 smooth.useSimpleCentreScheme();
109 // mark num_levels levels of nodes next to the boundary layer
110 QVector<bool> marked(m_Grid->GetNumberOfPoints(), false);
111 QList<vtkIdType> marked_nodes;
112 EG_FORALL_NODES(id_node, m_Grid) {
113 if (m_BoundaryLayerNode[id_node]) {
114 marked[id_node] = true;
115 marked_nodes << id_node;
118 for (int level = 0; level < num_levels; ++level) {
119 QList<vtkIdType> new_marked_nodes;
120 foreach (vtkIdType id_node, marked_nodes) {
121 for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
122 vtkIdType id_neigh = m_Part.n2nGG(id_node, i);
123 if (!marked[id_neigh]) {
124 marked[id_neigh] = true;
125 new_marked_nodes << id_neigh;
129 marked_nodes = new_marked_nodes;
132 QList<vtkIdType> cells;
133 EG_FORALL_CELLS(id_cell, m_Grid) {
134 if (isSurface(id_cell, m_Grid)) {
135 if (cell_code->GetValue(id_cell) == bc) {
136 cells << id_cell;
141 CadInterface *cad = GuiMainWindow::pointer()->getCadInterface(bc);
142 smooth.setCells(cells);
143 smooth.prepareCadInterface(cad);
145 QList<vtkIdType> bad_nodes;
146 int new_num_bad = m_Grid->GetNumberOfPoints();
147 int old_num_bad = 0;
148 updateNodeInfo();
149 do {
150 old_num_bad = new_num_bad;
151 scal_min = 1;
152 bad_nodes.clear();
153 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
154 if (!m_BoundaryLayerNode[id_node]){
155 bool node_has_bc = false;
156 bool all_bcs_adjacent = true;
157 for (int i = 0; i < m_Part.n2bcGSize(id_node); ++i) {
158 int bcn = m_Part.n2bcG(id_node, i);
159 if (bcn == bc) {
160 node_has_bc = true;
162 if (!m_LayerAdjacentBoundaryCodes.contains(bcn)) {
163 all_bcs_adjacent = false;
166 if (node_has_bc && all_bcs_adjacent && marked[id_node] && node_type->GetValue(id_node) != EG_FIXED_VERTEX) {
168 // one boundary condition and EG_SIMPLE_VERTEX
170 if (m_Part.n2bcGSize(id_node) == 1) {
171 vec3_t xs(0, 0, 0); // = smooth.smoothNode(id_node);
172 int count = 0;
173 for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
174 vtkIdType id_cell = m_Part.n2cGG(id_node, i);
175 if (cell_code->GetValue(id_cell) == bc) {
176 vec3_t x = cellCentre(m_Grid, id_cell);
177 xs += x;
178 ++count;
181 if (count > 0) {
182 xs *= 1.0/count;
183 if (node_type->GetValue(id_node) == EG_SIMPLE_VERTEX) {
184 xs = cad->snapNode(id_node, xs);
185 } else {
186 xs = cad->snapToEdge(xs);
188 if (!checkVector(xs)) {
189 EG_ERR_RETURN("error while correcting adjacent boundaries");
191 m_Grid->GetPoints()->SetPoint(id_node, xs.data());
195 // one boundary condition and EG_FEATURE_EDGE_VERTEX
197 else if (m_Part.n2bcGSize(id_node) == 1 && node_type->GetValue(id_node) == EG_FEATURE_EDGE_VERTEX) {
198 int count = 0;
199 vec3_t xs(0, 0, 0);
200 for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
201 vtkIdType id_neigh = m_Part.n2nGG(id_node, i);
202 if (node_type->GetValue(id_neigh) == EG_FEATURE_EDGE_VERTEX) {
203 vec3_t x;
204 m_Grid->GetPoint(id_neigh, x.data());
205 xs += x;
206 ++count;
209 if (count < 2) {
210 vec3_t x;
211 m_Grid->GetPoint(id_node, x.data());
212 xs += x;
213 ++count;
215 xs *= 1.0/count;
216 xs = cad->snapToEdge(xs);
217 m_Grid->GetPoints()->SetPoint(id_node, xs.data());
220 // two boundary conditions and no EG_FIXED_VERTEX
222 else if (m_Part.n2bcGSize(id_node) == 2 && node_type->GetValue(id_node) != EG_FIXED_VERTEX) {
223 int count = 0;
224 vec3_t xs(0, 0, 0);
225 for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
226 vtkIdType id_neigh = m_Part.n2nGG(id_node, i);
227 if (m_Part.n2bcGSize(id_neigh) > 1) {
228 vec3_t x;
229 m_Grid->GetPoint(id_neigh, x.data());
230 xs += x;
231 ++count;
234 if (count < 2) {
235 vec3_t x;
236 m_Grid->GetPoint(id_node, x.data());
237 xs += x;
238 ++count;
240 xs *= 1.0/count;
241 xs = cad->snapToEdge(xs);
242 m_Grid->GetPoints()->SetPoint(id_node, xs.data());
245 bool node_bad = false;
246 for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
247 vtkIdType id_cell = m_Part.n2cGG(id_node, i);
248 if (isSurface(id_cell, m_Grid)) {
249 if (cell_code->GetValue(id_cell) == bc) {
250 CadInterface *cad = GuiMainWindow::pointer()->getCadInterface(cell_code->GetValue(id_cell));
252 // Hier stinkts vielleicht ...
254 cad->snap(cellCentre(m_Grid, id_cell));
255 vec3_t n = cellNormal(m_Grid, id_cell);
256 n.normalise();
257 double scal = n*cad->getLastNormal();
258 scal_min = min(scal_min, scal);
259 if (scal < 0.5 && !node_bad) {
260 bad_nodes << id_node;
261 node_bad = true;
269 new_num_bad = bad_nodes.size();
270 //reduceSurface();
271 //cout << " " << bad_nodes.size() << " node defects" << endl;
272 ++count;
273 } while (scal_min < 0.5 && count < 20);
274 if (scal_min < 0.5) {
275 //writeGrid(grid, "adjacent_bc_failure");
276 //return false;
277 //EG_ERR_RETURN("failed to correct adjacent surfaces");
279 cout << " " << bad_nodes.size() << " node defects" << endl;
280 return bad_nodes;
283 void CreateBoundaryLayerShell::createLayerNodes(vtkIdType id_node)
285 vec3_t x1;
286 m_Grid->GetPoint(id_node, x1.data());
287 m_PrismaticGrid->GetPoints()->SetPoint(m_ShellNodeMap[id_node], x1.data());
288 vec3_t x2 = x1 + m_BoundaryLayerVectors[id_node];
290 double H = m_BoundaryLayerVectors[id_node].abs();
291 double h = H*(1.0 - m_StretchingRatio)/(1.0 - pow(m_StretchingRatio, m_NumLayers));
292 vec3_t dx = (1.0/H)*m_BoundaryLayerVectors[id_node];
293 vec3_t x = x1;
294 m_PrismaticGrid->GetPoints()->SetPoint(m_ShellNodeMap[id_node], x1.data());
295 for (int i = 1; i < m_NumLayers; ++i) {
296 x += h*dx;
297 h *= m_StretchingRatio;
298 m_PrismaticGrid->GetPoints()->SetPoint(i*m_ShellPart.getNumberOfNodes() + m_ShellNodeMap[id_node], x.data());
300 m_PrismaticGrid->GetPoints()->SetPoint(m_NumLayers*m_ShellPart.getNumberOfNodes() + m_ShellNodeMap[id_node], x2.data());
301 m_Grid->GetPoints()->SetPoint(id_node, x2.data());
304 void CreateBoundaryLayerShell::createPrismaticGrid()
306 QVector<vtkIdType> original_triangles, shell_triangles;
307 getSurfaceCells(m_BoundaryLayerCodes, original_triangles, m_OriginalGrid);
308 getSurfaceCells(m_BoundaryLayerCodes, shell_triangles, m_Grid);
310 MeshPartition part(m_Grid);
311 part.setCells(shell_triangles);
312 allocateGrid(m_PrismaticGrid, (m_NumLayers + 1)*part.getNumberOfCells(), (m_NumLayers + 1)*part.getNumberOfNodes());
315 m_ShellNodeMap.fill(-1, m_Grid->GetNumberOfPoints());
316 m_ShellPart.setGrid(m_Grid);
317 m_ShellPart.setCells(shell_triangles);
318 for (int i = 0; i < m_ShellPart.getNumberOfNodes(); ++i) {
319 m_ShellNodeMap[m_ShellPart.globalNode(i)] = i;
322 QVector<QSet<int> > n2bc(m_PrismaticGrid->GetNumberOfPoints());
324 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
325 if (m_ShellNodeMap[id_node] != -1) {
326 createLayerNodes(id_node);
327 for (int i = 0; i < m_Part.n2bcGSize(id_node); ++i) {
328 n2bc[m_ShellNodeMap[id_node]].insert(m_Part.n2bcG(id_node, i));
329 //n2bc[m_ShellNodeMap[id_node] + m_ShellPart.getNumberOfNodes()].insert(m_Part.n2bcG(id_node, i));
334 QList<QVector<vtkIdType> > adjacent_edges;
336 // create prismatic cells and prepare adjacent quad faces
338 foreach (vtkIdType id_cell, shell_triangles) {
339 vtkIdType num_pts, *pts;
340 m_Grid->GetCellPoints(id_cell, num_pts, pts);
341 vtkIdType tri_pts[3], pri_pts[6];
342 for (int i_pts = 0; i_pts < 3; ++i_pts) {
343 if (m_ShellNodeMap[pts[i_pts]] < 0) {
344 EG_BUG;
346 if (m_ShellNodeMap[pts[i_pts]] >= m_ShellPart.getNumberOfNodes()) {
347 EG_BUG;
349 QVector<vtkIdType> edge(4);
350 edge[1] = m_ShellNodeMap[pts[i_pts]];
351 edge[2] = m_ShellNodeMap[pts[0]];
352 if (i_pts < 2) {
353 edge[2] = m_ShellNodeMap[pts[i_pts+1]];
355 QSet<int> edge_codes = m_LayerAdjacentBoundaryCodes;
356 edge_codes.intersect(n2bc[edge[1]]);
357 edge_codes.intersect(n2bc[edge[2]]);
358 if (edge_codes.size() == 1) {
359 edge[0] = *edge_codes.begin();
360 edge[3] = id_cell;
361 adjacent_edges.append(edge);
363 tri_pts[i_pts] = m_ShellNodeMap[pts[i_pts]];
365 vtkIdType id_tri = m_PrismaticGrid->InsertNextCell(VTK_TRIANGLE, 3, tri_pts);
366 copyCellData(m_Grid, id_cell, m_PrismaticGrid, id_tri);
367 for (int i_layer = 0; i_layer < m_NumLayers; ++i_layer) {
368 for (int i_pts = 0; i_pts < 3; ++i_pts) {
369 pri_pts[i_pts] = m_ShellNodeMap[pts[i_pts]] + i_layer*m_ShellPart.getNumberOfNodes();
370 pri_pts[i_pts + 3] = m_ShellNodeMap[pts[i_pts]] + (i_layer + 1)*m_ShellPart.getNumberOfNodes();
372 vtkIdType id_pri = m_PrismaticGrid->InsertNextCell(VTK_WEDGE, 6, pri_pts);
376 // create quads on adjacent boundary faces
378 EG_VTKSP(vtkUnstructuredGrid, noquad_grid);
379 makeCopy(m_PrismaticGrid, noquad_grid);
380 allocateGrid(m_PrismaticGrid, m_PrismaticGrid->GetNumberOfCells() + m_NumLayers*adjacent_edges.size(), m_PrismaticGrid->GetNumberOfPoints());
381 makeCopyNoAlloc(noquad_grid, m_PrismaticGrid);
383 EG_VTKDCC(vtkIntArray, cell_code, m_PrismaticGrid, "cell_code");
385 foreach (QVector<vtkIdType> edge, adjacent_edges) {
386 vtkIdType qua_pts[4];
387 for (int i_layer = 0; i_layer < m_NumLayers; ++i_layer) {
388 qua_pts[0] = edge[2] + i_layer*m_ShellPart.getNumberOfNodes();
389 qua_pts[1] = edge[1] + i_layer*m_ShellPart.getNumberOfNodes();
390 qua_pts[2] = edge[1] + (i_layer + 1)*m_ShellPart.getNumberOfNodes();
391 qua_pts[3] = edge[2] + (i_layer + 1)*m_ShellPart.getNumberOfNodes();
392 vtkIdType id_qua = m_PrismaticGrid->InsertNextCell(VTK_QUAD, 4, qua_pts);
393 copyCellData(m_Grid, edge[3], m_PrismaticGrid, id_qua);
394 cell_code->SetValue(id_qua, edge[0]);
399 void CreateBoundaryLayerShell::reduceSurface()
401 RemovePoints remove_points;
402 MeshPartition part;
403 part.setGrid(m_Grid);
404 part.setAllCells();
405 remove_points.setMeshPartition(part);
406 remove_points.setBoundaryCodes(m_LayerAdjacentBoundaryCodes);
407 remove_points.setUpdatePSPOn();
408 remove_points.setThreshold(3);
409 QVector<bool> fix(m_Grid->GetNumberOfPoints(), false);
410 EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type");
411 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
412 // reset all node types to EG_SIMPLE_VERTEX to trigger a recalculation of the node types
413 node_type->SetValue(id_node, EG_SIMPLE_VERTEX);
414 for (int i = 0; i < part.n2cGSize(id_node); ++i) {
415 if (m_Grid->GetCellType(part.n2cGG(id_node, i)) == VTK_QUAD) {
416 fix[id_node] = true;
420 remove_points.fixNodes(fix);
421 remove_points();
424 void CreateBoundaryLayerShell::smoothSurface()
426 LaplaceSmoother smooth;
427 MeshPartition part;
428 part.setGrid(m_Grid);
429 part.setAllCells();
430 smooth.setMeshPartition(part);
431 smooth.setNumberOfIterations(2);
432 smooth.setBoundaryCodes(m_LayerAdjacentBoundaryCodes);
433 QVector<bool> fix(m_Grid->GetNumberOfPoints(), false);
434 EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type");
435 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
436 // reset all node types to EG_SIMPLE_VERTEX to trigger a recalculation of the node types
437 node_type->SetValue(id_node, EG_SIMPLE_VERTEX);
438 for (int i = 0; i < part.n2cGSize(id_node); ++i) {
439 if (m_Grid->GetCellType(part.n2cGG(id_node, i)) == VTK_QUAD) {
440 fix[id_node] = true;
445 for (int layer = 0; layer < 3; ++layer) {
446 QVector<bool> tmp = fix;
447 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
448 if (!tmp[id_node]) {
449 for (int i = 0; i < part.n2nGSize(id_node); ++i) {
450 fix[part.n2nGG(id_node, i)] = false;
456 smooth.fixNodes(fix);
457 smooth();
461 void CreateBoundaryLayerShell::operate()
463 prepare();
464 writeBoundaryLayerVectors("blayer");
465 createPrismaticGrid();
466 m_Success = true;
467 m_Part.trackGrid(m_Grid);
468 //return;
470 foreach (int bc, m_LayerAdjacentBoundaryCodes) {
472 QList<vtkIdType> bad_nodes;
473 int levels = 1;
475 QVector<vec3_t> x_save(m_Grid->GetNumberOfPoints());
477 do {
478 ++levels;
479 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
480 m_Grid->GetPoint(id_node, x_save[id_node].data());
482 QList<vtkIdType> last_bad = bad_nodes;
483 bad_nodes = correctAdjacentBC(bc, levels);
484 if (bad_nodes.size() >= last_bad.size() && last_bad.size() > 0) {
485 bad_nodes = last_bad;
486 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
487 m_Grid->GetPoints()->SetPoint(id_node, x_save[id_node].data());
489 break;
491 } while (bad_nodes.size() > 0);
493 if (bad_nodes.size() > 0) {
494 bool fixable = true;
496 QVector<bool> is_bad_cell(m_Grid->GetNumberOfCells(), false);
497 foreach (vtkIdType id_node, bad_nodes) {
498 if (m_Part.n2bcGSize(id_node) != 1) {
499 cout << "node " << id_node << " cannot be fixed." << endl;
500 fixable = false;
501 break;
503 for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
504 is_bad_cell[m_Part.n2cGG(id_node, i)] = true;
507 if (fixable) {
508 QList<vtkIdType> bad_cells;
509 EG_FORALL_CELLS (id_cell, m_Grid) {
510 if (is_bad_cell[id_cell]) {
511 bad_cells << id_cell;
515 DeleteCells del;
516 del.setGrid(m_Grid);
517 del.setCellsToDelete(bad_cells);
518 del();
519 StitchHoles stitch(bc);
520 stitch.setGrid(m_Grid);
521 stitch();
522 } else {
523 m_Success = false;
524 cout << "adjacent patch cannot be corrected!" << endl;
529 SwapTriangles swap;
530 swap.setGrid(m_Grid);
531 QSet<int> swap_codes = getAllBoundaryCodes(m_Grid);
532 swap_codes -= m_LayerAdjacentBoundaryCodes;
533 swap.setBoundaryCodes(swap_codes);
534 swap.setVerboseOff();
536 for (int iter = 0; iter < 5; ++iter) {
537 cout << "correcting adjacent boundaries\n" << " iteration: " << iter + 1 << endl;
538 swap();
539 smoothSurface();
540 reduceSurface();
541 swap();