src/libengrid/polymesh.cpp

   1 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   2 // +                                                                      +
   3 // + This file is part of enGrid.                                         +
   4 // +                                                                      +
   5 // + Copyright 2008-2014 enGits GmbH                                      +
   6 // +                                                                      +
   7 // + enGrid is free software: you can redistribute it and/or modify       +
   8 // + it under the terms of the GNU General Public License as published by +
   9 // + the Free Software Foundation, either version 3 of the License, or    +
  10 // + (at your option) any later version.                                  +
  11 // +                                                                      +
  12 // + enGrid is distributed in the hope that it will be useful,            +
  13 // + but WITHOUT ANY WARRANTY; without even the implied warranty of       +
  14 // + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the        +
  15 // + GNU General Public License for more details.                         +
  16 // +                                                                      +
  17 // + You should have received a copy of the GNU General Public License    +
  18 // + along with enGrid. If not, see <http://www.gnu.org/licenses/>.       +
  19 // +                                                                      +
  20 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  21
  22 #include "polymesh.h"
  23 #include "polymolecule.h"
  24 #include "guimainwindow.h"
  25
  26 #include <vtkXMLPolyDataWriter.h>
  27 #include <vtkTriangleFilter.h>
  28 #include <vtkCurvatures.h>
  29
  30 // ==========
  31 //   face_t
  32 // ==========
  33
  34 PolyMesh::face_t::face_t(int N, int o, int n, vec3_t rv, int b)
  35 {
  36   node.resize(N);
  37   owner = o;
  38   neighbour = n;
  39   bc = b;
  40   ref_vec = rv;
  41 }
  42
  43
  44 // ==========
  45 //   node_t
  46 // ==========
  47
  48 PolyMesh::node_t::node_t(const QVector<vtkIdType> &ids)
  49 {
  50   id = ids;
  51   qSort(id);
  52 }
  53
  54 PolyMesh::node_t::node_t(vtkIdType id1)
  55 {
  56   id.resize(1);
  57   id[0] = id1;
  58 }
  59
  60 PolyMesh::node_t::node_t(vtkIdType id1, vtkIdType id2)
  61 {
  62   id.resize(2);
  63   id[0] = id1;
  64   id[1] = id2;
  65   qSort(id);
  66 }
  67
  68 PolyMesh::node_t::node_t(vtkIdType id1, vtkIdType id2, vtkIdType id3)
  69 {
  70   id.resize(3);
  71   id[0] = id1;
  72   id[1] = id2;
  73   id[2] = id3;
  74   qSort(id);
  75 }
  76
  77 PolyMesh::node_t::node_t(vtkIdType id1, vtkIdType id2, vtkIdType id3, vtkIdType id4)
  78 {
  79   id.resize(4);
  80   id[0] = id1;
  81   id[1] = id2;
  82   id[2] = id3;
  83   id[3] = id4;
  84   qSort(id);
  85 }
  86
  87 bool PolyMesh::node_t::operator<(const PolyMesh::node_t &N) const
  88 {
  89   int num = min(id.size(), N.id.size());
  90   for (int i = 0; i < num; ++i) {
  91     if (id[i] < N.id[i]) {
  92       return true;
  93     }
  94     if (id[i] > N.id[i]) {
  95       return false;
  96     }
  97   }
  98   if (id.size() < N.id.size()) {
  99     return true;
 100   }
 101   return false;
 102 }
 103
 104 bool PolyMesh::node_t::operator>(const PolyMesh::node_t &N) const
 105 {
 106   int num = min(id.size(), N.id.size());
 107   for (int i = 0; i < num; ++i) {
 108     if (id[i] > N.id[i]) {
 109       return true;
 110     }
 111     if (id[i] < N.id[i]) {
 112       return false;
 113     }
 114   }
 115   if (id.size() > N.id.size()) {
 116     return true;
 117   }
 118   return false;
 119 }
 120
 121 bool PolyMesh::node_t::operator==(const PolyMesh::node_t &N) const
 122 {
 123   if (id.size() != N.id.size()) {
 124     return false;
 125   }
 126   for (int i = 0; i < id.size(); ++i) {
 127     if (id[i] != N.id[i]) {
 128       return false;
 129     }
 130   }
 131   return true;
 132 }
 133
 134
 135
 136
 137
 138 // ============
 139 //   PolyMesh
 140 // ============
 141
 142 PolyMesh::PolyMesh(vtkUnstructuredGrid *grid, bool dualise, double pull_in, bool optimise, bool split_faces, bool split_cells)
 143 {
 144   m_CreateDualMesh = false;
 145   m_OptimiseConvexity = false;
 146   m_SplitCells = false;
 147   m_SplitFaces = split_faces;
 148   if (dualise) {
 149     for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
 150       if (isVolume(id_cell, grid) && grid->GetCellType(id_cell) != VTK_POLYHEDRON) {
 151         m_CreateDualMesh = true;
 152         break;
 153       }
 154     }
 155     if (m_CreateDualMesh) {
 156       m_OptimiseConvexity = optimise;
 157
 158     }
 159   }
 160
 161   m_AttractorWeight = 0.0;
 162   m_PullInFactor = pull_in;
 163   m_OptimiseConvexity = optimise && m_CreateDualMesh;
 164   m_Grid = grid;
 165   m_Part.setGrid(m_Grid);
 166   m_Part.setAllCells();
 167   findPolyCells();
 168   createNodesAndFaces();
 169   checkFaceOrientation();
 170   buildPoint2Face();
 171   buildPCell2Face();
 172   computePoints();
 173
 174   if (m_OptimiseConvexity) {
 175     for (int iter = 0; iter < 2; ++iter) {
 176       int num_bad = 0;
 177       int i_improve = 0;
 178       for (int i = 0; i < numCells(); ++i) {
 179
 180         // check if any of the faces is a boundary face
 181         bool no_boundary = true;
 182         for (int j = 0; j < numFacesOfPCell(i); ++j) {
 183           int bc = boundaryCode(pcell2Face(i, j));
 184           if (bc != 0) {
 185             no_boundary = false;
 186             break;
 187           }
 188         }
 189
 190         if (no_boundary) {
 191           PolyMolecule pm(this, i);
 192           if (!pm.allPositive()) {
 193             ++i_improve;
 194             pm.optimise();
 195             if (pm.minPyramidVolume() < 0) {
 196               ++num_bad;
 197             }
 198           }
 199         }
 200       }
 201       cout << i_improve << " cells out of " << numCells() << " were optimised." << endl;
 202       if (num_bad == 0) {
 203         break;
 204       }
 205     }
 206   }
 207
 208   if (m_SplitFaces) {
 209     splitConcaveFaces();
 210   }
 211
 212   if (m_SplitCells) {
 213     for (int iter = 0; iter < 0; ++iter) {
 214       int num_improved = 0;
 215       int num_cells = numCells();
 216       for (int i = 0; i < num_cells; ++i) {
 217         PolyMolecule pm(this, i);
 218         if (!pm.allPositive()) {
 219           ++num_improved;
 220           pm.centreSplit();
 221         }
 222       }
 223       buildPoint2Face();
 224       buildPCell2Face();
 225       cout << num_improved << " cells out of " << num_cells << " were split." << endl;
 226       if (num_improved == 0) {
 227         break;
 228       }
 229     }
 230   }
 231   sortFaces();
 232   buildPoint2Face();
 233   buildPCell2Face();
 234 }
 235
 236 void PolyMesh::merge(PolyMesh *poly)
 237 {
 238   foreach (face_t face, poly->m_Faces) {
 239     for (int i = 0; i < face.node.size(); ++i) {
 240       face.node[i] += m_Points.size();
 241     }
 242     face.owner += m_NumPolyCells;
 243     face.neighbour += m_NumPolyCells;
 244     m_Faces.append(face);
 245   }
 246   m_Points += poly->m_Points;
 247   m_NumPolyCells += poly->m_NumPolyCells;
 248   sortFaces();
 249   collectBoundaryConditions();
 250 }
 251
 252 void PolyMesh::triangulateBadFaces()
 253 {
 254   m_IsBadCell.fill(false, numCells());
 255   int num_bad = 0;
 256   for (int i = 0; i < numCells(); ++i) {
 257     PolyMolecule pm(this, i);
 258     if (!pm.allPositive()) {
 259       ++num_bad;
 260       m_IsBadCell[i] = true;
 261     }
 262   }
 263   QList<int> bad_faces;
 264   for (int i = 0; i < numFaces(); ++i) {
 265     if (m_IsBadCell[owner(i)]) {
 266       bad_faces.append(i);
 267     } else if (neighbour(i) != -1) {
 268       if (m_IsBadCell[neighbour(i)]) {
 269         bad_faces.append(i);
 270       }
 271     }
 272   }
 273
 274   foreach (int i_face, bad_faces) {
 275     face_t face = m_Faces[i_face];
 276     QVector<vec3_t> x(face.node.size());
 277     for (int i = 0; i < x.size(); ++i) {
 278       x[i] = nodeVector(face.node[i]);
 279     }
 280     EG_VTKSP(vtkPolyData, poly);
 281     createPolyData(x, poly);
 282     EG_VTKSP(vtkTriangleFilter, tri);
 283     tri->SetInputData(poly);
 284     tri->Update();
 285     if (tri->GetOutput()->GetNumberOfPoints() > face.node.size()) {
 286       EG_BUG;
 287     }
 288     QVector<face_t> new_faces(tri->GetOutput()->GetNumberOfCells(), face);
 289     for (int i = 0; i < new_faces.size(); ++i) {
 290       new_faces[i].node.resize(3);
 291       EG_GET_CELL(i, tri->GetOutput());
 292       for (int j = 0; j < 3; ++j) {
 293         new_faces[i].node[j] = face.node[pts[j]];
 294       }
 295     }
 296     m_Faces[i_face] = new_faces[0];
 297     for (int i = 1; i < new_faces.size(); ++i) {
 298       m_Faces.append(new_faces[i]);
 299     }
 300   }
 301
 302   sortFaces();
 303   buildPoint2Face();
 304   buildPCell2Face();
 305
 306   cout << num_bad << " cells out of " << numCells() << " are concave" << endl;
 307   cout << bad_faces.size() << " faces have been triangulated" << endl;
 308 }
 309
 310 /*
 311 void PolyMesh::splitConcaveFaces()
 312 {
 313   for (int i_cell = 0; i_cell < numCells(); ++i_cell) {
 314     if (m_IsBadCell[i_cell]) {
 315       PolyMolecule pm(this, i_cell);
 316       EG_VTKSP(vtkPolyData, poly);
 317       pm.createPolyData(poly);
 318       EG_VTKSP(vtkCurvatures, curv);
 319       curv->SetCurvatureTypeToMinimum();
 320       curv->SetInput(poly);
 321       curv->Update();
 322       EG_VTKDCN(vtkDoubleArray, curvature, curv->GetOutput(), "Minimum_Curvature");
 323       vtkIdType id_ref = 0;
 324       vec3_t x_ref;
 325       double curv_min = 1e99;
 326       EG_FORALL_NODES(id_node, poly) {
 327         if (curvature->GetValue(id_node) < curv_min) {
 328           curv_min = curvature->GetValue(id_node);
 329           id_ref = id_node;
 330           poly->GetPoint(id_ref, x_ref.data());
 331         }
 332       }
 333       vec3_t n_ref;
 334       EG_FORALL_CELLS(id_face, poly) {
 335         EG_GET_CELL(id_face, poly);
 336         bool found = false;
 337         if (type_cell != VTK_TRIANGLE) {
 338           EG_BUG;
 339         }
 340         for (int i = 0; i < num_pts; ++i) {
 341           if (pts[i] == id_ref) {
 342             found = true;
 343             break;
 344           }
 345         }
 346         if (found) {
 347           QVector<vec3_t> x(num_pts);
 348           for (int i = 0; i < num_pts; ++i) {
 349             poly->GetPoint(pts[i], x[i].data());
 350           }
 351           n_ref = (x[1] - x[0]);
 352           n_ref = n_ref.cross(x[2] - x[0]);
 353           n_ref.normalise();
 354           break;
 355         }
 356       }
 357
 358
 359       EG_VTKSP(vtkDoubleArray, div);
 360       div->SetNumberOfComponents(1);
 361       div->SetNumberOfTuples(poly->GetNumberOfPoints());
 362       div->SetName("div");
 363
 364       bool done = false;
 365       vtkIdType id_node0 = id_ref;
 366       vtkIdType id_node1 = id_ref;
 367       QVector<QSet<vtkIdType> > n2n;
 368       QVector<QSet<vtkIdType> > n2c;
 369       QVector<QVector<vtkIdType> > c2c;
 370       createPolyDataN2N(poly, n2n);
 371       createPolyDataN2C(poly, n2c);
 372       createPolyDataC2C(poly, c2c);
 373
 374       vec3_t n_node(0,0,0);
 375       foreach (vtkIdType id_face, n2c[id_ref]) {
 376         EG_GET_CELL(id_face, poly);
 377         QVector<vec3_t> x(3);
 378         for (int i = 0; i < 3; ++i) {
 379           poly->GetPoint(pts[i], x[i].data());
 380         }
 381         n_node += triNormal(x[0], x[1], x[2]);
 382       }
 383       n_node.normalise();
 384
 385       QVector<bool> div_node(poly->GetNumberOfPoints(), false);
 386       EG_FORALL_NODES(id_node, poly) {
 387         div->SetValue(id_node, 0);
 388       }
 389       div_node[id_ref] = true;
 390       double dv = 1.0;
 391       div->SetValue(id_ref, dv);
 392       dv += 0.2;
 393       int count = 0;
 394       double c_min = 1e99;
 395       double L = 1e99;
 396       foreach (vtkIdType id, n2n[id_node0]) {
 397         vec3_t x;
 398         poly->GetPoint(id, x.data());
 399         L = min(L, (x - x_ref).abs());
 400         double c = curvature->GetValue(id);
 401         if (c < c_min) {
 402           c_min = c;
 403           id_node1 = id;
 404         }
 405       }
 406       div_node[id_node1] = true;
 407       div->SetValue(id_node1, dv);
 408       dv += 0.2;
 409       while (!done) {
 410         vtkIdType id_node2 = id_node1;
 411         double d_min = 1e99;
 412         foreach (vtkIdType id, n2n[id_node1]) {
 413           if (id == id_ref && id_node0 != id_ref) {
 414             done = true;
 415             break;
 416           }
 417         }
 418         if (!done) {
 419           foreach (vtkIdType id, n2n[id_node1]) {
 420             if (!div_node[id]) {
 421               vec3_t x;
 422               poly->GetPoint(id, x.data());
 423               if ((x - x_ref)*n_node < 0.25*L) {
 424                 double d = (x - x_ref)*n_ref;
 425                 if (fabs(d) < d_min && d < 0.1*L) {
 426                   d_min = d;
 427                   id_node2 = id;
 428                 }
 429               }
 430             }
 431           }
 432           div_node[id_node2] = true;
 433           div->SetValue(id_node2, dv);
 434           dv += 0.2;
 435           id_node0 = id_node1;
 436           id_node1 = id_node2;
 437           ++count;
 438           if (count == 100) {
 439             EG_BUG;
 440           }
 441         }
 442       }
 443       poly->GetPointData()->AddArray(div);
 444       EG_VTKSP(vtkXMLPolyDataWriter, vtp);
 445       vtp->SetFileName(qPrintable(GuiMainWindow::pointer()->getCwd() + "/triangulated_cell.vtp"));
 446       vtp->SetInput(poly);
 447       vtp->SetDataModeToAscii();
 448       vtp->Write();
 449       EG_BUG;
 450     }
 451   }
 452 }
 453 */
 454
 455 void PolyMesh::getFacesOfEdgeInsideCell(vtkIdType id_cell, vtkIdType id_node1, vtkIdType id_node2, int &face1, int &face2)
 456 {
 457   vtkIdType num_pts, *pts;
 458   m_Grid->GetCellPoints(id_cell, num_pts, pts);
 459   face1 = -1;
 460   face2 = -1;
 461
 462   if (m_Grid->GetCellType(id_cell) == VTK_TETRA) {
 463     if        (id_node1 == pts[0]) {
 464       if      (id_node2 == pts[1]) { face1 = 0; face2 = 1; }
 465       else if (id_node2 == pts[2]) { face1 = 2; face2 = 0; }
 466       else if (id_node2 == pts[3]) { face1 = 1; face2 = 2; }
 467     } else if (id_node1 == pts[1]) {
 468       if      (id_node2 == pts[0]) { face1 = 1; face2 = 0; }
 469       else if (id_node2 == pts[2]) { face1 = 0; face2 = 3; }
 470       else if (id_node2 == pts[3]) { face1 = 3; face2 = 1; }
 471     } else if (id_node1 == pts[2]) {
 472       if      (id_node2 == pts[0]) { face1 = 0; face2 = 2; }
 473       else if (id_node2 == pts[1]) { face1 = 3; face2 = 0; }
 474       else if (id_node2 == pts[3]) { face1 = 2; face2 = 3; }
 475     } else if (id_node1 == pts[3]) {
 476       if      (id_node2 == pts[0]) { face1 = 2; face2 = 1; }
 477       else if (id_node2 == pts[1]) { face1 = 1; face2 = 3; }
 478       else if (id_node2 == pts[2]) { face1 = 3; face2 = 2; }
 479     }
 480   } else if (m_Grid->GetCellType(id_cell) == VTK_PYRAMID) {
 481     if        (id_node1 == pts[0]) {
 482       if      (id_node2 == pts[1]) { face1 = 0; face2 = 1; }
 483       else if (id_node2 == pts[3]) { face1 = 4; face2 = 0; }
 484       else if (id_node2 == pts[4]) { face1 = 1; face2 = 4; }
 485     } else if (id_node1 == pts[1]) {
 486       if      (id_node2 == pts[0]) { face1 = 1; face2 = 0; }
 487       else if (id_node2 == pts[2]) { face1 = 0; face2 = 2; }
 488       else if (id_node2 == pts[4]) { face1 = 2; face2 = 1; }
 489     } else if (id_node1 == pts[2]) {
 490       if      (id_node2 == pts[1]) { face1 = 2; face2 = 0; }
 491       else if (id_node2 == pts[3]) { face1 = 0; face2 = 3; }
 492       else if (id_node2 == pts[4]) { face1 = 3; face2 = 2; }
 493     } else if (id_node1 == pts[3]) {
 494       if      (id_node2 == pts[0]) { face1 = 0; face2 = 4; }
 495       else if (id_node2 == pts[2]) { face1 = 3; face2 = 0; }
 496       else if (id_node2 == pts[4]) { face1 = 4; face2 = 3; }
 497     } else if (id_node1 == pts[4]) {
 498       if      (id_node2 == pts[0]) { face1 = 4; face2 = 1; }
 499       else if (id_node2 == pts[1]) { face1 = 1; face2 = 2; }
 500       else if (id_node2 == pts[2]) { face1 = 2; face2 = 3; }
 501       else if (id_node2 == pts[3]) { face1 = 3; face2 = 4; }
 502     }
 503   } else {
 504     EG_BUG;
 505   }
 506
 507   if (face1 == -1 || face2 == -1) {
 508     EG_BUG;
 509   }
 510 }
 511
 512 void PolyMesh::getSortedEdgeCells(vtkIdType id_node1, vtkIdType id_node2, QList<vtkIdType> &cells, bool &is_loop)
 513 {
 514   cells.clear();
 515   vtkIdType id_start = -1;
 516   for (int i = 0; i < m_Part.n2cGSize(id_node1); ++i) {
 517     vtkIdType id_cell = m_Part.n2cGG(id_node1, i);
 518     if (isVolume(id_cell, m_Grid)) {
 519       if (m_Cell2PCell[id_cell] == -1) {
 520         QVector<vtkIdType> pts;
 521         getPointsOfCell(m_Grid, id_cell, pts);
 522         if (pts.contains(id_node2)) {
 523           id_start = id_cell;
 524         }
 525       }
 526     }
 527   }
 528   if (id_start == -1) {
 529     EG_BUG;
 530   }
 531   vtkIdType id_cell_old = id_start;
 532   vtkIdType id_cell_new = id_start;
 533   cells.append(id_cell_new);
 534   bool finished = false;
 535   bool prepend = false;
 536   while (!finished) {
 537     do {
 538       int f1, f2;
 539       getFacesOfEdgeInsideCell(id_cell_new, id_node1, id_node2, f1, f2);
 540       vtkIdType id_neigh1 = m_Part.c2cGG(id_cell_new, f1);
 541       vtkIdType id_neigh2 = m_Part.c2cGG(id_cell_new, f2);
 542       if (id_neigh1 == -1) EG_BUG;
 543       if (id_neigh2 == -1) EG_BUG;
 544       if (id_neigh1 == id_cell_old) {
 545         id_cell_old = id_cell_new;
 546         id_cell_new = id_neigh2;
 547       } else {
 548         id_cell_old = id_cell_new;
 549         id_cell_new = id_neigh1;
 550       }
 551       if (prepend) {
 552         cells.prepend(id_cell_new);
 553       } else {
 554         cells.append(id_cell_new);
 555       }
 556     } while (id_cell_new != id_start  &&  !isSurface(id_cell_new, m_Grid)  &&  m_Cell2PCell[id_cell_new] == -1);
 557     if ((isSurface(id_cell_new, m_Grid) || m_Cell2PCell[id_cell_new] != -1) && !prepend) {
 558       id_cell_new = cells[0];
 559       id_cell_old = cells[1];
 560       prepend = true;
 561     } else {
 562       finished = true;
 563     }
 564   }
 565
 566   // remove last cell for loops
 567   if (cells.size() > 1 && cells.first() == cells.last() && m_Cell2PCell[cells.first()] == -1) {
 568     cells.pop_back();
 569     is_loop = true;
 570   } else {
 571     is_loop = false;
 572   }
 573 }
 574
 575 bool PolyMesh::isDualFace(vtkIdType id_face)
 576 {
 577   vtkIdType id_cell = m_Part.getVolumeCell(id_face);
 578   if (m_Cell2PCell[id_cell] == -1) {
 579     return true;
 580   }
 581   return false;
 582 }
 583
 584 void PolyMesh::getSortedPointFaces(vtkIdType id_node, int bc, QList<vtkIdType> &faces, bool &is_loop)
 585 {
 586   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
 587   faces.clear();
 588   vtkIdType id_start = -1;
 589   for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
 590     vtkIdType id_cell = m_Part.n2cGG(id_node, i);
 591     if (isSurface(id_cell, m_Grid)) {
 592       if (cell_code->GetValue(id_cell) == bc && isDualFace(id_cell)) {
 593         id_start = id_cell;
 594         break;
 595       }
 596     }
 597   }
 598   if (id_start == -1) {
 599     return;
 600   }
 601   vtkIdType id_face_old = id_start;
 602   vtkIdType id_face_new = id_start;
 603   faces.append(id_face_new);
 604   bool finished = false;
 605   bool prepend = false;
 606   while (!finished) {
 607     do {
 608       QList<vtkIdType> id_neigh;
 609       for (int i = 0; i < m_Part.c2cGSize(id_face_new); ++i) {
 610         if (cellContainsNode(m_Grid, m_Part.c2cGG(id_face_new, i), id_node)) {
 611           id_neigh.append(m_Part.c2cGG(id_face_new, i));
 612         }
 613       }
 614       if (id_neigh[0] == id_face_old) {
 615         id_face_old = id_face_new;
 616         id_face_new = id_neigh[1];
 617       } else {
 618         id_face_old = id_face_new;
 619         id_face_new = id_neigh[0];
 620       }
 621       if (cell_code->GetValue(id_face_new) == bc  &&  isDualFace(id_face_new)) {
 622         if (prepend) {
 623           faces.prepend(id_face_new);
 624         } else {
 625           faces.append(id_face_new);
 626         }
 627       }
 628     } while (id_face_new != id_start  &&  cell_code->GetValue(id_face_new) == bc  &&  isDualFace(id_face_new));
 629     if ((cell_code->GetValue(id_face_new) != bc || !isDualFace(id_face_new))  &&  !prepend) {
 630       id_face_old = id_face_new;
 631       id_face_new = faces[0];
 632       if (faces.size() > 1) {
 633         id_face_old = faces[1];
 634       }
 635       prepend = true;
 636     } else {
 637       finished = true;
 638     }
 639   }
 640
 641   // remove last face for loops
 642   if (faces.size() > 1 && faces.first() == faces.last()) {
 643     faces.pop_back();
 644     is_loop = true;
 645   } else {
 646     is_loop = false;
 647   }
 648
 649 }
 650
 651 void PolyMesh::findPolyCells()
 652 {
 653   m_Cell2PCell.fill(-1, m_Grid->GetNumberOfCells());
 654   m_Node2PCell.fill(-1, m_Grid->GetNumberOfPoints());
 655   m_NumPolyCells = 0;
 656   if (m_CreateDualMesh) {
 657     for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 658       if (!isHexCoreNode(id_node)) {
 659         bool tetra_or_pyramid = false;
 660         for (int j = 0; j < m_Part.n2cGSize(id_node); ++j) {
 661           vtkIdType id_cell = m_Part.n2cGG(id_node, j);
 662           vtkIdType type_cell = m_Grid->GetCellType(id_cell);
 663           if (type_cell == VTK_TETRA || type_cell == VTK_PYRAMID) {
 664             tetra_or_pyramid = true;
 665           }
 666         }
 667         if (tetra_or_pyramid) {
 668           m_Node2PCell[id_node] = m_NumPolyCells;
 669           ++m_NumPolyCells;
 670         }
 671       }
 672     }
 673   }
 674   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 675     vtkIdType type_cell = m_Grid->GetCellType(id_cell);
 676     if (isVolume(id_cell, m_Grid)) {
 677       if (type_cell == VTK_WEDGE || type_cell == VTK_HEXAHEDRON || type_cell == VTK_POLYHEDRON || !m_CreateDualMesh) {
 678         m_Cell2PCell[id_cell] = m_NumPolyCells;
 679         ++m_NumPolyCells;
 680       }
 681     }
 682   }
 683   m_CellCentre.resize(m_NumPolyCells);
 684   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 685     if (m_Node2PCell[id_node] != -1) {
 686       m_Grid->GetPoint(id_node, m_CellCentre[m_Node2PCell[id_node]].data());
 687     }
 688   }
 689   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 690     if (m_Cell2PCell[id_cell] != -1) {
 691       m_CellCentre[m_Cell2PCell[id_cell]] = cellCentre(m_Grid, id_cell);
 692     }
 693   }
 694 }
 695
 696 void PolyMesh::createFace(QList<node_t> nodes, int owner, int neighbour, vec3_t ref_vec, int bc)
 697 {
 698   if (owner > neighbour && neighbour != -1) {
 699     swap(owner, neighbour);
 700     ref_vec *= -1;
 701   }
 702   face_t face(nodes.size(), owner, neighbour, ref_vec, bc);
 703   for (int i = 0; i < nodes.size(); ++i) {
 704     int idx = m_Nodes.insert(nodes[i]);
 705     face.node[i] = idx;
 706   }
 707   m_Faces.append(face);
 708 }
 709
 710 void PolyMesh::createCornerFace(vtkIdType id_cell, int i_face, vtkIdType id_node)
 711 {
 712   QList<vtkIdType> edge_nodes;
 713   QVector<vtkIdType> face_nodes;
 714   getFaceOfCell(m_Grid, id_cell, i_face, face_nodes);
 715   vtkIdType num_pts, *pts;
 716   m_Grid->GetCellPoints(id_cell, num_pts, pts);
 717   for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
 718     vtkIdType id_neigh = m_Part.n2nGG(id_node, i);
 719     if (face_nodes.contains(id_neigh)) {
 720       edge_nodes.append(id_neigh);
 721     }
 722   }
 723   if (edge_nodes.size() != 2) {
 724     EG_BUG;
 725   }
 726   int owner = m_Cell2PCell[id_cell];
 727   if (owner == -1) {
 728     EG_BUG;
 729   }
 730   int neighbour = m_Node2PCell[id_node];
 731   int bc = 0;
 732   if (neighbour == -1) {
 733     EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
 734     vtkIdType id_face = m_Part.c2cGG(id_cell, i_face);
 735     if (id_face == -1) {
 736       EG_BUG;
 737     }
 738     if (!isSurface(id_cell, m_Grid)) {
 739       EG_BUG;
 740     }
 741     bc = cell_code->GetValue(id_face);
 742   }
 743   QList<node_t> nodes;
 744   nodes.append(node_t(id_node));
 745   nodes.append(node_t(id_node, edge_nodes[0]));
 746   nodes.append(node_t(face_nodes));
 747   nodes.append(node_t(id_node, edge_nodes[1]));
 748   vec3_t n = getNormalOfCell(m_Grid, id_cell, i_face);
 749   n.normalise();
 750   createFace(nodes, owner, neighbour, n, bc);
 751 }
 752
 753 void PolyMesh::createEdgeFace(vtkIdType id_node1, vtkIdType id_node2)
 754 {
 755   // check if additional edge node needs to be created
 756   // (transition from boundary layer to far-field)
 757   // try to find a boundary layer cell which contains both nodes
 758   bool add_edge_node = false;
 759   for (int i = 0; i < m_Part.n2cGSize(id_node1); ++i) {
 760     vtkIdType id_cell = m_Part.n2cGG(id_node1, i);
 761     if (m_Cell2PCell[id_cell] != -1) {
 762       if (cellContainsNode(m_Grid, id_cell, id_node1) && cellContainsNode(m_Grid, id_cell, id_node2)) {
 763         add_edge_node = true;
 764         break;
 765       }
 766     }
 767   }
 768   if (!add_edge_node) {
 769     for (int i = 0; i < m_Part.n2cGSize(id_node2); ++i) {
 770       vtkIdType id_cell = m_Part.n2cGG(id_node2, i);
 771       if (m_Cell2PCell[id_cell] != -1) {
 772         if (cellContainsNode(m_Grid, id_cell, id_node1) && cellContainsNode(m_Grid, id_cell, id_node2)) {
 773           add_edge_node = true;
 774           break;
 775         }
 776       }
 777     }
 778   }
 779
 780   QList<vtkIdType> cells;
 781   bool loop = false;
 782   getSortedEdgeCells(id_node1, id_node2, cells, loop);
 783   int owner     = m_Node2PCell[id_node1];
 784   int neighbour = m_Node2PCell[id_node2];
 785   if (owner == -1) {
 786     EG_BUG;
 787   }
 788   if (neighbour == -1) {
 789     EG_BUG;
 790   }
 791   if (owner > neighbour) {
 792     swap(id_node1, id_node2);
 793     swap(owner, neighbour);
 794   }
 795   QList<node_t> nodes;
 796   for (int i_cells = 0; i_cells < cells.size(); ++i_cells) {
 797     node_t node;
 798     QVector<vtkIdType> pts1;
 799     getPointsOfCell(m_Grid, cells[i_cells], pts1);
 800     if (m_Cell2PCell[cells[i_cells]] != -1) {
 801       int i_cells2 = 0;
 802       if (i_cells == 0) {
 803         i_cells2 = 1;
 804       } else if (i_cells == cells.size() - 1) {
 805         i_cells2 = cells.size() - 2;
 806       } else {
 807         EG_BUG;
 808       }
 809       QVector<vtkIdType> pts2;
 810       getPointsOfCell(m_Grid, cells[i_cells2], pts2);
 811       QSet<vtkIdType> p1, p2;
 812       for (int i_pts1 = 0; i_pts1 < pts1.size(); ++i_pts1) {
 813         p1.insert(pts1[i_pts1]);
 814       }
 815       for (int i_pts2 = 0; i_pts2 < pts2.size(); ++i_pts2) {
 816         p2.insert(pts2[i_pts2]);
 817       }
 818       QSet<vtkIdType> face_nodes = p1.intersect(p2);
 819       node.id.resize(face_nodes.size());
 820       int i = 0;
 821       foreach (vtkIdType id_node, face_nodes) {
 822         node.id[i] = id_node;
 823         ++i;
 824       }
 825     } else {
 826       node.id.resize(pts1.size());
 827       for (int i_pts1 = 0; i_pts1 < pts1.size(); ++i_pts1) {
 828         node.id[i_pts1] = pts1[i_pts1];
 829       }
 830     }
 831     qSort(node.id);
 832     nodes.append(node);
 833   }
 834   if (!loop) {
 835     EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
 836     if (cell_code->GetValue(cells.first()) != cell_code->GetValue(cells.last()) || add_edge_node) {
 837       nodes.append(node_t(id_node1, id_node2));
 838     }
 839   }
 840   vec3_t x1, x2;
 841   m_Grid->GetPoint(id_node1, x1.data());
 842   m_Grid->GetPoint(id_node2, x2.data());
 843   createFace(nodes, owner, neighbour, x2 - x1, 0);
 844 }
 845
 846 void PolyMesh::createFaceFace(vtkIdType id_cell, int i_face)
 847 {
 848   if (m_Cell2PCell[id_cell] == -1) {
 849     EG_BUG;
 850   }
 851   int owner = m_Cell2PCell[id_cell];
 852   int neighbour = -1;
 853   vtkIdType id_neigh = m_Part.c2cGG(id_cell, i_face);
 854   int bc = 0;
 855   if (id_neigh != -1) {
 856     if (isVolume(id_neigh, m_Grid)) {
 857       if (m_Cell2PCell[id_neigh] == -1) {
 858         EG_BUG;
 859       }
 860       neighbour = m_Cell2PCell[id_neigh];
 861     } else {
 862       EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
 863       bc = cell_code->GetValue(id_neigh);
 864     }
 865   }
 866   QVector<vtkIdType> tmp_node_ids;
 867   getFaceOfCell(m_Grid, id_cell, i_face, tmp_node_ids);
 868   vec3_t n = getNormalOfCell(m_Grid, id_cell, i_face);
 869   n.normalise();
 870   QVector<vtkIdType> node_ids(tmp_node_ids.size() + 1);
 871   for (int i = 0; i < tmp_node_ids.size(); ++i) {
 872     node_ids[i] = tmp_node_ids[i];
 873   }
 874   node_ids[node_ids.size() - 1] = node_ids[0];
 875   QList<node_t> nodes;
 876   for (int i = 0; i < node_ids.size() - 1; ++i) {
 877     nodes.append(node_t(node_ids[i]));
 878     if (m_Node2PCell[node_ids[i]] != -1 && m_Node2PCell[node_ids[i+1]] != -1) {
 879       nodes.append(node_t(node_ids[i], node_ids[i+1]));
 880     }
 881   }
 882   createFace(nodes, owner, neighbour, n, bc);
 883 }
 884
 885 void PolyMesh::createPointFace(vtkIdType id_node, int bc)
 886 {
 887   bool is_loop;
 888   QList<vtkIdType> faces;
 889   getSortedPointFaces(id_node, bc, faces, is_loop);
 890   if (faces.size() == 0) {
 891     return;
 892   }
 893   QList<node_t> nodes;
 894   vec3_t n(0,0,0);
 895   if (faces.size() == 0) {
 896     EG_BUG;
 897   }
 898   foreach (vtkIdType id_face, faces) {
 899     node_t node;
 900     vtkIdType num_pts, *pts;
 901     m_Grid->GetCellPoints(id_face, num_pts, pts);
 902     node.id.resize(num_pts);
 903     for (int i_pts = 0; i_pts < num_pts; ++i_pts) {
 904       node.id[i_pts] = pts[i_pts];
 905     }
 906     qSort(node.id);
 907     nodes.append(node);
 908     n += GeometryTools::cellNormal(m_Grid, id_face);
 909   }
 910   n.normalise();
 911   if (!is_loop) {
 912     bool prepend = false;
 913     vtkIdType id_face = faces.last();
 914     while (id_face != -1) {
 915       bool found = false;
 916       vtkIdType num_pts, *pts;
 917       m_Grid->GetCellPoints(id_face, num_pts, pts);
 918       QList<vtkIdType> id_neigh_node;
 919       vtkIdType id_neigh = -1;
 920       for (int i = 0; i < m_Part.c2cGSize(id_face); ++i) {
 921         id_neigh = m_Part.c2cGG(id_face, i);
 922         if (id_neigh == -1) {
 923           EG_BUG;
 924         }
 925         if (!isSurface(id_neigh, m_Grid)) {
 926           EG_BUG;
 927         }
 928         if (cellContainsNode(m_Grid, id_neigh, id_node)) {
 929           if (!faces.contains(id_neigh)) {
 930             if (found) {
 931               EG_BUG;
 932             }
 933             for (int j = 0; j < num_pts; ++j) {
 934               if (pts[j] != id_node && cellContainsNode(m_Grid, id_neigh, pts[j])) {
 935                 id_neigh_node.append(pts[j]);
 936               }
 937             }
 938           }
 939         }
 940       }
 941
 942       if (id_neigh_node.size() == 0) {
 943         EG_BUG;
 944       }
 945
 946       if (id_neigh_node.size() == 1) {
 947         if (prepend) {
 948           nodes.prepend(node_t(id_node, id_neigh_node[0]));
 949           id_face = -1;
 950         } else {
 951           nodes.append(node_t(id_node, id_neigh_node[0]));
 952           prepend = true;
 953           id_face = faces.first();
 954         }
 955       } else {
 956         if (id_neigh_node.size() > 2) {
 957           EG_BUG;
 958         }
 959         if (faces.size() != 1) {
 960           EG_BUG;
 961         }
 962         nodes.prepend(node_t(id_node, id_neigh_node[0]));
 963         nodes.append(node_t(id_node, id_neigh_node[1]));
 964         id_face = -1;
 965       }
 966     }
 967
 968     // check if this is a transition node (boundary layer -> far-field)
 969     bool dual_cell_found = false;
 970     bool cell_cell_found = false;
 971     for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
 972       if (m_Cell2PCell[m_Part.n2cGG(id_node,i)] == -1) {
 973         cell_cell_found = true;
 974       } else {
 975         dual_cell_found = true;
 976       }
 977     }
 978
 979     if (m_Part.n2bcGSize(id_node) > 2 || (dual_cell_found && cell_cell_found)) {
 980       nodes.prepend(node_t(id_node));
 981     }
 982   }
 983   int owner     = m_Node2PCell[id_node];
 984   int neighbour = -1;
 985   createFace(nodes, owner, neighbour, n, bc);
 986 }
 987
 988 void PolyMesh::computePoints()
 989 {
 990   QVector<node_t> nodes;
 991   m_Nodes.getQVector(nodes);
 992   m_Points.resize(nodes.size());
 993   m_PointWeights.fill(1.0, m_Grid->GetNumberOfPoints());
 994
 995   // find transition nodes
 996   QVector<bool> is_transition_node(m_Grid->GetNumberOfPoints(), false);
 997   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 998     if (m_Node2PCell[id_node] != -1) {
 999       for (int i_cell = 0; i_cell < m_Part.n2cGSize(id_node); ++i_cell) {
1000         vtkIdType id_cell = m_Part.n2cGG(id_node, i_cell);
1001         if (m_Cell2PCell[id_cell] != -1) {
1002           is_transition_node[id_node] = true;
1003           break;
1004         }
1005       }
1006     }
1007   }
1008
1009   // compute weights (attraction for convex corners)
1010   // mark convex edge nodes
1011   QVector<bool> is_convex_node(m_Grid->GetNumberOfPoints(), false);
1012   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
1013     if (m_Cell2PCell[id_cell] != -1) {
1014       vec3_t xc1 = cellCentre(m_Grid, id_cell);
1015       int N = 0;
1016       vec3_t n_out(0,0,0);
1017       for (int i_face = 0; i_face < m_Part.c2cGSize(id_cell); ++i_face) {
1018         vtkIdType id_neigh = m_Part.c2cGG(id_cell, i_face);
1019         if (id_neigh != -1) {
1020           if (m_Cell2PCell[id_neigh] == -1 && !isSurface(id_neigh, m_Grid)) {
1021             ++N;
1022             n_out += getNormalOfCell(m_Grid, id_cell, i_face);
1023           }
1024         }
1025       }
1026       if (N > 0) {
1027         n_out.normalise();
1028         for (int i_face = 0; i_face < m_Part.c2cGSize(id_cell); ++i_face) {
1029           vtkIdType id_neigh = m_Part.c2cGG(id_cell, i_face);
1030           if (id_neigh != -1) {
1031             if (m_Cell2PCell[id_neigh] != -1) {
1032               vec3_t xc2 = cellCentre(m_Grid, id_neigh);
1033               QVector<vtkIdType> face_nodes;
1034               getFaceOfCell(m_Grid, id_cell, i_face, face_nodes);
1035               if (face_nodes.size() == 0) {
1036                 EG_BUG;
1037               }
1038               vec3_t xf(0,0,0);
1039               foreach (vtkIdType id_node, face_nodes) {
1040                 vec3_t x;
1041                 m_Grid->GetPoint(id_node, x.data());
1042                 xf += x;
1043               }
1044               xf *= 1.0/face_nodes.size();
1045               vec3_t v1 = xc1 - xf;
1046               vec3_t v2 = xc2 - xf;
1047               v1.normalise();
1048               v2.normalise();
1049               if ((v1+v2)*n_out < 0) {
1050                 double w = 1 + m_AttractorWeight*(1.0 + v1*v2);
1051                 foreach (vtkIdType id_node, face_nodes) {
1052                   if (m_Node2PCell[id_node] != -1) {
1053                     m_PointWeights[id_node] = max(m_PointWeights[id_node], w);
1054                     if (v1*v2 > 0) {
1055                       //is_convex_node[id_node] = true;
1056                     }
1057                   }
1058                 }
1059               }
1060             }
1061           }
1062         }
1063       }
1064     }
1065   }
1066
1067   // mapping for simple nodes
1068   QVector<int> prime2dual(m_Grid->GetNumberOfPoints(), -1);
1069
1070   // compute the point locations
1071   for (int i = 0; i < nodes.size(); ++i) {
1072     if (nodes[i].id.size() == 0) {
1073       EG_BUG;
1074     }
1075     m_Points[i] = vec3_t(0,0,0);
1076     double weight_sum = 0.0;
1077     foreach (vtkIdType id, nodes[i].id) {
1078       vec3_t x;
1079       m_Grid->GetPoint(id, x.data());
1080       weight_sum += m_PointWeights[id];
1081       m_Points[i] += m_PointWeights[id]*x;
1082     }
1083     //m_Points[i] *= 1.0/nodes[i].id.size();
1084     m_Points[i] *= 1.0/weight_sum;
1085     if (nodes[i].id.size() == 1) {
1086       prime2dual[nodes[i].id[0]] = i;
1087     }
1088   }
1089
1090   // correct transition nodes (pull-in)
1091   for (vtkIdType id_node1 = 0; id_node1 < m_Grid->GetNumberOfPoints(); ++id_node1) {
1092     if (is_transition_node[id_node1]) {
1093       if (prime2dual[id_node1] == -1) {
1094         EG_BUG;
1095       }
1096       vtkIdType id_node2 = -1;
1097       bool pull_in = false;
1098       for (int i_neigh = 0; i_neigh < m_Part.n2nGSize(id_node1); ++i_neigh) {
1099         vtkIdType id_neigh = m_Part.n2nGG(id_node1, i_neigh);
1100         if (m_Node2PCell[id_neigh] == -1 && !is_transition_node[id_neigh]) {
1101           if (id_node2 == -1) {
1102             pull_in = true;
1103             id_node2 = id_neigh;
1104           } else {
1105             pull_in = false;
1106           }
1107         }
1108       }
1109       if (pull_in) {
1110         double w = m_PullInFactor;
1111         m_Points[prime2dual[id_node1]] = w*m_Points[prime2dual[id_node2]] + (1-w)*m_Points[prime2dual[id_node1]];
1112       }
1113     }
1114   }
1115
1116 }
1117
1118 void PolyMesh::splitConcaveFaces()
1119 {
1120   QList<int> delete_faces;
1121   do {
1122     delete_faces.clear();
1123     QList<face_t> new_faces;
1124     for (int i_face = 0; i_face < m_Faces.size(); ++i_face) {
1125       face_t face = m_Faces[i_face];
1126       int num_nodes = face.node.size();
1127       if (num_nodes >= 4) {
1128         QVector<vec3_t> x_face(num_nodes);
1129         for (int i = 0; i < num_nodes; ++i) {
1130           x_face[i] = nodeVector(face.node[i]);
1131         }
1132         vec3_t xc, n;
1133         GeometryTools::planeFit(x_face, xc, n);
1134         QVector<vec3_t> x(num_nodes + 2);
1135         for (int i = 0; i < num_nodes; ++i) {
1136           x[i+1] = x_face[i];
1137         }
1138         x.first() = x_face.last();
1139         x.last() = x_face.first();
1140         double L_max = 0.1;
1141         int i1 = -1;
1142         vec3_t v;
1143         for (int i = 1; i <= num_nodes; ++i) {
1144           vec3_t x_neigh = 0.5*(x[i-1] + x[i+1]);
1145           double Lnorm = (x[i-1] - x[i+1]).abs();
1146           double L = ((x_neigh - xc).abs() - (x[i] - xc).abs())/Lnorm;
1147           if (L > L_max) {
1148             i1 = i;
1149             L_max = L;
1150             v = x[i] - x_neigh;
1151           }
1152         }
1153         if (i1 != -1) {
1154           int i2 = -1;
1155           double alpha_min = 1e99;
1156           for (int i = 1; i <= num_nodes; ++i) {
1157             if ((i < i1 - 1 || i > i1 + 1) && (i < i1 + num_nodes -1 || i > i1 + num_nodes + 1)) {
1158               double alpha = GeometryTools::angle(x[i] - x[i1], v);
1159               if (alpha < alpha_min) {
1160                 i2 = i;
1161                 alpha_min = alpha;
1162               }
1163             }
1164           }
1165           if (i2 != -1) {
1166             --i1;
1167             --i2;
1168             if (i1 > i2) {
1169               swap(i1, i2);
1170             }
1171             delete_faces.append(i_face);
1172             QVector<int> nodes1;
1173             {
1174               int i = 0;
1175               while (i < num_nodes) {
1176                 nodes1.append(face.node[i]);
1177                 if (i == i1) {
1178                   i = i2;
1179                 } else {
1180                   ++i;
1181                 }
1182               }
1183             }
1184             QVector<int> nodes2;
1185             for (int i = i1; i <= i2; ++i) {
1186               nodes2.append(face.node[i]);
1187             }
1188             face_t face1(nodes1.size(), face.owner, face.neighbour, face.ref_vec, face.bc);
1189             face_t face2(nodes2.size(), face.owner, face.neighbour, face.ref_vec, face.bc);
1190             face1.node = nodes1;
1191             face2.node = nodes2;
1192             new_faces.append(face1);
1193             new_faces.append(face2);
1194           }
1195         }
1196       }
1197     }
1198     {
1199       int offset = 0;
1200       foreach (int i_face, delete_faces) {
1201         m_Faces.removeAt(i_face - offset);
1202         ++offset;
1203       }
1204     }
1205     m_Faces.append(new_faces);
1206     cout << delete_faces.size() << " concave faces have been split." << endl;
1207     delete_faces.clear();
1208   } while (delete_faces.size() > 0);
1209   sortFaces();
1210   buildPoint2Face();
1211   buildPCell2Face();
1212 }
1213
1214 void PolyMesh::collectBoundaryConditions()
1215 {
1216   QSet<int> bcs;
1217   foreach (face_t face, m_Faces) {
1218     if (face.bc != 0) {
1219       bcs.insert(face.bc);
1220     }
1221   }
1222   m_BCs.resize(bcs.size());
1223   qCopy(bcs.begin(), bcs.end(), m_BCs.begin());
1224 }
1225
1226 void PolyMesh::createNodesAndFaces()
1227 {
1228   m_Nodes.resize(m_Grid->GetNumberOfPoints());
1229   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
1230
1231   int num_pcells = 0;
1232   int num_dcells = 0;
1233   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
1234     if (m_Cell2PCell[id_cell] != -1) {
1235       ++num_pcells;
1236     }
1237   }
1238   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
1239     if (m_Node2PCell[id_node] != -1) {
1240       ++num_dcells;
1241     }
1242   }
1243   // create all remaining elements (prisms, hexes, and polyhedra)
1244   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
1245     if (m_Cell2PCell[id_cell] != -1) {
1246       for (int i_face = 0; i_face < m_Part.c2cGSize(id_cell); ++i_face) {
1247         vtkIdType id_neigh = m_Part.c2cGG(id_cell, i_face);
1248         if (id_neigh == -1) {
1249           EG_BUG;
1250         }
1251         bool create_corner_faces = false;
1252         if (m_Cell2PCell[id_neigh] == -1 && !isSurface(id_neigh, m_Grid)) {
1253           create_corner_faces = true;
1254         }
1255         if (create_corner_faces) {
1256           QVector<vtkIdType> face_nodes;
1257           getFaceOfCell(m_Grid, id_cell, i_face, face_nodes);
1258           foreach (vtkIdType id_node, face_nodes) {
1259             if (m_Node2PCell[id_node] == -1) {
1260               EG_BUG;
1261             }
1262             createCornerFace(id_cell, i_face, id_node);
1263           }
1264         } else {
1265           if (id_neigh > id_cell || isSurface(id_neigh, m_Grid)) {
1266             createFaceFace(id_cell, i_face);
1267           }
1268         }
1269       }
1270     }
1271   }
1272
1273   // create all dual cells
1274   if (m_CreateDualMesh) {
1275     for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
1276       QSet<vtkIdType> c1;
1277       for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
1278         vtkIdType id_cell = m_Part.n2cGG(id_node, i);
1279         if (m_Cell2PCell[id_cell] == -1) {
1280           c1.insert(id_cell);
1281         }
1282       }
1283       if (m_Node2PCell[id_node] != -1) {
1284         for (int i_neigh = 0; i_neigh < m_Part.n2nGSize(id_node); ++i_neigh) {
1285           vtkIdType id_neigh = m_Part.n2nGG(id_node, i_neigh);
1286           if (m_Node2PCell[id_neigh] != -1 && id_neigh > id_node) {
1287
1288             // check if any of the adjacent cells (id_node <-> id_neigh) needs to be "dualised"
1289             QSet<vtkIdType> c2;
1290             for (int i = 0; i < m_Part.n2cGSize(id_neigh); ++i) {
1291               vtkIdType id_cell = m_Part.n2cGG(id_neigh, i);
1292               if (m_Cell2PCell[id_cell] == -1) {
1293                 c2.insert(id_cell);
1294               }
1295             }
1296             c2.intersect(c1);
1297             if (c2.size() > 0) {
1298               createEdgeFace(id_node, id_neigh);
1299             }
1300           }
1301         }
1302         QSet<int> bcs;
1303         for (int i_cell = 0; i_cell < m_Part.n2cGSize(id_node); ++i_cell) {
1304           vtkIdType id_cell = m_Part.n2cGG(id_node, i_cell);
1305           if (isSurface(id_cell, m_Grid)) {
1306             bcs.insert(cell_code->GetValue(id_cell));
1307           }
1308         }
1309         foreach (int bc, bcs) {
1310           createPointFace(id_node, bc);
1311         }
1312       }
1313     }
1314   }
1315
1316   computePoints();
1317
1318   //splitConcaveFaces();
1319   sortFaces();
1320   collectBoundaryConditions();
1321 }
1322
1323 vec3_t PolyMesh::faceNormal(int i)
1324 {
1325   int N = m_Faces[i].node.size();
1326   QVector<vec3_t> x(N + 1);
1327   vec3_t xc(0,0,0);
1328   for (int j = 0; j < N; ++j) {
1329     x[j] = m_Points[m_Faces[i].node[j]];
1330     xc += x[j];
1331   }
1332   x[N] = x[0];
1333   xc *= 1.0/N;
1334   vec3_t n(0,0,0);
1335   for (int j = 0; j < N; ++j) {
1336     vec3_t u = x[j] - xc;
1337     vec3_t v = x[j+1] - xc;
1338     n += 0.5*u.cross(v);
1339   }
1340   return n;
1341 }
1342
1343 void PolyMesh::checkFaceOrientation()
1344 {
1345   for (int i = 0; i < m_Faces.size(); ++i) {
1346     int N = m_Faces[i].node.size();
1347     vec3_t n = faceNormal(i);
1348     n.normalise();
1349     if (n*m_Faces[i].ref_vec < 0) {
1350       QVector<int> old_node = m_Faces[i].node;
1351       for (int j = 0; j < N; ++j) {
1352         m_Faces[i].node[j] = old_node[N-1-j];
1353       }
1354     }
1355   }
1356 }
1357
1358 void PolyMesh::buildPoint2Face()
1359 {
1360   m_Point2Face.clear();
1361   m_Point2Face.resize(totalNumNodes());
1362   for (int face = 0; face < numFaces(); ++face) {
1363     foreach (int node, m_Faces[face].node) {
1364       m_Point2Face[node].append(face);
1365     }
1366   }
1367 }
1368
1369 void PolyMesh::buildPCell2Face()
1370 {
1371   m_PCell2Face.clear();
1372   m_PCell2Face.resize(m_NumPolyCells);
1373   for (int i = 0; i < m_Faces.size(); ++i) {
1374     if (owner(i) >= m_NumPolyCells) {
1375       EG_BUG;
1376     }
1377     m_PCell2Face[owner(i)].append(i);
1378     if (neighbour(i) >= 0) {
1379       if (neighbour(i) >= m_NumPolyCells) {
1380         EG_BUG;
1381       }
1382       m_PCell2Face[neighbour(i)].append(i);
1383     }
1384   }
1385 }
1386
1387 void PolyMesh::invertFace(int i)
1388 {
1389   invertQContainer(m_Faces[i].node);
1390   m_Faces[i].ref_vec *= -1;
1391 }
1392
1393 void PolyMesh::sortFaces()
1394 {
1395   // compute hashes
1396   int max_bc = -1;
1397   foreach (face_t face, m_Faces) {
1398     max_bc = max(max_bc, face.bc);
1399   }
1400   int hash_stride = 1;
1401   foreach (face_t face, m_Faces) {
1402     hash_stride = max(hash_stride, face.owner);
1403   }
1404   QVector<QList<face_t> > sort_lists(hash_stride*(max_bc + 1));
1405   foreach (face_t face, m_Faces) {
1406     int i = face.bc*hash_stride + face.owner;
1407     sort_lists[i].append(face);
1408   }
1409   m_Faces.clear();
1410   for (int i = 0; i < sort_lists.size(); ++i) {
1411     qSort(sort_lists[i]);
1412     m_Faces += sort_lists[i];
1413   }
1414 }