2 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
4 // + This file is part of enGrid. +
6 // + Copyright 2008,2009 Oliver Gloth +
8 // + enGrid is free software: you can redistribute it and/or modify +
9 // + it under the terms of the GNU General Public License as published by +
10 // + the Free Software Foundation, either version 3 of the License, or +
11 // + (at your option) any later version. +
13 // + enGrid is distributed in the hope that it will be useful, +
14 // + but WITHOUT ANY WARRANTY; without even the implied warranty of +
15 // + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +
16 // + GNU General Public License for more details. +
18 // + You should have received a copy of the GNU General Public License +
19 // + along with enGrid. If not, see <http://www.gnu.org/licenses/>. +
21 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
25 uint
qHash(PolyMesh::node_t N
)
28 if (N
.type
== PolyMesh::node
) {
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
];
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
;
46 bool PolyMesh::node_t::operator==(const node_t N
) const
50 poly
->getEdge(idx
,subidx
,nds
);
53 poly
->getFace(idx
,subidx
,nds
);
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
)) {
59 QList
<vtkIdType
> nds1
, nds2
;
61 poly
->getEdge(idx
,subidx
,nds1
);
62 poly
->getEdge(N
.idx
,N
.subidx
,nds2
);
65 poly
->getFace(idx
,subidx
,nds1
);
66 poly
->getFace(N
.idx
,N
.subidx
,nds2
);
68 foreach (vtkIdType p
, nds1
) {
69 if (!nds2
.contains(p
)) {
79 ostream
& operator<<(ostream
&s
, PolyMesh::node_t N
)
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
) {
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
] << ',';
95 if (N
.type
== PolyMesh::edge
) {
96 N
.poly
->getEdge(N
.idx
,N
.subidx
,nds
);
98 N
.poly
->getFace(N
.idx
,N
.subidx
,nds
);
100 QList
<vtkIdType
>::iterator i
= nds
.begin();
103 while (i
!= nds
.end()) {
113 ostream
& operator<<(ostream
&s
, PolyMesh::face_t F
)
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
];
122 void PolyMesh::face_t::checkOrientation()
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];
130 swap(owner
,neighbour
);
134 bool PolyMesh::face_t::operator<(const face_t
&F
) const
139 } else if (bc
== F
.bc
) {
140 if (owner
< F
.owner
) {
142 } else if (owner
== F
.owner
) {
143 if (neighbour
< F
.neighbour
) {
152 vec3_t
PolyMesh::face_t::normal()
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
);
159 x
*= 1.0/node
.size();
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
));
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
));
172 void PolyMesh::getFace(vtkIdType idx
, int subidx
, QList
<vtkIdType
> &nodes
)
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
) {
180 nodes
.append(pts
[2]);
181 nodes
.append(pts
[1]);
182 nodes
.append(pts
[0]);
185 nodes
.append(pts
[1]);
186 nodes
.append(pts
[3]);
187 nodes
.append(pts
[0]);
190 nodes
.append(pts
[3]);
191 nodes
.append(pts
[2]);
192 nodes
.append(pts
[0]);
195 nodes
.append(pts
[2]);
196 nodes
.append(pts
[3]);
197 nodes
.append(pts
[1]);
200 if (type_cell
== VTK_WEDGE
) {
202 nodes
.append(pts
[0]);
203 nodes
.append(pts
[1]);
204 nodes
.append(pts
[2]);
207 nodes
.append(pts
[3]);
208 nodes
.append(pts
[5]);
209 nodes
.append(pts
[4]);
212 nodes
.append(pts
[3]);
213 nodes
.append(pts
[4]);
214 nodes
.append(pts
[1]);
215 nodes
.append(pts
[0]);
218 nodes
.append(pts
[1]);
219 nodes
.append(pts
[4]);
220 nodes
.append(pts
[2]);
221 nodes
.append(pts
[5]);
224 nodes
.append(pts
[0]);
225 nodes
.append(pts
[2]);
226 nodes
.append(pts
[5]);
227 nodes
.append(pts
[3]);
230 if (type_cell
== VTK_HEXAHEDRON
) {
232 nodes
.append(pts
[0]);
233 nodes
.append(pts
[3]);
234 nodes
.append(pts
[2]);
235 nodes
.append(pts
[1]);
238 nodes
.append(pts
[4]);
239 nodes
.append(pts
[5]);
240 nodes
.append(pts
[6]);
241 nodes
.append(pts
[7]);
244 nodes
.append(pts
[0]);
245 nodes
.append(pts
[1]);
246 nodes
.append(pts
[5]);
247 nodes
.append(pts
[4]);
250 nodes
.append(pts
[3]);
251 nodes
.append(pts
[7]);
252 nodes
.append(pts
[6]);
253 nodes
.append(pts
[2]);
256 nodes
.append(pts
[0]);
257 nodes
.append(pts
[4]);
258 nodes
.append(pts
[7]);
259 nodes
.append(pts
[3]);
262 nodes
.append(pts
[1]);
263 nodes
.append(pts
[2]);
264 nodes
.append(pts
[6]);
265 nodes
.append(pts
[5]);
270 void PolyMesh::getEdge(vtkIdType idx
, int subidx
, QList
<vtkIdType
> &nodes
)
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
) {
278 nodes
.append(pts
[0]);
279 nodes
.append(pts
[1]);
282 nodes
.append(pts
[0]);
283 nodes
.append(pts
[2]);
286 nodes
.append(pts
[0]);
287 nodes
.append(pts
[3]);
290 nodes
.append(pts
[1]);
291 nodes
.append(pts
[2]);
294 nodes
.append(pts
[1]);
295 nodes
.append(pts
[3]);
298 nodes
.append(pts
[2]);
299 nodes
.append(pts
[3]);
302 if (type_cell
== VTK_WEDGE
) {
304 nodes
.append(pts
[0]);
305 nodes
.append(pts
[1]);
308 nodes
.append(pts
[0]);
309 nodes
.append(pts
[2]);
312 nodes
.append(pts
[0]);
313 nodes
.append(pts
[3]);
316 nodes
.append(pts
[1]);
317 nodes
.append(pts
[2]);
320 nodes
.append(pts
[1]);
321 nodes
.append(pts
[4]);
324 nodes
.append(pts
[2]);
325 nodes
.append(pts
[5]);
328 nodes
.append(pts
[3]);
329 nodes
.append(pts
[4]);
332 nodes
.append(pts
[3]);
333 nodes
.append(pts
[5]);
336 nodes
.append(pts
[4]);
337 nodes
.append(pts
[5]);
340 if (type_cell
== VTK_HEXAHEDRON
) {
342 nodes
.append(pts
[0]);
343 nodes
.append(pts
[1]);
346 nodes
.append(pts
[0]);
347 nodes
.append(pts
[3]);
350 nodes
.append(pts
[0]);
351 nodes
.append(pts
[4]);
354 nodes
.append(pts
[1]);
355 nodes
.append(pts
[2]);
358 nodes
.append(pts
[1]);
359 nodes
.append(pts
[5]);
362 nodes
.append(pts
[2]);
363 nodes
.append(pts
[3]);
366 nodes
.append(pts
[2]);
367 nodes
.append(pts
[6]);
370 nodes
.append(pts
[3]);
371 nodes
.append(pts
[7]);
374 nodes
.append(pts
[4]);
375 nodes
.append(pts
[5]);
378 nodes
.append(pts
[4]);
379 nodes
.append(pts
[7]);
382 nodes
.append(pts
[5]);
383 nodes
.append(pts
[6]);
386 nodes
.append(pts
[6]);
387 nodes
.append(pts
[7]);
392 PolyMesh::PolyMesh(vtkUnstructuredGrid
*a_grid
, bool dual_mesh
)
397 weight
.fill(1.0,grid
->GetNumberOfPoints());
403 int PolyMesh::pcIdxNode(vtkIdType id_node
)
405 if (node2pc
[id_node
] != -1) {
406 return node2pc
[id_node
];
408 node2pc
[id_node
] = id_newpc
;
410 return node2pc
[id_node
];
413 int PolyMesh::pcIdxCell(vtkIdType id_cell
)
415 if (cell2pc
[id_cell
] != -1) {
416 return cell2pc
[id_cell
];
418 cell2pc
[id_cell
] = id_newpc
;
420 return cell2pc
[id_cell
];
423 void PolyMesh::pass1Tetras()
425 EG_VTKDCC(vtkIntArray
, bc
, grid
, "cell_code");
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
);
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]);
444 F
.checkOrientation();
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]);
456 F
.checkOrientation();
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]);
468 F
.checkOrientation();
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]);
480 F
.checkOrientation();
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]);
492 F
.checkOrientation();
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]);
504 F
.checkOrientation();
509 vtkIdType id_bcell
= c2c
[id_cell
][0];
510 if (id_bcell
== -1) EG_BUG
;
511 if (isSurface(id_bcell
, grid
)) {
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]);
520 F
.bc
= bc
->GetValue(id_bcell
);
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]);
530 F
.bc
= bc
->GetValue(id_bcell
);
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]);
540 F
.bc
= bc
->GetValue(id_bcell
);
548 vtkIdType id_bcell
= c2c
[id_cell
][1];
549 if (id_bcell
== -1) EG_BUG
;
550 if (isSurface(id_bcell
, grid
)) {
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]);
559 F
.bc
= bc
->GetValue(id_bcell
);
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]);
569 F
.bc
= bc
->GetValue(id_bcell
);
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]);
579 F
.bc
= bc
->GetValue(id_bcell
);
587 vtkIdType id_bcell
= c2c
[id_cell
][2];
588 if (id_bcell
== -1) EG_BUG
;
589 if (isSurface(id_bcell
, grid
)) {
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]);
598 F
.bc
= bc
->GetValue(id_bcell
);
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]);
608 F
.bc
= bc
->GetValue(id_bcell
);
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]);
618 F
.bc
= bc
->GetValue(id_bcell
);
626 vtkIdType id_bcell
= c2c
[id_cell
][3];
627 if (id_bcell
== -1) EG_BUG
;
628 if (isSurface(id_bcell
, grid
)) {
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]);
637 F
.bc
= bc
->GetValue(id_bcell
);
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]);
647 F
.bc
= bc
->GetValue(id_bcell
);
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]);
657 F
.bc
= bc
->GetValue(id_bcell
);
668 void PolyMesh::pass1Prisms()
670 EG_VTKDCC(vtkIntArray
, bc
, grid
, "cell_code");
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
);
681 id_ncell
= c2c
[id_cell
][0];
683 if (id_ncell
== -1) EG_BUG
;
684 if (isSurface(id_ncell
, grid
)) {
685 F
.bc
= bc
->GetValue(id_ncell
);
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
);
694 if (grid
->GetCellType(id_ncell
) == VTK_WEDGE
) {
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
);
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]);
714 F
.checkOrientation();
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]);
725 F
.checkOrientation();
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]);
736 F
.checkOrientation();
742 id_ncell
= c2c
[id_cell
][1];
744 if (id_ncell
== -1) EG_BUG
;
745 if (isSurface(id_ncell
, grid
)) {
746 F
.bc
= bc
->GetValue(id_ncell
);
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
);
755 if (grid
->GetCellType(id_ncell
) == VTK_WEDGE
) {
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
);
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]);
775 F
.checkOrientation();
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]);
786 F
.checkOrientation();
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]);
797 F
.checkOrientation();
803 id_ncell
= c2c
[id_cell
][2];
805 if (id_ncell
== -1) EG_BUG
;
806 if (isSurface(id_ncell
, grid
)) {
807 F
.bc
= bc
->GetValue(id_ncell
);
811 F
.owner
= pcIdxCell(id_cell
);
812 if (!f0_split
&& !f1_split
) {
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);
819 if (f0_split
&& !f1_split
) {
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);
827 if (!f0_split
&& f1_split
) {
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);
835 if (f0_split
&& f1_split
) {
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);
844 if (id_ncell
!= -1) {
845 F
.neighbour
= pcIdxCell(id_ncell
);
846 if (F
.owner
< F
.neighbour
) face_list
.append(F
);
853 id_ncell
= c2c
[id_cell
][3];
855 if (id_ncell
== -1) EG_BUG
;
856 if (isSurface(id_ncell
, grid
)) {
857 F
.bc
= bc
->GetValue(id_ncell
);
860 F
.owner
= pcIdxCell(id_cell
);
861 if (!f0_split
&& !f1_split
) {
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);
868 if (f0_split
&& !f1_split
) {
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);
876 if (!f0_split
&& f1_split
) {
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);
884 if (f0_split
&& f1_split
) {
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);
894 if (id_ncell
!= -1) {
895 F
.neighbour
= pcIdxCell(id_ncell
);
896 if (F
.owner
< F
.neighbour
) face_list
.append(F
);
903 id_ncell
= c2c
[id_cell
][4];
905 if (id_ncell
== -1) EG_BUG
;
906 if (isSurface(id_ncell
, grid
)) {
907 F
.bc
= bc
->GetValue(id_ncell
);
910 F
.owner
= pcIdxCell(id_cell
);
911 if (!f0_split
&& !f1_split
) {
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
);
919 if (f0_split
&& !f1_split
) {
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);
927 if (!f0_split
&& f1_split
) {
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);
935 if (f0_split
&& f1_split
) {
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);
944 if (id_ncell
!= -1) {
945 F
.neighbour
= pcIdxCell(id_ncell
);
946 if (F
.owner
< F
.neighbour
) face_list
.append(F
);
956 void PolyMesh::pass1Hexas()
958 EG_VTKDCC(vtkIntArray
, bc
, grid
, "cell_code");
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
);
967 id_ncell
= c2c
[id_cell
][0];
969 if (id_ncell
== -1) EG_BUG
;
970 if (isSurface(id_ncell
, grid
)) {
971 F
.bc
= bc
->GetValue(id_ncell
);
974 F
.owner
= pcIdxCell(id_cell
);
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
);
989 id_ncell
= c2c
[id_cell
][1];
991 if (id_ncell
== -1) EG_BUG
;
992 if (isSurface(id_ncell
, grid
)) {
993 F
.bc
= bc
->GetValue(id_ncell
);
996 F
.owner
= pcIdxCell(id_cell
);
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
);
1007 face_list
.append(F
);
1011 id_ncell
= c2c
[id_cell
][2];
1013 if (id_ncell
== -1) EG_BUG
;
1014 if (isSurface(id_ncell
, grid
)) {
1015 F
.bc
= bc
->GetValue(id_ncell
);
1018 F
.owner
= pcIdxCell(id_cell
);
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
);
1029 face_list
.append(F
);
1033 id_ncell
= c2c
[id_cell
][3];
1035 if (id_ncell
== -1) EG_BUG
;
1036 if (isSurface(id_ncell
, grid
)) {
1037 F
.bc
= bc
->GetValue(id_ncell
);
1040 F
.owner
= pcIdxCell(id_cell
);
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
);
1051 face_list
.append(F
);
1055 id_ncell
= c2c
[id_cell
][4];
1057 if (id_ncell
== -1) EG_BUG
;
1058 if (isSurface(id_ncell
, grid
)) {
1059 F
.bc
= bc
->GetValue(id_ncell
);
1062 F
.owner
= pcIdxCell(id_cell
);
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
);
1073 face_list
.append(F
);
1077 id_ncell
= c2c
[id_cell
][5];
1079 if (id_ncell
== -1) EG_BUG
;
1080 if (isSurface(id_ncell
, grid
)) {
1081 F
.bc
= bc
->GetValue(id_ncell
);
1084 F
.owner
= pcIdxCell(id_cell
);
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
);
1095 face_list
.append(F
);
1102 void PolyMesh::sortFaces()
1105 foreach (face_t F
, face_list
) {
1106 if (!bcs
.contains(F
.bc
)) {
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
) {
1115 ordered_face_list
.append(F
);
1119 faces
.resize(ordered_face_list
.size());
1120 qCopy(ordered_face_list
.begin(), ordered_face_list
.end(), faces
.begin());
1123 void PolyMesh::pass1()
1125 cout
<< "creating polyhedral grid" << endl
;
1126 cell2pc
.fill(-1,grid
->GetNumberOfCells());
1127 node2pc
.fill(-1,grid
->GetNumberOfPoints());
1130 getAllCells(cells
, grid
);
1131 createCellToCell(cells
, c2c
, grid
);
1132 cout
<< "pass 1 for tetras" << endl
;
1134 cout
<< "pass 1 for prisms" << endl
;
1136 cout
<< "pass 1 for hexas" << endl
;
1138 cout
<< "sorting faces" << endl
;
1142 void PolyMesh::createNodes()
1144 QHash
<node_t
,int> nodemap
;
1147 foreach (face_t F
, faces
) {
1148 foreach (node_t N
, F
.node
) {
1149 if (!nodemap
.contains(N
)) {
1150 nodemap
.insert(N
, id_node
);
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
);
1165 nodes
[i
] = vec3_t(N
.x
,N
.y
,N
.z
);
1176 void PolyMesh::computeNodes()
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
) {
1184 grid
->GetPoint(N
.idx
, x
.data());
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
);
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]];
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]];
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]];
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]];
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]];
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]];
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]];
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]];
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]];
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]];
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]];
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]];
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]];
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]];
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]];
1290 vec3_t x
= 1.0/(w1
+w2
)*(w1
*x1
+w2
*x2
);
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
) {
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]]);
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]]);
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]]);
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]]);
1336 vec3_t x
= 1.0/(w1
+w2
+w3
)*(w1
*x1
+w2
*x2
+w3
*x3
);
1342 if (type_cell
== VTK_WEDGE
) {
1343 vec3_t x1
, x2
, x3
, x4
;
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]]);
1353 vec3_t x
= 1.0/(w1
+w2
+w3
)*(w1
*x1
+w2
*x2
+w3
*x3
);
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]]);
1367 vec3_t x
= 1.0/(w1
+w2
+w3
)*(w1
*x1
+w2
*x2
+w3
*x3
);
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());
1379 vec3_t x
= 0.25*(x1
+x2
+x3
+x4
);
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());
1391 vec3_t x
= 0.25*(x1
+x2
+x3
+x4
);
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());
1403 vec3_t x
= 0.25*(x1
+x2
+x3
+x4
);
1411 if (N
.type
== cell
) {
1412 vtkIdType N_pts
, *pts
;
1413 grid
->GetCellPoints(N
.idx
, N_pts
, pts
);
1417 for (int j
= 0; j
< N_pts
; ++j
) {
1419 grid
->GetPoint(pts
[j
], x
.data());
1434 PolyMesh::face_t
PolyMesh::combineFaces(QList
<face_t
> faces
)
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];
1446 node_t centre
, before
, after
, match_node
;
1448 cout
<< "\n\ncombining faces:\n";
1449 for (int i
= 0; i
< src_faces
.size(); ++i
) cout
<< src_faces
[i
] << endl
;
1453 for (i
= 0; i
< F_cur
.node
.size(); ++i
) {
1454 if (F_cur
[i
].type
== node
) {
1459 if (i
== F_cur
.node
.size()) {
1460 for (i
= 0; i
< F_cur
.node
.size(); ++i
) {
1461 if (F_cur
[i
].type
== edge
) {
1466 if (i
== F_cur
.node
.size()) EG_BUG
;
1468 before
= F_cur
[i
-1];
1471 if (dbg
) cout
<< "check" << endl
;
1477 for (i
= 0; i
< src_faces
.size(); ++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;
1484 if (j
!= src_faces
[i
].node
.size()) break;
1487 if (i
!= src_faces
.size()) {
1489 F_cur
= src_faces
[i
];
1491 before
= F_cur
[j
-1];
1496 } while ((i
!= src_faces
.size()) && (loops
< src_faces
.size()));
1498 if (dbg
) cout
<< "check" << endl
;
1499 dst_faces
[0] = F_cur
;
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;
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
;
1512 for (i
= 0; i
< src_faces
.size(); ++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;
1519 if (j
!= src_faces
[i
].node
.size()) break;
1522 if (i
== src_faces
.size()) EG_BUG
;
1523 dst_faces
[i_dst
] = src_faces
[i
];
1525 //cout << "face " << i_dst << ": " << dst_faces[i_dst] << endl;
1528 F_cur
= src_faces
[i
];
1529 match_node
= F_cur
[j
-1];
1531 if (dbg
) cout
<< "check" << endl
;
1532 QList
<node_t
> new_nodes
;
1536 for (i
= 0; i
< dst_faces
.size(); ++i
) {
1538 for (j
= 0; j
< dst_faces
[i
].node
.size(); ++j
) {
1539 if (dst_faces
[i
][j
] == centre
) break;
1541 if (j
== dst_faces
[i
].node
.size()) EG_BUG
;
1542 //cout << "j=" << j << endl;
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;
1550 if (!(dst_faces
[i
-1][k
] == new_nodes
.first())) {
1551 new_nodes
.append(dst_faces
[i
-1][k
]);
1552 new_nodes
.append(centre
);
1555 if (dbg
) cout
<< "check" << endl
;
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());
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;
1571 qCopy(new_nodes.begin(), new_nodes.end(), F_new.node.begin());
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
;
1582 void PolyMesh::pass2()
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());
1590 cout
<< "computing node coordinates" << endl
;
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
) {
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;
1601 vtkIdType N_pts
, *pts
;
1602 grid
->GetCellPoints(id_cell
, N_pts
, pts
);
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;
1610 if (!node_changed
[pts
[1]]) {
1611 x
[1] = 0.5*(x
[1]+x
[4]);
1612 node_changed
[pts
[1]] = true;
1614 if (!node_changed
[pts
[2]]) {
1615 x
[2] = 0.5*(x
[2]+x
[5]);
1616 node_changed
[pts
[2]] = true;
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;
1624 if (!node_changed
[pts
[4]]) {
1625 x
[4] = 0.5*(x
[1]+x
[4]);
1626 node_changed
[pts
[4]] = true;
1628 if (!node_changed
[pts
[5]]) {
1629 x
[5] = 0.5*(x
[2]+x
[5]);
1630 node_changed
[pts
[5]] = true;
1633 for (int i
= 0; i
< 6; ++i
) grid
->GetPoints()->SetPoint(pts
[i
], x
[i
].data());
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
) {
1643 grid
->GetPoint(N
.idx
, x
.data());
1652 for (vtkIdType id_node
= 0; id_node
< grid
->GetNumberOfPoints(); ++id_node
) {
1653 grid
->GetPoints()->SetPoint(id_node
, nodex_save
[id_node
].data());
1656 if (faces
.size() < 2) return;
1657 cout
<< "sorting faces" << endl
;
1658 qStableSort(faces
.begin(), faces
.end());
1659 QList
<face_t
> new_faces
;
1661 cout
<< "combining faces (" << faces
.size() << ")" << endl
;
1662 while (i
< faces
.size()) {
1663 //cout << i << endl;
1664 //if (i > 2230000) dbg = true;
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
;
1673 bool combine
= false;
1674 if (combine_faces
.size() > 1) {
1675 QList
<face_t
>::iterator j
= combine_faces
.begin();
1678 double mincosa
= 1.0;
1680 while (j
!= combine_faces
.end()) {
1681 if (j
->bc
!= bc
) combine
= false;
1687 j
= combine_faces
.begin();
1688 while (j
!= combine_faces
.end()) {
1689 vec3_t nj
= j
->normal();
1692 mincosa
= min(cosa
,mincosa
);
1695 if ((mincosa
< cos(45.0*M_PI
/180)) && (bc
!= 0)) {
1701 new_faces
.append(combineFaces(combine_faces
));
1703 foreach (face_t F
, combine_faces
) {
1704 new_faces
.append(F
);
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());
1714 void PolyMesh::pass3()
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;
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];
1730 faces[i].node[j] = N;
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) {
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) {
1747 faces[i].node[j] = N;
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();
1758 for (int j = 0; j < faces[i].node.size(); ++j) {
1759 node_t N = faces[i][j];
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;
1768 cout << "maximal edge weighting : " << w_max << endl;
1771 cout
<< "pass 3" << endl
;
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
;
1782 is_boundary
[N
.idx
] = true;
1786 foreach (face_t F
, faces
) {
1787 foreach (node_t N
, F
.node
) {
1788 if (is_boundary
[N
.idx
]) {
1790 ++face_count
[N
.idx
];
1793 ++face_count
[N
.idx
];
1797 QVector
<int> old2new(nodes
.size());
1801 for (int i
= 0; i
< nodes
.size(); ++i
) {
1804 if (face_count
[i
] <= fc
) {
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
];
1822 F
.node
.resize(new_node
.size());
1823 qCopy(new_node
.begin(), new_node
.end(), F
.node
.begin());
1826 QVector
<vec3_t
> new_nodes(nodes
.size()-N
);
1827 for (int i
= 0; i
< nodes
.size(); ++i
) {
1829 new_nodes
[old2new
[i
]] = nodes
[i
];
1832 nodes
.resize(new_nodes
.size());
1833 qCopy(new_nodes
.begin(), new_nodes
.end(), nodes
.begin());
1835 cout
<< "done" << endl
;