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