src/operation.cpp

   1 //
   2 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   3 // +                                                                      +
   4 // + This file is part of enGrid.                                         +
   5 // +                                                                      +
   6 // + Copyright 2008,2009 Oliver Gloth                                     +
   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 "operation.h"
  24 #include "guimainwindow.h"
  25 #include "vtkTriangleFilter.h"
  26 #include "vtkInformation.h"
  27 #include "vtkInformationVector.h"
  28 #include "vtkObjectFactory.h"
  29 #include "vtkPointData.h"
  30 #include "vtkPolyData.h"
  31 #include "vtkPolygon.h"
  32 #include "vtkStreamingDemandDrivenPipeline.h"
  33 #include "vtkCellArray.h"
  34 #include "vtkCellData.h"
  35 #include "vtkCellLocator.h"
  36 #include "vtkFloatArray.h"
  37 #include "vtkMath.h"
  38 #include <vtkCharArray.h>
  39 #include "egvtkobject.h"
  40
  41 #include "geometrytools.h"
  42 using namespace GeometryTools;
  43
  44 #include <QApplication>
  45
  46 QSet<Operation*> Operation::garbage_operations;
  47
  48 void Operation::collectGarbage()
  49 {
  50   QSet<Operation*> delete_operations;
  51   foreach (Operation *op, garbage_operations) {
  52     if (!op->getThread().isRunning()) {
  53       delete_operations.insert(op);
  54       cout << "deleting Operation " << op << endl;
  55       delete op;
  56     };
  57   };
  58   foreach (Operation *op, delete_operations) {
  59     garbage_operations.remove(op);
  60   };
  61 };
  62
  63 Operation::Operation()
  64 {
  65   grid = NULL;
  66   gui = false;
  67   m_quicksave = false;
  68   m_resetoperationcounter = false;
  69   err = NULL;
  70   autoset = true;
  71
  72   //default values for determining node types and for smoothing operations
  73   Convergence=0;
  74   NumberOfIterations=20;
  75   RelaxationFactor=0.01;
  76   FeatureEdgeSmoothing=1;//0 by default in VTK, but we need 1 to avoid the "potatoe effect" ^^
  77   FeatureAngle=45;
  78   EdgeAngle=15;
  79   BoundarySmoothing=1;
  80   GenerateErrorScalars=0;
  81   GenerateErrorVectors=0;
  82 };
  83
  84 Operation::~Operation()
  85 {
  86   if (err) {
  87     err->display();
  88     delete err;
  89   };
  90 };
  91
  92 void Operation::del()
  93 {
  94   garbage_operations.insert(this);
  95 };
  96
  97 void OperationThread::run()
  98 {
  99   try {
 100     GuiMainWindow::lock();
 101     GuiMainWindow::pointer()->setBusy();
 102     op->operate();
 103   } catch (Error err) {
 104     op->err = new Error();
 105     *(op->err) = err;
 106   };
 107   GuiMainWindow::unlock();
 108   GuiMainWindow::pointer()->setIdle();
 109 };
 110
 111 void Operation::operator()()
 112 {
 113   if (gui) {
 114     if (GuiMainWindow::tryLock()) {
 115       checkGrid();
 116       thread.setOperation(this);
 117       GuiMainWindow::unlock();
 118       thread.start(QThread::LowPriority);
 119     } else {
 120       QMessageBox::warning(NULL, "not permitted", "Operation is not permitted while background process is running!");
 121     };
 122   } else {
 123     checkGrid();
 124     operate();
 125     if(m_resetoperationcounter) GuiMainWindow::pointer()->ResetOperationCounter();
 126     if(m_quicksave) GuiMainWindow::pointer()->QuickSave();
 127   };
 128 };
 129
 130 void Operation::setAllCells()
 131 {
 132   QVector<vtkIdType> all_cells;
 133   getAllCells(all_cells, grid);
 134   setCells(all_cells);
 135 };
 136
 137 void Operation::setAllVolumeCells()
 138 {
 139   QVector<vtkIdType> cells;
 140   getAllVolumeCells(cells, grid);
 141   setCells(cells);
 142 };
 143
 144 void Operation::setAllSurfaceCells()
 145 {
 146   QVector<vtkIdType> cells;
 147   getAllSurfaceCells(cells, grid);
 148   setCells(cells);
 149 };
 150
 151 void Operation::initMapping()
 152 {
 153   nodes_map.resize(nodes.size());
 154   for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
 155     nodes_map[i_nodes] = nodes[i_nodes];
 156   };
 157   cells_map.resize(cells.size());
 158   for (int i_cells = 0; i_cells < cells.size(); ++i_cells) {
 159     cells_map[i_cells] = cells[i_cells];
 160   };
 161 };
 162
 163 void Operation::checkGrid()
 164 {
 165   if (grid == NULL) {
 166     grid = GuiMainWindow::pointer()->getGrid();
 167   };
 168   if ((cells.size() == 0) && autoset) {
 169     setAllCells();
 170   };
 171 };
 172
 173 void Operation::updateActors()
 174 {
 175   mainWindow()->updateActors();
 176 };
 177
 178 GuiMainWindow* Operation::mainWindow()
 179 {
 180   return GuiMainWindow::pointer();
 181 };
 182
 183 void Operation::populateBoundaryCodes(QListWidget *lw)
 184 {
 185   QSet<int> bcs;
 186   mainWindow()->getAllBoundaryCodes(bcs);
 187   foreach(int bc, bcs) {
 188     QListWidgetItem *lwi = new QListWidgetItem(lw);
 189     lwi->setCheckState(Qt::Unchecked);
 190     QString text = "";
 191     QTextStream ts(&text);
 192     ts << bc;
 193     lwi->setText(text);
 194     lwi->setFlags(Qt::ItemIsUserCheckable | Qt::ItemIsEnabled);
 195   };
 196 };
 197
 198 //TODO: Get the EG_BUG error again and figure out where it came from
 199 stencil_t Operation::getStencil(vtkIdType id_cell1, int j1)
 200 {
 201   stencil_t S;
 202   S.valid = true;
 203   S.id_cell1 = id_cell1;
 204   if (c2c[_cells[id_cell1]][j1] != -1) {
 205     S.id_cell2 = cells[c2c[_cells[id_cell1]][j1]];
 206     if (grid->GetCellType(S.id_cell2) != VTK_TRIANGLE) {
 207       EG_BUG;
 208     };
 209     vtkIdType N1, N2, *pts1, *pts2;
 210     grid->GetCellPoints(S.id_cell1, N1, pts1);
 211     grid->GetCellPoints(S.id_cell2, N2, pts2);
 212     if      (j1 == 0) { S.p[0] = pts1[2]; S.p[1] = pts1[0]; S.p[3] = pts1[1]; }
 213     else if (j1 == 1) { S.p[0] = pts1[0]; S.p[1] = pts1[1]; S.p[3] = pts1[2]; }
 214     else if (j1 == 2) { S.p[0] = pts1[1]; S.p[1] = pts1[2]; S.p[3] = pts1[0]; };
 215     bool p2 = false;
 216     if (c2c[_cells[S.id_cell2]][0] != -1) {
 217       if (cells[c2c[_cells[S.id_cell2]][0]] == S.id_cell1) {
 218         S.p[2] = pts2[2];
 219         p2 = true;
 220       };
 221     };
 222     if (c2c[_cells[S.id_cell2]][1] != -1) {
 223       if (cells[c2c[_cells[S.id_cell2]][1]] == S.id_cell1) {
 224         S.p[2] = pts2[0];
 225         p2 = true;
 226       };
 227     };
 228     if (c2c[_cells[S.id_cell2]][2] != -1) {
 229       if (cells[c2c[_cells[S.id_cell2]][2]] == S.id_cell1) {
 230         S.p[2] = pts2[1];
 231         p2 = true;
 232       };
 233     };
 234     if (!p2) {
 235       cout<<"S.id_cell1="<<S.id_cell1<<endl;
 236       cout<<"S.id_cell2="<<S.id_cell2<<endl;
 237       EG_BUG;
 238     };
 239   } else {
 240     S.valid = false;
 241     S.id_cell2 = -1;
 242     vtkIdType N1, *pts1;
 243     grid->GetCellPoints(S.id_cell1, N1, pts1);
 244     if      (j1 == 0) { S.p[0] = pts1[2]; S.p[1] = pts1[0]; S.p[3] = pts1[1]; }
 245     else if (j1 == 1) { S.p[0] = pts1[0]; S.p[1] = pts1[1]; S.p[3] = pts1[2]; }
 246     else if (j1 == 2) { S.p[0] = pts1[1]; S.p[1] = pts1[2]; S.p[3] = pts1[0]; };
 247   };
 248   return S;
 249 };
 250
 251 ostream& operator<<(ostream &out, stencil_t S)
 252 {
 253   out<<"S.id_cell1="<<S.id_cell1<<" ";
 254   out<<"S.id_cell2="<<S.id_cell2<<" ";
 255   out<<"S.valid="<<S.valid<<" ";
 256   out<<"[";
 257   for(int i=0;i<4;i++){
 258     out<<S.p[i];
 259     if(i!=3) out<<",";
 260   }
 261   out<<"]";
 262   return(out);
 263 }
 264
 265 //////////////////////////////////////////////
 266 double CurrentVertexAvgDist(vtkIdType a_vertex,QVector< QSet< int > > &n2n,vtkUnstructuredGrid *a_grid)
 267 {
 268   double total_dist=0;
 269   double avg_dist=0;
 270   int N=n2n[a_vertex].size();
 271   vec3_t C;
 272   a_grid->GetPoint(a_vertex, C.data());
 273   foreach(int i,n2n[a_vertex])
 274   {
 275     vec3_t M;
 276     a_grid->GetPoint(i, M.data());
 277     total_dist+=(M-C).abs();
 278   }
 279   avg_dist=total_dist/(double)N;
 280   return(avg_dist);
 281 }
 282
 283 double CurrentMeshDensity(vtkIdType a_vertex,QVector< QSet< int > > &n2n,vtkUnstructuredGrid *a_grid)
 284 {
 285   double total_dist=0;
 286   double avg_dist=0;
 287   int N=n2n[a_vertex].size();
 288   vec3_t C;
 289   a_grid->GetPoint(a_vertex, C.data());
 290   foreach(int i,n2n[a_vertex])
 291   {
 292     vec3_t M;
 293     a_grid->GetPoint(i, M.data());
 294     total_dist+=(M-C).abs();
 295   }
 296   avg_dist=total_dist/(double)N;
 297   double avg_density=1./avg_dist;
 298   return(avg_density);
 299 }
 300
 301 double DesiredVertexAvgDist(vtkIdType a_vertex,QVector< QSet< int > > &n2n,vtkUnstructuredGrid *a_grid)
 302 {
 303   double total_dist=0;
 304   double avg_dist=0;
 305   EG_VTKDCN(vtkDoubleArray, node_meshdensity, a_grid, "node_meshdensity");
 306   int N=n2n[a_vertex].size();
 307   foreach(int i,n2n[a_vertex])
 308   {
 309     total_dist+=1./node_meshdensity->GetValue(i);
 310   }
 311   avg_dist=total_dist/(double)N;
 312   return(avg_dist);
 313 }
 314
 315 double DesiredMeshDensity(vtkIdType a_vertex,QVector< QSet< int > > &n2n,vtkUnstructuredGrid *a_grid)
 316 {
 317   double total_density=0;
 318   double avg_density=0;
 319   EG_VTKDCN(vtkDoubleArray, node_meshdensity, a_grid, "node_meshdensity");
 320   int N=n2n[a_vertex].size();
 321   foreach(int i,n2n[a_vertex])
 322   {
 323     total_density+=node_meshdensity->GetValue(i);
 324   }
 325   avg_density=total_density/(double)N;
 326   return(avg_density);
 327 }
 328
 329 ///////////////////////////////////////////
 330 vtkIdType Operation::getClosestNode(vtkIdType a_id_node,vtkUnstructuredGrid* a_grid)
 331 {
 332   vec3_t C;
 333   a_grid->GetPoint(a_id_node,C.data());
 334   vtkIdType id_minlen=-1;
 335   double minlen=-1;
 336   foreach(vtkIdType neighbour,n2n[a_id_node])
 337   {
 338     vec3_t M;
 339     a_grid->GetPoint(neighbour,M.data());
 340     double len=(M-C).abs();
 341     if(minlen<0 or len<minlen)
 342     {
 343       minlen=len;
 344       id_minlen=neighbour;
 345     }
 346   }
 347   return(id_minlen);
 348 }
 349
 350 vtkIdType Operation::getFarthestNode(vtkIdType a_id_node,vtkUnstructuredGrid* a_grid)
 351 {
 352   vec3_t C;
 353   a_grid->GetPoint(a_id_node,C.data());
 354   vtkIdType id_maxlen=-1;
 355   double maxlen=-1;
 356   foreach(vtkIdType neighbour,n2n[a_id_node])
 357   {
 358     vec3_t M;
 359     a_grid->GetPoint(neighbour,M.data());
 360     double len=(M-C).abs();
 361     if(maxlen<0 or len>maxlen)
 362     {
 363       maxlen=len;
 364       id_maxlen=neighbour;
 365     }
 366   }
 367   return(id_maxlen);
 368 }
 369
 370 bool Operation::SwapCells(vtkUnstructuredGrid* a_grid, stencil_t S)
 371 {
 372   bool swap = false;
 373   if (S.valid) {
 374     vec3_t x3[4], x3_0(0,0,0);
 375     vec2_t x[4];
 376     for (int k = 0; k < 4; ++k) {
 377       a_grid->GetPoints()->GetPoint(S.p[k], x3[k].data());
 378       x3_0 += x3[k];
 379     };
 380     vec3_t n1 = triNormal(x3[0], x3[1], x3[3]);
 381     vec3_t n2 = triNormal(x3[1], x3[2], x3[3]);
 382     n1.normalise();
 383     n2.normalise();
 384     if ( (n1*n2) > 0.8) {
 385       vec3_t n = n1 + n2;
 386       n.normalise();
 387       vec3_t ex = orthogonalVector(n);
 388       vec3_t ey = ex.cross(n);
 389       for (int k = 0; k < 4; ++k) {
 390         x[k] = vec2_t(x3[k]*ex, x3[k]*ey);
 391       };
 392       vec2_t r1, r2, r3, u1, u2, u3;
 393       r1 = 0.5*(x[0] + x[1]); u1 = turnLeft(x[1] - x[0]);
 394       r2 = 0.5*(x[1] + x[2]); u2 = turnLeft(x[2] - x[1]);
 395       r3 = 0.5*(x[1] + x[3]); u3 = turnLeft(x[3] - x[1]);
 396       double k, l;
 397       vec2_t xm1, xm2;
 398       bool ok = true;
 399       if (intersection(k, l, r1, u1, r3, u3)) {
 400         xm1 = r1 + k*u1;
 401         if (intersection(k, l, r2, u2, r3, u3)) {
 402           xm2 = r2 + k*u2;
 403         } else {
 404           ok = false;
 405         };
 406       } else {
 407         ok = false;
 408         swap = true;
 409       };
 410       if (ok) {
 411         if ((xm1 - x[2]).abs() < (xm1 - x[0]).abs()) {
 412           swap = true;
 413         };
 414         if ((xm2 - x[0]).abs() < (xm2 - x[2]).abs()) {
 415           swap = true;
 416         };
 417       };
 418     };
 419   };
 420   if (swap) {
 421     vtkIdType new_pts1[3], new_pts2[3];
 422     new_pts1[0] = S.p[1];
 423     new_pts1[1] = S.p[2];
 424     new_pts1[2] = S.p[0];
 425     new_pts2[0] = S.p[2];
 426     new_pts2[1] = S.p[3];
 427     new_pts2[2] = S.p[0];
 428     a_grid->ReplaceCell(S.id_cell1, 3, new_pts1);
 429     a_grid->ReplaceCell(S.id_cell2, 3, new_pts2);
 430   };
 431   return(swap);
 432 }
 433
 434 void Operation::quad2triangle(vtkUnstructuredGrid* src,vtkIdType quadcell)
 435 {
 436   vtkIdType type_cell = grid->GetCellType(quadcell);
 437   cout<<"It's a "<<type_cell<<endl;
 438   if(type_cell==VTK_QUAD)
 439   {
 440     cout_grid(cout,src,true,true,true,true);
 441     EG_VTKSP(vtkUnstructuredGrid, dst);
 442     //src grid info
 443     int N_points=src->GetNumberOfPoints();
 444     int N_cells=src->GetNumberOfCells();
 445     allocateGrid(dst,N_cells+1,N_points);
 446
 447     for (vtkIdType id_node = 0; id_node < src->GetNumberOfPoints(); ++id_node) {
 448       vec3_t x;
 449       src->GetPoints()->GetPoint(id_node, x.data());
 450       dst->GetPoints()->SetPoint(id_node, x.data());
 451       copyNodeData(src, id_node, dst, id_node);
 452     };
 453     for (vtkIdType id_cell = 0; id_cell < src->GetNumberOfCells(); ++id_cell) {
 454       vtkIdType N_pts, *pts;
 455       vtkIdType type_cell = src->GetCellType(id_cell);
 456       src->GetCellPoints(id_cell, N_pts, pts);
 457       vtkIdType id_new_cell;
 458       vtkIdType id_new_cell1;
 459       vtkIdType id_new_cell2;
 460       if(id_cell!=quadcell)
 461       {
 462         id_new_cell = dst->InsertNextCell(type_cell, N_pts, pts);
 463         copyCellData(src, id_cell, dst, id_new_cell);
 464       }
 465       else
 466       {
 467         vtkIdType triangle1[3];
 468         vtkIdType triangle2[3];
 469         triangle1[0]=pts[1];
 470         triangle1[1]=pts[3];
 471         triangle1[2]=pts[0];
 472         triangle2[0]=pts[3];
 473         triangle2[1]=pts[1];
 474         triangle2[2]=pts[2];
 475         id_new_cell1 = dst->InsertNextCell(VTK_TRIANGLE, 3, triangle1);
 476         copyCellData(src, id_cell, dst, id_new_cell1);
 477         id_new_cell2 = dst->InsertNextCell(VTK_TRIANGLE, 3, triangle2);
 478         copyCellData(src, id_cell, dst, id_new_cell2);
 479         stencil_t S;
 480         S.id_cell1=id_new_cell1;
 481         S.id_cell2=id_new_cell2;
 482         S.p[0]=pts[0];
 483         S.p[1]=pts[1];
 484         S.p[2]=pts[2];
 485         S.p[3]=pts[3];
 486         S.valid=true;
 487         SwapCells(dst,S);
 488       }
 489     };
 490     cout_grid(cout,dst,true,true,true,true);
 491     makeCopy(dst, src);
 492   }//end of if quad
 493 }
 494
 495 void Operation::quad2triangle(vtkUnstructuredGrid* src,vtkIdType quadcell,vtkIdType MovingPoint)
 496 {
 497   vtkIdType type_cell = grid->GetCellType(quadcell);
 498   cout<<"It's a "<<type_cell<<endl;
 499   if(type_cell==VTK_QUAD)
 500   {
 501     cout_grid(cout,src,true,true,true,true);
 502     EG_VTKSP(vtkUnstructuredGrid, dst);
 503     //src grid info
 504     int N_points=src->GetNumberOfPoints();
 505     int N_cells=src->GetNumberOfCells();
 506     allocateGrid(dst,N_cells+1,N_points);
 507
 508     for (vtkIdType id_node = 0; id_node < src->GetNumberOfPoints(); ++id_node) {
 509       vec3_t x;
 510       src->GetPoints()->GetPoint(id_node, x.data());
 511       dst->GetPoints()->SetPoint(id_node, x.data());
 512       copyNodeData(src, id_node, dst, id_node);
 513     };
 514     for (vtkIdType id_cell = 0; id_cell < src->GetNumberOfCells(); ++id_cell) {
 515       vtkIdType N_pts, *pts;
 516       src->GetCellPoints(id_cell, N_pts, pts);
 517       vtkIdType type_cell = src->GetCellType(id_cell);
 518       vtkIdType id_new_cell;
 519       vtkIdType id_new_cell1;
 520       vtkIdType id_new_cell2;
 521       if(id_cell!=quadcell)
 522       {
 523         id_new_cell = dst->InsertNextCell(type_cell, N_pts, pts);
 524         copyCellData(src, id_cell, dst, id_new_cell);
 525       }
 526       else
 527       {
 528         vtkIdType triangle1[3];
 529         vtkIdType triangle2[3];
 530         if(MovingPoint==pts[0] || MovingPoint==pts[2])
 531         {
 532           triangle1[0]=pts[1];
 533           triangle1[1]=pts[3];
 534           triangle1[2]=pts[0];
 535           triangle2[0]=pts[3];
 536           triangle2[1]=pts[1];
 537           triangle2[2]=pts[2];
 538         }
 539         else
 540         {
 541           triangle1[0]=pts[2];
 542           triangle1[1]=pts[0];
 543           triangle1[2]=pts[1];
 544           triangle2[0]=pts[0];
 545           triangle2[1]=pts[2];
 546           triangle2[2]=pts[3];
 547         }
 548         id_new_cell1 = dst->InsertNextCell(VTK_TRIANGLE, 3, triangle1);
 549         copyCellData(src, id_cell, dst, id_new_cell1);
 550         id_new_cell2 = dst->InsertNextCell(VTK_TRIANGLE, 3, triangle2);
 551         copyCellData(src, id_cell, dst, id_new_cell2);
 552       }
 553     };
 554     cout_grid(cout,dst,true,true,true,true);
 555     makeCopy(dst, src);
 556   }//end of if quad
 557 }
 558
 559 int Operation::NumberOfCommonPoints(vtkIdType node1, vtkIdType node2, bool& IsTetra)
 560 {
 561 //   QVector< QSet< int > >     n2n
 562   QSet <int> node1_neighbours=n2n[node1];
 563   QSet <int> node2_neighbours=n2n[node2];
 564   QSet <int> intersection=node1_neighbours.intersect(node2_neighbours);
 565   int N=intersection.size();
 566   IsTetra=false;
 567   if(N==2)
 568   {
 569     QSet<int>::const_iterator p1=intersection.begin();
 570     QSet<int>::const_iterator p2=p1+1;
 571     vtkIdType intersection1=_nodes[*p1];
 572     vtkIdType intersection2=_nodes[*p2];
 573     if(n2n[intersection1].contains(intersection2))//if there's an edge between intersection1 and intersection2
 574     {
 575       //check if (node1,intersection1,intersection2) and (node2,intersection1,intersection2) are defined as cells!
 576   //     QVector< QSet< int > >         n2c
 577       QSet< int > S1=n2c[intersection1];
 578       QSet< int > S2=n2c[intersection2];
 579       QSet< int > Si=S1.intersect(S2);
 580       int counter=0;
 581       foreach(vtkIdType C,Si){
 582         vtkIdType N_pts, *pts;
 583         grid->GetCellPoints(C, N_pts, pts);
 584         for(int i=0;i<N_pts;i++)
 585         {
 586           if(pts[i]==node1 || pts[i]==node2) counter++;
 587         }
 588       }
 589       if(counter>=2) IsTetra=true;
 590     }
 591   }
 592   return(N);
 593 }
 594
 595 // vtkIdType Operation::FindSnapPoint(vtkUnstructuredGrid *src, vtkIdType DeadNode)
 596 // {
 597 //   return(0);
 598 // }
 599
 600 bool Operation::DeletePoint(vtkUnstructuredGrid *src, vtkIdType DeadNode)
 601 {
 602   EG_VTKSP(vtkUnstructuredGrid, dst);
 603
 604     //src grid info
 605   int N_points=src->GetNumberOfPoints();
 606   int N_cells=src->GetNumberOfCells();
 607   int N_newpoints=-1;
 608   int N_newcells=0;
 609
 610   if(DeadNode<0 || DeadNode>=N_points)
 611   {
 612     cout<<"Warning: Point out of range: DeadNode="<<DeadNode<<" N_points="<<N_points<<endl;
 613     return(false);
 614   }
 615
 616   QSet <vtkIdType> DeadCells;
 617   QSet <vtkIdType> MutatedCells;
 618   QSet <vtkIdType> MutilatedCells;
 619
 620   vtkIdType SnapPoint=-1;
 621   //Find closest point to DeadNode
 622 //   vtkIdType SnapPoint = getClosestNode(DeadNode,src);//DeadNode moves to SnapPoint
 623
 624   foreach(vtkIdType PSP, n2n[DeadNode])
 625   {
 626     bool IsValidSnapPoint=true;
 627
 628     cout<<"====>PSP="<<PSP<<endl;
 629     bool IsTetra=true;
 630     if(NumberOfCommonPoints(DeadNode,PSP,IsTetra)>2)//common point check
 631     {
 632       cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
 633       IsValidSnapPoint=false;
 634     }
 635     if(IsTetra)//tetra check
 636     {
 637       cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
 638       IsValidSnapPoint=false;
 639     }
 640
 641     //count number of points and cells to remove + analyse cell transformations
 642     N_newpoints=-1;
 643     N_newcells=0;
 644     DeadCells.clear();
 645     MutatedCells.clear();
 646     MutilatedCells.clear();
 647     foreach(vtkIdType C, n2c[DeadNode])//loop through potentially dead cells
 648     {
 649       cout<<"C="<<C<<endl;
 650       //get points around cell
 651       vtkIdType N_pts, *pts;
 652       src->GetCellPoints(C, N_pts, pts);
 653
 654       bool ContainsSnapPoint=false;
 655       bool invincible=false;
 656       for(int i=0;i<N_pts;i++)
 657       {
 658         cout<<"pts["<<i<<"]="<<pts[i]<<" and PSP="<<PSP<<endl;
 659         if(pts[i]==PSP) {ContainsSnapPoint=true;}
 660         if(pts[i]!=DeadNode && pts[i]!=PSP &&  n2c[pts[i]].size()<=1) invincible=true;
 661       }
 662       if(ContainsSnapPoint)
 663       {
 664         if(N_pts==3)//dead cell
 665         {
 666           if(invincible)//Check that empty lines aren't left behind when a cell is killed
 667           {
 668             cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
 669             IsValidSnapPoint=false;
 670           }
 671           else
 672           {
 673             DeadCells.insert(C);
 674             N_newcells-=1;
 675             cout<<"cell "<<C<<" has been pwned!"<<endl;
 676           }
 677         }
 678   /*      else if(N_pts==4)//mutilated cell
 679         {
 680           MutilatedCells.insert(C);
 681           cout<<"cell "<<C<<" has lost a limb!"<<endl;
 682         }*/
 683         else
 684         {
 685           cout<<"RED ALERT: Xenomorph detected!"<<endl;
 686           EG_BUG;
 687         }
 688       }
 689       else
 690       {
 691         vtkIdType src_N_pts, *src_pts;
 692         src->GetCellPoints(C, src_N_pts, src_pts);
 693
 694         if(src_N_pts!=3)
 695         {
 696           cout<<"RED ALERT: Xenomorph detected!"<<endl;
 697           EG_BUG;
 698         }
 699
 700         vtkIdType OldTriangle[3];
 701         vtkIdType NewTriangle[3];
 702
 703         for(int i=0;i<src_N_pts;i++)
 704         {
 705           OldTriangle[i]=src_pts[i];
 706           NewTriangle[i]=( (src_pts[i]==DeadNode) ? PSP : src_pts[i] );
 707         }
 708         vec3_t Old_N= triNormal(src, OldTriangle[0], OldTriangle[1], OldTriangle[2]);
 709         vec3_t New_N= triNormal(src, NewTriangle[0], NewTriangle[1], NewTriangle[2]);
 710         double OldArea=Old_N.abs();
 711         double NewArea=New_N.abs();
 712         double scal=Old_N*New_N;
 713         double cross=(Old_N.cross(New_N)).abs();//double-cross on Nar Shadaa B-)
 714
 715         cout<<"OldArea="<<OldArea<<endl;
 716         cout<<"NewArea="<<NewArea<<endl;
 717         cout<<"scal="<<scal<<endl;
 718         cout<<"cross="<<cross<<endl;
 719
 720         if(Old_N*New_N<Old_N*Old_N*1./100.)//area + inversion check
 721         {
 722           cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
 723           IsValidSnapPoint=false;
 724         }
 725
 726   /*      if(NewArea<OldArea*1./100.)
 727         {
 728           cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
 729           IsValidSnapPoint=false;
 730         }
 731
 732         if(abs(cross)>10e-4)
 733         {
 734           cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
 735           IsValidSnapPoint=false;
 736         }*/
 737
 738         //mutated cell
 739         MutatedCells.insert(C);
 740         cout<<"cell "<<C<<" has been infected!"<<endl;
 741       }
 742     }
 743
 744     if(N_cells+N_newcells<=0)//survivor check
 745     {
 746       cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
 747       IsValidSnapPoint=false;
 748     }
 749 /*    if(EmptyVolume(DeadNode,PSP))//simplified volume check
 750     {
 751       cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
 752       IsValidSnapPoint=false;
 753     }*/
 754     if(IsValidSnapPoint) {SnapPoint=PSP; break;}
 755   }//end of loop through potential SnapPoints
 756
 757   cout<<"===>SNAPPOINT="<<SnapPoint<<endl;
 758   if(SnapPoint<0) {cout<<"Sorry no possible SnapPoint found."<<endl; return(false);}
 759
 760   //allocate
 761   cout<<"N_points="<<N_points<<endl;
 762   cout<<"N_cells="<<N_cells<<endl;
 763   cout<<"N_newpoints="<<N_newpoints<<endl;
 764   cout<<"N_newcells="<<N_newcells<<endl;
 765   allocateGrid(dst,N_cells+N_newcells,N_points+N_newpoints);
 766
 767   //vector used to redefine the new point IDs
 768   QVector <vtkIdType> OffSet(N_points);
 769
 770   //copy undead points
 771   vtkIdType dst_id_node=0;
 772   for (vtkIdType src_id_node = 0; src_id_node < N_points; src_id_node++) {//loop through src points
 773     if(src_id_node!=DeadNode)//if the node isn't dead, copy it
 774     {
 775       vec3_t x;
 776       src->GetPoints()->GetPoint(src_id_node, x.data());
 777       dst->GetPoints()->SetPoint(dst_id_node, x.data());
 778       copyNodeData(src, src_id_node, dst, dst_id_node);
 779       OffSet[src_id_node]=src_id_node-dst_id_node;
 780       dst_id_node++;
 781     }
 782   };
 783
 784   cout<<"DeadCells="<<DeadCells<<endl;
 785   cout<<"MutatedCells="<<MutatedCells<<endl;
 786   cout<<"MutilatedCells="<<MutilatedCells<<endl;
 787
 788   //Copy undead cells
 789   for (vtkIdType id_cell = 0; id_cell < src->GetNumberOfCells(); ++id_cell) {//loop through src cells
 790     if(!DeadCells.contains(id_cell))//if the cell isn't dead
 791     {
 792       vtkIdType src_N_pts, *src_pts;
 793       vtkIdType dst_N_pts, *dst_pts;
 794       src->GetCellPoints(id_cell, src_N_pts, src_pts);
 795
 796       vtkIdType type_cell = src->GetCellType(id_cell);
 797       cout<<"-->id_cell="<<id_cell<<endl;
 798       for(int i=0;i<src_N_pts;i++) cout<<"src_pts["<<i<<"]="<<src_pts[i]<<endl;
 799 //       src->GetCellPoints(id_cell, dst_N_pts, dst_pts);
 800       dst_N_pts=src_N_pts;
 801       dst_pts=new vtkIdType[dst_N_pts];
 802       if(MutatedCells.contains(id_cell))//mutated cell
 803       {
 804         cout<<"processing mutated cell "<<id_cell<<endl;
 805         for(int i=0;i<src_N_pts;i++)
 806         {
 807           if(src_pts[i]==DeadNode) {
 808             cout<<"SnapPoint="<<SnapPoint<<endl;
 809             cout<<"OffSet[SnapPoint]="<<OffSet[SnapPoint]<<endl;
 810             cout<<"src_pts["<<i<<"]="<<src_pts[i]<<endl;
 811             dst_pts[i]=SnapPoint-OffSet[SnapPoint];
 812           }
 813           else dst_pts[i]=src_pts[i]-OffSet[src_pts[i]];
 814         }
 815         cout<<"--->dst_pts:"<<endl;
 816         for(int i=0;i<dst_N_pts;i++) cout<<"dst_pts["<<i<<"]="<<dst_pts[i]<<endl;
 817
 818       }
 819       else if(MutilatedCells.contains(id_cell))//mutilated cell
 820       {
 821         cout<<"processing mutilated cell "<<id_cell<<endl;
 822
 823         if(type_cell==VTK_QUAD) {
 824           type_cell=VTK_TRIANGLE;
 825           dst_N_pts-=1;
 826         }
 827         else {cout<<"FATAL ERROR: Unknown mutilated cell detected! It is not a quad! Potential xenomorph infestation!"<<endl;EG_BUG;}
 828         //merge points
 829         int j=0;
 830         for(int i=0;i<src_N_pts;i++)
 831         {
 832           if(src_pts[i]==SnapPoint) { dst_pts[j]=SnapPoint-OffSet[SnapPoint];j++; }//SnapPoint
 833           else if(src_pts[i]!=DeadNode) { dst_pts[j]=src_pts[i]-OffSet[src_pts[i]];j++; }//pre-snap/dead + post-snap/dead
 834           //do nothing in case of DeadNode
 835         }
 836       }
 837       else//normal cell
 838       {
 839         cout<<"processing normal cell "<<id_cell<<endl;
 840         for(int i=0;i<src_N_pts;i++)
 841         {
 842           dst_pts[i]=src_pts[i]-OffSet[src_pts[i]];
 843         }
 844       }
 845       //copy the cell
 846       vtkIdType id_new_cell = dst->InsertNextCell(type_cell, dst_N_pts, dst_pts);
 847       copyCellData(src, id_cell, dst, id_new_cell);
 848       cout<<"===Copying cell "<<id_cell<<" to "<<id_new_cell<<"==="<<endl;
 849       cout<<"src_pts:"<<endl;
 850       for(int i=0;i<src_N_pts;i++) cout<<"src_pts["<<i<<"]="<<src_pts[i]<<endl;
 851       cout<<"dst_pts:"<<endl;
 852       for(int i=0;i<dst_N_pts;i++) cout<<"dst_pts["<<i<<"]="<<dst_pts[i]<<endl;
 853       cout<<"OffSet="<<OffSet<<endl;
 854       cout<<"===Copying cell end==="<<endl;
 855       delete dst_pts;
 856     }
 857   };
 858 //   cout_grid(cout,dst,true,true,true,true);
 859   makeCopy(dst, src);
 860   return(true);
 861 }
 862
 863 bool Operation::EmptyVolume(vtkIdType DeadNode, vtkIdType PSP)
 864 {
 865   c2c[DeadNode];
 866   c2c[PSP];
 867   return(true);
 868 }
 869
 870 vec3_t Operation::GetCenter(vtkIdType cellId, double& R)
 871 {
 872   vtkIdType *pts, Npts;
 873   grid->GetCellPoints(cellId, Npts, pts);
 874
 875   vec3_t x(0,0,0);
 876   for (vtkIdType i = 0; i < Npts; ++i) {
 877     vec3_t xp;
 878     grid->GetPoints()->GetPoint(pts[i], xp.data());
 879     x += double(1)/Npts * xp;
 880   };
 881
 882   R = 1e99;
 883   for (vtkIdType i = 0; i < Npts; ++i) {
 884     vec3_t xp;
 885     grid->GetPoints()->GetPoint(pts[i], xp.data());
 886     R = min(R, 0.25*(xp-x).abs());
 887   };
 888
 889   return(x);
 890 }
 891
 892 bool Operation::getNeighbours(vtkIdType Boss, QVector <vtkIdType>& Peons, int BC)
 893 {
 894 //   QVector <vtkIdType> Peons;
 895
 896   QSet <int> S1=n2c[Boss];
 897   cout<<"S1="<<S1<<endl;
 898   foreach(vtkIdType PN,n2n[Boss])
 899   {
 900     cout<<"PN="<<PN<<endl;
 901     QSet <int> S2=n2c[PN];
 902     cout<<"S2="<<S2<<endl;
 903     QSet <int> Si=S2.intersect(S1);
 904     cout<<"PN="<<PN<<" Si="<<Si<<endl;
 905     if(Si.size()<2)//only one common cell
 906     {
 907       Peons.push_back(PN);
 908     }
 909     else
 910     {
 911       QSet <int> bc_set;
 912       foreach(vtkIdType C,Si)
 913       {
 914         EG_VTKDCC(vtkIntArray, cell_code, grid, "cell_code");
 915         int bc=cell_code->GetValue(C);
 916         cout<<"C="<<C<<" bc="<<bc<<endl;
 917         bc_set.insert(bc);
 918       }
 919       if(bc_set.size()>1)//2 different boundary codes
 920       {
 921         Peons.push_back(PN);
 922       }
 923     }
 924   }
 925   if(Peons.size()==2)
 926   {
 927 /*    Peon1=Peons[0];
 928     Peon2=Peons[1];*/
 929     return(true);
 930   }
 931   else
 932   {
 933     int N=n2n[Boss].size();
 934     QVector <vtkIdType> neighbours(N);
 935     qCopy(n2n[Boss].begin(), n2n[Boss].end(), neighbours.begin());
 936
 937     double alphamin_value;
 938     vtkIdType alphamin_i;
 939     vtkIdType alphamin_j;
 940     bool first=true;
 941
 942     for(int i=0;i<N;i++)
 943     {
 944       for(int j=i+1;j<N;j++)
 945       {
 946         double alpha=deviation(grid,neighbours[i],Boss,neighbours[j]);
 947 //         cout<<"alpha("<<neighbours[i]<<","<<Boss<<","<<neighbours[j]<<")="<<alpha<<endl;
 948         if(first) {
 949           alphamin_value=alpha;
 950           alphamin_i=i;
 951           alphamin_j=j;
 952           first=false;
 953         }
 954         else
 955         {
 956           if(alpha<alphamin_value)
 957           {
 958             alphamin_value=alpha;
 959             alphamin_i=i;
 960             alphamin_j=j;
 961           }
 962         }
 963       }
 964     }
 965 //     cout<<"alphamin_value="<<alphamin_value<<endl;
 966
 967     Peons.resize(2);
 968     Peons[0]=neighbours[alphamin_i];
 969     Peons[1]=neighbours[alphamin_j];
 970     return(true);
 971 /*    cout<<"FATAL ERROR: number of neighbours != 2"<<endl;
 972     EG_BUG;*/
 973   }
 974   return(false);
 975 }
 976
 977 int Operation::UpdateMeshDensity()
 978 {
 979   if(DebugLevel>0) cout<<"===UpdateMeshDensity START==="<<endl;
 980
 981   getAllSurfaceCells(cells,grid);
 982   EG_VTKDCC(vtkIntArray, cell_code, grid, "cell_code");
 983   EG_VTKDCN(vtkDoubleArray, node_meshdensity, grid, "node_meshdensity");
 984   getNodesFromCells(cells, nodes, grid);
 985   setGrid(grid);
 986   setCells(cells);
 987
 988   if(DebugLevel>5) cout<<"cells.size()="<<cells.size()<<endl;
 989
 990   EG_VTKDCN(vtkDoubleArray, node_meshdensity_current, grid, "node_meshdensity_current");
 991   foreach(vtkIdType node,nodes)
 992   {
 993     double L=CurrentVertexAvgDist(node,n2n,grid);
 994     double D=1./L;
 995     node_meshdensity_current->SetValue(node, D);
 996   }
 997   if(DebugLevel>0) cout<<"===UpdateMeshDensity END==="<<endl;
 998   return(0);
 999 }
