?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
limited volume meshing to boundary layer only
[engrid-github.git] / src / libengrid / egvtkobject.cpp
blobf1c1d7446ea1cdcbd3ae31f0b37aebe34cc8eb8c
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 #include "egvtkobject.h"
22 #include "engrid.h"
23 #include "guimainwindow.h"
24
25 #include <vtkCellData.h>
26 #include <vtkPointData.h>
27 #include <vtkCellLinks.h>
28 #include <vtkCellType.h>
29 #include <vtkIdList.h>
30 #include <vtkCell.h>
31 #include <vtkSignedCharArray.h>
32 #include <vtkCellArray.h>
33 #include <vtkSmartPointer.h>
34 #include <vtkType.h>
35
36 int EgVtkObject::DebugLevel;
37
38 void EgVtkObject::computeNormals
39 (
40 QVector<vec3_t> &cell_normals,
41 QVector<vec3_t> &node_normals,
42 QVector<vtkIdType> &cells,
43 QVector<vtkIdType> &nodes,
44 vtkUnstructuredGrid *grid
45 )
46 {
47 using namespace GeometryTools;
48
49 cell_normals.resize(cells.size());
50 node_normals.fill(vec3_t(0,0,0), nodes.size());
51 QVector<int> g2s;
52 createNodeMapping(nodes, g2s, grid);
53 for (int i_cell = 0; i_cell < cells.count(); ++i_cell) {
54 vtkIdType id_cell = cells[i_cell];
55 EG_GET_CELL(id_cell, grid);
56 cell_normals[i_cell] = cellNormal(grid, id_cell);
57 cell_normals[i_cell].normalise();
58 for (int i_pts = 0; i_pts < num_pts; ++i_pts) {
59 if (g2s[pts[i_pts]] != -1) {
60 node_normals[g2s[pts[i_pts]]] += cell_normals[i_cell];
61 }
62 }
63 }
64 for (int i_node = 0; i_node < nodes.count(); ++i_node) {
65 node_normals[i_node].normalise();
66 //cout << node_normals[i_node] << endl;
67 }
68 }
69
70 void EgVtkObject::createNodeMapping
71 (
72 QVector<vtkIdType> &nodes,
73 QVector<int> &_nodes,
74 vtkUnstructuredGrid *grid
75 )
76 {
77 _nodes.fill(-1,grid->GetNumberOfPoints());
78 for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
79 _nodes[nodes[i_nodes]] = i_nodes;
80 }
81 }
82
83 void EgVtkObject::createCellMapping
84 (
85 QVector<vtkIdType> &cells,
86 QVector<int> &_cells,
87 vtkUnstructuredGrid *grid
88 )
89 {
90 _cells.fill(-1,grid->GetNumberOfCells());
91 for (int i_cells = 0; i_cells < cells.size(); ++i_cells) {
92 _cells[cells[i_cells]] = i_cells;
93 }
94 }
95
96 void EgVtkObject::createNodeToBcMapping
97 (
98 QVector<QSet<int> > &bcs,
99 vtkUnstructuredGrid *grid
100 )
101 {
102 EG_BUG;
103 bcs.fill(QSet<int>(), grid->GetNumberOfPoints());
104 grid->BuildLinks();
105 auto abstract_links = grid->GetCellLinks();
106 auto cell_links = dynamic_cast<vtkCellLinks*>(abstract_links);
107 EG_VTKDCC(vtkIntArray, cell_code, grid, "cell_code");
108 for (vtkIdType nodeId = 0; nodeId < grid->GetNumberOfPoints(); ++nodeId) {
109 int Ncells = cell_links->GetNcells(nodeId);
110 for (int i = 0; i < Ncells; ++i) {
111 vtkIdType id_cell = cell_links->GetCells(nodeId)[i];
112 vtkIdType ct = grid->GetCellType(id_cell);
113 if ((ct == VTK_TRIANGLE) || (ct = VTK_QUAD)) {
114 if (cell_code->GetValue(id_cell) > 0) {
115 bcs[nodeId].insert(cell_code->GetValue(id_cell));
116 }
117 }
118 }
119 }
120 }
121
122 void EgVtkObject::createNodeToCell
123 (
124 QVector<vtkIdType> &cells,
125 QVector<vtkIdType> &nodes,
126 QVector<int> &_nodes,
127 QVector<QSet<int> > &n2c,
128 vtkUnstructuredGrid *grid
129 )
130 {
131 n2c.fill(QSet<int>(), nodes.size());
132 for (vtkIdType i_cells = 0; i_cells < cells.size(); ++i_cells) {
133 QList<vtkIdType> pts;
134 getPointsOfCell(grid, cells[i_cells], pts);
135 foreach (vtkIdType id_node, pts) {
136 n2c[_nodes[id_node]].insert(i_cells);
137 }
138 }
139 }
140
141 void EgVtkObject::createNodeToCell
142 (
143 QVector<vtkIdType> &cells,
144 QVector<vtkIdType> &nodes,
145 QVector<int> &_nodes,
146 QVector<QVector<int> > &n2c,
147 vtkUnstructuredGrid *grid
148 )
149 {
150 n2c.fill(QVector<int>(), nodes.size());
151 QVector<int> count(nodes.size(),0);
152 for (vtkIdType i_cells = 0; i_cells < cells.size(); ++i_cells) {
153 QList<vtkIdType> pts;
154 getPointsOfCell(grid, cells[i_cells], pts);
155 foreach (vtkIdType id_node, pts) {
156 ++count[_nodes[id_node]];
157 }
158 }
159 for (int i = 0; i < nodes.size(); ++i) {
160 n2c[i].resize(count[i]);
161 count[i] = 0;
162 }
163 for (vtkIdType i_cells = 0; i_cells < cells.size(); ++i_cells) {
164 QList<vtkIdType> pts;
165 getPointsOfCell(grid, cells[i_cells], pts);
166 foreach (vtkIdType id_node, pts) {
167 int i_nodes = _nodes[id_node];
168 n2c[i_nodes][count[i_nodes]] = i_cells;
169 ++count[i_nodes];
170 }
171 }
172 }
173
174 void EgVtkObject::addToN2N(QVector<QSet<int> > &n2n, int n1, int n2)
175 {
176 n2n[n1].insert(n2);
177 n2n[n2].insert(n1);
178 }
179
180 void EgVtkObject::createNodeToNode(QVector<vtkIdType> &cells, QVector<vtkIdType> &nodes, QVector<int> &_nodes, QVector<QSet<int> > &n2n, vtkUnstructuredGrid *grid)
181 {
182 n2n.fill(QSet<int>(), nodes.size());
183 foreach (vtkIdType id_cell, cells) {
184 if (isSurface(id_cell, grid)) {
185 EG_GET_CELL(id_cell, grid);
186 QVector<int> n(num_pts);
187 for (int i = 0; i < num_pts; ++i) {
188 n[i] = _nodes[pts[i]];
189 }
190 for (int i = 0; i < num_pts - 1; ++i) {
191 addToN2N(n2n, n[i], n[i+1]);
192 }
193 addToN2N(n2n, n.last(), n.first());
194 } else {
195 vtkIdType num_faces = 0;
196 vtkIdType type_cell = grid->GetCellType(id_cell);
197
198 if (type_cell == VTK_TETRA) {
199 num_faces = 4;
200 } else if (type_cell == VTK_WEDGE) {
201 num_faces = 5;
202 } else if (type_cell == VTK_PYRAMID) {
203 num_faces = 5;
204 } else if (type_cell == VTK_HEXAHEDRON) {
205 num_faces = 6;
206 } else if (type_cell == VTK_POLYHEDRON) {
207 vtkSmartPointer<vtkIdList> stream = vtkSmartPointer<vtkIdList>::New();
208 grid->GetFaceStream(id_cell, stream);
209 num_faces = stream->GetId(0);
210 }
211
212 for (int i = 0; i < num_faces; ++i) {
213 QVector<vtkIdType> face;
214 getFaceOfCell(grid, id_cell, i, face);
215 QVector<int> n(face.size());
216 for (int j = 0; j < face.size(); ++j) {
217 n[j] = _nodes[face[j]];
218 }
219 for (int j = 0; j < face.size() - 1; ++j) {
220 addToN2N(n2n, n[j], n[j+1]);
221 }
222 addToN2N(n2n, n.last(), n.first());
223 }
224 }
225
226
227
228 /*
229 vtkIdType *pts;
230 vtkIdType Npts;
231 grid->GetCellPoints(id_cell, Npts, pts);
232 vector<int> n(Npts);
233 for (int i = 0; i < Npts; ++i) {
234 n[i] = _nodes[pts[i]];
235 }
236 vtkIdType type_cell = grid->GetCellType(id_cell);
237 if (type_cell == VTK_TRIANGLE) {
238 addToN2N(n2n, n[0], n[1]);
239 addToN2N(n2n, n[1], n[2]);
240 addToN2N(n2n, n[2], n[0]);
241 } else if (type_cell == VTK_QUAD) {
242 addToN2N(n2n, n[0], n[1]);
243 addToN2N(n2n, n[1], n[2]);
244 addToN2N(n2n, n[2], n[3]);
245 addToN2N(n2n, n[3], n[0]);
246 } else if (type_cell == VTK_POLYGON) {
247 EG_BUG;
248 } else if (type_cell == VTK_TETRA) {
249 addToN2N(n2n, n[0], n[1]);
250 addToN2N(n2n, n[0], n[2]);
251 addToN2N(n2n, n[0], n[3]);
252 addToN2N(n2n, n[1], n[2]);
253 addToN2N(n2n, n[1], n[3]);
254 addToN2N(n2n, n[2], n[3]);
255 } else if (type_cell == VTK_PYRAMID) {
256 addToN2N(n2n, n[0], n[1]);
257 addToN2N(n2n, n[0], n[3]);
258 addToN2N(n2n, n[0], n[4]);
259 addToN2N(n2n, n[1], n[2]);
260 addToN2N(n2n, n[1], n[4]);
261 addToN2N(n2n, n[2], n[3]);
262 addToN2N(n2n, n[2], n[4]);
263 addToN2N(n2n, n[3], n[4]);
264 } else if (type_cell == VTK_WEDGE) {
265 addToN2N(n2n, n[0], n[1]);
266 addToN2N(n2n, n[0], n[2]);
267 addToN2N(n2n, n[0], n[3]);
268 addToN2N(n2n, n[1], n[2]);
269 addToN2N(n2n, n[1], n[4]);
270 addToN2N(n2n, n[2], n[5]);
271 addToN2N(n2n, n[3], n[4]);
272 addToN2N(n2n, n[3], n[5]);
273 addToN2N(n2n, n[4], n[5]);
274 } else if (type_cell == VTK_HEXAHEDRON) {
275 addToN2N(n2n, n[0], n[1]);
276 addToN2N(n2n, n[0], n[3]);
277 addToN2N(n2n, n[0], n[4]);
278 addToN2N(n2n, n[1], n[2]);
279 addToN2N(n2n, n[1], n[5]);
280 addToN2N(n2n, n[2], n[3]);
281 addToN2N(n2n, n[2], n[6]);
282 addToN2N(n2n, n[3], n[7]);
283 addToN2N(n2n, n[4], n[5]);
284 addToN2N(n2n, n[4], n[7]);
285 addToN2N(n2n, n[5], n[6]);
286 addToN2N(n2n, n[6], n[7]);
287 } else if (type_cell == VTK_POLYHEDRON) {
288 EG_BUG;
289 }
290 */
291 }
292 }
293
294 void EgVtkObject::createNodeToNode
295 (
296 QVector<vtkIdType> &cells,
297 QVector<vtkIdType> &nodes,
298 QVector<int> &_nodes,
299 QVector<QVector<int> > &n2n,
300 vtkUnstructuredGrid *grid
301 )
302 {
303 QVector<QSet<int> > n2n_set;
304 createNodeToNode(cells, nodes, _nodes, n2n_set, grid);
305 n2n.resize(n2n_set.size());
306 for (int i = 0; i < n2n.size(); ++i) {
307 n2n[i].resize(n2n_set[i].size());
308 qCopy(n2n_set[i].begin(), n2n_set[i].end(), n2n[i].begin());
309 }
310 }
311
312 void EgVtkObject::getAllCells
313 (
314 QVector<vtkIdType> &cells,
315 vtkUnstructuredGrid *grid
316 )
317 {
318 int N = 0;
319 cells.resize(grid->GetNumberOfCells());
320 for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
321 cells[N] = id_cell;
322 ++N;
323 }
324 }
325
326 void EgVtkObject::getAllCellsOfType
327 (
328 vtkIdType type,
329 QVector<vtkIdType> &cells,
330 vtkUnstructuredGrid *grid
331 )
332 {
333 int N = 0;
334 for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
335 if (grid->GetCellType(id_cell) == type) {
336 ++N;
337 }
338 }
339 cells.resize(N);
340 N = 0;
341 for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
342 if (grid->GetCellType(id_cell) == type) {
343 cells[N] = id_cell;
344 ++N;
345 }
346 }
347 }
348
349
350 void EgVtkObject::getAllVolumeCells
351 (
352 QVector<vtkIdType> &cells,
353 vtkUnstructuredGrid *grid
354 )
355 {
356 int N = 0;
357 for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
358 if (isVolume(id_cell, grid)) {
359 ++N;
360 }
361 }
362 cells.resize(N);
363 N = 0;
364 for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
365 if (isVolume(id_cell, grid)) {
366 cells[N] = id_cell;
367 ++N;
368 }
369 }
370 }
371
372 void EgVtkObject::getAllSurfaceCells
373 (
374 QVector<vtkIdType> &cells,
375 vtkUnstructuredGrid *grid
376 )
377 {
378 int N = 0;
379 for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
380 if (isSurface(id_cell, grid)) {
381 ++N;
382 }
383 }
384 cells.resize(N);
385 N = 0;
386 for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
387 if (isSurface(id_cell, grid)) {
388 cells[N] = id_cell;
389 ++N;
390 }
391 }
392 }
393
394 void EgVtkObject::addToC2C(vtkIdType id_cell, QVector<int> &_cells, QVector<QVector<int> > &c2c, int j, vtkIdList *nds, vtkIdList *cls, vtkUnstructuredGrid *grid)
395 {
396 c2c[_cells[id_cell]][j] = -1;
397 grid->GetCellNeighbors(id_cell, nds, cls);
398 if (isSurface(id_cell, grid)) {
399 for (int i = 0; i < cls->GetNumberOfIds(); ++i) {
400 if (cls->GetId(i) != id_cell) {
401 if (_cells[cls->GetId(i)] != -1) {
402 if (isSurface(cls->GetId(i), grid)) {
403 c2c[_cells[id_cell]][j] = _cells[cls->GetId(i)];
404 }
405 }
406 }
407 }
408 } else {
409 for (int i = 0; i < cls->GetNumberOfIds(); ++i) {
410 if (cls->GetId(i) != id_cell) {
411 if (_cells[cls->GetId(i)] != -1) {
412 if (isVolume(cls->GetId(i), grid) || c2c[_cells[id_cell]][j] == -1) {
413 c2c[_cells[id_cell]][j] = _cells[cls->GetId(i)];
414 }
415 }
416 }
417 }
418 }
419 }
420
421
422 void EgVtkObject::createCellToCell(QVector<vtkIdType> &cells, QVector<QVector<int> > &c2c, vtkUnstructuredGrid *grid)
423 {
424 // GetCellNeighbors(vtkIdType id_cell, vtkIdList *ptIds, vtkIdList *id_cells)
425 grid->BuildLinks();
426 QVector<int> _cells;
427 createCellMapping(cells, _cells, grid);
428 c2c.fill(QVector<int>(), cells.size());
429 EG_VTKSP(vtkIdList, nds);
430 EG_VTKSP(vtkIdList, cls);
431 for (int i = 0; i < cells.size(); ++i) {
432 vtkIdType id_cell = cells[i];
433 //vtkIdType *pts;
434 //vtkIdType Npts;
435 //grid->GetCellPoints(id_cell, Npts, pts);
436 QList<vtkIdType> pts;
437 getPointsOfCell(grid, id_cell, pts);
438 if (grid->GetCellType(id_cell) == VTK_TRIANGLE) {
439 c2c[i].resize(3);
440 nds->Reset();
441 nds->InsertNextId(pts[0]);
442 nds->InsertNextId(pts[1]);
443 addToC2C(id_cell, _cells, c2c, 0, nds, cls, grid);
444 nds->Reset();
445 nds->InsertNextId(pts[1]);
446 nds->InsertNextId(pts[2]);
447 addToC2C(id_cell, _cells, c2c, 1, nds, cls, grid);
448 nds->Reset();
449 nds->InsertNextId(pts[2]);
450 nds->InsertNextId(pts[0]);
451 addToC2C(id_cell, _cells, c2c, 2, nds, cls, grid);
452 } else if (grid->GetCellType(id_cell) == VTK_QUAD) {
453 c2c[i].resize(4);
454 nds->Reset();
455 nds->InsertNextId(pts[0]);
456 nds->InsertNextId(pts[1]);
457 addToC2C(id_cell, _cells, c2c, 0, nds, cls, grid);
458 nds->Reset();
459 nds->InsertNextId(pts[1]);
460 nds->InsertNextId(pts[2]);
461 addToC2C(id_cell, _cells, c2c, 1, nds, cls, grid);
462 nds->Reset();
463 nds->InsertNextId(pts[2]);
464 nds->InsertNextId(pts[3]);
465 addToC2C(id_cell, _cells, c2c, 2, nds, cls, grid);
466 nds->Reset();
467 nds->InsertNextId(pts[3]);
468 nds->InsertNextId(pts[0]);
469 addToC2C(id_cell, _cells, c2c, 3, nds, cls, grid);
470 } else if (grid->GetCellType(id_cell) == VTK_POLYGON) {
471 c2c[i].resize(pts.size());
472 pts.append(pts[0]);
473 for (int i = 0; i < pts.size() - 1; ++i) {
474 nds->Reset();
475 nds->InsertNextId(pts[i]);
476 nds->InsertNextId(pts[i+1]);
477 addToC2C(id_cell, _cells, c2c, i, nds, cls, grid);
478 }
479 } else if (grid->GetCellType(id_cell) == VTK_TETRA) {
480 c2c[i].resize(4);
481 nds->Reset();
482 nds->InsertNextId(pts[0]);
483 nds->InsertNextId(pts[1]);
484 nds->InsertNextId(pts[2]);
485 addToC2C(id_cell, _cells, c2c, 0, nds, cls, grid);
486 nds->Reset();
487 nds->InsertNextId(pts[0]);
488 nds->InsertNextId(pts[1]);
489 nds->InsertNextId(pts[3]);
490 addToC2C(id_cell, _cells, c2c, 1, nds, cls, grid);
491 nds->Reset();
492 nds->InsertNextId(pts[0]);
493 nds->InsertNextId(pts[3]);
494 nds->InsertNextId(pts[2]);
495 addToC2C(id_cell, _cells, c2c, 2, nds, cls, grid);
496 nds->Reset();
497 nds->InsertNextId(pts[1]);
498 nds->InsertNextId(pts[2]);
499 nds->InsertNextId(pts[3]);
500 addToC2C(id_cell, _cells, c2c, 3, nds, cls, grid);
501 } else if (grid->GetCellType(id_cell) == VTK_PYRAMID) {
502 c2c[i].resize(5);
503 nds->Reset();
504 nds->InsertNextId(pts[0]);
505 nds->InsertNextId(pts[1]);
506 nds->InsertNextId(pts[2]);
507 nds->InsertNextId(pts[3]);
508 addToC2C(id_cell, _cells, c2c, 0, nds, cls, grid);
509 nds->Reset();
510 nds->InsertNextId(pts[0]);
511 nds->InsertNextId(pts[1]);
512 nds->InsertNextId(pts[4]);
513 addToC2C(id_cell, _cells, c2c, 1, nds, cls, grid);
514 nds->Reset();
515 nds->InsertNextId(pts[1]);
516 nds->InsertNextId(pts[2]);
517 nds->InsertNextId(pts[4]);
518 addToC2C(id_cell, _cells, c2c, 2, nds, cls, grid);
519 nds->Reset();
520 nds->InsertNextId(pts[2]);
521 nds->InsertNextId(pts[3]);
522 nds->InsertNextId(pts[4]);
523 addToC2C(id_cell, _cells, c2c, 3, nds, cls, grid);
524 nds->Reset();
525 nds->InsertNextId(pts[3]);
526 nds->InsertNextId(pts[0]);
527 nds->InsertNextId(pts[4]);
528 addToC2C(id_cell, _cells, c2c, 4, nds, cls, grid);
529 } else if (grid->GetCellType(id_cell) == VTK_WEDGE) {
530 c2c[i].resize(5);
531 nds->Reset();
532 nds->InsertNextId(pts[0]);
533 nds->InsertNextId(pts[1]);
534 nds->InsertNextId(pts[2]);
535 addToC2C(id_cell, _cells, c2c, 0, nds, cls, grid);
536 nds->Reset();
537 nds->InsertNextId(pts[3]);
538 nds->InsertNextId(pts[4]);
539 nds->InsertNextId(pts[5]);
540 addToC2C(id_cell, _cells, c2c, 1, nds, cls, grid);
541 nds->Reset();
542 nds->InsertNextId(pts[0]);
543 nds->InsertNextId(pts[1]);
544 nds->InsertNextId(pts[4]);
545 nds->InsertNextId(pts[3]);
546 addToC2C(id_cell, _cells, c2c, 2, nds, cls, grid);
547 nds->Reset();
548 nds->InsertNextId(pts[1]);
549 nds->InsertNextId(pts[4]);
550 nds->InsertNextId(pts[5]);
551 nds->InsertNextId(pts[2]);
552 addToC2C(id_cell, _cells, c2c, 3, nds, cls, grid);
553 nds->Reset();
554 nds->InsertNextId(pts[0]);
555 nds->InsertNextId(pts[2]);
556 nds->InsertNextId(pts[5]);
557 nds->InsertNextId(pts[3]);
558 addToC2C(id_cell, _cells, c2c, 4, nds, cls, grid);
559 } else if (grid->GetCellType(id_cell) == VTK_HEXAHEDRON) {
560 c2c[i].resize(6);
561 nds->Reset();
562 nds->InsertNextId(pts[0]);
563 nds->InsertNextId(pts[3]);
564 nds->InsertNextId(pts[2]);
565 nds->InsertNextId(pts[1]);
566 addToC2C(id_cell, _cells, c2c, 0, nds, cls, grid);
567 nds->Reset();
568 nds->InsertNextId(pts[4]);
569 nds->InsertNextId(pts[5]);
570 nds->InsertNextId(pts[6]);
571 nds->InsertNextId(pts[7]);
572 addToC2C(id_cell, _cells, c2c, 1, nds, cls, grid);
573 nds->Reset();
574 nds->InsertNextId(pts[0]);
575 nds->InsertNextId(pts[1]);
576 nds->InsertNextId(pts[5]);
577 nds->InsertNextId(pts[4]);
578 addToC2C(id_cell, _cells, c2c, 2, nds, cls, grid);
579 nds->Reset();
580 nds->InsertNextId(pts[3]);
581 nds->InsertNextId(pts[7]);
582 nds->InsertNextId(pts[6]);
583 nds->InsertNextId(pts[2]);
584 addToC2C(id_cell, _cells, c2c, 3, nds, cls, grid);
585 nds->Reset();
586 nds->InsertNextId(pts[0]);
587 nds->InsertNextId(pts[4]);
588 nds->InsertNextId(pts[7]);
589 nds->InsertNextId(pts[3]);
590 addToC2C(id_cell, _cells, c2c, 4, nds, cls, grid);
591 nds->Reset();
592 nds->InsertNextId(pts[1]);
593 nds->InsertNextId(pts[2]);
594 nds->InsertNextId(pts[6]);
595 nds->InsertNextId(pts[5]);
596 addToC2C(id_cell, _cells, c2c, 5, nds, cls, grid);
597 } else if (grid->GetCellType(id_cell) == VTK_POLYHEDRON) {
598 EG_VTKSP(vtkIdList, stream);
599 grid->GetFaceStream(id_cell, stream);
600 int num_faces = stream->GetId(0);
601 c2c[i].resize(num_faces);
602 int id = 1;
603 for (int i_face = 0; i_face < num_faces; ++i_face) {
604 nds->Reset();
605 int num_pts = stream->GetId(id++);
606 for (int i_pt = 0; i_pt < num_pts; ++i_pt) {
607 nds->InsertNextId(stream->GetId(id++));
608 }
609 addToC2C(id_cell, _cells, c2c, i_face, nds, cls, grid);
610 }
611 }
612 }
613 }
614
615 bool EgVtkObject::isVolume(vtkIdType id_cell, vtkUnstructuredGrid *grid)
616 {
617 bool isVol = false;
618 if (grid->GetCellType(id_cell) == VTK_TETRA) isVol = true;
619 else if (grid->GetCellType(id_cell) == VTK_PYRAMID) isVol = true;
620 else if (grid->GetCellType(id_cell) == VTK_WEDGE) isVol = true;
621 else if (grid->GetCellType(id_cell) == VTK_HEXAHEDRON) isVol = true;
622 else if (grid->GetCellType(id_cell) == VTK_POLYHEDRON) isVol = true;
623 return isVol;
624 }
625
626 bool EgVtkObject::isSurface(vtkIdType id_cell, vtkUnstructuredGrid *grid)
627 {
628 bool isSurf = false;
629 if (grid->GetCellType(id_cell) == VTK_TRIANGLE) isSurf = true;
630 else if (grid->GetCellType(id_cell) == VTK_QUAD) isSurf = true;
631 else if (grid->GetCellType(id_cell) == VTK_POLYGON) isSurf = true;
632 return isSurf;
633 }
634
635 void EgVtkObject::updateCellIndex(vtkUnstructuredGrid *grid)
636 {
637 if (!grid->GetCellData()->GetArray("cell_index")) {
638 EG_VTKSP(vtkLongArray_t, cell_index);
639 cell_index->SetName("cell_index");
640 cell_index->SetNumberOfValues(grid->GetNumberOfCells());
641 grid->GetCellData()->AddArray(cell_index);
642 }
643 EG_VTKDCC(vtkLongArray_t, cell_index, grid, "cell_index");
644 for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
645 cell_index->SetValue(id_cell, id_cell);
646 }
647 }
648
649 void EgVtkObject::updateNodeIndex(vtkUnstructuredGrid *grid)
650 {
651 if (!grid->GetPointData()->GetArray("node_index")) {
652 EG_VTKSP(vtkLongArray_t, node_index);
653 node_index->SetName("node_index");
654 node_index->SetNumberOfValues(grid->GetNumberOfPoints());
655 grid->GetPointData()->AddArray(node_index);
656 }
657 EG_VTKDCN(vtkLongArray_t, node_index, grid, "node_index");
658 for (vtkIdType pointId = 0; pointId < grid->GetNumberOfPoints(); ++pointId) {
659 node_index->SetValue(pointId, pointId);
660 }
661 }
662
663 void EgVtkObject::addToPolyData
664 (
665 QVector<vtkIdType> &cells,
666 vtkPolyData *pdata,
667 vtkUnstructuredGrid *grid
668 )
669 {
670 updateCellIndex(grid);
671 updateNodeIndex(grid);
672 QVector<vtkIdType> nodes;
673 QVector<int> _nodes;
674 getNodesFromCells(cells, nodes, grid);
675 createNodeMapping(nodes, _nodes, grid);
676 EG_VTKSP(vtkDoubleArray, pcoords);
677 pcoords->SetNumberOfComponents(3);
678 pcoords->SetNumberOfTuples(nodes.size());
679 EG_VTKSP(vtkPoints, points);
680 points->SetData(pcoords);
681 pdata->SetPoints(points);
682 pdata->Allocate(cells.size());
683 if (!pdata->GetCellData()->GetArray("cell_index")) {
684 EG_VTKSP(vtkLongArray_t, cell_index);
685 cell_index->SetName("cell_index");
686 //cell_index->SetNumberOfValues(cells.size());
687 pdata->GetCellData()->AddArray(cell_index);
688 }
689 if (!pdata->GetPointData()->GetArray("node_index")) {
690 EG_VTKSP(vtkLongArray_t, node_index);
691 node_index->SetName("node_index");
692 //node_index->SetNumberOfValues(nodes.size());
693 pdata->GetPointData()->AddArray(node_index);
694 }
695 EG_VTKDCC(vtkLongArray_t, pd_cell_index, pdata, "cell_index");
696 EG_VTKDCN(vtkLongArray_t, pd_node_index, pdata, "node_index");
697 pd_cell_index->SetNumberOfValues(cells.size());
698 pd_node_index->SetNumberOfValues(nodes.size());
699 for (int i_cell = 0; i_cell < cells.size(); ++i_cell) {
700 vtkIdType id_cell = cells[i_cell];
701 vtkIdType cell_type = grid->GetCellType(id_cell);
702 if ((cell_type != VTK_TRIANGLE) && (cell_type != VTK_QUAD)) {
703 EG_ERR_RETURN("unsupported cell type for this operation");
704 }
705 EG_GET_CELL(id_cell, grid);
706 vtkSmartPointer<vtkIdList> new_pts = vtkSmartPointer<vtkIdList>::New();
707 new_pts->SetNumberOfIds(num_pts);
708 for (int i = 0; i < num_pts; ++i) {
709 new_pts->SetId(i, _nodes[pts[i]]);
710 }
711 vtkIdType new_cell_id = pdata->InsertNextCell(cell_type, new_pts);
712 pd_cell_index->SetValue(new_cell_id, id_cell);
713 }
714 for (int i_node = 0; i_node < nodes.size(); ++i_node) {
715 vec3_t x;
716 grid->GetPoints()->GetPoint(nodes[i_node], x.data());
717 pdata->GetPoints()->SetPoint(i_node, x.data());
718 pd_node_index->SetValue(i_node, nodes[i_node]);
719 }
720 }
721
722 #define EGVTKOBJECT_COPYCELLDATA(FIELD,TYPE) \
723 { \
724 if (old_grid->GetCellData()->GetArray(FIELD)) { \
725 EG_VTKDCC(TYPE, var1, old_grid, FIELD); \
726 EG_VTKDCC(TYPE, var2, new_grid, FIELD); \
727 var2->SetValue(newId, var1->GetValue(oldId)); \
728 } \
729 }
730
731 void EgVtkObject::copyCellData
732 (
733 vtkUnstructuredGrid *old_grid,
734 vtkIdType oldId,
735 vtkUnstructuredGrid *new_grid,
736 vtkIdType newId
737 )
738 {
739 EGVTKOBJECT_COPYCELLDATA("vtk_type", vtkIntArray);
740 EGVTKOBJECT_COPYCELLDATA("cell_code", vtkIntArray);
741 EGVTKOBJECT_COPYCELLDATA("cell_orgdir", vtkIntArray);
742 EGVTKOBJECT_COPYCELLDATA("cell_curdir", vtkIntArray);
743 EGVTKOBJECT_COPYCELLDATA("cell_voldir", vtkIntArray);
744 EGVTKOBJECT_COPYCELLDATA("cell_index", vtkLongArray_t);
745 EGVTKOBJECT_COPYCELLDATA("cell_subres", vtkDoubleArray);
746 }
747
748 #define EGVTKOBJECT_COPYNODEDATA(FIELD,TYPE) \
749 { \
750 if (old_grid->GetPointData()->GetArray(FIELD)) { \
751 EG_VTKDCN(TYPE, var1, old_grid, FIELD); \
752 EG_VTKDCN(TYPE, var2, new_grid, FIELD); \
753 var2->SetValue(newId, var1->GetValue(oldId)); \
754 } \
755 }
756
757 void EgVtkObject::copyNodeData
758 (
759 vtkUnstructuredGrid *old_grid,
760 vtkIdType oldId,
761 vtkUnstructuredGrid *new_grid,
762 vtkIdType newId
763 )
764 {
765 EGVTKOBJECT_COPYNODEDATA("node_status", vtkIntArray);
766 EGVTKOBJECT_COPYNODEDATA("node_layer", vtkIntArray);
767 EGVTKOBJECT_COPYNODEDATA("node_index", vtkLongArray_t);
768 EGVTKOBJECT_COPYNODEDATA("node_specified_density", vtkIntArray);
769 EGVTKOBJECT_COPYNODEDATA("node_meshdensity_desired", vtkDoubleArray);
770 //EGVTKOBJECT_COPYNODEDATA("node_meshdensity_current", vtkDoubleArray);
771 EGVTKOBJECT_COPYNODEDATA("node_type", vtkCharArray_t);
772 EGVTKOBJECT_COPYNODEDATA("node_pindex", vtkLongArray_t);
773 }
774
775 #define EGVTKOBJECT_CREATECELLFIELD(FIELD,TYPE,OW) \
776 if (!grid->GetCellData()->GetArray(FIELD)) { \
777 EG_VTKSP(TYPE, var); \
778 var->SetName(FIELD); \
779 var->SetNumberOfValues(Ncells); \
780 grid->GetCellData()->AddArray(var); \
781 for (int i = 0; i < grid->GetNumberOfCells(); ++i) { \
782 var->SetValue(i,0); \
783 } \
784 } else if (OW) { \
785 EG_VTKDCC(TYPE, var, grid, FIELD); \
786 var->SetNumberOfValues(Ncells); \
787 for (int i = 0; i < grid->GetNumberOfCells(); ++i) { \
788 var->SetValue(i,0); \
789 } \
790 }
791
792 #define EGVTKOBJECT_CREATENODEFIELD(FIELD,TYPE,OW) \
793 if (!grid->GetPointData()->GetArray(FIELD)) { \
794 EG_VTKSP(TYPE, var); \
795 var->SetName(FIELD); \
796 var->SetNumberOfValues(Nnodes); \
797 grid->GetPointData()->AddArray(var); \
798 for (int i = 0; i < grid->GetNumberOfPoints(); ++i) { \
799 var->SetValue(i,0); \
800 } \
801 } else if (OW) { \
802 EG_VTKDCN(TYPE, var, grid, FIELD); \
803 var->SetNumberOfValues(Nnodes); \
804 for (int i = 0; i < grid->GetNumberOfPoints(); ++i) { \
805 var->SetValue(i,0); \
806 } \
807 }
808
809 void EgVtkObject::createBasicFields(vtkUnstructuredGrid *grid, vtkIdType num_cells, vtkIdType num_nodes, bool overwrite)
810 {
811 if (num_cells == -1) {
812 num_cells = grid->GetNumberOfCells();
813 }
814 if (num_nodes == -1) {
815 num_nodes = grid->GetNumberOfPoints();
816 }
817 createBasicNodeFields(grid, num_nodes, overwrite);
818 createBasicCellFields(grid, num_cells, overwrite);
819 }
820
821 void EgVtkObject::createBasicCellFields(vtkUnstructuredGrid *grid, vtkIdType Ncells, bool overwrite)
822 {
823 EGVTKOBJECT_CREATECELLFIELD("vtk_type" , vtkIntArray, overwrite);
824 EGVTKOBJECT_CREATECELLFIELD("cell_code", vtkIntArray, overwrite);
825 EGVTKOBJECT_CREATECELLFIELD("cell_index", vtkLongArray_t, overwrite);
826 EGVTKOBJECT_CREATECELLFIELD("cell_orgdir", vtkIntArray, overwrite); // original orientation
827 EGVTKOBJECT_CREATECELLFIELD("cell_curdir", vtkIntArray, overwrite); // current orientation
828 EGVTKOBJECT_CREATECELLFIELD("cell_voldir", vtkIntArray, overwrite); // volume orientation -- only valid for a single (i.e. the current) volume
829 EGVTKOBJECT_CREATECELLFIELD("cell_mesh_quality", vtkDoubleArray, overwrite); // generic field to store different quality measures
830 EGVTKOBJECT_CREATECELLFIELD("cell_subres", vtkDoubleArray, overwrite); // mesh size for sub-cell resolution (DrNUM)
831 }
832
833 void EgVtkObject::createBasicNodeFields(vtkUnstructuredGrid *grid, vtkIdType Nnodes, bool overwrite)
834 {
835 EGVTKOBJECT_CREATENODEFIELD("node_status", vtkIntArray, overwrite);
836 EGVTKOBJECT_CREATENODEFIELD("node_layer", vtkIntArray, overwrite);
837 EGVTKOBJECT_CREATENODEFIELD("node_index", vtkLongArray_t, overwrite);
838 EGVTKOBJECT_CREATENODEFIELD("node_specified_density", vtkIntArray, overwrite); //density index from table
839 EGVTKOBJECT_CREATENODEFIELD("node_meshdensity_desired", vtkDoubleArray, overwrite); //what we want
840 EGVTKOBJECT_CREATENODEFIELD("node_type", vtkCharArray_t, overwrite); //node type
841 EGVTKOBJECT_CREATENODEFIELD("node_type_counter", vtkIntArray, overwrite); // counter field to delay node type demotion
842 EGVTKOBJECT_CREATENODEFIELD("node_pindex", vtkLongArray_t, overwrite);
843 EGVTKOBJECT_CREATENODEFIELD("node_mesh_quality", vtkDoubleArray, overwrite); // generic field to store different quality measures
844 }
845
846 void EgVtkObject::allocateGrid(vtkUnstructuredGrid *grid, vtkIdType Ncells, vtkIdType Nnodes, bool create_fields)
847 {
848 EG_VTKSP(vtkPoints, points);
849 points->SetNumberOfPoints(Nnodes);
850 grid->SetPoints(points);
851 grid->Allocate(0, 0);
852 grid->Squeeze();
853 grid->Reset();
854 grid->Allocate(Ncells,max(vtkIdType(1),Ncells/10));
855 if (create_fields) {
856 createBasicFields(grid, Ncells, Nnodes, true);
857 }
858 }
859
860 vec3_t EgVtkObject::cellCentre(vtkUnstructuredGrid *grid, vtkIdType id_cell)
861 {
862 vec3_t x,xc(0,0,0);
863 EG_GET_CELL(id_cell, grid);
864 double f = 1.0/num_pts;
865 for (int i_pts = 0; i_pts < num_pts; ++i_pts) {
866 grid->GetPoint(pts[i_pts], x.data());
867 xc += f*x;
868 }
869 return xc;
870 }
871
872 double EgVtkObject::faceAngle(vtkUnstructuredGrid *grid, vtkIdType id_face1, vtkIdType id_face2)
873 {
874 QList<vtkIdType> edge_nodes;
875 sharedNodesOfCells(grid, id_face1, id_face2, edge_nodes);
876 if (edge_nodes.size() != 2) {
877 EG_BUG;
878 }
879 vec3_t x1, x2;
880 grid->GetPoint(edge_nodes[0], x1.data());
881 grid->GetPoint(edge_nodes[1], x2.data());
882 vec3_t xf1 = cellCentre(grid, id_face1);
883 vec3_t xf2 = cellCentre(grid, id_face2);
884 vec3_t xe = 0.5*(x1 + x2);
885 vec3_t xfe = 0.5*(xf1 + xf2);
886 vec3_t n1 = cellNormal(grid, id_face1);
887 vec3_t n2 = cellNormal(grid, id_face2);
888 vec3_t ve = xe - xfe;
889 double alpha = angle(n1, n2);
890 if ((n1 + n2)*ve < 0) {
891 return -alpha;
892 }
893 return alpha;
894 }
895
896 void EgVtkObject::getRestCells(vtkUnstructuredGrid *grid,
897 const QVector<vtkIdType> &cells,
898 QVector<vtkIdType> &rest_cells)
899 {
900 QVector<bool> is_in_cells(grid->GetNumberOfCells(), false);
901 foreach (vtkIdType id_cell, cells) {
902 is_in_cells[id_cell] = true;
903 }
904 rest_cells.resize(grid->GetNumberOfCells() - cells.size());
905 int i_rest_cells = 0;
906 for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
907 if (!is_in_cells[id_cell]) {
908 rest_cells[i_rest_cells] = id_cell;
909 ++i_rest_cells;
910 }
911 }
912 }
913
914 void EgVtkObject::makeCopy(vtkUnstructuredGrid *src, vtkUnstructuredGrid *dst, bool copy_data)
915 {
916 allocateGrid(dst, src->GetNumberOfCells(), src->GetNumberOfPoints(), copy_data);
917 for (vtkIdType id_node = 0; id_node < src->GetNumberOfPoints(); ++id_node) {
918 vec3_t x;
919 src->GetPoints()->GetPoint(id_node, x.data());
920 dst->GetPoints()->SetPoint(id_node, x.data());
921 if (copy_data) {
922 copyNodeData(src, id_node, dst, id_node);
923 }
924 }
925 for (vtkIdType id_cell = 0; id_cell < src->GetNumberOfCells(); ++id_cell) {
926 vtkIdType id_new_cell = copyCell(src, id_cell, dst);
927 /*
928 if (id_cell == 2995) {
929 vtkUnstructuredGrid *G = dst;
930 EG_VTKSP(vtkIdList, stream);
931 vtkIdType type_cell = G->GetCellType(id_cell);
932 if (type_cell == VTK_POLYHEDRON) {
933 G->GetFaceStream(id_cell, stream);
934 } else {
935 G->GetCellPoints(id_cell, stream);
936 }
937 QList<int> ids;
938 vtkIdList2QContainer(stream, ids);
939 vtkIdType *pointer = stream->GetPointer(18);
940 cerr << "break" << endl;
941 }
942 */
943 if (copy_data) {
944 copyCellData(src, id_cell, dst, id_new_cell);
945 }
946 }
947 /*
948 if (!dbg) return;
949
950 {
951 vtkIdType id_cell = 2995;
952 vtkUnstructuredGrid *G = src;
953 EG_VTKSP(vtkIdList, stream);
954 vtkIdType type_cell = G->GetCellType(id_cell);
955 if (type_cell == VTK_POLYHEDRON) {
956 G->GetFaceStream(id_cell, stream);
957 } else {
958 G->GetCellPoints(id_cell, stream);
959 }
960 QList<int> ids;
961 vtkIdList2QContainer(stream, ids);
962 vtkIdType *pointer = stream->GetPointer(18);
963 cerr << "break" << endl;
964 }
965 {
966 vtkIdType id_cell = 2995;
967 vtkUnstructuredGrid *G = dst;
968 EG_VTKSP(vtkIdList, stream);
969 vtkIdType type_cell = G->GetCellType(id_cell);
970 if (type_cell == VTK_POLYHEDRON) {
971 G->GetFaceStream(id_cell, stream);
972 } else {
973 G->GetCellPoints(id_cell, stream);
974 }
975 QList<int> ids;
976 vtkIdList2QContainer(stream, ids);
977 vtkIdType *pointer = stream->GetPointer(18);
978 cerr << "break" << endl;
979 }
980 */
981 }
982
983 void EgVtkObject::makeCopyNoAlloc(vtkUnstructuredGrid *src, vtkUnstructuredGrid *dst)
984 {
985 for (vtkIdType id_node = 0; id_node < src->GetNumberOfPoints(); ++id_node) {
986 vec3_t x;
987 src->GetPoints()->GetPoint(id_node, x.data());
988 dst->GetPoints()->SetPoint(id_node, x.data());
989 copyNodeData(src, id_node, dst, id_node);
990 }
991 for (vtkIdType id_cell = 0; id_cell < src->GetNumberOfCells(); ++id_cell) {
992 vtkSmartPointer<vtkIdList> stream = vtkSmartPointer<vtkIdList>::New();
993 src->GetFaceStream(id_cell, stream);
994 auto type_cell = src->GetCellType(id_cell);
995 vtkIdType id_new_cell = dst->InsertNextCell(type_cell, stream);
996 copyCellData(src, id_cell, dst, id_new_cell);
997 }
998 }
999
1000 void EgVtkObject::reorientateFace(vtkUnstructuredGrid *grid, vtkIdType id_face)
1001 {
1002 EG_VTKDCC(vtkIntArray, cell_curdir, grid, "cell_curdir");
1003 EG_GET_CELL(id_face, grid);
1004 QVector<vtkIdType> new_pts(num_pts);
1005 for (int i = 0; i < num_pts; ++i) {
1006 new_pts[i] = pts[num_pts - i - 1];
1007 }
1008 if (cell_curdir->GetValue(id_face) == 0) {
1009 cell_curdir->SetValue(id_face, 1);
1010 } else {
1011 cell_curdir->SetValue(id_face, 0);
1012 }
1013 grid->ReplaceCell(id_face, num_pts, new_pts.data());
1014 }
1015
1016 void EgVtkObject::resetOrientation(vtkUnstructuredGrid *grid)
1017 {
1018 EG_VTKDCC(vtkIntArray, cell_orgdir, grid, "cell_orgdir");
1019 EG_VTKDCC(vtkIntArray, cell_curdir, grid, "cell_curdir");
1020 EG_VTKDCC(vtkIntArray, cell_voldir, grid, "cell_voldir");
1021 QVector<vtkIdType> faces;
1022 getAllSurfaceCells(faces, grid);
1023 foreach (vtkIdType id_face, faces) {
1024 if (cell_curdir->GetValue(id_face) != cell_orgdir->GetValue(id_face)) {
1025 reorientateFace(grid, id_face);
1026 cell_curdir->SetValue(id_face, cell_orgdir->GetValue(id_face));
1027 }
1028 cell_voldir->SetValue(id_face, 0);
1029 }
1030 }
1031
1032 vtkIdType EgVtkObject::findVolumeCell(vtkUnstructuredGrid *grid, vtkIdType id_surf, g2l_t _nodes, l2g_t cells, g2l_t _cells, l2l_t n2c)
1033 {
1034 EG_GET_CELL(id_surf, grid);
1035 QVector<QSet<int> > inters(num_pts-1);
1036 qcontIntersection(n2c[_nodes[pts[0]]], n2c[_nodes[pts[1]]], inters[0]);
1037 int i_pts = 2;
1038 while (i_pts < num_pts) {
1039 qcontIntersection(inters[i_pts-2], n2c[_nodes[pts[i_pts]]], inters[i_pts-1]);
1040 ++i_pts;
1041 }
1042 if (inters[num_pts-2].size() == 0) {
1043 return -1;
1044 } else if (inters[num_pts-2].size() > 2) {
1045 EG_BUG;
1046 }
1047 vtkIdType id_vol = -1;
1048 foreach (int i_cells, inters[num_pts-2]) {
1049 if (cells[i_cells] != id_surf) {
1050 id_vol = cells[i_cells];
1051 }
1052 }
1053 return id_vol;
1054 }
1055
1056 void EgVtkObject::setBoundaryCodes(const QSet<int> &bcs)
1057 {
1058 m_BoundaryCodes = bcs;
1059 }
1060
1061 QSet<int> EgVtkObject::getBoundaryCodes()
1062 {
1063 return m_BoundaryCodes;
1064 }
1065
1066 void EgVtkObject::createIndices(vtkUnstructuredGrid *grid)
1067 {
1068 if (!grid->GetCellData()->GetArray("cell_index")) {
1069 EG_VTKSP(vtkLongArray_t, var);
1070 var->SetName("cell_index");
1071 var->SetNumberOfValues(grid->GetNumberOfCells());
1072 grid->GetCellData()->AddArray(var);
1073 } else {
1074 EG_VTKDCC(vtkLongArray_t, var, grid, "cell_index");
1075 var->SetNumberOfValues(grid->GetNumberOfCells());
1076 }
1077 EG_VTKDCC(vtkLongArray_t, cell_index, grid, "cell_index");
1078 for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
1079 cell_index->SetValue(id_cell, id_cell);
1080 }
1081
1082 if (!grid->GetCellData()->GetArray("vtk_type")) {
1083 EG_VTKSP(vtkIntArray, var);
1084 var->SetName("vtk_type");
1085 var->SetNumberOfValues(grid->GetNumberOfCells());
1086 grid->GetCellData()->AddArray(var);
1087 } else {
1088 EG_VTKDCC(vtkIntArray, var, grid, "vtk_type");
1089 var->SetNumberOfValues(grid->GetNumberOfCells());
1090 }
1091 EG_VTKDCC(vtkIntArray, vtk_type, grid, "vtk_type");
1092 for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
1093 vtk_type->SetValue(id_cell, grid->GetCellType(id_cell));
1094 }
1095
1096 if (!grid->GetCellData()->GetArray("node_index")) {
1097 EG_VTKSP(vtkLongArray_t, var);
1098 var->SetName("node_index");
1099 var->SetNumberOfValues(grid->GetNumberOfPoints());
1100 grid->GetPointData()->AddArray(var);
1101 } else {
1102 EG_VTKDCC(vtkLongArray_t, var, grid, "node_index");
1103 var->SetNumberOfValues(grid->GetNumberOfPoints());
1104 }
1105 EG_VTKDCN(vtkLongArray_t, node_index, grid, "node_index");
1106 for (vtkIdType id_node = 0; id_node < grid->GetNumberOfPoints(); ++id_node) {
1107 node_index->SetValue(id_node, id_node);
1108 }
1109 }
1110
1111 BoundaryCondition EgVtkObject::getBC(int bc)
1112 {
1113 return GuiMainWindow::pointer()->getBC(bc);
1114 }
1115
1116 int EgVtkObject::getSet(QString group, QString key, int value, int& variable)
1117 {
1118 QSettings *qset = GuiMainWindow::settings();
1119 QString typed_key = "int/" + key;
1120 if(group!=QObject::tr("General")) qset->beginGroup(group);
1121 //if key=value pair not found in settings file, write it
1122 if (!qset->contains(typed_key)) qset->setValue(typed_key,value);
1123 //read key value from settings file and assign it to variable
1124 variable = (qset->value(typed_key,variable)).toInt();
1125 if(group!=QObject::tr("General")) qset->endGroup();
1126 return(variable);
1127 }
1128
1129 double EgVtkObject::getSet(QString group, QString key, double value, double& variable)
1130 {
1131 QSettings *qset = GuiMainWindow::settings();
1132 QString typed_key = "double/" + key;
1133 if(group!=QObject::tr("General")) qset->beginGroup(group);
1134 //if key=value pair not found in settings file, write it
1135 if (!qset->contains(typed_key)) qset->setValue(typed_key,value);
1136 //read key value from settings file and assign it to variable
1137 variable = (qset->value(typed_key,variable)).toDouble();
1138 if(group!=QObject::tr("General")) qset->endGroup();
1139 return(variable);
1140 }
1141
1142 bool EgVtkObject::getSet(QString group, QString key, bool value, bool& variable)
1143 {
1144 QSettings *qset = GuiMainWindow::settings();
1145 QString typed_key = "bool/" + key;
1146 if(group!=QObject::tr("General")) qset->beginGroup(group);
1147 Qt::CheckState state = (Qt::CheckState) ( value ? 2 : 0 );
1148 //if key=value pair not found in settings file, write it
1149 if (!qset->contains(typed_key)) qset->setValue(typed_key,state);
1150 //read key value from settings file and assign it to variable
1151 variable = (qset->value(typed_key,variable)).toBool();
1152 if(group!=QObject::tr("General")) qset->endGroup();
1153 return(variable);
1154 }
1155
1156 QString EgVtkObject::getSet(QString group, QString key, QString value, QString& variable)
1157 {
1158 QSettings *qset = GuiMainWindow::settings();
1159 QString typed_key;
1160 typed_key = QObject::tr("QString/") + key;
1161 if (group != QObject::tr("General")) qset->beginGroup(group);
1162 //if key=value pair not found in settings file, write it
1163 if (!qset->contains(typed_key)) qset->setValue(typed_key, value);
1164 //read key value from settings file and assign it to variable
1165 variable = (qset->value(typed_key)).toString();
1166 if (group != QObject::tr("General")) qset->endGroup();
1167 return(variable);
1168 }
1169
1170 QString EgVtkObject::getSet(QString group, QString key, QString value, QString& variable, int type)
1171 {
1172 QSettings *qset = GuiMainWindow::settings();
1173 QString typed_key;
1174 if (type == 0) {
1175 typed_key = QObject::tr("QString/") + key;
1176 }
1177 else if (type == 1) {
1178 typed_key = QObject::tr("Filename/") + key;
1179 }
1180 else {
1181 typed_key = QObject::tr("Directory/") + key;
1182 }
1183 if (group != QObject::tr("General")) qset->beginGroup(group);
1184 //if key=value pair not found in settings file, write it
1185 if (!qset->contains(typed_key)) qset->setValue(typed_key, value);
1186 //read key value from settings file and assign it to variable
1187 variable = (qset->value(typed_key)).toString();
1188 if (group != QObject::tr("General")) qset->endGroup();
1189 return(variable);
1190 }
1191
1192 void EgVtkObject::writeGrid(vtkUnstructuredGrid *grid, QString name)
1193 {
1194 QVector<vtkIdType> cells;
1195 getAllCells(cells, grid);
1196 name = GuiMainWindow::pointer()->getCwd() + "/" + name + ".vtu";
1197 writeCells(grid, cells, name);
1198 }
1199
1200 void EgVtkObject::getAllNodeDataNames(QVector<QString> &field_names, vtkUnstructuredGrid *grid)
1201 {
1202 int N = grid->GetPointData()->GetNumberOfArrays();
1203 field_names.resize(N);
1204 for (int i = 0; i < N; ++i) {
1205 field_names[i] = grid->GetPointData()->GetArrayName(i);
1206 }
1207 }
1208
1209 void EgVtkObject::getAllCellDataNames(QVector<QString> &field_names, vtkUnstructuredGrid *grid)
1210 {
1211 int N = grid->GetCellData()->GetNumberOfArrays();
1212 field_names.resize(N);
1213 for (int i = 0; i < N; ++i) {
1214 field_names[i] = grid->GetCellData()->GetArrayName(i);
1215 }
1216 }
1217
1218 QString EgVtkObject::stripFromExtension(QString file_name)
1219 {
1220 int i = file_name.size() - 1;
1221 while ((i > 0) && (file_name[i] != '.') && (file_name[i] != '/') && (file_name[i] != '\\')) {
1222 --i;
1223 }
1224 if (file_name[i] == '.') {
1225 return file_name.left(i);
1226 }
1227 return file_name;
1228 }
1229
1230 QString EgVtkObject::getExtension(QString file_name)
1231 {
1232 int i = file_name.size();
1233 while ((i > 0) && (file_name[i] != '.') && (file_name[i] != '/') && (file_name[i] != '\\')) {
1234 --i;
1235 }
1236 if (file_name[i] == '.') {
1237 return (file_name.right(file_name.size() - i - 1)).toLower();
1238 }
1239 return "";
1240 }
1241
1242 ///////////////////////////////////////////
1243
1244 void EgVtkObject::getFaceOfCell(vtkUnstructuredGrid *grid, vtkIdType id_cell, int i_face, QVector<vtkIdType> &ids)
1245 {
1246 vtkIdType type_cell = grid->GetCellType(id_cell);
1247 ids.clear();
1248 EG_VTKSP(vtkIdList, stream);
1249 QVector<vtkIdType> pts;
1250 if (type_cell == VTK_POLYHEDRON) {
1251 grid->GetFaceStream(id_cell, stream);
1252 } else {
1253 grid->GetCellPoints(id_cell, stream);
1254 }
1255 vtkIdList2QContainer(stream, pts);
1256
1257 if (type_cell == VTK_TETRA) {
1258 if (i_face == 0) { ids.resize(3); ids[0] = pts[2]; ids[1] = pts[1]; ids[2] = pts[0]; }
1259 else if (i_face == 1) { ids.resize(3); ids[0] = pts[1]; ids[1] = pts[3]; ids[2] = pts[0]; }
1260 else if (i_face == 2) { ids.resize(3); ids[0] = pts[3]; ids[1] = pts[2]; ids[2] = pts[0]; }
1261 else if (i_face == 3) { ids.resize(3); ids[0] = pts[2]; ids[1] = pts[3]; ids[2] = pts[1]; }
1262
1263 } else if (type_cell == VTK_PYRAMID) {
1264 if (i_face == 0) { ids.resize(4); ids[0] = pts[0]; ids[1] = pts[3]; ids[2] = pts[2]; ids[3] = pts[1]; }
1265 else if (i_face == 1) { ids.resize(3); ids[0] = pts[0]; ids[1] = pts[1]; ids[2] = pts[4]; }
1266 else if (i_face == 2) { ids.resize(3); ids[0] = pts[1]; ids[1] = pts[2]; ids[2] = pts[4]; }
1267 else if (i_face == 3) { ids.resize(3); ids[0] = pts[2]; ids[1] = pts[3]; ids[2] = pts[4]; }
1268 else if (i_face == 4) { ids.resize(3); ids[0] = pts[3]; ids[1] = pts[0]; ids[2] = pts[4]; }
1269
1270 } else if (type_cell == VTK_WEDGE) {
1271 if (i_face == 0) { ids.resize(3); ids[0] = pts[0]; ids[1] = pts[1]; ids[2] = pts[2]; }
1272 else if (i_face == 1) { ids.resize(3); ids[0] = pts[3]; ids[1] = pts[5]; ids[2] = pts[4]; }
1273 else if (i_face == 2) { ids.resize(4); ids[0] = pts[3]; ids[1] = pts[4]; ids[2] = pts[1]; ids[3] = pts[0]; }
1274 else if (i_face == 3) { ids.resize(4); ids[0] = pts[1]; ids[1] = pts[4]; ids[2] = pts[5]; ids[3] = pts[2]; }
1275 else if (i_face == 4) { ids.resize(4); ids[0] = pts[0]; ids[1] = pts[2]; ids[2] = pts[5]; ids[3] = pts[3]; }
1276
1277 } else if (type_cell == VTK_HEXAHEDRON) {
1278 if (i_face == 0) { ids.resize(4); ids[0] = pts[0]; ids[1] = pts[3]; ids[2] = pts[2]; ids[3] = pts[1]; }
1279 else if (i_face == 1) { ids.resize(4); ids[0] = pts[4]; ids[1] = pts[5]; ids[2] = pts[6]; ids[3] = pts[7]; }
1280 else if (i_face == 2) { ids.resize(4); ids[0] = pts[0]; ids[1] = pts[1]; ids[2] = pts[5]; ids[3] = pts[4]; }
1281 else if (i_face == 3) { ids.resize(4); ids[0] = pts[3]; ids[1] = pts[7]; ids[2] = pts[6]; ids[3] = pts[2]; }
1282 else if (i_face == 4) { ids.resize(4); ids[0] = pts[0]; ids[1] = pts[4]; ids[2] = pts[7]; ids[3] = pts[3]; }
1283 else if (i_face == 5) { ids.resize(4); ids[0] = pts[1]; ids[1] = pts[2]; ids[2] = pts[6]; ids[3] = pts[5]; }
1284
1285 } else if (type_cell == VTK_POLYHEDRON) {
1286 vtkIdType num_faces = pts[0];
1287 if (i_face >= num_faces) {
1288 EG_BUG;
1289 }
1290 int id = 1;
1291 for (int i = 0; i < i_face; ++i) {
1292 id += pts[id] + 1;
1293 }
1294 ids.resize(pts[id]);
1295 for (int i = 0; i < ids.size(); ++i) {
1296 ids[i] = pts[id + 1 + i];
1297 }
1298
1299 } else {
1300 EG_BUG; // not implemented
1301 }
1302 }
1303
1304 vec3_t EgVtkObject::getNormalOfCell(vtkUnstructuredGrid *grid, vtkIdType id_cell, int i_face)
1305 {
1306 QVector<vtkIdType> ids;
1307 getFaceOfCell(grid, id_cell, i_face, ids);
1308 if (ids.size() == 0) {
1309 EG_BUG;
1310 }
1311 QVector<vec3_t> x(ids.size() + 1);
1312 vec3_t xc(0,0,0);
1313 for (int i = 0; i < ids.size(); ++i) {
1314 grid->GetPoint(ids[i], x[i].data());
1315 xc += x[i];
1316 }
1317 x[ids.size()] = x[0];
1318 xc *= 1.0/ids.size();
1319 vec3_t n(0,0,0);
1320 for (int i = 0; i < ids.size(); ++i) {
1321 vec3_t u = x[i] - xc;
1322 vec3_t v = x[i+1] - xc;
1323 n += 0.5*u.cross(v);
1324 }
1325 return n;
1326 }
1327
1328 vec3_t EgVtkObject::getCentreOfCellFace(vtkUnstructuredGrid *grid, vtkIdType id_cell, int i_face)
1329 {
1330 QVector<vtkIdType> ids;
1331 getFaceOfCell(grid, id_cell, i_face, ids);
1332 if (ids.size() == 0) {
1333 EG_BUG;
1334 }
1335 vec3_t xc(0,0,0), x;
1336 for (int i = 0; i < ids.size(); ++i) {
1337 grid->GetPoint(ids[i], x.data());
1338 xc += x;
1339 }
1340 xc *= 1.0/ids.size();
1341 return xc;
1342 }
1343
1344 void EgVtkObject::getEdgeOfCell(vtkUnstructuredGrid *grid, vtkIdType id_cell, int i_edge, QVector<vtkIdType> &ids)
1345 {
1346 ids.clear();
1347 EG_GET_CELL(id_cell, grid);
1348 if (type_cell == VTK_TETRA) {
1349 ids.resize(2);
1350 if (i_edge == 0) { ids[0] = pts[0]; ids[1] = pts[1]; }
1351 else if (i_edge == 1) { ids[0] = pts[0]; ids[1] = pts[2]; }
1352 else if (i_edge == 2) { ids[0] = pts[0]; ids[1] = pts[3]; }
1353 else if (i_edge == 3) { ids[0] = pts[1]; ids[1] = pts[2]; }
1354 else if (i_edge == 4) { ids[0] = pts[1]; ids[1] = pts[3]; }
1355 else if (i_edge == 5) { ids[0] = pts[2]; ids[1] = pts[3]; }
1356 } else {
1357 EG_BUG; // not implemented
1358 }
1359 }
1360
1361 bool EgVtkObject::saveGrid(vtkUnstructuredGrid* a_grid, QString file_name)
1362 {
1363 file_name += ".vtu";
1364 addVtkTypeInfo(a_grid);
1365 createIndices(a_grid);
1366 EG_VTKSP(vtkXMLUnstructuredGridWriter,vtu);
1367 vtu->SetFileName(qPrintable(file_name));
1368 vtu->SetDataModeToBinary();
1369 vtu->SetInputData(a_grid);
1370 vtu->Write();
1371 if(vtu->GetErrorCode()) {
1372 return false;
1373 }
1374 else {
1375 return true;
1376 }
1377 }
1378
1379 void EgVtkObject::addVtkTypeInfo(vtkUnstructuredGrid* a_grid)
1380 {
1381 EG_VTKSP(vtkIntArray, vtk_type);
1382 vtk_type->SetName("vtk_type");
1383 vtk_type->SetNumberOfValues(a_grid->GetNumberOfCells());
1384 //EG_VTKDCC(vtkDoubleArray, cell_VA, a_grid, "cell_VA");
1385 for (vtkIdType cellId = 0; cellId < a_grid->GetNumberOfCells(); ++cellId) {
1386 vtk_type->SetValue(cellId, a_grid->GetCellType(cellId));
1387 //cell_VA->SetValue(cellId, GeometryTools::cellVA(a_grid, cellId, true));
1388 }
1389 a_grid->GetCellData()->AddArray(vtk_type);
1390 }
1391
1392 vtkIdType EgVtkObject::addGrid(vtkUnstructuredGrid *main_grid, vtkUnstructuredGrid *grid_to_add, vtkIdType offset)
1393 {
1394 EG_VTKSP(vtkUnstructuredGrid, main_copy);
1395 makeCopy(main_grid, main_copy);
1396 allocateGrid(main_grid, main_copy->GetNumberOfCells() + grid_to_add->GetNumberOfCells(), main_copy->GetNumberOfPoints() + grid_to_add->GetNumberOfPoints());
1397 makeCopyNoAlloc(main_copy, main_grid);
1398 for (vtkIdType id_node = 0; id_node < grid_to_add->GetNumberOfPoints(); ++id_node) {
1399 vec3_t x;
1400 grid_to_add->GetPoints()->GetPoint(id_node, x.data());
1401 main_grid->GetPoints()->SetPoint(offset + id_node, x.data());
1402 copyNodeData(grid_to_add, id_node, main_grid, offset + id_node);
1403 }
1404 for (vtkIdType id_cell = 0; id_cell < grid_to_add->GetNumberOfCells(); ++id_cell) {
1405 /*
1406 vtkIdType N_pts, *pts;
1407 vtkIdType type_cell = grid_to_add->GetCellType(id_cell);
1408 grid_to_add->GetCellPoints(id_cell, N_pts, pts);
1409 QVector <vtkIdType> new_pts(N_pts);
1410 for(int i=0;i<N_pts;i++) new_pts[i] = offset + pts[i];
1411 */
1412 vtkIdType id_new_cell = copyCell(grid_to_add, id_cell, main_grid, offset);
1413 copyCellData(grid_to_add, id_cell, main_grid, id_new_cell);
1414 }
1415 return( offset + grid_to_add->GetNumberOfPoints() );
1416 }
1417
1418 QSet<int> EgVtkObject::getAllBoundaryCodes(vtkUnstructuredGrid *grid)
1419 {
1420 EG_VTKDCC(vtkIntArray, cell_code, grid, "cell_code");
1421 QSet<int> bcs;
1422 for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
1423 if (isSurface(id_cell, grid)) {
1424 bcs.insert(cell_code->GetValue(id_cell));
1425 }
1426 }
1427 return bcs;
1428 }
1429
1430 bool EgVtkObject::cellContainsNode(vtkUnstructuredGrid *grid, vtkIdType id_cell, vtkIdType id_node)
1431 {
1432 /*
1433 vtkIdType num_pts, *pts;
1434 grid->GetCellPoints(id_cell, num_pts, pts);
1435 for (int i = 0; i < num_pts; ++i) {
1436 if (pts[i] == id_node) {
1437 return true;
1438 }
1439 }
1440 return false;
1441 */
1442
1443 QVector<vtkIdType> pts;
1444 getPointsOfCell(grid, id_cell, pts);
1445 if (pts.contains(id_node)) {
1446 return true;
1447 }
1448 return false;
1449
1450 }
1451
1452 void EgVtkObject::createPolyDataC2C(vtkPolyData *poly_data, QVector<QVector<vtkIdType> > &c2c)
1453 {
1454 c2c.clear();
1455 c2c.resize(poly_data->GetNumberOfCells());
1456 QVector<QSet<vtkIdType> > n2c;
1457 createPolyDataN2C(poly_data, n2c);
1458 EG_FORALL_CELLS(id_poly, poly_data) {
1459 EG_GET_CELL(id_poly, poly_data);
1460 QSet<vtkIdType> ids;
1461 for (int i = 0; i < num_pts; ++i) {
1462 vtkIdType p1 = pts[i];
1463 vtkIdType p2 = pts[0];
1464 if (i < num_pts - 1) {
1465 p2 = pts[i+1];
1466 }
1467 QSet<vtkIdType> shared_polys = n2c[p1];
1468 shared_polys.intersect(n2c[p2]);
1469 if (shared_polys.size() == 0) {
1470 EG_BUG;
1471 }
1472 if (shared_polys.size() > 2) {
1473 EG_BUG;
1474 }
1475 foreach (vtkIdType id_neigh, shared_polys) {
1476 if (id_neigh != id_poly) {
1477 ids.insert(id_neigh);
1478 }
1479 }
1480 }
1481 {
1482 c2c[id_poly].resize(ids.size());
1483 int i = 0;
1484 foreach (vtkIdType id_neigh, ids) {
1485 c2c[id_poly][i] = id_neigh;
1486 ++i;
1487 }
1488 }
1489 }
1490 }
1491
1492 void EgVtkObject::createPolyDataN2C(vtkPolyData *poly_data, QVector<QSet<vtkIdType> > &n2c)
1493 {
1494 n2c.clear();
1495 n2c.resize(poly_data->GetNumberOfPoints());
1496 EG_FORALL_CELLS(id_face, poly_data) {
1497 EG_GET_CELL(id_face, poly_data);
1498 for (int i = 0; i < num_pts; ++i) {
1499 n2c[pts[i]].insert(id_face);
1500 }
1501 }
1502 }
1503
1504 void EgVtkObject::createPolyDataN2N(vtkPolyData *poly_data, QVector<QSet<vtkIdType> > &n2n)
1505 {
1506 n2n.clear();
1507 n2n.resize(poly_data->GetNumberOfPoints());
1508 EG_FORALL_CELLS(id_face, poly_data) {
1509 EG_GET_CELL(id_face, poly_data);
1510 n2n[pts[0]].insert(pts[num_pts - 1]);
1511 n2n[pts[num_pts-1]].insert(pts[0]);
1512 for (int i = 1; i < num_pts; ++i) {
1513 n2n[pts[i]].insert(pts[i-1]);
1514 n2n[pts[i-1]].insert(pts[i]);
1515 }
1516 }
1517 }
1518
1519 void EgVtkObject::checkGridConsitency(vtkUnstructuredGrid *grid)
1520 {
1521 vtkIdType num_cells = grid->GetNumberOfCells();
1522 vtkIdType num_nodes = grid->GetNumberOfPoints();
1523 for (vtkIdType id_cell = 0; id_cell < num_cells; ++id_cell) {
1524 QList<vtkIdType> pts;
1525 getPointsOfCell(grid, id_cell, pts);
1526 foreach (vtkIdType id_node, pts) {
1527 if (id_node >= num_nodes) {
1528 EG_BUG;
1529 }
1530 }
1531 }
1532 }
1533
1534 QString EgVtkObject::getXmlSection(QString name)
1535 {
1536 return GuiMainWindow::pointer()->getXmlSection(name);
1537 }
1538
1539 void EgVtkObject::getSurfaceCells(int bc, QVector<vtkIdType> &cells, vtkUnstructuredGrid *grid)
1540 {
1541 int N = 0;
1542 EG_VTKDCC(vtkIntArray, cell_code, grid, "cell_code");
1543 for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
1544 if (isSurface(id_cell, grid)) {
1545 if (bc == cell_code->GetValue(id_cell)) {
1546 ++N;
1547 }
1548 }
1549 }
1550 cells.resize(N);
1551 N = 0;
1552 for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
1553 if (isSurface(id_cell, grid)) {
1554 if (bc == cell_code->GetValue(id_cell)) {
1555 cells[N] = id_cell;
1556 ++N;
1557 }
1558 }
1559 }
1560 }
1561
This is a mirror repository of the new location for engrid@github