?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
DeleteSetOfPoints function working
[engrid.git] / polymesh.cpp
blobadfe0bf4cf8b78087c699a9425b601776e4502d9
1 //
2 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3 // + +
4 // + This file is part of enGrid. +
5 // + +
6 // + Copyright 2008,2009 Oliver Gloth +
7 // + +
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. +
12 // + +
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. +
17 // + +
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/>. +
20 // + +
21 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
22 //
23 #include "polymesh.h"
24
25 uint qHash(PolyMesh::node_t N)
26 {
27 uint h = 0;
28 if (N.type == PolyMesh::node) {
29 h = N.idx;
30 } else {
31 if (N.type == PolyMesh::cell) {
32 vtkIdType N_pts, *pts;
33 N.poly->grid->GetCellPoints(N.idx,N_pts,pts);
34 for (int i = 0; i < N_pts; ++i) h += pts[i];
35 } else {
36 QList<vtkIdType> nds;
37 if (N.type == PolyMesh::edge) N.poly->getEdge(N.idx,N.subidx,nds);
38 if (N.type == PolyMesh::face) N.poly->getFace(N.idx,N.subidx,nds);
39 foreach (vtkIdType p, nds) h += p;
40 };
41 };
42 return h;
43 };
44
45
46 bool PolyMesh::node_t::operator==(const node_t N) const
47 {
48 QList<vtkIdType> nds;
49 if (type == edge) {
50 poly->getEdge(idx,subidx,nds);
51 };
52 if (type == face) {
53 poly->getFace(idx,subidx,nds);
54 };
55 bool match = ((type==N.type) && (idx==N.idx) && (subidx==N.subidx));
56 if (!match && (type == N.type) && (idx != N.idx)) {
57 if ((type == edge) || (type == face)) {
58 match = true;
59 QList<vtkIdType> nds1, nds2;
60 if (type == edge) {
61 poly->getEdge(idx,subidx,nds1);
62 poly->getEdge(N.idx,N.subidx,nds2);
63 };
64 if (type == face) {
65 poly->getFace(idx,subidx,nds1);
66 poly->getFace(N.idx,N.subidx,nds2);
67 };
68 foreach (vtkIdType p, nds1) {
69 if (!nds2.contains(p)) {
70 match = false;
71 break;
72 };
73 };
74 };
75 };
76 return match;
77 };
78
79 ostream& operator<<(ostream &s, PolyMesh::node_t N)
80 {
81 if (N.type == PolyMesh::node) s << "node: ";
82 if (N.type == PolyMesh::edge) s << "edge: ";
83 if (N.type == PolyMesh::face) s << "face: ";
84 if (N.type == PolyMesh::cell) s << "cell: ";
85 s << N.idx << ',' << N.subidx;
86 if (N.type != PolyMesh::node) {
87 s << "(";
88 if (N.type == PolyMesh::cell) {
89 vtkIdType N_pts, *pts;
90 N.poly->grid->GetCellPoints(N.idx,N_pts,pts);
91 for (int i = 0; i < N_pts-1; ++i) s << pts[i] << ',';
92 s << pts[N_pts-1];
93 } else {
94 QList<vtkIdType> nds;
95 if (N.type == PolyMesh::edge) {
96 N.poly->getEdge(N.idx,N.subidx,nds);
97 } else {
98 N.poly->getFace(N.idx,N.subidx,nds);
99 };
100 QList<vtkIdType>::iterator i = nds.begin();
101 s << *i;
102 ++i;
103 while (i != nds.end()) {
104 s << ',' << *i;
105 ++i;
106 };
107 };
108 s << ")";
109 };
110 return s;
111 };
112
113 ostream& operator<<(ostream &s, PolyMesh::face_t F)
114 {
115 s << "bc: " << F.bc << " owner: " << F.owner << " neighbour: " << F.neighbour;
116 for (int i = 0; i < F.node.size(); ++i) {
117 s << "\n " << F.node[i];
118 };
119 return s;
120 };
121
122 void PolyMesh::face_t::checkOrientation()
123 {
124 if (owner > neighbour) {
125 QVector<node_t> node_copy(node.size());
126 qCopy(node.begin(),node.end(),node_copy.begin());
127 for (int i = 0; i < node.size(); ++i) {
128 node[i] = node_copy[node.size()-i-1];
129 };
130 swap(owner,neighbour);
131 };
132 };
133
134 bool PolyMesh::face_t::operator<(const face_t &F) const
135 {
136 bool less = false;
137 if (bc < F.bc) {
138 less = true;
139 } else if (bc == F.bc) {
140 if (owner < F.owner) {
141 less = true;
142 } else if (owner == F.owner) {
143 if (neighbour < F.neighbour) {
144 less = true;
145 };
146 };
147 };
148
149 return less;
150 };
151
152 vec3_t PolyMesh::face_t::normal()
153 {
154 vec3_t x(0,0,0);
155 if (node.size() < 3) EG_BUG;
156 for (int i = 0; i < node.size(); ++i) {
157 x += vec3_t(node[i].x,node[i].y,node[i].z);
158 };
159 x *= 1.0/node.size();
160 vec3_t n(0,0,0);
161 for (int i = 0; i < node.size()-1; ++i) {
162 vec3_t a = vec3_t(node[i].x,node[i].y,node[i].z)-x;
163 vec3_t b = vec3_t(node[i+1].x,node[i+1].y,node[i+1].z)-x;
164 n += 0.5*(a.cross(b));
165 };
166 vec3_t a = vec3_t(node[node.size()-1].x,node[node.size()-1].y,node[node.size()-1].z) - x;
167 vec3_t b = vec3_t(node[0].x,node[0].y,node[0].z) - x;
168 n += 0.5*(a.cross(b));
169 return n;
170 };
171
172 void PolyMesh::getFace(vtkIdType idx, int subidx, QList<vtkIdType> &nodes)
173 {
174 nodes.clear();
175 vtkIdType N_pts, *pts;
176 grid->GetCellPoints(idx, N_pts, pts);
177 vtkIdType type_cell = grid->GetCellType(idx);
178 if (type_cell == VTK_TETRA) {
179 if (subidx == 0) {
180 nodes.append(pts[2]);
181 nodes.append(pts[1]);
182 nodes.append(pts[0]);
183 };
184 if (subidx == 1) {
185 nodes.append(pts[1]);
186 nodes.append(pts[3]);
187 nodes.append(pts[0]);
188 };
189 if (subidx == 2) {
190 nodes.append(pts[3]);
191 nodes.append(pts[2]);
192 nodes.append(pts[0]);
193 };
194 if (subidx == 3) {
195 nodes.append(pts[2]);
196 nodes.append(pts[3]);
197 nodes.append(pts[1]);
198 };
199 };
200 if (type_cell == VTK_WEDGE) {
201 if (subidx == 0) {
202 nodes.append(pts[0]);
203 nodes.append(pts[1]);
204 nodes.append(pts[2]);
205 };
206 if (subidx == 1) {
207 nodes.append(pts[3]);
208 nodes.append(pts[5]);
209 nodes.append(pts[4]);
210 };
211 if (subidx == 2) {
212 nodes.append(pts[3]);
213 nodes.append(pts[4]);
214 nodes.append(pts[1]);
215 nodes.append(pts[0]);
216 };
217 if (subidx == 3) {
218 nodes.append(pts[1]);
219 nodes.append(pts[4]);
220 nodes.append(pts[2]);
221 nodes.append(pts[5]);
222 };
223 if (subidx == 4) {
224 nodes.append(pts[0]);
225 nodes.append(pts[2]);
226 nodes.append(pts[5]);
227 nodes.append(pts[3]);
228 };
229 };
230 if (type_cell == VTK_HEXAHEDRON) {
231 if (subidx == 0) {
232 nodes.append(pts[0]);
233 nodes.append(pts[3]);
234 nodes.append(pts[2]);
235 nodes.append(pts[1]);
236 };
237 if (subidx == 1) {
238 nodes.append(pts[4]);
239 nodes.append(pts[5]);
240 nodes.append(pts[6]);
241 nodes.append(pts[7]);
242 };
243 if (subidx == 2) {
244 nodes.append(pts[0]);
245 nodes.append(pts[1]);
246 nodes.append(pts[5]);
247 nodes.append(pts[4]);
248 };
249 if (subidx == 3) {
250 nodes.append(pts[3]);
251 nodes.append(pts[7]);
252 nodes.append(pts[6]);
253 nodes.append(pts[2]);
254 };
255 if (subidx == 4) {
256 nodes.append(pts[0]);
257 nodes.append(pts[4]);
258 nodes.append(pts[7]);
259 nodes.append(pts[3]);
260 };
261 if (subidx == 5) {
262 nodes.append(pts[1]);
263 nodes.append(pts[2]);
264 nodes.append(pts[6]);
265 nodes.append(pts[5]);
266 };
267 };
268 };
269
270 void PolyMesh::getEdge(vtkIdType idx, int subidx, QList<vtkIdType> &nodes)
271 {
272 nodes.clear();
273 vtkIdType N_pts, *pts;
274 grid->GetCellPoints(idx, N_pts, pts);
275 vtkIdType type_cell = grid->GetCellType(idx);
276 if (type_cell == VTK_TETRA) {
277 if (subidx == 0) {
278 nodes.append(pts[0]);
279 nodes.append(pts[1]);
280 };
281 if (subidx == 1) {
282 nodes.append(pts[0]);
283 nodes.append(pts[2]);
284 };
285 if (subidx == 2) {
286 nodes.append(pts[0]);
287 nodes.append(pts[3]);
288 };
289 if (subidx == 3) {
290 nodes.append(pts[1]);
291 nodes.append(pts[2]);
292 };
293 if (subidx == 4) {
294 nodes.append(pts[1]);
295 nodes.append(pts[3]);
296 };
297 if (subidx == 5) {
298 nodes.append(pts[2]);
299 nodes.append(pts[3]);
300 };
301 };
302 if (type_cell == VTK_WEDGE) {
303 if (subidx == 0) {
304 nodes.append(pts[0]);
305 nodes.append(pts[1]);
306 };
307 if (subidx == 1) {
308 nodes.append(pts[0]);
309 nodes.append(pts[2]);
310 };
311 if (subidx == 2) {
312 nodes.append(pts[0]);
313 nodes.append(pts[3]);
314 };
315 if (subidx == 3) {
316 nodes.append(pts[1]);
317 nodes.append(pts[2]);
318 };
319 if (subidx == 4) {
320 nodes.append(pts[1]);
321 nodes.append(pts[4]);
322 };
323 if (subidx == 5) {
324 nodes.append(pts[2]);
325 nodes.append(pts[5]);
326 };
327 if (subidx == 6) {
328 nodes.append(pts[3]);
329 nodes.append(pts[4]);
330 };
331 if (subidx == 7) {
332 nodes.append(pts[3]);
333 nodes.append(pts[5]);
334 };
335 if (subidx == 8) {
336 nodes.append(pts[4]);
337 nodes.append(pts[5]);
338 };
339 };
340 if (type_cell == VTK_HEXAHEDRON) {
341 if (subidx == 0) {
342 nodes.append(pts[0]);
343 nodes.append(pts[1]);
344 };
345 if (subidx == 1) {
346 nodes.append(pts[0]);
347 nodes.append(pts[3]);
348 };
349 if (subidx == 2) {
350 nodes.append(pts[0]);
351 nodes.append(pts[4]);
352 };
353 if (subidx == 3) {
354 nodes.append(pts[1]);
355 nodes.append(pts[2]);
356 };
357 if (subidx == 4) {
358 nodes.append(pts[1]);
359 nodes.append(pts[5]);
360 };
361 if (subidx == 5) {
362 nodes.append(pts[2]);
363 nodes.append(pts[3]);
364 };
365 if (subidx == 6) {
366 nodes.append(pts[2]);
367 nodes.append(pts[6]);
368 };
369 if (subidx == 7) {
370 nodes.append(pts[3]);
371 nodes.append(pts[7]);
372 };
373 if (subidx == 8) {
374 nodes.append(pts[4]);
375 nodes.append(pts[5]);
376 };
377 if (subidx == 9) {
378 nodes.append(pts[4]);
379 nodes.append(pts[7]);
380 };
381 if (subidx == 10) {
382 nodes.append(pts[5]);
383 nodes.append(pts[6]);
384 };
385 if (subidx == 11) {
386 nodes.append(pts[6]);
387 nodes.append(pts[7]);
388 };
389 };
390 };
391
392 PolyMesh::PolyMesh(vtkUnstructuredGrid *a_grid, bool dual_mesh)
393 {
394 dbg = false;
395 dual = dual_mesh;
396 grid = a_grid;
397 weight.fill(1.0,grid->GetNumberOfPoints());
398 pass1();
399 pass2();
400 pass3();
401 };
402
403 int PolyMesh::pcIdxNode(vtkIdType id_node)
404 {
405 if (node2pc[id_node] != -1) {
406 return node2pc[id_node];
407 };
408 node2pc[id_node] = id_newpc;
409 ++id_newpc;
410 return node2pc[id_node];
411 };
412
413 int PolyMesh::pcIdxCell(vtkIdType id_cell)
414 {
415 if (cell2pc[id_cell] != -1) {
416 return cell2pc[id_cell];
417 };
418 cell2pc[id_cell] = id_newpc;
419 ++id_newpc;
420 return cell2pc[id_cell];
421 };
422
423 void PolyMesh::pass1Tetras()
424 {
425 EG_VTKDCC(vtkIntArray, bc, grid, "cell_code");
426 face_t F;
427 for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
428 vtkIdType type_cell = grid->GetCellType(id_cell);
429 if (type_cell == VTK_TETRA) {
430 vtkIdType N_pts, *pts;
431 grid->GetCellPoints(id_cell, N_pts, pts);
432
433 if (dual) {
434
435 // edge 0 -> 1
436 F.node.resize(4);
437 F.node[0] = node_t(this, face, id_cell, 0);
438 F.node[1] = node_t(this, cell, id_cell, 0);
439 F.node[2] = node_t(this, face, id_cell, 1);
440 F.node[3] = node_t(this, edge, id_cell, 0);
441 F.owner = pcIdxNode(pts[0]);
442 F.neighbour = pcIdxNode(pts[1]);
443 F.bc = 0;
444 F.checkOrientation();
445 face_list.append(F);
446
447 // edge 0 -> 2
448 F.node.resize(4);
449 F.node[0] = node_t(this, face, id_cell, 2);
450 F.node[1] = node_t(this, cell, id_cell, 0);
451 F.node[2] = node_t(this, face, id_cell, 0);
452 F.node[3] = node_t(this, edge, id_cell, 1);
453 F.owner = pcIdxNode(pts[0]);
454 F.neighbour = pcIdxNode(pts[2]);
455 F.bc = 0;
456 F.checkOrientation();
457 face_list.append(F);
458
459 // edge 0 -> 3
460 F.node.resize(4);
461 F.node[0] = node_t(this, face, id_cell, 1);
462 F.node[1] = node_t(this, cell, id_cell, 0);
463 F.node[2] = node_t(this, face, id_cell, 2);
464 F.node[3] = node_t(this, edge, id_cell, 2);
465 F.owner = pcIdxNode(pts[0]);
466 F.neighbour = pcIdxNode(pts[3]);
467 F.bc = 0;
468 F.checkOrientation();
469 face_list.append(F);
470
471 // edge 1 -> 2
472 F.node.resize(4);
473 F.node[0] = node_t(this, face, id_cell, 0);
474 F.node[1] = node_t(this, cell, id_cell, 0);
475 F.node[2] = node_t(this, face, id_cell, 3);
476 F.node[3] = node_t(this, edge, id_cell, 3);
477 F.owner = pcIdxNode(pts[1]);
478 F.neighbour = pcIdxNode(pts[2]);
479 F.bc = 0;
480 F.checkOrientation();
481 face_list.append(F);
482
483 // edge 1 -> 3
484 F.node.resize(4);
485 F.node[0] = node_t(this, face, id_cell, 3);
486 F.node[1] = node_t(this, cell, id_cell, 0);
487 F.node[2] = node_t(this, face, id_cell, 1);
488 F.node[3] = node_t(this, edge, id_cell, 4);
489 F.owner = pcIdxNode(pts[1]);
490 F.neighbour = pcIdxNode(pts[3]);
491 F.bc = 0;
492 F.checkOrientation();
493 face_list.append(F);
494
495 // edge 2 -> 3
496 F.node.resize(4);
497 F.node[0] = node_t(this, face, id_cell, 2);
498 F.node[1] = node_t(this, cell, id_cell, 0);
499 F.node[2] = node_t(this, face, id_cell, 3);
500 F.node[3] = node_t(this, edge, id_cell, 5);
501 F.owner = pcIdxNode(pts[2]);
502 F.neighbour = pcIdxNode(pts[3]);
503 F.bc = 0;
504 F.checkOrientation();
505 face_list.append(F);
506
507 // boundary face 0
508 {
509 vtkIdType id_bcell = c2c[id_cell][0];
510 if (id_bcell == -1) EG_BUG;
511 if (isSurface(id_bcell, grid)) {
512
513 F.node.resize(4);
514 F.node[0] = node_t(this, edge, id_cell, 1);
515 F.node[1] = node_t(this, face, id_cell, 0);
516 F.node[2] = node_t(this, edge, id_cell, 0);
517 F.node[3] = node_t(this, node, pts[0], 0);
518 F.owner = pcIdxNode(pts[0]);
519 F.neighbour = -1;
520 F.bc = bc->GetValue(id_bcell);
521 face_list.append(F);
522
523 F.node.resize(4);
524 F.node[0] = node_t(this, edge, id_cell, 0);
525 F.node[1] = node_t(this, face, id_cell, 0);
526 F.node[2] = node_t(this, edge, id_cell, 3);
527 F.node[3] = node_t(this, node, pts[1], 0);
528 F.owner = pcIdxNode(pts[1]);
529 F.neighbour = -1;
530 F.bc = bc->GetValue(id_bcell);
531 face_list.append(F);
532
533 F.node.resize(4);
534 F.node[0] = node_t(this, edge, id_cell, 3);
535 F.node[1] = node_t(this, face, id_cell, 0);
536 F.node[2] = node_t(this, edge, id_cell, 1);
537 F.node[3] = node_t(this, node, pts[2], 0);
538 F.owner = pcIdxNode(pts[2]);
539 F.neighbour = -1;
540 F.bc = bc->GetValue(id_bcell);
541 face_list.append(F);
542
543 };
544 };
545
546 // boundary face 1
547 {
548 vtkIdType id_bcell = c2c[id_cell][1];
549 if (id_bcell == -1) EG_BUG;
550 if (isSurface(id_bcell, grid)) {
551
552 F.node.resize(4);
553 F.node[0] = node_t(this, edge, id_cell, 0);
554 F.node[1] = node_t(this, face, id_cell, 1);
555 F.node[2] = node_t(this, edge, id_cell, 2);
556 F.node[3] = node_t(this, node, pts[0], 0);
557 F.owner = pcIdxNode(pts[0]);
558 F.neighbour = -1;
559 F.bc = bc->GetValue(id_bcell);
560 face_list.append(F);
561
562 F.node.resize(4);
563 F.node[0] = node_t(this, edge, id_cell, 4);
564 F.node[1] = node_t(this, face, id_cell, 1);
565 F.node[2] = node_t(this, edge, id_cell, 0);
566 F.node[3] = node_t(this, node, pts[1], 0);
567 F.owner = pcIdxNode(pts[1]);
568 F.neighbour = -1;
569 F.bc = bc->GetValue(id_bcell);
570 face_list.append(F);
571
572 F.node.resize(4);
573 F.node[0] = node_t(this, edge, id_cell, 2);
574 F.node[1] = node_t(this, face, id_cell, 1);
575 F.node[2] = node_t(this, edge, id_cell, 4);
576 F.node[3] = node_t(this, node, pts[3], 0);
577 F.owner = pcIdxNode(pts[3]);
578 F.neighbour = -1;
579 F.bc = bc->GetValue(id_bcell);
580 face_list.append(F);
581
582 };
583 };
584
585 // boundary face 2
586 {
587 vtkIdType id_bcell = c2c[id_cell][2];
588 if (id_bcell == -1) EG_BUG;
589 if (isSurface(id_bcell, grid)) {
590
591 F.node.resize(4);
592 F.node[0] = node_t(this, edge, id_cell, 2);
593 F.node[1] = node_t(this, face, id_cell, 2);
594 F.node[2] = node_t(this, edge, id_cell, 1);
595 F.node[3] = node_t(this, node, pts[0], 0);
596 F.owner = pcIdxNode(pts[0]);
597 F.neighbour = -1;
598 F.bc = bc->GetValue(id_bcell);
599 face_list.append(F);
600
601 F.node.resize(4);
602 F.node[0] = node_t(this, edge, id_cell, 1);
603 F.node[1] = node_t(this, face, id_cell, 2);
604 F.node[2] = node_t(this, edge, id_cell, 5);
605 F.node[3] = node_t(this, node, pts[2], 0);
606 F.owner = pcIdxNode(pts[2]);
607 F.neighbour = -1;
608 F.bc = bc->GetValue(id_bcell);
609 face_list.append(F);
610
611 F.node.resize(4);
612 F.node[0] = node_t(this, edge, id_cell, 5);
613 F.node[1] = node_t(this, face, id_cell, 2);
614 F.node[2] = node_t(this, edge, id_cell, 2);
615 F.node[3] = node_t(this, node, pts[3], 0);
616 F.owner = pcIdxNode(pts[3]);
617 F.neighbour = -1;
618 F.bc = bc->GetValue(id_bcell);
619 face_list.append(F);
620
621 };
622 };
623
624 // boundary face 3
625 {
626 vtkIdType id_bcell = c2c[id_cell][3];
627 if (id_bcell == -1) EG_BUG;
628 if (isSurface(id_bcell, grid)) {
629
630 F.node.resize(4);
631 F.node[0] = node_t(this, edge, id_cell, 3);
632 F.node[1] = node_t(this, face, id_cell, 3);
633 F.node[2] = node_t(this, edge, id_cell, 4);
634 F.node[3] = node_t(this, node, pts[1], 0);
635 F.owner = pcIdxNode(pts[1]);
636 F.neighbour = -1;
637 F.bc = bc->GetValue(id_bcell);
638 face_list.append(F);
639
640 F.node.resize(4);
641 F.node[0] = node_t(this, edge, id_cell, 5);
642 F.node[1] = node_t(this, face, id_cell, 3);
643 F.node[2] = node_t(this, edge, id_cell, 3);
644 F.node[3] = node_t(this, node, pts[2], 0);
645 F.owner = pcIdxNode(pts[2]);
646 F.neighbour = -1;
647 F.bc = bc->GetValue(id_bcell);
648 face_list.append(F);
649
650 F.node.resize(4);
651 F.node[0] = node_t(this, edge, id_cell, 4);
652 F.node[1] = node_t(this, face, id_cell, 3);
653 F.node[2] = node_t(this, edge, id_cell, 5);
654 F.node[3] = node_t(this, node, pts[3], 0);
655 F.owner = pcIdxNode(pts[3]);
656 F.neighbour = -1;
657 F.bc = bc->GetValue(id_bcell);
658 face_list.append(F);
659
660 };
661 };
662 } else {
663 };
664 };
665 };
666 };
667
668 void PolyMesh::pass1Prisms()
669 {
670 EG_VTKDCC(vtkIntArray, bc, grid, "cell_code");
671 face_t F;
672 for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
673 vtkIdType type_cell = grid->GetCellType(id_cell);
674 if (type_cell == VTK_WEDGE) {
675 vtkIdType N_pts, *pts, id_ncell;
676 bool f0_split = false;
677 bool f1_split = false;
678 grid->GetCellPoints(id_cell, N_pts, pts);
679
680 // face 0
681 id_ncell = c2c[id_cell][0];
682 F.bc = 0;
683 if (id_ncell == -1) EG_BUG;
684 if (isSurface(id_ncell, grid)) {
685 F.bc = bc->GetValue(id_ncell);
686 F.node.resize(3);
687 F.node[0] = node_t(this, node, pts[0], 0);
688 F.node[1] = node_t(this, node, pts[1], 0);
689 F.node[2] = node_t(this, node, pts[2], 0);
690 F.owner = pcIdxCell(id_cell);
691 F.neighbour = -1;
692 face_list.append(F);
693 } else {
694 if (grid->GetCellType(id_ncell) == VTK_WEDGE) {
695 F.bc = 0;
696 F.node.resize(3);
697 F.node[0] = node_t(this, node, pts[0], 0);
698 F.node[1] = node_t(this, node, pts[1], 0);
699 F.node[2] = node_t(this, node, pts[2], 0);
700 F.owner = pcIdxCell(id_cell);
701 F.neighbour = pcIdxCell(id_ncell);
702 if (F.owner < F.neighbour) face_list.append(F);
703 } else {
704 f0_split = true;
705
706 F.node.resize(4);
707 F.node[0] = node_t(this, node, pts[0], 0);
708 F.node[1] = node_t(this, edge, id_cell, 0);
709 F.node[2] = node_t(this, face, id_cell, 0);
710 F.node[3] = node_t(this, edge, id_cell, 1);
711 F.owner = pcIdxCell(id_cell);
712 F.neighbour = pcIdxNode(pts[0]);
713 F.bc = 0;
714 F.checkOrientation();
715 face_list.append(F);
716
717 F.node.resize(4);
718 F.node[0] = node_t(this, node, pts[1], 0);
719 F.node[1] = node_t(this, edge, id_cell, 3);
720 F.node[2] = node_t(this, face, id_cell, 0);
721 F.node[3] = node_t(this, edge, id_cell, 0);
722 F.owner = pcIdxCell(id_cell);
723 F.neighbour = pcIdxNode(pts[1]);
724 F.bc = 0;
725 F.checkOrientation();
726 face_list.append(F);
727
728 F.node.resize(4);
729 F.node[0] = node_t(this, node, pts[2], 0);
730 F.node[1] = node_t(this, edge, id_cell, 1);
731 F.node[2] = node_t(this, face, id_cell, 0);
732 F.node[3] = node_t(this, edge, id_cell, 3);
733 F.owner = pcIdxCell(id_cell);
734 F.neighbour = pcIdxNode(pts[2]);
735 F.bc = 0;
736 F.checkOrientation();
737 face_list.append(F);
738 };
739 };
740
741 // face 1
742 id_ncell = c2c[id_cell][1];
743 F.bc = 0;
744 if (id_ncell == -1) EG_BUG;
745 if (isSurface(id_ncell, grid)) {
746 F.bc = bc->GetValue(id_ncell);
747 F.node.resize(3);
748 F.node[0] = node_t(this, node, pts[3], 0);
749 F.node[1] = node_t(this, node, pts[5], 0);
750 F.node[2] = node_t(this, node, pts[4], 0);
751 F.owner = pcIdxCell(id_cell);
752 F.neighbour = -1;
753 face_list.append(F);
754 } else {
755 if (grid->GetCellType(id_ncell) == VTK_WEDGE) {
756 F.bc = 0;
757 F.node.resize(3);
758 F.node[0] = node_t(this, node, pts[3], 0);
759 F.node[1] = node_t(this, node, pts[5], 0);
760 F.node[2] = node_t(this, node, pts[4], 0);
761 F.owner = pcIdxCell(id_cell);
762 F.neighbour = pcIdxCell(id_ncell);
763 if (F.owner < F.neighbour) face_list.append(F);
764 } else {
765 f1_split = true;
766
767 F.node.resize(4);
768 F.node[0] = node_t(this, node, pts[3], 0);
769 F.node[1] = node_t(this, edge, id_cell, 7);
770 F.node[2] = node_t(this, face, id_cell, 1);
771 F.node[3] = node_t(this, edge, id_cell, 6);
772 F.owner = pcIdxCell(id_cell);
773 F.neighbour = pcIdxNode(pts[3]);
774 F.bc = 0;
775 F.checkOrientation();
776 face_list.append(F);
777
778 F.node.resize(4);
779 F.node[0] = node_t(this, node, pts[4], 0);
780 F.node[1] = node_t(this, edge, id_cell, 6);
781 F.node[2] = node_t(this, face, id_cell, 1);
782 F.node[3] = node_t(this, edge, id_cell, 8);
783 F.owner = pcIdxCell(id_cell);
784 F.neighbour = pcIdxNode(pts[4]);
785 F.bc = 0;
786 F.checkOrientation();
787 face_list.append(F);
788
789 F.node.resize(4);
790 F.node[0] = node_t(this, node, pts[5], 0);
791 F.node[1] = node_t(this, edge, id_cell, 8);
792 F.node[2] = node_t(this, face, id_cell, 1);
793 F.node[3] = node_t(this, edge, id_cell, 7);
794 F.owner = pcIdxCell(id_cell);
795 F.neighbour = pcIdxNode(pts[5]);
796 F.bc = 0;
797 F.checkOrientation();
798 face_list.append(F);
799 };
800 };
801
802 // face 2
803 id_ncell = c2c[id_cell][2];
804 F.bc = 0;
805 if (id_ncell == -1) EG_BUG;
806 if (isSurface(id_ncell, grid)) {
807 F.bc = bc->GetValue(id_ncell);
808 id_ncell = -1;
809 };
810 F.node.resize(4);
811 F.owner = pcIdxCell(id_cell);
812 if (!f0_split && !f1_split) {
813 F.node.resize(4);
814 F.node[0] = node_t(this, node, pts[1], 0);
815 F.node[1] = node_t(this, node, pts[0], 0);
816 F.node[2] = node_t(this, node, pts[3], 0);
817 F.node[3] = node_t(this, node, pts[4], 0);
818 };
819 if (f0_split && !f1_split) {
820 F.node.resize(5);
821 F.node[0] = node_t(this, node, pts[1], 0);
822 F.node[1] = node_t(this, edge, id_cell, 0);
823 F.node[2] = node_t(this, node, pts[0], 0);
824 F.node[3] = node_t(this, node, pts[3], 0);
825 F.node[4] = node_t(this, node, pts[4], 0);
826 };
827 if (!f0_split && f1_split) {
828 F.node.resize(5);
829 F.node[0] = node_t(this, node, pts[1], 0);
830 F.node[1] = node_t(this, node, pts[0], 0);
831 F.node[2] = node_t(this, node, pts[3], 0);
832 F.node[3] = node_t(this, edge, id_cell, 6);
833 F.node[4] = node_t(this, node, pts[4], 0);
834 };
835 if (f0_split && f1_split) {
836 F.node.resize(6);
837 F.node[0] = node_t(this, node, pts[1], 0);
838 F.node[1] = node_t(this, edge, id_cell, 0);
839 F.node[2] = node_t(this, node, pts[0], 0);
840 F.node[3] = node_t(this, node, pts[3], 0);
841 F.node[4] = node_t(this, edge, id_cell, 6);
842 F.node[5] = node_t(this, node, pts[4], 0);
843 };
844 if (id_ncell != -1) {
845 F.neighbour = pcIdxCell(id_ncell);
846 if (F.owner < F.neighbour) face_list.append(F);
847 } else {
848 F.neighbour = -1;
849 face_list.append(F);
850 };
851
852 // face 3
853 id_ncell = c2c[id_cell][3];
854 F.bc = 0;
855 if (id_ncell == -1) EG_BUG;
856 if (isSurface(id_ncell, grid)) {
857 F.bc = bc->GetValue(id_ncell);
858 id_ncell = -1;
859 };
860 F.owner = pcIdxCell(id_cell);
861 if (!f0_split && !f1_split) {
862 F.node.resize(4);
863 F.node[0] = node_t(this, node, pts[2], 0);
864 F.node[1] = node_t(this, node, pts[1], 0);
865 F.node[2] = node_t(this, node, pts[4], 0);
866 F.node[3] = node_t(this, node, pts[5], 0);
867 };
868 if (f0_split && !f1_split) {
869 F.node.resize(5);
870 F.node[0] = node_t(this, node, pts[2], 0);
871 F.node[1] = node_t(this, edge, id_cell, 3);
872 F.node[2] = node_t(this, node, pts[1], 0);
873 F.node[3] = node_t(this, node, pts[4], 0);
874 F.node[4] = node_t(this, node, pts[5], 0);
875 };
876 if (!f0_split && f1_split) {
877 F.node.resize(5);
878 F.node[0] = node_t(this, node, pts[2], 0);
879 F.node[1] = node_t(this, node, pts[1], 0);
880 F.node[2] = node_t(this, node, pts[4], 0);
881 F.node[3] = node_t(this, edge, id_cell, 8);
882 F.node[4] = node_t(this, node, pts[5], 0);
883 };
884 if (f0_split && f1_split) {
885 F.node.resize(6);
886 F.node[0] = node_t(this, node, pts[2], 0);
887 F.node[1] = node_t(this, edge, id_cell, 3);
888 F.node[2] = node_t(this, node, pts[1], 0);
889 F.node[3] = node_t(this, node, pts[4], 0);
890 F.node[4] = node_t(this, edge, id_cell, 8);
891 F.node[5] = node_t(this, node, pts[5], 0);
892 };
893
894 if (id_ncell != -1) {
895 F.neighbour = pcIdxCell(id_ncell);
896 if (F.owner < F.neighbour) face_list.append(F);
897 } else {
898 F.neighbour = -1;
899 face_list.append(F);
900 };
901
902 // face 4
903 id_ncell = c2c[id_cell][4];
904 F.bc = 0;
905 if (id_ncell == -1) EG_BUG;
906 if (isSurface(id_ncell, grid)) {
907 F.bc = bc->GetValue(id_ncell);
908 id_ncell = -1;
909 };
910 F.owner = pcIdxCell(id_cell);
911 if (!f0_split && !f1_split) {
912 F.node.resize(4);
913 F.node[0] = node_t(this, node, pts[0], 0);
914 F.node[1] = node_t(this, node, pts[2], 0);
915 F.node[2] = node_t(this, node, pts[5], 0);
916 F.node[3] = node_t(this, node, pts[3], 0);
917 F.owner = pcIdxCell(id_cell);
918 };
919 if (f0_split && !f1_split) {
920 F.node.resize(5);
921 F.node[0] = node_t(this, node, pts[0], 0);
922 F.node[1] = node_t(this, edge, id_cell, 1);
923 F.node[2] = node_t(this, node, pts[2], 0);
924 F.node[3] = node_t(this, node, pts[5], 0);
925 F.node[4] = node_t(this, node, pts[3], 0);
926 };
927 if (!f0_split && f1_split) {
928 F.node.resize(5);
929 F.node[0] = node_t(this, node, pts[0], 0);
930 F.node[1] = node_t(this, node, pts[2], 0);
931 F.node[2] = node_t(this, node, pts[5], 0);
932 F.node[3] = node_t(this, edge, id_cell, 7);
933 F.node[4] = node_t(this, node, pts[3], 0);
934 };
935 if (f0_split && f1_split) {
936 F.node.resize(6);
937 F.node[0] = node_t(this, node, pts[0], 0);
938 F.node[1] = node_t(this, edge, id_cell, 1);
939 F.node[2] = node_t(this, node, pts[2], 0);
940 F.node[3] = node_t(this, node, pts[5], 0);
941 F.node[4] = node_t(this, edge, id_cell, 7);
942 F.node[5] = node_t(this, node, pts[3], 0);
943 };
944 if (id_ncell != -1) {
945 F.neighbour = pcIdxCell(id_ncell);
946 if (F.owner < F.neighbour) face_list.append(F);
947 } else {
948 F.neighbour = -1;
949 face_list.append(F);
950 };
951
952 };
953 };
954 };
955
956 void PolyMesh::pass1Hexas()
957 {
958 EG_VTKDCC(vtkIntArray, bc, grid, "cell_code");
959 face_t F;
960 for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
961 vtkIdType type_cell = grid->GetCellType(id_cell);
962 if (type_cell == VTK_HEXAHEDRON) {
963 vtkIdType N_pts, *pts, id_ncell;
964 grid->GetCellPoints(id_cell, N_pts, pts);
965
966 // face 0
967 id_ncell = c2c[id_cell][0];
968 F.bc = 0;
969 if (id_ncell == -1) EG_BUG;
970 if (isSurface(id_ncell, grid)) {
971 F.bc = bc->GetValue(id_ncell);
972 id_ncell = -1;
973 };
974 F.owner = pcIdxCell(id_cell);
975 F.node.resize(4);
976 F.node[0] = node_t(this, node, pts[0], 0);
977 F.node[1] = node_t(this, node, pts[3], 0);
978 F.node[2] = node_t(this, node, pts[2], 0);
979 F.node[3] = node_t(this, node, pts[1], 0);
980 if (id_ncell != -1) {
981 F.neighbour = pcIdxCell(id_ncell);
982 if (F.owner < F.neighbour) face_list.append(F);
983 } else {
984 F.neighbour = -1;
985 face_list.append(F);
986 };
987
988 // face 1
989 id_ncell = c2c[id_cell][1];
990 F.bc = 0;
991 if (id_ncell == -1) EG_BUG;
992 if (isSurface(id_ncell, grid)) {
993 F.bc = bc->GetValue(id_ncell);
994 id_ncell = -1;
995 };
996 F.owner = pcIdxCell(id_cell);
997 F.node.resize(4);
998 F.node[0] = node_t(this, node, pts[4], 0);
999 F.node[1] = node_t(this, node, pts[5], 0);
1000 F.node[2] = node_t(this, node, pts[6], 0);
1001 F.node[3] = node_t(this, node, pts[7], 0);
1002 if (id_ncell != -1) {
1003 F.neighbour = pcIdxCell(id_ncell);
1004 if (F.owner < F.neighbour) face_list.append(F);
1005 } else {
1006 F.neighbour = -1;
1007 face_list.append(F);
1008 };
1009
1010 // face 2
1011 id_ncell = c2c[id_cell][2];
1012 F.bc = 0;
1013 if (id_ncell == -1) EG_BUG;
1014 if (isSurface(id_ncell, grid)) {
1015 F.bc = bc->GetValue(id_ncell);
1016 id_ncell = -1;
1017 };
1018 F.owner = pcIdxCell(id_cell);
1019 F.node.resize(4);
1020 F.node[0] = node_t(this, node, pts[0], 0);
1021 F.node[1] = node_t(this, node, pts[1], 0);
1022 F.node[2] = node_t(this, node, pts[5], 0);
1023 F.node[3] = node_t(this, node, pts[4], 0);
1024 if (id_ncell != -1) {
1025 F.neighbour = pcIdxCell(id_ncell);
1026 if (F.owner < F.neighbour) face_list.append(F);
1027 } else {
1028 F.neighbour = -1;
1029 face_list.append(F);
1030 };
1031
1032 // face 3
1033 id_ncell = c2c[id_cell][3];
1034 F.bc = 0;
1035 if (id_ncell == -1) EG_BUG;
1036 if (isSurface(id_ncell, grid)) {
1037 F.bc = bc->GetValue(id_ncell);
1038 id_ncell = -1;
1039 };
1040 F.owner = pcIdxCell(id_cell);
1041 F.node.resize(4);
1042 F.node[0] = node_t(this, node, pts[3], 0);
1043 F.node[1] = node_t(this, node, pts[7], 0);
1044 F.node[2] = node_t(this, node, pts[6], 0);
1045 F.node[3] = node_t(this, node, pts[2], 0);
1046 if (id_ncell != -1) {
1047 F.neighbour = pcIdxCell(id_ncell);
1048 if (F.owner < F.neighbour) face_list.append(F);
1049 } else {
1050 F.neighbour = -1;
1051 face_list.append(F);
1052 };
1053
1054 // face 4
1055 id_ncell = c2c[id_cell][4];
1056 F.bc = 0;
1057 if (id_ncell == -1) EG_BUG;
1058 if (isSurface(id_ncell, grid)) {
1059 F.bc = bc->GetValue(id_ncell);
1060 id_ncell = -1;
1061 };
1062 F.owner = pcIdxCell(id_cell);
1063 F.node.resize(4);
1064 F.node[0] = node_t(this, node, pts[0], 0);
1065 F.node[1] = node_t(this, node, pts[4], 0);
1066 F.node[2] = node_t(this, node, pts[7], 0);
1067 F.node[3] = node_t(this, node, pts[3], 0);
1068 if (id_ncell != -1) {
1069 F.neighbour = pcIdxCell(id_ncell);
1070 if (F.owner < F.neighbour) face_list.append(F);
1071 } else {
1072 F.neighbour = -1;
1073 face_list.append(F);
1074 };
1075
1076 // face 5
1077 id_ncell = c2c[id_cell][5];
1078 F.bc = 0;
1079 if (id_ncell == -1) EG_BUG;
1080 if (isSurface(id_ncell, grid)) {
1081 F.bc = bc->GetValue(id_ncell);
1082 id_ncell = -1;
1083 };
1084 F.owner = pcIdxCell(id_cell);
1085 F.node.resize(4);
1086 F.node[0] = node_t(this, node, pts[1], 0);
1087 F.node[1] = node_t(this, node, pts[2], 0);
1088 F.node[2] = node_t(this, node, pts[6], 0);
1089 F.node[3] = node_t(this, node, pts[5], 0);
1090 if (id_ncell != -1) {
1091 F.neighbour = pcIdxCell(id_ncell);
1092 if (F.owner < F.neighbour) face_list.append(F);
1093 } else {
1094 F.neighbour = -1;
1095 face_list.append(F);
1096 };
1097
1098 };
1099 };
1100 };
1101
1102 void PolyMesh::sortFaces()
1103 {
1104 bcs.clear();
1105 foreach (face_t F, face_list) {
1106 if (!bcs.contains(F.bc)) {
1107 bcs.append(F.bc);
1108 };
1109 };
1110 qSort(bcs.begin(),bcs.end());
1111 QList<face_t> ordered_face_list;
1112 foreach (int bc, bcs) {
1113 foreach (face_t F, face_list) {
1114 if (F.bc == bc) {
1115 ordered_face_list.append(F);
1116 };
1117 };
1118 };
1119 faces.resize(ordered_face_list.size());
1120 qCopy(ordered_face_list.begin(), ordered_face_list.end(), faces.begin());
1121 };
1122
1123 void PolyMesh::pass1()
1124 {
1125 cout << "creating polyhedral grid" << endl;
1126 cell2pc.fill(-1,grid->GetNumberOfCells());
1127 node2pc.fill(-1,grid->GetNumberOfPoints());
1128 face_list.clear();
1129 id_newpc = 0;
1130 getAllCells(cells, grid);
1131 createCellToCell(cells, c2c, grid);
1132 cout << "pass 1 for tetras" << endl;
1133 pass1Tetras();
1134 cout << "pass 1 for prisms" << endl;
1135 pass1Prisms();
1136 cout << "pass 1 for hexas" << endl;
1137 pass1Hexas();
1138 cout << "sorting faces" << endl;
1139 sortFaces();
1140 };
1141
1142 void PolyMesh::createNodes()
1143 {
1144 QHash<node_t,int> nodemap;
1145 {
1146 int id_node = 0;
1147 foreach (face_t F, faces) {
1148 foreach (node_t N, F.node) {
1149 if (!nodemap.contains(N)) {
1150 nodemap.insert(N, id_node);
1151 ++id_node;
1152 };
1153 };
1154 };
1155 };
1156 nodes.resize(nodemap.size());
1157 QVector<bool> nodeset(nodemap.size(), false);
1158 for (int i_face = 0; i_face < faces.size(); ++i_face) {
1159 face_t F = faces[i_face];
1160 for (int i_node = 0; i_node < F.node.size(); ++i_node) {
1161 node_t N = F.node[i_node];
1162 int i = nodemap.value(N);
1163 if (!nodeset[i]) {
1164 nodeset[i] = true;
1165 nodes[i] = vec3_t(N.x,N.y,N.z);
1166 };
1167 N.type = node;
1168 N.idx = i;
1169 N.subidx = 0;
1170 F.node[i_node] = N;
1171 };
1172 faces[i_face] = F;
1173 };
1174 };
1175
1176 void PolyMesh::computeNodes()
1177 {
1178 for (int i_face = 0; i_face < faces.size(); ++i_face) {
1179 face_t F = faces[i_face];
1180 for (int i_node = 0; i_node < F.node.size(); ++i_node) {
1181 node_t N = F.node[i_node];
1182 if (N.type == node) {
1183 vec3_t x;
1184 grid->GetPoint(N.idx, x.data());
1185 N.x = x[0];
1186 N.y = x[1];
1187 N.z = x[2];
1188 };
1189 if (N.type == edge) {
1190 vtkIdType N_pts, *pts;
1191 vtkIdType type_cell = grid->GetCellType(N.idx);
1192 grid->GetCellPoints(N.idx, N_pts, pts);
1193 vec3_t x1, x2;
1194 double w1=1, w2=1;
1195 if (type_cell == VTK_TETRA) {
1196 if (N.subidx == 0) {
1197 grid->GetPoint(pts[0], x1.data());
1198 grid->GetPoint(pts[1], x2.data());
1199 w1 = weight[pts[0]];
1200 w2 = weight[pts[1]];
1201 };
1202 if (N.subidx == 1) {
1203 grid->GetPoint(pts[0], x1.data());
1204 grid->GetPoint(pts[2], x2.data());
1205 w1 = weight[pts[0]];
1206 w2 = weight[pts[2]];
1207 };
1208 if (N.subidx == 2) {
1209 grid->GetPoint(pts[0], x1.data());
1210 grid->GetPoint(pts[3], x2.data());
1211 w1 = weight[pts[0]];
1212 w2 = weight[pts[3]];
1213 };
1214 if (N.subidx == 3) {
1215 grid->GetPoint(pts[1], x1.data());
1216 grid->GetPoint(pts[2], x2.data());
1217 w1 = weight[pts[1]];
1218 w2 = weight[pts[2]];
1219 };
1220 if (N.subidx == 4) {
1221 grid->GetPoint(pts[1], x1.data());
1222 grid->GetPoint(pts[3], x2.data());
1223 w1 = weight[pts[1]];
1224 w2 = weight[pts[3]];
1225 };
1226 if (N.subidx == 5) {
1227 grid->GetPoint(pts[2], x1.data());
1228 grid->GetPoint(pts[3], x2.data());
1229 w1 = weight[pts[2]];
1230 w2 = weight[pts[3]];
1231 };
1232 };
1233 if (type_cell == VTK_WEDGE) {
1234 if (N.subidx == 0) {
1235 grid->GetPoint(pts[0], x1.data());
1236 grid->GetPoint(pts[1], x2.data());
1237 w1 = weight[pts[0]];
1238 w2 = weight[pts[1]];
1239 };
1240 if (N.subidx == 1) {
1241 grid->GetPoint(pts[0], x1.data());
1242 grid->GetPoint(pts[2], x2.data());
1243 w1 = weight[pts[0]];
1244 w2 = weight[pts[2]];
1245 };
1246 if (N.subidx == 2) {
1247 grid->GetPoint(pts[0], x1.data());
1248 grid->GetPoint(pts[3], x2.data());
1249 w1 = weight[pts[0]];
1250 w2 = weight[pts[3]];
1251 };
1252 if (N.subidx == 3) {
1253 grid->GetPoint(pts[1], x1.data());
1254 grid->GetPoint(pts[2], x2.data());
1255 w1 = weight[pts[1]];
1256 w2 = weight[pts[2]];
1257 };
1258 if (N.subidx == 4) {
1259 grid->GetPoint(pts[1], x1.data());
1260 grid->GetPoint(pts[4], x2.data());
1261 w1 = weight[pts[1]];
1262 w2 = weight[pts[4]];
1263 };
1264 if (N.subidx == 5) {
1265 grid->GetPoint(pts[2], x1.data());
1266 grid->GetPoint(pts[5], x2.data());
1267 w1 = weight[pts[2]];
1268 w2 = weight[pts[5]];
1269 };
1270 if (N.subidx == 6) {
1271 grid->GetPoint(pts[3], x1.data());
1272 grid->GetPoint(pts[4], x2.data());
1273 w1 = weight[pts[3]];
1274 w2 = weight[pts[4]];
1275 };
1276 if (N.subidx == 7) {
1277 grid->GetPoint(pts[3], x1.data());
1278 grid->GetPoint(pts[5], x2.data());
1279 w1 = weight[pts[3]];
1280 w2 = weight[pts[5]];
1281 };
1282 if (N.subidx == 8) {
1283 grid->GetPoint(pts[4], x1.data());
1284 grid->GetPoint(pts[5], x2.data());
1285 w1 = weight[pts[4]];
1286 w2 = weight[pts[5]];
1287 };
1288 };
1289 {
1290 vec3_t x = 1.0/(w1+w2)*(w1*x1+w2*x2);
1291 N.x = x[0];
1292 N.y = x[1];
1293 N.z = x[2];
1294 };
1295 };
1296 if (N.type == face) {
1297 vtkIdType N_pts, *pts;
1298 vtkIdType type_cell = grid->GetCellType(N.idx);
1299 grid->GetCellPoints(N.idx, N_pts, pts);
1300 if (type_cell == VTK_TETRA) {
1301 vec3_t x1, x2, x3;
1302 double w1=1, w2=1, w3=1;
1303 if (N.subidx == 0) {
1304 grid->GetPoint(pts[0], x1.data());
1305 grid->GetPoint(pts[1], x2.data());
1306 grid->GetPoint(pts[2], x3.data());
1307 w1 = faceW(weight[pts[0]]);
1308 w2 = faceW(weight[pts[1]]);
1309 w3 = faceW(weight[pts[2]]);
1310 };
1311 if (N.subidx == 1) {
1312 grid->GetPoint(pts[0], x1.data());
1313 grid->GetPoint(pts[1], x2.data());
1314 grid->GetPoint(pts[3], x3.data());
1315 w1 = faceW(weight[pts[0]]);
1316 w2 = faceW(weight[pts[1]]);
1317 w3 = faceW(weight[pts[3]]);
1318 };
1319 if (N.subidx == 2) {
1320 grid->GetPoint(pts[0], x1.data());
1321 grid->GetPoint(pts[2], x2.data());
1322 grid->GetPoint(pts[3], x3.data());
1323 w1 = faceW(weight[pts[0]]);
1324 w2 = faceW(weight[pts[2]]);
1325 w3 = faceW(weight[pts[3]]);
1326 };
1327 if (N.subidx == 3) {
1328 grid->GetPoint(pts[1], x1.data());
1329 grid->GetPoint(pts[2], x2.data());
1330 grid->GetPoint(pts[3], x3.data());
1331 w1 = faceW(weight[pts[1]]);
1332 w2 = faceW(weight[pts[2]]);
1333 w3 = faceW(weight[pts[3]]);
1334 };
1335 {
1336 vec3_t x = 1.0/(w1+w2+w3)*(w1*x1+w2*x2+w3*x3);
1337 N.x = x[0];
1338 N.y = x[1];
1339 N.z = x[2];
1340 };
1341 };
1342 if (type_cell == VTK_WEDGE) {
1343 vec3_t x1, x2, x3, x4;
1344 double w1, w2, w3;
1345 if (N.subidx == 0) {
1346 grid->GetPoint(pts[0], x1.data());
1347 grid->GetPoint(pts[1], x2.data());
1348 grid->GetPoint(pts[2], x3.data());
1349 w1 = faceW(weight[pts[0]]);
1350 w2 = faceW(weight[pts[1]]);
1351 w3 = faceW(weight[pts[2]]);
1352 {
1353 vec3_t x = 1.0/(w1+w2+w3)*(w1*x1+w2*x2+w3*x3);
1354 N.x = x[0];
1355 N.y = x[1];
1356 N.z = x[2];
1357 };
1358 };
1359 if (N.subidx == 1) {
1360 grid->GetPoint(pts[3], x1.data());
1361 grid->GetPoint(pts[4], x2.data());
1362 grid->GetPoint(pts[5], x3.data());
1363 w1 = faceW(weight[pts[3]]);
1364 w2 = faceW(weight[pts[4]]);
1365 w3 = faceW(weight[pts[5]]);
1366 {
1367 vec3_t x = 1.0/(w1+w2+w3)*(w1*x1+w2*x2+w3*x3);
1368 N.x = x[0];
1369 N.y = x[1];
1370 N.z = x[2];
1371 };
1372 };
1373 if (N.subidx == 2) {
1374 grid->GetPoint(pts[0], x1.data());
1375 grid->GetPoint(pts[1], x2.data());
1376 grid->GetPoint(pts[3], x3.data());
1377 grid->GetPoint(pts[4], x4.data());
1378 {
1379 vec3_t x = 0.25*(x1+x2+x3+x4);
1380 N.x = x[0];
1381 N.y = x[1];
1382 N.z = x[2];
1383 };
1384 };
1385 if (N.subidx == 3) {
1386 grid->GetPoint(pts[1], x1.data());
1387 grid->GetPoint(pts[2], x2.data());
1388 grid->GetPoint(pts[4], x3.data());
1389 grid->GetPoint(pts[5], x4.data());
1390 {
1391 vec3_t x = 0.25*(x1+x2+x3+x4);
1392 N.x = x[0];
1393 N.y = x[1];
1394 N.z = x[2];
1395 };
1396 };
1397 if (N.subidx == 4) {
1398 grid->GetPoint(pts[0], x1.data());
1399 grid->GetPoint(pts[2], x2.data());
1400 grid->GetPoint(pts[3], x3.data());
1401 grid->GetPoint(pts[5], x4.data());
1402 {
1403 vec3_t x = 0.25*(x1+x2+x3+x4);
1404 N.x = x[0];
1405 N.y = x[1];
1406 N.z = x[2];
1407 };
1408 };
1409 };
1410 };
1411 if (N.type == cell) {
1412 vtkIdType N_pts, *pts;
1413 grid->GetCellPoints(N.idx, N_pts, pts);
1414 N.x = 0;
1415 N.y = 0;
1416 N.z = 0;
1417 for (int j = 0; j < N_pts; ++j) {
1418 vec3_t x;
1419 grid->GetPoint(pts[j], x.data());
1420 N.x += x[0];
1421 N.y += x[1];
1422 N.z += x[2];
1423 };
1424 N.x *= 1.0/N_pts;
1425 N.y *= 1.0/N_pts;
1426 N.z *= 1.0/N_pts;
1427 };
1428 F.node[i_node] = N;
1429 };
1430 faces[i_face] = F;
1431 };
1432 };
1433
1434 PolyMesh::face_t PolyMesh::combineFaces(QList<face_t> faces)
1435 {
1436 if (faces.size() < 2) EG_BUG;
1437 QVector<face_t> src_faces(faces.size());
1438 QVector<face_t> dst_faces(faces.size());
1439 qCopy(faces.begin(), faces.end(), src_faces.begin());
1440 QVector<bool> src_used(src_faces.size(),false);
1441 int N = src_faces.size()-1;
1442 face_t F_cur = src_faces[0];
1443 int i_cur = 0;
1444 src_used[0] = true;
1445 //int dir = 0;
1446 node_t centre, before, after, match_node;
1447 if (dbg) {
1448 cout << "\n\ncombining faces:\n";
1449 for (int i = 0; i < src_faces.size(); ++i) cout << src_faces[i] << endl;
1450 };
1451 {
1452 int i;
1453 for (i = 0; i < F_cur.node.size(); ++i) {
1454 if (F_cur[i].type == node) {
1455 centre = F_cur[i];
1456 break;
1457 };
1458 };
1459 if (i == F_cur.node.size()) {
1460 for (i = 0; i < F_cur.node.size(); ++i) {
1461 if (F_cur[i].type == edge) {
1462 centre = F_cur[i];
1463 break;
1464 };
1465 };
1466 if (i == F_cur.node.size()) EG_BUG;
1467 };
1468 before = F_cur[i-1];
1469 after = F_cur[i+1];
1470 };
1471 if (dbg) cout << "check" << endl;
1472 {
1473 int i, j = 0;
1474 match_node = after;
1475 int loops = 0;
1476 do {
1477 for (i = 0; i < src_faces.size(); ++i) {
1478 if (!src_used[i]) {
1479 for (j = 0; j < src_faces[i].node.size(); ++j) {
1480 if (src_faces[i][j] == centre) {
1481 if (src_faces[i][j-1] == match_node) break;
1482 };
1483 };
1484 if (j != src_faces[i].node.size()) break;
1485 };
1486 };
1487 if (i != src_faces.size()) {
1488 ++loops;
1489 F_cur = src_faces[i];
1490 i_cur = i;
1491 before = F_cur[j-1];
1492 after = F_cur[j+1];
1493 match_node = after;
1494 src_used[i] = true;
1495 };
1496 } while ((i != src_faces.size()) && (loops < src_faces.size()));
1497 };
1498 if (dbg) cout << "check" << endl;
1499 dst_faces[0] = F_cur;
1500 int i_dst = 1;
1501 for (int i = 0; i < src_faces.size(); ++i) {
1502 if (i == i_cur) src_used[i] = true;
1503 else src_used[i] = false;
1504 };
1505 match_node = before;
1506 //cout << "centre: " << centre << endl;
1507 //cout << "match_node: " << match_node << endl;
1508 //cout << "face 0: " << F_cur << endl;
1509 if (dbg) cout << "check" << endl;
1510 while (N > 0) {
1511 int i, j = 0;
1512 for (i = 0; i < src_faces.size(); ++i) {
1513 if (!src_used[i]) {
1514 for (j = 0; j < src_faces[i].node.size(); ++j) {
1515 if (src_faces[i][j] == centre) {
1516 if (src_faces[i][j+1] == match_node) break;
1517 };
1518 };
1519 if (j != src_faces[i].node.size()) break;
1520 };
1521 };
1522 if (i == src_faces.size()) EG_BUG;
1523 dst_faces[i_dst] = src_faces[i];
1524 src_used[i] = true;
1525 //cout << "face " << i_dst << ": " << dst_faces[i_dst] << endl;
1526 ++i_dst;
1527 --N;
1528 F_cur = src_faces[i];
1529 match_node = F_cur[j-1];
1530 };
1531 if (dbg) cout << "check" << endl;
1532 QList<node_t> new_nodes;
1533 {
1534 int i;
1535 int k = 0;
1536 for (i = 0; i < dst_faces.size(); ++i) {
1537 int j;
1538 for (j = 0; j < dst_faces[i].node.size(); ++j) {
1539 if (dst_faces[i][j] == centre) break;
1540 };
1541 if (j == dst_faces[i].node.size()) EG_BUG;
1542 //cout << "j=" << j << endl;
1543 k = j + 1;
1544 while (!(dst_faces[i][k+1] == dst_faces[i][j])) {
1545 new_nodes.append(dst_faces[i][k]);
1546 //cout << "appending " << dst_faces[i][k] << endl;
1547 k += 1;
1548 };
1549 };
1550 if (!(dst_faces[i-1][k] == new_nodes.first())) {
1551 new_nodes.append(dst_faces[i-1][k]);
1552 new_nodes.append(centre);
1553 };
1554 };
1555 if (dbg) cout << "check" << endl;
1556 face_t F_new;
1557 F_new.owner = src_faces[0].owner;
1558 F_new.neighbour = src_faces[0].neighbour;
1559 F_new.bc = src_faces[0].bc;
1560 F_new.node.resize(new_nodes.size());
1561 /*
1562 if (dir < 0) {
1563 QList<node_t>::iterator iter = new_nodes.begin();
1564 int i = new_nodes.size()-1;
1565 while (iter != new_nodes.end()) {
1566 F_new.node[i] = *iter;
1567 ++iter;
1568 --i;
1569 };
1570 } else {
1571 qCopy(new_nodes.begin(), new_nodes.end(), F_new.node.begin());
1572 };
1573 */
1574 if (dbg) cout << "check" << endl;
1575 qCopy(new_nodes.begin(), new_nodes.end(), F_new.node.begin());
1576 //cout << "dir: " << dir << endl;
1577 //cout << "new face: " << F_new << endl;
1578 if (dbg) cout << "check" << endl;
1579 return F_new;
1580 };
1581
1582 void PolyMesh::pass2()
1583 {
1584 cout << "pass 2" << endl;
1585 QVector<vec3_t> nodex_save(grid->GetNumberOfPoints());
1586 QVector<bool> node_changed(grid->GetNumberOfPoints(),false);
1587 for (vtkIdType id_node = 0; id_node < grid->GetNumberOfPoints(); ++id_node) {
1588 grid->GetPoint(id_node, nodex_save[id_node].data());
1589 };
1590 cout << "computing node coordinates" << endl;
1591 computeNodes();
1592
1593 cout << "adjusting outer layer prisms" << endl;
1594 for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
1595 if (grid->GetCellType(id_cell) == VTK_WEDGE) {
1596 bool ok = true;
1597 if (grid->GetCellType(c2c[id_cell][2]) == VTK_HEXAHEDRON) ok = false;
1598 if (grid->GetCellType(c2c[id_cell][3]) == VTK_HEXAHEDRON) ok = false;
1599 if (grid->GetCellType(c2c[id_cell][4]) == VTK_HEXAHEDRON) ok = false;
1600 if (ok) {
1601 vtkIdType N_pts, *pts;
1602 grid->GetCellPoints(id_cell, N_pts, pts);
1603 vec3_t x[6];
1604 for (int i = 0; i < 6; ++i) grid->GetPoint(pts[i], x[i].data());
1605 if (grid->GetCellType(c2c[id_cell][0]) == VTK_TETRA) {
1606 if (!node_changed[pts[0]]) {
1607 x[0] = 0.5*(x[0]+x[3]);
1608 node_changed[pts[0]] = true;
1609 };
1610 if (!node_changed[pts[1]]) {
1611 x[1] = 0.5*(x[1]+x[4]);
1612 node_changed[pts[1]] = true;
1613 };
1614 if (!node_changed[pts[2]]) {
1615 x[2] = 0.5*(x[2]+x[5]);
1616 node_changed[pts[2]] = true;
1617 };
1618 };
1619 if (grid->GetCellType(c2c[id_cell][1]) == VTK_TETRA) {
1620 if (!node_changed[pts[3]]) {
1621 x[3] = 0.5*(x[0]+x[3]);
1622 node_changed[pts[3]] = true;
1623 };
1624 if (!node_changed[pts[4]]) {
1625 x[4] = 0.5*(x[1]+x[4]);
1626 node_changed[pts[4]] = true;
1627 };
1628 if (!node_changed[pts[5]]) {
1629 x[5] = 0.5*(x[2]+x[5]);
1630 node_changed[pts[5]] = true;
1631 };
1632 };
1633 for (int i = 0; i < 6; ++i) grid->GetPoints()->SetPoint(pts[i], x[i].data());
1634 };
1635 };
1636 };
1637 for (int i_face = 0; i_face < faces.size(); ++i_face) {
1638 face_t F = faces[i_face];
1639 for (int i_node = 0; i_node < F.node.size(); ++i_node) {
1640 node_t N = F.node[i_node];
1641 if (N.type == node) {
1642 vec3_t x;
1643 grid->GetPoint(N.idx, x.data());
1644 N.x = x[0];
1645 N.y = x[1];
1646 N.z = x[2];
1647 };
1648 F.node[i_node] = N;
1649 };
1650 faces[i_face] = F;
1651 };
1652 for (vtkIdType id_node = 0; id_node < grid->GetNumberOfPoints(); ++id_node) {
1653 grid->GetPoints()->SetPoint(id_node, nodex_save[id_node].data());
1654 };
1655
1656 if (faces.size() < 2) return;
1657 cout << "sorting faces" << endl;
1658 qStableSort(faces.begin(), faces.end());
1659 QList<face_t> new_faces;
1660 int i = 0;
1661 cout << "combining faces (" << faces.size() << ")" << endl;
1662 while (i < faces.size()) {
1663 //cout << i << endl;
1664 //if (i > 2230000) dbg = true;
1665 //else dbg = false;
1666 face_t F = faces[i];
1667 QList<face_t> combine_faces;
1668 while ((i < faces.size()) && (faces[i].owner == F.owner) && (faces[i].neighbour == F.neighbour)) {
1669 combine_faces.append(faces[i]);
1670 if (dbg) cout << i << endl;
1671 ++i;
1672 };
1673 bool combine = false;
1674 if (combine_faces.size() > 1) {
1675 QList<face_t>::iterator j = combine_faces.begin();
1676 int bc = j->bc;
1677 combine = true;
1678 double mincosa = 1.0;
1679 vec3_t n(0,0,0);
1680 while (j != combine_faces.end()) {
1681 if (j->bc != bc) combine = false;
1682 n += j->normal();
1683 ++j;
1684 };
1685 if (combine) {
1686 n.normalise();
1687 j = combine_faces.begin();
1688 while (j != combine_faces.end()) {
1689 vec3_t nj = j->normal();
1690 nj.normalise();
1691 double cosa = n*nj;
1692 mincosa = min(cosa,mincosa);
1693 ++j;
1694 };
1695 if ((mincosa < cos(45.0*M_PI/180)) && (bc != 0)) {
1696 combine = false;
1697 };
1698 };
1699 };
1700 if (combine) {
1701 new_faces.append(combineFaces(combine_faces));
1702 } else {
1703 foreach (face_t F, combine_faces) {
1704 new_faces.append(F);
1705 };
1706 };
1707 };
1708 cout << "sorting faces" << endl;
1709 faces.resize(new_faces.size());
1710 qCopy(new_faces.begin(), new_faces.end(), faces.begin());
1711 qStableSort(faces.begin(), faces.end());
1712 };
1713
1714 void PolyMesh::pass3()
1715 {
1716 /*
1717 QVector<vec3_t> n(grid->GetNumberOfPoints(),vec3_t(0,0,0));
1718 for (int i = 0; i < faces.size(); ++i) {
1719 for (int j = 0; j < faces[i].node.size(); ++j) {
1720 node_t N = faces[i][j];
1721 N.n = vec3_t(0,0,0);
1722 faces[i].node[j] = N;
1723 };
1724 };
1725 for (int i = 0; i < faces.size(); ++i) {
1726 vec3_t ni = faces[i].normal();
1727 for (int j = 0; j < faces[i].node.size(); ++j) {
1728 node_t N = faces[i][j];
1729 N.n += ni;
1730 faces[i].node[j] = N;
1731 };
1732 };
1733 for (int i = 0; i < faces.size(); ++i) {
1734 for (int j = 0; j < faces[i].node.size(); ++j) {
1735 node_t N = faces[i][j];
1736 if (N.type == node) {
1737 n[N.idx] += N.n;
1738 };
1739 };
1740 };
1741 for (int i = 0; i < faces.size(); ++i) {
1742 for (int j = 0; j < faces[i].node.size(); ++j) {
1743 node_t N = faces[i][j];
1744 if (N.type == node) {
1745 N.n = n[N.idx];
1746 };
1747 faces[i].node[j] = N;
1748 };
1749 };
1750 for (int i = 0; i < n.size(); ++i) n[i].normalise();
1751 double w_flat = 1.0;
1752 double w_sharp = 5.0;
1753 double w_max = w_flat;
1754 for (int i = 0; i < weight.size(); ++i) weight[i] = w_flat;
1755 for (int i = 0; i < faces.size(); ++i) {
1756 vec3_t ni = faces[i].normal();
1757 ni.normalise();
1758 for (int j = 0; j < faces[i].node.size(); ++j) {
1759 node_t N = faces[i][j];
1760 N.n.normalise();
1761 double f = min(1.0,(1.0-N.n*ni));
1762 N.w = w_flat + f*(w_sharp - w_flat);
1763 w_max = max(w_max,N.w);
1764 if (N.type == node) weight[N.idx] = N.w;
1765 faces[i].node[j] = N;
1766 };
1767 };
1768 cout << "maximal edge weighting : " << w_max << endl;
1769 computeNodes();
1770 */
1771 cout << "pass 3" << endl;
1772 createNodes();
1773 return;
1774 {
1775 QVector<bool> is_boundary(nodes.size(),false);
1776 QVector<bool> del_node(nodes.size(),false);
1777 QVector<int> face_count(nodes.size(),0);
1778 foreach (face_t F, faces) {
1779 foreach (node_t N, F.node) {
1780 if (N.type != node) EG_BUG;
1781 if (F.bc != 0) {
1782 is_boundary[N.idx] = true;
1783 };
1784 };
1785 };
1786 foreach (face_t F, faces) {
1787 foreach (node_t N, F.node) {
1788 if (is_boundary[N.idx]) {
1789 if (F.bc != 0) {
1790 ++face_count[N.idx];
1791 };
1792 } else {
1793 ++face_count[N.idx];
1794 };
1795 };
1796 };
1797 QVector<int> old2new(nodes.size());
1798 int N = 0;
1799 {
1800 int j = 0;
1801 for (int i = 0; i < nodes.size(); ++i) {
1802 old2new[i] = j;
1803 int fc = 2;
1804 if (face_count[i] <= fc) {
1805 del_node[i] = true;
1806 ++N;
1807 } else {
1808 ++j;
1809 };
1810 };
1811 };
1812 cout << "removing " << N << " nodes" << endl;
1813 for (int i = 0; i < faces.size(); ++i) {
1814 face_t F = faces[i];
1815 QList<node_t> new_node;
1816 foreach (node_t N, F.node) {
1817 if (!del_node[N.idx]) {
1818 N.idx = old2new[N.idx];
1819 new_node.append(N);
1820 };
1821 };
1822 F.node.resize(new_node.size());
1823 qCopy(new_node.begin(), new_node.end(), F.node.begin());
1824 faces[i] = F;
1825 };
1826 QVector<vec3_t> new_nodes(nodes.size()-N);
1827 for (int i = 0; i < nodes.size(); ++i) {
1828 if (!del_node[i]) {
1829 new_nodes[old2new[i]] = nodes[i];
1830 };
1831 };
1832 nodes.resize(new_nodes.size());
1833 qCopy(new_nodes.begin(), new_nodes.end(), nodes.begin());
1834 };
1835 cout << "done" << endl;
1836 };
Clone of the official engrid repository