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