2 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
4 // + This file is part of enGrid. +
6 // + Copyright 2008,2009 Oliver Gloth +
8 // + enGrid is free software: you can redistribute it and/or modify +
9 // + it under the terms of the GNU General Public License as published by +
10 // + the Free Software Foundation, either version 3 of the License, or +
11 // + (at your option) any later version. +
13 // + enGrid is distributed in the hope that it will be useful, +
14 // + but WITHOUT ANY WARRANTY; without even the implied warranty of +
15 // + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +
16 // + GNU General Public License for more details. +
18 // + You should have received a copy of the GNU General Public License +
19 // + along with enGrid. If not, see <http://www.gnu.org/licenses/>. +
21 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
23 #include "guidivideboundarylayer.h"
24 #include "math/linsolve.h"
26 void GuiDivideBoundaryLayer::before()
28 populateBoundaryCodes(ui
.listWidget
);
31 void GuiDivideBoundaryLayer::findBoundaryLayer()
34 insert_cell
.fill(true,grid
->GetNumberOfCells());
37 for (int i_cells
= 0; i_cells
< cells
.size(); ++i_cells
) {
38 if (grid
->GetCellType(cells
[i_cells
]) == VTK_WEDGE
) {
40 vtkIdType N_pts
, *pts
;
41 grid
->GetCellPoints(cells
[i_cells
],N_pts
,pts
);
42 for (int j
= 0; j
< 3; ++j
) {
43 pairs
.insert(QPair
<vtkIdType
,vtkIdType
>(pts
[j
],pts
[j
+3]));
45 for (int j
= 2; j
< 5; ++j
) {
46 if (c2c
[i_cells
][j
] != -1) {
47 vtkIdType type_ncell
= grid
->GetCellType(cells
[c2c
[i_cells
][j
]]);
48 if ((type_ncell
!= VTK_WEDGE
) && (type_ncell
!= VTK_QUAD
)) {
49 EG_ERR_RETURN("unable to identify boundary layer");
55 for (int j
= 0; j
< 2; ++j
) {
56 if (c2c
[i_cells
][j
] != -1) {
57 vtkIdType type_ncell
= grid
->GetCellType(cells
[c2c
[i_cells
][j
]]);
58 if (type_ncell
== VTK_WEDGE
) {
59 EG_ERR_RETURN("the boundary layer seems to have been split already");
66 if (grid
->GetCellType(cells
[i_cells
]) == VTK_QUAD
) {
71 EG_ERR_RETURN("unable to identify boundary layer");
73 is_blayer_node
.clear();
74 is_blayer_node
.fill(false, grid
->GetNumberOfPoints());
76 QPair
<vtkIdType
,vtkIdType
> P
;
78 is_blayer_node
[P
.second
] = true;
82 void GuiDivideBoundaryLayer::findBoundaryLayer1()
87 getSelectedItems(ui
.listWidget
, bcs
);
88 QVector
<vtkIdType
> scells
;
89 getSurfaceCells(bcs
, scells
, grid
);
95 insert_cell
.fill(true,grid
->GetNumberOfCells());
97 foreach (int i_scells
, scells
) {
99 vtkIdType id_scell
= scells
[i_scells
];
100 int i_cells
= findVolumeCell(grid
,id_scell
,_nodes
,cells
,_cells
,n2c
);
101 vtkIdType id_cell
= cells
[i_cells
];
102 vtkIdType type_cell
= grid
->GetCellType(id_cell
);
103 vtkIdType type_scell
= grid
->GetCellType(id_scell
);
105 vtkIdType N_pts
, *pts
;
106 grid
->GetCellPoints(id_cell
,N_pts
,pts
);
108 insert_cell
[id_cell
] = false;
110 if (type_cell
== VTK_WEDGE
) {
115 if (type_scell
!= VTK_QUAD
) {
116 EG_ERR_RETURN("unable to identify boundary layer");
118 if (cells
[c2c
[i_cells
][0]] == id_scell
) {
119 for (int j
= 0; j
< 3; ++j
) {
120 pairs
.insert(QPair
<vtkIdType
,vtkIdType
>(pts
[j
],pts
[j
+3]));
122 } else if (cells
[c2c
[i_cells
][1]] == id_scell
) {
123 for (int j
= 0; j
< 3; ++j
) {
124 pairs
.insert(QPair
<vtkIdType
,vtkIdType
>(pts
[j
+3],pts
[j
]));
127 EG_ERR_RETURN("unable to identify boundary layer");
129 for (int j
= 2; j
< 5; ++j
) {
130 if (c2c
[i_cells
][j
] != -1) {
131 vtkIdType id_ncell
= cells
[c2c
[i_cells
][j
]];
132 vtkIdType type_ncell
= grid
->GetCellType(id_ncell
);
133 if ((type_ncell
!= VTK_WEDGE
) && (type_ncell
!= VTK_HEXAHEDRON
) && (type_ncell
!= VTK_QUAD
)) {
134 EG_ERR_RETURN("unable to identify boundary layer");
136 if (type_ncell
== VTK_QUAD
) {
137 insert_cell
[id_ncell
] = false;
143 } else if (type_cell
== VTK_HEXAHEDRON
) {
148 if (type_scell
!= VTK_QUAD
) {
149 EG_ERR_RETURN("unable to identify boundary layer");
151 QVector
<bool> check_neigh(6,true);
152 if (cells
[c2c
[i_cells
][0]] == id_scell
) {
153 pairs
.insert(QPair
<vtkIdType
,vtkIdType
>(pts
[0],pts
[4]));
154 pairs
.insert(QPair
<vtkIdType
,vtkIdType
>(pts
[1],pts
[5]));
155 pairs
.insert(QPair
<vtkIdType
,vtkIdType
>(pts
[2],pts
[6]));
156 pairs
.insert(QPair
<vtkIdType
,vtkIdType
>(pts
[3],pts
[7]));
157 check_neigh
[0] = false;
158 check_neigh
[1] = false;
159 } else if (cells
[c2c
[i_cells
][1]] == id_scell
) {
160 pairs
.insert(QPair
<vtkIdType
,vtkIdType
>(pts
[4],pts
[0]));
161 pairs
.insert(QPair
<vtkIdType
,vtkIdType
>(pts
[5],pts
[1]));
162 pairs
.insert(QPair
<vtkIdType
,vtkIdType
>(pts
[6],pts
[2]));
163 pairs
.insert(QPair
<vtkIdType
,vtkIdType
>(pts
[7],pts
[3]));
164 check_neigh
[0] = false;
165 check_neigh
[1] = false;
166 } else if (cells
[c2c
[i_cells
][2]] == id_scell
) {
167 pairs
.insert(QPair
<vtkIdType
,vtkIdType
>(pts
[0],pts
[3]));
168 pairs
.insert(QPair
<vtkIdType
,vtkIdType
>(pts
[1],pts
[2]));
169 pairs
.insert(QPair
<vtkIdType
,vtkIdType
>(pts
[5],pts
[6]));
170 pairs
.insert(QPair
<vtkIdType
,vtkIdType
>(pts
[4],pts
[7]));
171 check_neigh
[2] = false;
172 check_neigh
[3] = false;
173 } else if (cells
[c2c
[i_cells
][3]] == id_scell
) {
174 pairs
.insert(QPair
<vtkIdType
,vtkIdType
>(pts
[3],pts
[0]));
175 pairs
.insert(QPair
<vtkIdType
,vtkIdType
>(pts
[2],pts
[1]));
176 pairs
.insert(QPair
<vtkIdType
,vtkIdType
>(pts
[6],pts
[5]));
177 pairs
.insert(QPair
<vtkIdType
,vtkIdType
>(pts
[7],pts
[4]));
178 check_neigh
[2] = false;
179 check_neigh
[3] = false;
180 } else if (cells
[c2c
[i_cells
][4]] == id_scell
) {
181 pairs
.insert(QPair
<vtkIdType
,vtkIdType
>(pts
[0],pts
[1]));
182 pairs
.insert(QPair
<vtkIdType
,vtkIdType
>(pts
[4],pts
[5]));
183 pairs
.insert(QPair
<vtkIdType
,vtkIdType
>(pts
[7],pts
[6]));
184 pairs
.insert(QPair
<vtkIdType
,vtkIdType
>(pts
[3],pts
[2]));
185 check_neigh
[4] = false;
186 check_neigh
[5] = false;
187 } else if (cells
[c2c
[i_cells
][5]] == id_scell
) {
188 pairs
.insert(QPair
<vtkIdType
,vtkIdType
>(pts
[1],pts
[0]));
189 pairs
.insert(QPair
<vtkIdType
,vtkIdType
>(pts
[5],pts
[4]));
190 pairs
.insert(QPair
<vtkIdType
,vtkIdType
>(pts
[6],pts
[7]));
191 pairs
.insert(QPair
<vtkIdType
,vtkIdType
>(pts
[2],pts
[3]));
192 check_neigh
[4] = false;
193 check_neigh
[5] = false;
195 for (int j
= 0; j
< 6; ++j
) {
196 if (check_neigh
[j
]) {
197 if (c2c
[i_cells
][j
] != -1) {
198 vtkIdType id_ncell
= cells
[c2c
[i_cells
][j
]];
199 vtkIdType type_ncell
= grid
->GetCellType(id_ncell
);
200 if ((type_ncell
!= VTK_WEDGE
) && (type_ncell
!= VTK_HEXAHEDRON
) && (type_ncell
!= VTK_QUAD
)) {
201 EG_ERR_RETURN("unable to identify boundary layer");
203 if (type_ncell
== VTK_QUAD
) {
204 insert_cell
[id_ncell
] = false;
212 } else if (type_cell
== VTK_QUAD
) {
219 EG_ERR_RETURN("unable to identify boundary layer");
227 for (int i_cells = 0; i_cells < cells.size(); ++i_cells) {
228 if (grid->GetCellType(cells[i_cells]) == VTK_WEDGE) {
230 vtkIdType N_pts, *pts;
231 grid->GetCellPoints(cells[i_cells],N_pts,pts);
232 for (int j = 0; j < 3; ++j) {
233 pairs.insert(QPair<vtkIdType,vtkIdType>(pts[j],pts[j+3]));
235 for (int j = 2; j < 5; ++j) {
236 if (c2c[i_cells][j] != -1) {
237 vtkIdType type_ncell = grid->GetCellType(cells[c2c[i_cells][j]]);
238 if ((type_ncell != VTK_WEDGE) && (type_ncell != VTK_QUAD)) {
239 EG_ERR_RETURN("unable to identify boundary layer");
245 for (int j = 0; j < 2; ++j) {
246 if (c2c[i_cells][j] != -1) {
247 vtkIdType type_ncell = grid->GetCellType(cells[c2c[i_cells][j]]);
248 if (type_ncell == VTK_WEDGE) {
249 EG_ERR_RETURN("the boundary layer seems to have been split already");
256 if (grid->GetCellType(cells[i_cells]) == VTK_QUAD) {
262 if (N_prisms
+ N_hexes
== 0) {
263 EG_ERR_RETURN("unable to identify boundary layer");
266 is_blayer_node
.clear();
267 is_blayer_node
.fill(false, grid
->GetNumberOfPoints());
269 QPair
<vtkIdType
,vtkIdType
> P
;
271 is_blayer_node
[P
.second
] = true;
275 void GuiDivideBoundaryLayer::bisectF(double &f1
, double &f2
)
279 if (y
[y
.size()-1] > 1) f2
= f
;
285 //begin GROSSER MURKS
287 void GuiDivideBoundaryLayer::computeF()
290 double f2
= 100;//1.0/h;
292 y
[1] = min(0.33,y
[1]);
294 //cout << f1 << ',' << f2 << endl;
296 //while (fabs(1-y[y.size()-1]) > 1e-6);
297 err
= fabs((y
[y
.size()-1]-1)/(y
[y
.size()-1]-y
[y
.size()-2]));
298 } while (err
> 0.01);
302 void GuiDivideBoundaryLayer::computeY()
305 for (int i
= 2; i
< y
.size(); ++i
) {
306 y
[i
] = y
[i
-1] + C
*(y
[i
-1]-y
[i
-2]);
313 void GuiDivideBoundaryLayer::createEdges(vtkUnstructuredGrid
*new_grid
)
315 edges
.fill(QVector
<vtkIdType
>(N_layers
+1), pairs
.size());
316 old2edge
.fill(-1, grid
->GetNumberOfPoints());
318 vtkIdType id_new_node
= grid
->GetNumberOfPoints();
319 QPair
<vtkIdType
,vtkIdType
> P
;
324 edges
[N
][0] = P
.first
;
325 edges
[N
][N_layers
] = P
.second
;
326 old2edge
[P
.first
] = N
;
327 old2edge
[P
.second
] = N
;
330 grid
->GetPoint(P
.first
, x1
.data());
331 grid
->GetPoint(P
.second
, x2
.data());
333 if (!ui
.checkBoxH
->isChecked() || !y_computed
) {
334 y
.resize(N_layers
+ 1);
336 for (int i
= 0; i
< x
.size(); ++i
) {
337 x
[i
] = i
*1.0/(x
.size() - 1);
340 if (ui
.checkBoxH
->isChecked()) {
351 ymin
= min(ymin
,y
[1]*n
.abs());
352 ymax
= max(ymax
,y
[1]*n
.abs());
353 for (int i
= 1; i
< N_layers
; ++i
) {
354 vec3_t x
= x1
+ y
[i
]*n
;
355 max_step
= max(max_step
,(y
[i
+1]-y
[i
])/(y
[i
]-y
[i
-1]));
356 new_grid
->GetPoints()->SetPoint(id_new_node
, x
.data());
357 edges
[N
][i
] = id_new_node
;
363 cout
<< "maximal increment : " << max_step
<< endl
;
364 cout
<< "min(y) : " << ymin
<< endl
;
365 cout
<< "max(y) : " << ymax
<< endl
;
369 void GuiDivideBoundaryLayer::operate1()
372 N_layers
= ui
.spinBoxLayers
->value();
373 h
= ui
.lineEditH
->text().toDouble();
374 F
= ui
.doubleSpinBoxF
->value();
375 cout
<< "dividing boundary layer into " << N_layers
<< " layers:" << endl
;
376 findBoundaryLayer1();
378 EG_VTKSP(vtkUnstructuredGrid
,new_grid
);
379 int N_new_cells
= grid
->GetNumberOfCells() + (N_prisms
+ N_hexes
+ N_quads
)*(N_layers
-1);
380 int N_new_nodes
= grid
->GetNumberOfPoints() + pairs
.size()*(N_layers
-1);
381 allocateGrid(new_grid
, N_new_cells
, N_new_nodes
);
383 // copy existing mesh without prisms and adjacent cells
384 vtkIdType id_new_node
= 0;
385 for (vtkIdType id_node
= 0; id_node
< grid
->GetNumberOfPoints(); ++id_node
) {
387 grid
->GetPoint(id_node
, x
.data());
388 new_grid
->GetPoints()->SetPoint(id_new_node
, x
.data());
389 copyNodeData(grid
, id_node
, new_grid
, id_new_node
);
392 vtkIdType id_new_cell
;
393 for (vtkIdType id_cell
= 0; id_cell
< grid
->GetNumberOfCells(); ++id_cell
) {
394 vtkIdType N_pts
, *pts
;
395 grid
->GetCellPoints(id_cell
, N_pts
, pts
);
396 //bool insert_cell = true;
397 //if (grid->GetCellType(id_cell) == VTK_WEDGE) insert_cell = false;
398 //if (grid->GetCellType(id_cell) == VTK_QUAD) insert_cell = false;
399 if (insert_cell
[id_cell
]) {
400 id_new_cell
= new_grid
->InsertNextCell(grid
->GetCellType(id_cell
), N_pts
, pts
);
401 copyCellData(grid
, id_cell
, new_grid
, id_new_cell
);
405 // create divided boundary layer
406 createEdges(new_grid
);
408 for (vtkIdType id_cell
= 0; id_cell
< grid
->GetNumberOfCells(); ++id_cell
) {
409 if (!insert_cell
[id_cell
]) {
410 if (grid
->GetCellType(id_cell
) == VTK_WEDGE
) {
411 vtkIdType N_pts
, *pts
;
412 grid
->GetCellPoints(id_cell
, N_pts
, pts
);
413 for (int i
= 0; i
< N_layers
; ++i
) {
415 p
[0] = edges
[old2edge
[pts
[0]]][i
];
416 p
[1] = edges
[old2edge
[pts
[1]]][i
];
417 p
[2] = edges
[old2edge
[pts
[2]]][i
];
418 p
[3] = edges
[old2edge
[pts
[0]]][i
+1];
419 p
[4] = edges
[old2edge
[pts
[1]]][i
+1];
420 p
[5] = edges
[old2edge
[pts
[2]]][i
+1];
421 id_new_cell
= new_grid
->InsertNextCell(VTK_WEDGE
, 6, p
);
422 copyCellData(grid
, id_cell
, new_grid
, id_new_cell
);
425 if (grid
->GetCellType(id_cell
) == VTK_HEXAHEDRON
) {
428 vtkIdType N_pts
, *pts
;
429 grid
->GetCellPoints(id_cell
, N_pts
, pts
);
430 for (int i
= 0; i
< N_layers
; ++i
) {
432 p
[0] = edges
[old2edge
[pts
[0]]][i
];
433 p
[1] = edges
[old2edge
[pts
[1]]][i
];
434 p
[2] = edges
[old2edge
[pts
[2]]][i
];
435 p
[3] = edges
[old2edge
[pts
[0]]][i
+1];
436 p
[4] = edges
[old2edge
[pts
[1]]][i
+1];
437 p
[5] = edges
[old2edge
[pts
[2]]][i
+1];
438 id_new_cell
= new_grid
->InsertNextCell(VTK_WEDGE
, 6, p
);
439 copyCellData(grid
, id_cell
, new_grid
, id_new_cell
);
444 if (grid
->GetCellType(id_cell
) == VTK_QUAD
) {
445 vtkIdType N_pts
, *pts
;
446 grid
->GetCellPoints(id_cell
, N_pts
, pts
);
447 if ((old2edge
[pts
[0]] != -1) && (old2edge
[pts
[1]] != -1) && (old2edge
[pts
[2]] != -1) && (old2edge
[pts
[3]] != -1)) {
448 for (int i
= 0; i
< N_layers
; ++i
) {
450 p
[0] = edges
[old2edge
[pts
[0]]][i
];
451 p
[1] = edges
[old2edge
[pts
[1]]][i
];
452 p
[2] = edges
[old2edge
[pts
[1]]][i
+1];
453 p
[3] = edges
[old2edge
[pts
[0]]][i
+1];
454 id_new_cell
= new_grid
->InsertNextCell(VTK_QUAD
, 4, p
);
455 copyCellData(grid
, id_cell
, new_grid
, id_new_cell
);
462 makeCopy(new_grid
, grid
);
465 void GuiDivideBoundaryLayer::operate()
468 N_layers
= ui
.spinBoxLayers
->value();
469 h
= ui
.lineEditH
->text().toDouble();
470 F
= ui
.doubleSpinBoxF
->value();
471 cout
<< "dividing boundary layer into " << N_layers
<< " layers:" << endl
;
474 EG_VTKSP(vtkUnstructuredGrid
,new_grid
);
475 allocateGrid(new_grid
, grid
->GetNumberOfCells() + (N_prisms
+ N_quads
)*(N_layers
-1), grid
->GetNumberOfPoints() + pairs
.size()*(N_layers
-1));
477 // copy existing mesh without prisms and adjacent cells
478 vtkIdType id_new_node
= 0;
479 for (vtkIdType id_node
= 0; id_node
< grid
->GetNumberOfPoints(); ++id_node
) {
481 grid
->GetPoint(id_node
, x
.data());
482 new_grid
->GetPoints()->SetPoint(id_new_node
, x
.data());
483 copyNodeData(grid
, id_node
, new_grid
, id_new_node
);
486 vtkIdType id_new_cell
;
487 for (vtkIdType id_cell
= 0; id_cell
< grid
->GetNumberOfCells(); ++id_cell
) {
488 vtkIdType N_pts
, *pts
;
489 grid
->GetCellPoints(id_cell
, N_pts
, pts
);
490 bool insert_cell
= true;
491 if (grid
->GetCellType(id_cell
) == VTK_WEDGE
) insert_cell
= false;
492 if (grid
->GetCellType(id_cell
) == VTK_QUAD
) insert_cell
= false;
494 id_new_cell
= new_grid
->InsertNextCell(grid
->GetCellType(id_cell
), N_pts
, pts
);
495 copyCellData(grid
, id_cell
, new_grid
, id_new_cell
);
499 // create divided boundary layer
500 createEdges(new_grid
);
502 for (vtkIdType id_cell
= 0; id_cell
< grid
->GetNumberOfCells(); ++id_cell
) {
503 if (grid
->GetCellType(id_cell
) == VTK_WEDGE
) {
504 vtkIdType N_pts
, *pts
;
505 grid
->GetCellPoints(id_cell
, N_pts
, pts
);
506 for (int i
= 0; i
< N_layers
; ++i
) {
508 p
[0] = edges
[old2edge
[pts
[0]]][i
];
509 p
[1] = edges
[old2edge
[pts
[1]]][i
];
510 p
[2] = edges
[old2edge
[pts
[2]]][i
];
511 p
[3] = edges
[old2edge
[pts
[0]]][i
+1];
512 p
[4] = edges
[old2edge
[pts
[1]]][i
+1];
513 p
[5] = edges
[old2edge
[pts
[2]]][i
+1];
514 id_new_cell
= new_grid
->InsertNextCell(VTK_WEDGE
, 6, p
);
515 copyCellData(grid
, id_cell
, new_grid
, id_new_cell
);
518 if (grid
->GetCellType(id_cell
) == VTK_QUAD
) {
519 vtkIdType N_pts
, *pts
;
520 grid
->GetCellPoints(id_cell
, N_pts
, pts
);
521 if ((old2edge
[pts
[0]] != -1) && (old2edge
[pts
[1]] != -1) && (old2edge
[pts
[2]] != -1) && (old2edge
[pts
[3]] != -1)) {
522 for (int i
= 0; i
< N_layers
; ++i
) {
524 p
[0] = edges
[old2edge
[pts
[0]]][i
];
525 p
[1] = edges
[old2edge
[pts
[1]]][i
];
526 p
[2] = edges
[old2edge
[pts
[1]]][i
+1];
527 p
[3] = edges
[old2edge
[pts
[0]]][i
+1];
528 id_new_cell
= new_grid
->InsertNextCell(VTK_QUAD
, 4, p
);
529 copyCellData(grid
, id_cell
, new_grid
, id_new_cell
);
535 makeCopy(new_grid
, grid
);