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 "seedsimpleprismaticlayer.h"
24 #include <vtkIdList.h>
26 SeedSimplePrismaticLayer::SeedSimplePrismaticLayer()
32 void SeedSimplePrismaticLayer::setLayerCells(const QVector
<vtkIdType
> &cells
)
34 layer_cells
.resize(cells
.size());
35 qCopy(cells
.begin(), cells
.end(), layer_cells
.begin());
38 void SeedSimplePrismaticLayer::getLayerCells(QVector
<vtkIdType
> &cells
)
40 cells
.resize(layer_cells
.size());
41 qCopy(layer_cells
.begin(), layer_cells
.end(), cells
.begin());
44 void SeedSimplePrismaticLayer::prepareLayer()
47 EG_VTKSP(vtkIdList
, nds
);
48 EG_VTKSP(vtkIdList
, cls
);
49 EG_VTKDCN(vtkIntArray
, node_layer
, grid
, "node_layer");
50 vol_cells
.fill(-1, layer_cells
.size());
51 faces
.resize(layer_cells
.size());
53 QSet
<vtkIdType
> new_points
;
54 for (int i_layer_cell
= 0; i_layer_cell
< layer_cells
.size(); ++i_layer_cell
) {
55 vtkIdType id_cell
= layer_cells
[i_layer_cell
];
56 vtkIdType type_cell
= grid
->GetCellType(id_cell
);
59 grid
->GetCellPoints(id_cell
, Npts
, pts
);
60 if (type_cell
== VTK_TRIANGLE
) {
62 faces
[i_layer_cell
].resize(Npts
);
63 for (int i_pts
= 0; i_pts
< Npts
; ++i_pts
) {
64 faces
[i_layer_cell
][Npts
- 1 - i_pts
] = pts
[i_pts
];
65 nds
->InsertNextId(pts
[i_pts
]);
66 new_points
.insert(pts
[i_pts
]);
68 grid
->GetCellNeighbors(id_cell
, nds
, cls
);
69 for (int i_cls
= 0; i_cls
< cls
->GetNumberOfIds(); ++i_cls
) {
70 if (isVolume(cls
->GetId(i_cls
), grid
)) {
71 if (cls
->GetId(i_cls
) != id_cell
) {
72 vol_cells
[i_layer_cell
] = cls
->GetId(i_cls
);
77 } else if (type_cell
== VTK_WEDGE
) {
79 faces
[i_layer_cell
].resize(3);
80 for (int i_pts
= 0; i_pts
< 3; ++i_pts
) {
81 faces
[i_layer_cell
][2 - i_pts
] = pts
[i_pts
+3];
82 //faces[i_layer_cell][i_pts] = pts[i_pts+3];
83 nds
->InsertNextId(pts
[i_pts
+3]);
84 new_points
.insert(pts
[i_pts
+3]);
86 grid
->GetCellNeighbors(id_cell
, nds
, cls
);
87 for (int i_cls
= 0; i_cls
< cls
->GetNumberOfIds(); ++i_cls
) {
88 if (isVolume(cls
->GetId(i_cls
), grid
)) {
89 if (cls
->GetId(i_cls
) != id_cell
) {
90 vol_cells
[i_layer_cell
] = cls
->GetId(i_cls
);
99 if (faces
.size() > 0) {
100 old_layer
= node_layer
->GetValue(faces
[0][0]);
101 new_layer
= old_layer
+ 1;
103 N_new_cells
+= countBoundaryElements();
104 N_new_points
= new_points
.size();
107 int SeedSimplePrismaticLayer::countBoundaryElements()
110 QVector
<QSet
<int> > n2f(grid
->GetNumberOfPoints());
111 for (int i_faces
= 0; i_faces
< faces
.size(); ++i_faces
) {
112 for (int j_faces
= 0; j_faces
< faces
[i_faces
].size(); ++j_faces
) {
113 n2f
[faces
[i_faces
][j_faces
]].insert(i_faces
);
116 EG_VTKDCN(vtkIntArray
, node_layer
, grid
, "node_layer");
117 EG_VTKDCN(vtkIntArray
, node_status
, grid
, "node_status");
118 for (int i_faces
= 0; i_faces
< faces
.size(); ++i_faces
) {
119 for (int j_faces
= 0; j_faces
< faces
[i_faces
].size(); ++j_faces
) {
120 vtkIdType p1
= faces
[i_faces
][j_faces
];
121 vtkIdType p2
= faces
[i_faces
][0];
122 if (j_faces
< faces
[i_faces
].size() - 1) {
123 p2
= faces
[i_faces
][j_faces
+1];
125 bool consider_edge
= false;
126 if ((node_status
->GetValue(p1
) & 2) && (node_status
->GetValue(p2
) & 2)) {
127 consider_edge
= true;
129 if ((node_layer
->GetValue(p1
) == 0) && (node_layer
->GetValue(p2
) == 0)) {
130 consider_edge
= true;
133 QSet
<int> faces_on_edge
;
134 setIntersection(n2f
[p1
], n2f
[p2
], faces_on_edge
);
135 if (faces_on_edge
.size() == 0) {
137 } else if (faces_on_edge
.size() == 1) {
139 } else if (faces_on_edge
.size() > 2) {
148 void SeedSimplePrismaticLayer::createBoundaryElements(vtkUnstructuredGrid
*new_grid
)
150 QVector
<QSet
<int> > n2f(grid
->GetNumberOfPoints());
151 for (int i_faces
= 0; i_faces
< faces
.size(); ++i_faces
) {
152 for (int j_faces
= 0; j_faces
< faces
[i_faces
].size(); ++j_faces
) {
153 n2f
[faces
[i_faces
][j_faces
]].insert(i_faces
);
156 QVector
<vtkIdType
> bcells
;
157 QVector
<vtkIdType
> bnodes
;
158 QVector
<int> _bnodes
;
159 QVector
<QSet
<int> > bn2bc
;
160 getAllSurfaceCells(bcells
, new_grid
);
161 getNodesFromCells(bcells
, bnodes
, new_grid
);
162 createNodeMapping(bnodes
, _bnodes
, new_grid
);
163 createNodeToCell(bcells
, bnodes
, _bnodes
, bn2bc
, new_grid
);
164 EG_VTKDCC(vtkIntArray
, cell_code
, new_grid
, "cell_code");
165 EG_VTKDCN(vtkIntArray
, node_layer
, new_grid
, "node_layer");
166 EG_VTKDCN(vtkIntArray
, node_status
, new_grid
, "node_status");
167 for (int i_faces
= 0; i_faces
< faces
.size(); ++i_faces
) {
168 for (int j_faces
= 0; j_faces
< faces
[i_faces
].size(); ++j_faces
) {
169 vtkIdType p1
= faces
[i_faces
][j_faces
];
170 vtkIdType p2
= faces
[i_faces
][0];
171 if (j_faces
< faces
[i_faces
].size() - 1) {
172 p2
= faces
[i_faces
][j_faces
+1];
174 bool consider_edge
= false;
175 if ((node_status
->GetValue(p1
) & 2) && (node_status
->GetValue(p2
) & 2)) {
176 consider_edge
= true;
178 if ((node_layer
->GetValue(p1
) == 0) && (node_layer
->GetValue(p2
) == 0)) {
179 consider_edge
= true;
182 QSet
<int> faces_on_edge
;
183 setIntersection(n2f
[p1
], n2f
[p2
], faces_on_edge
);
184 if (faces_on_edge
.size() == 0) {
186 } else if (faces_on_edge
.size() == 1) {
187 vtkIdType pts
[4] = { p1
, p2
, old2new
[p2
], old2new
[p1
] };
188 vtkIdType id_new_bcell
= new_grid
->InsertNextCell(VTK_QUAD
,4, pts
);
190 QSet
<int> bc1
, bc2
, bc3
;
191 if (_bnodes
[old2new
[p1
]] != -1) {
192 foreach (int i_bcells
, bn2bc
[_bnodes
[old2new
[p1
]]]) {
193 if (!boundary_codes
.contains(cell_code
->GetValue(bcells
[i_bcells
]))) {
194 bc1
.insert(cell_code
->GetValue(bcells
[i_bcells
]));
198 if (_bnodes
[old2new
[p2
]] != -1) {
199 foreach (int i_bcells
, bn2bc
[_bnodes
[old2new
[p2
]]]) {
200 if (!boundary_codes
.contains(cell_code
->GetValue(bcells
[i_bcells
]))) {
201 bc2
.insert(cell_code
->GetValue(bcells
[i_bcells
]));
205 setIntersection(bc1
, bc2
, bc3
);
206 if (bc3
.size() == 1) {
212 cell_code
->SetValue(id_new_bcell
, bc
);
213 node_status
->SetValue(p1
, node_status
->GetValue(p1
) | 2);
214 node_status
->SetValue(old2new
[p1
], node_status
->GetValue(p1
) | 2);
215 node_status
->SetValue(p2
, node_status
->GetValue(p2
) | 2);
216 node_status
->SetValue(old2new
[p2
], node_status
->GetValue(p2
) | 2);
217 } else if (faces_on_edge
.size() > 2) {
225 void SeedSimplePrismaticLayer::operate()
227 cout
<< "seeding prismatic layer" << endl
;
229 EG_VTKSP(vtkUnstructuredGrid
, new_grid
);
230 allocateGrid(new_grid
, grid
->GetNumberOfCells() + N_new_cells
, grid
->GetNumberOfPoints() + N_new_points
);
231 vtkIdType id_new_node
= 0;
232 EG_VTKDCN(vtkIntArray
, node_status_old
, grid
, "node_status");
233 EG_VTKDCN(vtkIntArray
, node_status_new
, new_grid
, "node_status");
234 EG_VTKDCN(vtkIntArray
, node_layer_old
, grid
, "node_layer");
235 EG_VTKDCN(vtkIntArray
, node_layer_new
, new_grid
, "node_layer");
236 EG_VTKDCC(vtkIntArray
, bc
, grid
, "cell_code");
238 QVector
<QSet
<int> > n2f(grid
->GetNumberOfPoints());
239 for (int i_faces
= 0; i_faces
< faces
.size(); ++i_faces
) {
240 for (int j_faces
= 0; j_faces
< faces
[i_faces
].size(); ++j_faces
) {
241 n2f
[faces
[i_faces
][j_faces
]].insert(i_faces
);
245 // copy old grid nodes to the new grid
246 old2new
.resize(grid
->GetNumberOfPoints());
247 for (vtkIdType id_node
= 0; id_node
< grid
->GetNumberOfPoints(); ++id_node
) {
249 grid
->GetPoints()->GetPoint(id_node
, x
.data());
250 new_grid
->GetPoints()->SetPoint(id_new_node
, x
.data());
251 copyNodeData(grid
, id_node
, new_grid
, id_new_node
);
252 old2new
[id_node
] = id_new_node
;
256 QSet
<vtkIdType
> split_nodes
;
257 for (int i_faces
= 0; i_faces
< faces
.size(); ++i_faces
) {
258 for (int j_faces
= 0; j_faces
< faces
[i_faces
].size(); ++j_faces
) {
259 split_nodes
.insert(faces
[i_faces
][j_faces
]);
262 qDebug() << split_nodes
.size();
263 QVector
<bool> is_split_node(grid
->GetNumberOfPoints(), false);
264 foreach (vtkIdType id_node
, split_nodes
) {
266 grid
->GetPoints()->GetPoint(id_node
, x
.data());
269 if (n2f
[id_node
].size() > 0) {
272 foreach (int i_faces
, n2f
[id_node
]) {
274 if (faces
[i_faces
][0] == id_node
) {
275 grid
->GetPoint(faces
[i_faces
][1],a
.data());
276 grid
->GetPoint(faces
[i_faces
][2],b
.data());
278 if (faces
[i_faces
][1] == id_node
) {
279 grid
->GetPoint(faces
[i_faces
][2],a
.data());
280 grid
->GetPoint(faces
[i_faces
][0],b
.data());
282 if (faces
[i_faces
][2] == id_node
) {
283 grid
->GetPoint(faces
[i_faces
][0],a
.data());
284 grid
->GetPoint(faces
[i_faces
][1],b
.data());
292 vec3_t nf
= a
.cross(b
);
294 double alpha
= acos(a
*b
);
297 for (int i_boundary_correction
= 0; i_boundary_correction
< 20; ++i_boundary_correction
) {
298 foreach (vtkIdType id_cell
, n2c
[_nodes
[id_node
]]) {
299 if (isSurface(id_cell
, grid
)) {
300 if (!boundary_codes
.contains(bc
->GetValue(id_cell
))) {
301 double A
= GeometryTools::cellVA(grid
, id_cell
);
303 vec3_t nf
= GeometryTools::cellNormal(grid
, id_cell
);
317 new_grid
->GetPoints()->SetPoint(id_new_node
, x
.data());
318 old2new
[id_node
] = id_new_node
;
319 node_status_new
->SetValue(id_new_node
, node_status_old
->GetValue(id_node
));
320 node_layer_new
->SetValue(id_new_node
, node_layer_old
->GetValue(id_node
) + 1);
322 is_split_node
[id_node
] = true;
324 QVector
<bool> needs_correction(grid
->GetNumberOfCells(), false);
325 for (vtkIdType id_cell
= 0; id_cell
< grid
->GetNumberOfCells(); ++id_cell
) {
326 vtkIdType type_cell
= grid
->GetCellType(id_cell
);
327 if ((type_cell
== VTK_TRIANGLE
) || (type_cell
== VTK_TETRA
)) {
328 vtkIdType
*pts
, N_pts
;
329 grid
->GetCellPoints(id_cell
, N_pts
, pts
);
331 for (int i_pts
= 0; i_pts
< N_pts
; ++i_pts
) {
332 if (is_split_node
[pts
[i_pts
]]) {
334 if (node_layer_old
->GetValue(pts
[i_pts
]) != old_layer
) {
340 if (type_cell
== VTK_TETRA
) {
341 needs_correction
[id_cell
] = true;
345 for (int i_pts
= 0; i_pts
< N_pts
; ++i_pts
) {
346 if (node_layer_old
->GetValue(pts
[i_pts
]) == old_layer
- 1) {
352 needs_correction
[id_cell
] = true;
354 cout
<< "dodgy face: " << id_cell
<< " ";
355 for (int i_pts
= 0; i_pts
< N_pts
; ++i_pts
) {
356 cout
<< "(" << pts
[i_pts
] << "," << node_layer_old
->GetValue(pts
[i_pts
]) << ") ";
359 for (int i_pts
= 0; i_pts
< N_pts
; ++i_pts
) {
360 if (is_split_node
[pts
[i_pts
]]) {
362 grid
->GetPoint(pts
[i_pts
], x
.data());
363 cout
<< "split node: " << pts
[i_pts
] << " " << x
<< endl
;
371 foreach (vtkIdType id_cell
, layer_cells
) {
372 needs_correction
[id_cell
] = false;
375 QVector
<bool> is_in_vol_cells(grid
->GetNumberOfCells(), false);
376 for (int i_faces
= 0; i_faces
< faces
.size(); ++i_faces
) {
377 vtkIdType id_vol_cell
= vol_cells
[i_faces
];
378 if (id_vol_cell
>= 0) {
379 is_in_vol_cells
[id_vol_cell
] = true;
380 vtkIdType type_vol_cell
= grid
->GetCellType(id_vol_cell
);
381 if (type_vol_cell
== VTK_TETRA
) {
382 vtkIdType
*pts
, N_pts
, p
[6];
383 grid
->GetCellPoints(id_vol_cell
, N_pts
, pts
);
384 p
[0] = faces
[i_faces
][0];
385 p
[1] = faces
[i_faces
][1];
386 p
[2] = faces
[i_faces
][2];
387 p
[3] = old2new
[p
[0]];
388 p
[4] = old2new
[p
[1]];
389 p
[5] = old2new
[p
[2]];
391 vtkIdType pts
[6] = { p
[2], p
[1], p
[0], p
[5], p
[4], p
[3] };
392 layer_cells
[i_faces
] = new_grid
->InsertNextCell(VTK_WEDGE
, 6, pts
);
395 qDebug() << type_vol_cell
;
401 // copy old cells to the new grid
402 for (vtkIdType id_cell
= 0; id_cell
< grid
->GetNumberOfCells(); ++id_cell
) {
403 vtkIdType
*pts
, N_pts
;
404 grid
->GetCellPoints(id_cell
, N_pts
, pts
);
405 vtkIdType type_cell
= grid
->GetCellType(id_cell
);
406 if (needs_correction
[id_cell
]) {
407 for (int i_pts
= 0; i_pts
< N_pts
; ++i_pts
) {
408 pts
[i_pts
] = old2new
[pts
[i_pts
]];
411 vtkIdType id_new_cell
= new_grid
->InsertNextCell(type_cell
, N_pts
, pts
);
412 copyCellData(grid
, id_cell
, new_grid
, id_new_cell
);
415 createBoundaryElements(new_grid
);
416 UpdateCellIndex(new_grid
);
417 grid
->DeepCopy(new_grid
);
418 cout
<< "done." << endl
;