src/libengrid/removepoints.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 #include "removepoints.h"
  22
  23 #include "checksurfaceintegrity.h"
  24
  25 #include "geometrytools.h"
  26 using namespace GeometryTools;
  27
  28 #include <cmath>
  29 using namespace std;
  30
  31 #include <iostream>
  32 using namespace std;
  33
  34 RemovePoints::RemovePoints() : SurfaceOperation() {
  35   setQuickSave(true);
  36   getSet("surface meshing", "point removal threshold", 2, m_Threshold);
  37   m_ProtectFeatureEdges = false;
  38   m_PerformGeometricChecks = true;
  39   m_UpdatePSP = false;
  40 }
  41
  42 void RemovePoints::markFeatureEdges()
  43 {
  44   EG_VTKDCN(vtkDoubleArray, characteristic_length_desired, m_Grid, "node_meshdensity_desired");
  45   m_IsFeatureNode.fill(false, m_Part.getNumberOfNodes());
  46   if(m_ProtectFeatureEdges) {
  47     for (int i_nodes = 0; i_nodes < m_Part.getNumberOfNodes(); ++i_nodes) {
  48       if (!m_IsFeatureNode[i_nodes]) {
  49         vtkIdType id_node1 = m_Part.globalNode(i_nodes);
  50         for (int j = 0; j < m_Part.n2nLSize(i_nodes); ++j) {
  51           vtkIdType id_node2 = m_Part.n2nLG(i_nodes, j);
  52           QSet<vtkIdType> edge_cells;
  53           int N = getEdgeCells(id_node1, id_node2, edge_cells);
  54           if (N != 2) {
  55             m_IsFeatureNode[i_nodes] = true;
  56             m_IsFeatureNode[m_Part.localNode(id_node2)] = true;
  57           } else {
  58             QSet<vtkIdType>::iterator iter = edge_cells.begin();
  59             vtkIdType id_cell1 = *iter;
  60             ++iter;
  61             vtkIdType id_cell2 = *iter;
  62             vec3_t n1 = cellNormal(m_Grid, id_cell1);
  63             vec3_t n2 = cellNormal(m_Grid, id_cell2);
  64             if (angle(n1, n2) >= m_FeatureAngle) {
  65               m_IsFeatureNode[i_nodes] = true;
  66               m_IsFeatureNode[m_Part.localNode(id_node2)] = true;
  67             }
  68           }
  69         }
  70       }
  71     }
  72     int N = 0;
  73     for (int i = 0; i < m_IsFeatureNode.size(); ++i) {
  74       if(m_IsFeatureNode[i]) {
  75         ++N;
  76       }
  77     }
  78     cout << N << " nodes on feature edges (angle >= " << GeometryTools::rad2deg(m_FeatureAngle) << "deg)" << endl;
  79   }
  80 }
  81
  82 void RemovePoints::fixNodes(const QVector<bool> &fixnodes)
  83 {
  84   if (fixnodes.size() != m_Grid->GetNumberOfPoints()) {
  85     EG_BUG;
  86   }
  87   m_Fixed = fixnodes;
  88 }
  89
  90 bool RemovePoints::checkEdge(vtkIdType id_node1, vtkIdType id_node2)
  91 {
  92   EG_VTKDCN(vtkDoubleArray, cl, m_Grid, "node_meshdensity_desired");
  93   double cl1 = cl->GetValue(id_node1);
  94   double cl2 = cl->GetValue(id_node2);
  95   vec3_t x1, x2;
  96   m_Grid->GetPoint(id_node1, x1.data());
  97   m_Grid->GetPoint(id_node2, x2.data());
  98   double L = (x1 - x2).abs();
  99   double cl_crit = max(cl1, cl2) / m_Threshold;
 100   if (L < cl_crit) {
 101     return true;
 102   }
 103   return false;
 104 }
 105
 106 void RemovePoints::operate()
 107 {
 108   if (m_Fixed.size() != m_Grid->GetNumberOfPoints()) {
 109     m_Fixed.fill(false, m_Grid->GetNumberOfPoints());
 110   }
 111
 112   MeshPartition full_partition(m_Grid, true);
 113   full_partition.setAllCells();
 114
 115   int N1 = m_Grid->GetNumberOfPoints();
 116
 117   QVector<vtkIdType> selected_cells;
 118   getSurfaceCells(m_BoundaryCodes, selected_cells, m_Grid);
 119   QVector<vtkIdType> selected_nodes;
 120   getNodesFromCells(selected_cells, selected_nodes, m_Grid);
 121
 122   setAllSurfaceCells();
 123   computeNormals();
 124   updateNodeInfo();
 125
 126   EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type");
 127   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
 128   EG_VTKDCN(vtkDoubleArray, characteristic_length_desired, m_Grid, "node_meshdensity_desired");
 129
 130   // global values
 131   QVector <vtkIdType> all_deadcells;
 132   QVector <vtkIdType> all_mutatedcells;
 133   int num_newpoints = 0;
 134   int num_newcells = 0;
 135
 136   QVector<bool> marked_nodes(m_Grid->GetNumberOfPoints(), false);
 137
 138   QVector <vtkIdType> deadnode_vector;
 139   QVector <vtkIdType> snappoint_vector;
 140
 141   foreach (vtkIdType id_node, selected_nodes) {
 142     if (node_type->GetValue(id_node) != EG_FIXED_VERTEX && !m_Fixed[id_node]) {
 143
 144       // check that node is only surrounded by triangles
 145       bool tri_only = true;
 146       for (int j = 0; j < m_Part.n2cGSize(id_node); ++j) {
 147         vtkIdType id_cell = m_Part.n2cGG(id_node, j);
 148         if(m_Grid->GetCellType(id_cell) != VTK_TRIANGLE) {
 149           tri_only = false;
 150           break;
 151         }
 152       }
 153
 154       if (!marked_nodes[id_node] && tri_only) {
 155         for (int j = 0; j < m_Part.n2nGSize(id_node); ++j) {
 156           vtkIdType id_neigh = m_Part.n2nGG(id_node, j);
 157           if (checkEdge(id_node, id_neigh)) {
 158             QVector <vtkIdType> dead_cells;
 159             QVector <vtkIdType> mutated_cells;
 160             int l_num_newpoints = 0;
 161             int l_num_newcells = 0;
 162             if (isSnapPoint(id_node, id_neigh, dead_cells, mutated_cells, l_num_newpoints, l_num_newcells, marked_nodes)) {
 163
 164               // add deadnode/snappoint pair
 165               deadnode_vector.push_back(id_node);
 166               snappoint_vector.push_back(id_neigh);
 167
 168               // update global values
 169               num_newpoints += l_num_newpoints;
 170               num_newcells  += l_num_newcells;
 171               all_deadcells += dead_cells;
 172               all_mutatedcells += mutated_cells;
 173
 174               // mark neighbour nodes
 175               for (int k = 0; k < m_Part.n2nGSize(id_node); ++k) {
 176                 marked_nodes[m_Part.n2nGG(id_node, k)] = true;
 177               }
 178               for (int k = 0; k < m_Part.n2nGSize(id_neigh); ++k) {
 179                 marked_nodes[m_Part.n2nGG(id_neigh, k)] = true;
 180               }
 181
 182               break;
 183             }
 184           }
 185         }
 186       }
 187     }
 188   }
 189
 190   //delete
 191   if (num_newpoints != -deadnode_vector.size()) {
 192     EG_BUG;
 193   }
 194   if (num_newcells != -all_deadcells.size()) {
 195     EG_BUG;
 196   }
 197
 198   deleteSetOfPoints(deadnode_vector, snappoint_vector, all_deadcells, all_mutatedcells);
 199
 200   int N2 = m_Grid->GetNumberOfPoints();
 201   m_NumRemoved = N1 - N2;
 202 }
 203
 204 /// \todo finish this function and optimize it.
 205 bool RemovePoints::checkForDestroyedVolumes(vtkIdType id_node1, vtkIdType id_node2, int& N_common_points)
 206 {
 207   if (id_node1 == id_node2) {
 208     EG_BUG;
 209   }
 210
 211   l2l_t  n2n   = getPartN2N();
 212   g2l_t _nodes = getPartLocalNodes();
 213   l2g_t nodes  = getPartNodes();
 214
 215   QVector<int> node1_neighbours = n2n[_nodes[id_node1]];
 216   QVector<int> node2_neighbours = n2n[_nodes[id_node2]];
 217   QVector<int> common_points;
 218   qcontIntersection(node1_neighbours, node2_neighbours, common_points);
 219   // set N_common_points
 220   N_common_points = common_points.size();
 221
 222   // TEST 0: TOPOLOGICAL: DeadNode, PSP and any common point must belong to a cell.
 223   for(int i = 0; i < N_common_points; i++) {
 224     int i_common_point_1 = common_points[i];
 225     vtkIdType id_common_point_1 = nodes[i_common_point_1];
 226     if(!isCell(id_node1, id_node2, id_common_point_1)) {
 227       if(DebugLevel > 100) {
 228         qDebug() << "test 0 failed";
 229         qDebug() << "id_node1, id_node2, id_common_point_1=" << id_node1 << id_node2 << id_common_point_1;
 230       }
 231       return true;
 232     }
 233     // TEST 1: TOPOLOGICAL: Moving DeadNode to PSP must not lay any cell on another cell.
 234     //                      => For any pair of common points (cp1,cp2), (cp1,cp2,DeadNode)+(cp1,cp2,PSP)
 235     //                         must not be cells at the same time!
 236     for(int j = i + 1; j < common_points.size(); j++) {
 237       int i_common_point_2 = common_points[j];
 238       vtkIdType id_common_point_2 = nodes[i_common_point_2];
 239       if(isCell(id_common_point_1, id_common_point_2, id_node1) && isCell(id_common_point_1, id_common_point_2, id_node2)) {
 240         if(DebugLevel > 100) {
 241           qDebug() << "test 1 failed";
 242           qDebug() << "id_common_point_1, id_common_point_2, id_node1=" << id_common_point_1 << id_common_point_2 << id_node1;
 243           qDebug() << "id_common_point_1, id_common_point_2, id_node2=" << id_common_point_1 << id_common_point_2 << id_node2;
 244         }
 245         return true;
 246       }
 247     }
 248   }
 249
 250   QSet<vtkIdType> all_faces;
 251   for (int i = 0; i < m_Part.n2nGSize(id_node1); ++i) {
 252     for (int j = 0; j < m_Part.n2cGSize(m_Part.n2nGG(id_node1, i)); ++j) {
 253       all_faces.insert(m_Part.n2cGG(m_Part.n2nGG(id_node1, i), j));
 254     }
 255   }
 256   for (int i = 0; i < m_Part.n2nGSize(id_node2); ++i) {
 257     for (int j = 0; j < m_Part.n2cGSize(m_Part.n2nGG(id_node2, i)); ++j) {
 258       all_faces.insert(m_Part.n2cGG(m_Part.n2nGG(id_node2, i), j));
 259     }
 260   }
 261   QSet<vtkIdType> near_faces;
 262   for (int i = 0; i < m_Part.n2cGSize(id_node1); ++i) {
 263     near_faces.insert(m_Part.n2cGG(id_node1, i));
 264   }
 265   for (int i = 0; i < m_Part.n2cGSize(id_node2); ++i) {
 266     near_faces.insert(m_Part.n2cGG(id_node2, i));
 267   }
 268   QSet<vtkIdType> far_faces = all_faces - near_faces;
 269   bool tetra = true;
 270   foreach (vtkIdType id_cell, far_faces) {
 271     vtkIdType N_pts, *pts;
 272     m_Grid->GetCellPoints(id_cell, N_pts, pts);
 273     for (int i = 0; i < N_pts; ++i) {
 274       if (!m_Part.hasNeighNode(pts[i], id_node1) && !m_Part.hasNeighNode(pts[i], id_node2)) {
 275         tetra = false;
 276         break;
 277       }
 278     }
 279     if (!tetra) {
 280       break;
 281     }
 282   }
 283   if (tetra) {
 284     return true;
 285   }
 286
 287
 288
 289   //FIX THIS!!!!
 290   /*
 291   // check if DeadNode, PSP and common points form a tetrahedron.
 292   if ( n2n[_nodes[intersection1]].contains( _nodes[intersection2] ) ) { //if there's an edge between intersection1 and intersection2
 293     //check if (node1,intersection1,intersection2) and (node2,intersection1,intersection2) are defined as cells!
 294     QVector<int> S1 = n2c[_nodes[intersection1]];
 295     QVector<int> S2 = n2c[_nodes[intersection2]];
 296     QVector<int> Si;
 297     qcontIntersection( S1, S2, Si );
 298     int counter = 0;
 299     foreach( int i_cell, Si ) {
 300       vtkIdType num_pts, *pts;
 301       m_Grid->GetCellPoints( cells[i_cell], num_pts, pts );
 302       for ( int i = 0; i < num_pts; ++i ) {
 303         if ( pts[i] == id_node1 || pts[i] == id_node2 ) counter++;
 304       }
 305     }
 306     if ( counter >= 2 ) {
 307       IsTetra = true;
 308     }
 309   }
 310 */
 311   return false;
 312 }
 313
 314 int RemovePoints::numberOfCommonPoints(vtkIdType id_node1, vtkIdType id_node2, bool& IsTetra)
 315 {
 316   l2l_t  n2n   = getPartN2N();
 317   l2l_t  n2c   = getPartN2C();
 318   g2l_t _nodes = getPartLocalNodes();
 319   l2g_t nodes  = getPartNodes();
 320   l2g_t cells = getPartCells();
 321
 322   QVector<int> node1_neighbours = n2n[_nodes[id_node1]];
 323   QVector<int> node2_neighbours = n2n[_nodes[id_node2]];
 324   QVector<int> intersection;
 325   qcontIntersection(node1_neighbours, node2_neighbours, intersection);
 326   int N = intersection.size();
 327   IsTetra = false;
 328   if(N == 2) {
 329     vtkIdType intersection1 = nodes[intersection[0]];
 330     vtkIdType intersection2 = nodes[intersection[1]];
 331
 332     // test if id_node1, id_node2 and intersection* form a cell
 333     QVector <vtkIdType> EdgeCells_1i;
 334     QVector <vtkIdType> EdgeCells_2i;
 335     QVector <vtkIdType> inter;
 336     int Ncells;
 337
 338     // intersection1
 339     Ncells = getEdgeCells(id_node1, intersection1, EdgeCells_1i);
 340     if(N != 2) {
 341       qWarning() << "Ncells=" << Ncells;
 342       EG_BUG;
 343     }
 344     Ncells = getEdgeCells(id_node2, intersection1, EdgeCells_2i);
 345     if(Ncells != 2) {
 346       qWarning() << "Ncells=" << Ncells;
 347       EG_BUG;
 348     }
 349     qcontIntersection(EdgeCells_1i, EdgeCells_2i, inter);
 350     if(inter.size() <= 0) EG_BUG;   // (id_node1, id_node2, intersection1) is not a cell
 351
 352     // intersection2
 353     Ncells = getEdgeCells(id_node1, intersection2, EdgeCells_1i);
 354     if(Ncells != 2) {
 355       qWarning() << "Ncells=" << Ncells;
 356       EG_BUG;
 357     }
 358     Ncells = getEdgeCells(id_node2, intersection2, EdgeCells_2i);
 359     if(Ncells != 2) {
 360       qWarning() << "Ncells=" << Ncells;
 361       EG_BUG;
 362     }
 363     qcontIntersection(EdgeCells_1i, EdgeCells_2i, inter);
 364     if(inter.size() <= 0) EG_BUG;   // (id_node1, id_node2, intersection2) is not a cell
 365
 366     // check if DeadNode, PSP and common points form a tetrahedron.
 367     if(n2n[_nodes[intersection1]].contains(_nodes[intersection2])) {      //if there's an edge between intersection1 and intersection2
 368       //check if (node1,intersection1,intersection2) and (node2,intersection1,intersection2) are defined as cells!
 369       QVector<int> S1 = n2c[_nodes[intersection1]];
 370       QVector<int> S2 = n2c[_nodes[intersection2]];
 371       QVector<int> Si;
 372       qcontIntersection(S1, S2, Si);
 373       int counter = 0;
 374       foreach(int i_cell, Si) {
 375         vtkIdType num_pts, *pts;
 376         m_Grid->GetCellPoints(cells[i_cell], num_pts, pts);
 377         for(int i = 0; i < num_pts; ++i) {
 378           if(pts[i] == id_node1 || pts[i] == id_node2) counter++;
 379         }
 380       }
 381       if(counter >= 2) {
 382         IsTetra = true;
 383       }
 384     }
 385   }
 386   return(N);
 387 }
 388
 389 bool RemovePoints::flippedCell2(vtkIdType id_node, vec3_t x_new) {
 390   /*
 391   for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
 392
 393     vtkIdType id_cell = m_Part.n2cGG(id_node, i);
 394
 395     vtkIdType N_pts, *pts;
 396     m_Grid->GetCellPoints(id_cell, N_pts, pts);
 397     if(N_pts!=3) EG_BUG;
 398
 399     int i_pts=0;
 400     for(i_pts=0; i_pts<N_pts; i_pts++) {
 401       if(pts[i_pts]==id_node) break;
 402     }
 403     if(pts[i_pts]!=id_node) EG_BUG;
 404
 405     vec3_t x1, x2, x_old;
 406     m_Grid->GetPoint(pts[(i_pts+1)%N_pts],x1.data());
 407     m_Grid->GetPoint(pts[(i_pts+2)%N_pts],x2.data());
 408
 409     vec3_t old_cell_normal = GeometryTools::triNormal(x_old, x1, x2);
 410     vec3_t new_cell_normal = GeometryTools::triNormal(x_new, x1, x2);
 411
 412     if(old_cell_normal.abs2()==0) EG_BUG;
 413     if(old_cell_normal.abs2()==0) EG_BUG;
 414
 415     GeometryTools::cellNormal(m_Grid, );
 416     cell_normals.normalise();
 417
 418     vtkIdType *pts;
 419     vtkIdType npts;
 420     vec3_t n(0,0,0);
 421     grid->GetCellPoints(i, npts, pts);
 422     if (npts == 3) {
 423       return triNormal(grid,pts[0],pts[1],pts[2]);
 424     } else if (npts == 4) {
 425       return quadNormal(grid,pts[0],pts[1],pts[2],pts[3]);
 426     } else {
 427       EG_BUG;
 428     }
 429     return n;
 430
 431   }
 432   */
 433   return true;
 434 }
 435
 436 /// \todo adapt for multiple volumes
 437 bool RemovePoints::flippedCell(vtkIdType id_node, vec3_t x_new, vtkIdType id_cell) {
 438   if( m_Grid->GetCellType(id_cell) == VTK_WEDGE ) EG_BUG;
 439
 440   vec3_t x_old;
 441   m_Grid->GetPoint(id_node, x_old.data());
 442
 443   vec3_t n(0, 0, 0);
 444   bool move = true;
 445   QVector<vec3_t> cell_normals(m_Part.n2cGSize(id_node));
 446   double A_max = 0;
 447   for(int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
 448     double A = fabs(GeometryTools::cellVA(m_Grid, m_Part.n2cGG(id_node, i)));
 449     A_max = max(A, A_max);
 450     cell_normals[i] = GeometryTools::cellNormal(m_Grid, m_Part.n2cGG(id_node, i));
 451     cell_normals[i].normalise();
 452   }
 453   int N = 0;
 454   for(int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
 455     double A = fabs(GeometryTools::cellVA(m_Grid, m_Part.n2cGG(id_node, i)));
 456     if(A > 0.01 * A_max) {
 457       n += cell_normals[i];
 458       ++N;
 459     }
 460   }
 461   if(N == 0) {
 462     move = false;
 463   } else {
 464     n.normalise();
 465     double L_max = 0;
 466     for(int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
 467       vec3_t xn;
 468       m_Grid->GetPoint(m_Part.n2nGG(id_node, i), xn.data());
 469       double L = (xn - x_old).abs();
 470       L_max = max(L, L_max);
 471     }
 472     vec3_t x_summit = x_old + L_max * n;
 473     vec3_t x[3];
 474     vtkIdType N_pts, *pts;
 475     m_Grid->GetCellPoints(id_cell, N_pts, pts);
 476     if(N_pts != 3) {
 477       EG_BUG;
 478     }
 479     for(int j = 0; j < N_pts; ++j) {
 480       m_Grid->GetPoint(pts[j], x[j].data());
 481     }
 482     if(GeometryTools::tetraVol(x[0], x[1], x[2], x_summit, false) <= 0) {
 483       move = false;
 484     }
 485   }
 486
 487   return !move;
 488 }
 489
 490 /** This function tries to find a valid snappoint for DeadNode and returns its ID if it finds one, otherwise it returns -1.
 491     If a valid snappoint is found, the corresponding dead and mutated cells are returned via DeadCells and MutatedCells.
 492
 493  DEFINITIONS:
 494  Normal cell: nothing has changed
 495  Dead cell: the cell does not exist anymore
 496  Mutated cell: the cell's form has changed
 497
 498  Basic algorithm:\n
 499  foreach(potential snap point of DeadNode) {\n
 500    bool IsValidSnapPoint = true;\n
 501    some tests; if any fails: IsValidSnapPoint = false; continue;\n
 502    // reset output variables\n
 503    num_newpoints = -1;\n
 504    num_newcells = 0;\n
 505    DeadCells.clear();\n
 506    MutatedCells.clear();\n
 507    foreach(neighbour cell of DeadNode) {\n
 508      more tests; if any fails: IsValidSnapPoint = false; continue;\n
 509      fill DeadCells + MutatedCells;\n
 510    }\n
 511    even more tests; if any fails: IsValidSnapPoint = false; continue;\n
 512    if(IsValidSnapPoint) {\n
 513      SnapPoint = PSP;\n
 514      break;\n
 515    }\n
 516  }\n
 517
 518  \todo Clean up this function
 519  */
 520 bool RemovePoints::isSnapPoint(vtkIdType id_node1, vtkIdType id_node2,
 521                                QVector<vtkIdType>& dead_cells, QVector<vtkIdType>& mutated_cells,
 522                                int& num_newpoints, int& num_newcells,
 523                                const QVector<bool>& marked_nodes)
 524 {
 525   EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type");
 526   if (node_type->GetValue(id_node1) == EG_FIXED_VERTEX) {
 527     EG_BUG;
 528   }
 529
 530   QVector <vtkIdType> psp_vector = getPotentialSnapPoints(id_node1);
 531   if (!psp_vector.contains(id_node2)) {
 532     return false;
 533   }
 534
 535   vec3_t x1;
 536   m_Grid->GetPoint(id_node1, x1.data());
 537
 538   // TEST 1 : TOPOLOGICAL : Is the node already marked?
 539   if (marked_nodes[id_node2]) {
 540     return false;
 541   }
 542
 543   // TEST 2: TOPOLOGICAL: do not cut off feature corners
 544   {
 545     QSet<vtkIdType> common_nodes, n2n2;
 546     m_Part.getGlobalN2N(id_node1, common_nodes);
 547     m_Part.getGlobalN2N(id_node2, n2n2);
 548     common_nodes.intersect(n2n2);
 549     foreach (vtkIdType id_neigh, common_nodes) {
 550       if (node_type->GetValue(id_neigh) == EG_FEATURE_CORNER_VERTEX) {
 551         return false;
 552       }
 553     }
 554   }
 555
 556   // TEST 3: TOPOLOGICAL: Moving id_node1 to id_node2 must not lay any cell on another cell.
 557   int num_common_points = 0;
 558   if (checkForDestroyedVolumes(id_node1, id_node2, num_common_points)) {
 559     return false;
 560   }
 561
 562   // TEST 4: normal irregularity
 563   if (normalIrregularity(id_node1) > normalIrregularity(id_node2)) {
 564     //return false;
 565   }
 566
 567   //count number of points and cells to remove + analyse cell transformations
 568   num_newpoints = -1;
 569   num_newcells = 0;
 570   dead_cells.clear();
 571   mutated_cells.clear();
 572   //foreach (int i_cell, n2c[_nodes[DeadNode]]) {   //loop through potentially dead cells
 573   for (int i = 0; i < m_Part.n2cGSize(id_node1); ++i) {
 574     vtkIdType id_cell = m_Part.n2cGG(id_node1, i);
 575
 576     EG_GET_CELL(id_cell, m_Grid);
 577     if (type_cell == VTK_WEDGE) {
 578       EG_BUG;
 579     }
 580     if(num_pts != 3) {
 581       return false;
 582     }
 583
 584     bool contains_snap_point = false;
 585     bool invincible = false; // a point with only one cell is declared invincible.
 586     for (int j = 0; j < num_pts; ++j) {
 587       if (pts[j] == id_node2) {
 588         contains_snap_point = true;
 589       }
 590       if (pts[j] != id_node1 && pts[j] != id_node2 && m_Part.n2cGSize(pts[j]) <= 1) {
 591         invincible = true;
 592       }
 593     }
 594
 595     if (contains_snap_point) {    // potential dead cell
 596       if (invincible) {
 597         // TEST 4: TOPOLOGICAL: Check that empty lines aren't left behind when a cell is killed
 598         return false;
 599       } else {
 600         dead_cells.push_back(id_cell);
 601         num_newcells -= 1;
 602       }
 603     } else { // if the cell does not contain the SnapPoint (potential mutated cell)
 604
 605       vtkIdType old_triangle[3];
 606       vtkIdType new_triangle[3];
 607
 608       for (int j = 0; j < num_pts; ++j) {
 609         old_triangle[j] = pts[j];
 610         new_triangle[j] = ((pts[j] == id_node1) ? id_node2 : pts[j]);
 611       }
 612       vec3_t n_old = triNormal(m_Grid, old_triangle[0], old_triangle[1], old_triangle[2]);
 613       vec3_t n_new = triNormal(m_Grid, new_triangle[0], new_triangle[1], new_triangle[2]);
 614       double A_old = n_old.abs();
 615       double A_new = n_new.abs();
 616       n_old.normalise();
 617       n_new.normalise();
 618
 619       // TEST 5: GEOMETRICAL: area + inversion check
 620       if (m_PerformGeometricChecks) {
 621         //if(Old_N * New_N < 0.1 || New_N * New_N < Old_N * Old_N * 1. / 100.) {
 622         if (n_old * n_new < 0.2 || A_new < 0.1*A_old) {
 623           return false;
 624         }
 625       }
 626
 627       mutated_cells.push_back(id_cell);
 628     } // end of if the cell does not contain the SnapPoint (potential mutated cell)
 629   }
 630
 631   // TEST 6: TOPOLOGICAL: survivor check
 632   if (m_Grid->GetNumberOfCells() + num_newcells <= 0) {
 633     return false;
 634   }
 635
 636   return true;
 637 }
 638
 639 bool RemovePoints::deleteSetOfPoints(const QVector<vtkIdType>& deadnode_vector,
 640                                      const QVector<vtkIdType>& snappoint_vector,
 641                                      const QVector<vtkIdType>& all_deadcells,
 642                                      const QVector<vtkIdType>& all_mutatedcells) {
 643   int initial_num_points = m_Grid->GetNumberOfPoints();
 644
 645   //src grid info
 646   int num_points = m_Grid->GetNumberOfPoints();
 647   int num_cells = m_Grid->GetNumberOfCells();
 648
 649   int num_newcells = -all_deadcells.size();
 650   int num_newpoints = -deadnode_vector.size();
 651
 652   //  if ( num_newcells != 2*num_newpoints ) {
 653   //    EG_BUG;
 654   //  }
 655
 656   //allocate
 657   EG_VTKSP(vtkUnstructuredGrid, dst);
 658   allocateGrid(dst, num_cells + num_newcells, num_points + num_newpoints);
 659
 660   //vector used to redefine the new point IDs
 661   QVector <vtkIdType> OffSet(num_points);
 662
 663   //copy undead points
 664   QVector<bool> is_deadnode(m_Grid->GetNumberOfPoints(), false);
 665   QVector<int> glob2dead(m_Grid->GetNumberOfPoints(), -1);
 666   for(int i_deadnodes = 0; i_deadnodes < deadnode_vector.size(); ++i_deadnodes) {
 667     vtkIdType id_node = deadnode_vector[i_deadnodes];
 668     if(id_node > m_Grid->GetNumberOfPoints()) {
 669       EG_BUG;
 670     }
 671     is_deadnode[id_node] = true;
 672     glob2dead[id_node] = i_deadnodes;
 673   }
 674   vtkIdType dst_id_node = 0;
 675   for(vtkIdType src_id_node = 0; src_id_node < num_points; ++src_id_node) { //loop through src points
 676     OffSet[src_id_node] = src_id_node - dst_id_node;
 677     if(!is_deadnode[src_id_node]) {  //if the node isn't dead, copy it
 678       vec3_t x;
 679       m_Grid->GetPoints()->GetPoint(src_id_node, x.data());
 680       dst->GetPoints()->SetPoint(dst_id_node, x.data());
 681       copyNodeData(m_Grid, src_id_node, dst, dst_id_node);
 682       dst_id_node++;
 683     } else {
 684       if(DebugLevel > 0) {
 685         cout << "dead node encountered: src_id_node=" << src_id_node << " dst_id_node=" << dst_id_node << endl;
 686       }
 687     }
 688   }
 689
 690   //Copy undead cells
 691
 692   // Fill is_deadcell
 693   QVector<bool> is_deadcell(m_Grid->GetNumberOfCells(), false);
 694   foreach(vtkIdType id_cell, all_deadcells) {
 695     if( m_Grid->GetCellType(id_cell) == VTK_WEDGE ) EG_BUG;
 696     is_deadcell[id_cell] = true;
 697   }
 698
 699   // Fill is_mutatedcell
 700   QVector<bool> is_mutatedcell(m_Grid->GetNumberOfCells(), false);
 701   foreach(vtkIdType id_cell, all_mutatedcells) {
 702     if( m_Grid->GetCellType(id_cell) == VTK_WEDGE ) EG_BUG;
 703     is_mutatedcell[id_cell] = true;
 704   }
 705
 706   for(vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) { //loop through src cells
 707     //    if( m_Grid->GetCellType(id_cell) == VTK_WEDGE ) continue;
 708     //    if(isVolume(id_cell, m_Grid)) continue;
 709
 710     if(!is_deadcell[id_cell]) {  //if the cell isn't dead
 711       vtkIdType src_num_pts, *src_pts;
 712       m_Grid->GetCellPoints(id_cell, src_num_pts, src_pts);
 713       vtkIdType type_cell = m_Grid->GetCellType(id_cell);
 714
 715       vtkIdType dst_num_pts = src_num_pts;
 716       QVector<vtkIdType> dst_pts(dst_num_pts);
 717
 718       if(is_mutatedcell[id_cell]) {  //mutated cell
 719         if(dst_num_pts != 3) {
 720           // Not fully supported yet
 721           qWarning() << "all_mutatedcells=" << all_mutatedcells;
 722           qWarning() << "A non-triangle cell was mutated!";
 723           EG_BUG;
 724         }
 725         int num_deadnode = 0;
 726         for(int i = 0; i < src_num_pts; i++) {
 727           int DeadIndex = glob2dead[src_pts[i]];
 728           if(DeadIndex != -1) { // It is a dead node.
 729             dst_pts[i] = snappoint_vector[DeadIndex] - OffSet[snappoint_vector[DeadIndex]]; // dead node
 730             num_deadnode++;
 731           } else {
 732             dst_pts[i] = src_pts[i] - OffSet[src_pts[i]]; // not a dead node
 733           }
 734         }
 735         if(num_deadnode != 1) {
 736           qWarning() << "FATAL ERROR: Mutated cell has more than one dead node!";
 737           qWarning() << "num_deadnode=" << num_deadnode;
 738           qWarning() << "type_cell=" << type_cell << " VTK_TRIANGLE=" << VTK_TRIANGLE << " VTK_QUAD=" << VTK_QUAD;
 739           EG_BUG;
 740         }
 741       } else { //normal cell
 742
 743         if(DebugLevel > 10) {
 744           cout << "processing normal cell " << id_cell << endl;
 745         }
 746
 747         if(isVolume(id_cell, m_Grid)) {
 748           int num_deadnode = 0;
 749           for(int i = 0; i < src_num_pts; i++) {
 750             int DeadIndex = glob2dead[src_pts[i]];
 751             if(DeadIndex != -1) { // It is a dead node.
 752               dst_pts[i] = snappoint_vector[DeadIndex] - OffSet[snappoint_vector[DeadIndex]]; // dead node
 753               num_deadnode++;
 754             } else {
 755               dst_pts[i] = src_pts[i] - OffSet[src_pts[i]]; // not a dead node
 756             }
 757           }
 758           if(num_deadnode > 1) {
 759             qWarning() << "FATAL ERROR: Mutated cell has more than one dead node!";
 760             qWarning() << "num_deadnode=" << num_deadnode;
 761             qWarning() << "type_cell=" << type_cell << " VTK_TRIANGLE=" << VTK_TRIANGLE << " VTK_QUAD=" << VTK_QUAD;
 762             for(int k = 0; k < src_num_pts; k++) {
 763               int DeadIndex = glob2dead[src_pts[k]];
 764               qWarning()<<"k="<<k<<" DeadIndex="<<"glob2dead["<<src_pts[k]<<"]="<<DeadIndex;
 765             }
 766             EG_BUG;
 767           }
 768         }
 769         else {
 770           for(int j = 0; j < src_num_pts; j++) {
 771             if(is_deadnode[src_pts[j]]) {
 772               qWarning() << "FATAL ERROR: Normal cell contains a dead node!";
 773               qWarning() << "is_deadnode=" << is_deadnode;
 774               qWarning() << "src_pts[" << j << "]=" << src_pts[j];
 775               qWarning() << "type_cell=" << type_cell << " VTK_TRIANGLE=" << VTK_TRIANGLE << " VTK_QUAD=" << VTK_QUAD;
 776               saveGrid(m_Grid, "crash");
 777               EG_BUG;
 778             }
 779             dst_pts[j] = src_pts[j] - OffSet[src_pts[j]];
 780           }
 781         }
 782
 783       } // end of normal cell processing
 784
 785       // copy the cell
 786       //\todo adapt type_cell in the case of mutilated cells!
 787       vtkIdType id_new_cell = dst->InsertNextCell(type_cell, dst_num_pts, dst_pts.data());
 788       copyCellData(m_Grid, id_cell, dst, id_new_cell);
 789     } //end of undead cell processing
 790   } //end of loop through src cells
 791
 792   // update m_Grid
 793   makeCopy(dst, m_Grid);
 794
 795   int final_num_points = m_Grid->GetNumberOfPoints();
 796   if(initial_num_points - final_num_points != deadnode_vector.size()) {
 797     EG_BUG;
 798   }
 799
 800   return(true);
 801 }
 802 //End of DeleteSetOfPoints