2 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
4 // + This file is part of enGrid. +
6 // + Copyright 2008,2009 Oliver Gloth +
8 // + enGrid is free software: you can redistribute it and/or modify +
9 // + it under the terms of the GNU General Public License as published by +
10 // + the Free Software Foundation, either version 3 of the License, or +
11 // + (at your option) any later version. +
13 // + enGrid is distributed in the hope that it will be useful, +
14 // + but WITHOUT ANY WARRANTY; without even the implied warranty of +
15 // + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +
16 // + GNU General Public License for more details. +
18 // + You should have received a copy of the GNU General Public License +
19 // + along with enGrid. If not, see <http://www.gnu.org/licenses/>. +
21 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
23 #include "createvolumemesh.h"
24 #include "deletetetras.h"
25 #include <vtkXMLUnstructuredGridWriter.h>
27 CreateVolumeMesh::CreateVolumeMesh()
33 void CreateVolumeMesh::setTraceCells(const QVector
<vtkIdType
> &cells
)
35 trace_cells
.resize(cells
.size());
36 qCopy(cells
.begin(), cells
.end(), trace_cells
.begin());
39 void CreateVolumeMesh::getTraceCells(QVector
<vtkIdType
> &cells
)
41 cells
.resize(trace_cells
.size());
42 qCopy(trace_cells
.begin(), trace_cells
.end(), cells
.begin());
45 void CreateVolumeMesh::prepare()
47 using namespace nglib
;
51 QVector
<vtkIdType
> cells
, nodes
;
52 QVector
<int> _cells
, _nodes
;
53 QVector
<QVector
< int > > c2c
;
54 QVector
<QSet
<int> > n2c
;
55 getAllCells(cells
, grid
);
56 createCellMapping(cells
, _cells
, grid
);
57 getNodesFromCells(cells
, nodes
, grid
);
58 createNodeMapping(nodes
, _nodes
, grid
);
59 createCellToCell(cells
, c2c
, grid
);
60 createNodeToCell(cells
, nodes
, _nodes
, n2c
, grid
);
61 QList
<QVector
<vtkIdType
> > ex_tri
;
62 QList
<vec3_t
> centres
;
69 foreach (vtkIdType id_cell
, cells
) {
70 vtkIdType type_cell
= grid
->GetCellType(id_cell
);
71 vtkIdType
*pts
, N_pts
;
72 grid
->GetCellPoints(id_cell
, N_pts
, pts
);
73 QVector
<vtkIdType
> T(3);
74 if (type_cell
== VTK_TRIANGLE
) {
75 if (findVolumeCell(grid
, id_cell
, _nodes
, cells
, _cells
, n2c
) == -1) {
82 } else if (type_cell
== VTK_QUAD
) {
83 if (findVolumeCell(grid
, id_cell
, _nodes
, cells
, _cells
, n2c
) == -1) {
95 } else if (type_cell
== VTK_WEDGE
) {
96 if (c2c
[id_cell
][0] == -1) {
103 if (c2c
[id_cell
][1] == -1) {
110 if (c2c
[id_cell
][2] == -1) {
121 for (int i
= 0; i
< 6; ++i
) grid
->GetPoint(pts
[i
], x
[i
].data());
122 vec3_t xc
= 0.25*(x
[0]+x
[1]+x
[3]+x
[4]);
125 conn1
.append(cellCentre(grid
, id_cell
));
126 conn2
.append(cellCentre(grid
, cells
[c2c
[id_cell
][2]]));
128 if (c2c
[id_cell
][3] == -1) {
139 for (int i
= 0; i
< 6; ++i
) grid
->GetPoint(pts
[i
], x
[i
].data());
140 vec3_t xc
= 0.25*(x
[4]+x
[1]+x
[2]+x
[5]);
143 conn1
.append(cellCentre(grid
, id_cell
));
144 conn2
.append(cellCentre(grid
, cells
[c2c
[id_cell
][3]]));
146 if (c2c
[id_cell
][4] == -1) {
157 for (int i
= 0; i
< 6; ++i
) grid
->GetPoint(pts
[i
], x
[i
].data());
158 vec3_t xc
= 0.25*(x
[0]+x
[3]+x
[2]+x
[5]);
161 conn1
.append(cellCentre(grid
, id_cell
));
162 conn2
.append(cellCentre(grid
, cells
[c2c
[id_cell
][4]]));
168 cout
<< "*********************************************************************" << endl
;
169 cout
<< "prism quads : " << N1
<< endl
;
170 cout
<< "prism triangles : " << N2
<< endl
;
171 cout
<< "stray quads : " << N3
<< endl
;
172 cout
<< "stray triangles : " << N4
<< endl
;
173 cout
<< "*********************************************************************" << endl
;
174 tri
.resize(ex_tri
.size());
175 qCopy(ex_tri
.begin(), ex_tri
.end(), tri
.begin());
176 add_to_ng
.fill(false, grid
->GetNumberOfPoints());
177 foreach (QVector
<vtkIdType
> T
, tri
) {
178 add_to_ng
[T
[0]] = true;
179 add_to_ng
[T
[1]] = true;
180 add_to_ng
[T
[2]] = true;
182 num_nodes_to_add
= 0;
184 for (vtkIdType id_node
= 0; id_node
< grid
->GetNumberOfPoints(); ++id_node
) {
185 if (add_to_ng
[id_node
]) {
193 QVector
<vtkIdType
> old2tri(grid
->GetNumberOfPoints(),-1);
195 for (vtkIdType id_node
= 0; id_node
< grid
->GetNumberOfPoints(); ++id_node
) {
196 if (add_to_ng
[id_node
]) {
197 old2tri
[id_node
] = N
;
201 EG_VTKSP(vtkUnstructuredGrid
,tri_grid
);
202 allocateGrid(tri_grid
, tri
.size(), N
);
203 for (vtkIdType id_node
= 0; id_node
< grid
->GetNumberOfPoints(); ++id_node
) {
204 if (add_to_ng
[id_node
]) {
206 grid
->GetPoint(id_node
, x
.data());
207 tri_grid
->GetPoints()->SetPoint(old2tri
[id_node
], x
.data());
208 copyNodeData(grid
, id_node
, tri_grid
, old2tri
[id_node
]);
211 foreach (QVector
<vtkIdType
> T
, tri
) {
213 pts
[0] = old2tri
[T
[0]];
214 pts
[1] = old2tri
[T
[1]];
215 pts
[2] = old2tri
[T
[2]];
216 tri_grid
->InsertNextCell(VTK_TRIANGLE
, 3, pts
);
218 EG_VTKSP(vtkXMLUnstructuredGridWriter
,vtu
);
219 vtu
->SetFileName("triangles.vtu");
220 vtu
->SetDataModeToBinary();
221 vtu
->SetInput(tri_grid
);
225 EG_VTKSP(vtkUnstructuredGrid
,centre_grid
);
226 allocateGrid(centre_grid
, centres
.size(), centres
.size());
229 foreach (vec3_t x
, centres
) {
230 centre_grid
->GetPoints()->SetPoint(N
, x
.data());
232 centre_grid
->InsertNextCell(VTK_VERTEX
, 1, pts
);
235 EG_VTKSP(vtkXMLUnstructuredGridWriter
,vtu
);
236 vtu
->SetFileName("centres.vtu");
237 vtu
->SetDataModeToBinary();
238 vtu
->SetInput(centre_grid
);
242 EG_VTKSP(vtkUnstructuredGrid
,conn_grid
);
243 allocateGrid(conn_grid
, conn1
.size(), 2*conn1
.size());
246 foreach (vec3_t x
, conn1
) {
247 conn_grid
->GetPoints()->SetPoint(N
, x
.data());
250 foreach (vec3_t x
, conn2
) {
251 conn_grid
->GetPoints()->SetPoint(N
, x
.data());
254 for (int i
= 0; i
< conn1
.size(); ++i
) {
256 pts
[1] = i
+conn1
.size();
257 conn_grid
->InsertNextCell(VTK_LINE
, 2, pts
);
259 EG_VTKSP(vtkXMLUnstructuredGridWriter
,vtu
);
260 vtu
->SetFileName("connections.vtu");
261 vtu
->SetDataModeToBinary();
262 vtu
->SetInput(conn_grid
);
266 EG_VTKSP(vtkXMLUnstructuredGridWriter
,vtu
);
268 vtu
->SetFileName("prevol.vtu");
269 vtu
->SetDataModeToBinary();
277 void CreateVolumeMesh::computeMeshDensity()
279 using namespace nglib
;
280 QVector
<vtkIdType
> cells
;
281 QVector
<vtkIdType
> nodes
;
283 QVector
<QVector
<int> > c2c
;
284 QVector
<QSet
<int> > n2n
;
285 getAllCellsOfType(VTK_TETRA
, cells
, grid
);
286 getNodesFromCells(cells
, nodes
, grid
);
287 createNodeMapping(nodes
, _nodes
, grid
);
288 createCellToCell(cells
, c2c
, grid
);
289 createNodeToNode(cells
, nodes
, _nodes
, n2n
, grid
);
290 QVector
<bool> fixed(nodes
.size(), false);
291 for (int i_cells
= 0; i_cells
< cells
.size(); ++i_cells
) {
292 vtkIdType id_cell
= cells
[i_cells
];
293 vtkIdType N_pts
, *pts
;
294 grid
->GetCellPoints(id_cell
, N_pts
, pts
);
295 if (c2c
[i_cells
][0] == -1) {
296 fixed
[_nodes
[pts
[0]]] = true;
297 fixed
[_nodes
[pts
[1]]] = true;
298 fixed
[_nodes
[pts
[2]]] = true;
300 if (c2c
[i_cells
][1] == -1) {
301 fixed
[_nodes
[pts
[0]]] = true;
302 fixed
[_nodes
[pts
[1]]] = true;
303 fixed
[_nodes
[pts
[3]]] = true;
305 if (c2c
[i_cells
][2] == -1) {
306 fixed
[_nodes
[pts
[0]]] = true;
307 fixed
[_nodes
[pts
[2]]] = true;
308 fixed
[_nodes
[pts
[3]]] = true;
310 if (c2c
[i_cells
][3] == -1) {
311 fixed
[_nodes
[pts
[1]]] = true;
312 fixed
[_nodes
[pts
[2]]] = true;
313 fixed
[_nodes
[pts
[3]]] = true;
316 QVector
<double> H(nodes
.size(), 0.0);
319 bool first_node
= true;
321 for (int i_nodes
= 0; i_nodes
< nodes
.size(); ++i_nodes
) {
322 if (fixed
[i_nodes
]) {
326 grid
->GetPoint(nodes
[i_nodes
], xi
.data());
332 for (int k
= 0; k
< 3; ++k
) {
333 X1
[k
] = min(xi
[k
], X1
[k
]);
334 X2
[k
] = max(xi
[k
], X2
[k
]);
337 foreach (int j_nodes
, n2n
[i_nodes
]) {
338 if (fixed
[j_nodes
]) {
340 grid
->GetPoint(nodes
[j_nodes
], xj
.data());
341 H
[i_nodes
] += (xi
-xj
).abs();
349 H_min
= min(H
[i_nodes
], H_min
);
354 cout
<< "relaxing mesh size : " << num
.toAscii().data() << "% done" << endl
;
355 if (N_non_fixed
> 0) {
359 for (int i_nodes
= 0; i_nodes
< nodes
.size(); ++i_nodes
) {
360 if (!fixed
[i_nodes
]) {
361 double H0
= H
[i_nodes
];
364 foreach (int j_nodes
, n2n
[i_nodes
]) {
365 H
[i_nodes
] += H
[j_nodes
];
372 DH_max
= max(H
[i_nodes
] - H0
, DH_max
);
376 double e
= min(1.0,max(0.0,-log10(DH_max
/H_min
)/3));
377 new_num
.setNum(100*e
,'f',0);
378 if (new_num
!= num
) {
380 cout
<< "relaxing mesh size : " << num
.toAscii().data() << "% done" << endl
;
382 } while (DH_max
> 1e-3*H_min
);
383 for (int i_nodes
= 0; i_nodes
< nodes
.size(); ++i_nodes
) {
385 grid
->GetPoint(nodes
[i_nodes
], x1
.data());
387 foreach (int j_nodes
, n2n
[i_nodes
]) {
389 grid
->GetPoint(nodes
[j_nodes
], xj
.data());
390 for (int k
= 0; k
< 3; ++k
) {
391 x1
[k
] = min(xj
[k
], x1
[k
]);
392 x2
[k
] = max(xj
[k
], x2
[k
]);
404 void CreateVolumeMesh::operate()
406 using namespace nglib
;
408 Ng_Meshing_Parameters mp
;
410 mp
.fineness
= fineness
;
412 Ng_Mesh
*mesh
= Ng_NewMesh();
413 computeMeshDensity();
415 QVector
<vtkIdType
> ng2eg(num_nodes_to_add
+1);
416 QVector
<vtkIdType
> eg2ng(grid
->GetNumberOfPoints());
419 for (vtkIdType id_node
= 0; id_node
< grid
->GetNumberOfPoints(); ++id_node
) {
420 if (add_to_ng
[id_node
]) {
422 grid
->GetPoints()->GetPoint(id_node
, x
.data());
423 Ng_AddPoint(mesh
, x
.data());
431 foreach (QVector
<vtkIdType
> T
, tri
) {
433 for (int i
= 0; i
< 3; ++i
) {
434 trig
[i
] = eg2ng
[T
[i
]];
436 Ng_AddSurfaceElement(mesh
, NG_TRIG
, trig
);
440 foreach (box_t B
, boxes
) Ng_RestrictMeshSizeBox(mesh
, B
.x1
.data(), B
.x2
.data(), B
.h
);
441 res
= Ng_GenerateVolumeMesh (mesh
, &mp
);
442 } catch (netgen::NgException ng_err
) {
444 QString msg
= "Netgen stopped with the following error:\n";
445 msg
+= ng_err
.What().c_str();
446 err
.setType(Error::ExitOperation
);
451 int Npoints_ng
= Ng_GetNP(mesh
);
452 int Ncells_ng
= Ng_GetNE(mesh
);
453 int Nscells_ng
= Ng_GetNSE(mesh
);
454 EG_VTKSP(vtkUnstructuredGrid
,vol_grid
);
455 allocateGrid(vol_grid
, grid
->GetNumberOfCells() + Ncells_ng
, Npoints_ng
+ num_old_nodes
);
456 vtkIdType new_point
= 0;
458 // copy existing points
459 QVector
<vtkIdType
> old2new(grid
->GetNumberOfPoints(), -1);
460 for (vtkIdType id_point
= 0; id_point
< grid
->GetNumberOfPoints(); ++id_point
) {
462 grid
->GetPoints()->GetPoint(id_point
, x
.data());
463 vol_grid
->GetPoints()->SetPoint(new_point
, x
.data());
464 copyNodeData(grid
, id_point
, vol_grid
, new_point
);
465 old2new
[id_point
] = new_point
;
469 // mark all surface nodes coming from NETGEN
470 QVector
<bool> ng_surf_node(Npoints_ng
+ 1, false);
471 for (int i
= 1; i
<= Nscells_ng
; ++i
) {
473 Ng_Surface_Element_Type ng_type
;
474 ng_type
= Ng_GetSurfaceElement(mesh
, i
, pts
);
476 if (ng_type
== NG_TRIG
) {
478 } else if (ng_type
== NG_QUAD
) {
483 for (int j
= 0; j
< N
; ++j
) {
484 ng_surf_node
[pts
[j
]] = true;
488 // add new points from NETGEN
489 QVector
<vtkIdType
> ng2new(Npoints_ng
+1, -1);
490 for (int i
= 1; i
<= Npoints_ng
; ++i
) {
491 if (!ng_surf_node
[i
]) {
493 Ng_GetPoint(mesh
, i
, x
.data());
494 vol_grid
->GetPoints()->SetPoint(new_point
, x
.data());
495 ng2new
[i
] = new_point
;
500 // copy existing cells
501 QVector
<vtkIdType
> old2new_cell(grid
->GetNumberOfCells(), -1);
502 for (vtkIdType id_cell
= 0; id_cell
< grid
->GetNumberOfCells(); ++id_cell
) {
503 bool ins_cell
= false;
504 vtkIdType type_cell
= grid
->GetCellType(id_cell
);
505 if (type_cell
== VTK_TRIANGLE
) ins_cell
= true;
506 if (type_cell
== VTK_QUAD
) ins_cell
= true;
507 if (type_cell
== VTK_WEDGE
) ins_cell
= true;
509 vtkIdType N_pts
, *pts
;
510 grid
->GetCellPoints(id_cell
, N_pts
, pts
);
511 for (int i
= 0; i
< N_pts
; ++i
) {
512 pts
[i
] = old2new
[pts
[i
]];
517 vtkIdType id_new
= vol_grid
->InsertNextCell(type_cell
, N_pts
, pts
);
518 copyCellData(grid
, id_cell
, vol_grid
, id_new
);
519 old2new_cell
[id_cell
] = id_new
;
524 vtkIdType id_new_cell
;
525 for (vtkIdType cellId
= 0; cellId
< Ncells_ng
; ++cellId
) {
527 vtkIdType new_pts
[4];
528 for (int i
= 0; i
< 8; ++i
) {
531 Ng_Volume_Element_Type ng_type
;
532 ng_type
= Ng_GetVolumeElement(mesh
, cellId
+ 1, pts
);
533 if (ng_type
!= NG_TET
) {
536 for (int i
= 0; i
< 4; ++i
) {
537 if (!ng_surf_node
[pts
[i
]]) {
538 new_pts
[i
] = ng2new
[pts
[i
]];
540 new_pts
[i
] = ng2eg
[pts
[i
]];
543 if (ng_type
== NG_TET
) {
549 id_new_cell
= vol_grid
->InsertNextCell(VTK_TETRA
, 4, tet
);
550 } else if (ng_type
== NG_PYRAMID
) {
551 EG_ERR_RETURN("pyramids cannot be handled yet");
552 } else if (ng_type
== NG_PRISM
) {
553 EG_ERR_RETURN("prisms cannot be handled yet");
555 EG_ERR_RETURN("bug encountered");
558 makeCopy(vol_grid
, grid
);
559 for (int i
= 0; i
< trace_cells
.size(); ++i
) {
560 if (old2new_cell
[trace_cells
[i
]] == -1) {
563 trace_cells
[i
] = old2new_cell
[trace_cells
[i
]];
567 QString msg
= "NETGEN did not succeed.\nPlease check if the surface mesh is oriented correctly";
568 msg
+= " (red edges on the outside of the domain)";
569 err
.setType(Error::ExitOperation
);
575 cout
<< "\n\nNETGEN call finished" << endl
;