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
  28 MeshPartition::MeshPartition()
  29 {
  30   m_Grid = NULL;
  31   m_TrackGrid = false;
  32   resetTimeStamps();
  33 }
  34
  35 MeshPartition::MeshPartition(vtkUnstructuredGrid *grid, bool use_all_cells)
  36 {
  37   setGrid(grid, use_all_cells);
  38 }
  39
  40 MeshPartition::MeshPartition(QString volume_name)
  41 {
  42   m_TrackGrid = false;
  43   resetTimeStamps();
  44   setVolume(volume_name);
  45 }
  46
  47 void MeshPartition::setGrid(vtkUnstructuredGrid *grid, bool use_all_cells)
  48 {
  49   m_TrackGrid = false;
  50   resetTimeStamps();
  51   m_Grid = grid;
  52   if (use_all_cells) {
  53     QVector<vtkIdType> cls(grid->GetNumberOfCells());
  54     for (vtkIdType id_cell = 0; id_cell < cls.size(); ++id_cell) {
  55       cls[id_cell] = id_cell;
  56     }
  57     setCells(cls);
  58   }
  59 }
  60
  61
  62 void MeshPartition::resetTimeStamps()
  63 {
  64   m_CellsStamp = 0;
  65   m_LCellsStamp = 0;
  66   m_NodesStamp = 0;
  67   m_LNodesStamp = 0;
  68   m_N2NStamp = 0;
  69   m_N2CStamp = 0;
  70   m_N2BCStamp = 0;
  71   m_C2CStamp = 0;
  72 }
  73
  74 void MeshPartition::trackGrid(vtkUnstructuredGrid *grid)
  75 {
  76   setGrid(grid);
  77   setAllCells();
  78   m_GridMTime = m_Grid->GetMTime();
  79   m_TrackGrid = true;
  80 }
  81
  82 void MeshPartition::setVolume(QString volume_name)
  83 {
  84   m_Grid = GuiMainWindow::pointer()->getGrid();
  85   resetOrientation(m_Grid);
  86   VolumeDefinition V = GuiMainWindow::pointer()->getVol(volume_name);
  87   QList<vtkIdType> cls;
  88   EG_VTKDCC(vtkIntArray, cell_code,   m_Grid, "cell_code");
  89   EG_VTKDCC(vtkIntArray, cell_orgdir, m_Grid, "cell_orgdir");
  90   EG_VTKDCC(vtkIntArray, cell_curdir, m_Grid, "cell_curdir");
  91   EG_VTKDCC(vtkIntArray, cell_voldir, m_Grid, "cell_voldir");
  92   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
  93     if (isSurface(id_cell, m_Grid)) {
  94       int bc = cell_code->GetValue(id_cell);
  95       cell_voldir->SetValue(id_cell, 0);
  96       if (V.getSign(bc) != 0) {
  97         cls.append(id_cell);
  98         if (V.getSign(bc) == -1) {
  99           cell_voldir->SetValue(id_cell, 1);
 100         }
 101       }
 102     } else {
 103       if (cell_code->GetValue(id_cell) == V.getVC()) {
 104         cls.append(id_cell);
 105       }
 106     }
 107   }
 108   setCells(cls);
 109 }
 110
 111 void MeshPartition::setVolumeOrientation()
 112 {
 113   EG_VTKDCC(vtkIntArray, cell_curdir, m_Grid, "cell_curdir");
 114   EG_VTKDCC(vtkIntArray, cell_voldir, m_Grid, "cell_voldir");
 115   foreach (vtkIdType id_cell, m_Cells) {
 116     if (isSurface(id_cell, m_Grid)) {
 117       if (cell_curdir->GetValue(id_cell) != cell_voldir->GetValue(id_cell)) {
 118         reorientateFace(m_Grid, id_cell);
 119       }
 120     }
 121   }
 122 }
 123
 124 void MeshPartition::setOriginalOrientation()
 125 {
 126   EG_VTKDCC(vtkIntArray, cell_curdir, m_Grid, "cell_curdir");
 127   EG_VTKDCC(vtkIntArray, cell_orgdir, m_Grid, "cell_orgdir");
 128   foreach (vtkIdType id_cell, m_Cells) {
 129     if (isSurface(id_cell, m_Grid)) {
 130       if (cell_curdir->GetValue(id_cell) != cell_orgdir->GetValue(id_cell)) {
 131         reorientateFace(m_Grid, id_cell);
 132       }
 133     }
 134   }
 135 }
 136
 137 void MeshPartition::setRemainder(const MeshPartition& part)
 138 {
 139   setGrid(part.getGrid());
 140   QVector<vtkIdType> rcells;
 141   getRestCells(m_Grid, part.m_Cells, rcells);
 142   setCells(rcells);
 143 }
 144
 145 void MeshPartition::extractToVtkGrid(vtkUnstructuredGrid *new_grid)
 146 {
 147   checkLNodes();
 148   allocateGrid(new_grid, m_Cells.size(), m_Nodes.size());
 149   for (int i_nodes = 0; i_nodes < m_Nodes.size(); ++i_nodes) {
 150     vec3_t x;
 151     m_Grid->GetPoints()->GetPoint(m_Nodes[i_nodes], x.data());
 152     new_grid->GetPoints()->SetPoint(i_nodes, x.data());
 153     copyNodeData(m_Grid, m_Nodes[i_nodes], new_grid, i_nodes);
 154   }
 155   foreach (vtkIdType id_cell, m_Cells) {
 156     /*
 157     vtkIdType N_pts, *pts;
 158     vtkIdType type_cell = m_Grid->GetCellType(id_cell);
 159     m_Grid->GetCellPoints(id_cell, N_pts, pts);
 160     QVector<vtkIdType> new_pts(N_pts);
 161     for (int i = 0; i < N_pts; ++i) {
 162       new_pts[i] = m_LNodes[pts[i]];
 163     }
 164     // update for polyhedral cells here
 165     vtkIdType id_new_cell = new_grid->InsertNextCell(type_cell, N_pts, new_pts.data());
 166     */
 167     vtkIdType id_new_cell = copyCell(m_Grid, id_cell, new_grid, m_LNodes);
 168     copyCellData(m_Grid, id_cell, new_grid, id_new_cell);
 169   }
 170 }
 171
 172 void MeshPartition::addPartition(const MeshPartition& part, double tol)
 173 {
 174   if (m_Grid == part.m_Grid) {
 175     QVector<bool> cell_marked(m_Grid->GetNumberOfCells(), false);
 176     foreach (vtkIdType id_cell, m_Cells) {
 177       cell_marked[id_cell] = true;
 178     }
 179     foreach (vtkIdType id_cell, part.m_Cells) {
 180       cell_marked[id_cell] = true;
 181     }
 182     QList<vtkIdType> new_cells;
 183     for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 184       if (cell_marked[id_cell]) {
 185         new_cells.append(id_cell);
 186       }
 187     }
 188     setCells(new_cells);
 189   } else {
 190     if (tol < 0) {
 191       tol *= -min(getSmallestEdgeLength(), part.getSmallestEdgeLength());
 192     }
 193     tol = max(tol, 1e-30);
 194     EG_VTKSP(vtkUnstructuredGrid, new_grid);
 195     EG_VTKSP(vtkKdTreePointLocator,loc);
 196     loc->SetDataSet(m_Grid);
 197     loc->BuildLocator();
 198     QVector<vtkIdType> pnode2node(part.m_Grid->GetNumberOfPoints());
 199     vtkIdType N = m_Grid->GetNumberOfPoints();
 200     Timer T(10);
 201     for (vtkIdType id_pnode = 0; id_pnode < part.m_Grid->GetNumberOfPoints(); ++id_pnode) {
 202       vec3_t xp, x;
 203       part.m_Grid->GetPoint(id_pnode, xp.data());
 204       vtkIdType id_node = loc->FindClosestPoint(xp.data());
 205       m_Grid->GetPoint(id_node, x.data());
 206       if ((x - xp).abs() < tol) {
 207         pnode2node[id_pnode] = id_node;
 208       } else {
 209         pnode2node[id_pnode] = N;
 210         ++N;
 211       }
 212     }
 213     allocateGrid(new_grid, m_Grid->GetNumberOfCells() + part.m_Cells.size(), N);
 214     for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 215       vec3_t x;
 216       m_Grid->GetPoint(id_node, x.data());
 217       new_grid->GetPoints()->SetPoint(id_node, x.data());
 218       copyNodeData(m_Grid, id_node, new_grid, id_node);
 219     }
 220     QVector<vtkIdType> part_nodes;
 221     getNodesFromCells(part.m_Cells, part_nodes, part.m_Grid);
 222     foreach (vtkIdType id_pnode, part_nodes) {
 223       vec3_t x;
 224       part.m_Grid->GetPoint(id_pnode, x.data());
 225       new_grid->GetPoints()->SetPoint(pnode2node[id_pnode], x.data());
 226       copyNodeData(part.m_Grid, id_pnode, new_grid, pnode2node[id_pnode]);
 227     }
 228     QList<vtkIdType> new_cells;
 229     for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 230       vtkIdType id_new_cell = copyCell(m_Grid, id_cell, new_grid);
 231       copyCellData(m_Grid, id_cell, new_grid, id_new_cell);
 232       new_cells.append(id_new_cell);
 233     }
 234     foreach (vtkIdType id_pcell, part.m_Cells) {
 235       vtkIdType id_new_cell = copyCell(part.m_Grid, id_pcell, new_grid, pnode2node);
 236       copyCellData(part.m_Grid, id_pcell, new_grid, id_new_cell);
 237       new_cells.append(id_new_cell);
 238     }
 239     makeCopy(new_grid, m_Grid);
 240   }
 241 }
 242
 243 double MeshPartition::getSmallestEdgeLength() const
 244 {
 245   double L = 1e99;
 246   foreach (vtkIdType id_cell, m_Cells) {
 247     vtkIdType type_cell = m_Grid->GetCellType(id_cell);
 248     vtkIdType N_pts, *pts;
 249     m_Grid->GetCellPoints(id_cell, N_pts, pts);
 250     QVector<vec3_t> x(N_pts);
 251     for (int i = 0; i < N_pts; ++i) {
 252       m_Grid->GetPoint(pts[i], x[i].data());
 253     }
 254     if        (type_cell == VTK_TRIANGLE) {
 255       L = min(L, (x[0] - x[1]).abs());
 256       L = min(L, (x[1] - x[2]).abs());
 257       L = min(L, (x[2] - x[0]).abs());
 258     } else if (type_cell == VTK_QUAD) {
 259       L = min(L, (x[0] - x[1]).abs());
 260       L = min(L, (x[1] - x[2]).abs());
 261       L = min(L, (x[2] - x[3]).abs());
 262       L = min(L, (x[3] - x[0]).abs());
 263     } else if (type_cell == VTK_TETRA) {
 264       L = min(L, (x[0] - x[1]).abs());
 265       L = min(L, (x[1] - x[2]).abs());
 266       L = min(L, (x[2] - x[0]).abs());
 267       L = min(L, (x[3] - x[0]).abs());
 268       L = min(L, (x[3] - x[1]).abs());
 269       L = min(L, (x[3] - x[2]).abs());
 270     } else if (type_cell == VTK_PYRAMID) {
 271       L = min(L, (x[0] - x[1]).abs());
 272       L = min(L, (x[1] - x[2]).abs());
 273       L = min(L, (x[2] - x[3]).abs());
 274       L = min(L, (x[3] - x[0]).abs());
 275       L = min(L, (x[4] - x[0]).abs());
 276       L = min(L, (x[4] - x[1]).abs());
 277       L = min(L, (x[4] - x[2]).abs());
 278       L = min(L, (x[4] - x[3]).abs());
 279     } else if (type_cell == VTK_WEDGE) {
 280       L = min(L, (x[0] - x[1]).abs());
 281       L = min(L, (x[1] - x[2]).abs());
 282       L = min(L, (x[2] - x[0]).abs());
 283       L = min(L, (x[3] - x[4]).abs());
 284       L = min(L, (x[4] - x[5]).abs());
 285       L = min(L, (x[5] - x[3]).abs());
 286       L = min(L, (x[0] - x[3]).abs());
 287       L = min(L, (x[1] - x[4]).abs());
 288       L = min(L, (x[2] - x[5]).abs());
 289     } else if (type_cell == VTK_HEXAHEDRON) {
 290       L = min(L, (x[0] - x[1]).abs());
 291       L = min(L, (x[1] - x[2]).abs());
 292       L = min(L, (x[2] - x[3]).abs());
 293       L = min(L, (x[3] - x[0]).abs());
 294       L = min(L, (x[4] - x[5]).abs());
 295       L = min(L, (x[5] - x[6]).abs());
 296       L = min(L, (x[6] - x[7]).abs());
 297       L = min(L, (x[7] - x[4]).abs());
 298       L = min(L, (x[0] - x[4]).abs());
 299       L = min(L, (x[1] - x[5]).abs());
 300       L = min(L, (x[2] - x[6]).abs());
 301       L = min(L, (x[3] - x[7]).abs());
 302     }
 303   }
 304   return L;
 305 }
 306
 307 bool MeshPartition::hasNeighNode(vtkIdType id_node, vtkIdType id_neigh)
 308 {
 309   for (int i = 0; i < n2nGSize(id_node); ++i) {
 310     if (n2nGG(id_node, i) == id_neigh) {
 311       return true;
 312     }
 313   }
 314   return false;
 315 }
 316
 317 void MeshPartition::createNodeToBC()
 318 {
 319   EG_VTKDCC(vtkIntArray, cell_code,   m_Grid, "cell_code");
 320   m_N2BC.resize(m_Nodes.size());
 321   for (int i_node = 0; i_node < m_Nodes.size(); ++i_node) {
 322     QSet<int> bcs;
 323     for (int j = 0; j < n2cLSize(i_node); ++j) {
 324       vtkIdType id_cell = n2cLG(i_node, j);
 325       if (isSurface(id_cell, m_Grid)) {
 326         int bc = cell_code->GetValue(n2cLG(i_node, j));
 327         if (bc != 0) {
 328           bcs.insert(cell_code->GetValue(n2cLG(i_node, j)));
 329         }
 330       }
 331     }
 332     m_N2BC[i_node].resize(bcs.size());
 333     int i = 0;
 334     foreach (int bc, bcs) {
 335       m_N2BC[i_node][i++] = bc;
 336     }
 337   }
 338 }
 339
 340 bool MeshPartition::hasBC(vtkIdType id_node, int bc)
 341 {
 342   bool found = false;
 343   for (int j = 0; j < n2bcGSize(id_node); ++j) {
 344     if (n2bcG(id_node, j) == bc) {
 345       found = true;
 346       break;
 347     }
 348   }
 349   return found;
 350 }
 351
 352 vtkIdType MeshPartition::getVolumeCell(vtkIdType id_face)
 353 {
 354   checkLNodes();
 355   checkLCells();
 356   checkN2C();
 357   return findVolumeCell(m_Grid, id_face, m_LNodes, m_Cells, m_LCells, m_N2C);
 358 }
 359
 360 vec3_t MeshPartition::globalNormal(vtkIdType id_node)
 361 {
 362   vec3_t normal(0,0,0);
 363   for (int i = 0; i < n2cGSize(id_node); ++i) {
 364     vtkIdType id_cell = n2cGG(id_node, i);
 365     if (isSurface(id_cell, m_Grid)) {
 366       vtkIdType N_pts, *pts;
 367       m_Grid->GetCellPoints(id_cell, N_pts, pts);
 368       vec3_t a, b, c;
 369       for (int j = 0; j < N_pts; ++j) {
 370         if (pts[j] == id_node) {
 371           m_Grid->GetPoint(pts[j], a.data());
 372           if (j > 0) {
 373             m_Grid->GetPoint(pts[j-1], b.data());
 374           } else {
 375             m_Grid->GetPoint(pts[N_pts-1], b.data());
 376           }
 377           if (j < N_pts - 1) {
 378             m_Grid->GetPoint(pts[j+1], c.data());
 379           } else {
 380             m_Grid->GetPoint(pts[0], c.data());
 381           }
 382         }
 383       }
 384       vec3_t u = b - a;
 385       vec3_t v = c - a;
 386       double alpha = GeometryTools::angle(u, v);
 387       vec3_t n = u.cross(v);
 388       n.normalise();
 389       if (checkVector(n)) {
 390         normal -= alpha*n;
 391       }
 392     }
 393   }
 394   normal.normalise();
 395   return normal;
 396 }
 397
 398 double MeshPartition::getAverageSurfaceEdgeLength(vtkIdType id_node)
 399 {
 400   QSet<vtkIdType> surface_neighbours;
 401   for (int i = 0; i < n2cGSize(id_node); ++i) {
 402     vtkIdType id_cell = n2cGG(id_node, i);
 403     if (isSurface(id_cell, m_Grid)) {
 404       EG_GET_CELL(id_cell, m_Grid);
 405       for (int j = 0; j < num_pts; ++j) {
 406         if (pts[j] != id_node) {
 407           surface_neighbours.insert(pts[j]);
 408         }
 409       }
 410     }
 411   }
 412   double L = 0;
 413   if (surface_neighbours.size() > 0) {
 414     vec3_t x, xn;
 415     m_Grid->GetPoint(id_node, x.data());
 416     foreach (vtkIdType id_neigh, surface_neighbours) {
 417       m_Grid->GetPoint(id_neigh, xn.data());
 418       L += (x - xn).abs();
 419     }
 420     L /= surface_neighbours.size();
 421   }
 422   return L;
 423 }
 424
 425 void MeshPartition::computeMinAndMaxSurfaceStencilEdgeLengths(vtkIdType id_node, double &l_min, double &l_max)
 426 {
 427   l_min = 1e99;
 428   l_max = 0;
 429   for (int i = 0; i < n2cGSize(id_node); ++i) {
 430     vtkIdType id_cell = n2cGG(id_node, i);
 431     if (isSurface(id_cell, m_Grid)) {
 432       vtkIdType *pts, num_pts;
 433       m_Grid->GetCellPoints(id_cell, num_pts, pts);
 434       vec3_t x1, x2;
 435       m_Grid->GetPoint(pts[num_pts - 1], x1.data());
 436       for (int j = 0; j < num_pts; ++j) {
 437         m_Grid->GetPoint(pts[j], x2.data());
 438         double L = (x1 - x2).abs();
 439         l_min = min(L, l_min);
 440         l_max = max(L, l_max);
 441         x1 = x2;
 442       }
 443     }
 444   }
 445 }
 446
 447 double MeshPartition::getMinSurfaceStencilEdgeLength(vtkIdType id_node)
 448 {
 449   double l_min, l_max;
 450   computeMinAndMaxSurfaceStencilEdgeLengths(id_node, l_min, l_max);
 451   return l_min;
 452 }
 453
 454 double MeshPartition::getMaxSurfaceStencilEdgeLength(vtkIdType id_node)
 455 {
 456   double l_min, l_max;
 457   computeMinAndMaxSurfaceStencilEdgeLengths(id_node, l_min, l_max);
 458   return l_max;
 459 }
 460
 461 int MeshPartition::getNumberOfFeatureNeighbours(vtkIdType id_node)
 462 {
 463   EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type");
 464   int N = 0;
 465   for (int i = 0; i < n2nGSize(id_node); ++i) {
 466     char type = node_type->GetValue(n2nGG(id_node, i));
 467     if (type == EG_FEATURE_EDGE_VERTEX || type == EG_FEATURE_CORNER_VERTEX) {
 468       ++N;
 469     }
 470   }
 471   return N;
 472 }
 473
 474 int MeshPartition::getEdgeType(vtkIdType id_node1, vtkIdType id_node2)
 475 {
 476   QList <vtkIdType> edge_faces;
 477   getEdgeFaces(id_node1, id_node2, edge_faces);
 478   int edge = 0;
 479   if (edge_faces.size() == 1) {
 480     edge = 2;
 481   } else if (edge_faces.size() >= 2) {
 482     EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
 483     QSet<int> bcs;
 484     foreach (vtkIdType id_face, edge_faces) {
 485       bcs.insert(cell_code->GetValue(id_face));
 486     }
 487     if (bcs.size() > 1) {
 488       edge = 2;
 489     } else {
 490       edge = 1;
 491     }
 492   }
 493   return edge;
 494 }
 495
 496 int MeshPartition::computeTopoDistance(vtkIdType id_node1, vtkIdType id_node2, int max_dist, int restriction_type)
 497 {
 498   int dist = 0;
 499   QSet<vtkIdType> candidates;
 500   candidates.insert(id_node1);
 501   while (dist < max_dist && !candidates.contains(id_node2)) {
 502     foreach (vtkIdType id_node, candidates) {
 503       for (int i = 0; i < n2nGSize(id_node); ++i) {
 504         vtkIdType id_neigh = n2nGG(id_node,i);
 505         if (!candidates.contains(id_neigh)) {
 506           bool insert = true;
 507           if (restriction_type > 0) {
 508             insert = getEdgeType(id_node, id_neigh) >= restriction_type;
 509           }
 510           if (insert) {
 511             candidates.insert(id_neigh);
 512           }
 513         }
 514       }
 515     }
 516     ++dist;
 517   }
 518   return dist;
 519 }
 520
 521 void MeshPartition::getCommonNodes(vtkIdType id_cell1, vtkIdType id_cell2, QVector<vtkIdType> &common_nodes)
 522 {
 523   common_nodes.clear();
 524   QSet<vtkIdType> nodes1, nodes2;
 525   vtkIdType num_pts, *pts;
 526   m_Grid->GetCellPoints(id_cell1, num_pts, pts);
 527   for (int i = 0; i < num_pts; ++i) {
 528     nodes1.insert(pts[i]);
 529   }
 530   m_Grid->GetCellPoints(id_cell2, num_pts, pts);
 531   for (int i = 0; i < num_pts; ++i) {
 532     nodes2.insert(pts[i]);
 533   }
 534   nodes1.intersect(nodes2);
 535   common_nodes.resize(nodes1.size());
 536   qCopy(nodes1.begin(), nodes1.end(), common_nodes.begin());
 537 }
 538
 539 bool MeshPartition::isFeatureEdge(vtkIdType id_node1, vtkIdType id_node2, double feature_angle)
 540 {
 541   bool is_feature_edge = false;
 542   QVector<int> bcs;
 543   QVector<vtkIdType> nodes(2);
 544   nodes[0] = id_node1;
 545   nodes[1] = id_node2;
 546   getCommonBcs(nodes, bcs);
 547   if (bcs.size() == 1) {
 548     QList<vtkIdType> faces;
 549     getEdgeFaces(id_node1, id_node2, faces);
 550     if (faces.size() == 2) {
 551       vec3_t n1 = cellNormal(m_Grid, faces[0]);
 552       vec3_t n2 = cellNormal(m_Grid, faces[1]);
 553       if (GeometryTools::angle(n1, n2) > feature_angle) {
 554         EG_VTKDCN(vtkDoubleArray, cl, m_Grid, "node_meshdensity_desired");
 555         CadInterface *cad = GuiMainWindow::pointer()->getCadInterface(bcs.first());
 556         vec3_t x1, x2;
 557         m_Grid->GetPoint(id_node1, x1.data());
 558         m_Grid->GetPoint(id_node2, x2.data());
 559         vec3_t x = 0.5*(x1 + x2);
 560         double h = min(cl->GetValue(id_node1), cl->GetValue(id_node2));
 561         vec3_t xs = cad->snapToEdge(x);
 562         if ((xs - x).abs() < 0.2*h) {
 563           is_feature_edge = true;
 564         }
 565       }
 566     }
 567   }
 568   return is_feature_edge;
 569 }
 570
 571 bool MeshPartition::isConvexNode(vtkIdType id_node)
 572 {
 573   int n_faces = n2cGSize(id_node);
 574   if (n_faces <= 1) {
 575     return false;
 576   }
 577
 578   vec3_t x1, cell_centers(0,0,0), cell_normals(0,0,0);
 579   m_Grid->GetPoint(id_node, x1.data());
 580   for (int i = 0; i < n_faces; ++i) {
 581     vtkIdType id_cell = n2cGG(id_node, i);
 582     cell_centers += cellCentre(m_Grid, id_cell);
 583     cell_normals += cellNormal(m_Grid, id_cell);
 584   }
 585   cell_centers *= 1.0/n_faces;
 586   if ((x1 - cell_centers)*cell_normals.normalise() > 0) {
 587     return true;
 588   }
 589   return false;
 590 }
 591
 592
 593 bool MeshPartition::isConvexNode(vtkIdType id_node, QVector<int> bl_codes)
 594 {
 595   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
 596   int n_faces = n2cGSize(id_node);
 597   if (n_faces <= 1) {
 598     return false;
 599   }
 600
 601   vec3_t x1, cell_centers(0,0,0), cell_normals(0,0,0);
 602   m_Grid->GetPoint(id_node, x1.data());
 603   int n_used = 0;
 604   for (int i = 0; i < n_faces; ++i) {
 605     vtkIdType id_cell = n2cGG(id_node, i);
 606     if (bl_codes.contains(cell_code->GetValue(id_cell))) {
 607       cell_centers += cellCentre(m_Grid, id_cell);
 608       cell_normals += cellNormal(m_Grid, id_cell);
 609       ++n_used;
 610     }
 611   }
 612   if (n_used > 0) {
 613     cell_centers *= 1.0/n_used;
 614     if ((x1 - cell_centers)*cell_normals.normalise() > 0) {
 615       return true;
 616     }
 617   }
 618   return false;
 619 }