fixed edge display for volume cells
[engrid-github.git] / src / libengrid / octree.cpp
blobab0dd307ebff24e200711466fdf66a932c7644ae
1 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 // + +
3 // + This file is part of enGrid. +
4 // + +
5 // + Copyright 2008-2014 enGits GmbH +
6 // + +
7 // + enGrid is free software: you can redistribute it and/or modify +
8 // + it under the terms of the GNU General Public License as published by +
9 // + the Free Software Foundation, either version 3 of the License, or +
10 // + (at your option) any later version. +
11 // + +
12 // + enGrid is distributed in the hope that it will be useful, +
13 // + but WITHOUT ANY WARRANTY; without even the implied warranty of +
14 // + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +
15 // + GNU General Public License for more details. +
16 // + +
17 // + You should have received a copy of the GNU General Public License +
18 // + along with enGrid. If not, see <http://www.gnu.org/licenses/>. +
19 // + +
20 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
22 #include "octree.h"
23 #include "geometrytools.h"
24 #include "deletestraynodes.h"
25 #include "guimainwindow.h"
27 OctreeCell::OctreeCell()
29 for (int i = 0; i < 8; ++i) {
30 m_Node[i] = -1;
31 m_Child[i] = -1;
33 for (int i = 0; i < 6; ++i) {
34 m_Neighbour[i] = -1;
36 m_Parent = -1;
37 m_Level = 0;
40 int OctreeCell::findNode(int i_node)
42 int i_local = -1;
43 for (int i = 0; i < 8; ++i) {
44 if (m_Node[i] == i_node) {
45 i_local = i;
46 break;
49 return i_local;
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);
57 int edge_node = -1;
58 if (hasChildren()) {
59 if (n1 > n2) {
60 swap(n1, n2);
62 if (n1 == 0) {
63 if (n2 == 1 || n2 == 2 || n2 == 4) {
64 edge_node = octree->m_Cells[m_Child[n1]].m_Node[n2];
66 } else if (n1 == 1) {
67 if (n2 == 3 || n2 == 5) {
68 edge_node = octree->m_Cells[m_Child[n1]].m_Node[n2];
70 } else if (n1 == 2) {
71 if (n2 == 3 || n2 == 6) {
72 edge_node = octree->m_Cells[m_Child[n1]].m_Node[n2];
74 } else if (n1 == 3) {
75 if (n2 == 7) {
76 edge_node = octree->m_Cells[m_Child[n1]].m_Node[n2];
78 } else if (n1 == 4) {
79 if (n2 == 5 || n2 == 6) {
80 edge_node = octree->m_Cells[m_Child[n1]].m_Node[n2];
82 } else if (n1 == 5) {
83 if (n2 == 7) {
84 edge_node = octree->m_Cells[m_Child[n1]].m_Node[n2];
86 } else if (n1 == 6) {
87 if (n2 == 7) {
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) {
101 break;
106 return edge_node;
109 void OctreeCell::getFaceNodes(int i, Octree* octree, QVector<int>& face_nodes, bool reverse)
111 //face_nodes.resize(4);
112 QList<int> nodes;
113 int edge_node;
114 if (i == 0) {
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);
123 } else if (i == 1) {
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);
132 } else if (i == 2) {
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);
141 } else if (i == 3) {
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);
150 } else if (i == 4) {
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);
159 } else if (i == 5) {
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());
171 if (reverse) {
172 int i = face_nodes.size() - 1;
173 foreach(int node, nodes) {
174 face_nodes[i] = node;
175 --i;
177 } else {
178 int i = 0;
179 foreach(int node, nodes) {
180 face_nodes[i] = node;
181 ++i;
186 void OctreeCell::getFaceNodes(int i, Octree* octree, QVector<QVector<int> >& face_nodes, bool reverse)
188 if (hasChildren()) {
189 face_nodes.resize(4);
190 if (i == 0) {
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);
195 } else if (i == 1) {
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);
200 } else if (i == 2) {
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);
205 } else if (i == 3) {
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);
210 } else if (i == 4) {
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);
215 } else if (i == 5) {
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);
221 } else {
222 face_nodes.resize(1);
223 getFaceNodes(i, octree, face_nodes[0], reverse);
227 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
229 Octree::Octree()
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));
233 m_Nodes.resize(8);
234 m_Cells.resize(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) {
251 markToRefine(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;
263 int max_level = 0;
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);
267 buildNode2Cell();
269 for (int i_node = 0; i_node < getNumNodes(); ++i_node) {
270 QSet<int> canditates;
271 double tol = 1e99;
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]);
283 if (tol > 1e98) {
284 EG_BUG;
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;
292 } else {
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]) {
312 EG_BUG;
315 int N = 0;
316 for (int i_nodes = 0; i_nodes < m_Nodes.size(); ++i_nodes) {
317 if (is_dup[i_nodes]) {
318 ++N;
323 void Octree::mergeNodes_compactNodes()
325 QVector<int> offset(m_Nodes.size());
326 int last_offset = 0;
327 double max_dist = 0;
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) {
332 EG_BUG;
334 ++last_offset;
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]];
342 } else {
343 m_SameNodes[i] -= offset[i];
347 for (int i = 0; i < m_Nodes.size(); ++i) {
348 if (m_SameNodes[i] > i) {
349 EG_BUG;
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()
377 int Nerr = 0;
378 int Nno = 0;
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) {
385 ++i_back_faces;
386 } else {
387 --i_back_faces;
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;
393 ++Nerr;
396 } else {
397 cout << "no neighbour: " << i_cells << ',' << i_faces << ',' << getFaceCentre(i_cells, i_faces) << endl;
398 ++Nno;
402 cout << Nerr << " errors and " << Nno << " faces without neighbour" << endl;
405 void Octree::buildNode2Cell()
407 m_Node2Cell.clear();
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) {
415 int max_level = 1;
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)
431 m_Origin = 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);
452 m_Corner1 = corner1;
453 m_Corner2 = corner2;
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()
493 int N1;
494 int N2;
495 int count = 0;
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;
501 buildNode2Cell();
502 do {
503 N1 = 0;
504 N2 = 0;
505 for (int i_cells = 0; i_cells < m_Cells.size(); ++i_cells) {
506 if (m_ToRefine[i_cells]) {
507 ++N1;
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)) {
512 markToRefine(neigh);
513 ++N2;
520 ++count;
521 if (count > 100) {
522 EG_BUG;
524 } while (N2 > 0);
525 N2 = m_Cells.size();
526 if (N2 + 8*N1 > m_MaxCells) {
527 QString num;
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";
533 EG_ERR_RETURN(msg);
535 m_Cells.insert(N2, 8*N1, OctreeCell());
536 int Nrefine = N1;
537 int new_node = m_Nodes.size();
538 m_Nodes.insert(m_Nodes.size(), N1*19, OctreeNode());
539 N1 = N2;
540 for (int i_cells = 0; i_cells < N1; ++i_cells) {
541 if (m_ToRefine[i_cells]) {
542 if (hasChildren(i_cells)) {
543 EG_BUG;
545 int nn[8];
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];
554 // create new nodes
555 int ne[12];
556 ne[0] = new_node;
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);
558 ne[1] = new_node;
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);
560 ne[2] = new_node;
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);
562 ne[3] = new_node;
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);
564 ne[4] = new_node;
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);
566 ne[5] = new_node;
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);
568 ne[6] = new_node;
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);
570 ne[7] = new_node;
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);
572 ne[8] = new_node;
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);
574 ne[9] = new_node;
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);
576 ne[10] = new_node;
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);
578 ne[11] = new_node;
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);
580 int nf[6];
581 nf[0] = new_node;
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);
584 nf[1] = new_node;
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);
587 nf[2] = new_node;
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);
590 nf[3] = new_node;
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);
593 nf[4] = new_node;
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);
596 nf[5] = new_node;
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);
599 int nv = new_node;
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;
607 ++N2;
610 // child 0
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;
619 // child 1
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];
628 // child 2
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];
637 // child 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];
646 // child 4
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];
655 // child 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];
664 // child 6
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];
673 // child 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];
683 // - - -
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];
687 // + - -
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];
691 // - + -
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];
695 // + + -
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];
699 // - - +
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];
703 // + - +
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];
707 // - + +
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];
711 // + + +
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];
722 if (neigh != -1) {
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];
728 } else {
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];
738 if (neigh != -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];
744 } else {
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];
754 if (neigh != -1) {
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];
760 } else {
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];
770 if (neigh != -1) {
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];
776 } else {
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];
786 if (neigh != -1) {
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];
792 } else {
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];
802 if (neigh != -1) {
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];
808 } else {
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());
820 mergeNodes();
821 buildNode2Cell();
822 return Nrefine;
825 void Octree::toVtkGridHangingNodes(vtkUnstructuredGrid *grid, bool)
827 int N = 0;
828 for (int i = 0; i < m_Cells.size(); ++i) {
829 if (m_Cells[i].m_Child[0] == -1) {
830 ++N;
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) {
842 vtkIdType Npts = 8;
843 vtkIdType pts[8];
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()) {
854 EG_BUG;
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;
871 EG_BUG;
872 return -1;
875 void Octree::toVtkGridPolyhedral(vtkUnstructuredGrid *grid, bool create_fields)
877 if (!m_SmoothTransition) {
878 EG_BUG;
880 buildNode2Cell();
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
885 ++num_new_cells;
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);
910 } else {
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) {
918 EG_BUG;
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;
928 break;
932 if (simple_hex_cell) {
933 vtkIdType Npts = 8;
934 vtkIdType pts[8];
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()) {
945 EG_BUG;
948 vtkIdType id_cell = grid->InsertNextCell(VTK_HEXAHEDRON, Npts, pts);
949 } else {
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);
956 vtkIdType i_id = 0;
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) {
973 EG_BUG;
975 buildNode2Cell();
976 int num_new_nodes = 0;
977 int num_pyramids = 0;
978 int num_hexes = 0;
979 int num_tetras = 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);
997 } else {
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) {
1005 EG_BUG;
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;
1015 break;
1019 if (simple_hex_cell) {
1020 ++num_hexes;
1021 } else {
1022 ++num_new_nodes;
1023 foreach (QVector<int> face, all_faces) {
1024 if (face.size() > 4) {
1025 num_tetras += face.size();
1026 ++num_new_nodes;
1027 } else {
1028 ++num_pyramids;
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());
1040 ++id_node;
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);
1062 } else {
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;
1072 } else {
1073 foreach (QVector<int> face, all_faces) {
1074 if (face.size() > 4) {
1075 simple_hex_cell = false;
1076 break;
1080 if (simple_hex_cell) {
1081 vtkIdType pts[8];
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);
1091 } else {
1092 grid->GetPoints()->SetPoint(id_node, getCellCentre(i_cells).data());
1093 int id_cell_centre = id_node;
1094 ++id_node;
1095 foreach (QVector<int> face, all_faces) {
1096 if (face.size() == 4) {
1097 vtkIdType pts[5];
1098 for (int i = 0; i < 4; ++i) {
1099 pts[i] = face[i];
1101 pts[4] = id_cell_centre;
1102 grid->InsertNextCell(VTK_PYRAMID, 5, pts);
1103 } else {
1104 vec3_t xf(0,0,0);
1105 vec3_t xc(0,0,0);
1106 for (int i = 0; i < face.size(); ++i) {
1107 xf += m_Nodes[face[i]].getPosition();
1109 xf *= 1.0/face.size();
1110 double Atot = 0;
1111 double f13 = 1.0/3.0;
1112 for (int i = 0; i < face.size(); ++i) {
1113 vec3_t x1, x2;
1114 grid->GetPoint(face[i], x1.data());
1115 if (i < face.size()-1) {
1116 grid->GetPoint(face[i+1], x2.data());
1117 } else {
1118 grid->GetPoint(face[0], x2.data());
1120 double A = GeometryTools::triArea(xf, x1, x2);
1121 Atot += A;
1122 xc += f13*A*(x1 + x2 + xf);
1124 xc *= 1.0/Atot;
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;
1144 ++id_node;
1145 node_t N;
1146 N.x = xc;
1147 N.id = id_face_centre;
1148 face_nodes[i_cells].append(N);
1150 for (int i = 0; i < face.size(); ++i) {
1151 vtkIdType pts[4];
1152 pts[0] = face[i];
1153 if (i < face.size()-1) {
1154 pts[1] = face[i+1];
1155 } else {
1156 pts[1] = face[0];
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();
1171 del_stray();
1175 vec3_t Octree::getCellCentre(int cell)
1177 vec3_t x(0,0,0);
1178 for (int i = 0; i < 8; ++i) {
1179 x += m_Nodes[m_Cells[cell].m_Node[i]].m_Position;
1181 x *= 1.0/8.0;
1182 return x;
1185 vec3_t Octree::getFaceCentre(int i_cells, int i_faces)
1187 vec3_t x(0,0,0);
1188 const OctreeCell& cell = m_Cells[i_cells];
1189 if (i_faces == 0) {
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;
1220 x *= 0.25;
1221 return x;
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");
1231 int i_cells = 0;
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]) {
1236 break;
1238 ++i_cells;
1240 if (getLevel(i_cells) > 0) {
1241 EG_BUG;
1243 while (hasChildren(i_cells)) {
1244 vec3_t xc = getCellCentre(i_cells);
1245 if (x[0] > xc[0]) {
1246 if (x[1] > xc[1]) {
1247 if (x[2] > xc[2]) {
1248 i_cells = m_Cells[i_cells].m_Child[7];
1249 } else {
1250 i_cells = m_Cells[i_cells].m_Child[3];
1252 } else {
1253 if (x[2] > xc[2]) {
1254 i_cells = m_Cells[i_cells].m_Child[5];
1255 } else {
1256 i_cells = m_Cells[i_cells].m_Child[1];
1259 } else {
1260 if (x[1] > xc[1]) {
1261 if (x[2] > xc[2]) {
1262 i_cells = m_Cells[i_cells].m_Child[6];
1263 } else {
1264 i_cells = m_Cells[i_cells].m_Child[2];
1266 } else {
1267 if (x[2] > xc[2]) {
1268 i_cells = m_Cells[i_cells].m_Child[4];
1269 } else {
1270 i_cells = m_Cells[i_cells].m_Child[0];
1275 return i_cells;
1278 bool Octree::intersectsFace(int cell, int face, vec3_t x1, vec3_t x2, double scale, double &k, double tol)
1280 vec3_t a, b, c;
1281 vec3_t x_centre = getCellCentre(cell);
1282 if (face == 0) {
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);
1310 vec3_t g1 = b-a;
1311 vec3_t g2 = c-a;
1312 double g1abs = g1.abs();
1313 double g2abs = g2.abs();
1314 g1.normalise();
1315 g2.normalise();
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;
1321 double xg1 = x*g1;
1322 double xg2 = x*g2;
1323 xg1 /= g1abs;
1324 xg2 /= g2abs;
1325 if (fabs(x*n) > 1e-4) {
1326 EG_BUG;
1328 if ((xg1 > 0 - tol) && (xg1 < 1 + tol) && (xg2 > 0 - tol) && (xg2 < 1 + tol)) {
1329 intersects = true;
1332 return intersects;
1335 void Octree::resetRefineMarks()
1337 m_ToRefine.fill(false, m_Cells.size());
1340 void Octree::getEdges(int cell, QVector<SortedPair<int> >& edges)
1342 edges.resize(12);
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)
1359 bool inside = true;
1360 for (int i = 0; i < 3; ++i) {
1361 if ((x[i] < m_Corner1[i]) || (x[i] > m_Corner2[i])) {
1362 inside = false;
1363 break;
1366 return inside;
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]) {
1379 EG_BUG;
1381 if (x1[i] >= x2[i]) {
1382 EG_BUG;
1386 bool inside = true;
1387 for (int i = 0; i < 3; ++i) {
1388 if ((x[i] < x1[i]) || (x[i] > x2[i])) {
1389 inside = false;
1390 break;
1393 return inside;
1396 bool Octree::triangleIntersectsCell(int cell, QVector<vec3_t> tri, double scale)
1398 if (tri.size() != 3) {
1399 EG_BUG;
1402 // any node inside cell?
1403 foreach (vec3_t x, tri) {
1404 if (isInsideCell(cell, x, (scale - 1))) {
1405 return true;
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);
1415 vec3_t xi, ri;
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)) {
1420 return true;
1424 // any edge of triangle intersects any face of cell?
1425 for (int face = 0; face < 6; ++face) {
1426 double k;
1427 if (intersectsFace(cell, face, tri[0], tri[1], scale, k)) {
1428 return true;
1430 if (intersectsFace(cell, face, tri[1], tri[2], scale, k)) {
1431 return true;
1433 if (intersectsFace(cell, face, tri[2], tri[0], scale, k)) {
1434 return true;
1438 return false;
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);
1447 } else {
1448 finest_children << cell;
1452 void Octree::getNeighbourRegion(int cell, QList<int> &neighbour_cells)
1454 QSet<int> 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)) {
1460 cells.insert(cell);
1464 foreach (int cell, cells) {
1465 neighbour_cells.append(cell);