?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
various improvements to speed up surface meshing of BRL-CAD geometries
[engrid.git] / src / libengrid / polymolecule.cpp
blob0798f4337f87e13119f53edce9b78dc921c21aeb
1 //
2 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3 // + +
4 // + This file is part of enGrid. +
5 // + +
6 // + Copyright 2008-2013 enGits GmbH +
7 // + +
8 // + enGrid is free software: you can redistribute it and/or modify +
9 // + it under the terms of the GNU General Public License as published by +
10 // + the Free Software Foundation, either version 3 of the License, or +
11 // + (at your option) any later version. +
12 // + +
13 // + enGrid is distributed in the hope that it will be useful, +
14 // + but WITHOUT ANY WARRANTY; without even the implied warranty of +
15 // + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +
16 // + GNU General Public License for more details. +
17 // + +
18 // + You should have received a copy of the GNU General Public License +
19 // + along with enGrid. If not, see <http://www.gnu.org/licenses/>. +
20 // + +
21 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
22 //
23
24 #include "polymolecule.h"
25 #include "guimainwindow.h"
26
27 #include <vtkUnstructuredGridWriter.h>
28 #include <vtkGeometryFilter.h>
29 #include <vtkXMLPolyDataWriter.h>
30 #include <vtkDelaunay3D.h>
31 #include <vtkHull.h>
32 #include <vtkTriangleFilter.h>
33
34 PolyMolecule::PolyMolecule()
35 {
36 m_AllPositive = false;
37 m_PolyMesh = NULL;
38 m_SplitCell1 = NULL;
39 m_SplitCell2 = NULL;
40 }
41
42 PolyMolecule::PolyMolecule(PolyMesh *poly_mesh, int i_pcell)
43 {
44 m_PCell = i_pcell;
45 QVector<QList<int> > faces(poly_mesh->numFacesOfPCell(i_pcell));
46 for (int i = 0; i < poly_mesh->numFacesOfPCell(i_pcell); ++i) {
47 int i_face = poly_mesh->pcell2Face(i_pcell, i);
48 int num_nodes = poly_mesh->numNodes(i_face);
49 for (int j = 0; j < num_nodes; ++j) {
50 int node = poly_mesh->nodeIndex(i_face, j);
51 if (poly_mesh->owner(i_face) == i_pcell) {
52 faces[i].append(node);
53 } else {
54 faces[i].prepend(node);
55 }
56 }
57 }
58 m_AllPositive = false;
59 m_SplitCell1 = NULL;
60 m_SplitCell2 = NULL;
61 init(poly_mesh, faces);
62 }
63
64 void PolyMolecule::writeVtkFile(QString file_name)
65 {
66 EG_VTKSP(vtkUnstructuredGrid, grid);
67 allocateGrid(grid, m_Faces.size(), m_Nodes.size());
68
69 EG_VTKSP(vtkDoubleArray, normals);
70 normals->SetNumberOfComponents(3);
71 normals->SetNumberOfTuples(m_Nodes.size());
72 normals->SetName("N");
73
74 EG_VTKSP(vtkDoubleArray, badness);
75 badness->SetNumberOfComponents(1);
76 badness->SetNumberOfTuples(m_Nodes.size());
77 badness->SetName("badness");
78
79 for (int i = 0; i < m_Nodes.size(); ++i) {
80 vec3_t x = getXNode(i);
81 grid->GetPoints()->SetPoint(i, x.data());
82 normals->SetTuple(i, m_NodeNormals[i].data());
83 badness->SetValue(i, double(m_N2BadFaces[i].size()));
84 }
85 grid->GetPointData()->AddArray(normals);
86 grid->GetPointData()->AddArray(badness);
87
88 for (int i = 0; i < m_Faces.size(); ++i) {
89 QVector<vtkIdType> pts(m_Faces[i].size());
90 for (int j = 0; j < m_Faces[i].size(); ++j) {
91 pts[j] = m_Faces[i][j];
92 }
93 grid->InsertNextCell(VTK_POLYGON, m_Faces[i].size(), pts.data());
94 }
95 EG_VTKSP(vtkXMLUnstructuredGridWriter, vtu);
96 file_name = GuiMainWindow::pointer()->getCwd() + "/" + file_name + ".vtu";
97 vtu->SetFileName(qPrintable(file_name));
98 vtu->SetDataModeToBinary();
99 vtu->SetInput(grid);
100 vtu->Write();
101 }
102
103 void PolyMolecule::createPolyData(vtkPolyData *poly_data)
104 {
105 int N = m_Nodes.size();
106 EG_VTKSP(vtkPoints, points);
107 points->SetNumberOfPoints(N);
108 for (int i = 0; i < N; ++i) {
109 vec3_t x = getXNode(i);
110 points->SetPoint(i, x.data());
111 }
112 poly_data->Allocate(m_Faces.size());
113 poly_data->SetPoints(points);
114 for (int i = 0; i < m_Faces.size(); ++i) {
115 QVector<vtkIdType> pts(m_Faces[i].size());
116 for (int j = 0; j < m_Faces[i].size(); ++j) {
117 pts[j] = m_Faces[i][j];
118 }
119 if (pts.size() == 3) {
120 poly_data->InsertNextCell(VTK_TRIANGLE, pts.size(), pts.data());
121 } else if (pts.size() == 4) {
122 poly_data->InsertNextCell(VTK_QUAD, pts.size(), pts.data());
123 } else {
124 poly_data->InsertNextCell(VTK_POLYGON, pts.size(), pts.data());
125 }
126 }
127 }
128
129 void PolyMolecule::computeCentreOfGravity()
130 {
131 m_AllPositive = true;
132 // first guess of centre of gravity
133 vec3_t x_centre(0,0,0);
134 for (int i = 0; i < m_Faces.size(); ++i) {
135 x_centre += getXFace(i);
136 }
137 m_N2BadFaces.resize(m_Nodes.size());
138 for (int i = 0; i < m_Nodes.size(); ++i) {
139 m_N2BadFaces[i].clear();
140 }
141 x_centre *= 1.0/m_Faces.size();
142 for (int iter = 0; iter < 2; ++iter) {
143 double V_total = 0;
144 m_CentreOfGravity = vec3_t(0,0,0);
145 m_MaxPyramidVolume = -1e99;
146 m_MinPyramidVolume = 1e99;
147 for (int i = 0; i < m_Faces.size(); ++i) {
148 double V_pyramid = 0;
149 vec3_t x_face = getXFace(i);
150 QVector<vec3_t> x(m_Faces[i].size() + 1);
151 for (int j = 0; j < m_Faces[i].size(); ++j) {
152 x[j] = getXNode(m_Faces[i][j]);
153 }
154 x.last() = x.first();
155 for (int j = 0; j < m_Faces[i].size(); ++j) {
156 double V = GeometryTools::tetraVol(x_face, x[j+1], x[j], x_centre, true);
157 if (V <= 0) {
158 m_AllPositive = false;
159 }
160 V_total += V;
161 V_pyramid += V;
162 m_CentreOfGravity += 0.25*V*(x_face + x[j+1] + x[j] + x_centre);
163 }
164 if (V_pyramid <= 0) {
165 for (int j = 0; j < m_Faces[i].size(); ++j) {
166 m_N2BadFaces[m_Faces[i][j]].insert(i);
167 }
168 }
169 m_MaxPyramidVolume = max(V_pyramid, m_MaxPyramidVolume);
170 m_MinPyramidVolume = min(V_pyramid, m_MinPyramidVolume);
171 }
172 m_CentreOfGravity *= 1.0/V_total;
173 x_centre = m_CentreOfGravity;
174 }
175 }
176
177 void PolyMolecule::buildNode2Node()
178 {
179 m_N2N.clear();
180 m_N2N.resize(m_Nodes.size());
181 for (int i = 0; i < m_Faces.size(); ++i) {
182 for (int j = 0; j < m_Faces[i].size() - 1; ++j) {
183 m_N2N[m_Faces[i][j]].insert(m_Faces[i][j+1]);
184 m_N2N[m_Faces[i][j+1]].insert(m_Faces[i][j]);
185 }
186 }
187 }
188
189 void PolyMolecule::computeNormals()
190 {
191 m_FaceNormals.fill(vec3_t(0,0,0), m_Faces.size());
192 m_NodeNormals.fill(vec3_t(0,0,0), m_Nodes.size());
193 for (int i = 0; i < m_Faces.size(); ++i) {
194 QVector<vec3_t> x(m_Faces[i].size() + 1);
195 vec3_t x_face(0,0,0);
196 for (int j = 0; j < m_Faces[i].size(); ++j) {
197 x[j] = getXNode(m_Faces[i][j]);
198 x_face += x[j];
199 }
200 x.last() = x.first();
201 x_face *= 1.0/m_Faces[i].size();
202 m_FaceNormals[i] = vec3_t(0,0,0);
203 for (int j = 0; j < m_Faces[i].size(); ++j) {
204 m_FaceNormals[i] += triNormal(x_face, x[j], x[j+1]);
205 }
206 for (int j = 0; j < m_Faces[i].size(); ++j) {
207 vec3_t n = m_FaceNormals[i];
208 //n.normalise();
209 m_NodeNormals[m_Faces[i][j]] += n;
210 }
211 }
212 for (int i = 0; i < m_Nodes.size(); ++i) {
213 m_NodeNormals[i].normalise();
214 }
215 }
216
217 vec3_t PolyMolecule::getXFace(int face)
218 {
219 vec3_t x(0,0,0);
220 foreach (int i, m_Faces[face]) {
221 x += getXNode(i);
222 }
223 if (m_Nodes.size() == 0) {
224 EG_BUG;
225 }
226 return (1.0/m_Faces[face].size())*x;
227 }
228
229 void PolyMolecule::smooth(bool delaunay, bool write)
230 {
231 computeNormals();
232 if (write) writeVtkFile("bad_cell");
233
234 // collect all "healthy" cells in the neighbourhood
235 // they will be checked for volume inversion caused by the smoothing of this cell
236 //
237 QSet<int> check_cells;
238 check_cells.insert(m_PCell);
239 for (int i = 0; i < m_PolyMesh->numFacesOfPCell(m_PCell); ++i) {
240 int i_face = m_PolyMesh->pcell2Face(m_PCell, i);
241 int i_cell = m_PolyMesh->owner(i_face);
242 PolyMolecule pm(m_PolyMesh, i_cell);
243 if (pm.minPyramidVolume() > 0) {
244 check_cells.insert(i_cell);
245 }
246 i_cell = m_PolyMesh->neighbour(i_face);
247 if (i_cell != -1) {
248 PolyMolecule pm(m_PolyMesh, i_cell);
249 if (pm.minPyramidVolume() > 0) {
250 check_cells.insert(i_cell);
251 }
252 }
253 }
254
255 // VTK pipeline to get convex hull
256 // The cell will be "inflated" into this convex hull.
257 //
258 QVector<QVector<vec3_t> > convex_hull;
259 {
260 EG_VTKSP(vtkPolyData, poly_data);
261 EG_VTKSP(vtkPoints, points);
262 points->SetNumberOfPoints(m_Nodes.size());
263 for (vtkIdType id_node = 0; id_node < m_Nodes.size(); ++id_node) {
264 vec3_t x = getXNode(id_node);
265 points->SetPoint(id_node, x.data());
266 }
267 poly_data->Allocate(m_Faces.size());
268 poly_data->SetPoints(points);
269 for (int i = 0; i < m_Faces.size(); ++i) {
270 QVector<vtkIdType> pts(m_Faces[i].size());
271 for (int j = 0; j < m_Faces[i].size(); ++j) {
272 pts[j] = m_Faces[i][j];
273 }
274 poly_data->InsertNextCell(VTK_POLYGON, pts.size(), pts.data());
275 }
276 EG_VTKSP(vtkGeometryFilter, surface);
277 if (delaunay) {
278 EG_VTKSP(vtkDelaunay3D, delaunay);
279 delaunay->SetOffset(100);
280 delaunay->SetInput(poly_data);
281 surface->SetInput(delaunay->GetOutput());
282 } else {
283 EG_VTKSP(vtkHull, hull);
284 hull->SetInput(poly_data);
285 hull->AddRecursiveSpherePlanes(5);
286 EG_VTKSP(vtkTriangleFilter, tri);
287 tri->SetInput(hull->GetOutput());
288 surface->SetInput(tri->GetOutput());
289 }
290 surface->Update();
291 if (surface->GetOutput()->GetNumberOfCells() < 4) {
292 return;
293 }
294 EG_VTKSP(vtkXMLPolyDataWriter, vtp);
295 if (write) {
296 QString file_name = GuiMainWindow::pointer()->getCwd() + "/" + "convex_hull.vtp";
297 vtp->SetFileName(qPrintable(file_name));
298 vtp->SetDataModeToAscii();
299 vtp->SetInput(surface->GetOutput());
300 vtp->Write();
301 }
302 convex_hull.fill(QVector<vec3_t>(3), surface->GetOutput()->GetNumberOfCells());
303 for (vtkIdType id_tri = 0; id_tri < surface->GetOutput()->GetNumberOfCells(); ++id_tri) {
304 EG_GET_CELL(id_tri, surface->GetOutput());
305 if (type_cell != VTK_TRIANGLE) {
306 EG_BUG;
307 }
308 for (int i = 0; i < 3; ++i) {
309 vec3_t x;
310 surface->GetOutput()->GetPoint(pts[i], x.data());
311 convex_hull[id_tri][i] = x;
312 }
313 }
314 }
315
316 // Save old nodes for under-relaxations
317 //
318 QVector<vec3_t> x_old(m_Nodes.size());
319 for (int i = 0; i < m_Nodes.size(); ++i) {
320 x_old[i] = getXNode(i);
321 }
322
323 // "inflate" cell
324 //
325 for (int iter = 1; iter <= 100; ++iter) {
326 for (int i = 0; i < m_Nodes.size(); ++i) {
327 vec3_t x = getXNode(i);
328 double L = (x - m_CentreOfGravity).abs();
329
330 vec3_t Dx_n = 0.1*L*m_NodeNormals[i];
331 vec3_t Dx_lp(0,0,0);
332 foreach (int neigh, m_N2N[i]) {
333 Dx_lp += 0.1*(getXNode(neigh) - x);
334 }
335 //Dx_lp -= (Dx_lp * m_NodeNormals[i])*m_NodeNormals[i];
336 vec3_t Dx = Dx_n + Dx_lp;
337
338 for (int j = 0; j < convex_hull.size(); ++j) {
339 vec3_t np = triNormal(convex_hull[j][0], convex_hull[j][1], convex_hull[j][2]);
340 np.normalise();
341 double d = GeometryTools::intersection(x, Dx, convex_hull[j][0], np);
342 if (d > -1e-3 && d < 1) {
343 Dx -= (Dx*np)*np;
344 }
345 }
346
347 for (int j = 0; j < convex_hull.size(); ++j) {
348 vec3_t np = triNormal(convex_hull[j][0], convex_hull[j][1], convex_hull[j][2]);
349 np.normalise();
350 double d = GeometryTools::intersection(x, Dx, convex_hull[j][0], np);
351 if (d > -1e-3 && d < 1) {
352 Dx *= 0;
353 }
354 }
355
356 x += Dx;
357 m_PolyMesh->setNodeVector(m_Nodes[i], x);
358 }
359 QString file_name;
360 file_name.setNum(iter);
361 file_name = file_name.rightJustified(4, '0');
362 file_name = "improved_cell." + file_name;
363 if (write) writeVtkFile(file_name);
364 computeCentreOfGravity();
365 computeNormals();
366 }
367
368 // get new state of the cell for under relaxation
369 //
370 QVector<vec3_t> x_new(m_Nodes.size());
371 for (int i = 0; i < m_Nodes.size(); ++i) {
372 x_new[i] = getXNode(i);
373 }
374
375 // check all formerly "healthy" cells and under relax as rquired
376 //
377 bool done = false;
378 double relax = 1;
379 while (!done) {
380 done = true;
381 foreach (int i_pcell, check_cells) {
382 PolyMolecule pm(m_PolyMesh, i_pcell);
383 if (!pm.allPositive()) {
384 done = false;
385 break;
386 }
387 }
388 if (!done) {
389 relax -= 0.01;
390 if (relax < 1e-3) {
391 //relax = 1.0;
392 done = true;
393 }
394 for (int i = 0; i < m_Nodes.size(); ++i) {
395 m_PolyMesh->setNodeVector(m_Nodes[i], x_old[i] + relax*(x_new[i] - x_old[i]));
396 }
397 }
398 }
399
400
401 if (write) writeVtkFile("improved_cell");
402 }
403
404 void PolyMolecule::split(bool write)
405 {
406 if (write) writeVtkFile("concave_cell");
407 }
408
409 void PolyMolecule::fix(bool write)
410 {
411 /*
412 if (m_MinPyramidVolume/m_MaxPyramidVolume < -0.1) {
413 split(true);
414 EG_BUG;
415 }
416 */
417
418 smooth(true, write);
419 if (m_MinPyramidVolume <= 0) {
420 smooth(false, write);
421 }
422 //EG_BUG;
423 }
424
425
426
Clone of the official engrid repository