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