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