src/libengrid/seedsimpleprismaticlayer.cpp

   1 //
   2 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   3 // +                                                                      +
   4 // + This file is part of enGrid.                                         +
   5 // +                                                                      +
   6 // + Copyright 2008-2013 enGits GmbH                                      +
   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   m_Grid->BuildLinks();
  47   EG_VTKSP(vtkIdList, nds);
  48   EG_VTKSP(vtkIdList, cls);
  49   EG_VTKDCN(vtkIntArray, node_layer,  m_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 = m_Grid->GetCellType(id_cell);
  57     vtkIdType *pts;
  58     vtkIdType  Npts;
  59     m_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       m_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), m_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       m_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), m_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(m_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,  m_Grid, "node_layer");
 117   EG_VTKDCN(vtkIntArray, node_status, m_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(m_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, m_Grid->GetNumberOfCells() + N_new_cells, m_Grid->GetNumberOfPoints() + N_new_points);
 273   vtkIdType id_new_node = 0;
 274   EG_VTKDCN(vtkIntArray,    node_status_old, m_Grid,   "node_status");
 275   EG_VTKDCN(vtkIntArray,    node_status_new, new_grid, "node_status");
 276   EG_VTKDCN(vtkIntArray,    node_layer_old,  m_Grid,   "node_layer");
 277   EG_VTKDCN(vtkIntArray,    node_layer_new,  new_grid, "node_layer");
 278   EG_VTKDCC(vtkIntArray,    bc,              m_Grid,   "cell_code");
 279   EG_VTKDCN(vtkDoubleArray, cl_old,          m_Grid,   "node_meshdensity_desired");
 280   EG_VTKDCN(vtkDoubleArray, cl_new,          new_grid, "node_meshdensity_desired");
 281
 282   l2l_t  n2c   = getPartN2C();
 283   g2l_t _nodes = getPartLocalNodes();
 284
 285   QVector<QSet<int> > n2f(m_Grid->GetNumberOfPoints());
 286   for (int i_faces = 0; i_faces < faces.size(); ++i_faces) {
 287     for (int j_faces = 0; j_faces < faces[i_faces].size(); ++j_faces) {
 288       n2f[faces[i_faces][j_faces]].insert(i_faces);
 289     }
 290   }
 291
 292   // copy old grid nodes to the new grid
 293   old2new.resize(m_Grid->GetNumberOfPoints());
 294   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 295     vec3_t x;
 296     m_Grid->GetPoints()->GetPoint(id_node, x.data());
 297     new_grid->GetPoints()->SetPoint(id_new_node, x.data());
 298     copyNodeData(m_Grid, id_node, new_grid, id_new_node);
 299     old2new[id_node] = id_new_node;
 300     ++id_new_node;
 301   }
 302
 303   QSet<vtkIdType> split_nodes;
 304   for (int i_faces = 0; i_faces < faces.size(); ++i_faces) {
 305     for (int j_faces = 0; j_faces < faces[i_faces].size(); ++j_faces) {
 306       split_nodes.insert(faces[i_faces][j_faces]);
 307     }
 308   }
 309   qDebug() << split_nodes.size();
 310   QVector<bool> is_split_node(m_Grid->GetNumberOfPoints(), false);
 311   foreach (vtkIdType id_node, split_nodes) {
 312     vec3_t x;
 313     m_Grid->GetPoints()->GetPoint(id_node, x.data());
 314     double h = cl_old->GetValue(id_node);
 315
 316     if (n2f[id_node].size() > 0) {
 317       vec3_t n(0,0,0);
 318       double L = 1e99;
 319       foreach (int i_faces, n2f[id_node]) {
 320         vec3_t a,b;
 321         if (faces[i_faces][0] == id_node) {
 322           m_Grid->GetPoint(faces[i_faces][1],a.data());
 323           m_Grid->GetPoint(faces[i_faces][2],b.data());
 324         }
 325         if (faces[i_faces][1] == id_node) {
 326           m_Grid->GetPoint(faces[i_faces][2],a.data());
 327           m_Grid->GetPoint(faces[i_faces][0],b.data());
 328         }
 329         if (faces[i_faces][2] == id_node) {
 330           m_Grid->GetPoint(faces[i_faces][0],a.data());
 331           m_Grid->GetPoint(faces[i_faces][1],b.data());
 332         }
 333         a -= x;
 334         b -= x;
 335         L = min(a.abs(),L);
 336         L = min(b.abs(),L);
 337         a.normalise();
 338         b.normalise();
 339         vec3_t nf = a.cross(b);
 340         nf.normalise();
 341         double alpha = acos(a*b);
 342         n += alpha*nf;
 343       }
 344       for (int i_boundary_correction = 0; i_boundary_correction < 20; ++i_boundary_correction) {
 345         foreach (vtkIdType id_cell, n2c[_nodes[id_node]]) {
 346           if (isSurface(id_cell, m_Grid)) {
 347             if (!m_BoundaryCodes.contains(bc->GetValue(id_cell))) {
 348               double A = GeometryTools::cellVA(m_Grid, id_cell);
 349               if (A > 1e-60) {
 350                 vec3_t nf = GeometryTools::cellNormal(m_Grid, id_cell);
 351                 nf.normalise();
 352                 n -= (nf*n)*nf;
 353               }
 354             }
 355           }
 356         }
 357       }
 358       n.normalise();
 359       x += 0.01*L*n;
 360     } else {
 361       EG_BUG;
 362     }
 363
 364     new_grid->GetPoints()->SetPoint(id_new_node, x.data());
 365     cl_new->SetValue(id_new_node, h);
 366     old2new[id_node] = id_new_node;
 367     node_status_new->SetValue(id_new_node, node_status_old->GetValue(id_node));
 368     node_layer_new->SetValue(id_new_node, node_layer_old->GetValue(id_node) + 1);
 369     ++id_new_node;
 370     is_split_node[id_node] = true;
 371   }
 372   QVector<bool> needs_correction(m_Grid->GetNumberOfCells(), false);
 373   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 374     vtkIdType type_cell = m_Grid->GetCellType(id_cell);
 375     if ((type_cell == VTK_TRIANGLE) || (type_cell == VTK_TETRA)) {
 376       vtkIdType *pts, N_pts;
 377       m_Grid->GetCellPoints(id_cell, N_pts, pts);
 378       bool split = false;
 379       for (int i_pts = 0; i_pts < N_pts; ++i_pts) {
 380         if (is_split_node[pts[i_pts]]) {
 381           split = true;
 382           if (node_layer_old->GetValue(pts[i_pts]) != old_layer) {
 383             EG_BUG;
 384           }
 385         }
 386       }
 387       if (split) {
 388         if (type_cell == VTK_TETRA) {
 389           needs_correction[id_cell] = true;
 390         } else {
 391           bool f = false;
 392           if (old_layer > 0) {
 393             for (int i_pts = 0; i_pts < N_pts; ++i_pts) {
 394               if (node_layer_old->GetValue(pts[i_pts]) == old_layer - 1) {
 395                 f = true;
 396               }
 397             }
 398           }
 399           if (!f) {
 400             needs_correction[id_cell] = true;
 401           } else {
 402             cout << "dodgy face: " << id_cell << "   ";
 403             for (int i_pts = 0; i_pts < N_pts; ++i_pts) {
 404               cout << "(" << pts[i_pts] << "," << node_layer_old->GetValue(pts[i_pts]) << ") ";
 405             }
 406             cout << endl;
 407             for (int i_pts = 0; i_pts < N_pts; ++i_pts) {
 408               if (is_split_node[pts[i_pts]]) {
 409                 vec3_t x;
 410                 m_Grid->GetPoint(pts[i_pts], x.data());
 411                 cout << "split node: " << pts[i_pts] << "  " << x << endl;
 412               }
 413             }
 414           }
 415         }
 416       }
 417     }
 418   }
 419   foreach (vtkIdType id_cell, layer_cells) {
 420     needs_correction[id_cell] = false;
 421   }
 422
 423   QVector<bool> is_in_vol_cells(m_Grid->GetNumberOfCells(), false);
 424   for (int i_faces = 0; i_faces < faces.size(); ++i_faces) {
 425     vtkIdType id_vol_cell = vol_cells[i_faces];
 426     if (id_vol_cell >= 0) {
 427       is_in_vol_cells[id_vol_cell] = true;
 428       vtkIdType type_vol_cell = m_Grid->GetCellType(id_vol_cell);
 429       if (type_vol_cell == VTK_TETRA) {
 430         vtkIdType *pts, N_pts, p[6];
 431         m_Grid->GetCellPoints(id_vol_cell, N_pts, pts);
 432         p[0] = faces[i_faces][0];
 433         p[1] = faces[i_faces][1];
 434         p[2] = faces[i_faces][2];
 435         p[3] = old2new[p[0]];
 436         p[4] = old2new[p[1]];
 437         p[5] = old2new[p[2]];
 438         {
 439           vtkIdType pts[6] = { p[2], p[1], p[0], p[5], p[4], p[3] };
 440           layer_cells[i_faces] = new_grid->InsertNextCell(VTK_WEDGE, 6, pts);
 441         }
 442       } else {
 443         qDebug() << type_vol_cell;
 444         EG_BUG;
 445       }
 446     }
 447   }
 448
 449 //   writeGrid(new_grid, "pre-cellcopy");
 450
 451   // copy old cells to the new grid
 452   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 453     vtkIdType *pts, N_pts;
 454     m_Grid->GetCellPoints(id_cell, N_pts, pts);
 455     vtkIdType type_cell = m_Grid->GetCellType(id_cell);
 456     if (needs_correction[id_cell]) {
 457       for (int i_pts = 0; i_pts < N_pts; ++i_pts) {
 458         pts[i_pts] = old2new[pts[i_pts]];
 459       }
 460     }
 461     vtkIdType id_new_cell = new_grid->InsertNextCell(type_cell, N_pts, pts);
 462     copyCellData(m_Grid, id_cell, new_grid, id_new_cell);
 463   }
 464
 465 //   writeGrid(new_grid, "pre-createBoundaryElements");
 466
 467   createBoundaryElements(new_grid);
 468   UpdateCellIndex(new_grid);
 469   m_Grid->DeepCopy(new_grid);
 470   cout << "done." << endl;
 471 }
 472