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     cout << "  tol=" << tol << endl;
 195     EG_VTKSP(vtkUnstructuredGrid, new_grid);
 196     EG_VTKSP(vtkKdTreePointLocator,loc);
 197     loc->SetDataSet(m_Grid);
 198     loc->BuildLocator();
 199     QVector<vtkIdType> pnode2node(part.m_Grid->GetNumberOfPoints());
 200     vtkIdType N = m_Grid->GetNumberOfPoints();
 201     Timer T(10);
 202     for (vtkIdType id_pnode = 0; id_pnode < part.m_Grid->GetNumberOfPoints(); ++id_pnode) {
 203       vec3_t xp, x;
 204       part.m_Grid->GetPoint(id_pnode, xp.data());
 205       vtkIdType id_node = loc->FindClosestPoint(xp.data());
 206       m_Grid->GetPoint(id_node, x.data());
 207       if ((x - xp).abs() < tol) {
 208         pnode2node[id_pnode] = id_node;
 209       } else {
 210         pnode2node[id_pnode] = N;
 211         ++N;
 212       }
 213     }
 214     allocateGrid(new_grid, m_Grid->GetNumberOfCells() + part.m_Cells.size(), N);
 215     for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 216       vec3_t x;
 217       m_Grid->GetPoint(id_node, x.data());
 218       new_grid->GetPoints()->SetPoint(id_node, x.data());
 219       copyNodeData(m_Grid, id_node, new_grid, id_node);
 220     }
 221     QVector<vtkIdType> part_nodes;
 222     getNodesFromCells(part.m_Cells, part_nodes, part.m_Grid);
 223     foreach (vtkIdType id_pnode, part_nodes) {
 224       vec3_t x;
 225       part.m_Grid->GetPoint(id_pnode, x.data());
 226       new_grid->GetPoints()->SetPoint(pnode2node[id_pnode], x.data());
 227       copyNodeData(part.m_Grid, id_pnode, new_grid, pnode2node[id_pnode]);
 228     }
 229     QList<vtkIdType> new_cells;
 230     for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 231       vtkIdType id_new_cell = copyCell(m_Grid, id_cell, new_grid);
 232       copyCellData(m_Grid, id_cell, new_grid, id_new_cell);
 233       new_cells.append(id_new_cell);
 234     }
 235     foreach (vtkIdType id_pcell, part.m_Cells) {
 236       vtkIdType id_new_cell = copyCell(part.m_Grid, id_pcell, new_grid, pnode2node);
 237       copyCellData(part.m_Grid, id_pcell, new_grid, id_new_cell);
 238       new_cells.append(id_new_cell);
 239     }
 240     makeCopy(new_grid, m_Grid);
 241   }
 242 }
 243
 244 double MeshPartition::getSmallestEdgeLength() const
 245 {
 246   double L = 1e99;
 247   foreach (vtkIdType id_cell, m_Cells) {
 248     vtkIdType type_cell = m_Grid->GetCellType(id_cell);
 249     vtkIdType N_pts, *pts;
 250     m_Grid->GetCellPoints(id_cell, N_pts, pts);
 251     QVector<vec3_t> x(N_pts);
 252     for (int i = 0; i < N_pts; ++i) {
 253       m_Grid->GetPoint(pts[i], x[i].data());
 254     }
 255     if        (type_cell == VTK_TRIANGLE) {
 256       L = min(L, (x[0] - x[1]).abs());
 257       L = min(L, (x[1] - x[2]).abs());
 258       L = min(L, (x[2] - x[0]).abs());
 259     } else if (type_cell == VTK_QUAD) {
 260       L = min(L, (x[0] - x[1]).abs());
 261       L = min(L, (x[1] - x[2]).abs());
 262       L = min(L, (x[2] - x[3]).abs());
 263       L = min(L, (x[3] - x[0]).abs());
 264     } else if (type_cell == VTK_TETRA) {
 265       L = min(L, (x[0] - x[1]).abs());
 266       L = min(L, (x[1] - x[2]).abs());
 267       L = min(L, (x[2] - x[0]).abs());
 268       L = min(L, (x[3] - x[0]).abs());
 269       L = min(L, (x[3] - x[1]).abs());
 270       L = min(L, (x[3] - x[2]).abs());
 271     } else if (type_cell == VTK_PYRAMID) {
 272       L = min(L, (x[0] - x[1]).abs());
 273       L = min(L, (x[1] - x[2]).abs());
 274       L = min(L, (x[2] - x[3]).abs());
 275       L = min(L, (x[3] - x[0]).abs());
 276       L = min(L, (x[4] - x[0]).abs());
 277       L = min(L, (x[4] - x[1]).abs());
 278       L = min(L, (x[4] - x[2]).abs());
 279       L = min(L, (x[4] - x[3]).abs());
 280     } else if (type_cell == VTK_WEDGE) {
 281       L = min(L, (x[0] - x[1]).abs());
 282       L = min(L, (x[1] - x[2]).abs());
 283       L = min(L, (x[2] - x[0]).abs());
 284       L = min(L, (x[3] - x[4]).abs());
 285       L = min(L, (x[4] - x[5]).abs());
 286       L = min(L, (x[5] - x[3]).abs());
 287       L = min(L, (x[0] - x[3]).abs());
 288       L = min(L, (x[1] - x[4]).abs());
 289       L = min(L, (x[2] - x[5]).abs());
 290     } else if (type_cell == VTK_HEXAHEDRON) {
 291       L = min(L, (x[0] - x[1]).abs());
 292       L = min(L, (x[1] - x[2]).abs());
 293       L = min(L, (x[2] - x[3]).abs());
 294       L = min(L, (x[3] - x[0]).abs());
 295       L = min(L, (x[4] - x[5]).abs());
 296       L = min(L, (x[5] - x[6]).abs());
 297       L = min(L, (x[6] - x[7]).abs());
 298       L = min(L, (x[7] - x[4]).abs());
 299       L = min(L, (x[0] - x[4]).abs());
 300       L = min(L, (x[1] - x[5]).abs());
 301       L = min(L, (x[2] - x[6]).abs());
 302       L = min(L, (x[3] - x[7]).abs());
 303     }
 304   }
 305   return L;
 306 }
 307
 308 bool MeshPartition::hasNeighNode(vtkIdType id_node, vtkIdType id_neigh)
 309 {
 310   for (int i = 0; i < n2nGSize(id_node); ++i) {
 311     if (n2nGG(id_node, i) == id_neigh) {
 312       return true;
 313     }
 314   }
 315   return false;
 316 }
 317
 318 void MeshPartition::createNodeToBC()
 319 {
 320   EG_VTKDCC(vtkIntArray, cell_code,   m_Grid, "cell_code");
 321   m_N2BC.resize(m_Nodes.size());
 322   for (int i_node = 0; i_node < m_Nodes.size(); ++i_node) {
 323     QSet<int> bcs;
 324     for (int j = 0; j < n2cLSize(i_node); ++j) {
 325       vtkIdType id_cell = n2cLG(i_node, j);
 326       if (isSurface(id_cell, m_Grid)) {
 327         int bc = cell_code->GetValue(n2cLG(i_node, j));
 328         if (bc != 0) {
 329           bcs.insert(cell_code->GetValue(n2cLG(i_node, j)));
 330         }
 331       }
 332     }
 333     m_N2BC[i_node].resize(bcs.size());
 334     int i = 0;
 335     foreach (int bc, bcs) {
 336       m_N2BC[i_node][i++] = bc;
 337     }
 338   }
 339 }
 340
 341 bool MeshPartition::hasBC(vtkIdType id_node, int bc)
 342 {
 343   bool found = false;
 344   for (int j = 0; j < n2bcGSize(id_node); ++j) {
 345     if (n2bcG(id_node, j) == bc) {
 346       found = true;
 347       break;
 348     }
 349   }
 350   return found;
 351 }
 352
 353 vtkIdType MeshPartition::getVolumeCell(vtkIdType id_face)
 354 {
 355   checkLNodes();
 356   checkLCells();
 357   checkN2C();
 358   return findVolumeCell(m_Grid, id_face, m_LNodes, m_Cells, m_LCells, m_N2C);
 359 }
 360
 361 vec3_t MeshPartition::globalNormal(vtkIdType id_node)
 362 {
 363   vec3_t normal(0,0,0);
 364   for (int i = 0; i < n2cGSize(id_node); ++i) {
 365     vtkIdType id_cell = n2cGG(id_node, i);
 366     if (isSurface(id_cell, m_Grid)) {
 367       vtkIdType N_pts, *pts;
 368       m_Grid->GetCellPoints(id_cell, N_pts, pts);
 369       vec3_t a, b, c;
 370       for (int j = 0; j < N_pts; ++j) {
 371         if (pts[j] == id_node) {
 372           m_Grid->GetPoint(pts[j], a.data());
 373           if (j > 0) {
 374             m_Grid->GetPoint(pts[j-1], b.data());
 375           } else {
 376             m_Grid->GetPoint(pts[N_pts-1], b.data());
 377           }
 378           if (j < N_pts - 1) {
 379             m_Grid->GetPoint(pts[j+1], c.data());
 380           } else {
 381             m_Grid->GetPoint(pts[0], c.data());
 382           }
 383         }
 384       }
 385       vec3_t u = b - a;
 386       vec3_t v = c - a;
 387       double alpha = GeometryTools::angle(u, v);
 388       vec3_t n = u.cross(v);
 389       n.normalise();
 390       if (checkVector(n)) {
 391         normal -= alpha*n;
 392       }
 393     }
 394   }
 395   normal.normalise();
 396   return normal;
 397 }
 398
 399 double MeshPartition::getAverageSurfaceEdgeLength(vtkIdType id_node)
 400 {
 401   QSet<vtkIdType> surface_neighbours;
 402   for (int i = 0; i < n2cGSize(id_node); ++i) {
 403     vtkIdType id_cell = n2cGG(id_node, i);
 404     if (isSurface(id_cell, m_Grid)) {
 405       EG_GET_CELL(id_cell, m_Grid);
 406       for (int j = 0; j < num_pts; ++j) {
 407         if (pts[j] != id_node) {
 408           surface_neighbours.insert(pts[j]);
 409         }
 410       }
 411     }
 412   }
 413   double L = 0;
 414   if (surface_neighbours.size() > 0) {
 415     vec3_t x, xn;
 416     m_Grid->GetPoint(id_node, x.data());
 417     foreach (vtkIdType id_neigh, surface_neighbours) {
 418       m_Grid->GetPoint(id_neigh, xn.data());
 419       L += (x - xn).abs();
 420     }
 421     L /= surface_neighbours.size();
 422   }
 423   return L;
 424 }
 425
 426 void MeshPartition::computeMinAndMaxSurfaceStencilEdgeLengths(vtkIdType id_node, double &l_min, double &l_max)
 427 {
 428   l_min = 1e99;
 429   l_max = 0;
 430   for (int i = 0; i < n2cGSize(id_node); ++i) {
 431     vtkIdType id_cell = n2cGG(id_node, i);
 432     if (isSurface(id_cell, m_Grid)) {
 433       vtkIdType *pts, num_pts;
 434       m_Grid->GetCellPoints(id_cell, num_pts, pts);
 435       vec3_t x1, x2;
 436       m_Grid->GetPoint(pts[num_pts - 1], x1.data());
 437       for (int j = 0; j < num_pts; ++j) {
 438         m_Grid->GetPoint(pts[j], x2.data());
 439         double L = (x1 - x2).abs();
 440         l_min = min(L, l_min);
 441         l_max = max(L, l_max);
 442         x1 = x2;
 443       }
 444     }
 445   }
 446 }
 447
 448 double MeshPartition::getMinSurfaceStencilEdgeLength(vtkIdType id_node)
 449 {
 450   double l_min, l_max;
 451   computeMinAndMaxSurfaceStencilEdgeLengths(id_node, l_min, l_max);
 452   return l_min;
 453 }
 454
 455 double MeshPartition::getMaxSurfaceStencilEdgeLength(vtkIdType id_node)
 456 {
 457   double l_min, l_max;
 458   computeMinAndMaxSurfaceStencilEdgeLengths(id_node, l_min, l_max);
 459   return l_max;
 460 }
 461
 462 int MeshPartition::getNumberOfFeatureNeighbours(vtkIdType id_node)
 463 {
 464   EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type");
 465   int N = 0;
 466   for (int i = 0; i < n2nGSize(id_node); ++i) {
 467     char type = node_type->GetValue(n2nGG(id_node, i));
 468     if (type == EG_FEATURE_EDGE_VERTEX || type == EG_FEATURE_CORNER_VERTEX) {
 469       ++N;
 470     }
 471   }
 472   return N;
 473 }
 474
 475 int MeshPartition::getEdgeType(vtkIdType id_node1, vtkIdType id_node2)
 476 {
 477   QList <vtkIdType> edge_faces;
 478   getEdgeFaces(id_node1, id_node2, edge_faces);
 479   int edge = 0;
 480   if (edge_faces.size() == 1) {
 481     edge = 2;
 482   } else if (edge_faces.size() >= 2) {
 483     EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
 484     QSet<int> bcs;
 485     foreach (vtkIdType id_face, edge_faces) {
 486       bcs.insert(cell_code->GetValue(id_face));
 487     }
 488     if (bcs.size() > 1) {
 489       edge = 2;
 490     } else {
 491       edge = 1;
 492     }
 493   }
 494   return edge;
 495 }
 496
 497 int MeshPartition::computeTopoDistance(vtkIdType id_node1, vtkIdType id_node2, int max_dist, int restriction_type)
 498 {
 499   int dist = 0;
 500   QSet<vtkIdType> candidates;
 501   candidates.insert(id_node1);
 502   while (dist < max_dist && !candidates.contains(id_node2)) {
 503     foreach (vtkIdType id_node, candidates) {
 504       for (int i = 0; i < n2nGSize(id_node); ++i) {
 505         vtkIdType id_neigh = n2nGG(id_node,i);
 506         if (!candidates.contains(id_neigh)) {
 507           if (getEdgeType(id_node, id_neigh) >= restriction_type) {
 508             candidates.insert(id_neigh);
 509           }
 510         }
 511       }
 512     }
 513     ++dist;
 514   }
 515   return dist;
 516 }
 517
 518 void MeshPartition::getCommonNodes(vtkIdType id_cell1, vtkIdType id_cell2, QVector<vtkIdType> &common_nodes)
 519 {
 520   common_nodes.clear();
 521   QSet<vtkIdType> nodes1, nodes2;
 522   vtkIdType num_pts, *pts;
 523   m_Grid->GetCellPoints(id_cell1, num_pts, pts);
 524   for (int i = 0; i < num_pts; ++i) {
 525     nodes1.insert(pts[i]);
 526   }
 527   m_Grid->GetCellPoints(id_cell2, num_pts, pts);
 528   for (int i = 0; i < num_pts; ++i) {
 529     nodes2.insert(pts[i]);
 530   }
 531   nodes1.intersect(nodes2);
 532   common_nodes.resize(nodes1.size());
 533   qCopy(nodes1.begin(), nodes1.end(), common_nodes.begin());
 534 }
 535
 536 bool MeshPartition::isFeatureEdge(vtkIdType id_node1, vtkIdType id_node2, double feature_angle)
 537 {
 538   bool is_feature_edge = false;
 539   QVector<int> bcs;
 540   QVector<vtkIdType> nodes(2);
 541   nodes[0] = id_node1;
 542   nodes[1] = id_node2;
 543   getCommonBcs(nodes, bcs);
 544   if (bcs.size() == 1) {
 545     QList<vtkIdType> faces;
 546     getEdgeFaces(id_node1, id_node2, faces);
 547     if (faces.size() == 2) {
 548       vec3_t n1 = cellNormal(m_Grid, faces[0]);
 549       vec3_t n2 = cellNormal(m_Grid, faces[1]);
 550       if (GeometryTools::angle(n1, n2) > feature_angle) {
 551         EG_VTKDCN(vtkDoubleArray, cl, m_Grid, "node_meshdensity_desired");
 552         CadInterface *cad = GuiMainWindow::pointer()->getCadInterface(bcs.first());
 553         vec3_t x1, x2;
 554         m_Grid->GetPoint(id_node1, x1.data());
 555         m_Grid->GetPoint(id_node2, x2.data());
 556         vec3_t x = 0.5*(x1 + x2);
 557         double h = min(cl->GetValue(id_node1), cl->GetValue(id_node2));
 558         vec3_t xs = cad->snapToEdge(x);
 559         if ((xs - x).abs() < 0.2*h) {
 560           is_feature_edge = true;
 561         }
 562       }
 563     }
 564   }
 565   return is_feature_edge;
 566 }