src/libengrid/meshpartition.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 "meshpartition.h"
  23 #include "guimainwindow.h"
  24
  25 #include <vtkKdTreePointLocator.h>
  26 #include <vtkClipDataSet.h>
  27 #include <vtkTriangleFilter.h>
  28 #include <vtkGeometryFilter.h>
  29 #include <vtkSTLWriter.h>
  30
  31 #include <vtkEgNormalExtrusion.h>
  32
  33 MeshPartition::MeshPartition()
  34 {
  35   m_Grid = NULL;
  36   m_TrackGrid = false;
  37   resetTimeStamps();
  38 }
  39
  40 MeshPartition::MeshPartition(vtkUnstructuredGrid *grid, bool use_all_cells)
  41 {
  42   setGrid(grid, use_all_cells);
  43 }
  44
  45 MeshPartition::MeshPartition(QString volume_name)
  46 {
  47   m_TrackGrid = false;
  48   resetTimeStamps();
  49   setVolume(volume_name);
  50 }
  51
  52 void MeshPartition::setGrid(vtkUnstructuredGrid *grid, bool use_all_cells)
  53 {
  54   m_TrackGrid = false;
  55   resetTimeStamps();
  56   m_Grid = grid;
  57   if (use_all_cells) {
  58     QVector<vtkIdType> cls(grid->GetNumberOfCells());
  59     for (vtkIdType id_cell = 0; id_cell < cls.size(); ++id_cell) {
  60       cls[id_cell] = id_cell;
  61     }
  62     setCells(cls);
  63   }
  64 }
  65
  66
  67 void MeshPartition::resetTimeStamps()
  68 {
  69   m_CellsStamp = 0;
  70   m_LCellsStamp = 0;
  71   m_NodesStamp = 0;
  72   m_LNodesStamp = 0;
  73   m_N2NStamp = 0;
  74   m_N2CStamp = 0;
  75   m_N2BCStamp = 0;
  76   m_C2CStamp = 0;
  77 }
  78
  79 void MeshPartition::trackGrid(vtkUnstructuredGrid *grid)
  80 {
  81   setGrid(grid);
  82   setAllCells();
  83   m_GridMTime = m_Grid->GetMTime();
  84   m_TrackGrid = true;
  85 }
  86
  87 void MeshPartition::setVolume(QString volume_name)
  88 {
  89   m_Grid = GuiMainWindow::pointer()->getGrid();
  90   resetOrientation(m_Grid);
  91   VolumeDefinition V = GuiMainWindow::pointer()->getVol(volume_name);
  92   QList<vtkIdType> cls;
  93   EG_VTKDCC(vtkIntArray, cell_code,   m_Grid, "cell_code");
  94   EG_VTKDCC(vtkIntArray, cell_orgdir, m_Grid, "cell_orgdir");
  95   EG_VTKDCC(vtkIntArray, cell_curdir, m_Grid, "cell_curdir");
  96   EG_VTKDCC(vtkIntArray, cell_voldir, m_Grid, "cell_voldir");
  97   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
  98     if (isSurface(id_cell, m_Grid)) {
  99       int bc = cell_code->GetValue(id_cell);
 100       cell_voldir->SetValue(id_cell, 0);
 101       if (V.getSign(bc) != 0) {
 102         cls.append(id_cell);
 103         if (V.getSign(bc) == -1) {
 104           cell_voldir->SetValue(id_cell, 1);
 105         }
 106       }
 107     } else {
 108       if (cell_code->GetValue(id_cell) == V.getVC()) {
 109         cls.append(id_cell);
 110       }
 111     }
 112   }
 113   setCells(cls);
 114 }
 115
 116 void MeshPartition::setVolumeOrientation()
 117 {
 118   EG_VTKDCC(vtkIntArray, cell_curdir, m_Grid, "cell_curdir");
 119   EG_VTKDCC(vtkIntArray, cell_voldir, m_Grid, "cell_voldir");
 120   foreach (vtkIdType id_cell, m_Cells) {
 121     if (isSurface(id_cell, m_Grid)) {
 122       if (cell_curdir->GetValue(id_cell) != cell_voldir->GetValue(id_cell)) {
 123         reorientateFace(m_Grid, id_cell);
 124       }
 125     }
 126   }
 127 }
 128
 129 void MeshPartition::setOriginalOrientation()
 130 {
 131   EG_VTKDCC(vtkIntArray, cell_curdir, m_Grid, "cell_curdir");
 132   EG_VTKDCC(vtkIntArray, cell_orgdir, m_Grid, "cell_orgdir");
 133   foreach (vtkIdType id_cell, m_Cells) {
 134     if (isSurface(id_cell, m_Grid)) {
 135       if (cell_curdir->GetValue(id_cell) != cell_orgdir->GetValue(id_cell)) {
 136         reorientateFace(m_Grid, id_cell);
 137       }
 138     }
 139   }
 140 }
 141
 142 void MeshPartition::setRemainder(const MeshPartition& part)
 143 {
 144   setGrid(part.getGrid());
 145   QVector<vtkIdType> rcells;
 146   getRestCells(m_Grid, part.m_Cells, rcells);
 147   setCells(rcells);
 148 }
 149
 150 void MeshPartition::extractToVtkGrid(vtkUnstructuredGrid *new_grid)
 151 {
 152   checkLNodes();
 153   allocateGrid(new_grid, m_Cells.size(), m_Nodes.size());
 154   for (int i_nodes = 0; i_nodes < m_Nodes.size(); ++i_nodes) {
 155     vec3_t x;
 156     m_Grid->GetPoints()->GetPoint(m_Nodes[i_nodes], x.data());
 157     new_grid->GetPoints()->SetPoint(i_nodes, x.data());
 158     copyNodeData(m_Grid, m_Nodes[i_nodes], new_grid, i_nodes);
 159   }
 160   foreach (vtkIdType id_cell, m_Cells) {
 161     /*
 162     vtkIdType N_pts, *pts;
 163     vtkIdType type_cell = m_Grid->GetCellType(id_cell);
 164     m_Grid->GetCellPoints(id_cell, N_pts, pts);
 165     QVector<vtkIdType> new_pts(N_pts);
 166     for (int i = 0; i < N_pts; ++i) {
 167       new_pts[i] = m_LNodes[pts[i]];
 168     }
 169     // update for polyhedral cells here
 170     vtkIdType id_new_cell = new_grid->InsertNextCell(type_cell, N_pts, new_pts.data());
 171     */
 172     vtkIdType id_new_cell = copyCell(m_Grid, id_cell, new_grid, m_LNodes);
 173     copyCellData(m_Grid, id_cell, new_grid, id_new_cell);
 174   }
 175 }
 176
 177 void MeshPartition::addPartition(const MeshPartition& part, double tol)
 178 {
 179   if (m_Grid == part.m_Grid) {
 180     QVector<bool> cell_marked(m_Grid->GetNumberOfCells(), false);
 181     foreach (vtkIdType id_cell, m_Cells) {
 182       cell_marked[id_cell] = true;
 183     }
 184     foreach (vtkIdType id_cell, part.m_Cells) {
 185       cell_marked[id_cell] = true;
 186     }
 187     QList<vtkIdType> new_cells;
 188     for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 189       if (cell_marked[id_cell]) {
 190         new_cells.append(id_cell);
 191       }
 192     }
 193     setCells(new_cells);
 194   } else {
 195     if (tol < 0) {
 196       tol *= -min(getSmallestEdgeLength(), part.getSmallestEdgeLength());
 197     }
 198     tol = max(tol, 1e-30);
 199     EG_VTKSP(vtkUnstructuredGrid, new_grid);
 200     EG_VTKSP(vtkKdTreePointLocator,loc);
 201     loc->SetDataSet(m_Grid);
 202     loc->BuildLocator();
 203     QVector<vtkIdType> pnode2node(part.m_Grid->GetNumberOfPoints());
 204     vtkIdType N = m_Grid->GetNumberOfPoints();
 205     Timer T(10);
 206     for (vtkIdType id_pnode = 0; id_pnode < part.m_Grid->GetNumberOfPoints(); ++id_pnode) {
 207       vec3_t xp, x;
 208       part.m_Grid->GetPoint(id_pnode, xp.data());
 209       vtkIdType id_node = loc->FindClosestPoint(xp.data());
 210       m_Grid->GetPoint(id_node, x.data());
 211       if ((x - xp).abs() < tol) {
 212         pnode2node[id_pnode] = id_node;
 213       } else {
 214         pnode2node[id_pnode] = N;
 215         ++N;
 216       }
 217     }
 218     allocateGrid(new_grid, m_Grid->GetNumberOfCells() + part.m_Cells.size(), N);
 219     for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 220       vec3_t x;
 221       m_Grid->GetPoint(id_node, x.data());
 222       new_grid->GetPoints()->SetPoint(id_node, x.data());
 223       copyNodeData(m_Grid, id_node, new_grid, id_node);
 224     }
 225     QVector<vtkIdType> part_nodes;
 226     getNodesFromCells(part.m_Cells, part_nodes, part.m_Grid);
 227     foreach (vtkIdType id_pnode, part_nodes) {
 228       vec3_t x;
 229       part.m_Grid->GetPoint(id_pnode, x.data());
 230       new_grid->GetPoints()->SetPoint(pnode2node[id_pnode], x.data());
 231       copyNodeData(part.m_Grid, id_pnode, new_grid, pnode2node[id_pnode]);
 232     }
 233     QList<vtkIdType> new_cells;
 234     for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 235       vtkIdType id_new_cell = copyCell(m_Grid, id_cell, new_grid);
 236       copyCellData(m_Grid, id_cell, new_grid, id_new_cell);
 237       new_cells.append(id_new_cell);
 238     }
 239     foreach (vtkIdType id_pcell, part.m_Cells) {
 240       vtkIdType id_new_cell = copyCell(part.m_Grid, id_pcell, new_grid, pnode2node);
 241       copyCellData(part.m_Grid, id_pcell, new_grid, id_new_cell);
 242       new_cells.append(id_new_cell);
 243     }
 244     makeCopy(new_grid, m_Grid);
 245   }
 246 }
 247
 248 void MeshPartition::concatenatePartition(const MeshPartition& part)
 249 {
 250   if (m_Grid == part.m_Grid) {
 251     QVector<bool> cell_marked(m_Grid->GetNumberOfCells(), false);
 252     foreach (vtkIdType id_cell, m_Cells) {
 253       cell_marked[id_cell] = true;
 254     }
 255     foreach (vtkIdType id_cell, part.m_Cells) {
 256       cell_marked[id_cell] = true;
 257     }
 258     QList<vtkIdType> new_cells;
 259     for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 260       if (cell_marked[id_cell]) {
 261         new_cells.append(id_cell);
 262       }
 263     }
 264     setCells(new_cells);
 265   } else {
 266     EG_VTKSP(vtkUnstructuredGrid, new_grid);
 267     QVector<vtkIdType> pnode2node(part.m_Grid->GetNumberOfPoints());
 268     vtkIdType N = m_Grid->GetNumberOfPoints();
 269     for (vtkIdType id_pnode = 0; id_pnode < part.m_Grid->GetNumberOfPoints(); ++id_pnode) {
 270       pnode2node[id_pnode] = N;
 271       ++N;
 272     }
 273     allocateGrid(new_grid, m_Grid->GetNumberOfCells() + part.m_Cells.size(), N);
 274     for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 275       vec3_t x;
 276       m_Grid->GetPoint(id_node, x.data());
 277       new_grid->GetPoints()->SetPoint(id_node, x.data());
 278       copyNodeData(m_Grid, id_node, new_grid, id_node);
 279     }
 280     QVector<vtkIdType> part_nodes;
 281     getNodesFromCells(part.m_Cells, part_nodes, part.m_Grid);
 282     foreach (vtkIdType id_pnode, part_nodes) {
 283       vec3_t x;
 284       part.m_Grid->GetPoint(id_pnode, x.data());
 285       new_grid->GetPoints()->SetPoint(pnode2node[id_pnode], x.data());
 286       copyNodeData(part.m_Grid, id_pnode, new_grid, pnode2node[id_pnode]);
 287     }
 288     QList<vtkIdType> new_cells;
 289     for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 290       vtkIdType id_new_cell = copyCell(m_Grid, id_cell, new_grid);
 291       copyCellData(m_Grid, id_cell, new_grid, id_new_cell);
 292       new_cells.append(id_new_cell);
 293     }
 294     foreach (vtkIdType id_pcell, part.m_Cells) {
 295       vtkIdType id_new_cell = copyCell(part.m_Grid, id_pcell, new_grid, pnode2node);
 296       copyCellData(part.m_Grid, id_pcell, new_grid, id_new_cell);
 297       new_cells.append(id_new_cell);
 298     }
 299     makeCopy(new_grid, m_Grid);
 300   }
 301 }
 302
 303 double MeshPartition::getSmallestEdgeLength() const
 304 {
 305   double L = 1e99;
 306   foreach (vtkIdType id_cell, m_Cells) {
 307     vtkIdType type_cell = m_Grid->GetCellType(id_cell);
 308     vtkIdType N_pts, *pts;
 309     m_Grid->GetCellPoints(id_cell, N_pts, pts);
 310     QVector<vec3_t> x(N_pts);
 311     for (int i = 0; i < N_pts; ++i) {
 312       m_Grid->GetPoint(pts[i], x[i].data());
 313     }
 314     if        (type_cell == VTK_TRIANGLE) {
 315       L = min(L, (x[0] - x[1]).abs());
 316       L = min(L, (x[1] - x[2]).abs());
 317       L = min(L, (x[2] - x[0]).abs());
 318     } else if (type_cell == VTK_QUAD) {
 319       L = min(L, (x[0] - x[1]).abs());
 320       L = min(L, (x[1] - x[2]).abs());
 321       L = min(L, (x[2] - x[3]).abs());
 322       L = min(L, (x[3] - x[0]).abs());
 323     } else if (type_cell == VTK_TETRA) {
 324       L = min(L, (x[0] - x[1]).abs());
 325       L = min(L, (x[1] - x[2]).abs());
 326       L = min(L, (x[2] - x[0]).abs());
 327       L = min(L, (x[3] - x[0]).abs());
 328       L = min(L, (x[3] - x[1]).abs());
 329       L = min(L, (x[3] - x[2]).abs());
 330     } else if (type_cell == VTK_PYRAMID) {
 331       L = min(L, (x[0] - x[1]).abs());
 332       L = min(L, (x[1] - x[2]).abs());
 333       L = min(L, (x[2] - x[3]).abs());
 334       L = min(L, (x[3] - x[0]).abs());
 335       L = min(L, (x[4] - x[0]).abs());
 336       L = min(L, (x[4] - x[1]).abs());
 337       L = min(L, (x[4] - x[2]).abs());
 338       L = min(L, (x[4] - x[3]).abs());
 339     } else if (type_cell == VTK_WEDGE) {
 340       L = min(L, (x[0] - x[1]).abs());
 341       L = min(L, (x[1] - x[2]).abs());
 342       L = min(L, (x[2] - x[0]).abs());
 343       L = min(L, (x[3] - x[4]).abs());
 344       L = min(L, (x[4] - x[5]).abs());
 345       L = min(L, (x[5] - x[3]).abs());
 346       L = min(L, (x[0] - x[3]).abs());
 347       L = min(L, (x[1] - x[4]).abs());
 348       L = min(L, (x[2] - x[5]).abs());
 349     } else if (type_cell == VTK_HEXAHEDRON) {
 350       L = min(L, (x[0] - x[1]).abs());
 351       L = min(L, (x[1] - x[2]).abs());
 352       L = min(L, (x[2] - x[3]).abs());
 353       L = min(L, (x[3] - x[0]).abs());
 354       L = min(L, (x[4] - x[5]).abs());
 355       L = min(L, (x[5] - x[6]).abs());
 356       L = min(L, (x[6] - x[7]).abs());
 357       L = min(L, (x[7] - x[4]).abs());
 358       L = min(L, (x[0] - x[4]).abs());
 359       L = min(L, (x[1] - x[5]).abs());
 360       L = min(L, (x[2] - x[6]).abs());
 361       L = min(L, (x[3] - x[7]).abs());
 362     }
 363   }
 364   return L;
 365 }
 366
 367 bool MeshPartition::hasNeighNode(vtkIdType id_node, vtkIdType id_neigh)
 368 {
 369   for (int i = 0; i < n2nGSize(id_node); ++i) {
 370     if (n2nGG(id_node, i) == id_neigh) {
 371       return true;
 372     }
 373   }
 374   return false;
 375 }
 376
 377 void MeshPartition::createNodeToBC()
 378 {
 379   EG_VTKDCC(vtkIntArray, cell_code,   m_Grid, "cell_code");
 380   m_N2BC.resize(m_Nodes.size());
 381   for (int i_node = 0; i_node < m_Nodes.size(); ++i_node) {
 382     QSet<int> bcs;
 383     for (int j = 0; j < n2cLSize(i_node); ++j) {
 384       vtkIdType id_cell = n2cLG(i_node, j);
 385       if (isSurface(id_cell, m_Grid)) {
 386         int bc = cell_code->GetValue(n2cLG(i_node, j));
 387         if (bc != 0) {
 388           bcs.insert(cell_code->GetValue(n2cLG(i_node, j)));
 389         }
 390       }
 391     }
 392     m_N2BC[i_node].resize(bcs.size());
 393     int i = 0;
 394     foreach (int bc, bcs) {
 395       m_N2BC[i_node][i++] = bc;
 396     }
 397   }
 398 }
 399
 400 bool MeshPartition::hasBC(vtkIdType id_node, int bc)
 401 {
 402   bool found = false;
 403   for (int j = 0; j < n2bcGSize(id_node); ++j) {
 404     if (n2bcG(id_node, j) == bc) {
 405       found = true;
 406       break;
 407     }
 408   }
 409   return found;
 410 }
 411
 412 vtkIdType MeshPartition::getVolumeCell(vtkIdType id_face)
 413 {
 414   checkLNodes();
 415   checkLCells();
 416   checkN2C();
 417   return findVolumeCell(m_Grid, id_face, m_LNodes, m_Cells, m_LCells, m_N2C);
 418 }
 419
 420 vec3_t MeshPartition::globalNormal(vtkIdType id_node)
 421 {
 422   vec3_t normal(0,0,0);
 423   for (int i = 0; i < n2cGSize(id_node); ++i) {
 424     vtkIdType id_cell = n2cGG(id_node, i);
 425     if (isSurface(id_cell, m_Grid)) {
 426       vtkIdType N_pts, *pts;
 427       m_Grid->GetCellPoints(id_cell, N_pts, pts);
 428       vec3_t a, b, c;
 429       for (int j = 0; j < N_pts; ++j) {
 430         if (pts[j] == id_node) {
 431           m_Grid->GetPoint(pts[j], a.data());
 432           if (j > 0) {
 433             m_Grid->GetPoint(pts[j-1], b.data());
 434           } else {
 435             m_Grid->GetPoint(pts[N_pts-1], b.data());
 436           }
 437           if (j < N_pts - 1) {
 438             m_Grid->GetPoint(pts[j+1], c.data());
 439           } else {
 440             m_Grid->GetPoint(pts[0], c.data());
 441           }
 442         }
 443       }
 444       vec3_t u = b - a;
 445       vec3_t v = c - a;
 446       double alpha = GeometryTools::angle(u, v);
 447       vec3_t n = u.cross(v);
 448       n.normalise();
 449       if (checkVector(n)) {
 450         normal -= alpha*n;
 451       }
 452     }
 453   }
 454   normal.normalise();
 455   return normal;
 456 }
 457
 458 double MeshPartition::getAverageSurfaceEdgeLength(vtkIdType id_node)
 459 {
 460   QSet<vtkIdType> surface_neighbours;
 461   for (int i = 0; i < n2cGSize(id_node); ++i) {
 462     vtkIdType id_cell = n2cGG(id_node, i);
 463     if (isSurface(id_cell, m_Grid)) {
 464       EG_GET_CELL(id_cell, m_Grid);
 465       for (int j = 0; j < num_pts; ++j) {
 466         if (pts[j] != id_node) {
 467           surface_neighbours.insert(pts[j]);
 468         }
 469       }
 470     }
 471   }
 472   double L = 0;
 473   if (surface_neighbours.size() > 0) {
 474     vec3_t x, xn;
 475     m_Grid->GetPoint(id_node, x.data());
 476     foreach (vtkIdType id_neigh, surface_neighbours) {
 477       m_Grid->GetPoint(id_neigh, xn.data());
 478       L += (x - xn).abs();
 479     }
 480     L /= surface_neighbours.size();
 481   }
 482   return L;
 483 }
 484
 485 void MeshPartition::computeMinAndMaxSurfaceStencilEdgeLengths(vtkIdType id_node, double &l_min, double &l_max)
 486 {
 487   l_min = 1e99;
 488   l_max = 0;
 489   for (int i = 0; i < n2cGSize(id_node); ++i) {
 490     vtkIdType id_cell = n2cGG(id_node, i);
 491     if (isSurface(id_cell, m_Grid)) {
 492       vtkIdType *pts, num_pts;
 493       m_Grid->GetCellPoints(id_cell, num_pts, pts);
 494       vec3_t x1, x2;
 495       m_Grid->GetPoint(pts[num_pts - 1], x1.data());
 496       for (int j = 0; j < num_pts; ++j) {
 497         m_Grid->GetPoint(pts[j], x2.data());
 498         double L = (x1 - x2).abs();
 499         l_min = min(L, l_min);
 500         l_max = max(L, l_max);
 501         x1 = x2;
 502       }
 503     }
 504   }
 505 }
 506
 507 double MeshPartition::getMinSurfaceStencilEdgeLength(vtkIdType id_node)
 508 {
 509   double l_min, l_max;
 510   computeMinAndMaxSurfaceStencilEdgeLengths(id_node, l_min, l_max);
 511   return l_min;
 512 }
 513
 514 double MeshPartition::getMaxSurfaceStencilEdgeLength(vtkIdType id_node)
 515 {
 516   double l_min, l_max;
 517   computeMinAndMaxSurfaceStencilEdgeLengths(id_node, l_min, l_max);
 518   return l_max;
 519 }
 520
 521 int MeshPartition::getNumberOfFeatureNeighbours(vtkIdType id_node)
 522 {
 523   EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type");
 524   int N = 0;
 525   for (int i = 0; i < n2nGSize(id_node); ++i) {
 526     char type = node_type->GetValue(n2nGG(id_node, i));
 527     if (type == EG_FEATURE_EDGE_VERTEX || type == EG_FEATURE_CORNER_VERTEX) {
 528       ++N;
 529     }
 530   }
 531   return N;
 532 }
 533
 534 int MeshPartition::getEdgeType(vtkIdType id_node1, vtkIdType id_node2)
 535 {
 536   QList <vtkIdType> edge_faces;
 537   getEdgeFaces(id_node1, id_node2, edge_faces);
 538   int edge = 0;
 539   if (edge_faces.size() == 1) {
 540     edge = 2;
 541   } else if (edge_faces.size() >= 2) {
 542     EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
 543     QSet<int> bcs;
 544     foreach (vtkIdType id_face, edge_faces) {
 545       bcs.insert(cell_code->GetValue(id_face));
 546     }
 547     if (bcs.size() > 1) {
 548       edge = 2;
 549     } else {
 550       edge = 1;
 551     }
 552   }
 553   return edge;
 554 }
 555
 556 int MeshPartition::computeTopoDistance(vtkIdType id_node1, vtkIdType id_node2, int max_dist, int restriction_type)
 557 {
 558   int dist = 0;
 559   QSet<vtkIdType> candidates;
 560   candidates.insert(id_node1);
 561   while (dist < max_dist && !candidates.contains(id_node2)) {
 562     foreach (vtkIdType id_node, candidates) {
 563       for (int i = 0; i < n2nGSize(id_node); ++i) {
 564         vtkIdType id_neigh = n2nGG(id_node,i);
 565         if (!candidates.contains(id_neigh)) {
 566           bool insert = true;
 567           if (restriction_type > 0) {
 568             insert = getEdgeType(id_node, id_neigh) >= restriction_type;
 569           }
 570           if (insert) {
 571             candidates.insert(id_neigh);
 572           }
 573         }
 574       }
 575     }
 576     ++dist;
 577   }
 578   return dist;
 579 }
 580
 581 void MeshPartition::getCommonNodes(vtkIdType id_cell1, vtkIdType id_cell2, QVector<vtkIdType> &common_nodes)
 582 {
 583   common_nodes.clear();
 584   QSet<vtkIdType> nodes1, nodes2;
 585   vtkIdType num_pts, *pts;
 586   m_Grid->GetCellPoints(id_cell1, num_pts, pts);
 587   for (int i = 0; i < num_pts; ++i) {
 588     nodes1.insert(pts[i]);
 589   }
 590   m_Grid->GetCellPoints(id_cell2, num_pts, pts);
 591   for (int i = 0; i < num_pts; ++i) {
 592     nodes2.insert(pts[i]);
 593   }
 594   nodes1.intersect(nodes2);
 595   common_nodes.resize(nodes1.size());
 596   qCopy(nodes1.begin(), nodes1.end(), common_nodes.begin());
 597 }
 598
 599 bool MeshPartition::isFeatureEdge(vtkIdType id_node1, vtkIdType id_node2, double feature_angle)
 600 {
 601   bool is_feature_edge = false;
 602   QVector<int> bcs;
 603   QVector<vtkIdType> nodes(2);
 604   nodes[0] = id_node1;
 605   nodes[1] = id_node2;
 606   getCommonBcs(nodes, bcs);
 607   if (bcs.size() == 1) {
 608     QList<vtkIdType> faces;
 609     getEdgeFaces(id_node1, id_node2, faces);
 610     if (faces.size() == 2) {
 611       vec3_t n1 = cellNormal(m_Grid, faces[0]);
 612       vec3_t n2 = cellNormal(m_Grid, faces[1]);
 613       if (GeometryTools::angle(n1, n2) > feature_angle) {
 614         EG_VTKDCN(vtkDoubleArray, cl, m_Grid, "node_meshdensity_desired");
 615         CadInterface *cad = GuiMainWindow::pointer()->getCadInterface(bcs.first());
 616         vec3_t x1, x2;
 617         m_Grid->GetPoint(id_node1, x1.data());
 618         m_Grid->GetPoint(id_node2, x2.data());
 619         vec3_t x = 0.5*(x1 + x2);
 620         double h = min(cl->GetValue(id_node1), cl->GetValue(id_node2));
 621         vec3_t xs = cad->snapToEdge(x);
 622         if ((xs - x).abs() < 0.2*h) {
 623           is_feature_edge = true;
 624         }
 625       }
 626     }
 627   }
 628   return is_feature_edge;
 629 }
 630
 631 bool MeshPartition::isConvexNode(vtkIdType id_node)
 632 {
 633   int n_faces = n2cGSize(id_node);
 634   if (n_faces <= 1) {
 635     return false;
 636   }
 637
 638   vec3_t x1, cell_centers(0,0,0), cell_normals(0,0,0);
 639   m_Grid->GetPoint(id_node, x1.data());
 640   for (int i = 0; i < n_faces; ++i) {
 641     vtkIdType id_cell = n2cGG(id_node, i);
 642     cell_centers += cellCentre(m_Grid, id_cell);
 643     cell_normals += cellNormal(m_Grid, id_cell);
 644   }
 645   cell_centers *= 1.0/n_faces;
 646   if ((x1 - cell_centers)*cell_normals.normalise() > 0) {
 647     return true;
 648   }
 649   return false;
 650 }
 651
 652
 653 bool MeshPartition::isConvexNode(vtkIdType id_node, QVector<int> bl_codes)
 654 {
 655   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
 656   int n_faces = n2cGSize(id_node);
 657   if (n_faces <= 1) {
 658     return false;
 659   }
 660
 661   vec3_t x1, cell_centers(0,0,0), cell_normals(0,0,0);
 662   m_Grid->GetPoint(id_node, x1.data());
 663   int n_used = 0;
 664   for (int i = 0; i < n_faces; ++i) {
 665     vtkIdType id_cell = n2cGG(id_node, i);
 666     if (bl_codes.contains(cell_code->GetValue(id_cell))) {
 667       cell_centers += cellCentre(m_Grid, id_cell);
 668       cell_normals += cellNormal(m_Grid, id_cell);
 669       ++n_used;
 670     }
 671   }
 672   if (n_used > 0) {
 673     cell_centers *= 1.0/n_used;
 674     if ((x1 - cell_centers)*cell_normals.normalise() > 0) {
 675       return true;
 676     }
 677   }
 678   return false;
 679 }
 680
 681 void MeshPartition::calcPlanarSurfaceMetrics(double &Dh, double &A, double &P, vec3_t &x, vec3_t &n)
 682 {
 683   A = 0;
 684   n = vec3_t(0, 0, 0);
 685   x = vec3_t(0, 0, 0);
 686   P = 0;
 687   foreach (vtkIdType id_cell, m_Cells) {
 688     if (isSurface(id_cell, m_Grid)) {
 689       vtkIdType num_pts, *pts;
 690       m_Grid->GetCellPoints(id_cell, num_pts, pts);
 691       QVector<vec3_t> x_pt(num_pts + 1);
 692       for (int i = 0; i < num_pts; ++i) {
 693         m_Grid->GetPoint(pts[i], x_pt[i].data());
 694       }
 695       x_pt[num_pts] = x_pt[0];
 696       double a = cellVA(m_Grid, id_cell);
 697       A += a;
 698       x += a*cellCentre(m_Grid, id_cell);
 699       n += cellNormal(m_Grid, id_cell);
 700       for (int i = 0; i < num_pts; ++i) {
 701         if (c2cGG(id_cell, i) == -1) {
 702           P += (x_pt[i+1] - x_pt[i]).abs();
 703         }
 704       }
 705     }
 706   }
 707   n.normalise();
 708   x *= 1.0/A;
 709   Dh = 4*A/P;
 710 }
 711
 712 bool MeshPartition::isPlanar(double tolerance_angle)
 713 {
 714   foreach (vtkIdType id_cell, m_Cells) {
 715     if (!isSurface(id_cell, m_Grid)) {
 716       return false;
 717     }
 718   }
 719   double Dh, A, P;
 720   vec3_t x, n0;
 721   calcPlanarSurfaceMetrics(Dh, A, P, n0, x);
 722   foreach (vtkIdType id_cell, m_Cells) {
 723     vec3_t n = cellNormal(m_Grid, id_cell);
 724     if (angle(n, n0) > tolerance_angle) {
 725       return false;
 726     }
 727   }
 728   return true;
 729 }
 730
 731 void MeshPartition::setBC(int bc)
 732 {
 733   QList<vtkIdType> cls;
 734   EG_VTKDCC(vtkIntArray, cell_code,   m_Grid, "cell_code");
 735   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 736     if (cell_code->GetValue(id_cell) == bc) {
 737       cls.append(id_cell);
 738     }
 739   }
 740   setCells(cls);
 741 }
 742
 743 void MeshPartition::duplicate()
 744 {
 745   vtkIdType old_num_cells = m_Grid->GetNumberOfCells();
 746   EG_VTKSP(vtkUnstructuredGrid, new_grid);
 747   extractToVtkGrid(new_grid);
 748   MeshPartition new_part(new_grid, true);
 749   concatenatePartition(new_part);
 750   QList<vtkIdType> cells;
 751   for (vtkIdType id_cell = old_num_cells; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 752     cells << id_cell;
 753   }
 754   setCells(cells);
 755 }
 756
 757 void MeshPartition::scale(double factor, vec3_t centre)
 758 {
 759   checkNodes();
 760   foreach (vtkIdType id_node, m_Nodes) {
 761     vec3_t x;
 762     m_Grid->GetPoint(id_node, x.data());
 763     x -= centre;
 764     x *= factor;
 765     x += centre;
 766     m_Grid->GetPoints()->SetPoint(id_node, x.data());
 767   }
 768 }
 769
 770 void MeshPartition::translate(vec3_t v)
 771 {
 772   checkNodes();
 773   foreach (vtkIdType id_node, m_Nodes) {
 774     vec3_t x;
 775     m_Grid->GetPoint(id_node, x.data());
 776     x += v;
 777     m_Grid->GetPoints()->SetPoint(id_node, x.data());
 778   }
 779 }
 780
 781 void MeshPartition::extrude(vec3_t dir, QList<double> h, BoundaryCondition extrude_bc,
 782                             bool force_bottom_bc, bool force_side_bc, bool force_top_bc,
 783                             BoundaryCondition bottom_bc, BoundaryCondition side_bc, BoundaryCondition top_bc)
 784 {
 785   checkNodes();
 786   if (onlySurfaceCells()) {
 787     QList<vtkIdType> new_cells;
 788     foreach (vtkIdType id_cell, m_Cells) {
 789       new_cells << id_cell;
 790     }
 791     vtkIdType old_num_cells = m_Grid->GetNumberOfCells();
 792     extrude_bc = GuiMainWindow::pointer()->getBC(extrude_bc);
 793     EG_VTKDCC(vtkIntArray, cell_code,   m_Grid, "cell_code");
 794     foreach (vtkIdType id_cell, m_Cells) {
 795       cell_code->SetValue(id_cell, extrude_bc.getCode());
 796     }
 797     QSet<int> bcs;
 798     bcs.insert(extrude_bc.getCode());
 799     EG_VTKSP(vtkEgNormalExtrusion, extr);
 800     QVector<double> y(h.size() + 1);
 801     {
 802       y[0] = 0;
 803       int i = 1;
 804       foreach (double dy, h) {
 805         y[i] = y[i-1] + dy;
 806         ++i;
 807       }
 808     }
 809     extr->SetLayers(y);
 810     extr->SetBoundaryCodes(bcs);
 811     if (force_bottom_bc) extr->SetCustomBottomBc (bottom_bc.getCode());
 812     if (force_side_bc)   extr->SetCustomSideBc   (side_bc.getCode());
 813     if (force_top_bc)    extr->SetCustomTopBc    (top_bc.getCode());
 814     dir.normalise();
 815     extr->SetNormal(dir);
 816     extr->SetFixed();
 817     EG_VTKSP(vtkUnstructuredGrid,ug);
 818     makeCopy(m_Grid, ug);
 819     extr->SetInputData(ug);
 820     extr->Update();
 821     makeCopy(extr->GetOutput(), m_Grid);
 822     setBC(extrude_bc.getCode());
 823   }
 824 }
 825
 826 bool MeshPartition::onlySurfaceCells()
 827 {
 828   foreach (vtkIdType id_cell, m_Cells) {
 829     if (!isSurface(id_cell, m_Grid)) {
 830       return false;
 831     }
 832   }
 833   return true;
 834 }
 835
 836 void MeshPartition::resetBC(QString bc_name, QString bc_type)
 837 {
 838   BoundaryCondition bc = GuiMainWindow::pointer()->getBC(BoundaryCondition(bc_name, bc_type));
 839   EG_VTKDCC(vtkIntArray, cell_code,   m_Grid, "cell_code");
 840   foreach (vtkIdType id_cell, m_Cells) {
 841     cell_code->SetValue(id_cell, bc.getCode());
 842   }
 843 }
 844
 845 void MeshPartition::writeSTL(QString file_name)
 846 {
 847   if (file_name.right(4).toLower() != ".stl") {
 848     file_name += ".stl";
 849   }
 850   EG_VTKSP(vtkUnstructuredGrid, grid);
 851   extractToVtkGrid(grid);
 852   EG_VTKSP(vtkGeometryFilter, geometry);
 853   geometry->SetInputData(grid);
 854   EG_VTKSP(vtkTriangleFilter, triangle);
 855   triangle->SetInputConnection(geometry->GetOutputPort());
 856   EG_VTKSP(vtkSTLWriter, stl);
 857   stl->SetInputConnection(  triangle->GetOutputPort());
 858   stl->SetFileName(qPrintable(file_name));
 859   stl->SetFileTypeToBinary();
 860   stl->Write();
 861 }