src/libengrid/eliminatesmallbranches.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 "eliminatesmallbranches.h"
  23 #include "createvolumemesh.h"
  24 #include "guimainwindow.h"
  25
  26 #include <QTime>
  27
  28
  29 EliminateSmallBranches::EliminateSmallBranches()
  30 {
  31   m_NumLayers = 0;
  32   m_NumFillLayers = 10;
  33 }
  34
  35 bool EliminateSmallBranches::needsToBeMarked(vtkIdType id_node, int layer)
  36 {
  37   if (m_IsSurfaceNode[id_node]) {
  38     return true;
  39   } else {
  40     if (layer < m_NumLayers) {
  41       for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
  42         if (needsToBeMarked(m_Part.n2nGG(id_node, i), layer+1)) {
  43           return true;
  44         }
  45       }
  46       return false;
  47     }
  48     return false;
  49   }
  50 }
  51
  52 void EliminateSmallBranches::unmarkNode(vtkIdType id_node, int layer)
  53 {
  54   if (layer <= m_NumLayers+2) {
  55     for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
  56       vtkIdType id_cell = m_Part.n2cGG(id_node, i);
  57       if (isVolume(id_cell, m_Grid)) {
  58         m_DeleteCell[id_cell] = false;
  59       }
  60     }
  61     for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
  62       unmarkNode(m_Part.n2nGG(id_node, i), layer+1);
  63     }
  64   }
  65 }
  66
  67 void EliminateSmallBranches::fill(vtkIdType id_cell)
  68 {
  69   if (!m_DeleteCell[id_cell] && !m_MainVolumeCell[id_cell]) {
  70     m_MainVolumeCell[id_cell] = true;
  71     for (int i = 0; i < m_Part.c2cGSize(id_cell); ++i) {
  72       vtkIdType id_neigh_cell = m_Part.c2cGG(id_cell, i);
  73       if (id_neigh_cell != -1) {
  74         m_FillCells.append(id_neigh_cell);
  75         //fill(id_neigh_cell);
  76       }
  77     }
  78   }
  79 }
  80
  81 void EliminateSmallBranches::fillFromLargestVolume()
  82 {
  83   l2g_t  cells = m_Part.getCells();
  84   cout << "filling from largest volume" << endl;
  85   double vol_max = 0;
  86   vtkIdType id_largest_cell = -1;
  87   foreach(vtkIdType id_cell, cells) {
  88     if (isVolume(id_cell, m_Grid)) {
  89       if (m_Grid->GetCellType(id_cell) != VTK_TETRA) {
  90         EG_BUG;
  91       }
  92       double vol = GeometryTools::cellVA(m_Grid, id_cell, true);
  93       if (vol > vol_max && !m_DeleteCell[id_cell]) {
  94         id_largest_cell = id_cell;
  95         vol_max = vol;
  96       }
  97     }
  98   }
  99   if (id_largest_cell == -1) {
 100     EG_BUG;
 101   }
 102   m_MainVolumeCell.fill(false, m_Grid->GetNumberOfCells());
 103   m_FillCells.clear();
 104   m_FillCells.append(id_largest_cell);
 105   while (m_FillCells.size() > 0) {
 106     QList<vtkIdType> fill_cells = m_FillCells;
 107     m_FillCells.clear();
 108     foreach (vtkIdType id_cell, fill_cells) {
 109       fill(id_cell);
 110     }
 111   }
 112 }
 113
 114
 115 void EliminateSmallBranches::fillLayers()
 116 {
 117   QVector<bool> main_volume_cell = m_MainVolumeCell;
 118   for (int i_layer = 0; i_layer < m_NumFillLayers; ++i_layer) {
 119     for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 120       if (m_MainVolumeCell[id_cell]) {
 121         for (int i = 0; i < m_Part.c2cGSize(id_cell); ++i) {
 122           vtkIdType id_neigh_cell = m_Part.c2cGG(id_cell, i);
 123           if (id_neigh_cell != -1) {
 124             if (isVolume(id_neigh_cell, m_Grid)) {
 125               main_volume_cell[id_neigh_cell] = true;
 126             }
 127           }
 128         }
 129       }
 130     }
 131     m_MainVolumeCell = main_volume_cell;
 132   }
 133 }
 134
 135 void EliminateSmallBranches::fillCraters()
 136 {
 137   cout << "trying to fill holes" << endl;
 138   QVector<bool> is_mainvol_node(m_Grid->GetNumberOfPoints(), false);
 139   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 140     if (m_MainVolumeCell[id_cell]) {
 141       vtkIdType N_pts, *pts;
 142       m_Grid->GetCellPoints(id_cell, N_pts, pts);
 143       for (int i = 0; i < N_pts; ++i) {
 144         is_mainvol_node[pts[i]] = true;
 145       }
 146     }
 147   }
 148   QVector<bool> is_mainvol_cell(m_MainVolumeCell.size(), true);
 149   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 150     if (!m_MainVolumeCell[id_cell]) {
 151       vtkIdType N_pts, *pts;
 152       m_Grid->GetCellPoints(id_cell, N_pts, pts);
 153       for (int i = 0; i < N_pts; ++i) {
 154         if (!is_mainvol_node[pts[i]]) {
 155           is_mainvol_cell[id_cell] = false;
 156           break;
 157         }
 158       }
 159     }
 160   }
 161   m_MainVolumeCell = is_mainvol_cell;
 162 }
 163
 164 void EliminateSmallBranches::fixNonManifold()
 165 {
 166   cout << "trying to fix non-manifold edges" << endl;
 167   int N = 1;
 168   int loop = 0;
 169   while (N > 0 && loop < 20) {
 170     N = 0;
 171     cout << loop + 1 << ". sweep" << endl;
 172     QVector<QVector<int> > num_faces(m_Grid->GetNumberOfPoints(), QVector<int>(0));
 173     for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 174       num_faces[id_node].fill(0, m_Part.n2nGSize(id_node));
 175     }
 176     for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 177       if (m_Grid->GetCellType(id_cell) == VTK_TETRA && m_MainVolumeCell[id_cell]) {
 178         for (int i_face = 0; i_face < 4; ++i_face) {
 179           vtkIdType id_neigh = m_Part.c2cGG(id_cell, i_face);
 180           bool proper_neigh = false;
 181           if (id_neigh != -1) {
 182             if (isVolume(id_neigh, m_Grid)) {
 183               if (m_MainVolumeCell[id_neigh]) {
 184                 proper_neigh = true;
 185               }
 186             }
 187           }
 188           if (!proper_neigh) {
 189             QVector<vtkIdType> nds;
 190             getFaceOfCell(m_Grid, id_cell, i_face, nds);
 191             for (int i_node1 = 0; i_node1 < 3; ++i_node1) {
 192               for (int i_node2 = 0; i_node2 < 3; ++i_node2) {
 193                 if (i_node1 != i_node2) {
 194                   int j12 = -1;
 195                   for (int j = 0; j < m_Part.n2nGSize(nds[i_node1]); ++j) {
 196                     if (m_Part.n2nGG(nds[i_node1], j) == nds[i_node2]) {
 197                       j12 = j;
 198                       break;
 199                     }
 200                   }
 201                   if (j12 == -1) {
 202                     EG_BUG;
 203                   }
 204                   ++num_faces[nds[i_node1]][j12];
 205                 }
 206               }
 207             }
 208           }
 209         }
 210       }
 211     }
 212     for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 213       if (m_Grid->GetCellType(id_cell) == VTK_TETRA) {
 214         if (!m_MainVolumeCell[id_cell]) {
 215           for (int i_edge = 0; i_edge < 6; ++i_edge) {
 216             QVector<vtkIdType> nds;
 217             getEdgeOfCell(m_Grid, id_cell, i_edge, nds);
 218             int j12 = 0;
 219             for (int j = 0; j < m_Part.n2nGSize(nds[0]); ++j) {
 220               if (m_Part.n2nGG(nds[0], j) == nds[1]) {
 221                 j12 = j;
 222                 break;
 223               }
 224             }
 225             if (num_faces[nds[0]][j12] != 2 && num_faces[nds[0]][j12] != 0) {
 226               m_MainVolumeCell[id_cell] = true;
 227               ++N;
 228               break;
 229             }
 230           }
 231         }
 232       }
 233     }
 234     cout << "found " << N << " cells" << endl;
 235     ++loop;
 236   }
 237 }
 238
 239 void EliminateSmallBranches::operate()
 240 {
 241   /*
 242   CreateVolumeMesh vol;
 243   vol.setGrid(m_Grid);
 244   vol();
 245   */
 246   //vol();
 247   setAllCells();
 248   l2g_t  cells = m_Part.getCells();
 249   g2l_t _cells = m_Part.getLocalCells();
 250   g2l_t _nodes = m_Part.getLocalNodes();
 251   l2l_t  n2c   = m_Part.getN2C();
 252   l2l_t  c2c   = m_Part.getC2C();
 253   QVector<vtkIdType> faces;
 254   getAllSurfaceCells(faces, m_Grid);
 255   m_DeleteCell.fill(false, m_Grid->GetNumberOfCells());
 256   m_IsSurfaceNode.fill(false, m_Grid->GetNumberOfPoints());
 257   foreach(vtkIdType id_face, faces) {
 258     vtkIdType N_pts, *pts;
 259     m_Grid->GetCellPoints(id_face, N_pts, pts);
 260     for (int i = 0; i < N_pts; ++i) {
 261       m_IsSurfaceNode[pts[i]] = true;
 262     }
 263     m_DeleteCell[id_face] = true;
 264   }
 265
 266   cout << "marking cells to be removed" << endl;
 267   foreach(vtkIdType id_cell, cells) {
 268     vtkIdType N_pts, *pts;
 269     m_Grid->GetCellPoints(id_cell, N_pts, pts);
 270     for (int i = 0; i < N_pts; ++i) {
 271       if (needsToBeMarked(pts[i])) {
 272         m_DeleteCell[id_cell] = true;
 273         break;
 274       }
 275     }
 276   }
 277   fillFromLargestVolume();
 278   fillLayers();
 279   for (int iter = 0; iter < 3; ++iter) {
 280     fillCraters();
 281     fixNonManifold();
 282   }
 283
 284   for (int i = 0; i < m_MainVolumeCell.size(); ++i) {
 285     if (m_MainVolumeCell[i]) {
 286       m_DeleteCell[i] = false;
 287     } else {
 288       m_DeleteCell[i] = true;
 289     }
 290     m_MainVolumeCell[i] = false;
 291   }
 292
 293   cout << "saving existing boundary faces" << endl;
 294   foreach (vtkIdType id_face, faces) {
 295     vtkIdType id_cell = findVolumeCell(m_Grid, id_face, _nodes, cells, _cells, n2c);
 296     if (id_cell != -1) {
 297       if (!m_DeleteCell[id_cell]) {
 298         m_DeleteCell[id_face] = false;
 299       }
 300     }
 301   }
 302
 303   cout << "counting new boundary faces" << endl;
 304   int num_new_faces = 0;
 305   foreach(vtkIdType id_cell, cells) {
 306     if (!m_DeleteCell[id_cell]) {
 307       for (int i = 0; i < m_Part.c2cGSize(id_cell); ++i) {
 308         if (m_DeleteCell[m_Part.c2cGG(id_cell, i)]) {
 309           ++num_new_faces;
 310         }
 311       }
 312     }
 313   }
 314
 315   cout << "creating reduced grid" << endl;
 316   QVector<vtkIdType> old2new(m_Grid->GetNumberOfPoints(), -1);
 317   vtkIdType num_new_nodes = 0;
 318   vtkIdType num_new_cells = 0;
 319   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 320     if (!m_DeleteCell[id_cell]) {
 321       ++num_new_cells;
 322       vtkIdType N_pts, *pts;
 323       m_Grid->GetCellPoints(id_cell, N_pts, pts);
 324       for (int i = 0; i < N_pts; ++i) {
 325         if (old2new[pts[i]] == -1) {
 326           old2new[pts[i]] = num_new_nodes;
 327           ++num_new_nodes;
 328         }
 329       }
 330     }
 331   }
 332   EG_VTKSP(vtkUnstructuredGrid, new_grid);
 333   allocateGrid(new_grid, num_new_cells + num_new_faces, num_new_nodes, true);
 334   EG_VTKDCC(vtkIntArray, bc, new_grid, "cell_code" );
 335   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 336     if (old2new[id_node] != -1) {
 337       vec3_t x;
 338       m_Grid->GetPoint(id_node, x.data());
 339       new_grid->GetPoints()->SetPoint(old2new[id_node], x.data());
 340       copyNodeData(m_Grid, id_node, new_grid, old2new[id_node]);
 341     }
 342   }
 343   EG_VTKDCC(vtkIntArray, cell_orgdir, new_grid, "cell_orgdir");
 344   EG_VTKDCC(vtkIntArray, cell_voldir, new_grid, "cell_voldir");
 345   EG_VTKDCC(vtkIntArray, cell_curdir, new_grid, "cell_curdir");
 346   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 347     if (!m_DeleteCell[id_cell]) {
 348       vtkIdType N_pts, *pts;
 349       m_Grid->GetCellPoints(id_cell, N_pts, pts);
 350       for (int i = 0; i < N_pts; ++i) {
 351         pts[i] = old2new[pts[i]];
 352         if (pts[i] == -1) {
 353           EG_BUG;
 354         }
 355       }
 356       vtkIdType id_new_cell = new_grid->InsertNextCell(m_Grid->GetCellType(id_cell), N_pts, pts);
 357       copyCellData(m_Grid, id_cell, new_grid, id_new_cell);
 358       if (m_Grid->GetCellType(id_cell) == VTK_TETRA) {
 359         for (int i = 0; i < m_Part.c2cGSize(id_cell); ++i) {
 360           if (m_DeleteCell[m_Part.c2cGG(id_cell, i)]) {
 361             vtkIdType face_pts[3];
 362             if (i == 0) {
 363               face_pts[0] = pts[2];
 364               face_pts[1] = pts[1];
 365               face_pts[2] = pts[0];
 366             } else if (i == 1) {
 367               face_pts[0] = pts[1];
 368               face_pts[1] = pts[3];
 369               face_pts[2] = pts[0];
 370             } else if (i == 2) {
 371               face_pts[0] = pts[3];
 372               face_pts[1] = pts[2];
 373               face_pts[2] = pts[0];
 374             } else if (i == 3) {
 375               face_pts[0] = pts[2];
 376               face_pts[1] = pts[3];
 377               face_pts[2] = pts[1];
 378             }
 379             vtkIdType id_new_cell = new_grid->InsertNextCell(VTK_TRIANGLE, 3, face_pts);
 380             copyCellData(m_Grid, id_cell, new_grid, id_new_cell);
 381             bc->SetValue(id_new_cell, 99);
 382             cell_orgdir->SetValue(id_new_cell, 0);
 383             cell_voldir->SetValue(id_new_cell, 0);
 384             cell_curdir->SetValue(id_new_cell, 0);
 385           }
 386         }
 387       }
 388     }
 389   }
 390   makeCopy(new_grid, m_Grid);
 391   GuiMainWindow::pointer()->updateBoundaryCodes(true);
 392 }