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