src/libengrid/removepoints.cpp

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