src/guidivideboundarylayer.cpp

   1 //
   2 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   3 // +                                                                      +
   4 // + This file is part of enGrid.                                         +
   5 // +                                                                      +
   6 // + Copyright 2008-2010 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 "guidivideboundarylayer.h"
  24 #include "math/linsolve.h"
  25
  26 #include "volumedefinition.h"
  27 #include "guimainwindow.h"
  28
  29 void GuiDivideBoundaryLayer::before()
  30 {
  31   populateBoundaryCodes(ui.listWidgetBC);
  32   populateVolumes(ui.listWidgetVC);
  33
  34   m_rest_grid = vtkUnstructuredGrid::New();
  35 }
  36
  37 bool GuiDivideBoundaryLayer::findBoundaryLayer()
  38 {
  39   l2g_t cells = getPartCells();
  40   l2l_t c2c   = getPartC2C();
  41
  42   pairs.clear();
  43   insert_cell.fill(true,m_Grid->GetNumberOfCells());
  44   N_prisms = 0;
  45   N_quads = 0;
  46   for (int i_cells = 0; i_cells < cells.size(); ++i_cells) {
  47     if (m_Grid->GetCellType(cells[i_cells]) == VTK_WEDGE) {
  48       ++N_prisms;
  49       vtkIdType N_pts, *pts;
  50       m_Grid->GetCellPoints(cells[i_cells],N_pts,pts);
  51       for (int j = 0; j < 3; ++j) {
  52         pairs.insert(QPair<vtkIdType,vtkIdType>(pts[j],pts[j+3]));
  53       }
  54       for (int j = 2; j < 5; ++j) {
  55         if (c2c[i_cells][j] != -1) {
  56           vtkIdType type_ncell = m_Grid->GetCellType(cells[c2c[i_cells][j]]);
  57           if ((type_ncell != VTK_WEDGE) && (type_ncell != VTK_QUAD)) {
  58             after();
  59             EG_ERR_RETURN("unable to identify boundary layer");
  60             return(false);
  61           }
  62         } else {
  63           vec3_t x;
  64           m_Grid->GetPoint(pts[0],x.data());
  65           cout << x << endl;
  66           EG_BUG;
  67         }
  68       }
  69       for (int j = 0; j < 2; ++j) {
  70         if (c2c[i_cells][j] != -1) {
  71           vtkIdType type_ncell = m_Grid->GetCellType(cells[c2c[i_cells][j]]);
  72           if (type_ncell == VTK_WEDGE) {
  73             after();
  74             EG_ERR_RETURN("the boundary layer seems to have been split already");
  75             return(false);
  76           }
  77         } else {
  78           EG_BUG;
  79         }
  80       }
  81     }
  82     if (m_Grid->GetCellType(cells[i_cells]) == VTK_QUAD) {
  83       ++N_quads;
  84     }
  85   }
  86   if (N_prisms == 0) {
  87     after();
  88     EG_ERR_RETURN("unable to identify boundary layer");
  89     return(false);
  90   }
  91   is_blayer_node.clear();
  92   is_blayer_node.fill(false, m_Grid->GetNumberOfPoints());
  93
  94   QPair<vtkIdType,vtkIdType> P;
  95   foreach (P, pairs) {
  96     is_blayer_node[P.second] = true;
  97   }
  98
  99   return(true);
 100 }
 101
 102 void GuiDivideBoundaryLayer::findBoundaryLayer1()
 103 {
 104   pairs.clear();
 105
 106   l2g_t  cells = getPartCells();
 107   g2l_t _nodes = getPartLocalNodes();
 108   g2l_t _cells = getPartLocalCells();
 109   l2l_t  n2c   = getPartN2C();
 110   l2l_t  c2c   = getPartC2C();
 111
 112   QSet<int> bcs;
 113   getSelectedItems(ui.listWidgetBC, bcs);
 114   QVector<vtkIdType> scells;
 115   getSurfaceCells(bcs, scells, m_Grid);
 116
 117   N_prisms = 0;
 118   N_quads = 0;
 119   N_hexes = 0;
 120
 121   insert_cell.fill(true,m_Grid->GetNumberOfCells());
 122
 123   foreach (int i_scells, scells) {
 124
 125     vtkIdType id_scell   = scells[i_scells];
 126     int       i_cells    = findVolumeCell(m_Grid,id_scell,_nodes,cells,_cells,n2c);
 127     vtkIdType id_cell    = cells[i_cells];
 128     vtkIdType type_cell  = m_Grid->GetCellType(id_cell);
 129     vtkIdType type_scell = m_Grid->GetCellType(id_scell);
 130
 131     vtkIdType N_pts, *pts;
 132     m_Grid->GetCellPoints(id_cell,N_pts,pts);
 133
 134     insert_cell[id_cell] = false;
 135
 136     if (type_cell == VTK_WEDGE) {
 137
 138       // prisms
 139       //
 140       ++N_prisms;
 141       if (type_scell != VTK_QUAD) {
 142         EG_ERR_RETURN("unable to identify boundary layer");
 143       }
 144       if (cells[c2c[i_cells][0]] == id_scell) {
 145         for (int j = 0; j < 3; ++j) {
 146           pairs.insert(QPair<vtkIdType,vtkIdType>(pts[j],pts[j+3]));
 147         }
 148       } else if (cells[c2c[i_cells][1]] == id_scell) {
 149         for (int j = 0; j < 3; ++j) {
 150           pairs.insert(QPair<vtkIdType,vtkIdType>(pts[j+3],pts[j]));
 151         }
 152       } else {
 153         EG_ERR_RETURN("unable to identify boundary layer");
 154       }
 155       for (int j = 2; j < 5; ++j) {
 156         if (c2c[i_cells][j] != -1) {
 157           vtkIdType id_ncell   = cells[c2c[i_cells][j]];
 158           vtkIdType type_ncell = m_Grid->GetCellType(id_ncell);
 159           if ((type_ncell != VTK_WEDGE) && (type_ncell != VTK_HEXAHEDRON) && (type_ncell != VTK_QUAD)) {
 160             EG_ERR_RETURN("unable to identify boundary layer");
 161           }
 162           if (type_ncell == VTK_QUAD) {
 163             insert_cell[id_ncell] = false;
 164           }
 165         } else {
 166           EG_BUG;
 167         }
 168       }
 169     } else if (type_cell == VTK_HEXAHEDRON) {
 170
 171       // hexes
 172       //
 173       ++N_hexes;
 174       if (type_scell != VTK_QUAD) {
 175         EG_ERR_RETURN("unable to identify boundary layer");
 176       }
 177       QVector<bool> check_neigh(6,true);
 178       if (cells[c2c[i_cells][0]] == id_scell) {
 179         pairs.insert(QPair<vtkIdType,vtkIdType>(pts[0],pts[4]));
 180         pairs.insert(QPair<vtkIdType,vtkIdType>(pts[1],pts[5]));
 181         pairs.insert(QPair<vtkIdType,vtkIdType>(pts[2],pts[6]));
 182         pairs.insert(QPair<vtkIdType,vtkIdType>(pts[3],pts[7]));
 183         check_neigh[0] = false;
 184         check_neigh[1] = false;
 185       } else if (cells[c2c[i_cells][1]] == id_scell) {
 186         pairs.insert(QPair<vtkIdType,vtkIdType>(pts[4],pts[0]));
 187         pairs.insert(QPair<vtkIdType,vtkIdType>(pts[5],pts[1]));
 188         pairs.insert(QPair<vtkIdType,vtkIdType>(pts[6],pts[2]));
 189         pairs.insert(QPair<vtkIdType,vtkIdType>(pts[7],pts[3]));
 190         check_neigh[0] = false;
 191         check_neigh[1] = false;
 192       } else if (cells[c2c[i_cells][2]] == id_scell) {
 193         pairs.insert(QPair<vtkIdType,vtkIdType>(pts[0],pts[3]));
 194         pairs.insert(QPair<vtkIdType,vtkIdType>(pts[1],pts[2]));
 195         pairs.insert(QPair<vtkIdType,vtkIdType>(pts[5],pts[6]));
 196         pairs.insert(QPair<vtkIdType,vtkIdType>(pts[4],pts[7]));
 197         check_neigh[2] = false;
 198         check_neigh[3] = false;
 199       } else if (cells[c2c[i_cells][3]] == id_scell) {
 200         pairs.insert(QPair<vtkIdType,vtkIdType>(pts[3],pts[0]));
 201         pairs.insert(QPair<vtkIdType,vtkIdType>(pts[2],pts[1]));
 202         pairs.insert(QPair<vtkIdType,vtkIdType>(pts[6],pts[5]));
 203         pairs.insert(QPair<vtkIdType,vtkIdType>(pts[7],pts[4]));
 204         check_neigh[2] = false;
 205         check_neigh[3] = false;
 206       } else if (cells[c2c[i_cells][4]] == id_scell) {
 207         pairs.insert(QPair<vtkIdType,vtkIdType>(pts[0],pts[1]));
 208         pairs.insert(QPair<vtkIdType,vtkIdType>(pts[4],pts[5]));
 209         pairs.insert(QPair<vtkIdType,vtkIdType>(pts[7],pts[6]));
 210         pairs.insert(QPair<vtkIdType,vtkIdType>(pts[3],pts[2]));
 211         check_neigh[4] = false;
 212         check_neigh[5] = false;
 213       } else if (cells[c2c[i_cells][5]] == id_scell) {
 214         pairs.insert(QPair<vtkIdType,vtkIdType>(pts[1],pts[0]));
 215         pairs.insert(QPair<vtkIdType,vtkIdType>(pts[5],pts[4]));
 216         pairs.insert(QPair<vtkIdType,vtkIdType>(pts[6],pts[7]));
 217         pairs.insert(QPair<vtkIdType,vtkIdType>(pts[2],pts[3]));
 218         check_neigh[4] = false;
 219         check_neigh[5] = false;
 220       }
 221       for (int j = 0; j < 6; ++j) {
 222         if (check_neigh[j]) {
 223           if (c2c[i_cells][j] != -1) {
 224             vtkIdType id_ncell   = cells[c2c[i_cells][j]];
 225             vtkIdType type_ncell = m_Grid->GetCellType(id_ncell);
 226             if ((type_ncell != VTK_WEDGE) && (type_ncell != VTK_HEXAHEDRON) && (type_ncell != VTK_QUAD)) {
 227               EG_ERR_RETURN("unable to identify boundary layer");
 228             }
 229             if (type_ncell == VTK_QUAD) {
 230               insert_cell[id_ncell] = false;
 231             }
 232           } else {
 233             EG_BUG;
 234           }
 235         }
 236       }
 237
 238     } else if (type_cell == VTK_QUAD) {
 239
 240       // quads
 241       //
 242       ++N_quads;
 243
 244     } else {
 245       EG_ERR_RETURN("unable to identify boundary layer");
 246     }
 247   }
 248
 249   /*
 250   N_prisms = 0;
 251   N_quads = 0;
 252   N_hexes = 0;
 253   for (int i_cells = 0; i_cells < cells.size(); ++i_cells) {
 254     if (m_Grid->GetCellType(cells[i_cells]) == VTK_WEDGE) {
 255       ++N_prisms;
 256       vtkIdType N_pts, *pts;
 257       m_Grid->GetCellPoints(cells[i_cells],N_pts,pts);
 258       for (int j = 0; j < 3; ++j) {
 259         pairs.insert(QPair<vtkIdType,vtkIdType>(pts[j],pts[j+3]));
 260       }
 261       for (int j = 2; j < 5; ++j) {
 262         if (c2c[i_cells][j] != -1) {
 263           vtkIdType type_ncell = m_Grid->GetCellType(cells[c2c[i_cells][j]]);
 264           if ((type_ncell != VTK_WEDGE) && (type_ncell != VTK_QUAD)) {
 265             EG_ERR_RETURN("unable to identify boundary layer");
 266           }
 267         } else {
 268           EG_BUG;
 269         }
 270       }
 271       for (int j = 0; j < 2; ++j) {
 272         if (c2c[i_cells][j] != -1) {
 273           vtkIdType type_ncell = m_Grid->GetCellType(cells[c2c[i_cells][j]]);
 274           if (type_ncell == VTK_WEDGE) {
 275             EG_ERR_RETURN("the boundary layer seems to have been split already");
 276           }
 277         } else {
 278           EG_BUG;
 279         }
 280       }
 281     }
 282     if (m_Grid->GetCellType(cells[i_cells]) == VTK_QUAD) {
 283       ++N_quads;
 284     }
 285   }
 286   */
 287
 288   if (N_prisms + N_hexes == 0) {
 289     EG_ERR_RETURN("unable to identify boundary layer");
 290   }
 291
 292   is_blayer_node.clear();
 293   is_blayer_node.fill(false, m_Grid->GetNumberOfPoints());
 294
 295   QPair<vtkIdType,vtkIdType> P;
 296   foreach (P, pairs) {
 297     is_blayer_node[P.second] = true;
 298   }
 299 }
 300
 301 void GuiDivideBoundaryLayer::bisectF(double &f1, double &f2)
 302 {
 303   f = 0.5*(f1+f2);
 304   computeY();
 305   if (y[y.size()-1] > 1) f2 = f;
 306   else                   f1 = f;
 307   f = 0.5*(f1+f2);
 308 }
 309
 310
 311 //begin GROSSER MURKS
 312
 313 void GuiDivideBoundaryLayer::computeF()
 314 {
 315   double f1 = 0;
 316   double f2 = 100;//1.0/h;
 317   double err = 0;
 318   y[1] = min(0.33,y[1]);
 319   do {
 320     //cout << f1 << ',' << f2 << endl;
 321     bisectF(f1,f2);
 322     //while (fabs(1-y[y.size()-1]) > 1e-6);
 323     err = fabs((y[y.size()-1]-1)/(y[y.size()-1]-y[y.size()-2]));
 324   } while (err > 0.01);
 325   f = 0.5*(f1+f2);
 326 }
 327
 328 void GuiDivideBoundaryLayer::computeY()
 329 {
 330   double C = F;
 331   for (int i = 2; i < y.size(); ++i) {
 332     y[i] = y[i-1] + C*(y[i-1]-y[i-2]);
 333     C *= f;
 334   }
 335 }
 336
 337 //end GROSSER MURKS
 338
 339 void GuiDivideBoundaryLayer::createEdges(vtkUnstructuredGrid *new_grid)
 340 {
 341   edges.fill(QVector<vtkIdType>(N_layers+1), pairs.size());
 342   old2edge.fill(-1, m_Grid->GetNumberOfPoints());
 343   int N = 0;
 344   vtkIdType id_new_node = m_Grid->GetNumberOfPoints();
 345   QPair<vtkIdType,vtkIdType> P;
 346   double max_step = 0;
 347   double ymax     = 0;
 348   double ymin     = 1e99;
 349   foreach (P, pairs) {
 350     edges[N][0] = P.first;
 351     edges[N][N_layers] = P.second;
 352     old2edge[P.first] = N;
 353     old2edge[P.second] = N;
 354
 355     vec3_t x1,x2;
 356     m_Grid->GetPoint(P.first, x1.data());
 357     m_Grid->GetPoint(P.second, x2.data());
 358     vec3_t n = x2-x1;
 359     if (!ui.checkBoxH->isChecked() || !y_computed) {
 360       y.resize(N_layers + 1);
 361       x.resize(y.size());
 362       for (int i = 0; i < x.size(); ++i) {
 363         x[i] = i*1.0/(x.size() - 1);
 364       }
 365       double Dy1;
 366       if (ui.checkBoxH->isChecked()) {
 367         Dy1 = h;
 368       } else {
 369         Dy1 = h/n.abs();
 370       }
 371       y[0]   = 0;
 372       y[1]   = Dy1;
 373       computeF();
 374       computeY();
 375       y_computed = true;
 376     }
 377     ymin = min(ymin,y[1]*n.abs());
 378     ymax = max(ymax,y[1]*n.abs());
 379     for (int i = 1; i < N_layers; ++i) {
 380       vec3_t x = x1 + y[i]*n;
 381       max_step = max(max_step,(y[i+1]-y[i])/(y[i]-y[i-1]));
 382       new_grid->GetPoints()->SetPoint(id_new_node, x.data());
 383       edges[N][i] = id_new_node;
 384       ++id_new_node;
 385     }
 386     ++N;
 387   }
 388   cout << LINE;
 389   cout << "maximal increment : " << max_step << endl;
 390   cout << "min(y) : " << ymin << endl;
 391   cout << "max(y) : " << ymax << endl;
 392   cout << LINE;
 393 }
 394
 395 void GuiDivideBoundaryLayer::operate1()
 396 {
 397   y_computed = false;
 398   N_layers = ui.spinBoxLayers->value();
 399   h = ui.lineEditH->text().toDouble();
 400   F = ui.doubleSpinBoxF->value();
 401   cout << "dividing boundary layer into " << N_layers << " layers:" << endl;
 402   findBoundaryLayer1();
 403
 404   EG_VTKSP(vtkUnstructuredGrid,new_grid);
 405   int N_new_cells = m_Grid->GetNumberOfCells() + (N_prisms + N_hexes + N_quads)*(N_layers-1);
 406   int N_new_nodes = m_Grid->GetNumberOfPoints() + pairs.size()*(N_layers-1);
 407   allocateGrid(new_grid, N_new_cells, N_new_nodes);
 408
 409   // copy existing mesh without prisms and adjacent cells
 410   vtkIdType id_new_node = 0;
 411   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 412     vec3_t x;
 413     m_Grid->GetPoint(id_node, x.data());
 414     new_grid->GetPoints()->SetPoint(id_new_node, x.data());
 415     copyNodeData(m_Grid, id_node, new_grid, id_new_node);
 416     ++id_new_node;
 417   }
 418   vtkIdType id_new_cell;
 419   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 420     vtkIdType N_pts, *pts;
 421     m_Grid->GetCellPoints(id_cell, N_pts, pts);
 422     //bool insert_cell = true;
 423     //if (m_Grid->GetCellType(id_cell) == VTK_WEDGE) insert_cell = false;
 424     //if (m_Grid->GetCellType(id_cell) == VTK_QUAD) insert_cell = false;
 425     if (insert_cell[id_cell]) {
 426       id_new_cell = new_grid->InsertNextCell(m_Grid->GetCellType(id_cell), N_pts, pts);
 427       copyCellData(m_Grid, id_cell, new_grid, id_new_cell);
 428     }
 429   }
 430
 431   // create divided boundary layer
 432   createEdges(new_grid);
 433
 434   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 435     if (!insert_cell[id_cell]) {
 436       if (m_Grid->GetCellType(id_cell) == VTK_WEDGE) {
 437         vtkIdType N_pts, *pts;
 438         m_Grid->GetCellPoints(id_cell, N_pts, pts);
 439         for (int i = 0; i < N_layers; ++i) {
 440           vtkIdType p[6];
 441           p[0] = edges[old2edge[pts[0]]][i];
 442           p[1] = edges[old2edge[pts[1]]][i];
 443           p[2] = edges[old2edge[pts[2]]][i];
 444           p[3] = edges[old2edge[pts[0]]][i+1];
 445           p[4] = edges[old2edge[pts[1]]][i+1];
 446           p[5] = edges[old2edge[pts[2]]][i+1];
 447           id_new_cell = new_grid->InsertNextCell(VTK_WEDGE, 6, p);
 448           copyCellData(m_Grid, id_cell, new_grid, id_new_cell);
 449         }
 450       }
 451       if (m_Grid->GetCellType(id_cell) == VTK_HEXAHEDRON) {
 452         EG_BUG;
 453         // MURKS!!!!
 454         vtkIdType N_pts, *pts;
 455         m_Grid->GetCellPoints(id_cell, N_pts, pts);
 456         for (int i = 0; i < N_layers; ++i) {
 457           vtkIdType p[6];
 458           p[0] = edges[old2edge[pts[0]]][i];
 459           p[1] = edges[old2edge[pts[1]]][i];
 460           p[2] = edges[old2edge[pts[2]]][i];
 461           p[3] = edges[old2edge[pts[0]]][i+1];
 462           p[4] = edges[old2edge[pts[1]]][i+1];
 463           p[5] = edges[old2edge[pts[2]]][i+1];
 464           id_new_cell = new_grid->InsertNextCell(VTK_WEDGE, 6, p);
 465           copyCellData(m_Grid, id_cell, new_grid, id_new_cell);
 466         }
 467
 468
 469       }
 470       if (m_Grid->GetCellType(id_cell) == VTK_QUAD) {
 471         vtkIdType N_pts, *pts;
 472         m_Grid->GetCellPoints(id_cell, N_pts, pts);
 473         if ((old2edge[pts[0]] != -1) && (old2edge[pts[1]] != -1) && (old2edge[pts[2]] != -1) && (old2edge[pts[3]] != -1)) {
 474           for (int i = 0; i < N_layers; ++i) {
 475             vtkIdType p[4];
 476             p[0] = edges[old2edge[pts[0]]][i];
 477             p[1] = edges[old2edge[pts[1]]][i];
 478             p[2] = edges[old2edge[pts[1]]][i+1];
 479             p[3] = edges[old2edge[pts[0]]][i+1];
 480             id_new_cell = new_grid->InsertNextCell(VTK_QUAD, 4, p);
 481             copyCellData(m_Grid, id_cell, new_grid, id_new_cell);
 482           }
 483         }
 484       }
 485     }
 486   }
 487
 488   makeCopy(new_grid, m_Grid);
 489 }
 490
 491 void GuiDivideBoundaryLayer::operate()
 492 {
 493   ///////////////////////////////////////////////////////////////
 494   // set m_Grid to selected volume
 495   getSelectedItems(ui.listWidgetBC, m_BoundaryCodes); // fill m_BoundaryCodes with values from listWidgetBC
 496   QString volume_name = getSelectedVolume(ui.listWidgetVC);
 497   VolumeDefinition V = GuiMainWindow::pointer()->getVol(volume_name);
 498   foreach (int bc, m_BoundaryCodes) {
 499     qDebug()<<"V.getSign("<<bc<<")="<<V.getSign(bc);
 500     if (V.getSign(bc) == 0) {
 501       QString msg;
 502       msg.setNum(bc);
 503       msg = "Boundary code " + msg + " is not part of the volume '" + volume_name +"'.";
 504       EG_ERR_RETURN(msg);
 505     }
 506   }
 507
 508 //   EG_VTKSP(vtkUnstructuredGrid, m_rest_grid);
 509   {
 510     EG_VTKSP(vtkUnstructuredGrid, vol_grid);
 511     MeshPartition volume(volume_name);
 512     MeshPartition rest(m_Grid);
 513     rest.setRemainder(volume);
 514     volume.setVolumeOrientation();
 515     volume.extractToVtkGrid(vol_grid);
 516     rest.extractToVtkGrid(m_rest_grid);
 517     makeCopy(vol_grid, m_Grid);
 518   }
 519   setAllCells();
 520
 521   //   writeGrid(m_Grid,"selected_volume");
 522   //   return;
 523   ///////////////////////////////////////////////////////////////
 524
 525   y_computed = false;
 526   N_layers = ui.spinBoxLayers->value();
 527   h = ui.lineEditH->text().toDouble();
 528   F = ui.doubleSpinBoxF->value();
 529   cout << "dividing boundary layer into " << N_layers << " layers:" << endl;
 530   if(findBoundaryLayer()) {
 531     EG_VTKSP(vtkUnstructuredGrid,new_grid);
 532     allocateGrid(new_grid, m_Grid->GetNumberOfCells() + (N_prisms + N_quads)*(N_layers-1), m_Grid->GetNumberOfPoints() + pairs.size()*(N_layers-1));
 533
 534
 535     EG_VTKDCC(vtkIntArray, old_orgdir, m_Grid, "cell_orgdir");
 536     EG_VTKDCC(vtkIntArray, old_voldir, m_Grid, "cell_voldir");
 537     EG_VTKDCC(vtkIntArray, old_curdir, m_Grid, "cell_curdir");
 538     EG_VTKDCC(vtkIntArray, new_orgdir, new_grid, "cell_orgdir");
 539     EG_VTKDCC(vtkIntArray, new_voldir, new_grid, "cell_voldir");
 540     EG_VTKDCC(vtkIntArray, new_curdir, new_grid, "cell_curdir");
 541
 542     int orgdir = -99;
 543     int curdir = -99;
 544     int voldir = -99;
 545
 546     // copy existing mesh without prisms and adjacent cells
 547     vtkIdType id_new_node = 0;
 548     for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 549       vec3_t x;
 550       m_Grid->GetPoint(id_node, x.data());
 551       new_grid->GetPoints()->SetPoint(id_new_node, x.data());
 552       copyNodeData(m_Grid, id_node, new_grid, id_new_node);
 553       ++id_new_node;
 554     }
 555     vtkIdType id_new_cell;
 556     for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 557       vtkIdType N_pts, *pts;
 558       m_Grid->GetCellPoints(id_cell, N_pts, pts);
 559       bool insert_cell = true;
 560       if (m_Grid->GetCellType(id_cell) == VTK_WEDGE) {
 561         insert_cell = false;
 562       }
 563       if (m_Grid->GetCellType(id_cell) == VTK_QUAD) {
 564         insert_cell = false;
 565         if (orgdir != -99 && old_orgdir->GetValue(id_cell) != orgdir) {
 566           EG_BUG;
 567         }
 568         if (voldir != -99 && old_voldir->GetValue(id_cell) != voldir) {
 569           EG_BUG;
 570         }
 571         if (curdir != -99 && old_curdir->GetValue(id_cell) != curdir) {
 572           EG_BUG;
 573         }
 574         orgdir = old_orgdir->GetValue(id_cell);
 575         voldir = old_voldir->GetValue(id_cell);
 576         curdir = old_curdir->GetValue(id_cell);
 577       }
 578       if (insert_cell) {
 579         id_new_cell = new_grid->InsertNextCell(m_Grid->GetCellType(id_cell), N_pts, pts);
 580         copyCellData(m_Grid, id_cell, new_grid, id_new_cell);
 581       }
 582     }
 583
 584   // create divided boundary layer
 585     createEdges(new_grid);
 586
 587     for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 588       if (m_Grid->GetCellType(id_cell) == VTK_WEDGE) {
 589         vtkIdType N_pts, *pts;
 590         m_Grid->GetCellPoints(id_cell, N_pts, pts);
 591         for (int i = 0; i < N_layers; ++i) {
 592           vtkIdType p[6];
 593           p[0] = edges[old2edge[pts[0]]][i];
 594           p[1] = edges[old2edge[pts[1]]][i];
 595           p[2] = edges[old2edge[pts[2]]][i];
 596           p[3] = edges[old2edge[pts[0]]][i+1];
 597           p[4] = edges[old2edge[pts[1]]][i+1];
 598           p[5] = edges[old2edge[pts[2]]][i+1];
 599           id_new_cell = new_grid->InsertNextCell(VTK_WEDGE, 6, p);
 600           copyCellData(m_Grid, id_cell, new_grid, id_new_cell);
 601         }
 602       }
 603       if (m_Grid->GetCellType(id_cell) == VTK_QUAD) {
 604         vtkIdType N_pts, *pts;
 605         m_Grid->GetCellPoints(id_cell, N_pts, pts);
 606         if ((old2edge[pts[0]] != -1) && (old2edge[pts[1]] != -1) && (old2edge[pts[2]] != -1) && (old2edge[pts[3]] != -1)) {
 607           for (int i = 0; i < N_layers; ++i) {
 608             vtkIdType p[4];
 609             p[0] = edges[old2edge[pts[0]]][i];
 610             p[1] = edges[old2edge[pts[1]]][i];
 611             p[2] = edges[old2edge[pts[1]]][i+1];
 612             p[3] = edges[old2edge[pts[0]]][i+1];
 613             id_new_cell = new_grid->InsertNextCell(VTK_QUAD, 4, p);
 614             copyCellData(m_Grid, id_cell, new_grid, id_new_cell);
 615             new_orgdir->SetValue(id_new_cell, orgdir);
 616             new_voldir->SetValue(id_new_cell, voldir);
 617             new_curdir->SetValue(id_new_cell, curdir);
 618           }
 619         }
 620       }
 621     }
 622
 623     makeCopy(new_grid, m_Grid);
 624   }
 625
 626   ///////////////////////////////////////////////////////////////
 627   after();
 628   ///////////////////////////////////////////////////////////////
 629 }
 630
 631 void GuiDivideBoundaryLayer::after()
 632 {
 633   // set m_Grid to modified selected volume + unselected volumes
 634   {
 635     MeshPartition volume(m_Grid, true);
 636     MeshPartition rest(m_rest_grid, true);
 637     volume.addPartition(rest);
 638   }
 639   resetOrientation(m_Grid);
 640   createIndices(m_Grid);
 641   m_rest_grid->Delete();
 642 }