src/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 "guimainwindow.h"
  26 #include <vtkXMLUnstructuredGridWriter.h>
  27
  28 CreateVolumeMesh::CreateVolumeMesh()
  29 {
  30   EG_TYPENAME;
  31   maxh     = 1e99;
  32   fineness = 0.0;
  33 }
  34
  35 void CreateVolumeMesh::setTraceCells(const QVector<vtkIdType> &cells)
  36 {
  37   trace_cells.resize(cells.size());
  38   qCopy(cells.begin(), cells.end(), trace_cells.begin());
  39 }
  40
  41 void CreateVolumeMesh::getTraceCells(QVector<vtkIdType> &cells)
  42 {
  43   cells.resize(trace_cells.size());
  44   qCopy(trace_cells.begin(), trace_cells.end(), cells.begin());
  45 }
  46
  47 void CreateVolumeMesh::prepare()
  48 {
  49   using namespace nglib;
  50   DeleteTetras del;
  51   del.setGrid(grid);
  52   del.setAllCells();
  53   del();
  54   QVector<vtkIdType> cells, nodes;
  55   QVector<int>       _cells, _nodes;
  56   QVector<QVector< int > > c2c;
  57   QVector<QVector<int> > n2c;
  58   getAllCells(cells, grid);
  59   createCellMapping(cells, _cells, grid);
  60   getNodesFromCells(cells, nodes, grid);
  61   createNodeMapping(nodes, _nodes, grid);
  62   createCellToCell(cells, c2c, grid);
  63   createNodeToCell(cells, nodes, _nodes, n2c, grid);
  64   QList<QVector<vtkIdType> > ex_tri;
  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       }
 121       if (c2c[id_cell][3] == -1) {
 122         T[0] = pts[4];
 123         T[1] = pts[1];
 124         T[2] = pts[2];
 125         ex_tri.append(T);
 126         T[0] = pts[4];
 127         T[1] = pts[2];
 128         T[2] = pts[5];
 129         ex_tri.append(T);
 130         ++N1;
 131       }
 132       if (c2c[id_cell][4] == -1) {
 133         T[0] = pts[0];
 134         T[1] = pts[3];
 135         T[2] = pts[2];
 136         ex_tri.append(T);
 137         T[0] = pts[3];
 138         T[1] = pts[5];
 139         T[2] = pts[2];
 140         ex_tri.append(T);
 141         ++N1;
 142       }
 143     } else {
 144       EG_BUG;
 145     }
 146   }
 147   cout << "*********************************************************************" << endl;
 148   cout << "prism quads     : " << N1 << endl;
 149   cout << "prism triangles : " << N2 << endl;
 150   cout << "stray quads     : " << N3 << endl;
 151   cout << "stray triangles : " << N4 << endl;
 152   cout << "*********************************************************************" << endl;
 153   tri.resize(ex_tri.size());
 154   qCopy(ex_tri.begin(), ex_tri.end(), tri.begin());
 155   add_to_ng.fill(false, grid->GetNumberOfPoints());
 156   foreach (QVector<vtkIdType> T, tri) {
 157     add_to_ng[T[0]] = true;
 158     add_to_ng[T[1]] = true;
 159     add_to_ng[T[2]] = true;
 160   }
 161   num_nodes_to_add = 0;
 162   num_old_nodes = 0;
 163   for (vtkIdType id_node = 0; id_node < grid->GetNumberOfPoints(); ++id_node) {
 164     if (add_to_ng[id_node]) {
 165       ++num_nodes_to_add;
 166     } else {
 167       ++num_old_nodes;
 168     }
 169   }
 170   old2tri.fill(-1, grid->GetNumberOfPoints());
 171   m_NumTriangles = 0;
 172   for (vtkIdType id_node = 0; id_node < grid->GetNumberOfPoints(); ++id_node) {
 173     if (add_to_ng[id_node]) {
 174       old2tri[id_node] = m_NumTriangles;
 175       ++m_NumTriangles;
 176     }
 177   }
 178 }
 179
 180 void CreateVolumeMesh::writeDebugInfo()
 181 {
 182   {
 183     EG_VTKSP(vtkUnstructuredGrid,tri_grid);
 184     allocateGrid(tri_grid, tri.size(), m_NumTriangles);
 185     for (vtkIdType id_node = 0; id_node < grid->GetNumberOfPoints(); ++id_node) {
 186       if (add_to_ng[id_node]) {
 187         vec3_t x;
 188         grid->GetPoint(id_node, x.data());
 189         tri_grid->GetPoints()->SetPoint(old2tri[id_node], x.data());
 190         copyNodeData(grid, id_node, tri_grid, old2tri[id_node]);
 191       }
 192     }
 193     foreach (QVector<vtkIdType> T, tri) {
 194       vtkIdType pts[3];
 195       pts[0] = old2tri[T[0]];
 196       pts[1] = old2tri[T[1]];
 197       pts[2] = old2tri[T[2]];
 198       tri_grid->InsertNextCell(VTK_TRIANGLE, 3, pts);
 199     }
 200     writeGrid(tri_grid, "triangles");
 201   }
 202   {
 203     writeGrid(grid, "last_grid");
 204   }
 205 }
 206
 207 void CreateVolumeMesh::computeMeshDensity()
 208 {
 209   using namespace nglib;
 210   QVector<vtkIdType>  cells;
 211   QVector<vtkIdType>  nodes;
 212   QVector<int>       _nodes;
 213   QVector<QVector<int> >  c2c;
 214   QVector<QSet<int> >     n2n;
 215   getAllCellsOfType(VTK_TETRA, cells, grid);
 216   getNodesFromCells(cells, nodes, grid);
 217   createNodeMapping(nodes, _nodes, grid);
 218   createCellToCell(cells, c2c, grid);
 219   createNodeToNode(cells, nodes, _nodes, n2n, grid);
 220   QVector<bool> fixed(nodes.size(), false);
 221   for (int i_cells = 0; i_cells < cells.size(); ++i_cells) {
 222     vtkIdType id_cell = cells[i_cells];
 223     vtkIdType N_pts, *pts;
 224     grid->GetCellPoints(id_cell, N_pts, pts);
 225     if (c2c[i_cells][0] == -1) {
 226       fixed[_nodes[pts[0]]] = true;
 227       fixed[_nodes[pts[1]]] = true;
 228       fixed[_nodes[pts[2]]] = true;
 229     }
 230     if (c2c[i_cells][1] == -1) {
 231       fixed[_nodes[pts[0]]] = true;
 232       fixed[_nodes[pts[1]]] = true;
 233       fixed[_nodes[pts[3]]] = true;
 234     }
 235     if (c2c[i_cells][2] == -1) {
 236       fixed[_nodes[pts[0]]] = true;
 237       fixed[_nodes[pts[2]]] = true;
 238       fixed[_nodes[pts[3]]] = true;
 239     }
 240     if (c2c[i_cells][3] == -1) {
 241       fixed[_nodes[pts[1]]] = true;
 242       fixed[_nodes[pts[2]]] = true;
 243       fixed[_nodes[pts[3]]] = true;
 244     }
 245   }
 246   QVector<double> H(nodes.size(), 0.0);
 247   double H_min = 1e99;
 248   vec3_t X1, X2;
 249   bool first_node = true;
 250   int N_non_fixed = 0;
 251   for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
 252     if (fixed[i_nodes]) {
 253       ++N_non_fixed;
 254       int N = 0;
 255       vec3_t xi;
 256       grid->GetPoint(nodes[i_nodes], xi.data());
 257       if (first_node) {
 258         X1 = xi;
 259         X2 = xi;
 260         first_node = false;
 261       } else {
 262         for (int k = 0; k < 3; ++k) {
 263           X1[k] = min(xi[k], X1[k]);
 264           X2[k] = max(xi[k], X2[k]);
 265         }
 266       }
 267       foreach (int j_nodes, n2n[i_nodes]) {
 268         if (fixed[j_nodes]) {
 269           vec3_t xj;
 270           grid->GetPoint(nodes[j_nodes], xj.data());
 271           H[i_nodes] += (xi-xj).abs();
 272           ++N;
 273         }
 274       }
 275       if (N < 2) {
 276         EG_BUG;
 277       }
 278       H[i_nodes] /= N;
 279       H_min = min(H[i_nodes], H_min);
 280     }
 281   }
 282   boxes.clear();
 283   QString num = "0";
 284   cout << "relaxing mesh size : " << qPrintable(num) << "% done" << endl;
 285   if (N_non_fixed > 0) {
 286     double DH_max;
 287     do {
 288       DH_max = 0;
 289       for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
 290         if (!fixed[i_nodes]) {
 291           double H0 = H[i_nodes];
 292           H[i_nodes] = 0.0;
 293           int N = 0;
 294           foreach (int j_nodes, n2n[i_nodes]) {
 295             H[i_nodes] += H[j_nodes];
 296             ++N;
 297           }
 298           if (N == 0) {
 299             EG_BUG;
 300           }
 301           H[i_nodes] /= N;
 302           DH_max = max(H[i_nodes] - H0, DH_max);
 303         }
 304       }
 305       QString new_num;
 306       double e = min(1.0,max(0.0,-log10(DH_max/H_min)/3));
 307       new_num.setNum(100*e,'f',0);
 308       if (new_num != num) {
 309         num = new_num;
 310         cout << "relaxing mesh size : " << qPrintable(num) << "% done" << endl;
 311       }
 312     } while (DH_max > 1e-3*H_min);
 313     for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
 314       vec3_t x1, x2;
 315       grid->GetPoint(nodes[i_nodes], x1.data());
 316       x2 = x1;
 317       foreach (int j_nodes, n2n[i_nodes]) {
 318         vec3_t xj;
 319         grid->GetPoint(nodes[j_nodes], xj.data());
 320         for (int k = 0; k < 3; ++k) {
 321           x1[k] = min(xj[k], x1[k]);
 322           x2[k] = max(xj[k], x2[k]);
 323         }
 324       }
 325       box_t B;
 326       B.x1 =x1;
 327       B.x2 =x2;
 328       B.h = H[i_nodes];
 329       boxes.append(B);
 330     }
 331   }
 332 }
 333
 334 void CreateVolumeMesh::operate()
 335 {
 336   using namespace nglib;
 337   if (grid->GetNumberOfCells() == 0) {
 338     EG_ERR_RETURN("The grid appears to be empty.");
 339   }
 340   nglib::Ng_Init();
 341   Ng_Meshing_Parameters mp;
 342   mp.maxh = maxh;
 343   mp.fineness = fineness;
 344   mp.secondorder = 0;
 345   Ng_Mesh *mesh = Ng_NewMesh();
 346   computeMeshDensity();
 347   prepare();
 348   QVector<vtkIdType> ng2eg(num_nodes_to_add+1);
 349   QVector<vtkIdType> eg2ng(grid->GetNumberOfPoints());
 350   {
 351     int N = 1;
 352     for (vtkIdType id_node = 0; id_node < grid->GetNumberOfPoints(); ++id_node) {
 353       if (add_to_ng[id_node]) {
 354         vec3_t x;
 355         grid->GetPoints()->GetPoint(id_node, x.data());
 356         Ng_AddPoint(mesh, x.data());
 357         ng2eg[N] = id_node;
 358         eg2ng[id_node] = N;
 359         ++N;
 360       }
 361     }
 362   }
 363
 364   foreach (QVector<vtkIdType> T, tri) {
 365     int trig[3];
 366     for (int i = 0; i < 3; ++i) {
 367       trig[i] = eg2ng[T[i]];
 368     }
 369     Ng_AddSurfaceElement(mesh, NG_TRIG, trig);
 370   }
 371   Ng_Result res;
 372   try {
 373     foreach (box_t B, boxes) Ng_RestrictMeshSizeBox(mesh, B.x1.data(), B.x2.data(), B.h);
 374     res = Ng_GenerateVolumeMesh (mesh, &mp);
 375   } catch (netgen::NgException ng_err) {
 376     writeDebugInfo();
 377     Error err;
 378     QString msg = "Netgen stopped with the following error:\n";
 379     msg += ng_err.What().c_str();
 380     msg += "\n\nDebug information has been saved to:\n" + mainWindow()->getCwd();
 381     err.setType(Error::ExitOperation);
 382     err.setText(msg);
 383     throw err;
 384   }
 385   if (res == NG_OK) {
 386     int Npoints_ng = Ng_GetNP(mesh);
 387     int Ncells_ng  = Ng_GetNE(mesh);
 388     int Nscells_ng = Ng_GetNSE(mesh);
 389     EG_VTKSP(vtkUnstructuredGrid,vol_grid);
 390     allocateGrid(vol_grid, grid->GetNumberOfCells() + Ncells_ng, Npoints_ng + num_old_nodes);
 391     vtkIdType new_point = 0;
 392
 393     // copy existing points
 394     QVector<vtkIdType> old2new(grid->GetNumberOfPoints(), -1);
 395     for (vtkIdType id_point = 0; id_point < grid->GetNumberOfPoints(); ++id_point) {
 396       vec3_t x;
 397       grid->GetPoints()->GetPoint(id_point, x.data());
 398       vol_grid->GetPoints()->SetPoint(new_point, x.data());
 399       copyNodeData(grid, id_point, vol_grid, new_point);
 400       old2new[id_point] = new_point;
 401       ++new_point;
 402     }
 403
 404     // mark all surface nodes coming from NETGEN
 405     QVector<bool> ng_surf_node(Npoints_ng + 1, false);
 406     for (int i = 1; i <= Nscells_ng; ++i) {
 407       int pts[8];
 408       Ng_Surface_Element_Type ng_type;
 409       ng_type = Ng_GetSurfaceElement(mesh, i, pts);
 410       int N = 0;
 411       if (ng_type == NG_TRIG) {
 412         N = 3;
 413       } else if (ng_type == NG_QUAD) {
 414         N = 4;
 415       } else {
 416         EG_BUG;
 417       }
 418       for (int j = 0; j < N; ++j) {
 419         ng_surf_node[pts[j]] = true;
 420       }
 421     }
 422
 423     // add new points from NETGEN
 424     QVector<vtkIdType> ng2new(Npoints_ng+1, -1);
 425     for (int i = 1; i <= Npoints_ng; ++i) {
 426       if (!ng_surf_node[i]) {
 427         vec3_t x;
 428         Ng_GetPoint(mesh, i, x.data());
 429         vol_grid->GetPoints()->SetPoint(new_point, x.data());
 430         ng2new[i] = new_point;
 431         ++new_point;
 432       }
 433     }
 434
 435     // copy existing cells
 436     QVector<vtkIdType> old2new_cell(grid->GetNumberOfCells(), -1);
 437     for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
 438       bool ins_cell = false;
 439       vtkIdType type_cell = grid->GetCellType(id_cell);
 440       if (type_cell == VTK_TRIANGLE) ins_cell = true;
 441       if (type_cell == VTK_QUAD)     ins_cell = true;
 442       if (type_cell == VTK_WEDGE)    ins_cell = true;
 443       if (ins_cell) {
 444         vtkIdType N_pts, *pts;
 445         grid->GetCellPoints(id_cell, N_pts, pts);
 446         for (int i = 0; i < N_pts; ++i) {
 447           pts[i] = old2new[pts[i]];
 448           if (pts[i] == -1) {
 449             EG_BUG;
 450           }
 451         }
 452         vtkIdType id_new = vol_grid->InsertNextCell(type_cell, N_pts, pts);
 453         copyCellData(grid, id_cell, vol_grid, id_new);
 454         old2new_cell[id_cell] = id_new;
 455       }
 456     }
 457
 458     // add new cells
 459     vtkIdType id_new_cell;
 460     for (vtkIdType cellId = 0; cellId < Ncells_ng; ++cellId) {
 461       int       pts[8];
 462       vtkIdType new_pts[4];
 463       for (int i = 0; i < 8; ++i) {
 464         pts[i] = 0;
 465       }
 466       Ng_Volume_Element_Type ng_type;
 467       ng_type = Ng_GetVolumeElement(mesh, cellId + 1, pts);
 468       if (ng_type != NG_TET) {
 469         EG_BUG;
 470       }
 471       for (int i = 0; i < 4; ++i) {
 472         if (!ng_surf_node[pts[i]]) {
 473           new_pts[i] = ng2new[pts[i]];
 474         } else {
 475           new_pts[i] = ng2eg[pts[i]];
 476         }
 477       }
 478       if (ng_type == NG_TET) {
 479         vtkIdType tet[4];
 480         tet[0] = new_pts[0];
 481         tet[1] = new_pts[1];
 482         tet[2] = new_pts[3];
 483         tet[3] = new_pts[2];
 484         id_new_cell = vol_grid->InsertNextCell(VTK_TETRA, 4, tet);
 485       } else if (ng_type == NG_PYRAMID) {
 486         EG_ERR_RETURN("pyramids cannot be handled yet");
 487       } else if (ng_type == NG_PRISM) {
 488         EG_ERR_RETURN("prisms cannot be handled yet");
 489       } else {
 490         EG_ERR_RETURN("bug encountered");
 491       }
 492     }
 493     makeCopy(vol_grid, grid);
 494     for (int i = 0; i < trace_cells.size(); ++i) {
 495       if (old2new_cell[trace_cells[i]] == -1) {
 496         EG_BUG;
 497       }
 498       trace_cells[i] = old2new_cell[trace_cells[i]];
 499     }
 500   } else {
 501     Error err;
 502     QString msg = "NETGEN did not succeed.\nPlease check if the surface mesh is oriented correctly";
 503     msg += " (red edges on the outside of the domain)";
 504     err.setType(Error::ExitOperation);
 505     err.setText(msg);
 506     throw err;
 507   }
 508   Ng_DeleteMesh(mesh);
 509   Ng_Exit();
 510   cout << "\n\nNETGEN call finished" << endl;
 511   cout << endl;
 512 }
 513