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 "engrid.h"
  24 #include "polymolecule.h"
  25 #include "guimainwindow.h"
  26
  27 #include <vtkXMLPolyDataWriter.h>
  28 #include <vtkTriangleFilter.h>
  29 #include <vtkCurvatures.h>
  30
  31 // ==========
  32 //   face_t
  33 // ==========
  34
  35 PolyMesh::face_t::face_t(int N, int o, int n, vec3_t rv, int b)
  36 {
  37   node.resize(N);
  38   owner = o;
  39   neighbour = n;
  40   bc = b;
  41   ref_vec = rv;
  42 }
  43
  44
  45 // ==========
  46 //   node_t
  47 // ==========
  48
  49 PolyMesh::node_t::node_t(const QVector<vtkIdType> &ids)
  50 {
  51   id = ids;
  52   qSort(id);
  53 }
  54
  55 PolyMesh::node_t::node_t(vtkIdType id1)
  56 {
  57   id.resize(1);
  58   id[0] = id1;
  59 }
  60
  61 PolyMesh::node_t::node_t(vtkIdType id1, vtkIdType id2)
  62 {
  63   id.resize(2);
  64   id[0] = id1;
  65   id[1] = id2;
  66   qSort(id);
  67 }
  68
  69 PolyMesh::node_t::node_t(vtkIdType id1, vtkIdType id2, vtkIdType id3)
  70 {
  71   id.resize(3);
  72   id[0] = id1;
  73   id[1] = id2;
  74   id[2] = id3;
  75   qSort(id);
  76 }
  77
  78 PolyMesh::node_t::node_t(vtkIdType id1, vtkIdType id2, vtkIdType id3, vtkIdType id4)
  79 {
  80   id.resize(4);
  81   id[0] = id1;
  82   id[1] = id2;
  83   id[2] = id3;
  84   id[3] = id4;
  85   qSort(id);
  86 }
  87
  88 bool PolyMesh::node_t::operator<(const PolyMesh::node_t &N) const
  89 {
  90   int num = min(id.size(), N.id.size());
  91   for (int i = 0; i < num; ++i) {
  92     if (id[i] < N.id[i]) {
  93       return true;
  94     }
  95     if (id[i] > N.id[i]) {
  96       return false;
  97     }
  98   }
  99   if (id.size() < N.id.size()) {
 100     return true;
 101   }
 102   return false;
 103 }
 104
 105 bool PolyMesh::node_t::operator>(const PolyMesh::node_t &N) const
 106 {
 107   int num = min(id.size(), N.id.size());
 108   for (int i = 0; i < num; ++i) {
 109     if (id[i] > N.id[i]) {
 110       return true;
 111     }
 112     if (id[i] < N.id[i]) {
 113       return false;
 114     }
 115   }
 116   if (id.size() > N.id.size()) {
 117     return true;
 118   }
 119   return false;
 120 }
 121
 122 bool PolyMesh::node_t::operator==(const PolyMesh::node_t &N) const
 123 {
 124   if (id.size() != N.id.size()) {
 125     return false;
 126   }
 127   for (int i = 0; i < id.size(); ++i) {
 128     if (id[i] != N.id[i]) {
 129       return false;
 130     }
 131   }
 132   return true;
 133 }
 134
 135
 136
 137
 138
 139 // ============
 140 //   PolyMesh
 141 // ============
 142
 143 PolyMesh::PolyMesh(vtkUnstructuredGrid *grid, bool dualise, double pull_in, bool optimise, bool split_faces, bool split_cells)
 144 {
 145   m_CreateDualMesh = false;
 146   m_OptimiseConvexity = false;
 147   m_SplitCells = false;
 148   m_SplitFaces = split_faces;
 149   if (dualise) {
 150     for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
 151       if (isVolume(id_cell, grid) && grid->GetCellType(id_cell) != VTK_POLYHEDRON) {
 152         m_CreateDualMesh = true;
 153         break;
 154       }
 155     }
 156     if (m_CreateDualMesh) {
 157       m_OptimiseConvexity = optimise;
 158
 159     }
 160   }
 161
 162   m_AttractorWeight = 0.0;
 163   m_PullInFactor = pull_in;
 164   m_OptimiseConvexity = optimise && m_CreateDualMesh;
 165   m_Grid = grid;
 166   m_Part.setGrid(m_Grid);
 167   m_Part.setAllCells();
 168   findPolyCells();
 169   createNodesAndFaces();
 170   checkFaceOrientation();
 171   buildPoint2Face();
 172   buildPCell2Face();
 173   computePoints();
 174
 175   if (m_OptimiseConvexity) {
 176     for (int iter = 0; iter < 2; ++iter) {
 177       int num_bad = 0;
 178       int i_improve = 0;
 179       for (int i = 0; i < numCells(); ++i) {
 180
 181         // check if any of the faces is a boundary face
 182         bool no_boundary = true;
 183         for (int j = 0; j < numFacesOfPCell(i); ++j) {
 184           int bc = boundaryCode(pcell2Face(i, j));
 185           if (bc != 0) {
 186             no_boundary = false;
 187             break;
 188           }
 189         }
 190
 191         if (no_boundary) {
 192           PolyMolecule pm(this, i);
 193           if (!pm.allPositive()) {
 194             ++i_improve;
 195             pm.optimise();
 196             if (pm.minPyramidVolume() < 0) {
 197               ++num_bad;
 198             }
 199           }
 200         }
 201       }
 202       cout << i_improve << " cells out of " << numCells() << " were optimised." << endl;
 203       if (num_bad == 0) {
 204         break;
 205       }
 206     }
 207   }
 208
 209   if (m_SplitFaces) {
 210     splitConcaveFaces();
 211   }
 212
 213   if (m_SplitCells) {
 214     for (int iter = 0; iter < 0; ++iter) {
 215       int num_improved = 0;
 216       int num_cells = numCells();
 217       for (int i = 0; i < num_cells; ++i) {
 218         PolyMolecule pm(this, i);
 219         if (!pm.allPositive()) {
 220           ++num_improved;
 221           pm.centreSplit();
 222         }
 223       }
 224       buildPoint2Face();
 225       buildPCell2Face();
 226       cout << num_improved << " cells out of " << num_cells << " were split." << endl;
 227       if (num_improved == 0) {
 228         break;
 229       }
 230     }
 231   }
 232   sortFaces();
 233   buildPoint2Face();
 234   buildPCell2Face();
 235 }
 236
 237 void PolyMesh::merge(PolyMesh *poly)
 238 {
 239   foreach (face_t face, poly->m_Faces) {
 240     for (int i = 0; i < face.node.size(); ++i) {
 241       face.node[i] += m_Points.size();
 242     }
 243     face.owner += m_NumPolyCells;
 244     face.neighbour += m_NumPolyCells;
 245     m_Faces.append(face);
 246   }
 247   m_Points += poly->m_Points;
 248   m_NumPolyCells += poly->m_NumPolyCells;
 249   sortFaces();
 250   collectBoundaryConditions();
 251 }
 252
 253 void PolyMesh::triangulateBadFaces()
 254 {
 255   m_IsBadCell.fill(false, numCells());
 256   int num_bad = 0;
 257   for (int i = 0; i < numCells(); ++i) {
 258     PolyMolecule pm(this, i);
 259     if (!pm.allPositive()) {
 260       ++num_bad;
 261       m_IsBadCell[i] = true;
 262     }
 263   }
 264   QList<int> bad_faces;
 265   for (int i = 0; i < numFaces(); ++i) {
 266     if (m_IsBadCell[owner(i)]) {
 267       bad_faces.append(i);
 268     } else if (neighbour(i) != -1) {
 269       if (m_IsBadCell[neighbour(i)]) {
 270         bad_faces.append(i);
 271       }
 272     }
 273   }
 274
 275   foreach (int i_face, bad_faces) {
 276     face_t face = m_Faces[i_face];
 277     QVector<vec3_t> x(face.node.size());
 278     for (int i = 0; i < x.size(); ++i) {
 279       x[i] = nodeVector(face.node[i]);
 280     }
 281     EG_VTKSP(vtkPolyData, poly);
 282     createPolyData(x, poly);
 283     EG_VTKSP(vtkTriangleFilter, tri);
 284     tri->SetInputData(poly);
 285     tri->Update();
 286     if (tri->GetOutput()->GetNumberOfPoints() > face.node.size()) {
 287       EG_BUG;
 288     }
 289     QVector<face_t> new_faces(tri->GetOutput()->GetNumberOfCells(), face);
 290     for (int i = 0; i < new_faces.size(); ++i) {
 291       new_faces[i].node.resize(3);
 292       EG_GET_CELL(i, tri->GetOutput());
 293       for (int j = 0; j < 3; ++j) {
 294         new_faces[i].node[j] = face.node[pts[j]];
 295       }
 296     }
 297     m_Faces[i_face] = new_faces[0];
 298     for (int i = 1; i < new_faces.size(); ++i) {
 299       m_Faces.append(new_faces[i]);
 300     }
 301   }
 302
 303   sortFaces();
 304   buildPoint2Face();
 305   buildPCell2Face();
 306
 307   cout << num_bad << " cells out of " << numCells() << " are concave" << endl;
 308   cout << bad_faces.size() << " faces have been triangulated" << endl;
 309 }
 310
 311 /*
 312 void PolyMesh::splitConcaveFaces()
 313 {
 314   for (int i_cell = 0; i_cell < numCells(); ++i_cell) {
 315     if (m_IsBadCell[i_cell]) {
 316       PolyMolecule pm(this, i_cell);
 317       EG_VTKSP(vtkPolyData, poly);
 318       pm.createPolyData(poly);
 319       EG_VTKSP(vtkCurvatures, curv);
 320       curv->SetCurvatureTypeToMinimum();
 321       curv->SetInput(poly);
 322       curv->Update();
 323       EG_VTKDCN(vtkDoubleArray, curvature, curv->GetOutput(), "Minimum_Curvature");
 324       vtkIdType id_ref = 0;
 325       vec3_t x_ref;
 326       double curv_min = 1e99;
 327       EG_FORALL_NODES(id_node, poly) {
 328         if (curvature->GetValue(id_node) < curv_min) {
 329           curv_min = curvature->GetValue(id_node);
 330           id_ref = id_node;
 331           poly->GetPoint(id_ref, x_ref.data());
 332         }
 333       }
 334       vec3_t n_ref;
 335       EG_FORALL_CELLS(id_face, poly) {
 336         EG_GET_CELL(id_face, poly);
 337         bool found = false;
 338         if (type_cell != VTK_TRIANGLE) {
 339           EG_BUG;
 340         }
 341         for (int i = 0; i < num_pts; ++i) {
 342           if (pts[i] == id_ref) {
 343             found = true;
 344             break;
 345           }
 346         }
 347         if (found) {
 348           QVector<vec3_t> x(num_pts);
 349           for (int i = 0; i < num_pts; ++i) {
 350             poly->GetPoint(pts[i], x[i].data());
 351           }
 352           n_ref = (x[1] - x[0]);
 353           n_ref = n_ref.cross(x[2] - x[0]);
 354           n_ref.normalise();
 355           break;
 356         }
 357       }
 358
 359
 360       EG_VTKSP(vtkDoubleArray, div);
 361       div->SetNumberOfComponents(1);
 362       div->SetNumberOfTuples(poly->GetNumberOfPoints());
 363       div->SetName("div");
 364
 365       bool done = false;
 366       vtkIdType id_node0 = id_ref;
 367       vtkIdType id_node1 = id_ref;
 368       QVector<QSet<vtkIdType> > n2n;
 369       QVector<QSet<vtkIdType> > n2c;
 370       QVector<QVector<vtkIdType> > c2c;
 371       createPolyDataN2N(poly, n2n);
 372       createPolyDataN2C(poly, n2c);
 373       createPolyDataC2C(poly, c2c);
 374
 375       vec3_t n_node(0,0,0);
 376       foreach (vtkIdType id_face, n2c[id_ref]) {
 377         EG_GET_CELL(id_face, poly);
 378         QVector<vec3_t> x(3);
 379         for (int i = 0; i < 3; ++i) {
 380           poly->GetPoint(pts[i], x[i].data());
 381         }
 382         n_node += triNormal(x[0], x[1], x[2]);
 383       }
 384       n_node.normalise();
 385
 386       QVector<bool> div_node(poly->GetNumberOfPoints(), false);
 387       EG_FORALL_NODES(id_node, poly) {
 388         div->SetValue(id_node, 0);
 389       }
 390       div_node[id_ref] = true;
 391       double dv = 1.0;
 392       div->SetValue(id_ref, dv);
 393       dv += 0.2;
 394       int count = 0;
 395       double c_min = 1e99;
 396       double L = 1e99;
 397       foreach (vtkIdType id, n2n[id_node0]) {
 398         vec3_t x;
 399         poly->GetPoint(id, x.data());
 400         L = min(L, (x - x_ref).abs());
 401         double c = curvature->GetValue(id);
 402         if (c < c_min) {
 403           c_min = c;
 404           id_node1 = id;
 405         }
 406       }
 407       div_node[id_node1] = true;
 408       div->SetValue(id_node1, dv);
 409       dv += 0.2;
 410       while (!done) {
 411         vtkIdType id_node2 = id_node1;
 412         double d_min = 1e99;
 413         foreach (vtkIdType id, n2n[id_node1]) {
 414           if (id == id_ref && id_node0 != id_ref) {
 415             done = true;
 416             break;
 417           }
 418         }
 419         if (!done) {
 420           foreach (vtkIdType id, n2n[id_node1]) {
 421             if (!div_node[id]) {
 422               vec3_t x;
 423               poly->GetPoint(id, x.data());
 424               if ((x - x_ref)*n_node < 0.25*L) {
 425                 double d = (x - x_ref)*n_ref;
 426                 if (fabs(d) < d_min && d < 0.1*L) {
 427                   d_min = d;
 428                   id_node2 = id;
 429                 }
 430               }
 431             }
 432           }
 433           div_node[id_node2] = true;
 434           div->SetValue(id_node2, dv);
 435           dv += 0.2;
 436           id_node0 = id_node1;
 437           id_node1 = id_node2;
 438           ++count;
 439           if (count == 100) {
 440             EG_BUG;
 441           }
 442         }
 443       }
 444       poly->GetPointData()->AddArray(div);
 445       EG_VTKSP(vtkXMLPolyDataWriter, vtp);
 446       vtp->SetFileName(qPrintable(GuiMainWindow::pointer()->getCwd() + "/triangulated_cell.vtp"));
 447       vtp->SetInput(poly);
 448       vtp->SetDataModeToAscii();
 449       vtp->Write();
 450       EG_BUG;
 451     }
 452   }
 453 }
 454 */
 455
 456 void PolyMesh::getFacesOfEdgeInsideCell(vtkIdType id_cell, vtkIdType id_node1, vtkIdType id_node2, int &face1, int &face2)
 457 {
 458   EG_GET_CELL(id_cell, m_Grid);
 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   EG_GET_CELL(id_cell, m_Grid);
 716   for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
 717     vtkIdType id_neigh = m_Part.n2nGG(id_node, i);
 718     if (face_nodes.contains(id_neigh)) {
 719       edge_nodes.append(id_neigh);
 720     }
 721   }
 722   if (edge_nodes.size() != 2) {
 723     EG_BUG;
 724   }
 725   int owner = m_Cell2PCell[id_cell];
 726   if (owner == -1) {
 727     EG_BUG;
 728   }
 729   int neighbour = m_Node2PCell[id_node];
 730   int bc = 0;
 731   if (neighbour == -1) {
 732     EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
 733     vtkIdType id_face = m_Part.c2cGG(id_cell, i_face);
 734     if (id_face == -1) {
 735       EG_BUG;
 736     }
 737     if (!isSurface(id_cell, m_Grid)) {
 738       EG_BUG;
 739     }
 740     bc = cell_code->GetValue(id_face);
 741   }
 742   QList<node_t> nodes;
 743   nodes.append(node_t(id_node));
 744   nodes.append(node_t(id_node, edge_nodes[0]));
 745   nodes.append(node_t(face_nodes));
 746   nodes.append(node_t(id_node, edge_nodes[1]));
 747   vec3_t n = getNormalOfCell(m_Grid, id_cell, i_face);
 748   n.normalise();
 749   createFace(nodes, owner, neighbour, n, bc);
 750 }
 751
 752 void PolyMesh::createEdgeFace(vtkIdType id_node1, vtkIdType id_node2)
 753 {
 754   // check if additional edge node needs to be created
 755   // (transition from boundary layer to far-field)
 756   // try to find a boundary layer cell which contains both nodes
 757   bool add_edge_node = false;
 758   for (int i = 0; i < m_Part.n2cGSize(id_node1); ++i) {
 759     vtkIdType id_cell = m_Part.n2cGG(id_node1, i);
 760     if (m_Cell2PCell[id_cell] != -1) {
 761       if (cellContainsNode(m_Grid, id_cell, id_node1) && cellContainsNode(m_Grid, id_cell, id_node2)) {
 762         add_edge_node = true;
 763         break;
 764       }
 765     }
 766   }
 767   if (!add_edge_node) {
 768     for (int i = 0; i < m_Part.n2cGSize(id_node2); ++i) {
 769       vtkIdType id_cell = m_Part.n2cGG(id_node2, i);
 770       if (m_Cell2PCell[id_cell] != -1) {
 771         if (cellContainsNode(m_Grid, id_cell, id_node1) && cellContainsNode(m_Grid, id_cell, id_node2)) {
 772           add_edge_node = true;
 773           break;
 774         }
 775       }
 776     }
 777   }
 778
 779   QList<vtkIdType> cells;
 780   bool loop = false;
 781   getSortedEdgeCells(id_node1, id_node2, cells, loop);
 782   int owner     = m_Node2PCell[id_node1];
 783   int neighbour = m_Node2PCell[id_node2];
 784   if (owner == -1) {
 785     EG_BUG;
 786   }
 787   if (neighbour == -1) {
 788     EG_BUG;
 789   }
 790   if (owner > neighbour) {
 791     swap(id_node1, id_node2);
 792     swap(owner, neighbour);
 793   }
 794   QList<node_t> nodes;
 795   for (int i_cells = 0; i_cells < cells.size(); ++i_cells) {
 796     node_t node;
 797     QVector<vtkIdType> pts1;
 798     getPointsOfCell(m_Grid, cells[i_cells], pts1);
 799     if (m_Cell2PCell[cells[i_cells]] != -1) {
 800       int i_cells2 = 0;
 801       if (i_cells == 0) {
 802         i_cells2 = 1;
 803       } else if (i_cells == cells.size() - 1) {
 804         i_cells2 = cells.size() - 2;
 805       } else {
 806         EG_BUG;
 807       }
 808       QVector<vtkIdType> pts2;
 809       getPointsOfCell(m_Grid, cells[i_cells2], pts2);
 810       QSet<vtkIdType> p1, p2;
 811       for (int i_pts1 = 0; i_pts1 < pts1.size(); ++i_pts1) {
 812         p1.insert(pts1[i_pts1]);
 813       }
 814       for (int i_pts2 = 0; i_pts2 < pts2.size(); ++i_pts2) {
 815         p2.insert(pts2[i_pts2]);
 816       }
 817       QSet<vtkIdType> face_nodes = p1.intersect(p2);
 818       node.id.resize(face_nodes.size());
 819       int i = 0;
 820       foreach (vtkIdType id_node, face_nodes) {
 821         node.id[i] = id_node;
 822         ++i;
 823       }
 824     } else {
 825       node.id.resize(pts1.size());
 826       for (int i_pts1 = 0; i_pts1 < pts1.size(); ++i_pts1) {
 827         node.id[i_pts1] = pts1[i_pts1];
 828       }
 829     }
 830     qSort(node.id);
 831     nodes.append(node);
 832   }
 833   if (!loop) {
 834     EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
 835     if (cell_code->GetValue(cells.first()) != cell_code->GetValue(cells.last()) || add_edge_node) {
 836       nodes.append(node_t(id_node1, id_node2));
 837     }
 838   }
 839   vec3_t x1, x2;
 840   m_Grid->GetPoint(id_node1, x1.data());
 841   m_Grid->GetPoint(id_node2, x2.data());
 842   createFace(nodes, owner, neighbour, x2 - x1, 0);
 843 }
 844
 845 void PolyMesh::createFaceFace(vtkIdType id_cell, int i_face)
 846 {
 847   if (m_Cell2PCell[id_cell] == -1) {
 848     EG_BUG;
 849   }
 850   int owner = m_Cell2PCell[id_cell];
 851   int neighbour = -1;
 852   vtkIdType id_neigh = m_Part.c2cGG(id_cell, i_face);
 853   int bc = 0;
 854   if (id_neigh != -1) {
 855     if (isVolume(id_neigh, m_Grid)) {
 856       if (m_Cell2PCell[id_neigh] == -1) {
 857         EG_BUG;
 858       }
 859       neighbour = m_Cell2PCell[id_neigh];
 860     } else {
 861       EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
 862       bc = cell_code->GetValue(id_neigh);
 863     }
 864   }
 865   QVector<vtkIdType> tmp_node_ids;
 866   getFaceOfCell(m_Grid, id_cell, i_face, tmp_node_ids);
 867   vec3_t n = getNormalOfCell(m_Grid, id_cell, i_face);
 868   n.normalise();
 869   QVector<vtkIdType> node_ids(tmp_node_ids.size() + 1);
 870   for (int i = 0; i < tmp_node_ids.size(); ++i) {
 871     node_ids[i] = tmp_node_ids[i];
 872   }
 873   node_ids[node_ids.size() - 1] = node_ids[0];
 874   QList<node_t> nodes;
 875   for (int i = 0; i < node_ids.size() - 1; ++i) {
 876     nodes.append(node_t(node_ids[i]));
 877     if (m_Node2PCell[node_ids[i]] != -1 && m_Node2PCell[node_ids[i+1]] != -1) {
 878       nodes.append(node_t(node_ids[i], node_ids[i+1]));
 879     }
 880   }
 881   createFace(nodes, owner, neighbour, n, bc);
 882 }
 883
 884 void PolyMesh::createPointFace(vtkIdType id_node, int bc)
 885 {
 886   bool is_loop;
 887   QList<vtkIdType> faces;
 888   getSortedPointFaces(id_node, bc, faces, is_loop);
 889   if (faces.size() == 0) {
 890     return;
 891   }
 892   QList<node_t> nodes;
 893   vec3_t n(0,0,0);
 894   if (faces.size() == 0) {
 895     EG_BUG;
 896   }
 897   foreach (vtkIdType id_face, faces) {
 898     node_t node;
 899     EG_GET_CELL(id_face, m_Grid);
 900     node.id.resize(num_pts);
 901     for (int i_pts = 0; i_pts < num_pts; ++i_pts) {
 902       node.id[i_pts] = pts[i_pts];
 903     }
 904     qSort(node.id);
 905     nodes.append(node);
 906     n += GeometryTools::cellNormal(m_Grid, id_face);
 907   }
 908   n.normalise();
 909   if (!is_loop) {
 910     bool prepend = false;
 911     vtkIdType id_face = faces.last();
 912     while (id_face != -1) {
 913       bool found = false;
 914       EG_GET_CELL(id_face, m_Grid);
 915       QList<vtkIdType> id_neigh_node;
 916       vtkIdType id_neigh = -1;
 917       for (int i = 0; i < m_Part.c2cGSize(id_face); ++i) {
 918         id_neigh = m_Part.c2cGG(id_face, i);
 919         if (id_neigh == -1) {
 920           EG_BUG;
 921         }
 922         if (!isSurface(id_neigh, m_Grid)) {
 923           EG_BUG;
 924         }
 925         if (cellContainsNode(m_Grid, id_neigh, id_node)) {
 926           if (!faces.contains(id_neigh)) {
 927             if (found) {
 928               EG_BUG;
 929             }
 930             for (int j = 0; j < num_pts; ++j) {
 931               if (pts[j] != id_node && cellContainsNode(m_Grid, id_neigh, pts[j])) {
 932                 id_neigh_node.append(pts[j]);
 933               }
 934             }
 935           }
 936         }
 937       }
 938
 939       if (id_neigh_node.size() == 0) {
 940         EG_BUG;
 941       }
 942
 943       if (id_neigh_node.size() == 1) {
 944         if (prepend) {
 945           nodes.prepend(node_t(id_node, id_neigh_node[0]));
 946           id_face = -1;
 947         } else {
 948           nodes.append(node_t(id_node, id_neigh_node[0]));
 949           prepend = true;
 950           id_face = faces.first();
 951         }
 952       } else {
 953         if (id_neigh_node.size() > 2) {
 954           EG_BUG;
 955         }
 956         if (faces.size() != 1) {
 957           EG_BUG;
 958         }
 959         nodes.prepend(node_t(id_node, id_neigh_node[0]));
 960         nodes.append(node_t(id_node, id_neigh_node[1]));
 961         id_face = -1;
 962       }
 963     }
 964
 965     // check if this is a transition node (boundary layer -> far-field)
 966     bool dual_cell_found = false;
 967     bool cell_cell_found = false;
 968     for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
 969       if (m_Cell2PCell[m_Part.n2cGG(id_node,i)] == -1) {
 970         cell_cell_found = true;
 971       } else {
 972         dual_cell_found = true;
 973       }
 974     }
 975
 976     if (m_Part.n2bcGSize(id_node) > 2 || (dual_cell_found && cell_cell_found)) {
 977       nodes.prepend(node_t(id_node));
 978     }
 979   }
 980   int owner     = m_Node2PCell[id_node];
 981   int neighbour = -1;
 982   createFace(nodes, owner, neighbour, n, bc);
 983 }
 984
 985 void PolyMesh::computePoints()
 986 {
 987   QVector<node_t> nodes;
 988   m_Nodes.getQVector(nodes);
 989   m_Points.resize(nodes.size());
 990   m_PointWeights.fill(1.0, m_Grid->GetNumberOfPoints());
 991
 992   // find transition nodes
 993   QVector<bool> is_transition_node(m_Grid->GetNumberOfPoints(), false);
 994   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 995     if (m_Node2PCell[id_node] != -1) {
 996       for (int i_cell = 0; i_cell < m_Part.n2cGSize(id_node); ++i_cell) {
 997         vtkIdType id_cell = m_Part.n2cGG(id_node, i_cell);
 998         if (m_Cell2PCell[id_cell] != -1) {
 999           is_transition_node[id_node] = true;
1000           break;
1001         }
1002       }
1003     }
1004   }
1005
1006   // compute weights (attraction for convex corners)
1007   // mark convex edge nodes
1008   QVector<bool> is_convex_node(m_Grid->GetNumberOfPoints(), false);
1009   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
1010     if (m_Cell2PCell[id_cell] != -1) {
1011       vec3_t xc1 = cellCentre(m_Grid, id_cell);
1012       int N = 0;
1013       vec3_t n_out(0,0,0);
1014       for (int i_face = 0; i_face < m_Part.c2cGSize(id_cell); ++i_face) {
1015         vtkIdType id_neigh = m_Part.c2cGG(id_cell, i_face);
1016         if (id_neigh != -1) {
1017           if (m_Cell2PCell[id_neigh] == -1 && !isSurface(id_neigh, m_Grid)) {
1018             ++N;
1019             n_out += getNormalOfCell(m_Grid, id_cell, i_face);
1020           }
1021         }
1022       }
1023       if (N > 0) {
1024         n_out.normalise();
1025         for (int i_face = 0; i_face < m_Part.c2cGSize(id_cell); ++i_face) {
1026           vtkIdType id_neigh = m_Part.c2cGG(id_cell, i_face);
1027           if (id_neigh != -1) {
1028             if (m_Cell2PCell[id_neigh] != -1) {
1029               vec3_t xc2 = cellCentre(m_Grid, id_neigh);
1030               QVector<vtkIdType> face_nodes;
1031               getFaceOfCell(m_Grid, id_cell, i_face, face_nodes);
1032               if (face_nodes.size() == 0) {
1033                 EG_BUG;
1034               }
1035               vec3_t xf(0,0,0);
1036               foreach (vtkIdType id_node, face_nodes) {
1037                 vec3_t x;
1038                 m_Grid->GetPoint(id_node, x.data());
1039                 xf += x;
1040               }
1041               xf *= 1.0/face_nodes.size();
1042               vec3_t v1 = xc1 - xf;
1043               vec3_t v2 = xc2 - xf;
1044               v1.normalise();
1045               v2.normalise();
1046               if ((v1+v2)*n_out < 0) {
1047                 double w = 1 + m_AttractorWeight*(1.0 + v1*v2);
1048                 foreach (vtkIdType id_node, face_nodes) {
1049                   if (m_Node2PCell[id_node] != -1) {
1050                     m_PointWeights[id_node] = max(m_PointWeights[id_node], w);
1051                     if (v1*v2 > 0) {
1052                       //is_convex_node[id_node] = true;
1053                     }
1054                   }
1055                 }
1056               }
1057             }
1058           }
1059         }
1060       }
1061     }
1062   }
1063
1064   // mapping for simple nodes
1065   QVector<int> prime2dual(m_Grid->GetNumberOfPoints(), -1);
1066
1067   // compute the point locations
1068   for (int i = 0; i < nodes.size(); ++i) {
1069     if (nodes[i].id.size() == 0) {
1070       EG_BUG;
1071     }
1072     m_Points[i] = vec3_t(0,0,0);
1073     double weight_sum = 0.0;
1074     foreach (vtkIdType id, nodes[i].id) {
1075       vec3_t x;
1076       m_Grid->GetPoint(id, x.data());
1077       weight_sum += m_PointWeights[id];
1078       m_Points[i] += m_PointWeights[id]*x;
1079     }
1080     //m_Points[i] *= 1.0/nodes[i].id.size();
1081     m_Points[i] *= 1.0/weight_sum;
1082     if (nodes[i].id.size() == 1) {
1083       prime2dual[nodes[i].id[0]] = i;
1084     }
1085   }
1086
1087   // correct transition nodes (pull-in)
1088   for (vtkIdType id_node1 = 0; id_node1 < m_Grid->GetNumberOfPoints(); ++id_node1) {
1089     if (is_transition_node[id_node1]) {
1090       if (prime2dual[id_node1] == -1) {
1091         EG_BUG;
1092       }
1093       vtkIdType id_node2 = -1;
1094       bool pull_in = false;
1095       for (int i_neigh = 0; i_neigh < m_Part.n2nGSize(id_node1); ++i_neigh) {
1096         vtkIdType id_neigh = m_Part.n2nGG(id_node1, i_neigh);
1097         if (m_Node2PCell[id_neigh] == -1 && !is_transition_node[id_neigh]) {
1098           if (id_node2 == -1) {
1099             pull_in = true;
1100             id_node2 = id_neigh;
1101           } else {
1102             pull_in = false;
1103           }
1104         }
1105       }
1106       if (pull_in) {
1107         double w = m_PullInFactor;
1108         m_Points[prime2dual[id_node1]] = w*m_Points[prime2dual[id_node2]] + (1-w)*m_Points[prime2dual[id_node1]];
1109       }
1110     }
1111   }
1112
1113 }
1114
1115 void PolyMesh::splitConcaveFaces()
1116 {
1117   QList<int> delete_faces;
1118   do {
1119     delete_faces.clear();
1120     QList<face_t> new_faces;
1121     for (int i_face = 0; i_face < m_Faces.size(); ++i_face) {
1122       face_t face = m_Faces[i_face];
1123       int num_nodes = face.node.size();
1124       if (num_nodes >= 4) {
1125         QVector<vec3_t> x_face(num_nodes);
1126         for (int i = 0; i < num_nodes; ++i) {
1127           x_face[i] = nodeVector(face.node[i]);
1128         }
1129         vec3_t xc, n;
1130         GeometryTools::planeFit(x_face, xc, n);
1131         QVector<vec3_t> x(num_nodes + 2);
1132         for (int i = 0; i < num_nodes; ++i) {
1133           x[i+1] = x_face[i];
1134         }
1135         x.first() = x_face.last();
1136         x.last() = x_face.first();
1137         double L_max = 0.1;
1138         int i1 = -1;
1139         vec3_t v;
1140         for (int i = 1; i <= num_nodes; ++i) {
1141           vec3_t x_neigh = 0.5*(x[i-1] + x[i+1]);
1142           double Lnorm = (x[i-1] - x[i+1]).abs();
1143           double L = ((x_neigh - xc).abs() - (x[i] - xc).abs())/Lnorm;
1144           if (L > L_max) {
1145             i1 = i;
1146             L_max = L;
1147             v = x[i] - x_neigh;
1148           }
1149         }
1150         if (i1 != -1) {
1151           int i2 = -1;
1152           double alpha_min = 1e99;
1153           for (int i = 1; i <= num_nodes; ++i) {
1154             if ((i < i1 - 1 || i > i1 + 1) && (i < i1 + num_nodes -1 || i > i1 + num_nodes + 1)) {
1155               double alpha = GeometryTools::angle(x[i] - x[i1], v);
1156               if (alpha < alpha_min) {
1157                 i2 = i;
1158                 alpha_min = alpha;
1159               }
1160             }
1161           }
1162           if (i2 != -1) {
1163             --i1;
1164             --i2;
1165             if (i1 > i2) {
1166               swap(i1, i2);
1167             }
1168             delete_faces.append(i_face);
1169             QVector<int> nodes1;
1170             {
1171               int i = 0;
1172               while (i < num_nodes) {
1173                 nodes1.append(face.node[i]);
1174                 if (i == i1) {
1175                   i = i2;
1176                 } else {
1177                   ++i;
1178                 }
1179               }
1180             }
1181             QVector<int> nodes2;
1182             for (int i = i1; i <= i2; ++i) {
1183               nodes2.append(face.node[i]);
1184             }
1185             face_t face1(nodes1.size(), face.owner, face.neighbour, face.ref_vec, face.bc);
1186             face_t face2(nodes2.size(), face.owner, face.neighbour, face.ref_vec, face.bc);
1187             face1.node = nodes1;
1188             face2.node = nodes2;
1189             new_faces.append(face1);
1190             new_faces.append(face2);
1191           }
1192         }
1193       }
1194     }
1195     {
1196       int offset = 0;
1197       foreach (int i_face, delete_faces) {
1198         m_Faces.removeAt(i_face - offset);
1199         ++offset;
1200       }
1201     }
1202     m_Faces.append(new_faces);
1203     cout << delete_faces.size() << " concave faces have been split." << endl;
1204     delete_faces.clear();
1205   } while (delete_faces.size() > 0);
1206   sortFaces();
1207   buildPoint2Face();
1208   buildPCell2Face();
1209 }
1210
1211 void PolyMesh::collectBoundaryConditions()
1212 {
1213   QSet<int> bcs;
1214   foreach (face_t face, m_Faces) {
1215     if (face.bc != 0) {
1216       bcs.insert(face.bc);
1217     }
1218   }
1219   m_BCs.resize(bcs.size());
1220   qCopy(bcs.begin(), bcs.end(), m_BCs.begin());
1221 }
1222
1223 void PolyMesh::createNodesAndFaces()
1224 {
1225   m_Nodes.resize(m_Grid->GetNumberOfPoints());
1226   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
1227
1228   int num_pcells = 0;
1229   int num_dcells = 0;
1230   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
1231     if (m_Cell2PCell[id_cell] != -1) {
1232       ++num_pcells;
1233     }
1234   }
1235   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
1236     if (m_Node2PCell[id_node] != -1) {
1237       ++num_dcells;
1238     }
1239   }
1240   // create all remaining elements (prisms, hexes, and polyhedra)
1241   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
1242     if (m_Cell2PCell[id_cell] != -1) {
1243       for (int i_face = 0; i_face < m_Part.c2cGSize(id_cell); ++i_face) {
1244         vtkIdType id_neigh = m_Part.c2cGG(id_cell, i_face);
1245         if (id_neigh == -1) {
1246           EG_BUG;
1247         }
1248         bool create_corner_faces = false;
1249         if (m_Cell2PCell[id_neigh] == -1 && !isSurface(id_neigh, m_Grid)) {
1250           create_corner_faces = true;
1251         }
1252         if (create_corner_faces) {
1253           QVector<vtkIdType> face_nodes;
1254           getFaceOfCell(m_Grid, id_cell, i_face, face_nodes);
1255           foreach (vtkIdType id_node, face_nodes) {
1256             if (m_Node2PCell[id_node] == -1) {
1257               EG_BUG;
1258             }
1259             createCornerFace(id_cell, i_face, id_node);
1260           }
1261         } else {
1262           if (id_neigh > id_cell || isSurface(id_neigh, m_Grid)) {
1263             createFaceFace(id_cell, i_face);
1264           }
1265         }
1266       }
1267     }
1268   }
1269
1270   // create all dual cells
1271   if (m_CreateDualMesh) {
1272     for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
1273       QSet<vtkIdType> c1;
1274       for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
1275         vtkIdType id_cell = m_Part.n2cGG(id_node, i);
1276         if (m_Cell2PCell[id_cell] == -1) {
1277           c1.insert(id_cell);
1278         }
1279       }
1280       if (m_Node2PCell[id_node] != -1) {
1281         for (int i_neigh = 0; i_neigh < m_Part.n2nGSize(id_node); ++i_neigh) {
1282           vtkIdType id_neigh = m_Part.n2nGG(id_node, i_neigh);
1283           if (m_Node2PCell[id_neigh] != -1 && id_neigh > id_node) {
1284
1285             // check if any of the adjacent cells (id_node <-> id_neigh) needs to be "dualised"
1286             QSet<vtkIdType> c2;
1287             for (int i = 0; i < m_Part.n2cGSize(id_neigh); ++i) {
1288               vtkIdType id_cell = m_Part.n2cGG(id_neigh, i);
1289               if (m_Cell2PCell[id_cell] == -1) {
1290                 c2.insert(id_cell);
1291               }
1292             }
1293             c2.intersect(c1);
1294             if (c2.size() > 0) {
1295               createEdgeFace(id_node, id_neigh);
1296             }
1297           }
1298         }
1299         QSet<int> bcs;
1300         for (int i_cell = 0; i_cell < m_Part.n2cGSize(id_node); ++i_cell) {
1301           vtkIdType id_cell = m_Part.n2cGG(id_node, i_cell);
1302           if (isSurface(id_cell, m_Grid)) {
1303             bcs.insert(cell_code->GetValue(id_cell));
1304           }
1305         }
1306         foreach (int bc, bcs) {
1307           createPointFace(id_node, bc);
1308         }
1309       }
1310     }
1311   }
1312
1313   computePoints();
1314
1315   //splitConcaveFaces();
1316   sortFaces();
1317   collectBoundaryConditions();
1318 }
1319
1320 vec3_t PolyMesh::faceNormal(int i)
1321 {
1322   int N = m_Faces[i].node.size();
1323   QVector<vec3_t> x(N + 1);
1324   vec3_t xc(0,0,0);
1325   for (int j = 0; j < N; ++j) {
1326     x[j] = m_Points[m_Faces[i].node[j]];
1327     xc += x[j];
1328   }
1329   x[N] = x[0];
1330   xc *= 1.0/N;
1331   vec3_t n(0,0,0);
1332   for (int j = 0; j < N; ++j) {
1333     vec3_t u = x[j] - xc;
1334     vec3_t v = x[j+1] - xc;
1335     n += 0.5*u.cross(v);
1336   }
1337   return n;
1338 }
1339
1340 void PolyMesh::checkFaceOrientation()
1341 {
1342   for (int i = 0; i < m_Faces.size(); ++i) {
1343     int N = m_Faces[i].node.size();
1344     vec3_t n = faceNormal(i);
1345     n.normalise();
1346     if (n*m_Faces[i].ref_vec < 0) {
1347       QVector<int> old_node = m_Faces[i].node;
1348       for (int j = 0; j < N; ++j) {
1349         m_Faces[i].node[j] = old_node[N-1-j];
1350       }
1351     }
1352   }
1353 }
1354
1355 void PolyMesh::buildPoint2Face()
1356 {
1357   m_Point2Face.clear();
1358   m_Point2Face.resize(totalNumNodes());
1359   for (int face = 0; face < numFaces(); ++face) {
1360     foreach (int node, m_Faces[face].node) {
1361       m_Point2Face[node].append(face);
1362     }
1363   }
1364 }
1365
1366 void PolyMesh::buildPCell2Face()
1367 {
1368   m_PCell2Face.clear();
1369   m_PCell2Face.resize(m_NumPolyCells);
1370   for (int i = 0; i < m_Faces.size(); ++i) {
1371     if (owner(i) >= m_NumPolyCells) {
1372       EG_BUG;
1373     }
1374     m_PCell2Face[owner(i)].append(i);
1375     if (neighbour(i) >= 0) {
1376       if (neighbour(i) >= m_NumPolyCells) {
1377         EG_BUG;
1378       }
1379       m_PCell2Face[neighbour(i)].append(i);
1380     }
1381   }
1382 }
1383
1384 void PolyMesh::invertFace(int i)
1385 {
1386   invertQContainer(m_Faces[i].node);
1387   m_Faces[i].ref_vec *= -1;
1388 }
1389
1390 void PolyMesh::sortFaces()
1391 {
1392   // compute hashes
1393   int max_bc = -1;
1394   foreach (face_t face, m_Faces) {
1395     max_bc = max(max_bc, face.bc);
1396   }
1397   if (max_bc < 0) {
1398     EG_ERR_RETURN("mesh is corrupted");
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) + 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 }