src/libengrid/polymesh.cpp

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