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