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