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