src/libengrid/createvolumemesh.cpp

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