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