?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
fairly stable now
[engrid-github.git] / src / libengrid / createboundarylayershell.cpp
blob6131d313db82e99b6f70c47d3e7ae45734d73df4
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 "createboundarylayershell.h"
23
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
32 CreateBoundaryLayerShell::CreateBoundaryLayerShell()
33 {
34 m_RestGrid = vtkSmartPointer<vtkUnstructuredGrid>::New();
35 m_OriginalGrid = vtkSmartPointer<vtkUnstructuredGrid>::New();
36 m_PrismaticGrid = vtkSmartPointer<vtkUnstructuredGrid>::New();
37 }
38
39 void CreateBoundaryLayerShell::prepare()
40 {
41 m_Part.trackGrid(m_Grid);
42
43 DeleteVolumeGrid delete_volume;
44 delete_volume.setGrid(m_Grid);
45 delete_volume.setAllCells();
46 delete_volume();
47
48 readSettings();
49 setAllCells();
50 getSurfaceCells(m_BoundaryLayerCodes, layer_cells, m_Grid);
51
52 // fill m_LayerAdjacentBoundaryCodes
53 EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
54 foreach (vtkIdType id_cell, layer_cells) {
55 for (int i = 0; i < m_Part.c2cGSize(id_cell); ++i) {
56 vtkIdType id_neigh = m_Part.c2cGG(id_cell, i);
57 int bc = cell_code->GetValue(id_neigh);
58 if (!m_BoundaryLayerCodes.contains(bc)) {
59 m_LayerAdjacentBoundaryCodes.insert(bc);
60 }
61 }
62 }
63
64 // compute normals and origins of adjacent planes
65 m_LayerAdjacentNormals.clear();
66 m_LayerAdjacentOrigins.clear();
67 foreach (int bc, m_LayerAdjacentBoundaryCodes) {
68 double L = EG_LARGE_REAL;
69 vec3_t n0(0, 0, 0);
70 vec3_t x0(0, 0, 0);
71 double total_area = 0;
72 for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
73 if (isSurface(id_cell, m_Grid) && cell_code->GetValue(id_cell) == bc) {
74 vec3_t n = cellNormal(m_Grid, id_cell);
75 double A = n.abs();
76 total_area += A;
77 n0 += n;
78 x0 += A*cellCentre(m_Grid, id_cell);
79 L = min(L, sqrt(4*A/sqrt(3.0)));
80 }
81 }
82 n0.normalise();
83 x0 *= 1.0/total_area;
84 /*
85 for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
86 if (isSurface(id_cell, m_Grid) && cell_code->GetValue(id_cell) == bc) {
87 vec3_t x = cellCentre(m_Grid, id_cell);
88 double l = fabs((x - x0)*n0);
89 if (l > 0.1*L) {
90 BoundaryCondition boundary_condition = GuiMainWindow::pointer()->getBC(bc);
91 QString err_msg = "The boundary \"" + boundary_condition.getName() + "\" is not planar.\n";
92 QString L_txt, l_txt;
93 L_txt.setNum(L);
94 l_txt.setNum(l);
95 err_msg += "L = " + L_txt + " , l = " + l_txt;
96 EG_ERR_RETURN(err_msg);
97 }
98 }
99 }*/
100 m_LayerAdjacentNormals[bc] = n0;
101 m_LayerAdjacentOrigins[bc] = x0;
102 }
103
104 computeBoundaryLayerVectors();
105 makeCopy(m_Grid, m_OriginalGrid);
106 }
107
108 void CreateBoundaryLayerShell::correctAdjacentBC(int bc, vtkUnstructuredGrid *grid)
109 {
110 EG_VTKDCC(vtkIntArray, cell_code, grid, "cell_code");
111 MeshPartition part(grid, true);
112 vec3_t n0 = m_LayerAdjacentNormals[bc];
113 vec3_t x0 = m_LayerAdjacentOrigins[bc];
114 double scal_min = -1;
115 int count = 0;
116 while (scal_min < 0.5 && count < 1000) {
117 scal_min = 1;
118 for (vtkIdType id_node = 0; id_node < grid->GetNumberOfPoints(); ++id_node) {
119 if (part.n2bcGSize(id_node) == 1) {
120 if (part.n2bcG(id_node, 0) == bc) {
121 vec3_t x(0,0,0);
122 int count = 0;
123 double A_total = 0;
124 for (int i = 0; i < part.n2cGSize(id_node); ++i) {
125 vtkIdType id_cell = part.n2cGG(id_node, i);
126 if (isSurface(id_cell, grid)) {
127 if (cell_code->GetValue(id_cell) == bc) {
128 vec3_t n = cellNormal(grid, id_cell);
129 double A = n.abs();
130 n.normalise();
131 if (n*n0 > 0.5) {
132 x += A*cellCentre(grid, id_cell);
133 A_total += A;
134 ++count;
135 }
136 }
137 }
138 }
139 if (count == 0) {
140 grid->GetPoint(id_node, x.data());
141 } else {
142 x *= 1.0/A_total;
143 }
144
145 x -= ((x - x0)*n0)*n0;
146 grid->GetPoints()->SetPoint(id_node, x.data());
147 for (int i = 0; i < part.n2cGSize(id_node); ++i) {
148 vtkIdType id_cell = part.n2cGG(id_node, i);
149 if (isSurface(id_cell, grid)) {
150 if (cell_code->GetValue(id_cell) == bc) {
151 vec3_t n = cellNormal(grid, id_cell);
152 n.normalise();
153 scal_min = min(scal_min, n*n0);
154 }
155 }
156 }
157 }
158 }
159 }
160 ++count;
161 }
162 if (scal_min < 0.5) {
163 writeGrid(grid, "adjacent_bc_failure");
164 EG_ERR_RETURN("failed to correct adjacent surfaces");
165 }
166 }
167
168 void CreateBoundaryLayerShell::createLayerNodes(vtkIdType id_node)
169 {
170 vec3_t x1;
171 m_Grid->GetPoint(id_node, x1.data());
172 m_PrismaticGrid->GetPoints()->SetPoint(m_ShellNodeMap[id_node], x1.data());
173 vec3_t x2 = x1 + m_BoundaryLayerVectors[id_node];
174
175 double H = m_BoundaryLayerVectors[id_node].abs();
176 double h = H*(1.0 - m_StretchingRatio)/(1.0 - pow(m_StretchingRatio, m_NumLayers));
177 vec3_t dx = (1.0/H)*m_BoundaryLayerVectors[id_node];
178 vec3_t x = x1;
179 m_PrismaticGrid->GetPoints()->SetPoint(m_ShellNodeMap[id_node], x1.data());
180 for (int i = 1; i < m_NumLayers; ++i) {
181 x += h*dx;
182 h *= m_StretchingRatio;
183 m_PrismaticGrid->GetPoints()->SetPoint(i*m_ShellPart.getNumberOfNodes() + m_ShellNodeMap[id_node], x.data());
184 }
185 m_PrismaticGrid->GetPoints()->SetPoint(m_NumLayers*m_ShellPart.getNumberOfNodes() + m_ShellNodeMap[id_node], x2.data());
186 m_Grid->GetPoints()->SetPoint(id_node, x2.data());
187 }
188
189 void CreateBoundaryLayerShell::createPrismaticGrid()
190 {
191 QVector<vtkIdType> original_triangles, shell_triangles;
192 getSurfaceCells(m_BoundaryLayerCodes, original_triangles, m_OriginalGrid);
193 getSurfaceCells(m_BoundaryLayerCodes, shell_triangles, m_Grid);
194 {
195 MeshPartition part(m_Grid);
196 part.setCells(shell_triangles);
197 allocateGrid(m_PrismaticGrid, (m_NumLayers + 1)*part.getNumberOfCells(), (m_NumLayers + 1)*part.getNumberOfNodes());
198 }
199
200 m_ShellNodeMap.fill(-1, m_Grid->GetNumberOfPoints());
201 m_ShellPart.setGrid(m_Grid);
202 m_ShellPart.setCells(shell_triangles);
203 for (int i = 0; i < m_ShellPart.getNumberOfNodes(); ++i) {
204 m_ShellNodeMap[m_ShellPart.globalNode(i)] = i;
205 }
206
207 QVector<QSet<int> > n2bc(m_PrismaticGrid->GetNumberOfPoints());
208
209 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
210 if (m_ShellNodeMap[id_node] != -1) {
211 createLayerNodes(id_node);
212 for (int i = 0; i < m_Part.n2bcGSize(id_node); ++i) {
213 n2bc[m_ShellNodeMap[id_node]].insert(m_Part.n2bcG(id_node, i));
214 //n2bc[m_ShellNodeMap[id_node] + m_ShellPart.getNumberOfNodes()].insert(m_Part.n2bcG(id_node, i));
215 }
216 }
217 }
218
219 QList<QVector<vtkIdType> > adjacent_edges;
220
221 // create prismatic cells and prepare adjacent quad faces
222 //
223 foreach (vtkIdType id_cell, shell_triangles) {
224 vtkIdType num_pts, *pts;
225 m_Grid->GetCellPoints(id_cell, num_pts, pts);
226 vtkIdType tri_pts[3], pri_pts[6];
227 for (int i_pts = 0; i_pts < 3; ++i_pts) {
228 if (m_ShellNodeMap[pts[i_pts]] < 0) {
229 EG_BUG;
230 }
231 if (m_ShellNodeMap[pts[i_pts]] >= m_ShellPart.getNumberOfNodes()) {
232 EG_BUG;
233 }
234 QVector<vtkIdType> edge(4);
235 edge[1] = m_ShellNodeMap[pts[i_pts]];
236 edge[2] = m_ShellNodeMap[pts[0]];
237 if (i_pts < 2) {
238 edge[2] = m_ShellNodeMap[pts[i_pts+1]];
239 }
240 QSet<int> edge_codes = m_LayerAdjacentBoundaryCodes;
241 edge_codes.intersect(n2bc[edge[1]]);
242 edge_codes.intersect(n2bc[edge[2]]);
243 if (edge_codes.size() == 1) {
244 edge[0] = *edge_codes.begin();
245 edge[3] = id_cell;
246 adjacent_edges.append(edge);
247 }
248 tri_pts[i_pts] = m_ShellNodeMap[pts[i_pts]];
249 }
250 vtkIdType id_tri = m_PrismaticGrid->InsertNextCell(VTK_TRIANGLE, 3, tri_pts);
251 copyCellData(m_Grid, id_cell, m_PrismaticGrid, id_tri);
252 for (int i_layer = 0; i_layer < m_NumLayers; ++i_layer) {
253 for (int i_pts = 0; i_pts < 3; ++i_pts) {
254 pri_pts[i_pts] = m_ShellNodeMap[pts[i_pts]] + i_layer*m_ShellPart.getNumberOfNodes();
255 pri_pts[i_pts + 3] = m_ShellNodeMap[pts[i_pts]] + (i_layer + 1)*m_ShellPart.getNumberOfNodes();
256 }
257 vtkIdType id_pri = m_PrismaticGrid->InsertNextCell(VTK_WEDGE, 6, pri_pts);
258 }
259 }
260
261 // create quads on adjacent boundary faces
262 //
263 EG_VTKSP(vtkUnstructuredGrid, noquad_grid);
264 makeCopy(m_PrismaticGrid, noquad_grid);
265 allocateGrid(m_PrismaticGrid, m_PrismaticGrid->GetNumberOfCells() + m_NumLayers*adjacent_edges.size(), m_PrismaticGrid->GetNumberOfPoints());
266 makeCopyNoAlloc(noquad_grid, m_PrismaticGrid);
267
268 EG_VTKDCC(vtkIntArray, cell_code, m_PrismaticGrid, "cell_code");
269
270 foreach (QVector<vtkIdType> edge, adjacent_edges) {
271 vtkIdType qua_pts[4];
272 for (int i_layer = 0; i_layer < m_NumLayers; ++i_layer) {
273 qua_pts[0] = edge[2] + i_layer*m_ShellPart.getNumberOfNodes();
274 qua_pts[1] = edge[1] + i_layer*m_ShellPart.getNumberOfNodes();
275 qua_pts[2] = edge[1] + (i_layer + 1)*m_ShellPart.getNumberOfNodes();
276 qua_pts[3] = edge[2] + (i_layer + 1)*m_ShellPart.getNumberOfNodes();
277 vtkIdType id_qua = m_PrismaticGrid->InsertNextCell(VTK_QUAD, 4, qua_pts);
278 copyCellData(m_Grid, edge[3], m_PrismaticGrid, id_qua);
279 cell_code->SetValue(id_qua, edge[0]);
280 }
281 }
282 }
283
284 void CreateBoundaryLayerShell::reduceSurface()
285 {
286 RemovePoints remove_points;
287 MeshPartition part;
288 part.setGrid(m_Grid);
289 part.setAllCells();
290 remove_points.setMeshPartition(part);
291 remove_points.setBoundaryCodes(m_LayerAdjacentBoundaryCodes);
292 remove_points.setUpdatePSPOn();
293 remove_points.setThreshold(3);
294 QVector<bool> fix(m_Grid->GetNumberOfPoints(), true);
295
296 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
297 for (int i = 0; i < part.n2cGSize(id_node); ++i) {
298 if (m_Grid->GetCellType(part.n2cGG(id_node, i)) == VTK_WEDGE) {
299 fix[id_node] = false;
300 }
301 }
302 }
303 for (int layer = 0; layer < 3; ++layer) {
304 QVector<bool> tmp = fix;
305 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
306 if (!tmp[id_node]) {
307 for (int i = 0; i < part.n2nGSize(id_node); ++i) {
308 fix[part.n2nGG(id_node, i)] = false;
309 }
310 }
311 }
312 }
313 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
314 for (int i = 0; i < part.n2cGSize(id_node); ++i) {
315 if (m_Grid->GetCellType(part.n2cGG(id_node, i)) == VTK_WEDGE) {
316 fix[id_node] = true;
317 }
318 }
319 }
320
321 remove_points.fixNodes(fix);
322 remove_points();
323 }
324
325 void CreateBoundaryLayerShell::smoothSurface()
326 {
327 LaplaceSmoother smooth;
328 MeshPartition part;
329 part.setGrid(m_Grid);
330 part.setAllCells();
331 smooth.setMeshPartition(part);
332 smooth.setNumberOfIterations(2);
333 smooth.setBoundaryCodes(m_LayerAdjacentBoundaryCodes);
334
335 QVector<bool> fix(m_Grid->GetNumberOfPoints(), true);
336 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
337 for (int i = 0; i < part.n2cGSize(id_node); ++i) {
338 if (m_Grid->GetCellType(part.n2cGG(id_node, i)) == VTK_WEDGE) {
339 fix[id_node] = false;
340 }
341 }
342 }
343 for (int layer = 0; layer < 3; ++layer) {
344 QVector<bool> tmp = fix;
345 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
346 if (!tmp[id_node]) {
347 for (int i = 0; i < part.n2nGSize(id_node); ++i) {
348 fix[part.n2nGG(id_node, i)] = false;
349 }
350 }
351 }
352 }
353
354 smooth.fixNodes(fix);
355 smooth();
356 }
357
358
359 void CreateBoundaryLayerShell::operate()
360 {
361 prepare();
362 writeBoundaryLayerVectors("blayer");
363 createPrismaticGrid();
364 foreach (int bc, m_LayerAdjacentBoundaryCodes) {
365 correctAdjacentBC(bc, m_Grid);
366 }
367 SwapTriangles swap;
368 swap.setGrid(m_Grid);
369 QSet<int> swap_codes = getAllBoundaryCodes(m_Grid);
370 swap_codes -= m_LayerAdjacentBoundaryCodes;
371 swap.setBoundaryCodes(swap_codes);
372 swap.setVerboseOff();
373
374 for (int iter = 0; iter < 1; ++iter) {
375 swap();
376 smoothSurface();
377 reduceSurface();
378 swap();
379 }
380
381 }
382
This is a mirror repository of the new location for engrid@github