1000
1001 // Special structure for marking vertices
1002 typedef struct _vtkMeshVertex
1003 {
1004   char      type;
1005   vtkIdList *edges; // connected edges (list of connected point ids)
1006 } vtkMeshVertex, *vtkMeshVertexPtr;
1007
1008 int Operation::UpdateNodeType_all()
1009 {
1010   cout<<"===UpdateNodeType_all START==="<<endl;
1011   if(DebugLevel>47) cout<<"this->FeatureAngle="<<this->FeatureAngle<<endl;
1012   if(DebugLevel>47) cout<<"this->EdgeAngle="<<this->EdgeAngle<<endl;
1013
1014   getAllSurfaceCells(cells,grid);
1015   if(DebugLevel>5) cout<<"cells.size()="<<cells.size()<<endl;
1016
1017   EG_VTKSP(vtkPolyData, pdata);
1018   //   addToPolyData(m_SelectedCells, pdata, grid);
1019   addToPolyData(cells, pdata, grid);
1020
1021   vtkPolyData* input=pdata;
1022
1023   vtkPolyData *source = 0;
1024
1025   vtkIdType numPts, numCells, i, numPolys;
1026   int j, k;
1027   vtkIdType npts = 0;
1028   vtkIdType *pts = 0;
1029   vtkIdType p1, p2;
1030   double x[3], y[3], deltaX[3], xNew[3], conv, maxDist, dist, factor;
1031   double x1[3], x2[3], x3[3], l1[3], l2[3];
1032   double CosFeatureAngle; //Cosine of angle between adjacent polys
1033   double CosEdgeAngle; // Cosine of angle between adjacent edges
1034   double closestPt[3], dist2, *w = NULL;
1035   int iterationNumber, abortExecute;
1036   vtkIdType numSimple=0, numBEdges=0, numFixed=0, numFEdges=0;
1037   vtkPolyData *inMesh, *Mesh;
1038   vtkPoints *inPts;
1039   vtkCellArray *inVerts, *inLines, *inPolys;
1040   vtkPoints *newPts;
1041   vtkMeshVertexPtr Verts;
1042   vtkCellLocator *cellLocator=NULL;
1043
1044     // Check input
1045     //
1046   numPts=input->GetNumberOfPoints();
1047   numCells=input->GetNumberOfCells();
1048   if (numPts < 1 || numCells < 1)
1049   {
1050     cout<<"No data to smooth!"<<endl;
1051     return 1;
1052   }
1053
1054   CosFeatureAngle =
1055     cos((double) vtkMath::DegreesToRadians() * this->FeatureAngle);
1056   CosEdgeAngle = cos((double) vtkMath::DegreesToRadians() * this->EdgeAngle);
1057
1058   if(DebugLevel>5) {
1059     cout<<"Smoothing " << numPts << " vertices, " << numCells
1060       << " cells with:\n"
1061       << "\tConvergence= " << this->Convergence << "\n"
1062       << "\tIterations= " << this->NumberOfIterations << "\n"
1063       << "\tRelaxation Factor= " << this->RelaxationFactor << "\n"
1064       << "\tEdge Angle= " << this->EdgeAngle << "\n"
1065       << "\tBoundary Smoothing " << (this->BoundarySmoothing ? "On\n" : "Off\n")
1066       << "\tFeature Edge Smoothing " << (this->FeatureEdgeSmoothing ? "On\n" : "Off\n")
1067       << "\tError Scalars " << (this->GenerateErrorScalars ? "On\n" : "Off\n")
1068       << "\tError Vectors " << (this->GenerateErrorVectors ? "On\n" : "Off\n")<<endl;
1069   }
1070     // Peform topological analysis. What we're gonna do is build a connectivity
1071     // array of connected vertices. The outcome will be one of three
1072     // classifications for a vertex: VTK_SIMPLE_VERTEX, VTK_FIXED_VERTEX. or
1073     // VTK_EDGE_VERTEX. Simple vertices are smoothed using all connected
1074     // vertices. FIXED vertices are never smoothed. Edge vertices are smoothed
1075     // using a subset of the attached vertices.
1076     //
1077   if(DebugLevel>5) cout<<"===>Analyze topology==="<<endl;
1078   if(DebugLevel>5) cout<<"Analyzing topology..."<<endl;
1079   if(DebugLevel>5) cout<<"0:numPts="<<numPts<<endl;
1080   Verts = new vtkMeshVertex[numPts];
1081   for (i=0; i<numPts; i++)
1082   {
1083     if(DebugLevel>5) cout<<"0:VTK_SIMPLE_VERTEX"<<endl;
1084     Verts[i].type = VTK_SIMPLE_VERTEX; //can smooth
1085     Verts[i].edges = NULL;
1086   }
1087
1088   inPts = input->GetPoints();
1089   conv = this->Convergence * input->GetLength();
1090
1091   if(DebugLevel>5) cout<<"==polygons and triangle strips=="<<endl;
1092     // now polygons and triangle strips-------------------------------
1093   inPolys=input->GetPolys();
1094   numPolys = inPolys->GetNumberOfCells();
1095
1096   if(DebugLevel>5) cout<<"numPolys="<<numPolys<<endl;
1097
1098   if ( numPolys > 0 )
1099   { //build cell structure
1100     vtkCellArray *polys;
1101     vtkIdType cellId;
1102     int numNei, nei, edge;
1103     vtkIdType numNeiPts;
1104     vtkIdType *neiPts;
1105     double normal[3], neiNormal[3];
1106     vtkIdList *neighbors;
1107
1108     neighbors = vtkIdList::New();
1109     neighbors->Allocate(VTK_CELL_SIZE);
1110
1111     inMesh = vtkPolyData::New();
1112     inMesh->SetPoints(inPts);
1113     inMesh->SetPolys(inPolys);
1114     Mesh = inMesh;
1115
1116     Mesh->BuildLinks(); //to do neighborhood searching
1117     polys = Mesh->GetPolys();
1118
1119     for (cellId=0, polys->InitTraversal(); polys->GetNextCell(npts,pts);
1120          cellId++)
1121     {
1122       if(DebugLevel>5) cout<<"->cellId="<<cellId<<endl;
1123       for (i=0; i < npts; i++)
1124       {
1125         if(DebugLevel>5) cout<<"-->i="<<i<<endl;
1126         p1 = pts[i];
1127         p2 = pts[(i+1)%npts];
1128
1129         if ( Verts[p1].edges == NULL )
1130         {
1131           Verts[p1].edges = vtkIdList::New();
1132           Verts[p1].edges->Allocate(16,6);
1133         }
1134         if ( Verts[p2].edges == NULL )
1135         {
1136           Verts[p2].edges = vtkIdList::New();
1137           Verts[p2].edges->Allocate(16,6);
1138         }
1139
1140         Mesh->GetCellEdgeNeighbors(cellId,p1,p2,neighbors);
1141         numNei = neighbors->GetNumberOfIds();
1142         if(DebugLevel>5) cout<<"-->numNei="<<numNei<<endl;
1143
1144         edge = VTK_SIMPLE_VERTEX;
1145         if ( numNei == 0 )
1146         {
1147           edge = VTK_BOUNDARY_EDGE_VERTEX;
1148         }
1149
1150         else if ( numNei >= 2 )
1151         {
1152             // check to make sure that this edge hasn't been marked already
1153           for (j=0; j < numNei; j++)
1154           {
1155             if ( neighbors->GetId(j) < cellId )
1156             {
1157               break;
1158             }
1159           }
1160           if ( j >= numNei )
1161           {
1162             edge = VTK_FEATURE_EDGE_VERTEX;
1163           }
1164         }
1165
1166         else if ( numNei == 1 && (nei=neighbors->GetId(0)) > cellId )
1167         {
1168           vtkPolygon::ComputeNormal(inPts,npts,pts,normal);
1169           Mesh->GetCellPoints(nei,numNeiPts,neiPts);
1170           vtkPolygon::ComputeNormal(inPts,numNeiPts,neiPts,neiNormal);
1171
1172           if ( this->FeatureEdgeSmoothing &&
1173                vtkMath::Dot(normal,neiNormal) <= CosFeatureAngle )
1174           {
1175             edge = VTK_FEATURE_EDGE_VERTEX;
1176           }
1177         }
1178         else // a visited edge; skip rest of analysis
1179         {
1180           continue;
1181         }
1182
1183         if ( edge && Verts[p1].type == VTK_SIMPLE_VERTEX )
1184         {
1185           Verts[p1].edges->Reset();
1186           Verts[p1].edges->InsertNextId(p2);
1187           Verts[p1].type = edge;
1188         }
1189         else if ( (edge && Verts[p1].type == VTK_BOUNDARY_EDGE_VERTEX) ||
1190                   (edge && Verts[p1].type == VTK_FEATURE_EDGE_VERTEX) ||
1191                   (!edge && Verts[p1].type == VTK_SIMPLE_VERTEX ) )
1192         {
1193           Verts[p1].edges->InsertNextId(p2);
1194           if ( Verts[p1].type && edge == VTK_BOUNDARY_EDGE_VERTEX )
1195           {
1196             Verts[p1].type = VTK_BOUNDARY_EDGE_VERTEX;
1197           }
1198         }
1199
1200         if ( edge && Verts[p2].type == VTK_SIMPLE_VERTEX )
1201         {
1202           Verts[p2].edges->Reset();
1203           Verts[p2].edges->InsertNextId(p1);
1204           Verts[p2].type = edge;
1205         }
1206         else if ( (edge && Verts[p2].type == VTK_BOUNDARY_EDGE_VERTEX ) ||
1207                   (edge && Verts[p2].type == VTK_FEATURE_EDGE_VERTEX) ||
1208                   (!edge && Verts[p2].type == VTK_SIMPLE_VERTEX ) )
1209         {
1210           Verts[p2].edges->InsertNextId(p1);
1211           if ( Verts[p2].type && edge == VTK_BOUNDARY_EDGE_VERTEX )
1212           {
1213             Verts[p2].type = VTK_BOUNDARY_EDGE_VERTEX;
1214           }
1215         }
1216       }
1217     }
1218
1219     inMesh->Delete();
1220
1221     neighbors->Delete();
1222   }//if strips or polys
1223
1224     //post-process edge vertices to make sure we can smooth them
1225   for (i=0; i<numPts; i++)
1226   {
1227     if ( Verts[i].type == VTK_SIMPLE_VERTEX )
1228     {
1229       numSimple++;
1230     }
1231
1232     else if ( Verts[i].type == VTK_FIXED_VERTEX )
1233     {
1234       numFixed++;
1235     }
1236
1237     else if ( Verts[i].type == VTK_FEATURE_EDGE_VERTEX ||
1238               Verts[i].type == VTK_BOUNDARY_EDGE_VERTEX )
1239     { //see how many edges; if two, what the angle is
1240
1241       if ( !this->BoundarySmoothing &&
1242            Verts[i].type == VTK_BOUNDARY_EDGE_VERTEX )
1243       {
1244         if(DebugLevel>5) cout<<"Verts[i].type = VTK_FIXED_VERTEX; 4"<<endl;
1245         Verts[i].type = VTK_FIXED_VERTEX;
1246         numBEdges++;
1247       }
1248
1249       else if ( (npts = Verts[i].edges->GetNumberOfIds()) != 2 )
1250       {
1251         if(DebugLevel>5) cout<<"Verts["<<i<<"].type = VTK_FIXED_VERTEX; 5"<<endl;
1252         Verts[i].type = VTK_FIXED_VERTEX;
1253         numFixed++;
1254       }
1255
1256       else //check angle between edges
1257       {
1258         inPts->GetPoint(Verts[i].edges->GetId(0),x1);
1259         inPts->GetPoint(i,x2);
1260         inPts->GetPoint(Verts[i].edges->GetId(1),x3);
1261
1262         for (k=0; k<3; k++)
1263         {
1264           l1[k] = x2[k] - x1[k];
1265           l2[k] = x3[k] - x2[k];
1266         }
1267         if ( vtkMath::Normalize(l1) >= 0.0 &&
1268              vtkMath::Normalize(l2) >= 0.0 &&
1269              vtkMath::Dot(l1,l2) < CosEdgeAngle)
1270         {
1271           if(DebugLevel>5) cout<<"Verts["<<i<<"].type = VTK_FIXED_VERTEX; 6"<<endl;
1272           Verts[i].type = VTK_FIXED_VERTEX;
1273           numFixed++;
1274         }
1275         else
1276         {
1277           if ( Verts[i].type == VTK_FEATURE_EDGE_VERTEX )
1278           {
1279             numFEdges++;
1280           }
1281           else
1282           {
1283             numBEdges++;
1284           }
1285         }
1286       }//if along edge
1287     }//if edge vertex
1288   }//for all points
1289
1290   if(DebugLevel>5) {
1291     cout<<"Found\n\t" << numSimple << " simple vertices\n\t"
1292       << numFEdges << " feature edge vertices\n\t"
1293       << numBEdges << " boundary edge vertices\n\t"
1294       << numFixed << " fixed vertices\n\t"<<endl;
1295     cout<<"1:numPts="<<numPts<<endl;
1296   }
1297
1298   for (i=0; i<numPts; i++)
1299   {
1300     if(DebugLevel>5) cout<<"Verts["<<i<<"].type="<<VertexType2Str(Verts[i].type)<<endl;
1301     if(Verts[i].edges != NULL && (npts = Verts[i].edges->GetNumberOfIds()) > 0)
1302     {
1303       for (j=0; j<npts; j++)
1304       {
1305         if(DebugLevel>5) cout<<"Verts["<<i<<"].edges->GetId("<<j<<")="<<Verts[i].edges->GetId(j)<<endl;
1306       }//for all connected points
1307     }
1308   }
1309
1310   //Copy node type info from Verts
1311   EG_VTKDCN(vtkCharArray, node_type, grid, "node_type");
1312   if(DebugLevel>5) cout<<"nodes.size()="<<nodes.size()<<endl;
1313   foreach(vtkIdType node,nodes)
1314   {
1315     if(DebugLevel>5) cout<<"Verts["<<node<<"].type="<<VertexType2Str(Verts[node].type)<<endl;
1316     node_type->SetValue(node,Verts[node].type);
1317   }
1318
1319   //free up connectivity storage
1320   for (int i=0; i<numPts; i++)
1321   {
1322     if ( Verts[i].edges != NULL )
1323     {
1324       Verts[i].edges->Delete();
1325       Verts[i].edges = NULL;
1326     }
1327   }
1328   delete [] Verts;
1329
1330   cout<<"===UpdateNodeType_all END==="<<endl;
1331   return(0);
1332 }
1333 //End of UpdateNodeType_all
1334
1335 int Operation::UpdateNodeType()
1336 {
1337   if(DebugLevel>47) cout<<"this->FeatureAngle="<<this->FeatureAngle<<endl;
1338   if(DebugLevel>47) cout<<"this->EdgeAngle="<<this->EdgeAngle<<endl;
1339   cout<<"===UpdateNodeType START==="<<endl;
1340
1341   getAllSurfaceCells(cells,grid);
1342   if(DebugLevel>5) cout<<"cells.size()="<<cells.size()<<endl;
1343
1344   EG_VTKSP(vtkPolyData, pdata);
1345   //   addToPolyData(m_SelectedCells, pdata, grid);
1346   addToPolyData(cells, pdata, grid);
1347
1348   vtkPolyData* input=pdata;
1349
1350   vtkPolyData *source = 0;
1351
1352   vtkIdType numPts, numCells, i, numPolys;
1353   int j, k;
1354   vtkIdType npts = 0;
1355   vtkIdType *pts = 0;
1356   vtkIdType p1, p2;
1357   double x[3], y[3], deltaX[3], xNew[3], conv, maxDist, dist, factor;
1358   double x1[3], x2[3], x3[3], l1[3], l2[3];
1359   double CosFeatureAngle; //Cosine of angle between adjacent polys
1360   double CosEdgeAngle; // Cosine of angle between adjacent edges
1361   double closestPt[3], dist2, *w = NULL;
1362   int iterationNumber, abortExecute;
1363   vtkIdType numSimple=0, numBEdges=0, numFixed=0, numFEdges=0;
1364   vtkPolyData *inMesh, *Mesh;
1365   vtkPoints *inPts;
1366   vtkCellArray *inVerts, *inLines, *inPolys;
1367   vtkPoints *newPts;
1368   vtkMeshVertexPtr Verts;
1369   vtkCellLocator *cellLocator=NULL;
1370
1371     // Check input
1372     //
1373   numPts=input->GetNumberOfPoints();
1374   numCells=input->GetNumberOfCells();
1375   if (numPts < 1 || numCells < 1)
1376   {
1377     cout<<"No data to smooth!"<<endl;
1378     return 1;
1379   }
1380
1381   CosFeatureAngle =
1382     cos((double) vtkMath::DegreesToRadians() * this->FeatureAngle);
1383   CosEdgeAngle = cos((double) vtkMath::DegreesToRadians() * this->EdgeAngle);
1384
1385   if(DebugLevel>5) {
1386     cout<<"Smoothing " << numPts << " vertices, " << numCells
1387       << " cells with:\n"
1388       << "\tConvergence= " << this->Convergence << "\n"
1389       << "\tIterations= " << this->NumberOfIterations << "\n"
1390       << "\tRelaxation Factor= " << this->RelaxationFactor << "\n"
1391       << "\tEdge Angle= " << this->EdgeAngle << "\n"
1392       << "\tBoundary Smoothing " << (this->BoundarySmoothing ? "On\n" : "Off\n")
1393       << "\tFeature Edge Smoothing " << (this->FeatureEdgeSmoothing ? "On\n" : "Off\n")
1394       << "\tError Scalars " << (this->GenerateErrorScalars ? "On\n" : "Off\n")
1395       << "\tError Vectors " << (this->GenerateErrorVectors ? "On\n" : "Off\n")<<endl;
1396   }
1397     // Peform topological analysis. What we're gonna do is build a connectivity
1398     // array of connected vertices. The outcome will be one of three
1399     // classifications for a vertex: VTK_SIMPLE_VERTEX, VTK_FIXED_VERTEX. or
1400     // VTK_EDGE_VERTEX. Simple vertices are smoothed using all connected
1401     // vertices. FIXED vertices are never smoothed. Edge vertices are smoothed
1402     // using a subset of the attached vertices.
1403     //
1404   if(DebugLevel>5) cout<<"===>Analyze topology==="<<endl;
1405   if(DebugLevel>5) cout<<"Analyzing topology..."<<endl;
1406   if(DebugLevel>5) cout<<"0:numPts="<<numPts<<endl;
1407   Verts = new vtkMeshVertex[numPts];
1408   for (i=0; i<numPts; i++)
1409   {
1410     if(DebugLevel>5) cout<<"0:VTK_SIMPLE_VERTEX"<<endl;
1411     Verts[i].type = VTK_SIMPLE_VERTEX; //can smooth
1412     Verts[i].edges = NULL;
1413   }
1414
1415   inPts = input->GetPoints();
1416   conv = this->Convergence * input->GetLength();
1417
1418   if(DebugLevel>5) cout<<"==polygons and triangle strips=="<<endl;
1419     // now polygons and triangle strips-------------------------------
1420   inPolys=input->GetPolys();
1421   numPolys = inPolys->GetNumberOfCells();
1422
1423   if(DebugLevel>5) cout<<"numPolys="<<numPolys<<endl;
1424
1425   if ( numPolys > 0 )
1426   { //build cell structure
1427     vtkCellArray *polys;
1428     vtkIdType cellId;
1429     int numNei, nei, edge;
1430     vtkIdType numNeiPts;
1431     vtkIdType *neiPts;
1432     double normal[3], neiNormal[3];
1433     vtkIdList *neighbors;
1434
1435     neighbors = vtkIdList::New();
1436     neighbors->Allocate(VTK_CELL_SIZE);
1437
1438     inMesh = vtkPolyData::New();
1439     inMesh->SetPoints(inPts);
1440     inMesh->SetPolys(inPolys);
1441     Mesh = inMesh;
1442
1443     Mesh->BuildLinks(); //to do neighborhood searching
1444     polys = Mesh->GetPolys();
1445
1446     for (cellId=0, polys->InitTraversal(); polys->GetNextCell(npts,pts);
1447          cellId++)
1448     {
1449       if(DebugLevel>5) cout<<"->cellId="<<cellId<<endl;
1450       for (i=0; i < npts; i++)
1451       {
1452         if(DebugLevel>5) cout<<"-->i="<<i<<endl;
1453         p1 = pts[i];
1454         p2 = pts[(i+1)%npts];
1455
1456         if ( Verts[p1].edges == NULL )
1457         {
1458           Verts[p1].edges = vtkIdList::New();
1459           Verts[p1].edges->Allocate(16,6);
1460         }
1461         if ( Verts[p2].edges == NULL )
1462         {
1463           Verts[p2].edges = vtkIdList::New();
1464           Verts[p2].edges->Allocate(16,6);
1465         }
1466
1467         Mesh->GetCellEdgeNeighbors(cellId,p1,p2,neighbors);
1468         numNei = neighbors->GetNumberOfIds();
1469         if(DebugLevel>5) cout<<"-->numNei="<<numNei<<endl;
1470
1471         edge = VTK_SIMPLE_VERTEX;
1472         if ( numNei == 0 )
1473         {
1474           edge = VTK_BOUNDARY_EDGE_VERTEX;
1475         }
1476
1477         else if ( numNei >= 2 )
1478         {
1479             // check to make sure that this edge hasn't been marked already
1480           for (j=0; j < numNei; j++)
1481           {
1482             if ( neighbors->GetId(j) < cellId )
1483             {
1484               break;
1485             }
1486           }
1487           if ( j >= numNei )
1488           {
1489             edge = VTK_FEATURE_EDGE_VERTEX;
1490           }
1491         }
1492
1493         else if ( numNei == 1 && (nei=neighbors->GetId(0)) > cellId )
1494         {
1495           vtkPolygon::ComputeNormal(inPts,npts,pts,normal);
1496           Mesh->GetCellPoints(nei,numNeiPts,neiPts);
1497           vtkPolygon::ComputeNormal(inPts,numNeiPts,neiPts,neiNormal);
1498
1499           if ( this->FeatureEdgeSmoothing &&
1500                vtkMath::Dot(normal,neiNormal) <= CosFeatureAngle )
1501           {
1502             edge = VTK_FEATURE_EDGE_VERTEX;
1503           }
1504         }
1505         else // a visited edge; skip rest of analysis
1506         {
1507           continue;
1508         }
1509
1510         if ( edge && Verts[p1].type == VTK_SIMPLE_VERTEX )
1511         {
1512           Verts[p1].edges->Reset();
1513           Verts[p1].edges->InsertNextId(p2);
1514           Verts[p1].type = edge;
1515         }
1516         else if ( (edge && Verts[p1].type == VTK_BOUNDARY_EDGE_VERTEX) ||
1517                   (edge && Verts[p1].type == VTK_FEATURE_EDGE_VERTEX) ||
1518                   (!edge && Verts[p1].type == VTK_SIMPLE_VERTEX ) )
1519         {
1520           Verts[p1].edges->InsertNextId(p2);
1521           if ( Verts[p1].type && edge == VTK_BOUNDARY_EDGE_VERTEX )
1522           {
1523             Verts[p1].type = VTK_BOUNDARY_EDGE_VERTEX;
1524           }
1525         }
1526
1527         if ( edge && Verts[p2].type == VTK_SIMPLE_VERTEX )
1528         {
1529           Verts[p2].edges->Reset();
1530           Verts[p2].edges->InsertNextId(p1);
1531           Verts[p2].type = edge;
1532         }
1533         else if ( (edge && Verts[p2].type == VTK_BOUNDARY_EDGE_VERTEX ) ||
1534                   (edge && Verts[p2].type == VTK_FEATURE_EDGE_VERTEX) ||
1535                   (!edge && Verts[p2].type == VTK_SIMPLE_VERTEX ) )
1536         {
1537           Verts[p2].edges->InsertNextId(p1);
1538           if ( Verts[p2].type && edge == VTK_BOUNDARY_EDGE_VERTEX )
1539           {
1540             Verts[p2].type = VTK_BOUNDARY_EDGE_VERTEX;
1541           }
1542         }
1543       }
1544     }
1545
1546     inMesh->Delete();
1547
1548     neighbors->Delete();
1549   }//if strips or polys
1550
1551     //post-process edge vertices to make sure we can smooth them
1552   for (i=0; i<numPts; i++)
1553   {
1554     if ( Verts[i].type == VTK_SIMPLE_VERTEX )
1555     {
1556       numSimple++;
1557     }
1558
1559     else if ( Verts[i].type == VTK_FIXED_VERTEX )
1560     {
1561       numFixed++;
1562     }
1563
1564     else if ( Verts[i].type == VTK_FEATURE_EDGE_VERTEX ||
1565               Verts[i].type == VTK_BOUNDARY_EDGE_VERTEX )
1566     { //see how many edges; if two, what the angle is
1567
1568       if ( !this->BoundarySmoothing &&
1569            Verts[i].type == VTK_BOUNDARY_EDGE_VERTEX )
1570       {
1571         if(DebugLevel>5) cout<<"Verts[i].type = VTK_FIXED_VERTEX; 4"<<endl;
1572         Verts[i].type = VTK_FIXED_VERTEX;
1573         numBEdges++;
1574       }
1575
1576       else if ( (npts = Verts[i].edges->GetNumberOfIds()) != 2 )
1577       {
1578         if(DebugLevel>5) cout<<"Verts["<<i<<"].type = VTK_FIXED_VERTEX; 5"<<endl;
1579         Verts[i].type = VTK_FIXED_VERTEX;
1580         numFixed++;
1581       }
1582
1583       else //check angle between edges
1584       {
1585         inPts->GetPoint(Verts[i].edges->GetId(0),x1);
1586         inPts->GetPoint(i,x2);
1587         inPts->GetPoint(Verts[i].edges->GetId(1),x3);
1588
1589         for (k=0; k<3; k++)
1590         {
1591           l1[k] = x2[k] - x1[k];
1592           l2[k] = x3[k] - x2[k];
1593         }
1594         if ( vtkMath::Normalize(l1) >= 0.0 &&
1595              vtkMath::Normalize(l2) >= 0.0 &&
1596              vtkMath::Dot(l1,l2) < CosEdgeAngle)
1597         {
1598           if(DebugLevel>5) cout<<"Verts["<<i<<"].type = VTK_FIXED_VERTEX; 6"<<endl;
1599           Verts[i].type = VTK_FIXED_VERTEX;
1600           numFixed++;
1601         }
1602         else
1603         {
1604           if ( Verts[i].type == VTK_FEATURE_EDGE_VERTEX )
1605           {
1606             numFEdges++;
1607           }
1608           else
1609           {
1610             numBEdges++;
1611           }
1612         }
1613       }//if along edge
1614     }//if edge vertex
1615   }//for all points
1616
1617   if(DebugLevel>5) {
1618     cout<<"Found\n\t" << numSimple << " simple vertices\n\t"
1619       << numFEdges << " feature edge vertices\n\t"
1620       << numBEdges << " boundary edge vertices\n\t"
1621       << numFixed << " fixed vertices\n\t"<<endl;
1622     cout<<"1:numPts="<<numPts<<endl;
1623   }
1624
1625   for (i=0; i<numPts; i++)
1626   {
1627     if(DebugLevel>5) cout<<"Verts["<<i<<"].type="<<VertexType2Str(Verts[i].type)<<endl;
1628     if(Verts[i].edges != NULL && (npts = Verts[i].edges->GetNumberOfIds()) > 0)
1629     {
1630       for (j=0; j<npts; j++)
1631       {
1632         if(DebugLevel>5) cout<<"Verts["<<i<<"].edges->GetId("<<j<<")="<<Verts[i].edges->GetId(j)<<endl;
1633       }//for all connected points
1634     }
1635   }
1636
1637   //Copy node type info from Verts
1638   EG_VTKDCN(vtkCharArray, node_type, grid, "node_type");
1639   if(DebugLevel>5) cout<<"nodes.size()="<<nodes.size()<<endl;
1640   foreach(vtkIdType node,nodes)
1641   {
1642     if(DebugLevel>5) cout<<"Verts["<<node<<"].type="<<VertexType2Str(Verts[node].type)<<endl;
1643     node_type->SetValue(node,Verts[node].type);
1644   }
1645
1646   //free up connectivity storage
1647   for (int i=0; i<numPts; i++)
1648   {
1649     if ( Verts[i].edges != NULL )
1650     {
1651       Verts[i].edges->Delete();
1652       Verts[i].edges = NULL;
1653     }
1654   }
1655   delete [] Verts;
1656
1657   cout<<"===UpdateNodeType END==="<<endl;
1658   return(0);
1659 }
1660 //End of UpdateNodeType
1661
1662 // DEFINITIONS:
1663 // Normal cell: nothing has changed
1664 // Dead cell: the cell does not exist anymore
1665 // Mutated cell: the cell's form has changed
1666 // Mutilated cell: the cell has less points than before
1667
1668 vtkIdType Operation::FindSnapPoint(vtkUnstructuredGrid *src, vtkIdType DeadNode,QSet <vtkIdType> & DeadCells,QSet <vtkIdType> & MutatedCells,QSet <vtkIdType> & MutilatedCells, int& N_newpoints, int& N_newcells)
1669 {
1670   getAllSurfaceCells(cells,src);
1671   getNodesFromCells(cells, nodes, src);
1672   setGrid(src);
1673   setCells(cells);
1674
1675   UpdateNodeType_all();
1676
1677   EG_VTKDCN(vtkCharArray, node_type, src, "node_type");
1678   if(node_type->GetValue(DeadNode)==VTK_FIXED_VERTEX)
1679   {
1680     cout<<"Sorry, unable to remove fixed vertex."<<endl;
1681     return(-1);
1682   }
1683
1684     //src grid info
1685   int N_points=src->GetNumberOfPoints();
1686   int N_cells=src->GetNumberOfCells();
1687   N_newpoints=-1;
1688   N_newcells=0;
1689
1690   vtkIdType SnapPoint=-1;
1691
1692   foreach(vtkIdType PSP, n2n[DeadNode])
1693   {
1694     bool IsValidSnapPoint=true;
1695
1696     if(DebugLevel>10) cout<<"====>PSP="<<PSP<<endl;
1697     bool IsTetra=true;
1698     if(NumberOfCommonPoints(DeadNode,PSP,IsTetra)>2)//common point check
1699     {
1700       if(DebugLevel>10) cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
1701       IsValidSnapPoint=false;
1702     }
1703     if(IsTetra)//tetra check
1704     {
1705       if(DebugLevel>10) cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
1706       IsValidSnapPoint=false;
1707     }
1708
1709     //count number of points and cells to remove + analyse cell transformations
1710     N_newpoints=-1;
1711     N_newcells=0;
1712     DeadCells.clear();
1713     MutatedCells.clear();
1714     MutilatedCells.clear();
1715     foreach(vtkIdType C, n2c[DeadNode])//loop through potentially dead cells
1716     {
1717       //get points around cell
1718       vtkIdType N_pts, *pts;
1719       src->GetCellPoints(C, N_pts, pts);
1720
1721       bool ContainsSnapPoint=false;
1722       bool invincible=false;
1723       for(int i=0;i<N_pts;i++)
1724       {
1725         if(DebugLevel>10) cout<<"pts["<<i<<"]="<<pts[i]<<" and PSP="<<PSP<<endl;
1726         if(pts[i]==PSP) {ContainsSnapPoint=true;}
1727         if(pts[i]!=DeadNode && pts[i]!=PSP &&  n2c[pts[i]].size()<=1) invincible=true;
1728       }
1729       if(ContainsSnapPoint)
1730       {
1731         if(N_pts==3)//dead cell
1732         {
1733           if(invincible)//Check that empty lines aren't left behind when a cell is killed
1734           {
1735             if(DebugLevel>10) cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
1736             IsValidSnapPoint=false;
1737           }
1738           else
1739           {
1740             DeadCells.insert(C);
1741             N_newcells-=1;
1742             if(DebugLevel>10) cout<<"cell "<<C<<" has been pwned!"<<endl;
1743           }
1744         }
1745         else
1746         {
1747           cout<<"RED ALERT: Xenomorph detected!"<<endl;
1748           EG_BUG;
1749         }
1750       }
1751       else
1752       {
1753         vtkIdType src_N_pts, *src_pts;
1754         src->GetCellPoints(C, src_N_pts, src_pts);
1755
1756         if(src_N_pts!=3)
1757         {
1758           cout<<"RED ALERT: Xenomorph detected!"<<endl;
1759           EG_BUG;
1760         }
1761
1762         vtkIdType OldTriangle[3];
1763         vtkIdType NewTriangle[3];
1764
1765         for(int i=0;i<src_N_pts;i++)
1766         {
1767           OldTriangle[i]=src_pts[i];
1768           NewTriangle[i]=( (src_pts[i]==DeadNode) ? PSP : src_pts[i] );
1769         }
1770         vec3_t Old_N= triNormal(src, OldTriangle[0], OldTriangle[1], OldTriangle[2]);
1771         vec3_t New_N= triNormal(src, NewTriangle[0], NewTriangle[1], NewTriangle[2]);
1772         double OldArea=Old_N.abs();
1773         double NewArea=New_N.abs();
1774         double scal=Old_N*New_N;
1775         double cross=(Old_N.cross(New_N)).abs();//double-cross on Nar Shadaa B-)
1776
1777         if(DebugLevel>10) {
1778           cout<<"OldArea="<<OldArea<<endl;
1779           cout<<"NewArea="<<NewArea<<endl;
1780           cout<<"scal="<<scal<<endl;
1781           cout<<"cross="<<cross<<endl;
1782         }
1783
1784         if(Old_N*New_N<Old_N*Old_N*1./100.)//area + inversion check
1785         {
1786           if(DebugLevel>10) cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
1787           IsValidSnapPoint=false;
1788         }
1789
1790         //mutated cell
1791         MutatedCells.insert(C);
1792         if(DebugLevel>10) cout<<"cell "<<C<<" has been infected!"<<endl;
1793       }
1794     }
1795
1796     if(N_cells+N_newcells<=0)//survivor check
1797     {
1798       if(DebugLevel>10) cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
1799       IsValidSnapPoint=false;
1800     }
1801
1802     if(node_type->GetValue(DeadNode)==VTK_BOUNDARY_EDGE_VERTEX && node_type->GetValue(PSP)==VTK_SIMPLE_VERTEX)
1803     {
1804       if(DebugLevel>10) cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
1805       IsValidSnapPoint=false;
1806     }
1807
1808     if(node_type->GetValue(DeadNode)==VTK_BOUNDARY_EDGE_VERTEX)
1809     {
1810       int BC=0;
1811       QVector <vtkIdType> Peons;
1812       getNeighbours(DeadNode, Peons, BC);
1813       if(!Peons.contains(PSP))
1814       {
1815         if(DebugLevel>0) cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
1816         IsValidSnapPoint=false;
1817       }
1818     }
1819
1820     if(node_type->GetValue(DeadNode)==VTK_FEATURE_EDGE_VERTEX && node_type->GetValue(PSP)==VTK_SIMPLE_VERTEX)
1821     {
1822       if(DebugLevel>10) cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
1823       IsValidSnapPoint=false;
1824     }
1825
1826     if(node_type->GetValue(DeadNode)==VTK_FEATURE_EDGE_VERTEX)
1827     {
1828       int BC=0;
1829       QVector <vtkIdType> Peons;
1830       getNeighbours(DeadNode, Peons, BC);
1831       if(!Peons.contains(PSP))
1832       {
1833         if(DebugLevel>0) cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
1834         IsValidSnapPoint=false;
1835       }
1836     }
1837
1838     if(IsValidSnapPoint) {SnapPoint=PSP; break;}
1839   }//end of loop through potential SnapPoints
1840
1841   if(DebugLevel>10)
1842   {
1843     cout<<"AT FINDSNAPPOINT EXIT"<<endl;
1844     cout<<"N_points="<<N_points<<endl;
1845     cout<<"N_cells="<<N_cells<<endl;
1846     cout<<"N_newpoints="<<N_newpoints<<endl;
1847     cout<<"N_newcells="<<N_newcells<<endl;
1848   }
1849   return(SnapPoint);
1850 }
1851 //End of FindSnapPoint
1852
1853 bool Operation::DeletePoint_2(vtkUnstructuredGrid *src, vtkIdType DeadNode, int& N_newpoints, int& N_newcells)
1854 {
1855   getAllSurfaceCells(cells,src);
1856 //   getNodesFromCells(cells, nodes, src);
1857   setGrid(src);
1858   setCells(cells);
1859   UpdateNodeType_all();
1860
1861   //src grid info
1862   int N_points=src->GetNumberOfPoints();
1863   int N_cells=src->GetNumberOfCells();
1864   N_newpoints=-1;
1865   N_newcells=0;
1866
1867   if(DeadNode<0 || DeadNode>=N_points)
1868   {
1869     cout<<"Warning: Point out of range: DeadNode="<<DeadNode<<" N_points="<<N_points<<endl;
1870     return(false);
1871   }
1872
1873   QSet <vtkIdType> DeadCells;
1874   QSet <vtkIdType> MutatedCells;
1875   QSet <vtkIdType> MutilatedCells;
1876
1877   if(DebugLevel>10) {
1878     cout<<"BEFORE FINDSNAPPOINT"<<endl;
1879     cout<<"N_points="<<N_points<<endl;
1880     cout<<"N_cells="<<N_cells<<endl;
1881     cout<<"N_newpoints="<<N_newpoints<<endl;
1882     cout<<"N_newcells="<<N_newcells<<endl;
1883   }
1884   vtkIdType SnapPoint=FindSnapPoint(src,DeadNode,DeadCells,MutatedCells,MutilatedCells, N_newpoints, N_newcells);
1885
1886   if(DebugLevel>0) cout<<"===>DeadNode="<<DeadNode<<" moving to SNAPPOINT="<<SnapPoint<<" DebugLevel="<<DebugLevel<<endl;
1887   if(SnapPoint<0) {cout<<"Sorry no possible SnapPoint found."<<endl; return(false);}
1888
1889   //allocate
1890   if(DebugLevel>10) {
1891     cout<<"BEFORE ALLOCATION"<<endl;
1892     cout<<"N_points="<<N_points<<endl;
1893     cout<<"N_cells="<<N_cells<<endl;
1894     cout<<"N_newpoints="<<N_newpoints<<endl;
1895     cout<<"N_newcells="<<N_newcells<<endl;
1896   }
1897   N_points=src->GetNumberOfPoints();
1898   N_cells=src->GetNumberOfCells();
1899
1900   if(DebugLevel>47) {
1901     cout<<"N_points="<<N_points<<endl;
1902     cout<<"N_cells="<<N_cells<<endl;
1903     cout<<"N_newpoints="<<N_newpoints<<endl;
1904     cout<<"N_newcells="<<N_newcells<<endl;
1905   }
1906   EG_VTKSP(vtkUnstructuredGrid,dst);
1907   allocateGrid(dst,N_cells+N_newcells,N_points+N_newpoints);
1908
1909   //vector used to redefine the new point IDs
1910   QVector <vtkIdType> OffSet(N_points);
1911
1912   //copy undead points
1913   vtkIdType dst_id_node=0;
1914   for (vtkIdType src_id_node = 0; src_id_node < N_points; src_id_node++) {//loop through src points
1915     if(src_id_node!=DeadNode)//if the node isn't dead, copy it
1916     {
1917       vec3_t x;
1918       src->GetPoints()->GetPoint(src_id_node, x.data());
1919       dst->GetPoints()->SetPoint(dst_id_node, x.data());
1920       copyNodeData(src, src_id_node, dst, dst_id_node);
1921       OffSet[src_id_node]=src_id_node-dst_id_node;
1922       dst_id_node++;
1923     }
1924     else
1925     {
1926       if(DebugLevel>0) cout<<"src_id_node="<<src_id_node<<" dst_id_node="<<dst_id_node<<endl;
1927     }
1928   };
1929   if(DebugLevel>10) {
1930     cout<<"DeadCells="<<DeadCells<<endl;
1931     cout<<"MutatedCells="<<MutatedCells<<endl;
1932     cout<<"MutilatedCells="<<MutilatedCells<<endl;
1933   }
1934   //Copy undead cells
1935   for (vtkIdType id_cell = 0; id_cell < src->GetNumberOfCells(); ++id_cell) {//loop through src cells
1936     if(!DeadCells.contains(id_cell))//if the cell isn't dead
1937     {
1938       vtkIdType src_N_pts, *src_pts;
1939       vtkIdType dst_N_pts, *dst_pts;
1940       src->GetCellPoints(id_cell, src_N_pts, src_pts);
1941
1942       vtkIdType type_cell = src->GetCellType(id_cell);
1943       if(DebugLevel>10) cout<<"-->id_cell="<<id_cell<<endl;
1944       if(DebugLevel>10) for(int i=0;i<src_N_pts;i++) cout<<"src_pts["<<i<<"]="<<src_pts[i]<<endl;
1945 //       src->GetCellPoints(id_cell, dst_N_pts, dst_pts);
1946       dst_N_pts=src_N_pts;
1947       dst_pts=new vtkIdType[dst_N_pts];
1948       if(MutatedCells.contains(id_cell))//mutated cell
1949       {
1950         if(DebugLevel>10) cout<<"processing mutated cell "<<id_cell<<endl;
1951         for(int i=0;i<src_N_pts;i++)
1952         {
1953           if(src_pts[i]==DeadNode) {
1954             if(DebugLevel>10) {
1955               cout<<"SnapPoint="<<SnapPoint<<endl;
1956               cout<<"OffSet[SnapPoint]="<<OffSet[SnapPoint]<<endl;
1957               cout<<"src_pts["<<i<<"]="<<src_pts[i]<<endl;
1958             }
1959             dst_pts[i]=SnapPoint-OffSet[SnapPoint];
1960           }
1961           else dst_pts[i]=src_pts[i]-OffSet[src_pts[i]];
1962         }
1963         if(DebugLevel>10) cout<<"--->dst_pts:"<<endl;
1964         if(DebugLevel>10) for(int i=0;i<dst_N_pts;i++) cout<<"dst_pts["<<i<<"]="<<dst_pts[i]<<endl;
1965
1966       }
1967       else if(MutilatedCells.contains(id_cell))//mutilated cell
1968       {
1969         if(DebugLevel>10) cout<<"processing mutilated cell "<<id_cell<<endl;
1970
1971         if(type_cell==VTK_QUAD) {
1972           type_cell=VTK_TRIANGLE;
1973           dst_N_pts-=1;
1974         }
1975         else {cout<<"FATAL ERROR: Unknown mutilated cell detected! It is not a quad! Potential xenomorph infestation!"<<endl;EG_BUG;}
1976         //merge points
1977         int j=0;
1978         for(int i=0;i<src_N_pts;i++)
1979         {
1980           if(src_pts[i]==SnapPoint) { dst_pts[j]=SnapPoint-OffSet[SnapPoint];j++; }//SnapPoint
1981           else if(src_pts[i]!=DeadNode) { dst_pts[j]=src_pts[i]-OffSet[src_pts[i]];j++; }//pre-snap/dead + post-snap/dead
1982           //do nothing in case of DeadNode
1983         }
1984       }
1985       else//normal cell
1986       {
1987         if(DebugLevel>10) cout<<"processing normal cell "<<id_cell<<endl;
1988         for(int i=0;i<src_N_pts;i++)
1989         {
1990           dst_pts[i]=src_pts[i]-OffSet[src_pts[i]];
1991         }
1992       }
1993       //copy the cell
1994       vtkIdType id_new_cell = dst->InsertNextCell(type_cell, dst_N_pts, dst_pts);
1995       copyCellData(src, id_cell, dst, id_new_cell);
1996       if(DebugLevel>10) {
1997         cout<<"===Copying cell "<<id_cell<<" to "<<id_new_cell<<"==="<<endl;
1998         cout<<"src_pts:"<<endl;
1999         for(int i=0;i<src_N_pts;i++) cout<<"src_pts["<<i<<"]="<<src_pts[i]<<endl;
2000         cout<<"dst_pts:"<<endl;
2001         for(int i=0;i<dst_N_pts;i++) cout<<"dst_pts["<<i<<"]="<<dst_pts[i]<<endl;
2002         cout<<"OffSet="<<OffSet<<endl;
2003         cout<<"===Copying cell end==="<<endl;
2004       }
2005       delete dst_pts;
2006     }
2007   };
2008 //   cout_grid(cout,dst,true,true,true,true);
2009   makeCopy(dst, src);
2010   return(true);
2011 }
2012 //End of DeletePoint_2
2013
2014 bool Operation::DeleteSetOfPoints(vtkUnstructuredGrid *src, QSet <vtkIdType> DeadNodes, int& N_newpoints, int& N_newcells)
2015 {
2016   QVector <vtkIdType> DeadNode_vector=Set2Vector(DeadNodes,false);
2017
2018   getAllSurfaceCells(cells,src);
2019 //   getNodesFromCells(cells, nodes, src);
2020   setGrid(src);
2021   setCells(cells);
2022   UpdateNodeType_all();
2023
2024   //src grid info
2025   int N_points=src->GetNumberOfPoints();
2026   int N_cells=src->GetNumberOfCells();
2027
2028   QSet <vtkIdType> DeadCells;
2029   QSet <vtkIdType> MutatedCells;
2030   QSet <vtkIdType> MutilatedCells;
2031   QVector <vtkIdType> SnapPoint(DeadNode_vector.size());
2032
2033     //counter init
2034   N_newpoints=0;
2035   N_newcells=0;
2036
2037   for(int i=0;i<DeadNode_vector.size();i++)
2038   {
2039     if(DeadNode_vector[i]<0 || DeadNode_vector[i]>=N_points)
2040     {
2041       cout<<"Warning: Point out of range: DeadNode_vector[i]="<<DeadNode_vector[i]<<" N_points="<<N_points<<endl;
2042       return(false);
2043     }
2044
2045     if(DebugLevel>10) {
2046       cout<<"BEFORE FINDSNAPPOINT"<<endl;
2047       cout<<"N_points="<<N_points<<endl;
2048       cout<<"N_cells="<<N_cells<<endl;
2049       cout<<"N_newpoints="<<N_newpoints<<endl;
2050       cout<<"N_newcells="<<N_newcells<<endl;
2051     }
2052
2053     //local values
2054     int l_N_newpoints;
2055     int l_N_newcells;
2056     QSet <vtkIdType> l_DeadCells;
2057     QSet <vtkIdType> l_MutatedCells;
2058     QSet <vtkIdType> l_MutilatedCells;
2059
2060     SnapPoint[i]=FindSnapPoint(src,DeadNode_vector[i], l_DeadCells, l_MutatedCells, l_MutilatedCells, l_N_newpoints, l_N_newcells);
2061     //global values
2062     N_newpoints+=l_N_newpoints;
2063     N_newcells+=l_N_newcells;
2064     DeadCells.unite(l_DeadCells);//DeadCells unite! Kill the living! :D
2065     MutatedCells.unite(l_MutatedCells);
2066     MutilatedCells.unite(l_MutilatedCells);
2067
2068     if(DebugLevel>0) cout<<"===>DeadNode_vector[i]="<<DeadNode_vector[i]<<" moving to SNAPPOINT="<<SnapPoint[i]<<" DebugLevel="<<DebugLevel<<endl;
2069     if(SnapPoint[i]<0) {cout<<"Sorry no possible SnapPoint found."<<endl; return(false);}
2070
2071   }
2072   //allocate
2073   if(DebugLevel>10) {
2074     cout<<"BEFORE ALLOCATION"<<endl;
2075     cout<<"N_points="<<N_points<<endl;
2076     cout<<"N_cells="<<N_cells<<endl;
2077     cout<<"N_newpoints="<<N_newpoints<<endl;
2078     cout<<"N_newcells="<<N_newcells<<endl;
2079   }
2080 //   N_points=src->GetNumberOfPoints();
2081 //   N_cells=src->GetNumberOfCells();
2082
2083   if(DebugLevel>47) {
2084     cout<<"N_points="<<N_points<<endl;
2085     cout<<"N_cells="<<N_cells<<endl;
2086     cout<<"N_newpoints="<<N_newpoints<<endl;
2087     cout<<"N_newcells="<<N_newcells<<endl;
2088   }
2089   EG_VTKSP(vtkUnstructuredGrid,dst);
2090   allocateGrid(dst,N_cells+N_newcells,N_points+N_newpoints);
2091
2092   //vector used to redefine the new point IDs
2093   QVector <vtkIdType> OffSet(N_points);
2094
2095   //copy undead points
2096   vtkIdType dst_id_node=0;
2097   for (vtkIdType src_id_node = 0; src_id_node < N_points; src_id_node++) {//loop through src points
2098     if(!DeadNode_vector.contains(src_id_node))//if the node isn't dead, copy it
2099     {
2100       vec3_t x;
2101       src->GetPoints()->GetPoint(src_id_node, x.data());
2102       dst->GetPoints()->SetPoint(dst_id_node, x.data());
2103       copyNodeData(src, src_id_node, dst, dst_id_node);
2104       OffSet[src_id_node]=src_id_node-dst_id_node;
2105       dst_id_node++;
2106     }
2107     else
2108     {
2109       if(DebugLevel>0) cout<<"src_id_node="<<src_id_node<<" dst_id_node="<<dst_id_node<<endl;
2110     }
2111   };
2112   if(DebugLevel>10) {
2113     cout<<"DeadCells="<<DeadCells<<endl;
2114     cout<<"MutatedCells="<<MutatedCells<<endl;
2115     cout<<"MutilatedCells="<<MutilatedCells<<endl;
2116   }
2117   //Copy undead cells
2118   for (vtkIdType id_cell = 0; id_cell < src->GetNumberOfCells(); ++id_cell) {//loop through src cells
2119     if(!DeadCells.contains(id_cell))//if the cell isn't dead
2120     {
2121       vtkIdType src_N_pts, *src_pts;
2122       vtkIdType dst_N_pts, *dst_pts;
2123       src->GetCellPoints(id_cell, src_N_pts, src_pts);
2124
2125       vtkIdType type_cell = src->GetCellType(id_cell);
2126       if(DebugLevel>10) cout<<"-->id_cell="<<id_cell<<endl;
2127       if(DebugLevel>10) for(int i=0;i<src_N_pts;i++) cout<<"src_pts["<<i<<"]="<<src_pts[i]<<endl;
2128 //       src->GetCellPoints(id_cell, dst_N_pts, dst_pts);
2129       dst_N_pts=src_N_pts;
2130       dst_pts=new vtkIdType[dst_N_pts];
2131       if(MutatedCells.contains(id_cell))//mutated cell
2132       {
2133         if(DebugLevel>10) cout<<"processing mutated cell "<<id_cell<<endl;
2134         for(int i=0;i<src_N_pts;i++)
2135         {
2136           int DeadIndex = DeadNode_vector.indexOf(src_pts[i]);
2137           if(DeadIndex!=-1) {
2138             if(DebugLevel>10) {
2139               cout<<"SnapPoint="<<SnapPoint[DeadIndex]<<endl;
2140               cout<<"OffSet[SnapPoint]="<<OffSet[SnapPoint[DeadIndex]]<<endl;
2141               cout<<"src_pts["<<i<<"]="<<src_pts[i]<<endl;
2142             }
2143             dst_pts[i]=SnapPoint[DeadIndex]-OffSet[SnapPoint[DeadIndex]];
2144           }
2145           else dst_pts[i]=src_pts[i]-OffSet[src_pts[i]];
2146         }
2147         if(DebugLevel>10) cout<<"--->dst_pts:"<<endl;
2148         if(DebugLevel>10) for(int i=0;i<dst_N_pts;i++) cout<<"dst_pts["<<i<<"]="<<dst_pts[i]<<endl;
2149
2150       }
2151       else if(MutilatedCells.contains(id_cell))//mutilated cell (ex: square becoming triangle)
2152       {
2153         cout<<"FATAL ERROR: Quads not supported yet."<<endl;EG_BUG;
2154
2155         if(DebugLevel>10) cout<<"processing mutilated cell "<<id_cell<<endl;
2156
2157         if(type_cell==VTK_QUAD) {
2158           type_cell=VTK_TRIANGLE;
2159           dst_N_pts-=1;
2160         }
2161         else {cout<<"FATAL ERROR: Unknown mutilated cell detected! It is not a quad! Potential xenomorph infestation!"<<endl;EG_BUG;}
2162         //merge points
2163         int j=0;
2164         for(int i=0;i<src_N_pts;i++)
2165         {
2166 /*          if(src_pts[i]==SnapPoint) { dst_pts[j]=SnapPoint-OffSet[SnapPoint];j++; }//SnapPoint
2167           else if(src_pts[i]!=DeadNode_vector[i]) { dst_pts[j]=src_pts[i]-OffSet[src_pts[i]];j++; }//pre-snap/dead + post-snap/dead*/
2168           //do nothing in case of DeadNode_vector[i]
2169         }
2170       }
2171       else//normal cell
2172       {
2173         if(DebugLevel>10) cout<<"processing normal cell "<<id_cell<<endl;
2174         for(int i=0;i<src_N_pts;i++)
2175         {
2176           dst_pts[i]=src_pts[i]-OffSet[src_pts[i]];
2177         }
2178       }
2179       //copy the cell
2180       vtkIdType id_new_cell = dst->InsertNextCell(type_cell, dst_N_pts, dst_pts);
2181       copyCellData(src, id_cell, dst, id_new_cell);
2182       if(DebugLevel>10) {
2183         cout<<"===Copying cell "<<id_cell<<" to "<<id_new_cell<<"==="<<endl;
2184         cout<<"src_pts:"<<endl;
2185         for(int i=0;i<src_N_pts;i++) cout<<"src_pts["<<i<<"]="<<src_pts[i]<<endl;
2186         cout<<"dst_pts:"<<endl;
2187         for(int i=0;i<dst_N_pts;i++) cout<<"dst_pts["<<i<<"]="<<dst_pts[i]<<endl;
2188         cout<<"OffSet="<<OffSet<<endl;
2189         cout<<"===Copying cell end==="<<endl;
2190       }
2191       delete dst_pts;
2192     }
2193   };
2194
2195 //   cout_grid(cout,dst,true,true,true,true);
2196   makeCopy(dst, src);
2197   return(true);
2198 }
2199 //End of DeleteSetOfPoints
2200
2201 void Operation::TxtSave(QString a_filename)
2202 {
2203   cout << a_filename.toAscii().data() << endl;
2204   ofstream file;
2205   file.open(a_filename.toAscii().data());
2206   cout_grid(file,grid,true,true,true,true);
2207   file.close();
2208 }
2209
2210 void Operation::DualSave(QString a_filename)
2211 {
2212   TxtSave(a_filename+".txt");
2213   GuiMainWindow::pointer()->QuickSave(a_filename+".vtu");
2214 }