createvolumemesh.cpp

   1 //
   2 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   3 // +                                                                      +
   4 // + This file is part of enGrid.                                         +
   5 // +                                                                      +
   6 // + Copyright 2008,2009 Oliver Gloth                                     +
   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 #include "createvolumemesh.h"
  24 #include "deletetetras.h"
  25 #include <vtkXMLUnstructuredGridWriter.h>
  26
  27 CreateVolumeMesh::CreateVolumeMesh()
  28 {
  29   maxh     = 1e99;
  30   fineness = 0.0;
  31 };
  32
  33 void CreateVolumeMesh::setTraceCells(const QVector<vtkIdType> &cells)
  34 {
  35   trace_cells.resize(cells.size());
  36   qCopy(cells.begin(), cells.end(), trace_cells.begin());
  37 };
  38
  39 void CreateVolumeMesh::getTraceCells(QVector<vtkIdType> &cells)
  40 {
  41   cells.resize(trace_cells.size());
  42   qCopy(trace_cells.begin(), trace_cells.end(), cells.begin());
  43 };
  44
  45 void CreateVolumeMesh::prepare()
  46 {
  47   using namespace nglib;
  48   DeleteTetras del;
  49   del.setGrid(grid);
  50   del();
  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;
  63   QList<vec3_t> conn1;
  64   QList<vec3_t> conn2;
  65   int N1 = 0;
  66   int N2 = 0;
  67   int N3 = 0;
  68   int N4 = 0;
  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) {
  76         T[0] = pts[0];
  77         T[1] = pts[1];
  78         T[2] = pts[2];
  79         ex_tri.append(T);
  80         ++N4;
  81       };
  82     } else if (type_cell == VTK_QUAD) {
  83       if (findVolumeCell(grid, id_cell, _nodes, cells, _cells, n2c) == -1) {
  84         EG_BUG;
  85         T[0] = pts[0];
  86         T[1] = pts[1];
  87         T[2] = pts[2];
  88         ex_tri.append(T);
  89         T[0] = pts[2];
  90         T[1] = pts[3];
  91         T[2] = pts[0];
  92         ex_tri.append(T);
  93         ++N3;
  94       };
  95     } else if (type_cell == VTK_WEDGE) {
  96       if (c2c[id_cell][0] == -1) {
  97         T[0] = pts[0];
  98         T[1] = pts[2];
  99         T[2] = pts[1];
 100         ex_tri.append(T);
 101         ++N2;
 102       };
 103       if (c2c[id_cell][1] == -1) {
 104         T[0] = pts[3];
 105         T[1] = pts[4];
 106         T[2] = pts[5];
 107         ex_tri.append(T);
 108         ++N2;
 109       };
 110       if (c2c[id_cell][2] == -1) {
 111         T[0] = pts[0];
 112         T[1] = pts[1];
 113         T[2] = pts[4];
 114         ex_tri.append(T);
 115         T[0] = pts[0];
 116         T[1] = pts[4];
 117         T[2] = pts[3];
 118         ex_tri.append(T);
 119         ++N1;
 120         vec3_t x[6];
 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]);
 123         centres.append(xc);
 124       } else {
 125         conn1.append(cellCentre(grid, id_cell));
 126         conn2.append(cellCentre(grid, cells[c2c[id_cell][2]]));
 127       };
 128       if (c2c[id_cell][3] == -1) {
 129         T[0] = pts[4];
 130         T[1] = pts[1];
 131         T[2] = pts[2];
 132         ex_tri.append(T);
 133         T[0] = pts[4];
 134         T[1] = pts[2];
 135         T[2] = pts[5];
 136         ex_tri.append(T);
 137         ++N1;
 138         vec3_t x[6];
 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]);
 141         centres.append(xc);
 142       } else {
 143         conn1.append(cellCentre(grid, id_cell));
 144         conn2.append(cellCentre(grid, cells[c2c[id_cell][3]]));
 145       };
 146       if (c2c[id_cell][4] == -1) {
 147         T[0] = pts[0];
 148         T[1] = pts[3];
 149         T[2] = pts[2];
 150         ex_tri.append(T);
 151         T[0] = pts[3];
 152         T[1] = pts[5];
 153         T[2] = pts[2];
 154         ex_tri.append(T);
 155         ++N1;
 156         vec3_t x[6];
 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]);
 159         centres.append(xc);
 160       } else {
 161         conn1.append(cellCentre(grid, id_cell));
 162         conn2.append(cellCentre(grid, cells[c2c[id_cell][4]]));
 163       };
 164     } else {
 165       EG_BUG;
 166     };
 167   };
 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;
 181   };
 182   num_nodes_to_add = 0;
 183   num_old_nodes = 0;
 184   for (vtkIdType id_node = 0; id_node < grid->GetNumberOfPoints(); ++id_node) {
 185     if (add_to_ng[id_node]) {
 186       ++num_nodes_to_add;
 187     } else {
 188       ++num_old_nodes;
 189     };
 190   };
 191
 192   {
 193     QVector<vtkIdType> old2tri(grid->GetNumberOfPoints(),-1);
 194     int N = 0;
 195     for (vtkIdType id_node = 0; id_node < grid->GetNumberOfPoints(); ++id_node) {
 196       if (add_to_ng[id_node]) {
 197         old2tri[id_node] = N;
 198         ++N;
 199       };
 200     };
 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]) {
 205         vec3_t x;
 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]);
 209       };
 210     };
 211     foreach (QVector<vtkIdType> T, tri) {
 212       vtkIdType pts[3];
 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);
 217     };
 218     EG_VTKSP(vtkXMLUnstructuredGridWriter,vtu);
 219     vtu->SetFileName("triangles.vtu");
 220     vtu->SetDataModeToBinary();
 221     vtu->SetInput(tri_grid);
 222     vtu->Write();
 223   };
 224   {
 225     EG_VTKSP(vtkUnstructuredGrid,centre_grid);
 226     allocateGrid(centre_grid, centres.size(), centres.size());
 227     vtkIdType N = 0;
 228     vtkIdType pts[1];
 229     foreach (vec3_t x, centres) {
 230       centre_grid->GetPoints()->SetPoint(N, x.data());
 231       pts[0] = N;
 232       centre_grid->InsertNextCell(VTK_VERTEX, 1, pts);
 233       ++N;
 234     };
 235     EG_VTKSP(vtkXMLUnstructuredGridWriter,vtu);
 236     vtu->SetFileName("centres.vtu");
 237     vtu->SetDataModeToBinary();
 238     vtu->SetInput(centre_grid);
 239     vtu->Write();
 240   };
 241   {
 242     EG_VTKSP(vtkUnstructuredGrid,conn_grid);
 243     allocateGrid(conn_grid, conn1.size(), 2*conn1.size());
 244     vtkIdType N = 0;
 245     vtkIdType pts[2];
 246     foreach (vec3_t x, conn1) {
 247       conn_grid->GetPoints()->SetPoint(N, x.data());
 248       ++N;
 249     };
 250     foreach (vec3_t x, conn2) {
 251       conn_grid->GetPoints()->SetPoint(N, x.data());
 252       ++N;
 253     };
 254     for (int i = 0; i < conn1.size(); ++i) {
 255       pts[0] = i;
 256       pts[1] = i+conn1.size();
 257       conn_grid->InsertNextCell(VTK_LINE, 2, pts);
 258     };
 259     EG_VTKSP(vtkXMLUnstructuredGridWriter,vtu);
 260     vtu->SetFileName("connections.vtu");
 261     vtu->SetDataModeToBinary();
 262     vtu->SetInput(conn_grid);
 263     vtu->Write();
 264   };
 265   {
 266     EG_VTKSP(vtkXMLUnstructuredGridWriter,vtu);
 267     createIndices(grid);
 268     vtu->SetFileName("prevol.vtu");
 269     vtu->SetDataModeToBinary();
 270     vtu->SetInput(grid);
 271     vtu->Write();
 272   };
 273
 274
 275 };
 276
 277 void CreateVolumeMesh::computeMeshDensity()
 278 {
 279   using namespace nglib;
 280   QVector<vtkIdType>  cells;
 281   QVector<vtkIdType>  nodes;
 282   QVector<int>       _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;
 299     };
 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;
 304     };
 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;
 309     };
 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;
 314     };
 315   };
 316   QVector<double> H(nodes.size(), 0.0);
 317   double H_min = 1e99;
 318   vec3_t X1, X2;
 319   bool first_node = true;
 320   int N_non_fixed = 0;
 321   for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
 322     if (fixed[i_nodes]) {
 323       ++N_non_fixed;
 324       int N = 0;
 325       vec3_t xi;
 326       grid->GetPoint(nodes[i_nodes], xi.data());
 327       if (first_node) {
 328         X1 = xi;
 329         X2 = xi;
 330         first_node = false;
 331       } else {
 332         for (int k = 0; k < 3; ++k) {
 333           X1[k] = min(xi[k], X1[k]);
 334           X2[k] = max(xi[k], X2[k]);
 335         };
 336       };
 337       foreach (int j_nodes, n2n[i_nodes]) {
 338         if (fixed[j_nodes]) {
 339           vec3_t xj;
 340           grid->GetPoint(nodes[j_nodes], xj.data());
 341           H[i_nodes] += (xi-xj).abs();
 342           ++N;
 343         };
 344       };
 345       if (N < 2) {
 346         EG_BUG;
 347       };
 348       H[i_nodes] /= N;
 349       H_min = min(H[i_nodes], H_min);
 350     };
 351   };
 352   boxes.clear();
 353   QString num = "0";
 354   cout << "relaxing mesh size : " << num.toAscii().data() << "% done" << endl;
 355   if (N_non_fixed > 0) {
 356     double DH_max;
 357     do {
 358       DH_max = 0;
 359       for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
 360         if (!fixed[i_nodes]) {
 361           double H0 = H[i_nodes];
 362           H[i_nodes] = 0.0;
 363           int N = 0;
 364           foreach (int j_nodes, n2n[i_nodes]) {
 365             H[i_nodes] += H[j_nodes];
 366             ++N;
 367           };
 368           if (N == 0) {
 369             EG_BUG;
 370           };
 371           H[i_nodes] /= N;
 372           DH_max = max(H[i_nodes] - H0, DH_max);
 373         };
 374       };
 375       QString new_num;
 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) {
 379         num = new_num;
 380         cout << "relaxing mesh size : " << num.toAscii().data() << "% done" << endl;
 381       };
 382     } while (DH_max > 1e-3*H_min);
 383     for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
 384       vec3_t x1, x2;
 385       grid->GetPoint(nodes[i_nodes], x1.data());
 386       x2 = x1;
 387       foreach (int j_nodes, n2n[i_nodes]) {
 388         vec3_t xj;
 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]);
 393         };
 394       };
 395       box_t B;
 396       B.x1 =x1;
 397       B.x2 =x2;
 398       B.h = H[i_nodes];
 399       boxes.append(B);
 400     };
 401   };
 402 };
 403
 404 void CreateVolumeMesh::operate()
 405 {
 406   using namespace nglib;
 407   Ng_Init();
 408   Ng_Meshing_Parameters mp;
 409   mp.maxh = maxh;
 410   mp.fineness = fineness;
 411   mp.secondorder = 0;
 412   Ng_Mesh *mesh = Ng_NewMesh();
 413   computeMeshDensity();
 414   prepare();
 415   QVector<vtkIdType> ng2eg(num_nodes_to_add+1);
 416   QVector<vtkIdType> eg2ng(grid->GetNumberOfPoints());
 417   {
 418     int N = 1;
 419     for (vtkIdType id_node = 0; id_node < grid->GetNumberOfPoints(); ++id_node) {
 420       if (add_to_ng[id_node]) {
 421         vec3_t x;
 422         grid->GetPoints()->GetPoint(id_node, x.data());
 423         Ng_AddPoint(mesh, x.data());
 424         ng2eg[N] = id_node;
 425         eg2ng[id_node] = N;
 426         ++N;
 427       };
 428     };
 429   };
 430
 431   foreach (QVector<vtkIdType> T, tri) {
 432     int trig[3];
 433     for (int i = 0; i < 3; ++i) {
 434       trig[i] = eg2ng[T[i]];
 435     };
 436     Ng_AddSurfaceElement(mesh, NG_TRIG, trig);
 437   };
 438   Ng_Result res;
 439   try {
 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) {
 443     Error err;
 444     QString msg = "Netgen stopped with the following error:\n";
 445     msg += ng_err.What().c_str();
 446     err.setType(Error::ExitOperation);
 447     err.setText(msg);
 448     throw err;
 449   };
 450   if (res == NG_OK) {
 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;
 457
 458     // copy existing points
 459     QVector<vtkIdType> old2new(grid->GetNumberOfPoints(), -1);
 460     for (vtkIdType id_point = 0; id_point < grid->GetNumberOfPoints(); ++id_point) {
 461       vec3_t x;
 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;
 466       ++new_point;
 467     };
 468
 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) {
 472       int pts[8];
 473       Ng_Surface_Element_Type ng_type;
 474       ng_type = Ng_GetSurfaceElement(mesh, i, pts);
 475       int N = 0;
 476       if (ng_type == NG_TRIG) {
 477         N = 3;
 478       } else if (ng_type == NG_QUAD) {
 479         N = 4;
 480       } else {
 481         EG_BUG;
 482       };
 483       for (int j = 0; j < N; ++j) {
 484         ng_surf_node[pts[j]] = true;
 485       };
 486     };
 487
 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]) {
 492         vec3_t x;
 493         Ng_GetPoint(mesh, i, x.data());
 494         vol_grid->GetPoints()->SetPoint(new_point, x.data());
 495         ng2new[i] = new_point;
 496         ++new_point;
 497       };
 498     };
 499
 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;
 508       if (ins_cell) {
 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]];
 513           if (pts[i] == -1) {
 514             EG_BUG;
 515           };
 516         };
 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;
 520       };
 521     };
 522
 523     // add new cells
 524     vtkIdType id_new_cell;
 525     for (vtkIdType cellId = 0; cellId < Ncells_ng; ++cellId) {
 526       int       pts[8];
 527       vtkIdType new_pts[4];
 528       for (int i = 0; i < 8; ++i) {
 529         pts[i] = 0;
 530       };
 531       Ng_Volume_Element_Type ng_type;
 532       ng_type = Ng_GetVolumeElement(mesh, cellId + 1, pts);
 533       if (ng_type != NG_TET) {
 534         EG_BUG;
 535       };
 536       for (int i = 0; i < 4; ++i) {
 537         if (!ng_surf_node[pts[i]]) {
 538           new_pts[i] = ng2new[pts[i]];
 539         } else {
 540           new_pts[i] = ng2eg[pts[i]];
 541         };
 542       };
 543       if (ng_type == NG_TET) {
 544         vtkIdType tet[4];
 545         tet[0] = new_pts[0];
 546         tet[1] = new_pts[1];
 547         tet[2] = new_pts[3];
 548         tet[3] = new_pts[2];
 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");
 554       } else {
 555         EG_ERR_RETURN("bug encountered");
 556       };
 557     };
 558     makeCopy(vol_grid, grid);
 559     for (int i = 0; i < trace_cells.size(); ++i) {
 560       if (old2new_cell[trace_cells[i]] == -1) {
 561         EG_BUG;
 562       };
 563       trace_cells[i] = old2new_cell[trace_cells[i]];
 564     };
 565   } else {
 566     Error err;
 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);
 570     err.setText(msg);
 571     throw err;
 572   };
 573   Ng_DeleteMesh(mesh);
 574   Ng_Exit();
 575   cout << "\n\nNETGEN call finished" << endl;
 576   cout << endl;
 577 };
 578