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