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