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         setIntersection(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   for (int i_faces = 0; i_faces < faces.size(); ++i_faces) {
 168     for (int j_faces = 0; j_faces < faces[i_faces].size(); ++j_faces) {
 169       vtkIdType p1 = faces[i_faces][j_faces];
 170       vtkIdType p2 = faces[i_faces][0];
 171       if (j_faces < faces[i_faces].size() - 1) {
 172         p2 = faces[i_faces][j_faces+1];
 173       };
 174       bool consider_edge = false;
 175       if ((node_status->GetValue(p1) & 2) && (node_status->GetValue(p2) & 2)) {
 176         consider_edge = true;
 177       };
 178       if ((node_layer->GetValue(p1) == 0) && (node_layer->GetValue(p2) == 0)) {
 179         consider_edge = true;
 180       };
 181       if (consider_edge) {
 182         QSet<int> faces_on_edge;
 183         setIntersection(n2f[p1], n2f[p2], faces_on_edge);
 184         if (faces_on_edge.size() == 0) {
 185           EG_BUG;
 186         } else if (faces_on_edge.size() == 1) {
 187           vtkIdType pts[4] = { p1, p2, old2new[p2], old2new[p1] };
 188           vtkIdType id_new_bcell = new_grid->InsertNextCell(VTK_QUAD ,4, pts);
 189           int bc = 9999;
 190           QSet<int> bc1, bc2, bc3;
 191           if (_bnodes[old2new[p1]] != -1) {
 192             foreach (int i_bcells, bn2bc[_bnodes[old2new[p1]]]) {
 193               if (!boundary_codes.contains(cell_code->GetValue(bcells[i_bcells]))) {
 194                 bc1.insert(cell_code->GetValue(bcells[i_bcells]));
 195               };
 196             };
 197           };
 198           if (_bnodes[old2new[p2]] != -1) {
 199             foreach (int i_bcells, bn2bc[_bnodes[old2new[p2]]]) {
 200               if (!boundary_codes.contains(cell_code->GetValue(bcells[i_bcells]))) {
 201                 bc2.insert(cell_code->GetValue(bcells[i_bcells]));
 202               };
 203             };
 204           };
 205           setIntersection(bc1, bc2, bc3);
 206           if (bc3.size() == 1) {
 207             bc = *(bc3.begin());
 208           } else {
 209             bc = 9999;
 210           //EG_BUG;
 211           };
 212           cell_code->SetValue(id_new_bcell, bc);
 213           node_status->SetValue(p1, node_status->GetValue(p1) | 2);
 214           node_status->SetValue(old2new[p1], node_status->GetValue(p1) | 2);
 215           node_status->SetValue(p2, node_status->GetValue(p2) | 2);
 216           node_status->SetValue(old2new[p2], node_status->GetValue(p2) | 2);
 217         } else if (faces_on_edge.size() > 2) {
 218           EG_BUG;
 219         };
 220       };
 221     };
 222   };
 223 };
 224
 225 void SeedSimplePrismaticLayer::operate()
 226 {
 227   cout << "seeding prismatic layer" << endl;
 228   prepareLayer();
 229   EG_VTKSP(vtkUnstructuredGrid, new_grid);
 230   allocateGrid(new_grid, grid->GetNumberOfCells() + N_new_cells, grid->GetNumberOfPoints() + N_new_points);
 231   vtkIdType id_new_node = 0;
 232   EG_VTKDCN(vtkIntArray, node_status_old, grid,     "node_status");
 233   EG_VTKDCN(vtkIntArray, node_status_new, new_grid, "node_status");
 234   EG_VTKDCN(vtkIntArray, node_layer_old,  grid,     "node_layer");
 235   EG_VTKDCN(vtkIntArray, node_layer_new,  new_grid, "node_layer");
 236   EG_VTKDCC(vtkIntArray, bc,              grid,     "cell_code");
 237
 238   QVector<QSet<int> > n2f(grid->GetNumberOfPoints());
 239   for (int i_faces = 0; i_faces < faces.size(); ++i_faces) {
 240     for (int j_faces = 0; j_faces < faces[i_faces].size(); ++j_faces) {
 241       n2f[faces[i_faces][j_faces]].insert(i_faces);
 242     };
 243   };
 244
 245   // copy old grid nodes to the new grid
 246   old2new.resize(grid->GetNumberOfPoints());
 247   for (vtkIdType id_node = 0; id_node < grid->GetNumberOfPoints(); ++id_node) {
 248     vec3_t x;
 249     grid->GetPoints()->GetPoint(id_node, x.data());
 250     new_grid->GetPoints()->SetPoint(id_new_node, x.data());
 251     copyNodeData(grid, id_node, new_grid, id_new_node);
 252     old2new[id_node] = id_new_node;
 253     ++id_new_node;
 254   };
 255
 256   QSet<vtkIdType> split_nodes;
 257   for (int i_faces = 0; i_faces < faces.size(); ++i_faces) {
 258     for (int j_faces = 0; j_faces < faces[i_faces].size(); ++j_faces) {
 259       split_nodes.insert(faces[i_faces][j_faces]);
 260     };
 261   };
 262   qDebug() << split_nodes.size();
 263   QVector<bool> is_split_node(grid->GetNumberOfPoints(), false);
 264   foreach (vtkIdType id_node, split_nodes) {
 265     vec3_t x;
 266     grid->GetPoints()->GetPoint(id_node, x.data());
 267
 268
 269     if (n2f[id_node].size() > 0) {
 270       vec3_t n(0,0,0);
 271       double L = 1e99;
 272       foreach (int i_faces, n2f[id_node]) {
 273         vec3_t a,b;
 274         if (faces[i_faces][0] == id_node) {
 275           grid->GetPoint(faces[i_faces][1],a.data());
 276           grid->GetPoint(faces[i_faces][2],b.data());
 277         };
 278         if (faces[i_faces][1] == id_node) {
 279           grid->GetPoint(faces[i_faces][2],a.data());
 280           grid->GetPoint(faces[i_faces][0],b.data());
 281         };
 282         if (faces[i_faces][2] == id_node) {
 283           grid->GetPoint(faces[i_faces][0],a.data());
 284           grid->GetPoint(faces[i_faces][1],b.data());
 285         };
 286         a -= x;
 287         b -= x;
 288         L = min(a.abs(),L);
 289         L = min(b.abs(),L);
 290         a.normalise();
 291         b.normalise();
 292         vec3_t nf = a.cross(b);
 293         nf.normalise();
 294         double alpha = acos(a*b);
 295         n += alpha*nf;
 296       };
 297       for (int i_boundary_correction = 0; i_boundary_correction < 20; ++i_boundary_correction) {
 298         foreach (vtkIdType id_cell, n2c[_nodes[id_node]]) {
 299           if (isSurface(id_cell, grid)) {
 300             if (!boundary_codes.contains(bc->GetValue(id_cell))) {
 301               double A = GeometryTools::cellVA(grid, id_cell);
 302               if (A > 1e-60) {
 303                 vec3_t nf = GeometryTools::cellNormal(grid, id_cell);
 304                 nf.normalise();
 305                 n -= (nf*n)*nf;
 306               };
 307             };
 308           };
 309         };
 310       };
 311       n.normalise();
 312       x += 0.01*L*n;
 313     } else {
 314       EG_BUG;
 315     };
 316
 317     new_grid->GetPoints()->SetPoint(id_new_node, x.data());
 318     old2new[id_node] = id_new_node;
 319     node_status_new->SetValue(id_new_node, node_status_old->GetValue(id_node));
 320     node_layer_new->SetValue(id_new_node, node_layer_old->GetValue(id_node) + 1);
 321     ++id_new_node;
 322     is_split_node[id_node] = true;
 323   };
 324   QVector<bool> needs_correction(grid->GetNumberOfCells(), false);
 325   for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
 326     vtkIdType type_cell = grid->GetCellType(id_cell);
 327     if ((type_cell == VTK_TRIANGLE) || (type_cell == VTK_TETRA)) {
 328       vtkIdType *pts, N_pts;
 329       grid->GetCellPoints(id_cell, N_pts, pts);
 330       bool split = false;
 331       for (int i_pts = 0; i_pts < N_pts; ++i_pts) {
 332         if (is_split_node[pts[i_pts]]) {
 333           split = true;
 334           if (node_layer_old->GetValue(pts[i_pts]) != old_layer) {
 335             EG_BUG;
 336           };
 337         };
 338       };
 339       if (split) {
 340         if (type_cell == VTK_TETRA) {
 341           needs_correction[id_cell] = true;
 342         } else {
 343           bool f = false;
 344           if (old_layer > 0) {
 345             for (int i_pts = 0; i_pts < N_pts; ++i_pts) {
 346               if (node_layer_old->GetValue(pts[i_pts]) == old_layer - 1) {
 347                 f = true;
 348               };
 349             };
 350           };
 351           if (!f) {
 352             needs_correction[id_cell] = true;
 353           } else {
 354             cout << "dodgy face: " << id_cell << "   ";
 355             for (int i_pts = 0; i_pts < N_pts; ++i_pts) {
 356               cout << "(" << pts[i_pts] << "," << node_layer_old->GetValue(pts[i_pts]) << ") ";
 357             };
 358             cout << endl;
 359             for (int i_pts = 0; i_pts < N_pts; ++i_pts) {
 360               if (is_split_node[pts[i_pts]]) {
 361                 vec3_t x;
 362                 grid->GetPoint(pts[i_pts], x.data());
 363                 cout << "split node: " << pts[i_pts] << "  " << x << endl;
 364               };
 365             };
 366           };
 367         };
 368       };
 369     };
 370   };
 371   foreach (vtkIdType id_cell, layer_cells) {
 372     needs_correction[id_cell] = false;
 373   };
 374
 375   QVector<bool> is_in_vol_cells(grid->GetNumberOfCells(), false);
 376   for (int i_faces = 0; i_faces < faces.size(); ++i_faces) {
 377     vtkIdType id_vol_cell = vol_cells[i_faces];
 378     if (id_vol_cell >= 0) {
 379       is_in_vol_cells[id_vol_cell] = true;
 380       vtkIdType type_vol_cell = grid->GetCellType(id_vol_cell);
 381       if (type_vol_cell == VTK_TETRA) {
 382         vtkIdType *pts, N_pts, p[6];
 383         grid->GetCellPoints(id_vol_cell, N_pts, pts);
 384         p[0] = faces[i_faces][0];
 385         p[1] = faces[i_faces][1];
 386         p[2] = faces[i_faces][2];
 387         p[3] = old2new[p[0]];
 388         p[4] = old2new[p[1]];
 389         p[5] = old2new[p[2]];
 390         {
 391           vtkIdType pts[6] = { p[2], p[1], p[0], p[5], p[4], p[3] };
 392           layer_cells[i_faces] = new_grid->InsertNextCell(VTK_WEDGE, 6, pts);
 393         };
 394       } else {
 395         qDebug() << type_vol_cell;
 396         EG_BUG;
 397       };
 398     };
 399   };
 400
 401   // copy old cells to the new grid
 402   for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
 403     vtkIdType *pts, N_pts;
 404     grid->GetCellPoints(id_cell, N_pts, pts);
 405     vtkIdType type_cell = grid->GetCellType(id_cell);
 406     if (needs_correction[id_cell]) {
 407       for (int i_pts = 0; i_pts < N_pts; ++i_pts) {
 408         pts[i_pts] = old2new[pts[i_pts]];
 409       };
 410     };
 411     vtkIdType id_new_cell = new_grid->InsertNextCell(type_cell, N_pts, pts);
 412     copyCellData(grid, id_cell, new_grid, id_new_cell);
 413   };
 414
 415   createBoundaryElements(new_grid);
 416   UpdateCellIndex(new_grid);
 417   grid->DeepCopy(new_grid);
 418   cout << "done." << endl;
 419 };
 420