src/eliminatesmallbranches.cpp

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