src/seedsimpleprismaticlayer.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 #include "seedsimpleprismaticlayer.h"
  24 #include <vtkIdList.h>
  25
  26 SeedSimplePrismaticLayer::SeedSimplePrismaticLayer()
  27 {
  28   layer_g  = 0.0;
  29   layer_dg = 1.0;
  30 }
  31
  32 void SeedSimplePrismaticLayer::setLayerCells(const QVector<vtkIdType> &cells)
  33 {
  34   layer_cells.resize(cells.size());
  35   qCopy(cells.begin(), cells.end(), layer_cells.begin());
  36 }
  37
  38 void SeedSimplePrismaticLayer::getLayerCells(QVector<vtkIdType> &cells)
  39 {
  40   cells.resize(layer_cells.size());
  41   qCopy(layer_cells.begin(), layer_cells.end(), cells.begin());
  42 }
  43
  44 void SeedSimplePrismaticLayer::prepareLayer()
  45 {
  46   grid->BuildLinks();
  47   EG_VTKSP(vtkIdList, nds);
  48   EG_VTKSP(vtkIdList, cls);
  49   EG_VTKDCN(vtkIntArray, node_layer,  grid, "node_layer");
  50   vol_cells.fill(-1, layer_cells.size());
  51   faces.resize(layer_cells.size());
  52   N_new_cells  = 0;
  53   QSet<vtkIdType> new_points;
  54   for (int i_layer_cell = 0; i_layer_cell < layer_cells.size(); ++i_layer_cell) {
  55     vtkIdType  id_cell   = layer_cells[i_layer_cell];
  56     vtkIdType  type_cell = grid->GetCellType(id_cell);
  57     vtkIdType *pts;
  58     vtkIdType  Npts;
  59     grid->GetCellPoints(id_cell, Npts, pts);
  60     if (type_cell == VTK_TRIANGLE) {
  61       nds->Reset();
  62       faces[i_layer_cell].resize(Npts);
  63       for (int i_pts = 0; i_pts < Npts; ++i_pts) {
  64         faces[i_layer_cell][Npts - 1 - i_pts] = pts[i_pts];
  65         nds->InsertNextId(pts[i_pts]);
  66         new_points.insert(pts[i_pts]);
  67       }
  68       grid->GetCellNeighbors(id_cell, nds, cls);
  69       for (int i_cls = 0; i_cls < cls->GetNumberOfIds(); ++i_cls) {
  70         if (isVolume(cls->GetId(i_cls), grid)) {
  71           if (cls->GetId(i_cls) != id_cell) {
  72             vol_cells[i_layer_cell] = cls->GetId(i_cls);
  73           }
  74         }
  75       }
  76       N_new_cells  += 1;
  77     } else if (type_cell == VTK_WEDGE) {
  78       nds->Reset();
  79       faces[i_layer_cell].resize(3);
  80       for (int i_pts = 0; i_pts < 3; ++i_pts) {
  81         faces[i_layer_cell][2 - i_pts] = pts[i_pts+3];
  82         //faces[i_layer_cell][i_pts] = pts[i_pts+3];
  83         nds->InsertNextId(pts[i_pts+3]);
  84         new_points.insert(pts[i_pts+3]);
  85       }
  86       grid->GetCellNeighbors(id_cell, nds, cls);
  87       for (int i_cls = 0; i_cls < cls->GetNumberOfIds(); ++i_cls) {
  88         if (isVolume(cls->GetId(i_cls), grid)) {
  89           if (cls->GetId(i_cls) != id_cell) {
  90             vol_cells[i_layer_cell] = cls->GetId(i_cls);
  91           }
  92         }
  93       }
  94       N_new_cells  += 1;
  95     }
  96   }
  97   new_layer = -1;
  98   old_layer = -1;
  99   if (faces.size() > 0) {
 100     old_layer = node_layer->GetValue(faces[0][0]);
 101     new_layer = old_layer + 1;
 102   }
 103   N_new_cells += countBoundaryElements();
 104   N_new_points = new_points.size();
 105 }
 106
 107 int SeedSimplePrismaticLayer::countBoundaryElements()
 108 {
 109   int N = 0;
 110   QVector<QSet<int> > n2f(grid->GetNumberOfPoints());
 111   for (int i_faces = 0; i_faces < faces.size(); ++i_faces) {
 112     for (int j_faces = 0; j_faces < faces[i_faces].size(); ++j_faces) {
 113       n2f[faces[i_faces][j_faces]].insert(i_faces);
 114     }
 115   }
 116   EG_VTKDCN(vtkIntArray, node_layer,  grid, "node_layer");
 117   EG_VTKDCN(vtkIntArray, node_status, grid, "node_status");
 118   for (int i_faces = 0; i_faces < faces.size(); ++i_faces) {
 119     for (int j_faces = 0; j_faces < faces[i_faces].size(); ++j_faces) {
 120       vtkIdType p1 = faces[i_faces][j_faces];
 121       vtkIdType p2 = faces[i_faces][0];
 122       if (j_faces < faces[i_faces].size() - 1) {
 123         p2 = faces[i_faces][j_faces+1];
 124       }
 125       bool consider_edge = false;
 126       if ((node_status->GetValue(p1) & 2) && (node_status->GetValue(p2) & 2)) {
 127         consider_edge = true;
 128       }
 129       if ((node_layer->GetValue(p1) == 0) && (node_layer->GetValue(p2) == 0)) {
 130         consider_edge = true;
 131       }
 132       if (consider_edge) {
 133         QSet<int> faces_on_edge;
 134         qcontIntersection(n2f[p1], n2f[p2], faces_on_edge);
 135         if (faces_on_edge.size() == 0) {
 136           EG_BUG;
 137         } else if (faces_on_edge.size() == 1) {
 138           ++N;
 139         } else if (faces_on_edge.size() > 2) {
 140           EG_BUG;
 141         }
 142       }
 143     }
 144   }
 145   return N;
 146 }
 147
 148 void SeedSimplePrismaticLayer::createBoundaryElements(vtkUnstructuredGrid *new_grid)
 149 {
 150   QVector<QSet<int> > n2f(grid->GetNumberOfPoints());
 151   for (int i_faces = 0; i_faces < faces.size(); ++i_faces) {
 152     for (int j_faces = 0; j_faces < faces[i_faces].size(); ++j_faces) {
 153       n2f[faces[i_faces][j_faces]].insert(i_faces);
 154     }
 155   }
 156   QVector<vtkIdType>  bcells;
 157   QVector<vtkIdType>  bnodes;
 158   QVector<int>        _bnodes;
 159   QVector<QSet<int> > bn2bc;
 160   getAllSurfaceCells(bcells, new_grid);
 161   getNodesFromCells(bcells, bnodes, new_grid);
 162   createNodeMapping(bnodes, _bnodes, new_grid);
 163   createNodeToCell(bcells, bnodes, _bnodes, bn2bc, new_grid);
 164   EG_VTKDCC(vtkIntArray, cell_code,   new_grid, "cell_code");
 165   EG_VTKDCN(vtkIntArray, node_layer,  new_grid, "node_layer");
 166   EG_VTKDCN(vtkIntArray, node_status, new_grid, "node_status");
 167   EG_VTKDCC(vtkIntArray, cell_orgdir, new_grid, "cell_orgdir");
 168   EG_VTKDCC(vtkIntArray, cell_voldir, new_grid, "cell_voldir");
 169   EG_VTKDCC(vtkIntArray, cell_curdir, new_grid, "cell_curdir");
 170   for (int i_faces = 0; i_faces < faces.size(); ++i_faces) {
 171     for (int j_faces = 0; j_faces < faces[i_faces].size(); ++j_faces) {
 172       vtkIdType p1 = faces[i_faces][j_faces];
 173       vtkIdType p2 = faces[i_faces][0];
 174       if (j_faces < faces[i_faces].size() - 1) {
 175         p2 = faces[i_faces][j_faces+1];
 176       }
 177       bool consider_edge = false;
 178       if ((node_status->GetValue(p1) & 2) && (node_status->GetValue(p2) & 2)) {
 179         consider_edge = true;
 180       }
 181       if ((node_layer->GetValue(p1) == 0) && (node_layer->GetValue(p2) == 0)) {
 182         consider_edge = true;
 183       }
 184       if (consider_edge) {
 185         QSet<int> faces_on_edge;
 186         qcontIntersection(n2f[p1], n2f[p2], faces_on_edge);
 187         if (faces_on_edge.size() == 0) {
 188           EG_BUG;
 189         } else if (faces_on_edge.size() == 1) {
 190           vtkIdType pts[4] = { p1, p2, old2new[p2], old2new[p1] };
 191           vtkIdType id_new_bcell = new_grid->InsertNextCell(VTK_QUAD ,4, pts);
 192           int bc = 9999;
 193           QSet<int> bc1, bc2, bc3;
 194           int org_dir = -99;
 195           int cur_dir = -99;
 196           int vol_dir = -99;
 197           if (_bnodes[old2new[p1]] != -1) {
 198             foreach (int i_bcells, bn2bc[_bnodes[old2new[p1]]]) {
 199               if (org_dir == -99) {
 200                 org_dir = cell_orgdir->GetValue(bcells[i_bcells]);
 201               } else if (cell_orgdir->GetValue(bcells[i_bcells]) != org_dir) {
 202                 EG_BUG;
 203               }
 204               if (cur_dir == -99) {
 205                 cur_dir = cell_curdir->GetValue(bcells[i_bcells]);
 206               } else if (cell_curdir->GetValue(bcells[i_bcells]) != cur_dir) {
 207                 EG_BUG;
 208               }
 209               if (vol_dir == -99) {
 210                 vol_dir = cell_voldir->GetValue(bcells[i_bcells]);
 211               } else if (cell_voldir->GetValue(bcells[i_bcells]) != vol_dir) {
 212                 EG_BUG;
 213               }
 214               if (!m_BoundaryCodes.contains(cell_code->GetValue(bcells[i_bcells]))) {
 215                 bc1.insert(cell_code->GetValue(bcells[i_bcells]));
 216               }
 217             }
 218           }
 219           if (_bnodes[old2new[p2]] != -1) {
 220             foreach (int i_bcells, bn2bc[_bnodes[old2new[p2]]]) {
 221               if (org_dir == -99) {
 222                 org_dir = cell_orgdir->GetValue(bcells[i_bcells]);
 223               } else if (cell_orgdir->GetValue(bcells[i_bcells]) != org_dir) {
 224                 EG_BUG;
 225               }
 226               if (cur_dir == -99) {
 227                 cur_dir = cell_curdir->GetValue(bcells[i_bcells]);
 228               } else if (cell_curdir->GetValue(bcells[i_bcells]) != cur_dir) {
 229                 EG_BUG;
 230               }
 231               if (vol_dir == -99) {
 232                 vol_dir = cell_voldir->GetValue(bcells[i_bcells]);
 233               } else if (cell_voldir->GetValue(bcells[i_bcells]) != vol_dir) {
 234                 EG_BUG;
 235               }
 236 /*              if (!boundary_codes.contains(cell_code->GetValue(bcells[i_bcells]))) {
 237                 bc1.insert(cell_code->GetValue(bcells[i_bcells]));
 238               }*/
 239               if (!m_BoundaryCodes.contains(cell_code->GetValue(bcells[i_bcells]))) {
 240                 bc2.insert(cell_code->GetValue(bcells[i_bcells]));
 241               }
 242             }
 243           }
 244           qcontIntersection(bc1, bc2, bc3);
 245           if (bc3.size() == 1) {
 246             bc = *(bc3.begin());
 247           } else {
 248             bc = 9999;
 249           //EG_BUG;
 250           }
 251           cell_code->SetValue(id_new_bcell, bc);
 252           cell_orgdir->SetValue(id_new_bcell, org_dir);
 253           cell_voldir->SetValue(id_new_bcell, vol_dir);
 254           cell_curdir->SetValue(id_new_bcell, cur_dir);
 255           node_status->SetValue(p1, node_status->GetValue(p1) | 2);
 256           node_status->SetValue(old2new[p1], node_status->GetValue(p1) | 2);
 257           node_status->SetValue(p2, node_status->GetValue(p2) | 2);
 258           node_status->SetValue(old2new[p2], node_status->GetValue(p2) | 2);
 259         } else if (faces_on_edge.size() > 2) {
 260           EG_BUG;
 261         }
 262       }
 263     }
 264   }
 265 }
 266
 267 void SeedSimplePrismaticLayer::operate()
 268 {
 269   cout << "seeding prismatic layer" << endl;
 270   prepareLayer();
 271   EG_VTKSP(vtkUnstructuredGrid, new_grid);
 272   allocateGrid(new_grid, grid->GetNumberOfCells() + N_new_cells, grid->GetNumberOfPoints() + N_new_points);
 273   vtkIdType id_new_node = 0;
 274   EG_VTKDCN(vtkIntArray, node_status_old, grid,     "node_status");
 275   EG_VTKDCN(vtkIntArray, node_status_new, new_grid, "node_status");
 276   EG_VTKDCN(vtkIntArray, node_layer_old,  grid,     "node_layer");
 277   EG_VTKDCN(vtkIntArray, node_layer_new,  new_grid, "node_layer");
 278   EG_VTKDCC(vtkIntArray, bc,              grid,     "cell_code");
 279
 280   l2l_t  n2c   = getPartN2C();
 281   g2l_t _nodes = getPartLocalNodes();
 282
 283   QVector<QSet<int> > n2f(grid->GetNumberOfPoints());
 284   for (int i_faces = 0; i_faces < faces.size(); ++i_faces) {
 285     for (int j_faces = 0; j_faces < faces[i_faces].size(); ++j_faces) {
 286       n2f[faces[i_faces][j_faces]].insert(i_faces);
 287     }
 288   }
 289
 290   // copy old grid nodes to the new grid
 291   old2new.resize(grid->GetNumberOfPoints());
 292   for (vtkIdType id_node = 0; id_node < grid->GetNumberOfPoints(); ++id_node) {
 293     vec3_t x;
 294     grid->GetPoints()->GetPoint(id_node, x.data());
 295     new_grid->GetPoints()->SetPoint(id_new_node, x.data());
 296     copyNodeData(grid, id_node, new_grid, id_new_node);
 297     old2new[id_node] = id_new_node;
 298     ++id_new_node;
 299   }
 300
 301   QSet<vtkIdType> split_nodes;
 302   for (int i_faces = 0; i_faces < faces.size(); ++i_faces) {
 303     for (int j_faces = 0; j_faces < faces[i_faces].size(); ++j_faces) {
 304       split_nodes.insert(faces[i_faces][j_faces]);
 305     }
 306   }
 307   qDebug() << split_nodes.size();
 308   QVector<bool> is_split_node(grid->GetNumberOfPoints(), false);
 309   foreach (vtkIdType id_node, split_nodes) {
 310     vec3_t x;
 311     grid->GetPoints()->GetPoint(id_node, x.data());
 312
 313
 314     if (n2f[id_node].size() > 0) {
 315       vec3_t n(0,0,0);
 316       double L = 1e99;
 317       foreach (int i_faces, n2f[id_node]) {
 318         vec3_t a,b;
 319         if (faces[i_faces][0] == id_node) {
 320           grid->GetPoint(faces[i_faces][1],a.data());
 321           grid->GetPoint(faces[i_faces][2],b.data());
 322         }
 323         if (faces[i_faces][1] == id_node) {
 324           grid->GetPoint(faces[i_faces][2],a.data());
 325           grid->GetPoint(faces[i_faces][0],b.data());
 326         }
 327         if (faces[i_faces][2] == id_node) {
 328           grid->GetPoint(faces[i_faces][0],a.data());
 329           grid->GetPoint(faces[i_faces][1],b.data());
 330         }
 331         a -= x;
 332         b -= x;
 333         L = min(a.abs(),L);
 334         L = min(b.abs(),L);
 335         a.normalise();
 336         b.normalise();
 337         vec3_t nf = a.cross(b);
 338         nf.normalise();
 339         double alpha = acos(a*b);
 340         n += alpha*nf;
 341       }
 342       for (int i_boundary_correction = 0; i_boundary_correction < 20; ++i_boundary_correction) {
 343         foreach (vtkIdType id_cell, n2c[_nodes[id_node]]) {
 344           if (isSurface(id_cell, grid)) {
 345             if (!m_BoundaryCodes.contains(bc->GetValue(id_cell))) {
 346               double A = GeometryTools::cellVA(grid, id_cell);
 347               if (A > 1e-60) {
 348                 vec3_t nf = GeometryTools::cellNormal(grid, id_cell);
 349                 nf.normalise();
 350                 n -= (nf*n)*nf;
 351               }
 352             }
 353           }
 354         }
 355       }
 356       n.normalise();
 357       x += 0.01*L*n;
 358     } else {
 359       EG_BUG;
 360     }
 361
 362     new_grid->GetPoints()->SetPoint(id_new_node, x.data());
 363     old2new[id_node] = id_new_node;
 364     node_status_new->SetValue(id_new_node, node_status_old->GetValue(id_node));
 365     node_layer_new->SetValue(id_new_node, node_layer_old->GetValue(id_node) + 1);
 366     ++id_new_node;
 367     is_split_node[id_node] = true;
 368   }
 369   QVector<bool> needs_correction(grid->GetNumberOfCells(), false);
 370   for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
 371     vtkIdType type_cell = grid->GetCellType(id_cell);
 372     if ((type_cell == VTK_TRIANGLE) || (type_cell == VTK_TETRA)) {
 373       vtkIdType *pts, N_pts;
 374       grid->GetCellPoints(id_cell, N_pts, pts);
 375       bool split = false;
 376       for (int i_pts = 0; i_pts < N_pts; ++i_pts) {
 377         if (is_split_node[pts[i_pts]]) {
 378           split = true;
 379           if (node_layer_old->GetValue(pts[i_pts]) != old_layer) {
 380             EG_BUG;
 381           }
 382         }
 383       }
 384       if (split) {
 385         if (type_cell == VTK_TETRA) {
 386           needs_correction[id_cell] = true;
 387         } else {
 388           bool f = false;
 389           if (old_layer > 0) {
 390             for (int i_pts = 0; i_pts < N_pts; ++i_pts) {
 391               if (node_layer_old->GetValue(pts[i_pts]) == old_layer - 1) {
 392                 f = true;
 393               }
 394             }
 395           }
 396           if (!f) {
 397             needs_correction[id_cell] = true;
 398           } else {
 399             cout << "dodgy face: " << id_cell << "   ";
 400             for (int i_pts = 0; i_pts < N_pts; ++i_pts) {
 401               cout << "(" << pts[i_pts] << "," << node_layer_old->GetValue(pts[i_pts]) << ") ";
 402             }
 403             cout << endl;
 404             for (int i_pts = 0; i_pts < N_pts; ++i_pts) {
 405               if (is_split_node[pts[i_pts]]) {
 406                 vec3_t x;
 407                 grid->GetPoint(pts[i_pts], x.data());
 408                 cout << "split node: " << pts[i_pts] << "  " << x << endl;
 409               }
 410             }
 411           }
 412         }
 413       }
 414     }
 415   }
 416   foreach (vtkIdType id_cell, layer_cells) {
 417     needs_correction[id_cell] = false;
 418   }
 419
 420   QVector<bool> is_in_vol_cells(grid->GetNumberOfCells(), false);
 421   for (int i_faces = 0; i_faces < faces.size(); ++i_faces) {
 422     vtkIdType id_vol_cell = vol_cells[i_faces];
 423     if (id_vol_cell >= 0) {
 424       is_in_vol_cells[id_vol_cell] = true;
 425       vtkIdType type_vol_cell = grid->GetCellType(id_vol_cell);
 426       if (type_vol_cell == VTK_TETRA) {
 427         vtkIdType *pts, N_pts, p[6];
 428         grid->GetCellPoints(id_vol_cell, N_pts, pts);
 429         p[0] = faces[i_faces][0];
 430         p[1] = faces[i_faces][1];
 431         p[2] = faces[i_faces][2];
 432         p[3] = old2new[p[0]];
 433         p[4] = old2new[p[1]];
 434         p[5] = old2new[p[2]];
 435         {
 436           vtkIdType pts[6] = { p[2], p[1], p[0], p[5], p[4], p[3] };
 437           layer_cells[i_faces] = new_grid->InsertNextCell(VTK_WEDGE, 6, pts);
 438         }
 439       } else {
 440         qDebug() << type_vol_cell;
 441         EG_BUG;
 442       }
 443     }
 444   }
 445
 446 //   writeGrid(new_grid, "pre-cellcopy");
 447
 448   // copy old cells to the new grid
 449   for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
 450     vtkIdType *pts, N_pts;
 451     grid->GetCellPoints(id_cell, N_pts, pts);
 452     vtkIdType type_cell = grid->GetCellType(id_cell);
 453     if (needs_correction[id_cell]) {
 454       for (int i_pts = 0; i_pts < N_pts; ++i_pts) {
 455         pts[i_pts] = old2new[pts[i_pts]];
 456       }
 457     }
 458     vtkIdType id_new_cell = new_grid->InsertNextCell(type_cell, N_pts, pts);
 459     copyCellData(grid, id_cell, new_grid, id_new_cell);
 460   }
 461
 462 //   writeGrid(new_grid, "pre-createBoundaryElements");
 463
 464   createBoundaryElements(new_grid);
 465   UpdateCellIndex(new_grid);
 466   grid->DeepCopy(new_grid);
 467   cout << "done." << endl;
 468 }
 469