src/libengrid/createhexcore.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 "createhexcore.h"
  23 #include "guimainwindow.h"
  24 #include "pointfinder.h"
  25
  26 CreateHexCore::CreateHexCore(vec3_t x1, vec3_t x2, vec3_t xi, int num_inital_refinement_levels)
  27 {
  28   m_X1 = x1;
  29   m_X2 = x2;
  30   m_Xi = xi;
  31   m_NumInitialRefinementLevels = num_inital_refinement_levels;
  32   m_NumBreakOutLayers = 1;
  33 }
  34
  35 void CreateHexCore::refineOctree()
  36 {
  37   m_Octree.resetRefineMarks();
  38   m_Octree.setSmoothTransitionOn();
  39
  40   EG_VTKDCN(vtkDoubleArray, cl, m_Grid, "node_meshdensity_desired");
  41   for (int i = 0; i < m_NumInitialRefinementLevels; ++i) {
  42     m_Octree.markAllToRefine();
  43     m_Octree.refineAll();
  44   }
  45   do {
  46     for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
  47       vtkIdType id_surf = -1;
  48       for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
  49         vtkIdType id_cell = m_Part.n2cGG(id_node, i);
  50         if (isSurface(id_cell, m_Grid)) {
  51           id_surf = id_node;
  52           break;
  53         } else {
  54           vtkIdType cell_type = m_Grid->GetCellType(id_cell);
  55           if (cell_type == VTK_WEDGE) {
  56             EG_GET_CELL(id_cell, m_Grid);
  57             if      (pts[3] == id_node) id_surf = pts[0];
  58             else if (pts[4] == id_node) id_surf = pts[1];
  59             else if (pts[5] == id_node) id_surf = pts[2];
  60           }
  61         }
  62       }
  63       if (id_surf != -1) {
  64         vec3_t x;
  65         m_Grid->GetPoint(id_node, x.data());
  66         int i_otcell = m_Octree.findCell(x);
  67         double h = m_Octree.getDx(i_otcell);
  68         h = max(h, m_Octree.getDy(i_otcell));
  69         h = max(h, m_Octree.getDz(i_otcell));
  70         if (h > cl->GetValue(id_surf)) {
  71           m_Octree.markToRefine(i_otcell);
  72         }
  73       }
  74     }
  75   } while (m_Octree.refineAll());
  76 }
  77
  78 void CreateHexCore::transferOctreeGrid()
  79 {
  80   QVector<bool> delete_node(m_Octree.getNumNodes(), false);
  81   QVector<bool> delete_cell(m_Octree.getNumCells(), false);
  82   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
  83     vec3_t x;
  84     m_Grid->GetPoint(id_node, x.data());
  85     int i_cell = m_Octree.findCell(x);
  86     if (m_Octree.hasChildren(i_cell)) {
  87       EG_BUG;
  88     }
  89     delete_cell[i_cell] = true;
  90     for (int i = 0; i < 8; ++i) {
  91       delete_node[m_Octree.getNode(i_cell, i)] = true;
  92     }
  93   }
  94
  95   for (int layer = 0; layer < m_NumBreakOutLayers; ++layer) {
  96     for (int i = 0; i < m_Octree.getNumCells(); ++i) {
  97       if (delete_cell[i] && !m_Octree.hasChildren(i)) {
  98         for (int j = 0; j < 8; ++j) {
  99           delete_node[m_Octree.getNode(i,j)] = true;
 100         }
 101       }
 102     }
 103     for (int i = 0; i < m_Octree.getNumCells(); ++i) {
 104       if (!m_Octree.hasChildren(i)) {
 105         for (int j = 0; j < 8; ++j) {
 106           if (delete_node[m_Octree.getNode(i,j)]) {
 107             delete_cell[i] = true;
 108             break;
 109           }
 110         }
 111       }
 112     }
 113   }
 114   EG_VTKSP(vtkUnstructuredGrid, otgrid);
 115   m_Octree.toVtkGridPolyhedral(otgrid, true);
 116   //m_Octree.toVtkGridConforming(otgrid, true);
 117   MeshPartition add_part(otgrid);
 118   QList<vtkIdType> add_cells;
 119   for (vtkIdType id_cell = 0; id_cell < otgrid->GetNumberOfCells(); ++id_cell) {
 120     int i_cell = m_Octree.findCell(cellCentre(otgrid, id_cell));
 121     if (!delete_cell[i_cell]) {
 122       add_cells.append(id_cell);
 123     }
 124   }
 125   add_part.setCells(add_cells);
 126   m_Part.addPartition(add_part);
 127   m_Part.setAllCells();
 128   deleteOutside(m_Grid);
 129 }
 130
 131 void CreateHexCore::deleteOutside(vtkUnstructuredGrid *grid)
 132 {
 133   MeshPartition part(grid, true);
 134   QVector<bool> is_inside(grid->GetNumberOfCells(), false);
 135   vtkIdType id_start = -1;
 136   double dmin = 1e99;
 137   for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
 138     if (isVolume(id_cell, grid)) {
 139       vec3_t x = cellCentre(grid, id_cell);
 140       double d = (x - m_Xi).abs();
 141       if (d < dmin) {
 142         dmin = d;
 143         id_start = id_cell;
 144       }
 145     }
 146   }
 147   if (id_start == -1) {
 148     EG_BUG;
 149   }
 150   is_inside[id_start] = true;
 151   bool added = true;
 152   while (added) {
 153     added = false;
 154     for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
 155       if (is_inside[id_cell]) {
 156         for (int j = 0; j < part.c2cGSize(id_cell); ++j) {
 157           vtkIdType id_neigh = part.c2cGG(id_cell, j);
 158           if (id_neigh >= 0) {
 159             if (!is_inside[id_neigh]) {
 160               is_inside[id_neigh] = true;
 161               added = true;
 162             }
 163           }
 164         }
 165       }
 166     }
 167   }
 168   int N = 0;
 169   for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
 170     if (isSurface(id_cell, grid)) {
 171       is_inside[id_cell] = true;
 172     }
 173     if (is_inside[id_cell]) {
 174       ++N;
 175     }
 176   }
 177   QVector<vtkIdType> cls(N);
 178   N = 0;
 179   for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
 180     if (is_inside[id_cell]) {
 181       cls[N] = id_cell;
 182       ++N;
 183     }
 184   }
 185   EG_VTKSP(vtkUnstructuredGrid, return_grid);
 186   makeCopy(grid, return_grid, cls);
 187   makeCopy(return_grid, grid);
 188 }
 189
 190 void CreateHexCore::createBoundaryFaces()
 191 {
 192   EG_VTKSP(vtkUnstructuredGrid, new_grid);
 193
 194   // find all polygons which need to be triangulated and collect the new triangles
 195   m_Part.setAllCells();
 196   QList<QVector<vtkIdType> > new_triangles;
 197   QVector<bool> adapt_cell(m_Grid->GetNumberOfCells(), false);
 198   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 199     vec3_t xc = cellCentre(m_Grid, id_cell);
 200     if (isVolume(id_cell, m_Grid)) {
 201       for (int i = 0; i < m_Part.c2cGSize(id_cell); ++i) {
 202         if (m_Part.c2cGG(id_cell, i) == -1) {
 203           QVector<vtkIdType> face;
 204           getFaceOfCell(m_Grid, id_cell, i, face);
 205           QVector<QVector<vtkIdType> > triangles;
 206           triangulatePolygon(m_Grid, face, triangles);
 207           foreach (QVector<vtkIdType> triangle, triangles) {
 208             vec3_t x1, x2, x3;
 209             m_Grid->GetPoint(triangle[0], x1.data());
 210             m_Grid->GetPoint(triangle[1], x2.data());
 211             m_Grid->GetPoint(triangle[2], x3.data());
 212             vec3_t xt = (1.0/3.0)*(x1 + x2 + x3);
 213             vec3_t nt = GeometryTools::triNormal(x1, x2, x3);
 214             if (nt*(xt - xc) < 0) {
 215               swap(triangle[0], triangle[1]);
 216             }
 217             new_triangles.append(triangle);
 218           }
 219           if (face.size() > 3) {
 220             adapt_cell[id_cell] = true;
 221           }
 222         }
 223       }
 224     }
 225   }
 226
 227   // allocate memory for the new grid
 228   allocateGrid(new_grid, m_Grid->GetNumberOfCells() + new_triangles.size(), m_Grid->GetNumberOfPoints());
 229
 230   // copy nodes
 231   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 232     vec3_t x;
 233     m_Grid->GetPoint(id_node, x.data());
 234     new_grid->GetPoints()->SetPoint(id_node, x.data());
 235   }
 236
 237   // copy existing cells and update if required
 238   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 239     vtkIdType id_new_cell;
 240     if (!adapt_cell[id_cell]) {
 241       id_new_cell = copyCell(m_Grid, id_cell, new_grid);
 242     } else {
 243       QList<QVector<vtkIdType> > faces_of_cell;
 244       for (int i = 0; i < m_Part.c2cGSize(id_cell); ++i) {
 245         QVector<vtkIdType> face;
 246         getFaceOfCell(m_Grid, id_cell, i, face);
 247         if (m_Part.c2cGG(id_cell, i) == -1) {
 248           QVector<QVector<vtkIdType> > triangles;
 249           triangulatePolygon(m_Grid, face, triangles);
 250           foreach (QVector<vtkIdType> triangle, triangles) {
 251             faces_of_cell.append(triangle);
 252           }
 253         } else {
 254           faces_of_cell.append(face);
 255         }
 256       }
 257       EG_VTKSP(vtkIdList, stream);
 258       int stream_size = 1;
 259       foreach (QVector<vtkIdType> face, faces_of_cell) {
 260         stream_size += face.size() + 1;
 261       }
 262       stream->SetNumberOfIds(stream_size);
 263       int id = 0;
 264       stream->SetId(id++, faces_of_cell.size());
 265       foreach (QVector<vtkIdType> face, faces_of_cell) {
 266         stream->SetId(id++, face.size());
 267         foreach (vtkIdType id_node, face) {
 268           stream->SetId(id++, id_node);
 269         }
 270       }
 271       id_new_cell = new_grid->InsertNextCell(VTK_POLYHEDRON, stream);
 272     }
 273     copyCellData(m_Grid, id_cell, new_grid, id_new_cell);
 274   }
 275
 276   // determine new boundary code
 277   EG_VTKDCC(vtkIntArray, cell_code, new_grid, "cell_code");
 278   QSet<int> bcs = GuiMainWindow::pointer()->getAllBoundaryCodes();
 279   int bc_new = 1;
 280   foreach (int bc, bcs) {
 281     bc_new = max(bc, bc_new);
 282   }
 283   ++bc_new;
 284
 285   // create triangular boundary faces
 286   EG_VTKDCC(vtkIntArray, cell_orgdir, new_grid, "cell_orgdir");
 287   EG_VTKDCC(vtkIntArray, cell_curdir, new_grid, "cell_curdir");
 288   EG_VTKDCC(vtkIntArray, cell_voldir, new_grid, "cell_voldir");
 289   GuiMainWindow::pointer()->setBC(bc_new, BoundaryCondition("HexCore", "unknown", bc_new));
 290   foreach (QVector<vtkIdType> face, new_triangles) {
 291     vtkIdType id_new_face;
 292     if (face.size() == 3) {
 293       id_new_face = new_grid->InsertNextCell(VTK_TRIANGLE, 3, face.data());
 294     } else if (face.size() == 4) {
 295       id_new_face = new_grid->InsertNextCell(VTK_QUAD, 4, face.data());
 296     } else {
 297       id_new_face = new_grid->InsertNextCell(VTK_POLYGON, face.size(), face.data());
 298     }
 299     cell_code->SetValue(id_new_face, bc_new);
 300     cell_orgdir->SetValue(id_new_face, 0);
 301     cell_curdir->SetValue(id_new_face, 0);
 302     cell_voldir->SetValue(id_new_face, 0);
 303   }
 304
 305   makeCopy(new_grid, m_Grid);
 306 }
 307
 308 void CreateHexCore::operate()
 309 {
 310   m_Octree.setBounds(m_X1, m_X2, 1, 1, 1);
 311   refineOctree();
 312   EG_VTKSP(vtkUnstructuredGrid, otgrid);
 313   transferOctreeGrid();
 314   createBoundaryFaces();
 315   updateCellIndex(m_Grid);
 316   GuiMainWindow::pointer()->updateBoundaryCodes(true);
 317 }