src/libengrid/createvolumemesh.cpp

   1 //
   2 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   3 // +                                                                      +
   4 // + This file is part of enGrid.                                         +
   5 // +                                                                      +
   6 // + Copyright 2008-2012 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     writeGrid(tri_grid, "triangles");
 203   }
 204   {
 205     writeGrid(m_Grid, "last_grid");
 206   }
 207 }
 208
 209 void CreateVolumeMesh::computeMeshDensity()
 210 {
 211   boxes.clear();
 212   QString buffer = GuiMainWindow::pointer()->getXmlSection("engrid/surface/settings").replace("\n", " ");
 213   if (!buffer.isEmpty()) {
 214     QTextStream in(&buffer, QIODevice::ReadOnly);
 215     in >> m_MaxEdgeLength;
 216     in >> m_MinEdgeLength;
 217     in >> m_GrowthFactor;
 218   } else {
 219     m_MaxEdgeLength = 1000.0;
 220     m_MinEdgeLength = 0.0;
 221     m_GrowthFactor = 1.5;
 222   }
 223   m_ELSManager.read();
 224   QVector<double> H(m_Grid->GetNumberOfPoints(), m_MaxEdgeLength);
 225
 226   QVector<bool> fixed(m_Grid->GetNumberOfPoints(), false);
 227   double H_min = 1e99;
 228   vtkIdType id_min = -1;
 229   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 230     bool volume_only = true;
 231     for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
 232       if (isSurface(m_Part.n2cGG(id_node, i), m_Grid)) {
 233         volume_only = false;
 234       }
 235       if (!volume_only) {
 236         fixed[id_node] = true;
 237         H[id_node] = 0;
 238         int N = 0;
 239         vec3_t xi;
 240         m_Grid->GetPoint(id_node, xi.data());
 241         for (int j = 0; j < m_Part.n2nGSize(id_node); ++j) {
 242           if (m_Part.n2nGG(id_node, j)) {
 243             vec3_t xj;
 244             m_Grid->GetPoint(m_Part.n2nGG(id_node, j), xj.data());
 245             H[id_node] += (xi-xj).abs();
 246             ++N;
 247           }
 248         }
 249         if (N < 2) {
 250           EG_BUG;
 251         }
 252         H[id_node] /= N;
 253                 if (H[id_node] < 0) {
 254               EG_BUG;
 255                 }
 256         if (H[id_node] < H_min) {
 257           id_min = id_node;
 258           H_min = H[id_node];
 259         }
 260       }
 261     }
 262   }
 263   if (id_min < 0) {
 264     EG_BUG;
 265   }
 266
 267   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 268     vec3_t x;
 269     m_Grid->GetPoint(id_node, x.data());
 270     double cl_src = m_ELSManager.minEdgeLength(x);
 271     if (cl_src > 0) {
 272       if (cl_src < H[id_node]) {
 273         H[id_node] = cl_src;
 274       }
 275     }
 276   }
 277
 278   QVector<bool> marked(m_Grid->GetNumberOfPoints(), false);
 279   marked[id_min] = true;
 280   bool done = false;
 281   while (!done) {
 282     done = true;
 283     for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 284       if (marked[id_node] && H[id_node] <= H_min) {
 285         vec3_t x1;
 286         m_Grid->GetPoint(id_node, x1.data());
 287         for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
 288           vtkIdType id_neigh = m_Part.n2nGG(id_node,i);
 289           if (!marked[id_neigh]) {
 290             vec3_t x2;
 291             m_Grid->GetPoint(id_neigh, x2.data());
 292             double dist = (x1 - x2).abs();
 293             double h = H[id_node] + (m_GrowthFactor - 1)*dist;
 294             if (h < 0) {
 295               EG_BUG;
 296             }
 297             H[id_neigh] = min(H[id_neigh], h);
 298             marked[id_neigh] = true;
 299             //H[id_neigh] += 1.0*H[id_node];
 300             done = false;
 301           }
 302         }
 303       }
 304     }
 305     H_min *= m_GrowthFactor;
 306   }
 307
 308   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 309     vec3_t x;
 310     m_Grid->GetPoint(id_node, x.data());
 311     double cl_src = m_ELSManager.minEdgeLength(x);
 312     if (cl_src > 0) {
 313       if (cl_src < H[id_node]) {
 314         H[id_node] = cl_src;
 315       }
 316     }
 317   }
 318
 319   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 320     vec3_t x1, x2;
 321     m_Grid->GetPoint(id_node, x1.data());
 322     x2 = x1;
 323     for (int j = 0; j < m_Part.n2nGSize(id_node); ++j) {
 324       vec3_t xj;
 325       m_Grid->GetPoint(m_Part.n2nGG(id_node, j), xj.data());
 326       for (int k = 0; k < 3; ++k) {
 327         x1[k] = min(xj[k], x1[k]);
 328         x2[k] = max(xj[k], x2[k]);
 329       }
 330     }
 331     box_t B;
 332     B.x1 =x1;
 333     B.x2 =x2;
 334     if (H[id_node] < 0) {
 335       EG_BUG;
 336     }
 337     B.h = H[id_node];
 338     boxes.append(B);
 339   }
 340 }
 341
 342
 343 void CreateVolumeMesh::operate()
 344 {
 345   using namespace nglib;
 346   setAllCells();
 347   if (m_Grid->GetNumberOfCells() == 0) {
 348     EG_ERR_RETURN("The grid appears to be empty.");
 349   }
 350   //nglib::Ng_Init();
 351   Ng_Init();
 352   Ng_Meshing_Parameters mp;
 353   mp.fineness = fineness;
 354   //mp.secondorder = 0;
 355   Ng_Mesh *mesh = Ng_NewMesh();
 356   computeMeshDensity();
 357   mp.maxh = m_MaxEdgeLength;
 358   mp.minh = m_MinEdgeLength;
 359   mp.grading = 1;
 360   prepare();
 361   QVector<vtkIdType> ng2eg(num_nodes_to_add+1);
 362   QVector<vtkIdType> eg2ng(m_Grid->GetNumberOfPoints());
 363   {
 364     int N = 1;
 365     for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 366       if (add_to_ng[id_node]) {
 367         vec3_t x;
 368         m_Grid->GetPoints()->GetPoint(id_node, x.data());
 369         Ng_AddPoint(mesh, x.data());
 370         ng2eg[N] = id_node;
 371         eg2ng[id_node] = N;
 372         ++N;
 373       }
 374     }
 375   }
 376
 377   foreach (QVector<vtkIdType> T, tri) {
 378     int trig[3];
 379     for (int i = 0; i < 3; ++i) {
 380       trig[i] = eg2ng[T[i]];
 381     }
 382     Ng_AddSurfaceElement(mesh, NG_TRIG, trig);
 383   }
 384   Ng_Result res;
 385   try {
 386     int box_counter = 0;
 387     foreach (box_t B, boxes) {
 388       Ng_RestrictMeshSizeBox(mesh, B.x1.data(), B.x2.data(), B.h);
 389       ++box_counter;
 390     }
 391     GuiMainWindow::pointer()->setSystemOutput();
 392     writeDebugInfo();
 393     res = Ng_GenerateVolumeMesh (mesh, &mp);
 394     GuiMainWindow::pointer()->setLogFileOutput();
 395   } catch (netgen::NgException ng_err) {
 396     GuiMainWindow::pointer()->setLogFileOutput();
 397     writeDebugInfo();
 398     Error err;
 399     QString msg = "Netgen stopped with the following error:\n";
 400     msg += ng_err.What().c_str();
 401     msg += "\n\nDebug information has been saved to:\n" + mainWindow()->getCwd();
 402     err.setType(Error::ExitOperation);
 403     err.setText(msg);
 404     throw err;
 405   }
 406   if (res == NG_OK) {
 407     int Npoints_ng = Ng_GetNP(mesh);
 408     int Ncells_ng  = Ng_GetNE(mesh);
 409     int Nscells_ng = Ng_GetNSE(mesh);
 410     EG_VTKSP(vtkUnstructuredGrid,vol_grid);
 411     allocateGrid(vol_grid, m_Grid->GetNumberOfCells() + Ncells_ng, Npoints_ng + num_old_nodes);
 412     vtkIdType new_point = 0;
 413
 414     // copy existing points
 415     QVector<vtkIdType> old2new(m_Grid->GetNumberOfPoints(), -1);
 416     for (vtkIdType id_point = 0; id_point < m_Grid->GetNumberOfPoints(); ++id_point) {
 417       vec3_t x;
 418       m_Grid->GetPoints()->GetPoint(id_point, x.data());
 419       vol_grid->GetPoints()->SetPoint(new_point, x.data());
 420       copyNodeData(m_Grid, id_point, vol_grid, new_point);
 421       old2new[id_point] = new_point;
 422       ++new_point;
 423     }
 424
 425     // mark all surface nodes coming from NETGEN
 426     QVector<bool> ng_surf_node(Npoints_ng + 1, false);
 427     for (int i = 1; i <= Nscells_ng; ++i) {
 428       int pts[8];
 429       Ng_Surface_Element_Type ng_type;
 430       ng_type = Ng_GetSurfaceElement(mesh, i, pts);
 431       int N = 0;
 432       if (ng_type == NG_TRIG) {
 433         N = 3;
 434       } else if (ng_type == NG_QUAD) {
 435         N = 4;
 436       } else {
 437         EG_BUG;
 438       }
 439       for (int j = 0; j < N; ++j) {
 440         ng_surf_node[pts[j]] = true;
 441       }
 442     }
 443
 444     // add new points from NETGEN
 445     QVector<vtkIdType> ng2new(Npoints_ng+1, -1);
 446     for (int i = 1; i <= Npoints_ng; ++i) {
 447       if (!ng_surf_node[i]) {
 448         vec3_t x;
 449         Ng_GetPoint(mesh, i, x.data());
 450         vol_grid->GetPoints()->SetPoint(new_point, x.data());
 451         ng2new[i] = new_point;
 452         ++new_point;
 453       }
 454     }
 455
 456     // copy existing cells
 457     QVector<vtkIdType> old2new_cell(m_Grid->GetNumberOfCells(), -1);
 458     for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 459       bool ins_cell = false;
 460       vtkIdType type_cell = m_Grid->GetCellType(id_cell);
 461       if (type_cell == VTK_TRIANGLE) ins_cell = true;
 462       if (type_cell == VTK_QUAD)     ins_cell = true;
 463       if (type_cell == VTK_WEDGE)    ins_cell = true;
 464       if (ins_cell) {
 465         vtkIdType N_pts, *pts;
 466         m_Grid->GetCellPoints(id_cell, N_pts, pts);
 467         for (int i = 0; i < N_pts; ++i) {
 468           pts[i] = old2new[pts[i]];
 469           if (pts[i] == -1) {
 470             EG_BUG;
 471           }
 472         }
 473         vtkIdType id_new = vol_grid->InsertNextCell(type_cell, N_pts, pts);
 474         copyCellData(m_Grid, id_cell, vol_grid, id_new);
 475         old2new_cell[id_cell] = id_new;
 476       }
 477     }
 478
 479     // add new cells
 480     vtkIdType id_new_cell;
 481     for (vtkIdType cellId = 0; cellId < Ncells_ng; ++cellId) {
 482       int       pts[8];
 483       vtkIdType new_pts[4];
 484       for (int i = 0; i < 8; ++i) {
 485         pts[i] = 0;
 486       }
 487       Ng_Volume_Element_Type ng_type;
 488       ng_type = Ng_GetVolumeElement(mesh, cellId + 1, pts);
 489       if (ng_type != NG_TET) {
 490         EG_BUG;
 491       }
 492       for (int i = 0; i < 4; ++i) {
 493         if (!ng_surf_node[pts[i]]) {
 494           new_pts[i] = ng2new[pts[i]];
 495         } else {
 496           new_pts[i] = ng2eg[pts[i]];
 497         }
 498       }
 499       if (ng_type == NG_TET) {
 500         vtkIdType tet[4];
 501         tet[0] = new_pts[0];
 502         tet[1] = new_pts[1];
 503         tet[2] = new_pts[3];
 504         tet[3] = new_pts[2];
 505         id_new_cell = vol_grid->InsertNextCell(VTK_TETRA, 4, tet);
 506       } else if (ng_type == NG_PYRAMID) {
 507         EG_ERR_RETURN("pyramids cannot be handled yet");
 508       } else if (ng_type == NG_PRISM) {
 509         EG_ERR_RETURN("prisms cannot be handled yet");
 510       } else {
 511         EG_ERR_RETURN("bug encountered");
 512       }
 513     }
 514     makeCopy(vol_grid, m_Grid);
 515     for (int i = 0; i < trace_cells.size(); ++i) {
 516       if (old2new_cell[trace_cells[i]] == -1) {
 517         EG_BUG;
 518       }
 519       trace_cells[i] = old2new_cell[trace_cells[i]];
 520     }
 521   } else {
 522     Error err;
 523     QString msg = "NETGEN did not succeed.\nPlease check if the surface mesh is oriented correctly";
 524     msg += " (red edges on the outside of the domain)";
 525     err.setType(Error::ExitOperation);
 526     err.setText(msg);
 527     throw err;
 528   }
 529   Ng_DeleteMesh(mesh);
 530   Ng_Exit();
 531   cout << "\n\nNETGEN call finished" << endl;
 532   cout << endl;
 533 }
 534