1 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3 // + This file is part of enGrid. +
5 // + Copyright 2008-2014 enGits GmbH +
7 // + enGrid is free software: you can redistribute it and/or modify +
8 // + it under the terms of the GNU General Public License as published by +
9 // + the Free Software Foundation, either version 3 of the License, or +
10 // + (at your option) any later version. +
12 // + enGrid is distributed in the hope that it will be useful, +
13 // + but WITHOUT ANY WARRANTY; without even the implied warranty of +
14 // + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +
15 // + GNU General Public License for more details. +
17 // + You should have received a copy of the GNU General Public License +
18 // + along with enGrid. If not, see <http://www.gnu.org/licenses/>. +
20 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
23 #include "geometrytools.h"
24 #include "deletestraynodes.h"
25 #include "guimainwindow.h"
27 OctreeCell::OctreeCell()
29 for (int i
= 0; i
< 8; ++i
) {
33 for (int i
= 0; i
< 6; ++i
) {
40 int OctreeCell::findNode(int i_node
)
43 for (int i
= 0; i
< 8; ++i
) {
44 if (m_Node
[i
] == i_node
) {
52 int OctreeCell::getEdgeNode(Octree
* octree
, int n1
, int n2
, bool this_cell_only
)
54 vec3_t x1
= octree
->getNodePosition(m_Node
[n1
]);
55 vec3_t x2
= octree
->getNodePosition(m_Node
[n2
]);
56 vec3_t xm
= 0.5*(x1
+ x2
);
63 if (n2
== 1 || n2
== 2 || n2
== 4) {
64 edge_node
= octree
->m_Cells
[m_Child
[n1
]].m_Node
[n2
];
67 if (n2
== 3 || n2
== 5) {
68 edge_node
= octree
->m_Cells
[m_Child
[n1
]].m_Node
[n2
];
71 if (n2
== 3 || n2
== 6) {
72 edge_node
= octree
->m_Cells
[m_Child
[n1
]].m_Node
[n2
];
76 edge_node
= octree
->m_Cells
[m_Child
[n1
]].m_Node
[n2
];
79 if (n2
== 5 || n2
== 6) {
80 edge_node
= octree
->m_Cells
[m_Child
[n1
]].m_Node
[n2
];
84 edge_node
= octree
->m_Cells
[m_Child
[n1
]].m_Node
[n2
];
88 edge_node
= octree
->m_Cells
[m_Child
[n1
]].m_Node
[n2
];
91 } else if (!this_cell_only
) {
92 QSet
<int> cells1
= QSet
<int>::fromList(octree
->m_Node2Cell
[m_Node
[n1
]]);
93 QSet
<int> cells2
= QSet
<int>::fromList(octree
->m_Node2Cell
[m_Node
[n2
]]);
94 QSet
<int> cells
= cells1
.intersect(cells2
);
95 foreach (int i_cell
, cells
) {
96 if (octree
->m_Cells
[i_cell
].m_Level
== m_Level
) {
97 int nn1
= octree
->m_Cells
[i_cell
].findNode(m_Node
[n1
]);
98 int nn2
= octree
->m_Cells
[i_cell
].findNode(m_Node
[n2
]);
99 edge_node
= octree
->m_Cells
[i_cell
].getEdgeNode(octree
, nn1
, nn2
, true);
100 if (edge_node
!= -1) {
109 void OctreeCell::getFaceNodes(int i
, Octree
* octree
, QVector
<int>& face_nodes
, bool reverse
)
111 //face_nodes.resize(4);
115 nodes
.push_back(m_Node
[0]);
116 edge_node
= getEdgeNode(octree
, 0,4); if (edge_node
!= -1) nodes
.push_back(edge_node
);
117 nodes
.push_back(m_Node
[4]);
118 edge_node
= getEdgeNode(octree
, 4,6); if (edge_node
!= -1) nodes
.push_back(edge_node
);
119 nodes
.push_back(m_Node
[6]);
120 edge_node
= getEdgeNode(octree
, 6,2); if (edge_node
!= -1) nodes
.push_back(edge_node
);
121 nodes
.push_back(m_Node
[2]);
122 edge_node
= getEdgeNode(octree
, 2,0); if (edge_node
!= -1) nodes
.push_back(edge_node
);
124 nodes
.push_back(m_Node
[1]);
125 edge_node
= getEdgeNode(octree
, 1,3); if (edge_node
!= -1) nodes
.push_back(edge_node
);
126 nodes
.push_back(m_Node
[3]);
127 edge_node
= getEdgeNode(octree
, 3,7); if (edge_node
!= -1) nodes
.push_back(edge_node
);
128 nodes
.push_back(m_Node
[7]);
129 edge_node
= getEdgeNode(octree
, 7,5); if (edge_node
!= -1) nodes
.push_back(edge_node
);
130 nodes
.push_back(m_Node
[5]);
131 edge_node
= getEdgeNode(octree
, 5,1); if (edge_node
!= -1) nodes
.push_back(edge_node
);
133 nodes
.push_back(m_Node
[0]);
134 edge_node
= getEdgeNode(octree
, 0,1); if (edge_node
!= -1) nodes
.push_back(edge_node
);
135 nodes
.push_back(m_Node
[1]);
136 edge_node
= getEdgeNode(octree
, 1,5); if (edge_node
!= -1) nodes
.push_back(edge_node
);
137 nodes
.push_back(m_Node
[5]);
138 edge_node
= getEdgeNode(octree
, 5,4); if (edge_node
!= -1) nodes
.push_back(edge_node
);
139 nodes
.push_back(m_Node
[4]);
140 edge_node
= getEdgeNode(octree
, 4,0); if (edge_node
!= -1) nodes
.push_back(edge_node
);
142 nodes
.push_back(m_Node
[3]);
143 edge_node
= getEdgeNode(octree
, 3,2); if (edge_node
!= -1) nodes
.push_back(edge_node
);
144 nodes
.push_back(m_Node
[2]);
145 edge_node
= getEdgeNode(octree
, 2,6); if (edge_node
!= -1) nodes
.push_back(edge_node
);
146 nodes
.push_back(m_Node
[6]);
147 edge_node
= getEdgeNode(octree
, 6,7); if (edge_node
!= -1) nodes
.push_back(edge_node
);
148 nodes
.push_back(m_Node
[7]);
149 edge_node
= getEdgeNode(octree
, 7,3); if (edge_node
!= -1) nodes
.push_back(edge_node
);
151 nodes
.push_back(m_Node
[0]);
152 edge_node
= getEdgeNode(octree
, 0,2); if (edge_node
!= -1) nodes
.push_back(edge_node
);
153 nodes
.push_back(m_Node
[2]);
154 edge_node
= getEdgeNode(octree
, 2,3); if (edge_node
!= -1) nodes
.push_back(edge_node
);
155 nodes
.push_back(m_Node
[3]);
156 edge_node
= getEdgeNode(octree
, 3,1); if (edge_node
!= -1) nodes
.push_back(edge_node
);
157 nodes
.push_back(m_Node
[1]);
158 edge_node
= getEdgeNode(octree
, 1,0); if (edge_node
!= -1) nodes
.push_back(edge_node
);
160 nodes
.push_back(m_Node
[4]);
161 edge_node
= getEdgeNode(octree
, 4,5); if (edge_node
!= -1) nodes
.push_back(edge_node
);
162 nodes
.push_back(m_Node
[5]);
163 edge_node
= getEdgeNode(octree
, 5,7); if (edge_node
!= -1) nodes
.push_back(edge_node
);
164 nodes
.push_back(m_Node
[7]);
165 edge_node
= getEdgeNode(octree
, 7,6); if (edge_node
!= -1) nodes
.push_back(edge_node
);
166 nodes
.push_back(m_Node
[6]);
167 edge_node
= getEdgeNode(octree
, 6,4); if (edge_node
!= -1) nodes
.push_back(edge_node
);
170 face_nodes
.resize(nodes
.size());
172 int i
= face_nodes
.size() - 1;
173 foreach(int node
, nodes
) {
174 face_nodes
[i
] = node
;
179 foreach(int node
, nodes
) {
180 face_nodes
[i
] = node
;
186 void OctreeCell::getFaceNodes(int i
, Octree
* octree
, QVector
<QVector
<int> >& face_nodes
, bool reverse
)
189 face_nodes
.resize(4);
191 octree
->m_Cells
[m_Child
[0]].getFaceNodes(0, octree
, face_nodes
[0], reverse
);
192 octree
->m_Cells
[m_Child
[2]].getFaceNodes(0, octree
, face_nodes
[1], reverse
);
193 octree
->m_Cells
[m_Child
[4]].getFaceNodes(0, octree
, face_nodes
[2], reverse
);
194 octree
->m_Cells
[m_Child
[6]].getFaceNodes(0, octree
, face_nodes
[3], reverse
);
196 octree
->m_Cells
[m_Child
[1]].getFaceNodes(1, octree
, face_nodes
[0], reverse
);
197 octree
->m_Cells
[m_Child
[3]].getFaceNodes(1, octree
, face_nodes
[1], reverse
);
198 octree
->m_Cells
[m_Child
[5]].getFaceNodes(1, octree
, face_nodes
[2], reverse
);
199 octree
->m_Cells
[m_Child
[7]].getFaceNodes(1, octree
, face_nodes
[3], reverse
);
201 octree
->m_Cells
[m_Child
[0]].getFaceNodes(2, octree
, face_nodes
[0], reverse
);
202 octree
->m_Cells
[m_Child
[1]].getFaceNodes(2, octree
, face_nodes
[1], reverse
);
203 octree
->m_Cells
[m_Child
[4]].getFaceNodes(2, octree
, face_nodes
[2], reverse
);
204 octree
->m_Cells
[m_Child
[5]].getFaceNodes(2, octree
, face_nodes
[3], reverse
);
206 octree
->m_Cells
[m_Child
[2]].getFaceNodes(3, octree
, face_nodes
[0], reverse
);
207 octree
->m_Cells
[m_Child
[3]].getFaceNodes(3, octree
, face_nodes
[1], reverse
);
208 octree
->m_Cells
[m_Child
[6]].getFaceNodes(3, octree
, face_nodes
[2], reverse
);
209 octree
->m_Cells
[m_Child
[7]].getFaceNodes(3, octree
, face_nodes
[3], reverse
);
211 octree
->m_Cells
[m_Child
[0]].getFaceNodes(4, octree
, face_nodes
[0], reverse
);
212 octree
->m_Cells
[m_Child
[1]].getFaceNodes(4, octree
, face_nodes
[1], reverse
);
213 octree
->m_Cells
[m_Child
[2]].getFaceNodes(4, octree
, face_nodes
[2], reverse
);
214 octree
->m_Cells
[m_Child
[3]].getFaceNodes(4, octree
, face_nodes
[3], reverse
);
216 octree
->m_Cells
[m_Child
[4]].getFaceNodes(5, octree
, face_nodes
[0], reverse
);
217 octree
->m_Cells
[m_Child
[5]].getFaceNodes(5, octree
, face_nodes
[1], reverse
);
218 octree
->m_Cells
[m_Child
[6]].getFaceNodes(5, octree
, face_nodes
[2], reverse
);
219 octree
->m_Cells
[m_Child
[7]].getFaceNodes(5, octree
, face_nodes
[3], reverse
);
222 face_nodes
.resize(1);
223 getFaceNodes(i
, octree
, face_nodes
[0], reverse
);
227 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
231 m_Origin
= vec3_t(0,0,0);
232 setBase(vec3_t(1,0,0), vec3_t(0,1,0), vec3_t(0,0,1));
235 m_ToRefine
.fill(false, 1);
236 setBounds(vec3_t(0,0,0), vec3_t(1,1,1));
237 setSmoothTransitionOn();
238 setMaxCells(1000000);
241 void Octree::markToRefine(int cell
)
243 if (!hasChildren(cell
)) {
244 m_ToRefine
[cell
] = true;
248 void Octree::markAllToRefine()
250 for (int cell
= 0; cell
< m_Cells
.size(); ++cell
) {
255 void Octree::mergeNodes_identifyDuplicates()
257 m_SameNodes
.resize(m_Nodes
.size());
258 for (int i_nodes
= 0; i_nodes
< m_Nodes
.size(); ++i_nodes
) {
259 m_SameNodes
[i_nodes
] = i_nodes
;
264 for (int i_cell
= 0; i_cell
< m_Cells
.size(); ++i_cell
) {
265 max_level
= max(max_level
, m_Cells
[i_cell
].m_Level
);
269 for (int i_node
= 0; i_node
< getNumNodes(); ++i_node
) {
270 QSet
<int> canditates
;
272 foreach (int i_cell
, m_Node2Cell
[i_node
]) {
273 tol
= min(tol
, 0.01*min(getDx(i_cell
), min(getDy(i_cell
), getDz(i_cell
))));
274 for (int i
= 0; i
< 8; ++i
) {
275 int j_node
= m_Cells
[i_cell
].m_Node
[i
];
276 foreach (int j_cell
, m_Node2Cell
[j_node
]) {
277 for (int j
= 0; j
< 8; ++j
) {
278 canditates
.insert(m_Cells
[j_cell
].m_Node
[j
]);
286 foreach (int j_node
, canditates
) {
287 if (i_node
!= j_node
) {
288 double distance
= (m_Nodes
[i_node
].m_Position
- m_Nodes
[j_node
].m_Position
).abs();
289 if (distance
< tol
) {
290 if (i_node
> j_node
) {
291 m_SameNodes
[i_node
] = j_node
;
293 m_SameNodes
[j_node
] = i_node
;
300 QVector
<bool> is_dup(m_Nodes
.size(), false);
301 for (int i_nodes
= 0; i_nodes
< m_Nodes
.size(); ++i_nodes
) {
303 int i_same_node
= i_nodes
;
304 while (i_same_node
!= m_SameNodes
[i_same_node
]) {
305 i_same_node
= m_SameNodes
[i_same_node
];
307 m_SameNodes
[i_nodes
] = i_same_node
;
308 if (m_SameNodes
[i_nodes
] != i_nodes
) {
309 is_dup
[m_SameNodes
[i_nodes
]] = true;
311 if (m_SameNodes
[m_SameNodes
[i_nodes
]] != m_SameNodes
[i_nodes
]) {
316 for (int i_nodes
= 0; i_nodes
< m_Nodes
.size(); ++i_nodes
) {
317 if (is_dup
[i_nodes
]) {
323 void Octree::mergeNodes_compactNodes()
325 QVector
<int> offset(m_Nodes
.size());
328 for (int i
= 0; i
< m_Nodes
.size(); ++i
) {
329 if (m_SameNodes
[i
] != i
) {
330 max_dist
= max(max_dist
, (m_Nodes
[i
].m_Position
- m_Nodes
[m_SameNodes
[i
]].m_Position
).abs());
331 if (m_SameNodes
[i
] > i
) {
336 offset
[i
] = last_offset
;
339 for (int i
= 0; i
< m_Nodes
.size(); ++i
) {
340 if (m_SameNodes
[i
] != i
) {
341 m_SameNodes
[i
] -= offset
[m_SameNodes
[i
]];
343 m_SameNodes
[i
] -= offset
[i
];
347 for (int i
= 0; i
< m_Nodes
.size(); ++i
) {
348 if (m_SameNodes
[i
] > i
) {
353 QVector
<int> copy_nodes(m_Nodes
.size() - last_offset
);
354 for (int i
= 0; i
< m_Nodes
.size(); ++i
) {
355 copy_nodes
[m_SameNodes
[i
]] = i
;
358 for (int i
= 0; i
< m_Nodes
.size() - last_offset
; ++i
) {
359 m_Nodes
[i
] = m_Nodes
[copy_nodes
[i
]];
362 m_Nodes
.remove(m_Nodes
.size() - last_offset
, last_offset
);
366 void Octree::mergeNodes_updateCells()
368 for (int i_cells
= 0; i_cells
< m_Cells
.size(); ++i_cells
) {
369 for (int j
= 0; j
< 8; ++j
) {
370 m_Cells
[i_cells
].m_Node
[j
] = m_SameNodes
[m_Cells
[i_cells
].m_Node
[j
]];
375 void Octree::checkNeighbours()
379 for (int i_cells
= 0; i_cells
< m_Cells
.size(); ++i_cells
) {
380 for (int i_faces
= 0; i_faces
< 6; ++i_faces
) {
381 int i_neigh_cells
= m_Cells
[i_cells
].m_Neighbour
[i_faces
];
382 if (i_neigh_cells
!= -1) {
383 int i_back_faces
= i_faces
;
384 if (i_faces
% 2 == 0) {
389 if (m_Cells
[i_cells
].m_Level
== m_Cells
[i_neigh_cells
].m_Level
) {
390 if (m_Cells
[i_neigh_cells
].m_Neighbour
[i_back_faces
] != i_cells
) {
391 cout
<< "neighbour error: " << i_cells
<< ',' << i_neigh_cells
<< ',' << i_faces
<< ',' << i_back_faces
<< ','
392 << getCellCentre(i_cells
) << ',' << getCellCentre(i_neigh_cells
) << endl
;
397 cout
<< "no neighbour: " << i_cells
<< ',' << i_faces
<< ',' << getFaceCentre(i_cells
, i_faces
) << endl
;
402 cout
<< Nerr
<< " errors and " << Nno
<< " faces without neighbour" << endl
;
405 void Octree::buildNode2Cell()
408 m_Node2Cell
.fill(QList
<int>(), getNumNodes());
409 for (int i_cell
= 0; i_cell
< getNumCells(); ++i_cell
) {
410 for (int j
= 0; j
< 8; ++j
) {
411 m_Node2Cell
[getNode(i_cell
, j
)].append(i_cell
);
414 for (int i_node
= 0; i_node
< getNumNodes(); ++i_node
) {
416 foreach (int i_cell
, m_Node2Cell
[i_node
]) {
417 max_level
= max(max_level
, m_Cells
[i_cell
].m_Level
);
422 void Octree::mergeNodes()
424 mergeNodes_identifyDuplicates();
425 mergeNodes_compactNodes();
426 mergeNodes_updateCells();
429 void Octree::setOrigin(vec3_t x0
)
434 void Octree::setBase(vec3_t g1
, vec3_t g2
, vec3_t g3
)
436 m_Base
.column(0, g1
);
437 m_Base
.column(1, g2
);
438 m_Base
.column(2, g3
);
439 m_InvBase
= m_Base
.inverse();
442 void Octree::setBounds(vec3_t corner1
, vec3_t corner2
, int num_i
, int num_j
, int num_k
)
444 vec3_t
g1(m_Base
[0][0], m_Base
[1][0], m_Base
[2][0]);
445 vec3_t
g2(m_Base
[0][1], m_Base
[1][1], m_Base
[2][1]);
446 vec3_t
g3(m_Base
[0][2], m_Base
[1][2], m_Base
[2][2]);
448 int num_cells
= num_i
*num_j
*num_k
;
449 int num_nodes
= (num_i
+ 1)*(num_j
+ 1)*(num_k
+ 1);
450 m_Cells
.resize(num_cells
);
451 m_Nodes
.resize(num_nodes
);
454 m_Dx
= fabs(corner1
[0] - corner2
[0])/num_i
;
455 m_Dy
= fabs(corner1
[1] - corner2
[1])/num_j
;
456 m_Dz
= fabs(corner1
[2] - corner2
[2])/num_k
;
458 for (int i
= 0; i
< num_i
+ 1; ++i
) {
459 for (int j
= 0; j
< num_j
+ 1; ++j
) {
460 for (int k
= 0; k
< num_k
+ 1; ++k
) {
461 int idx_node
= i
*(num_j
+ 1)*(num_k
+ 1) + j
*(num_k
+ 1) + k
;
462 m_Nodes
[idx_node
].setPosition(m_Corner1
+ i
*m_Dx
*g1
+ j
*m_Dy
*g2
+ k
*m_Dz
*g3
);
466 for (int i
= 0; i
< num_i
; ++i
) {
467 for (int j
= 0; j
< num_j
; ++j
) {
468 for (int k
= 0; k
< num_k
; ++k
) {
469 int idx_cell
= i
*num_j
*num_k
+ j
*num_k
+ k
;
470 int idx_node0
= (i
+ 0)*(num_j
+ 1)*(num_k
+ 1) + (j
+ 0)*(num_k
+ 1) + (k
+ 0);
471 int idx_node1
= (i
+ 1)*(num_j
+ 1)*(num_k
+ 1) + (j
+ 0)*(num_k
+ 1) + (k
+ 0);
472 int idx_node2
= (i
+ 0)*(num_j
+ 1)*(num_k
+ 1) + (j
+ 1)*(num_k
+ 1) + (k
+ 0);
473 int idx_node3
= (i
+ 1)*(num_j
+ 1)*(num_k
+ 1) + (j
+ 1)*(num_k
+ 1) + (k
+ 0);
474 int idx_node4
= (i
+ 0)*(num_j
+ 1)*(num_k
+ 1) + (j
+ 0)*(num_k
+ 1) + (k
+ 1);
475 int idx_node5
= (i
+ 1)*(num_j
+ 1)*(num_k
+ 1) + (j
+ 0)*(num_k
+ 1) + (k
+ 1);
476 int idx_node6
= (i
+ 0)*(num_j
+ 1)*(num_k
+ 1) + (j
+ 1)*(num_k
+ 1) + (k
+ 1);
477 int idx_node7
= (i
+ 1)*(num_j
+ 1)*(num_k
+ 1) + (j
+ 1)*(num_k
+ 1) + (k
+ 1);
478 m_Cells
[idx_cell
].m_Node
[0] = idx_node0
;
479 m_Cells
[idx_cell
].m_Node
[1] = idx_node1
;
480 m_Cells
[idx_cell
].m_Node
[2] = idx_node2
;
481 m_Cells
[idx_cell
].m_Node
[3] = idx_node3
;
482 m_Cells
[idx_cell
].m_Node
[4] = idx_node4
;
483 m_Cells
[idx_cell
].m_Node
[5] = idx_node5
;
484 m_Cells
[idx_cell
].m_Node
[6] = idx_node6
;
485 m_Cells
[idx_cell
].m_Node
[7] = idx_node7
;
491 int Octree::refineAll()
496 for (int i_cells
= 0; i_cells
< m_Cells
.size(); ++i_cells
) {
497 if (m_Cells
[i_cells
].m_Child
[0] != -1) {
498 m_ToRefine
[i_cells
] = false;
505 for (int i_cells
= 0; i_cells
< m_Cells
.size(); ++i_cells
) {
506 if (m_ToRefine
[i_cells
]) {
508 if (m_SmoothTransition
) {
509 for (int i
= 0; i
< 8; ++i
) {
510 foreach (int neigh
, m_Node2Cell
[m_Cells
[i_cells
].m_Node
[i
]]) {
511 if ((m_Cells
[neigh
].m_Level
< m_Cells
[i_cells
].m_Level
) && !m_ToRefine
[neigh
] && !hasChildren(neigh
)) {
526 if (N2
+ 8*N1
> m_MaxCells
) {
528 QString msg
= "maximal number of cells exceeded\n";
529 num
.setNum(N2
+ 8*N1
);
530 msg
+= num
+= " requested and ";
531 num
.setNum(m_MaxCells
);
532 msg
+= num
+ " allowed";
535 m_Cells
.insert(N2
, 8*N1
, OctreeCell());
537 int new_node
= m_Nodes
.size();
538 m_Nodes
.insert(m_Nodes
.size(), N1
*19, OctreeNode());
540 for (int i_cells
= 0; i_cells
< N1
; ++i_cells
) {
541 if (m_ToRefine
[i_cells
]) {
542 if (hasChildren(i_cells
)) {
546 nn
[0] = m_Cells
[i_cells
].m_Node
[0];
547 nn
[1] = m_Cells
[i_cells
].m_Node
[1];
548 nn
[2] = m_Cells
[i_cells
].m_Node
[2];
549 nn
[3] = m_Cells
[i_cells
].m_Node
[3];
550 nn
[4] = m_Cells
[i_cells
].m_Node
[4];
551 nn
[5] = m_Cells
[i_cells
].m_Node
[5];
552 nn
[6] = m_Cells
[i_cells
].m_Node
[6];
553 nn
[7] = m_Cells
[i_cells
].m_Node
[7];
557 m_Nodes
[new_node
++].m_Position
= 0.5*(m_Nodes
[m_Cells
[i_cells
].m_Node
[0]].m_Position
+ m_Nodes
[m_Cells
[i_cells
].m_Node
[2]].m_Position
);
559 m_Nodes
[new_node
++].m_Position
= 0.5*(m_Nodes
[m_Cells
[i_cells
].m_Node
[1]].m_Position
+ m_Nodes
[m_Cells
[i_cells
].m_Node
[3]].m_Position
);
561 m_Nodes
[new_node
++].m_Position
= 0.5*(m_Nodes
[m_Cells
[i_cells
].m_Node
[0]].m_Position
+ m_Nodes
[m_Cells
[i_cells
].m_Node
[1]].m_Position
);
563 m_Nodes
[new_node
++].m_Position
= 0.5*(m_Nodes
[m_Cells
[i_cells
].m_Node
[2]].m_Position
+ m_Nodes
[m_Cells
[i_cells
].m_Node
[3]].m_Position
);
565 m_Nodes
[new_node
++].m_Position
= 0.5*(m_Nodes
[m_Cells
[i_cells
].m_Node
[4]].m_Position
+ m_Nodes
[m_Cells
[i_cells
].m_Node
[6]].m_Position
);
567 m_Nodes
[new_node
++].m_Position
= 0.5*(m_Nodes
[m_Cells
[i_cells
].m_Node
[5]].m_Position
+ m_Nodes
[m_Cells
[i_cells
].m_Node
[7]].m_Position
);
569 m_Nodes
[new_node
++].m_Position
= 0.5*(m_Nodes
[m_Cells
[i_cells
].m_Node
[4]].m_Position
+ m_Nodes
[m_Cells
[i_cells
].m_Node
[5]].m_Position
);
571 m_Nodes
[new_node
++].m_Position
= 0.5*(m_Nodes
[m_Cells
[i_cells
].m_Node
[6]].m_Position
+ m_Nodes
[m_Cells
[i_cells
].m_Node
[7]].m_Position
);
573 m_Nodes
[new_node
++].m_Position
= 0.5*(m_Nodes
[m_Cells
[i_cells
].m_Node
[0]].m_Position
+ m_Nodes
[m_Cells
[i_cells
].m_Node
[4]].m_Position
);
575 m_Nodes
[new_node
++].m_Position
= 0.5*(m_Nodes
[m_Cells
[i_cells
].m_Node
[1]].m_Position
+ m_Nodes
[m_Cells
[i_cells
].m_Node
[5]].m_Position
);
577 m_Nodes
[new_node
++].m_Position
= 0.5*(m_Nodes
[m_Cells
[i_cells
].m_Node
[2]].m_Position
+ m_Nodes
[m_Cells
[i_cells
].m_Node
[6]].m_Position
);
579 m_Nodes
[new_node
++].m_Position
= 0.5*(m_Nodes
[m_Cells
[i_cells
].m_Node
[3]].m_Position
+ m_Nodes
[m_Cells
[i_cells
].m_Node
[7]].m_Position
);
582 m_Nodes
[new_node
++].m_Position
= 0.25*( m_Nodes
[ne
[0]].m_Position
+ m_Nodes
[ne
[4]].m_Position
583 + m_Nodes
[ne
[8]].m_Position
+ m_Nodes
[ne
[10]].m_Position
);
585 m_Nodes
[new_node
++].m_Position
= 0.25*( m_Nodes
[ne
[1]].m_Position
+ m_Nodes
[ne
[5]].m_Position
586 + m_Nodes
[ne
[9]].m_Position
+ m_Nodes
[ne
[11]].m_Position
);
588 m_Nodes
[new_node
++].m_Position
= 0.25*( m_Nodes
[ne
[2]].m_Position
+ m_Nodes
[ne
[6]].m_Position
589 + m_Nodes
[ne
[8]].m_Position
+ m_Nodes
[ne
[9]].m_Position
);
591 m_Nodes
[new_node
++].m_Position
= 0.25*( m_Nodes
[ne
[3]].m_Position
+ m_Nodes
[ne
[7]].m_Position
592 + m_Nodes
[ne
[10]].m_Position
+ m_Nodes
[ne
[11]].m_Position
);
594 m_Nodes
[new_node
++].m_Position
= 0.25*( m_Nodes
[ne
[0]].m_Position
+ m_Nodes
[ne
[1]].m_Position
595 + m_Nodes
[ne
[2]].m_Position
+ m_Nodes
[ne
[3]].m_Position
);
597 m_Nodes
[new_node
++].m_Position
= 0.25*( m_Nodes
[ne
[4]].m_Position
+ m_Nodes
[ne
[5]].m_Position
598 + m_Nodes
[ne
[6]].m_Position
+ m_Nodes
[ne
[7]].m_Position
);
600 m_Nodes
[new_node
++].m_Position
= 1.0/6.0*( m_Nodes
[nf
[0]].m_Position
+ m_Nodes
[nf
[1]].m_Position
+ m_Nodes
[nf
[2]].m_Position
601 + m_Nodes
[nf
[3]].m_Position
+ m_Nodes
[nf
[4]].m_Position
+ m_Nodes
[nf
[5]].m_Position
);
603 for (int child
= 0; child
< 8; ++child
) {
604 m_Cells
[i_cells
].m_Child
[child
] = N2
;
605 m_Cells
[N2
].m_Parent
= i_cells
;
606 m_Cells
[N2
].m_Level
= m_Cells
[i_cells
].m_Level
+ 1;
611 m_Cells
[m_Cells
[i_cells
].m_Child
[0]].m_Node
[0] = nn
[0];
612 m_Cells
[m_Cells
[i_cells
].m_Child
[0]].m_Node
[1] = ne
[2];
613 m_Cells
[m_Cells
[i_cells
].m_Child
[0]].m_Node
[2] = ne
[0];
614 m_Cells
[m_Cells
[i_cells
].m_Child
[0]].m_Node
[3] = nf
[4];
615 m_Cells
[m_Cells
[i_cells
].m_Child
[0]].m_Node
[4] = ne
[8];
616 m_Cells
[m_Cells
[i_cells
].m_Child
[0]].m_Node
[5] = nf
[2];
617 m_Cells
[m_Cells
[i_cells
].m_Child
[0]].m_Node
[6] = nf
[0];
618 m_Cells
[m_Cells
[i_cells
].m_Child
[0]].m_Node
[7] = nv
;
620 m_Cells
[m_Cells
[i_cells
].m_Child
[1]].m_Node
[0] = ne
[2];
621 m_Cells
[m_Cells
[i_cells
].m_Child
[1]].m_Node
[1] = nn
[1];
622 m_Cells
[m_Cells
[i_cells
].m_Child
[1]].m_Node
[2] = nf
[4];
623 m_Cells
[m_Cells
[i_cells
].m_Child
[1]].m_Node
[3] = ne
[1];
624 m_Cells
[m_Cells
[i_cells
].m_Child
[1]].m_Node
[4] = nf
[2];
625 m_Cells
[m_Cells
[i_cells
].m_Child
[1]].m_Node
[5] = ne
[9];
626 m_Cells
[m_Cells
[i_cells
].m_Child
[1]].m_Node
[6] = nv
;
627 m_Cells
[m_Cells
[i_cells
].m_Child
[1]].m_Node
[7] = nf
[1];
629 m_Cells
[m_Cells
[i_cells
].m_Child
[2]].m_Node
[0] = ne
[0];
630 m_Cells
[m_Cells
[i_cells
].m_Child
[2]].m_Node
[1] = nf
[4];
631 m_Cells
[m_Cells
[i_cells
].m_Child
[2]].m_Node
[2] = nn
[2];
632 m_Cells
[m_Cells
[i_cells
].m_Child
[2]].m_Node
[3] = ne
[3];
633 m_Cells
[m_Cells
[i_cells
].m_Child
[2]].m_Node
[4] = nf
[0];
634 m_Cells
[m_Cells
[i_cells
].m_Child
[2]].m_Node
[5] = nv
;
635 m_Cells
[m_Cells
[i_cells
].m_Child
[2]].m_Node
[6] = ne
[10];
636 m_Cells
[m_Cells
[i_cells
].m_Child
[2]].m_Node
[7] = nf
[3];
638 m_Cells
[m_Cells
[i_cells
].m_Child
[3]].m_Node
[0] = nf
[4];
639 m_Cells
[m_Cells
[i_cells
].m_Child
[3]].m_Node
[1] = ne
[1];
640 m_Cells
[m_Cells
[i_cells
].m_Child
[3]].m_Node
[2] = ne
[3];
641 m_Cells
[m_Cells
[i_cells
].m_Child
[3]].m_Node
[3] = nn
[3];
642 m_Cells
[m_Cells
[i_cells
].m_Child
[3]].m_Node
[4] = nv
;
643 m_Cells
[m_Cells
[i_cells
].m_Child
[3]].m_Node
[5] = nf
[1];
644 m_Cells
[m_Cells
[i_cells
].m_Child
[3]].m_Node
[6] = nf
[3];
645 m_Cells
[m_Cells
[i_cells
].m_Child
[3]].m_Node
[7] = ne
[11];
647 m_Cells
[m_Cells
[i_cells
].m_Child
[4]].m_Node
[0] = ne
[8];
648 m_Cells
[m_Cells
[i_cells
].m_Child
[4]].m_Node
[1] = nf
[2];
649 m_Cells
[m_Cells
[i_cells
].m_Child
[4]].m_Node
[2] = nf
[0];
650 m_Cells
[m_Cells
[i_cells
].m_Child
[4]].m_Node
[3] = nv
;
651 m_Cells
[m_Cells
[i_cells
].m_Child
[4]].m_Node
[4] = nn
[4];
652 m_Cells
[m_Cells
[i_cells
].m_Child
[4]].m_Node
[5] = ne
[6];
653 m_Cells
[m_Cells
[i_cells
].m_Child
[4]].m_Node
[6] = ne
[4];
654 m_Cells
[m_Cells
[i_cells
].m_Child
[4]].m_Node
[7] = nf
[5];
656 m_Cells
[m_Cells
[i_cells
].m_Child
[5]].m_Node
[0] = nf
[2];
657 m_Cells
[m_Cells
[i_cells
].m_Child
[5]].m_Node
[1] = ne
[9];
658 m_Cells
[m_Cells
[i_cells
].m_Child
[5]].m_Node
[2] = nv
;
659 m_Cells
[m_Cells
[i_cells
].m_Child
[5]].m_Node
[3] = nf
[1];
660 m_Cells
[m_Cells
[i_cells
].m_Child
[5]].m_Node
[4] = ne
[6];
661 m_Cells
[m_Cells
[i_cells
].m_Child
[5]].m_Node
[5] = nn
[5];
662 m_Cells
[m_Cells
[i_cells
].m_Child
[5]].m_Node
[6] = nf
[5];
663 m_Cells
[m_Cells
[i_cells
].m_Child
[5]].m_Node
[7] = ne
[5];
665 m_Cells
[m_Cells
[i_cells
].m_Child
[6]].m_Node
[0] = nf
[0];
666 m_Cells
[m_Cells
[i_cells
].m_Child
[6]].m_Node
[1] = nv
;
667 m_Cells
[m_Cells
[i_cells
].m_Child
[6]].m_Node
[2] = ne
[10];
668 m_Cells
[m_Cells
[i_cells
].m_Child
[6]].m_Node
[3] = nf
[3];
669 m_Cells
[m_Cells
[i_cells
].m_Child
[6]].m_Node
[4] = ne
[4];
670 m_Cells
[m_Cells
[i_cells
].m_Child
[6]].m_Node
[5] = nf
[5];
671 m_Cells
[m_Cells
[i_cells
].m_Child
[6]].m_Node
[6] = nn
[6];
672 m_Cells
[m_Cells
[i_cells
].m_Child
[6]].m_Node
[7] = ne
[7];
674 m_Cells
[m_Cells
[i_cells
].m_Child
[7]].m_Node
[0] = nv
;
675 m_Cells
[m_Cells
[i_cells
].m_Child
[7]].m_Node
[1] = nf
[1];
676 m_Cells
[m_Cells
[i_cells
].m_Child
[7]].m_Node
[2] = nf
[3];
677 m_Cells
[m_Cells
[i_cells
].m_Child
[7]].m_Node
[3] = ne
[11];
678 m_Cells
[m_Cells
[i_cells
].m_Child
[7]].m_Node
[4] = nf
[5];
679 m_Cells
[m_Cells
[i_cells
].m_Child
[7]].m_Node
[5] = ne
[5];
680 m_Cells
[m_Cells
[i_cells
].m_Child
[7]].m_Node
[6] = ne
[7];
681 m_Cells
[m_Cells
[i_cells
].m_Child
[7]].m_Node
[7] = nn
[7];
684 m_Cells
[m_Cells
[i_cells
].m_Child
[0]].m_Neighbour
[1] = m_Cells
[i_cells
].m_Child
[1];
685 m_Cells
[m_Cells
[i_cells
].m_Child
[0]].m_Neighbour
[3] = m_Cells
[i_cells
].m_Child
[2];
686 m_Cells
[m_Cells
[i_cells
].m_Child
[0]].m_Neighbour
[5] = m_Cells
[i_cells
].m_Child
[4];
688 m_Cells
[m_Cells
[i_cells
].m_Child
[1]].m_Neighbour
[0] = m_Cells
[i_cells
].m_Child
[0];
689 m_Cells
[m_Cells
[i_cells
].m_Child
[1]].m_Neighbour
[3] = m_Cells
[i_cells
].m_Child
[3];
690 m_Cells
[m_Cells
[i_cells
].m_Child
[1]].m_Neighbour
[5] = m_Cells
[i_cells
].m_Child
[5];
692 m_Cells
[m_Cells
[i_cells
].m_Child
[2]].m_Neighbour
[1] = m_Cells
[i_cells
].m_Child
[3];
693 m_Cells
[m_Cells
[i_cells
].m_Child
[2]].m_Neighbour
[2] = m_Cells
[i_cells
].m_Child
[0];
694 m_Cells
[m_Cells
[i_cells
].m_Child
[2]].m_Neighbour
[5] = m_Cells
[i_cells
].m_Child
[6];
696 m_Cells
[m_Cells
[i_cells
].m_Child
[3]].m_Neighbour
[0] = m_Cells
[i_cells
].m_Child
[2];
697 m_Cells
[m_Cells
[i_cells
].m_Child
[3]].m_Neighbour
[2] = m_Cells
[i_cells
].m_Child
[1];
698 m_Cells
[m_Cells
[i_cells
].m_Child
[3]].m_Neighbour
[5] = m_Cells
[i_cells
].m_Child
[7];
700 m_Cells
[m_Cells
[i_cells
].m_Child
[4]].m_Neighbour
[1] = m_Cells
[i_cells
].m_Child
[5];
701 m_Cells
[m_Cells
[i_cells
].m_Child
[4]].m_Neighbour
[3] = m_Cells
[i_cells
].m_Child
[6];
702 m_Cells
[m_Cells
[i_cells
].m_Child
[4]].m_Neighbour
[4] = m_Cells
[i_cells
].m_Child
[0];
704 m_Cells
[m_Cells
[i_cells
].m_Child
[5]].m_Neighbour
[0] = m_Cells
[i_cells
].m_Child
[4];
705 m_Cells
[m_Cells
[i_cells
].m_Child
[5]].m_Neighbour
[3] = m_Cells
[i_cells
].m_Child
[7];
706 m_Cells
[m_Cells
[i_cells
].m_Child
[5]].m_Neighbour
[4] = m_Cells
[i_cells
].m_Child
[1];
708 m_Cells
[m_Cells
[i_cells
].m_Child
[6]].m_Neighbour
[1] = m_Cells
[i_cells
].m_Child
[7];
709 m_Cells
[m_Cells
[i_cells
].m_Child
[6]].m_Neighbour
[2] = m_Cells
[i_cells
].m_Child
[4];
710 m_Cells
[m_Cells
[i_cells
].m_Child
[6]].m_Neighbour
[4] = m_Cells
[i_cells
].m_Child
[2];
712 m_Cells
[m_Cells
[i_cells
].m_Child
[7]].m_Neighbour
[0] = m_Cells
[i_cells
].m_Child
[6];
713 m_Cells
[m_Cells
[i_cells
].m_Child
[7]].m_Neighbour
[2] = m_Cells
[i_cells
].m_Child
[5];
714 m_Cells
[m_Cells
[i_cells
].m_Child
[7]].m_Neighbour
[4] = m_Cells
[i_cells
].m_Child
[3];
718 for (int i_cells
= 0; i_cells
< N1
; ++i_cells
) {
719 if (m_Cells
[i_cells
].m_Child
[0] >= N1
) {
721 int neigh
= m_Cells
[i_cells
].m_Neighbour
[0];
723 if ((m_Cells
[neigh
].m_Child
[0] != -1) && (m_Cells
[i_cells
].m_Level
== m_Cells
[neigh
].m_Level
)) {
724 m_Cells
[m_Cells
[i_cells
].m_Child
[0]].m_Neighbour
[0] = m_Cells
[neigh
].m_Child
[1];
725 m_Cells
[m_Cells
[i_cells
].m_Child
[2]].m_Neighbour
[0] = m_Cells
[neigh
].m_Child
[3];
726 m_Cells
[m_Cells
[i_cells
].m_Child
[4]].m_Neighbour
[0] = m_Cells
[neigh
].m_Child
[5];
727 m_Cells
[m_Cells
[i_cells
].m_Child
[6]].m_Neighbour
[0] = m_Cells
[neigh
].m_Child
[7];
729 m_Cells
[m_Cells
[i_cells
].m_Child
[0]].m_Neighbour
[0] = neigh
;
730 m_Cells
[m_Cells
[i_cells
].m_Child
[2]].m_Neighbour
[0] = neigh
;
731 m_Cells
[m_Cells
[i_cells
].m_Child
[4]].m_Neighbour
[0] = neigh
;
732 m_Cells
[m_Cells
[i_cells
].m_Child
[6]].m_Neighbour
[0] = neigh
;
737 int neigh
= m_Cells
[i_cells
].m_Neighbour
[1];
739 if (m_Cells
[neigh
].m_Child
[0] != -1) {
740 m_Cells
[m_Cells
[i_cells
].m_Child
[1]].m_Neighbour
[1] = m_Cells
[neigh
].m_Child
[0];
741 m_Cells
[m_Cells
[i_cells
].m_Child
[3]].m_Neighbour
[1] = m_Cells
[neigh
].m_Child
[2];
742 m_Cells
[m_Cells
[i_cells
].m_Child
[5]].m_Neighbour
[1] = m_Cells
[neigh
].m_Child
[4];
743 m_Cells
[m_Cells
[i_cells
].m_Child
[7]].m_Neighbour
[1] = m_Cells
[neigh
].m_Child
[6];
745 m_Cells
[m_Cells
[i_cells
].m_Child
[1]].m_Neighbour
[1] = neigh
;
746 m_Cells
[m_Cells
[i_cells
].m_Child
[3]].m_Neighbour
[1] = neigh
;
747 m_Cells
[m_Cells
[i_cells
].m_Child
[5]].m_Neighbour
[1] = neigh
;
748 m_Cells
[m_Cells
[i_cells
].m_Child
[7]].m_Neighbour
[1] = neigh
;
753 int neigh
= m_Cells
[i_cells
].m_Neighbour
[2];
755 if (m_Cells
[neigh
].m_Child
[0] != -1) {
756 m_Cells
[m_Cells
[i_cells
].m_Child
[0]].m_Neighbour
[2] = m_Cells
[neigh
].m_Child
[2];
757 m_Cells
[m_Cells
[i_cells
].m_Child
[1]].m_Neighbour
[2] = m_Cells
[neigh
].m_Child
[3];
758 m_Cells
[m_Cells
[i_cells
].m_Child
[4]].m_Neighbour
[2] = m_Cells
[neigh
].m_Child
[6];
759 m_Cells
[m_Cells
[i_cells
].m_Child
[5]].m_Neighbour
[2] = m_Cells
[neigh
].m_Child
[7];
761 m_Cells
[m_Cells
[i_cells
].m_Child
[0]].m_Neighbour
[2] = neigh
;
762 m_Cells
[m_Cells
[i_cells
].m_Child
[1]].m_Neighbour
[2] = neigh
;
763 m_Cells
[m_Cells
[i_cells
].m_Child
[4]].m_Neighbour
[2] = neigh
;
764 m_Cells
[m_Cells
[i_cells
].m_Child
[5]].m_Neighbour
[2] = neigh
;
769 int neigh
= m_Cells
[i_cells
].m_Neighbour
[3];
771 if (m_Cells
[neigh
].m_Child
[0] != -1) {
772 m_Cells
[m_Cells
[i_cells
].m_Child
[2]].m_Neighbour
[3] = m_Cells
[neigh
].m_Child
[0];
773 m_Cells
[m_Cells
[i_cells
].m_Child
[3]].m_Neighbour
[3] = m_Cells
[neigh
].m_Child
[1];
774 m_Cells
[m_Cells
[i_cells
].m_Child
[6]].m_Neighbour
[3] = m_Cells
[neigh
].m_Child
[4];
775 m_Cells
[m_Cells
[i_cells
].m_Child
[7]].m_Neighbour
[3] = m_Cells
[neigh
].m_Child
[5];
777 m_Cells
[m_Cells
[i_cells
].m_Child
[2]].m_Neighbour
[3] = neigh
;
778 m_Cells
[m_Cells
[i_cells
].m_Child
[3]].m_Neighbour
[3] = neigh
;
779 m_Cells
[m_Cells
[i_cells
].m_Child
[6]].m_Neighbour
[3] = neigh
;
780 m_Cells
[m_Cells
[i_cells
].m_Child
[7]].m_Neighbour
[3] = neigh
;
785 int neigh
= m_Cells
[i_cells
].m_Neighbour
[4];
787 if (m_Cells
[neigh
].m_Child
[0] != -1) {
788 m_Cells
[m_Cells
[i_cells
].m_Child
[0]].m_Neighbour
[4] = m_Cells
[neigh
].m_Child
[4];
789 m_Cells
[m_Cells
[i_cells
].m_Child
[1]].m_Neighbour
[4] = m_Cells
[neigh
].m_Child
[5];
790 m_Cells
[m_Cells
[i_cells
].m_Child
[2]].m_Neighbour
[4] = m_Cells
[neigh
].m_Child
[6];
791 m_Cells
[m_Cells
[i_cells
].m_Child
[3]].m_Neighbour
[4] = m_Cells
[neigh
].m_Child
[7];
793 m_Cells
[m_Cells
[i_cells
].m_Child
[0]].m_Neighbour
[4] = neigh
;
794 m_Cells
[m_Cells
[i_cells
].m_Child
[1]].m_Neighbour
[4] = neigh
;
795 m_Cells
[m_Cells
[i_cells
].m_Child
[2]].m_Neighbour
[4] = neigh
;
796 m_Cells
[m_Cells
[i_cells
].m_Child
[3]].m_Neighbour
[4] = neigh
;
801 int neigh
= m_Cells
[i_cells
].m_Neighbour
[5];
803 if (m_Cells
[neigh
].m_Child
[0] != -1) {
804 m_Cells
[m_Cells
[i_cells
].m_Child
[4]].m_Neighbour
[5] = m_Cells
[neigh
].m_Child
[0];
805 m_Cells
[m_Cells
[i_cells
].m_Child
[5]].m_Neighbour
[5] = m_Cells
[neigh
].m_Child
[1];
806 m_Cells
[m_Cells
[i_cells
].m_Child
[6]].m_Neighbour
[5] = m_Cells
[neigh
].m_Child
[2];
807 m_Cells
[m_Cells
[i_cells
].m_Child
[7]].m_Neighbour
[5] = m_Cells
[neigh
].m_Child
[3];
809 m_Cells
[m_Cells
[i_cells
].m_Child
[4]].m_Neighbour
[5] = neigh
;
810 m_Cells
[m_Cells
[i_cells
].m_Child
[5]].m_Neighbour
[5] = neigh
;
811 m_Cells
[m_Cells
[i_cells
].m_Child
[6]].m_Neighbour
[5] = neigh
;
812 m_Cells
[m_Cells
[i_cells
].m_Child
[7]].m_Neighbour
[5] = neigh
;
819 m_ToRefine
.fill(false, m_Cells
.size());
825 void Octree::toVtkGridHangingNodes(vtkUnstructuredGrid
*grid
, bool)
828 for (int i
= 0; i
< m_Cells
.size(); ++i
) {
829 if (m_Cells
[i
].m_Child
[0] == -1) {
833 allocateGrid(grid
, N
, m_Nodes
.size());
834 //allocateGrid(grid, N, m_Nodes.size(), create_fields);
835 for (int id_node
= 0; id_node
< m_Nodes
.size(); ++id_node
) {
836 grid
->GetPoints()->SetPoint(id_node
, m_Nodes
[id_node
].m_Position
.data());
838 EG_VTKDCC(vtkIntArray
, cell_index
, grid
, "cell_code");
839 for (int i_cells
= 0; i_cells
< m_Cells
.size(); ++i_cells
) {
840 OctreeCell cell
= m_Cells
[i_cells
];
841 if (cell
.m_Child
[0] == -1) {
844 pts
[0] = cell
.m_Node
[0];
845 pts
[1] = cell
.m_Node
[1];
846 pts
[2] = cell
.m_Node
[3];
847 pts
[3] = cell
.m_Node
[2];
848 pts
[4] = cell
.m_Node
[4];
849 pts
[5] = cell
.m_Node
[5];
850 pts
[6] = cell
.m_Node
[7];
851 pts
[7] = cell
.m_Node
[6];
852 for (int i
= 0; i
< 8; ++i
) {
853 if (pts
[i
] < 0 || pts
[i
] >= m_Nodes
.size()) {
857 vtkIdType id_cell
= grid
->InsertNextCell(VTK_HEXAHEDRON
, Npts
, pts
);
858 cell_index
->SetValue(id_cell
, i_cells
);
863 int Octree::opposingFace(int i
)
865 if (i
== 0) return 1;
866 if (i
== 1) return 0;
867 if (i
== 2) return 3;
868 if (i
== 3) return 2;
869 if (i
== 4) return 5;
870 if (i
== 5) return 4;
875 void Octree::toVtkGridPolyhedral(vtkUnstructuredGrid
*grid
, bool create_fields
)
877 if (!m_SmoothTransition
) {
881 vtkIdType num_new_cells
= 0;
882 for (int i_cell
= 0; i_cell
< getNumCells(); ++i_cell
) {
883 OctreeCell cell
= m_Cells
[i_cell
];
884 if (!cell
.hasChildren()) { // only use cells which do not have children
888 allocateGrid(grid
, num_new_cells
, m_Nodes
.size(), create_fields
);
889 for (vtkIdType id_node
= 0; id_node
< m_Nodes
.size(); ++id_node
) {
890 grid
->GetPoints()->SetPoint(id_node
, m_Nodes
[id_node
].getPosition().data());
893 for (int i_cell
= 0; i_cell
< getNumCells(); ++i_cell
) {
894 OctreeCell cell
= m_Cells
[i_cell
];
895 if (!cell
.hasChildren()) { // only use cells which do not have children
896 QList
<QVector
<int> > all_faces
;
897 QVector
<QVector
<int> > faces
;
898 for (int i
= 0; i
< 6; ++i
) {
899 bool use_neighbour_faces
= false;
900 if (cell
.getNeighbour(i
) != -1) {
901 OctreeCell neigh
= m_Cells
[cell
.getNeighbour(i
)];
902 if (neigh
.m_Level
== cell
.m_Level
) {
903 if (neigh
.hasChildren()) {
904 use_neighbour_faces
= true;
908 if (use_neighbour_faces
) {
909 m_Cells
[cell
.getNeighbour(i
)].getFaceNodes(opposingFace(i
), this, faces
, true);
911 cell
.getFaceNodes(i
, this, faces
, false);
913 foreach (QVector
<int> face
, faces
) {
914 all_faces
.push_back(face
);
917 if (all_faces
.size() < 6) {
920 bool simple_hex_cell
= true;
921 if (all_faces
.size() > 6) {
922 simple_hex_cell
= false;
925 foreach (QVector
<int> face
, all_faces
) {
926 if (face
.size() > 4) {
927 simple_hex_cell
= false;
932 if (simple_hex_cell
) {
935 pts
[0] = cell
.m_Node
[0];
936 pts
[1] = cell
.m_Node
[1];
937 pts
[2] = cell
.m_Node
[3];
938 pts
[3] = cell
.m_Node
[2];
939 pts
[4] = cell
.m_Node
[4];
940 pts
[5] = cell
.m_Node
[5];
941 pts
[6] = cell
.m_Node
[7];
942 pts
[7] = cell
.m_Node
[6];
943 for (int i
= 0; i
< 8; ++i
) {
944 if (pts
[i
] < 0 || pts
[i
] >= m_Nodes
.size()) {
948 vtkIdType id_cell
= grid
->InsertNextCell(VTK_HEXAHEDRON
, Npts
, pts
);
950 vtkIdType num_ids
= 1;
951 foreach (QVector
<int> face
, all_faces
) {
952 num_ids
+= 1 + face
.size();
954 EG_VTKSP(vtkIdList
, ptids
);
955 ptids
->SetNumberOfIds(num_ids
);
957 ptids
->SetId(i_id
++, all_faces
.size());
958 foreach (QVector
<int> face
, all_faces
) {
959 ptids
->SetId(i_id
++, face
.size());
960 foreach (int i
, face
) {
961 ptids
->SetId(i_id
++, i
);
964 grid
->InsertNextCell(VTK_POLYHEDRON
, ptids
);
970 void Octree::toVtkGridConforming(vtkUnstructuredGrid
* grid
, bool create_fields
)
972 if (!m_SmoothTransition
) {
976 int num_new_nodes
= 0;
977 int num_pyramids
= 0;
980 for (int i_cell
= 0; i_cell
< getNumCells(); ++i_cell
) {
981 OctreeCell cell
= m_Cells
[i_cell
];
982 if (!cell
.hasChildren()) { // only use cells which do not have children
983 QList
<QVector
<int> > all_faces
;
984 QVector
<QVector
<int> > faces
;
985 for (int i
= 0; i
< 6; ++i
) {
986 bool use_neighbour_faces
= false;
987 if (cell
.getNeighbour(i
) != -1) {
988 OctreeCell neigh
= m_Cells
[cell
.getNeighbour(i
)];
989 if (neigh
.m_Level
== cell
.m_Level
) {
990 if (neigh
.hasChildren()) {
991 use_neighbour_faces
= true;
995 if (use_neighbour_faces
) {
996 m_Cells
[cell
.getNeighbour(i
)].getFaceNodes(opposingFace(i
), this, faces
);
998 cell
.getFaceNodes(i
, this, faces
, true);
1000 foreach (QVector
<int> face
, faces
) {
1001 all_faces
.push_back(face
);
1004 if (all_faces
.size() < 6) {
1007 bool simple_hex_cell
= true;
1008 if (all_faces
.size() > 6) {
1009 simple_hex_cell
= false;
1012 foreach (QVector
<int> face
, all_faces
) {
1013 if (face
.size() > 4) {
1014 simple_hex_cell
= false;
1019 if (simple_hex_cell
) {
1023 foreach (QVector
<int> face
, all_faces
) {
1024 if (face
.size() > 4) {
1025 num_tetras
+= face
.size();
1035 allocateGrid(grid
, num_hexes
+ num_pyramids
+ num_tetras
, m_Nodes
.size() + num_new_nodes
, create_fields
);
1036 vtkIdType id_node
= 0;
1037 for (int i
= 0; i
< m_Nodes
.size(); ++i
) {
1038 vec3_t x
= m_Nodes
[i
].getPosition().data();
1039 grid
->GetPoints()->SetPoint(id_node
, m_Nodes
[i
].getPosition().data());
1043 QVector
<QList
<node_t
> > face_nodes(m_Cells
.size());
1045 for (int i_cells
= 0; i_cells
< m_Cells
.size(); ++i_cells
) {
1046 OctreeCell cell
= m_Cells
[i_cells
];
1047 if (!cell
.hasChildren()) {
1048 QList
<QVector
<int> > all_faces
;
1049 QVector
<QVector
<int> > faces
;
1050 for (int i
= 0; i
< 6; ++i
) {
1051 bool use_neighbour_faces
= false;
1052 if (cell
.getNeighbour(i
) != -1) {
1053 OctreeCell neigh
= m_Cells
[cell
.getNeighbour(i
)];
1054 if (neigh
.m_Level
== cell
.m_Level
) {
1055 if (neigh
.m_Child
[0] != -1) {
1056 use_neighbour_faces
= true;
1060 if (use_neighbour_faces
) {
1061 m_Cells
[cell
.getNeighbour(i
)].getFaceNodes(opposingFace(i
), this, faces
);
1063 cell
.getFaceNodes(i
, this, faces
, true);
1065 foreach (QVector
<int> face
, faces
) {
1066 all_faces
.push_back(face
);
1069 bool simple_hex_cell
= true;
1070 if (all_faces
.size() > 6) {
1071 simple_hex_cell
= false;
1073 foreach (QVector
<int> face
, all_faces
) {
1074 if (face
.size() > 4) {
1075 simple_hex_cell
= false;
1080 if (simple_hex_cell
) {
1082 pts
[0] = cell
.m_Node
[0];
1083 pts
[1] = cell
.m_Node
[1];
1084 pts
[2] = cell
.m_Node
[3];
1085 pts
[3] = cell
.m_Node
[2];
1086 pts
[4] = cell
.m_Node
[4];
1087 pts
[5] = cell
.m_Node
[5];
1088 pts
[6] = cell
.m_Node
[7];
1089 pts
[7] = cell
.m_Node
[6];
1090 grid
->InsertNextCell(VTK_HEXAHEDRON
, 8, pts
);
1092 grid
->GetPoints()->SetPoint(id_node
, getCellCentre(i_cells
).data());
1093 int id_cell_centre
= id_node
;
1095 foreach (QVector
<int> face
, all_faces
) {
1096 if (face
.size() == 4) {
1098 for (int i
= 0; i
< 4; ++i
) {
1101 pts
[4] = id_cell_centre
;
1102 grid
->InsertNextCell(VTK_PYRAMID
, 5, pts
);
1106 for (int i
= 0; i
< face
.size(); ++i
) {
1107 xf
+= m_Nodes
[face
[i
]].getPosition();
1109 xf
*= 1.0/face
.size();
1111 double f13
= 1.0/3.0;
1112 for (int i
= 0; i
< face
.size(); ++i
) {
1114 grid
->GetPoint(face
[i
], x1
.data());
1115 if (i
< face
.size()-1) {
1116 grid
->GetPoint(face
[i
+1], x2
.data());
1118 grid
->GetPoint(face
[0], x2
.data());
1120 double A
= GeometryTools::triArea(xf
, x1
, x2
);
1122 xc
+= f13
*A
*(x1
+ x2
+ xf
);
1126 int id_face_centre
= -1;
1127 double l_crit
= 0.001*min(getDx(i_cells
), min(getDy(i_cells
), getDz(i_cells
)));
1128 QList
<int> neighbours
;
1129 for (int i
= 0; i
< 6; ++i
) {
1130 if (cell
.getNeighbour(i
) != -1) {
1131 neighbours
.append(cell
.getNeighbour(i
));
1134 foreach (int i_neigh
, neighbours
) {
1135 foreach (node_t N
, face_nodes
[i_neigh
]) {
1136 if ((N
.x
- xc
).abs() < l_crit
) {
1137 id_face_centre
= N
.id
;
1141 if (id_face_centre
== -1) {
1142 grid
->GetPoints()->SetPoint(id_node
, xc
.data());
1143 id_face_centre
= id_node
;
1147 N
.id
= id_face_centre
;
1148 face_nodes
[i_cells
].append(N
);
1150 for (int i
= 0; i
< face
.size(); ++i
) {
1153 if (i
< face
.size()-1) {
1158 pts
[2] = id_face_centre
;
1159 pts
[3] = id_cell_centre
;
1160 grid
->InsertNextCell(VTK_TETRA
, 4, pts
);
1168 DeleteStrayNodes del_stray
;
1169 del_stray
.setGrid(grid
);
1170 del_stray
.setAllCells();
1175 vec3_t
Octree::getCellCentre(int cell
)
1178 for (int i
= 0; i
< 8; ++i
) {
1179 x
+= m_Nodes
[m_Cells
[cell
].m_Node
[i
]].m_Position
;
1185 vec3_t
Octree::getFaceCentre(int i_cells
, int i_faces
)
1188 const OctreeCell
& cell
= m_Cells
[i_cells
];
1190 x
+= m_Nodes
[cell
.m_Node
[0]].m_Position
;
1191 x
+= m_Nodes
[cell
.m_Node
[2]].m_Position
;
1192 x
+= m_Nodes
[cell
.m_Node
[4]].m_Position
;
1193 x
+= m_Nodes
[cell
.m_Node
[6]].m_Position
;
1194 } else if (i_faces
== 1) {
1195 x
+= m_Nodes
[cell
.m_Node
[1]].m_Position
;
1196 x
+= m_Nodes
[cell
.m_Node
[3]].m_Position
;
1197 x
+= m_Nodes
[cell
.m_Node
[5]].m_Position
;
1198 x
+= m_Nodes
[cell
.m_Node
[7]].m_Position
;
1199 } else if (i_faces
== 2) {
1200 x
+= m_Nodes
[cell
.m_Node
[0]].m_Position
;
1201 x
+= m_Nodes
[cell
.m_Node
[1]].m_Position
;
1202 x
+= m_Nodes
[cell
.m_Node
[4]].m_Position
;
1203 x
+= m_Nodes
[cell
.m_Node
[5]].m_Position
;
1204 } else if (i_faces
== 3) {
1205 x
+= m_Nodes
[cell
.m_Node
[2]].m_Position
;
1206 x
+= m_Nodes
[cell
.m_Node
[3]].m_Position
;
1207 x
+= m_Nodes
[cell
.m_Node
[6]].m_Position
;
1208 x
+= m_Nodes
[cell
.m_Node
[7]].m_Position
;
1209 } else if (i_faces
== 4) {
1210 x
+= m_Nodes
[cell
.m_Node
[0]].m_Position
;
1211 x
+= m_Nodes
[cell
.m_Node
[1]].m_Position
;
1212 x
+= m_Nodes
[cell
.m_Node
[2]].m_Position
;
1213 x
+= m_Nodes
[cell
.m_Node
[3]].m_Position
;
1214 } else if (i_faces
== 5) {
1215 x
+= m_Nodes
[cell
.m_Node
[4]].m_Position
;
1216 x
+= m_Nodes
[cell
.m_Node
[5]].m_Position
;
1217 x
+= m_Nodes
[cell
.m_Node
[6]].m_Position
;
1218 x
+= m_Nodes
[cell
.m_Node
[7]].m_Position
;
1224 int Octree::findCell(vec3_t x
)
1226 for (int i
= 0; i
< 3; ++i
) {
1227 if ((x
[i
] < m_Corner1
[i
]) || (x
[i
] > m_Corner2
[i
])) {
1228 EG_ERR_RETURN("node outside of octree mesh");
1232 while (getLevel(i_cells
) == 0) {
1233 vec3_t x1
= m_Nodes
[m_Cells
[i_cells
].m_Node
[0]].getPosition();
1234 vec3_t x2
= m_Nodes
[m_Cells
[i_cells
].m_Node
[7]].getPosition();
1235 if (x
[0] >= x1
[0] && x
[1] >= x1
[1] && x
[2] >= x1
[2] && x
[0] <= x2
[0] && x
[1] <= x2
[1] && x
[2] <= x2
[2]) {
1240 if (getLevel(i_cells
) > 0) {
1243 while (hasChildren(i_cells
)) {
1244 vec3_t xc
= getCellCentre(i_cells
);
1248 i_cells
= m_Cells
[i_cells
].m_Child
[7];
1250 i_cells
= m_Cells
[i_cells
].m_Child
[3];
1254 i_cells
= m_Cells
[i_cells
].m_Child
[5];
1256 i_cells
= m_Cells
[i_cells
].m_Child
[1];
1262 i_cells
= m_Cells
[i_cells
].m_Child
[6];
1264 i_cells
= m_Cells
[i_cells
].m_Child
[2];
1268 i_cells
= m_Cells
[i_cells
].m_Child
[4];
1270 i_cells
= m_Cells
[i_cells
].m_Child
[0];
1278 bool Octree::intersectsFace(int cell
, int face
, vec3_t x1
, vec3_t x2
, double scale
, double &k
, double tol
)
1281 vec3_t x_centre
= getCellCentre(cell
);
1283 a
= m_Nodes
[m_Cells
[cell
].m_Node
[0]].m_Position
;
1284 b
= m_Nodes
[m_Cells
[cell
].m_Node
[2]].m_Position
;
1285 c
= m_Nodes
[m_Cells
[cell
].m_Node
[4]].m_Position
;
1286 } else if (face
== 1) {
1287 a
= m_Nodes
[m_Cells
[cell
].m_Node
[1]].m_Position
;
1288 b
= m_Nodes
[m_Cells
[cell
].m_Node
[3]].m_Position
;
1289 c
= m_Nodes
[m_Cells
[cell
].m_Node
[5]].m_Position
;
1290 } else if (face
== 2) {
1291 a
= m_Nodes
[m_Cells
[cell
].m_Node
[0]].m_Position
;
1292 b
= m_Nodes
[m_Cells
[cell
].m_Node
[1]].m_Position
;
1293 c
= m_Nodes
[m_Cells
[cell
].m_Node
[4]].m_Position
;
1294 } else if (face
== 3) {
1295 a
= m_Nodes
[m_Cells
[cell
].m_Node
[2]].m_Position
;
1296 b
= m_Nodes
[m_Cells
[cell
].m_Node
[3]].m_Position
;
1297 c
= m_Nodes
[m_Cells
[cell
].m_Node
[6]].m_Position
;
1298 } else if (face
== 4) {
1299 a
= m_Nodes
[m_Cells
[cell
].m_Node
[0]].m_Position
;
1300 b
= m_Nodes
[m_Cells
[cell
].m_Node
[1]].m_Position
;
1301 c
= m_Nodes
[m_Cells
[cell
].m_Node
[2]].m_Position
;
1302 } else if (face
== 5) {
1303 a
= m_Nodes
[m_Cells
[cell
].m_Node
[4]].m_Position
;
1304 b
= m_Nodes
[m_Cells
[cell
].m_Node
[5]].m_Position
;
1305 c
= m_Nodes
[m_Cells
[cell
].m_Node
[6]].m_Position
;
1307 a
= x_centre
+ scale
*(a
- x_centre
);
1308 b
= x_centre
+ scale
*(b
- x_centre
);
1309 c
= x_centre
+ scale
*(c
- x_centre
);
1312 double g1abs
= g1
.abs();
1313 double g2abs
= g2
.abs();
1316 vec3_t n
= g1
.cross(g2
);
1317 k
= GeometryTools::intersection(x1
, x2
-x1
, a
, n
);
1318 bool intersects
= false;
1319 if ((k
> 0 - tol
) && (k
< 1 + tol
)) {
1320 vec3_t x
= x1
+ k
*(x2
-x1
) - a
;
1325 if (fabs(x
*n
) > 1e-4) {
1328 if ((xg1
> 0 - tol
) && (xg1
< 1 + tol
) && (xg2
> 0 - tol
) && (xg2
< 1 + tol
)) {
1335 void Octree::resetRefineMarks()
1337 m_ToRefine
.fill(false, m_Cells
.size());
1340 void Octree::getEdges(int cell
, QVector
<SortedPair
<int> >& edges
)
1343 edges
[0].v1
= getNode(cell
, 0); edges
[0].v2
= getNode(cell
, 1);
1344 edges
[1].v1
= getNode(cell
, 0); edges
[1].v2
= getNode(cell
, 2);
1345 edges
[2].v1
= getNode(cell
, 0); edges
[2].v2
= getNode(cell
, 4);
1346 edges
[3].v1
= getNode(cell
, 1); edges
[3].v2
= getNode(cell
, 3);
1347 edges
[4].v1
= getNode(cell
, 1); edges
[4].v2
= getNode(cell
, 5);
1348 edges
[5].v1
= getNode(cell
, 2); edges
[5].v2
= getNode(cell
, 3);
1349 edges
[6].v1
= getNode(cell
, 2); edges
[6].v2
= getNode(cell
, 6);
1350 edges
[7].v1
= getNode(cell
, 3); edges
[7].v2
= getNode(cell
, 7);
1351 edges
[8].v1
= getNode(cell
, 4); edges
[8].v2
= getNode(cell
, 5);
1352 edges
[9].v1
= getNode(cell
, 4); edges
[9].v2
= getNode(cell
, 6);
1353 edges
[10].v1
= getNode(cell
, 5); edges
[10].v2
= getNode(cell
, 7);
1354 edges
[11].v1
= getNode(cell
, 6); edges
[11].v2
= getNode(cell
, 7);
1357 bool Octree::isInsideBounds(vec3_t x
)
1360 for (int i
= 0; i
< 3; ++i
) {
1361 if ((x
[i
] < m_Corner1
[i
]) || (x
[i
] > m_Corner2
[i
])) {
1369 bool Octree::isInsideCell(int cell
, vec3_t x
, double overlap
)
1371 vec3_t x0_1
= getNodePosition(cell
, 0);
1372 vec3_t x0_2
= getNodePosition(cell
, 7);
1373 vec3_t dx
= x0_2
- x0_1
;
1374 vec3_t x1
= x0_1
- overlap
*dx
;
1375 vec3_t x2
= x0_2
+ overlap
*dx
;
1377 for (int i
= 0; i
< 3; ++i
) {
1378 if (x0_1
[i
] >= x0_2
[i
]) {
1381 if (x1
[i
] >= x2
[i
]) {
1387 for (int i
= 0; i
< 3; ++i
) {
1388 if ((x
[i
] < x1
[i
]) || (x
[i
] > x2
[i
])) {
1396 bool Octree::triangleIntersectsCell(int cell
, QVector
<vec3_t
> tri
, double scale
)
1398 if (tri
.size() != 3) {
1402 // any node inside cell?
1403 foreach (vec3_t x
, tri
) {
1404 if (isInsideCell(cell
, x
, (scale
- 1))) {
1409 // any edge of cell intersects triangle?
1410 QVector
<SortedPair
<int> > edges
;
1411 getEdges(cell
, edges
);
1412 foreach (SortedPair
<int> edge
, edges
) {
1413 vec3_t x1
= getNodePosition(edge
.v1
);
1414 vec3_t x2
= getNodePosition(edge
.v2
);
1416 vec3_t x_centre
= getCellCentre(cell
);
1417 x1
= x_centre
+ scale
*(x1
- x_centre
);
1418 x2
= x_centre
+ scale
*(x2
- x_centre
);
1419 if (GeometryTools::intersectEdgeAndTriangle(tri
[0], tri
[1], tri
[2], x1
, x2
, xi
, ri
)) {
1424 // any edge of triangle intersects any face of cell?
1425 for (int face
= 0; face
< 6; ++face
) {
1427 if (intersectsFace(cell
, face
, tri
[0], tri
[1], scale
, k
)) {
1430 if (intersectsFace(cell
, face
, tri
[1], tri
[2], scale
, k
)) {
1433 if (intersectsFace(cell
, face
, tri
[2], tri
[0], scale
, k
)) {
1441 void Octree::getFinestChildren(int cell
, QList
<int> &finest_children
)
1443 if (hasChildren(cell
)) {
1444 for (int i
= 0; i
< 8; ++i
) {
1445 getFinestChildren(m_Cells
[cell
].m_Child
[i
], finest_children
);
1448 finest_children
<< cell
;
1452 void Octree::getNeighbourRegion(int cell
, QList
<int> &neighbour_cells
)
1455 for (int i
= 0; i
< 8; ++i
) {
1456 int node
= getNode(cell
, i
);
1457 for (int j
= 0; j
< m_Node2Cell
[node
].size(); ++j
) {
1458 int cell
= m_Node2Cell
[node
][j
];
1459 if (!hasChildren(cell
)) {
1464 foreach (int cell
, cells
) {
1465 neighbour_cells
.append(cell
);