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