1 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3 // + This file is part of enGrid. +
5 // + Copyright 2008-2014 enGits GmbH +
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. +
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. +
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/>. +
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();
43 void CreateBoundaryLayerShell::prepare()
45 m_Part
.trackGrid(m_Grid
);
47 DeleteVolumeGrid delete_volume
;
48 delete_volume
.setGrid(m_Grid
);
49 delete_volume
.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
;
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
);
82 x0
+= A
*cellCentre(m_Grid
, id_cell
);
83 L
= min(L
, sqrt(4*A
/sqrt(3.0)));
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;
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
) {
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();
150 old_num_bad
= new_num_bad
;
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
);
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);
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
);
183 if (node_type
->GetValue(id_node
) == EG_SIMPLE_VERTEX
) {
184 xs
= cad
->snapNode(id_node
, xs
);
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
) {
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
) {
204 m_Grid
->GetPoint(id_neigh
, x
.data());
211 m_Grid
->GetPoint(id_node
, x
.data());
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
) {
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) {
229 m_Grid
->GetPoint(id_neigh
, x
.data());
236 m_Grid
->GetPoint(id_node
, x
.data());
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
);
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
;
269 new_num_bad
= bad_nodes
.size();
271 //cout << " " << bad_nodes.size() << " node defects" << endl;
273 } while (scal_min
< 0.5 && count
< 20);
274 if (scal_min
< 0.5) {
275 //writeGrid(grid, "adjacent_bc_failure");
277 //EG_ERR_RETURN("failed to correct adjacent surfaces");
279 cout
<< " " << bad_nodes
.size() << " node defects" << endl
;
283 void CreateBoundaryLayerShell::createLayerNodes(vtkIdType id_node
)
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
];
294 m_PrismaticGrid
->GetPoints()->SetPoint(m_ShellNodeMap
[id_node
], x1
.data());
295 for (int i
= 1; i
< m_NumLayers
; ++i
) {
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) {
346 if (m_ShellNodeMap
[pts
[i_pts
]] >= m_ShellPart
.getNumberOfNodes()) {
349 QVector
<vtkIdType
> edge(4);
350 edge
[1] = m_ShellNodeMap
[pts
[i_pts
]];
351 edge
[2] = m_ShellNodeMap
[pts
[0]];
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();
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
;
403 part
.setGrid(m_Grid
);
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
) {
420 remove_points
.fixNodes(fix
);
424 void CreateBoundaryLayerShell::smoothSurface()
426 LaplaceSmoother smooth
;
428 part
.setGrid(m_Grid
);
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
) {
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) {
449 for (int i = 0; i < part.n2nGSize(id_node); ++i) {
450 fix[part.n2nGG(id_node, i)] = false;
456 smooth
.fixNodes(fix
);
461 void CreateBoundaryLayerShell::operate()
464 writeBoundaryLayerVectors("blayer");
465 createPrismaticGrid();
467 m_Part
.trackGrid(m_Grid
);
470 foreach (int bc
, m_LayerAdjacentBoundaryCodes
) {
472 QList
<vtkIdType
> bad_nodes
;
475 QVector
<vec3_t
> x_save(m_Grid
->GetNumberOfPoints());
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());
491 } while (bad_nodes
.size() > 0);
493 if (bad_nodes
.size() > 0) {
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
;
503 for (int i
= 0; i
< m_Part
.n2cGSize(id_node
); ++i
) {
504 is_bad_cell
[m_Part
.n2cGG(id_node
, i
)] = true;
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
;
517 del
.setCellsToDelete(bad_cells
);
519 StitchHoles
stitch(bc
);
520 stitch
.setGrid(m_Grid
);
524 cout
<< "adjacent patch cannot be corrected!" << endl
;
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
;