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