?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
fixed edge display for volume cells
[engrid-github.git] / src / libengrid / polymolecule.cpp
blob25d35bb6fbbb5bdde1deb8f3425ab9a1fab29c52
1 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 // + +
3 // + This file is part of enGrid. +
4 // + +
5 // + Copyright 2008-2014 enGits GmbH +
6 // + +
7 // + enGrid is free software: you can redistribute it and/or modify +
8 // + it under the terms of the GNU General Public License as published by +
9 // + the Free Software Foundation, either version 3 of the License, or +
10 // + (at your option) any later version. +
11 // + +
12 // + enGrid is distributed in the hope that it will be useful, +
13 // + but WITHOUT ANY WARRANTY; without even the implied warranty of +
14 // + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +
15 // + GNU General Public License for more details. +
16 // + +
17 // + You should have received a copy of the GNU General Public License +
18 // + along with enGrid. If not, see <http://www.gnu.org/licenses/>. +
19 // + +
20 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
21
22 #include "polymolecule.h"
23 #include "guimainwindow.h"
24
25 #include <vtkUnstructuredGridWriter.h>
26 #include <vtkGeometryFilter.h>
27 #include <vtkXMLPolyDataWriter.h>
28 #include <vtkDelaunay3D.h>
29 #include <vtkHull.h>
30 #include <vtkTriangleFilter.h>
31
32 PolyMolecule::PolyMolecule()
33 {
34 m_AllPositive = false;
35 m_PolyMesh = NULL;
36 m_SplitCell1 = NULL;
37 m_SplitCell2 = NULL;
38 }
39
40 PolyMolecule::PolyMolecule(PolyMesh *poly_mesh, int i_pcell)
41 {
42 m_PCell = i_pcell;
43 QVector<QList<int> > faces(poly_mesh->numFacesOfPCell(i_pcell));
44 for (int i = 0; i < poly_mesh->numFacesOfPCell(i_pcell); ++i) {
45 int i_face = poly_mesh->pcell2Face(i_pcell, i);
46 m_FacesInPMesh.append(i_face);
47 int num_nodes = poly_mesh->numNodes(i_face);
48 for (int j = 0; j < num_nodes; ++j) {
49 int node = poly_mesh->nodeIndex(i_face, j);
50 if (poly_mesh->owner(i_face) == i_pcell) {
51 faces[i].append(node);
52 } else {
53 faces[i].prepend(node);
54 }
55 }
56 }
57 m_AllPositive = false;
58 m_SplitCell1 = NULL;
59 m_SplitCell2 = NULL;
60 init(poly_mesh, faces);
61 }
62
63 void PolyMolecule::writeVtkFile(QString file_name)
64 {
65 EG_VTKSP(vtkUnstructuredGrid, grid);
66 allocateGrid(grid, m_Faces.size(), m_Nodes.size());
67
68 EG_VTKSP(vtkDoubleArray, normals);
69 normals->SetNumberOfComponents(3);
70 normals->SetNumberOfTuples(m_Nodes.size());
71 normals->SetName("N");
72
73 EG_VTKSP(vtkDoubleArray, badness);
74 badness->SetNumberOfComponents(1);
75 badness->SetNumberOfTuples(m_Nodes.size());
76 badness->SetName("badness");
77
78 for (int i = 0; i < m_Nodes.size(); ++i) {
79 vec3_t x = getXNode(i);
80 grid->GetPoints()->SetPoint(i, x.data());
81 normals->SetTuple(i, m_NodeNormals[i].data());
82 badness->SetValue(i, double(m_N2BadFaces[i].size()));
83 }
84 grid->GetPointData()->AddArray(normals);
85 grid->GetPointData()->AddArray(badness);
86
87 for (int i = 0; i < m_Faces.size(); ++i) {
88 QVector<vtkIdType> pts(m_Faces[i].size());
89 for (int j = 0; j < m_Faces[i].size(); ++j) {
90 pts[j] = m_Faces[i][j];
91 }
92 grid->InsertNextCell(VTK_POLYGON, m_Faces[i].size(), pts.data());
93 }
94 EG_VTKSP(vtkXMLUnstructuredGridWriter, vtu);
95 file_name = GuiMainWindow::pointer()->getCwd() + "/" + file_name + ".vtu";
96 vtu->SetFileName(qPrintable(file_name));
97 vtu->SetDataModeToBinary();
98 vtu->SetInputData(grid);
99 vtu->Write();
100 }
101
102 void PolyMolecule::createPolyData(vtkPolyData *poly_data)
103 {
104 int N = m_Nodes.size();
105 EG_VTKSP(vtkPoints, points);
106 points->SetNumberOfPoints(N);
107 for (int i = 0; i < N; ++i) {
108 vec3_t x = getXNode(i);
109 points->SetPoint(i, x.data());
110 }
111 poly_data->Allocate(m_Faces.size());
112 poly_data->SetPoints(points);
113 for (int i = 0; i < m_Faces.size(); ++i) {
114 QVector<vtkIdType> pts(m_Faces[i].size());
115 for (int j = 0; j < m_Faces[i].size(); ++j) {
116 pts[j] = m_Faces[i][j];
117 }
118 if (pts.size() == 3) {
119 poly_data->InsertNextCell(VTK_TRIANGLE, pts.size(), pts.data());
120 } else if (pts.size() == 4) {
121 poly_data->InsertNextCell(VTK_QUAD, pts.size(), pts.data());
122 } else {
123 poly_data->InsertNextCell(VTK_POLYGON, pts.size(), pts.data());
124 }
125 }
126 }
127
128 void PolyMolecule::computeCentreOfGravity()
129 {
130 m_AllPositive = true;
131 // first guess of centre of gravity
132 vec3_t x_centre(0,0,0);
133 for (int i = 0; i < m_Faces.size(); ++i) {
134 x_centre += getXFace(i);
135 }
136 m_N2BadFaces.resize(m_Nodes.size());
137 for (int i = 0; i < m_Nodes.size(); ++i) {
138 m_N2BadFaces[i].clear();
139 }
140 x_centre *= 1.0/m_Faces.size();
141 for (int iter = 0; iter < 2; ++iter) {
142 double V_total = 0;
143 m_CentreOfGravity = vec3_t(0,0,0);
144 m_MaxPyramidVolume = -1e99;
145 m_MinPyramidVolume = 1e99;
146 for (int i = 0; i < m_Faces.size(); ++i) {
147 double V_pyramid = 0;
148 vec3_t x_face = getXFace(i);
149 QVector<vec3_t> x(m_Faces[i].size() + 1);
150 for (int j = 0; j < m_Faces[i].size(); ++j) {
151 x[j] = getXNode(m_Faces[i][j]);
152 }
153 x.last() = x.first();
154 for (int j = 0; j < m_Faces[i].size(); ++j) {
155 double V = GeometryTools::tetraVol(x_face, x[j+1], x[j], x_centre, true);
156 if (V <= 0) {
157 m_AllPositive = false;
158 }
159 V_total += V;
160 V_pyramid += V;
161 m_CentreOfGravity += 0.25*V*(x_face + x[j+1] + x[j] + x_centre);
162 }
163 if (V_pyramid <= 0) {
164 for (int j = 0; j < m_Faces[i].size(); ++j) {
165 m_N2BadFaces[m_Faces[i][j]].insert(i);
166 }
167 }
168 m_MaxPyramidVolume = max(V_pyramid, m_MaxPyramidVolume);
169 m_MinPyramidVolume = min(V_pyramid, m_MinPyramidVolume);
170 }
171 m_CentreOfGravity *= 1.0/V_total;
172 x_centre = m_CentreOfGravity;
173 }
174 }
175
176 void PolyMolecule::buildNode2Node()
177 {
178 m_N2N.clear();
179 m_N2N.resize(m_Nodes.size());
180 for (int i = 0; i < m_Faces.size(); ++i) {
181 for (int j = 0; j < m_Faces[i].size() - 1; ++j) {
182 m_N2N[m_Faces[i][j]].insert(m_Faces[i][j+1]);
183 m_N2N[m_Faces[i][j+1]].insert(m_Faces[i][j]);
184 }
185 }
186 }
187
188 void PolyMolecule::buildFace2Face()
189 {
190 m_F2F.clear();
191 m_F2F.resize(m_Faces.size());
192 for (int face = 0; face < m_Faces.size(); ++face) {
193 for (int i = 0; i < m_Faces[face].size(); ++i) {
194 int j = i + 1;
195 if (j >= m_Faces[face].size()) {
196 j = 0;
197 }
198 edge_t e = getEdge(m_Faces[face][i], m_Faces[face][j]);
199 if (e.face1 != face) {
200 m_F2F[face].insert(e.face1);
201 }
202 if (e.face2 != face) {
203 m_F2F[face].insert(e.face2);
204 }
205 }
206 }
207 }
208
209 void PolyMolecule::computeNormals()
210 {
211 m_FaceNormals.fill(vec3_t(0,0,0), m_Faces.size());
212 m_NodeNormals.fill(vec3_t(0,0,0), m_Nodes.size());
213 for (int i = 0; i < m_Faces.size(); ++i) {
214 QVector<vec3_t> x(m_Faces[i].size() + 1);
215 vec3_t x_face(0,0,0);
216 for (int j = 0; j < m_Faces[i].size(); ++j) {
217 x[j] = getXNode(m_Faces[i][j]);
218 x_face += x[j];
219 }
220 x.last() = x.first();
221 x_face *= 1.0/m_Faces[i].size();
222 m_FaceNormals[i] = vec3_t(0,0,0);
223 for (int j = 0; j < m_Faces[i].size(); ++j) {
224 m_FaceNormals[i] += triNormal(x_face, x[j], x[j+1]);
225 }
226 for (int j = 0; j < m_Faces[i].size(); ++j) {
227 vec3_t n = m_FaceNormals[i];
228 //n.normalise();
229 m_NodeNormals[m_Faces[i][j]] += n;
230 }
231 }
232 for (int i = 0; i < m_Nodes.size(); ++i) {
233 m_NodeNormals[i].normalise();
234 }
235 }
236
237 vec3_t PolyMolecule::getXFace(int face)
238 {
239 vec3_t x(0,0,0);
240 foreach (int i, m_Faces[face]) {
241 x += getXNode(i);
242 }
243 if (m_Nodes.size() == 0) {
244 EG_BUG;
245 }
246 return (1.0/m_Faces[face].size())*x;
247 }
248
249 void PolyMolecule::smooth(bool delaunay, bool write)
250 {
251 computeNormals();
252 if (write) {
253 writeVtkFile("bad_cell");
254 }
255
256 // collect all "healthy" cells in the neighbourhood
257 // they will be checked for volume inversion caused by the smoothing of this cell
258 //
259 QSet<int> check_cells;
260 check_cells.insert(m_PCell);
261 for (int i = 0; i < m_PolyMesh->numFacesOfPCell(m_PCell); ++i) {
262 int i_face = m_PolyMesh->pcell2Face(m_PCell, i);
263 int i_cell = m_PolyMesh->owner(i_face);
264 PolyMolecule pm(m_PolyMesh, i_cell);
265 if (pm.minPyramidVolume() > 0) {
266 check_cells.insert(i_cell);
267 }
268 i_cell = m_PolyMesh->neighbour(i_face);
269 if (i_cell != -1) {
270 PolyMolecule pm(m_PolyMesh, i_cell);
271 if (pm.minPyramidVolume() > 0) {
272 check_cells.insert(i_cell);
273 }
274 }
275 }
276
277 // VTK pipeline to get convex hull
278 // The cell will be "inflated" into this convex hull.
279 //
280 QVector<QVector<vec3_t> > convex_hull;
281 {
282 EG_VTKSP(vtkPolyData, poly_data);
283 EG_VTKSP(vtkPoints, points);
284 points->SetNumberOfPoints(m_Nodes.size());
285 for (vtkIdType id_node = 0; id_node < m_Nodes.size(); ++id_node) {
286 vec3_t x = getXNode(id_node);
287 points->SetPoint(id_node, x.data());
288 }
289 poly_data->Allocate(m_Faces.size());
290 poly_data->SetPoints(points);
291 for (int i = 0; i < m_Faces.size(); ++i) {
292 QVector<vtkIdType> pts(m_Faces[i].size());
293 for (int j = 0; j < m_Faces[i].size(); ++j) {
294 pts[j] = m_Faces[i][j];
295 }
296 poly_data->InsertNextCell(VTK_POLYGON, pts.size(), pts.data());
297 }
298 EG_VTKSP(vtkGeometryFilter, surface);
299 if (delaunay) {
300 EG_VTKSP(vtkDelaunay3D, delaunay);
301 delaunay->SetOffset(100);
302 delaunay->SetInputData(poly_data);
303 surface->SetInputConnection(delaunay->GetOutputPort());
304 } else {
305 EG_VTKSP(vtkHull, hull);
306 hull->SetInputData(poly_data);
307 hull->AddRecursiveSpherePlanes(5);
308 EG_VTKSP(vtkTriangleFilter, tri);
309 tri->SetInputConnection(hull->GetOutputPort());
310 surface->SetInputConnection(tri->GetOutputPort());
311 }
312 surface->Update();
313 if (surface->GetOutput()->GetNumberOfCells() < 4) {
314 return;
315 }
316 EG_VTKSP(vtkXMLPolyDataWriter, vtp);
317 if (write) {
318 QString file_name = GuiMainWindow::pointer()->getCwd() + "/" + "convex_hull.vtp";
319 vtp->SetFileName(qPrintable(file_name));
320 vtp->SetDataModeToAscii();
321 vtp->SetInputConnection(surface->GetOutputPort());
322 vtp->Write();
323 }
324 convex_hull.fill(QVector<vec3_t>(3), surface->GetOutput()->GetNumberOfCells());
325 for (vtkIdType id_tri = 0; id_tri < surface->GetOutput()->GetNumberOfCells(); ++id_tri) {
326 EG_GET_CELL(id_tri, surface->GetOutput());
327 if (type_cell != VTK_TRIANGLE) {
328 EG_BUG;
329 }
330 for (int i = 0; i < 3; ++i) {
331 vec3_t x;
332 surface->GetOutput()->GetPoint(pts[i], x.data());
333 convex_hull[id_tri][i] = x;
334 }
335 }
336 }
337
338 // Save old nodes for under-relaxations
339 //
340 QVector<vec3_t> x_old(m_Nodes.size());
341 for (int i = 0; i < m_Nodes.size(); ++i) {
342 x_old[i] = getXNode(i);
343 }
344
345 // "inflate" cell
346 //
347 for (int iter = 1; iter <= 100; ++iter) {
348 for (int i = 0; i < m_Nodes.size(); ++i) {
349 vec3_t x = getXNode(i);
350 double L = (x - m_CentreOfGravity).abs();
351
352 vec3_t Dx_n = 0.1*L*m_NodeNormals[i];
353 vec3_t Dx_lp(0,0,0);
354 foreach (int neigh, m_N2N[i]) {
355 Dx_lp += 0.1*(getXNode(neigh) - x);
356 }
357 //Dx_lp -= (Dx_lp * m_NodeNormals[i])*m_NodeNormals[i];
358 vec3_t Dx = Dx_n + Dx_lp;
359
360 for (int j = 0; j < convex_hull.size(); ++j) {
361 vec3_t np = triNormal(convex_hull[j][0], convex_hull[j][1], convex_hull[j][2]);
362 np.normalise();
363 double d = GeometryTools::intersection(x, Dx, convex_hull[j][0], np);
364 if (d > -1e-3 && d < 1) {
365 Dx -= (Dx*np)*np;
366 }
367 }
368
369 for (int j = 0; j < convex_hull.size(); ++j) {
370 vec3_t np = triNormal(convex_hull[j][0], convex_hull[j][1], convex_hull[j][2]);
371 np.normalise();
372 double d = GeometryTools::intersection(x, Dx, convex_hull[j][0], np);
373 if (d > -1e-3 && d < 1) {
374 Dx *= 0;
375 }
376 }
377
378 x += Dx;
379 m_PolyMesh->setNodeVector(m_Nodes[i], x);
380 }
381 QString file_name;
382 file_name.setNum(iter);
383 file_name = file_name.rightJustified(4, '0');
384 file_name = "improved_cell." + file_name;
385 if (write) writeVtkFile(file_name);
386 computeCentreOfGravity();
387 computeNormals();
388 }
389
390 // get new state of the cell for under relaxation
391 //
392 QVector<vec3_t> x_new(m_Nodes.size());
393 for (int i = 0; i < m_Nodes.size(); ++i) {
394 x_new[i] = getXNode(i);
395 }
396
397 // check all formerly "healthy" cells and under relax as rquired
398 //
399 bool done = false;
400 double relax = 1;
401 while (!done) {
402 done = true;
403 foreach (int i_pcell, check_cells) {
404 PolyMolecule pm(m_PolyMesh, i_pcell);
405 if (!pm.allPositive()) {
406 done = false;
407 break;
408 }
409 }
410 if (!done) {
411 relax -= 0.01;
412 if (relax < 1e-3) {
413 //relax = 1.0;
414 done = true;
415 }
416 for (int i = 0; i < m_Nodes.size(); ++i) {
417 m_PolyMesh->setNodeVector(m_Nodes[i], x_old[i] + relax*(x_new[i] - x_old[i]));
418 }
419 }
420 }
421
422
423 if (write) writeVtkFile("improved_cell");
424 }
425
426 PolyMolecule::edge_t PolyMolecule::getEdge(int node1, int node2)
427 {
428 edge_t e;
429 e.node1 = node1;
430 e.node2 = node2;
431 e.face1 = -1;
432 e.face2 = -1;
433 for (int i = 0; i < m_Faces.size(); ++i) {
434 for (int j = 0; j < m_Faces[i].size(); ++j) {
435 int n1 = m_Faces[i][j];
436 int n2 = m_Faces[i][0];
437 if (j < m_Faces[i].size() - 1) {
438 n2 = m_Faces[i][j+1];
439 }
440 if (n1 == node1 && n2 == node2) {
441 e.face1 = i;
442 }
443 if (n1 == node2 && n2 == node1) {
444 e.face2 = i;
445 }
446 }
447 }
448 if (e.face1 == -1 || e.face2 == -1) {
449 EG_BUG;
450 }
451 return e;
452 }
453
454 QList<PolyMolecule::edge_t> PolyMolecule::findConcaveEdges()
455 {
456 if (m_PCell == 32) {
457 cout << "break" << endl;
458 }
459 QList<edge_t> concave_edges;
460 for (int node1 = 0; node1 < m_Nodes.size(); ++node1) {
461 foreach (int node2, m_N2N[node1]) {
462 if (node2 > node1) {
463 edge_t e = getEdge(node1, node2);
464 vec3_t x1 = getXNode(node1);
465 vec3_t x2 = getXNode(node2);
466 vec3_t xe = 0.5*(getXNode(node1) + getXNode(node2));
467 vec3_t xf1 = getXFace(e.face1);
468 vec3_t xf2 = getXFace(e.face2);
469 vec3_t xc = m_CentreOfGravity;
470 vec3_t xn = 0.5*(xf1 + xf2);
471 double Lnorm = (xf2 - xf1).abs();
472 double L = ((xn - xc).abs() - (xe - xc).abs())/Lnorm;
473 if (L > 0.01) {
474 concave_edges.append(e);
475 }
476 }
477 }
478 }
479 return concave_edges;
480 }
481
482 void PolyMolecule::split(bool write)
483 {
484 if (write) {
485 writeVtkFile("concave_cell");
486 }
487 buildNode2Node();
488 buildFace2Face();
489 QList<edge_t> concave_edges = findConcaveEdges();
490 if (concave_edges.size() == 0) {
491 return;
492 }
493 QList<int> seed_faces;
494 foreach (edge_t e, concave_edges) {
495 if (!seed_faces.contains(e.face1)) {
496 seed_faces.append(e.face1);
497 }
498 if (!seed_faces.contains(e.face2)) {
499 seed_faces.append(e.face2);
500 }
501 }
502 if (seed_faces.size() != concave_edges.size() + 1) {
503 //EG_BUG;
504 }
505 QVector<int> old2new(m_Faces.size());
506 int num_faces = m_Faces.size();
507 for (int face = 0; face < m_Faces.size(); ++face) {
508 old2new[face] = seed_faces.indexOf(face);
509 if (old2new[face] != -1) {
510 --num_faces;
511 }
512 }
513 while (num_faces > 0) {
514 for (int face1 = 0; face1 < m_Faces.size(); ++face1) {
515 if (old2new[face1] != -1) {
516 foreach (int face2, m_F2F[face1]) {
517 if (old2new[face2] == -1) {
518 old2new[face2] = old2new[face1];
519 --num_faces;
520 }
521 }
522 }
523 }
524 }
525 QVector<int> new_pcell_index(seed_faces.size());
526 new_pcell_index[0] = m_PCell;
527 for (int i = 1; i < seed_faces.size(); ++i) {
528 new_pcell_index[i] = m_PolyMesh->numCells() + i - 1;
529 }
530 QVector<PolyMesh::face_t> new_faces(concave_edges.size());
531 QVector<vec3_t> new_cell_centre(seed_faces.size(), vec3_t(0,0,0));
532 {
533 QVector<int> count(seed_faces.size(), 0);
534 for (int face = 0; face < m_Faces.size(); ++face) {
535 new_cell_centre[old2new[face]] += getXFace(face);
536 ++count[old2new[face]];
537 }
538 for (int i = 0; i < new_cell_centre.size(); ++i) {
539 if (count[i] == 0) {
540 EG_BUG;
541 }
542 new_cell_centre[i] *= 1.0/count[i];
543 }
544 }
545 for (int i = 0; i < concave_edges.size(); ++i) {
546 edge_t e = concave_edges[i];
547 int start = e.node1;
548 int first = e.node2;
549 new_faces[i].node.append(m_Nodes[e.node1]);
550 int o2n1 = old2new[e.face1];
551 int o2n2 = old2new[e.face2];
552 new_faces[i].ref_vec = new_cell_centre[o2n2] - new_cell_centre[o2n1];
553 new_faces[i].owner = new_pcell_index[o2n1];
554 new_faces[i].neighbour = new_pcell_index[o2n2];
555 new_faces[i].bc = 0;
556 bool reversed = false;
557 while (e.node2 != start) {
558 if (reversed) {
559 new_faces[i].node.prepend(m_Nodes[e.node2]);
560 } else {
561 new_faces[i].node.append(m_Nodes[e.node2]);
562 }
563 bool found = false;
564 foreach (int n, m_N2N[e.node2]) {
565 if (n != e.node1) {
566 edge_t ne = getEdge(e.node2, n);
567 if (old2new[ne.face1] == o2n1 && old2new[ne.face2] == o2n2) {
568 e = ne;
569 found = true;
570 break;
571 }
572 }
573 }
574
575 // reached a dead-end here ...
576 if (!found && !reversed) {
577 reversed = true;
578 swap(o2n1, o2n2);
579 swap(start, first);
580 foreach (int n, m_N2N[first]) {
581 if (n != start) {
582 edge_t ne = getEdge(first, n);
583 if (old2new[ne.face1] == o2n1 && old2new[ne.face2] == o2n2) {
584 e = ne;
585 found = true;
586 break;
587 }
588 }
589 }
590 }
591 if (!found) {
592 break;
593 }
594 }
595 if (new_faces[i].owner > new_faces[i].neighbour && new_faces[i].neighbour != -1) {
596 swap(new_faces[i].owner, new_faces[i].neighbour);
597 new_faces[i].ref_vec *= -1;
598 invertQContainer(new_faces[i].node);
599 }
600 }
601 for (int face = 0; face < m_Faces.size(); ++face) {
602 int pface = m_FacesInPMesh[face];
603 //m_PolyMesh->m_Faces[pface].bc = 1;
604 if (m_PolyMesh->m_Faces[pface].owner == m_PCell) {
605 m_PolyMesh->m_Faces[pface].owner = new_pcell_index[old2new[face]];
606 }
607 if (m_PolyMesh->m_Faces[pface].neighbour == m_PCell) {
608 m_PolyMesh->m_Faces[pface].neighbour = new_pcell_index[old2new[face]];
609 }
610 if (m_PolyMesh->m_Faces[pface].owner > m_PolyMesh->m_Faces[pface].neighbour && m_PolyMesh->m_Faces[pface].neighbour != -1) {
611 swap(m_PolyMesh->m_Faces[pface].owner, m_PolyMesh->m_Faces[pface].neighbour);
612 m_PolyMesh->m_Faces[pface].ref_vec *= -1;
613 invertQContainer(m_PolyMesh->m_Faces[pface].node);
614 }
615 }
616 foreach (PolyMesh::face_t f, new_faces) {
617 if (f.owner == -1) {
618 EG_BUG;
619 }
620 if (f.neighbour == -1) {
621 EG_BUG;
622 }
623 if (f.owner > f.neighbour) {
624 swap(f.owner, f.neighbour);
625 invertQContainer(f.node);
626 }
627 //f.bc = 2;
628 m_PolyMesh->m_Faces.append(f);
629 }
630 m_PolyMesh->m_NumPolyCells += seed_faces.size() - 1;
631 }
632
633 void PolyMolecule::updatePMesh()
634 {
635 int offset = 0;
636 foreach (int face, m_FacesInPMesh) {
637 m_PolyMesh->m_Faces.removeAt(face - offset);
638 ++offset;
639 }
640 }
641
642 void PolyMolecule::optimise(bool write)
643 {
644 smooth(true, write);
645 if (m_MinPyramidVolume <= 0) {
646 smooth(false, write);
647 }
648 }
649
650 void PolyMolecule::updateFace(int face, int new_cell_index)
651 {
652 PolyMesh::face_t f = m_PolyMesh->m_Faces[m_FacesInPMesh[face]];
653 if (f.owner == m_PCell) {
654 f.owner = new_cell_index;
655 }
656 if (f.neighbour == m_PCell) {
657 f.neighbour = new_cell_index;
658 }
659 if (f.owner > f.neighbour && f.neighbour != -1) {
660 swap(f.owner, f.neighbour);
661 invertQContainer(f.node);
662 f.ref_vec *= -1;
663 }
664 m_PolyMesh->m_Faces[m_FacesInPMesh[face]] = f;
665 }
666
667 void PolyMolecule::centreSplit()
668 {
669 buildNode2Node();
670 buildFace2Face();
671 QList<edge_t> concave_edges = findConcaveEdges();
672 if (concave_edges.size() == 0) {
673 return;
674 }
675
676 // find common nodes in all edges and return if none exists
677 QSet<int> common_nodes;
678 common_nodes.insert(concave_edges.first().node1);
679 common_nodes.insert(concave_edges.first().node2);
680 foreach (edge_t e, concave_edges) {
681 QSet<int> edge_nodes;
682 edge_nodes.insert(e.node1);
683 edge_nodes.insert(e.node2);
684 common_nodes.intersect(edge_nodes);
685 }
686 if (common_nodes.size() == 0) {
687 return;
688 }
689
690 // start point of straight for divide point computation
691 vec3_t x1;
692 {
693 int count = 0;
694 foreach(int node, common_nodes) {
695 x1 += getXNode(node);
696 ++count;
697 }
698 x1 *= 1.0/count;
699 }
700
701 // end point of straight for divide point computation
702 computeNormals();
703 vec3_t x2;
704 {
705 vec3_t n(0,0,0);
706 QList<int> adjacent_faces;
707 foreach (edge_t e, concave_edges) {
708 if (!adjacent_faces.contains(e.face1)) {
709 adjacent_faces.append(e.face1);
710 }
711 if (!adjacent_faces.contains(e.face2)) {
712 adjacent_faces.append(e.face2);
713 }
714 }
715 foreach (int face, adjacent_faces) {
716 n += m_FaceNormals[face];
717 }
718 n.normalise();
719 bool x2_set = false;
720 double max_dist = 0;
721 for (int node = 0; node < m_Nodes.size(); ++node) {
722 vec3_t x = getXNode(node);
723 double dist = (x1 - x)*n;
724 if (dist > max_dist) {
725 x2 = x1 - dist*n;
726 max_dist = dist;
727 x2_set = true;
728 }
729 }
730 if (!x2_set) {
731 EG_BUG;
732 }
733 }
734
735 // tip point for all polygonial pyramids
736 vec3_t x_tip = 0.5*(x1 + x2);
737
738 // vector with new cell indices
739 QVector<int> cell_index(m_Faces.size());
740 cell_index[0] = m_PCell;
741 for (int i = 1; i < cell_index.size(); ++i) {
742 cell_index[i] = m_PolyMesh->m_NumPolyCells++;
743 }
744
745 // new node index
746 int new_node = m_PolyMesh->m_Points.size();
747 m_PolyMesh->m_Points.append(x_tip);
748
749 // re-map existing faces
750 for (int face = 0; face < m_FacesInPMesh.size(); ++face) {
751 updateFace(face, cell_index[face]);
752 }
753
754 // create new faces
755 for (int face = 0; face < m_Faces.size(); ++face) {
756 QList<int> nodes = m_Faces[face];
757 nodes.append(nodes.first());
758 for (int i_node = 0; i_node < nodes.size() - 1; ++i_node) {
759 int node1 = nodes[i_node];
760 int node2 = nodes[i_node + 1];
761 edge_t e = getEdge(node1, node2);
762 PolyMesh::face_t f;
763 if (e.face1 == face) {
764 f.neighbour = cell_index[e.face2];
765 } else {
766 f.neighbour = cell_index[e.face1];
767 }
768 if (f.owner < f.neighbour) {
769 f.bc = 0;
770 f.node.append(m_Nodes[node2]);
771 f.node.append(m_Nodes[node1]);
772 f.node.append(new_node);
773 f.owner = cell_index[face];
774 vec3_t x1 = m_PolyMesh->m_Points[f.node[0]];
775 vec3_t x2 = m_PolyMesh->m_Points[f.node[1]];
776 vec3_t x3 = m_PolyMesh->m_Points[f.node[2]];
777 f.ref_vec = GeometryTools::triNormal(x1, x2, x3);
778 /*
779 if (f.owner > f.neighbour) {
780 swap(f.owner, f.neighbour);
781 f.ref_vec *= -1;
782 invertQContainer(f.node);
783 }
784 */
785 m_PolyMesh->m_Faces.append(f);
786 }
787 }
788 }
789 }
790
791
This is a mirror repository of the new location for engrid@github