src/libengrid/tetgenoperation.cpp

   1 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   2 // +                                                                      +
   3 // + This file is part of enGrid.                                         +
   4 // +                                                                      +
   5 // + Copyright 2008-2014 enGits GmbH                                      +
   6 // +                                                                      +
   7 // + enGrid is free software: you can redistribute it and/or modify       +
   8 // + it under the terms of the GNU General Public License as published by +
   9 // + the Free Software Foundation, either version 3 of the License, or    +
  10 // + (at your option) any later version.                                  +
  11 // +                                                                      +
  12 // + enGrid is distributed in the hope that it will be useful,            +
  13 // + but WITHOUT ANY WARRANTY; without even the implied warranty of       +
  14 // + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the        +
  15 // + GNU General Public License for more details.                         +
  16 // +                                                                      +
  17 // + You should have received a copy of the GNU General Public License    +
  18 // + along with enGrid. If not, see <http://www.gnu.org/licenses/>.       +
  19 // +                                                                      +
  20 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  21 //
  22 #include "tetgenoperation.h"
  23 #include "engrid.h"
  24 #include "guimainwindow.h"
  25 #include "optimisenormalvector.h"
  26
  27 #include <QDir>
  28
  29 TetGenOperation::TetGenOperation()
  30 {
  31   // The following code will be used once the option to start TetGen in a
  32   // separate process has been implemented
  33
  34   /*
  35   QDir cwd(GuiMainWindow::pointer()->getCwd());
  36   m_TetGenPath = GuiMainWindow::pointer()->getCwd() + "tetgen";
  37   QDir tetgen_dir(m_TetGenPath);
  38   if (!tetgen_dir.exists()) {
  39     cwd.mkdir("tetgen");
  40   }
  41   */
  42 }
  43
  44 void TetGenOperation::copyToTetGen(tetgenio &tgio, vtkUnstructuredGrid *alt_grid)
  45 {
  46   vtkUnstructuredGrid *grid = alt_grid;
  47   if (!grid) {
  48     grid = m_Grid;
  49   }
  50   MeshPartition part(grid, true);
  51   EG_VTKDCC(vtkIntArray, cell_code, grid, "cell_code");
  52   EG_VTKDCC(vtkIntArray, cell_orgdir, grid, "cell_orgdir");
  53   EG_VTKDCC(vtkIntArray, cell_curdir, grid, "cell_curdir");
  54   EG_VTKDCC(vtkIntArray, cell_voldir, grid, "cell_voldir");
  55   tgio.initialize();
  56   tgio.numberofpoints = grid->GetNumberOfPoints();
  57   tgio.pointlist = new REAL [3*tgio.numberofpoints];
  58   for (vtkIdType id_node = 0; id_node < grid->GetNumberOfPoints(); ++id_node) {
  59     vec3_t x;
  60     grid->GetPoint(id_node, x.data());
  61     for (int i = 0; i < 3; ++i) {
  62       tgio.pointlist[3*id_node + i] = x[i];
  63     }
  64   }
  65   QList<vtkIdType> triangles, tetrahedra;
  66   for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
  67     if      (grid->GetCellType(id_cell) == VTK_TRIANGLE) triangles.append(id_cell);
  68     else if (grid->GetCellType(id_cell) == VTK_TETRA)    tetrahedra.append(id_cell);
  69     else {
  70       EG_BUG;
  71     }
  72   }
  73   tgio.numberoffacets = triangles.size();
  74   tgio.facetmarkerlist = new int [triangles.size()];
  75   tgio.facetlist = new tetgenio::facet [triangles.size()];
  76   {
  77     int i = 0;
  78     foreach (vtkIdType id_cell, triangles) {
  79       m_OrgDir = cell_orgdir->GetValue(id_cell);
  80       m_CurDir = cell_curdir->GetValue(id_cell);
  81       m_VolDir = cell_voldir->GetValue(id_cell);
  82       EG_GET_CELL(id_cell, grid);
  83       tgio.facetlist[i].numberofpolygons = 1;
  84       tgio.facetlist[i].polygonlist = new tetgenio::polygon[1];
  85       tgio.facetlist[i].numberofholes = 0;
  86       tgio.facetlist[i].holelist = NULL;
  87       tetgenio::polygon &poly = tgio.facetlist[i].polygonlist[0];
  88       poly.numberofvertices = 3;
  89       poly.vertexlist = new int[3];
  90       for (int j = 0; j < num_pts; ++j) {
  91         poly.vertexlist[j] = pts[j];
  92       }
  93       tgio.facetmarkerlist[i] = cell_code->GetValue(id_cell);
  94       ++i;
  95     }
  96   }
  97   {
  98     tgio.numberoftetrahedra = tetrahedra.size();
  99     tgio.tetrahedronlist = new int [tetrahedra.size()*4];
 100     int i = 0;
 101     foreach (vtkIdType id_cell, tetrahedra) {
 102       EG_GET_CELL(id_cell, grid);
 103       for (int j = 0; j < num_pts; ++j) {
 104         tgio.tetrahedronlist[i*num_pts + j] = pts[j];
 105       }
 106       ++i;
 107     }
 108
 109     QVector<bool> provide_mesh_resolution(grid->GetNumberOfPoints(), true);
 110     for (vtkIdType id_node = 0; id_node < grid->GetNumberOfPoints(); ++id_node) {
 111       provide_mesh_resolution[id_node] = false;
 112     }
 113     for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
 114       if (isSurface(id_cell, grid)) {
 115         EG_GET_CELL(id_cell, grid);
 116         for (int i = 0; i < num_pts; ++i) {
 117           provide_mesh_resolution[pts[i]] = true;
 118         }
 119       }
 120     }
 121     for (vtkIdType id_node = 0; id_node < grid->GetNumberOfPoints(); ++id_node) {
 122       for (int i = 0; i < part.n2nGSize(id_node); ++i) {
 123         provide_mesh_resolution[part.n2nGG(id_node,i)] = true;
 124       }
 125     }
 126
 127     tgio.numberofpointmtrs = 1;
 128     tgio.pointmtrlist = new REAL [tgio.numberofpoints];
 129     for (vtkIdType id_node = 0; id_node < grid->GetNumberOfPoints(); ++id_node) {
 130       tgio.pointmtrlist[id_node] = 0;
 131       if (provide_mesh_resolution[id_node]) {
 132         vec3_t x;
 133         grid->GetPoint(id_node, x.data());
 134         /// @todo look into this
 135         double l = min(m_MaximalEdgeLength, 0.66*m_ELSManager.minEdgeLength(x)); // apply empirical scaling to match surface resolution
 136         tgio.pointmtrlist[id_node] = l;
 137       }
 138     }
 139   }
 140 }
 141
 142 void TetGenOperation::copyFromTetGen(tetgenio &tgio)
 143 {
 144   // build search tree to correct surface orientation
 145   QVector<vtkIdType> tri_faces;
 146   getAllCellsOfType(VTK_TRIANGLE, tri_faces, m_Grid);
 147   QVector<Triangle> triangles(tri_faces.size());
 148   TriangleTree search_tree;
 149   {
 150     int i = 0;
 151     foreach (vtkIdType id_cell, tri_faces) {
 152       EG_GET_CELL(id_cell, m_Grid);
 153       vec3_t a, b, c;
 154       m_Grid->GetPoint(pts[0], a.data());
 155       m_Grid->GetPoint(pts[1], b.data());
 156       m_Grid->GetPoint(pts[2], c.data());
 157       triangles[i] = Triangle(Point(a[0], a[1], a[2]), Point(b[0], b[1], b[2]), Point(c[0], c[1], c[2]));
 158       ++i;
 159     }
 160   }
 161   search_tree.rebuild(triangles.begin(), triangles.end());
 162   search_tree.accelerate_distance_queries();
 163
 164   EG_VTKSP(vtkUnstructuredGrid, full_grid);
 165   allocateGrid(full_grid, tgio.numberoftetrahedra + tgio.numberoftrifaces, tgio.numberofpoints);
 166   EG_VTKDCC(vtkIntArray, cell_code, full_grid, "cell_code");
 167   EG_VTKDCC(vtkIntArray, cell_orgdir, full_grid, "cell_orgdir");
 168   EG_VTKDCC(vtkIntArray, cell_curdir, full_grid, "cell_curdir");
 169   EG_VTKDCC(vtkIntArray, cell_voldir, full_grid, "cell_voldir");
 170   for (vtkIdType id_node = 0; id_node < tgio.numberofpoints; ++id_node) {
 171     vec3_t x;
 172     for (int i = 0; i < 3; ++i) {
 173       x[i] = tgio.pointlist[3*id_node + i];
 174     }
 175     full_grid->GetPoints()->SetPoint(id_node, x.data());
 176   }
 177   for (int i = 0; i < tgio.numberoftetrahedra; ++i) {
 178     vtkIdType num_pts = 4, pts[4];
 179     for (int j = 0; j < num_pts; ++j) {
 180       pts[j] = tgio.tetrahedronlist[i*num_pts + j];
 181     }
 182     vtkIdType id_cell = full_grid->InsertNextCell(VTK_TETRA, num_pts, pts);
 183     cell_code->SetValue(id_cell, 0);
 184   }
 185   for (int i = 0; i < tgio.numberoftrifaces; ++i) {
 186     vtkIdType num_pts = 3, pts[3];
 187     for (int j = 0; j < num_pts; ++j) {
 188       //pts[2-j] = tgio.trifacelist[i*num_pts + j];
 189       pts[j] = tgio.trifacelist[i*num_pts + j];
 190     }
 191
 192     vec3_t x1, x2, x3;
 193     full_grid->GetPoint(pts[0], x1.data());
 194     full_grid->GetPoint(pts[1], x2.data());
 195     full_grid->GetPoint(pts[2], x3.data());
 196     vec3_t xm = (1.0/3.0)*(x1 + x2 + x3);
 197     vec3_t n1 = GeometryTools::triNormal(x1, x2, x3);
 198     TriangleTree::Point_and_primitive_id cp = search_tree.closest_point_and_primitive(Point(xm[0], xm[1], xm[2]));
 199
 200     Point  pa = cp.second->vertex(0);
 201     vec3_t xa(pa[0], pa[1], pa[2]);
 202     Point  pb = cp.second->vertex(1);
 203     vec3_t xb(pb[0], pb[1], pb[2]);
 204     Point  pc = cp.second->vertex(2);
 205     vec3_t xc(pc[0], pc[1], pc[2]);
 206
 207     vec3_t n2 = GeometryTools::triNormal(xa, xb, xc);
 208     if (n1*n2 < 0) {
 209       swap(pts[0], pts[1]);
 210     }
 211
 212     vtkIdType id_cell = full_grid->InsertNextCell(VTK_TRIANGLE, num_pts, pts);
 213     cell_code->SetValue(id_cell, tgio.trifacemarkerlist[i]);
 214     cell_orgdir->SetValue(id_cell, m_OrgDir);
 215     cell_curdir->SetValue(id_cell, m_CurDir);
 216     cell_voldir->SetValue(id_cell, m_VolDir);
 217   }
 218
 219   QVector<int> keep_cell(full_grid->GetNumberOfCells(), 0);
 220   MeshPartition full_part(full_grid, true);
 221   QList<vtkIdType> seed_cells;
 222   EG_FORALL_CELLS(id_cell, full_grid) {
 223     EG_GET_CELL(id_cell, full_grid);
 224     if (type_cell == VTK_TRIANGLE) {
 225       QVector<QSet<vtkIdType> > cells(3);
 226       for (int i = 0; i < 3; ++i) {
 227         for (int j = 0; j < full_part.n2cGSize(pts[i]); ++j) {
 228           cells[i].insert(full_part.n2cGG(pts[i], j));
 229         }
 230       }
 231       cells[0].intersect(cells[1]);
 232       cells[0].intersect(cells[2]);
 233       vec3_t x_face = cellCentre(full_grid, id_cell);
 234       vec3_t n_face = cellNormal(full_grid, id_cell);
 235       foreach (vtkIdType id_neigh, cells[0]) {
 236         if (full_grid->GetCellType(id_neigh) == VTK_TETRA) {
 237           vec3_t x_cell = cellCentre(full_grid, id_neigh);
 238           if ((x_cell - x_face)*n_face < 0) { // !!!
 239             keep_cell[id_neigh] = 1;
 240           } else {
 241             if (keep_cell[id_neigh] == 0) {
 242               keep_cell[id_neigh] = 2;
 243               seed_cells.append(id_neigh);
 244             }
 245           }
 246         }
 247       }
 248     }
 249   }
 250   while (seed_cells.size() > 0) {
 251     QList<vtkIdType> new_seed_cells;
 252     foreach (vtkIdType id_cell, seed_cells) {
 253       for (int i = 0; i < full_part.c2cGSize(id_cell); ++i) {
 254         vtkIdType id_neigh = full_part.c2cGG(id_cell, i);
 255         if (keep_cell[id_neigh] == 0) {
 256           new_seed_cells.append(id_neigh);
 257           keep_cell[id_neigh] = 2;
 258         }
 259       }
 260     }
 261     seed_cells = new_seed_cells;
 262   }
 263   QList<vtkIdType> keep_cells;
 264   EG_FORALL_CELLS(id_cell, full_grid) {
 265     if (keep_cell[id_cell] < 2) {
 266       keep_cells.append(id_cell);
 267     }
 268   }
 269   cout << "keeping " << keep_cells.size() << " cells" << endl;
 270   MeshPartition keep_part(full_grid);
 271   keep_part.setCells(keep_cells);
 272   keep_part.extractToVtkGrid(m_Grid);
 273   //makeCopy(full_grid, m_Grid);
 274 }
 275
 276 void TetGenOperation::readSettings()
 277 {
 278   QString buffer = GuiMainWindow::pointer()->getXmlSection("engrid/surface/settings");
 279   if(!buffer.isEmpty()) {
 280     QTextStream in(&buffer, QIODevice::ReadOnly);
 281     in >> m_MaximalEdgeLength;
 282     in >> m_MinimalEdgeLength;
 283     in >> m_GrowthFactor;
 284     in >> m_NodesPerQuarterCircle;
 285     int num_bcs;
 286     in >> num_bcs;
 287     int check_state;
 288     for (int i = 0; i < num_bcs; ++i) {
 289       in >> check_state;
 290     }
 291     in >> m_2dFeatureResolution;
 292     in >> m_3dFeatureResolution;
 293   }
 294
 295   m_ELSManager.read();
 296   m_ELSManager.readRules(m_Grid);
 297 }
 298
 299 QString TetGenOperation::qualityText()
 300 {
 301   using namespace GeometryTools;
 302   vec3_t x1(1.0/sqrt(3.0), 0, 0);
 303   vec3_t x2 = rotate(x1, vec3_t(0, 0, 1), deg2rad(60.0));
 304   vec3_t x3 = rotate(x2, vec3_t(0, 0, 1), deg2rad(60.0));
 305   vec3_t x4(0, 0, 1.0);
 306   vec3_t n1 = triNormal(x1, x2, x3);
 307   vec3_t n2 = triNormal(x1, x4, x2);
 308   vec3_t n3 = triNormal(x2, x4, x3);
 309   double alpha = angle(n1, n2);
 310   alpha = min(alpha, angle(n1, n3));
 311   alpha = min(alpha, angle(n2, n3));
 312   QString alpha_txt;
 313   alpha_txt.setNum(rad2deg(alpha));
 314   return QString("1.4") + "/" + alpha_txt;
 315 }
 316
 317 void TetGenOperation::tetgen(QString flags, vtkUnstructuredGrid *background_grid)
 318 {
 319   try {
 320     tetgenio in, out, background;
 321     readSettings();
 322     copyToTetGen(in);
 323     copyToTetGen(background, background_grid);
 324     char *char_flags = new char [flags.size() + 1];
 325     strcpy(char_flags, qPrintable(flags));
 326     tetrahedralize(char_flags, &in, &out, NULL, &background);
 327     copyFromTetGen(out);
 328     delete [] char_flags;
 329   } catch (int x) {
 330     switch (x) {
 331     case 1: // Out of memory.
 332       EG_ERR_RETURN("TetGen error:  Out of memory");
 333       break;
 334     case 2: // Encounter an internal error.
 335       EG_ERR_RETURN("TetGen error:  internal error");
 336       break;
 337     case 3:
 338       EG_ERR_RETURN("TetGen error:  A self-intersection was detected.");
 339       break;
 340     case 4:
 341       EG_ERR_RETURN("TetGen error:  A very small input feature size was detected.");
 342       break;
 343     case 5:
 344       EG_ERR_RETURN("TetGen error:  Two very close input facets were detected.");
 345       break;
 346     case 10:
 347       EG_ERR_RETURN("TetGen error:  An input error was detected.");
 348       EG_BUG;
 349       break;
 350     default:
 351       EG_ERR_RETURN("TetGen error:  unknown error");
 352       break;
 353     }
 354   }
 355 }
 356