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
  40 #include "geometrytools.h"
  41 using namespace GeometryTools;
  42
  43 #include <QApplication>
  44
  45 QSet<Operation*> Operation::garbage_operations;
  46
  47 void Operation::collectGarbage()
  48 {
  49   QSet<Operation*> delete_operations;
  50   foreach (Operation *op, garbage_operations) {
  51     if (!op->getThread().isRunning()) {
  52       delete_operations.insert(op);
  53       cout << "deleting Operation " << op << endl;
  54       delete op;
  55     };
  56   };
  57   foreach (Operation *op, delete_operations) {
  58     garbage_operations.remove(op);
  59   };
  60 };
  61
  62 Operation::Operation()
  63 {
  64   grid = NULL;
  65   gui = false;
  66   err = NULL;
  67   autoset = true;
  68 };
  69
  70 Operation::~Operation()
  71 {
  72   if (err) {
  73     err->display();
  74     delete err;
  75   };
  76 };
  77
  78 void Operation::del()
  79 {
  80   garbage_operations.insert(this);
  81 };
  82
  83 void OperationThread::run()
  84 {
  85   try {
  86     GuiMainWindow::lock();
  87     GuiMainWindow::pointer()->setBusy();
  88     op->operate();
  89   } catch (Error err) {
  90     op->err = new Error();
  91     *(op->err) = err;
  92   };
  93   GuiMainWindow::unlock();
  94   GuiMainWindow::pointer()->setIdle();
  95 };
  96
  97 void Operation::operator()()
  98 {
  99   if (gui) {
 100     if (GuiMainWindow::tryLock()) {
 101       checkGrid();
 102       thread.setOperation(this);
 103       GuiMainWindow::unlock();
 104       thread.start(QThread::LowPriority);
 105     } else {
 106       QMessageBox::warning(NULL, "not permitted", "Operation is not permitted while background process is running!");
 107     };
 108   } else {
 109     checkGrid();
 110     operate();
 111   };
 112   /*
 113   checkGrid();
 114   operate();
 115   */
 116 };
 117
 118 void Operation::setAllCells()
 119 {
 120   QVector<vtkIdType> all_cells;
 121   getAllCells(all_cells, grid);
 122   setCells(all_cells);
 123 };
 124
 125 void Operation::setAllVolumeCells()
 126 {
 127   QVector<vtkIdType> cells;
 128   getAllVolumeCells(cells, grid);
 129   setCells(cells);
 130 };
 131
 132 void Operation::setAllSurfaceCells()
 133 {
 134   QVector<vtkIdType> cells;
 135   getAllSurfaceCells(cells, grid);
 136   setCells(cells);
 137 };
 138
 139 void Operation::initMapping()
 140 {
 141   nodes_map.resize(nodes.size());
 142   for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
 143     nodes_map[i_nodes] = nodes[i_nodes];
 144   };
 145   cells_map.resize(cells.size());
 146   for (int i_cells = 0; i_cells < cells.size(); ++i_cells) {
 147     cells_map[i_cells] = cells[i_cells];
 148   };
 149 };
 150
 151 void Operation::checkGrid()
 152 {
 153   if (grid == NULL) {
 154     grid = GuiMainWindow::pointer()->getGrid();
 155   };
 156   if ((cells.size() == 0) && autoset) {
 157     setAllCells();
 158   };
 159 };
 160
 161 void Operation::updateActors()
 162 {
 163   mainWindow()->updateActors();
 164 };
 165
 166 GuiMainWindow* Operation::mainWindow()
 167 {
 168   return GuiMainWindow::pointer();
 169 };
 170
 171 void Operation::populateBoundaryCodes(QListWidget *lw)
 172 {
 173   QSet<int> bcs;
 174   mainWindow()->getAllBoundaryCodes(bcs);
 175   foreach(int bc, bcs) {
 176     QListWidgetItem *lwi = new QListWidgetItem(lw);
 177     lwi->setCheckState(Qt::Unchecked);
 178     QString text = "";
 179     QTextStream ts(&text);
 180     ts << bc;
 181     lwi->setText(text);
 182     lwi->setFlags(Qt::ItemIsUserCheckable | Qt::ItemIsEnabled);
 183   };
 184 };
 185
 186 stencil_t Operation::getStencil(vtkIdType id_cell1, int j1)
 187 {
 188   stencil_t S;
 189   S.valid = true;
 190   S.id_cell1 = id_cell1;
 191   if (c2c[_cells[id_cell1]][j1] != -1) {
 192     S.id_cell2 = cells[c2c[_cells[id_cell1]][j1]];
 193     if (grid->GetCellType(S.id_cell2) != VTK_TRIANGLE) {
 194       EG_BUG;
 195     };
 196     vtkIdType N1, N2, *pts1, *pts2;
 197     grid->GetCellPoints(S.id_cell1, N1, pts1);
 198     grid->GetCellPoints(S.id_cell2, N2, pts2);
 199     if      (j1 == 0) { S.p[0] = pts1[2]; S.p[1] = pts1[0]; S.p[3] = pts1[1]; }
 200     else if (j1 == 1) { S.p[0] = pts1[0]; S.p[1] = pts1[1]; S.p[3] = pts1[2]; }
 201     else if (j1 == 2) { S.p[0] = pts1[1]; S.p[1] = pts1[2]; S.p[3] = pts1[0]; };
 202     bool p2 = false;
 203     if (c2c[_cells[S.id_cell2]][0] != -1) {
 204       if (cells[c2c[_cells[S.id_cell2]][0]] == S.id_cell1) {
 205         S.p[2] = pts2[2];
 206         p2 = true;
 207       };
 208     };
 209     if (c2c[_cells[S.id_cell2]][1] != -1) {
 210       if (cells[c2c[_cells[S.id_cell2]][1]] == S.id_cell1) {
 211         S.p[2] = pts2[0];
 212         p2 = true;
 213       };
 214     };
 215     if (c2c[_cells[S.id_cell2]][2] != -1) {
 216       if (cells[c2c[_cells[S.id_cell2]][2]] == S.id_cell1) {
 217         S.p[2] = pts2[1];
 218         p2 = true;
 219       };
 220     };
 221     if (!p2) {
 222       EG_BUG;
 223     };
 224   } else {
 225     S.valid = false;
 226     S.id_cell2 = -1;
 227     vtkIdType N1, *pts1;
 228     grid->GetCellPoints(S.id_cell1, N1, pts1);
 229     if      (j1 == 0) { S.p[0] = pts1[2]; S.p[1] = pts1[0]; S.p[3] = pts1[1]; }
 230     else if (j1 == 1) { S.p[0] = pts1[0]; S.p[1] = pts1[1]; S.p[3] = pts1[2]; }
 231     else if (j1 == 2) { S.p[0] = pts1[1]; S.p[1] = pts1[2]; S.p[3] = pts1[0]; };
 232   };
 233   return S;
 234 };
 235
 236 ostream& operator<<(ostream &out, stencil_t S)
 237 {
 238   out<<"S.id_cell1="<<S.id_cell1<<" ";
 239   out<<"S.id_cell2="<<S.id_cell2<<" ";
 240   out<<"S.valid="<<S.valid<<" ";
 241   out<<"[";
 242   for(int i=0;i<4;i++){
 243     out<<S.p[i];
 244     if(i!=3) out<<",";
 245   }
 246   out<<"]";
 247   return(out);
 248 }
 249
 250 //////////////////////////////////////////////
 251 double CurrentVertexAvgDist(vtkIdType a_vertex,QVector< QSet< int > > &n2n,vtkUnstructuredGrid *a_grid)
 252 {
 253   double total_dist=0;
 254   double avg_dist=0;
 255   int N=n2n[a_vertex].size();
 256   vec3_t C;
 257   a_grid->GetPoint(a_vertex, C.data());
 258   foreach(int i,n2n[a_vertex])
 259   {
 260     vec3_t M;
 261     a_grid->GetPoint(i, M.data());
 262     total_dist+=(M-C).abs();
 263   }
 264   avg_dist=total_dist/(double)N;
 265   return(avg_dist);
 266 }
 267
 268 double CurrentMeshDensity(vtkIdType a_vertex,QVector< QSet< int > > &n2n,vtkUnstructuredGrid *a_grid)
 269 {
 270   double total_dist=0;
 271   double avg_dist=0;
 272   int N=n2n[a_vertex].size();
 273   vec3_t C;
 274   a_grid->GetPoint(a_vertex, C.data());
 275   foreach(int i,n2n[a_vertex])
 276   {
 277     vec3_t M;
 278     a_grid->GetPoint(i, M.data());
 279     total_dist+=(M-C).abs();
 280   }
 281   avg_dist=total_dist/(double)N;
 282   double avg_density=1./avg_dist;
 283   return(avg_density);
 284 }
 285
 286 double DesiredVertexAvgDist(vtkIdType a_vertex,QVector< QSet< int > > &n2n,vtkUnstructuredGrid *a_grid)
 287 {
 288   double total_dist=0;
 289   double avg_dist=0;
 290   EG_VTKDCN(vtkDoubleArray, node_meshdensity, a_grid, "node_meshdensity");
 291   int N=n2n[a_vertex].size();
 292   foreach(int i,n2n[a_vertex])
 293   {
 294     total_dist+=1./node_meshdensity->GetValue(i);
 295   }
 296   avg_dist=total_dist/(double)N;
 297   return(avg_dist);
 298 }
 299
 300 double DesiredMeshDensity(vtkIdType a_vertex,QVector< QSet< int > > &n2n,vtkUnstructuredGrid *a_grid)
 301 {
 302   double total_density=0;
 303   double avg_density=0;
 304   EG_VTKDCN(vtkDoubleArray, node_meshdensity, a_grid, "node_meshdensity");
 305   int N=n2n[a_vertex].size();
 306   foreach(int i,n2n[a_vertex])
 307   {
 308     total_density+=node_meshdensity->GetValue(i);
 309   }
 310   avg_density=total_density/(double)N;
 311   return(avg_density);
 312 }
 313
 314 ///////////////////////////////////////////
 315 vtkIdType Operation::getClosestNode(vtkIdType a_id_node,vtkUnstructuredGrid* a_grid)
 316 {
 317   vec3_t C;
 318   a_grid->GetPoint(a_id_node,C.data());
 319   vtkIdType id_minlen=-1;
 320   double minlen=-1;
 321   foreach(vtkIdType neighbour,n2n[a_id_node])
 322   {
 323     vec3_t M;
 324     a_grid->GetPoint(neighbour,M.data());
 325     double len=(M-C).abs();
 326     if(minlen<0 or len<minlen)
 327     {
 328       minlen=len;
 329       id_minlen=neighbour;
 330     }
 331   }
 332   return(id_minlen);
 333 }
 334
 335 vtkIdType Operation::getFarthestNode(vtkIdType a_id_node,vtkUnstructuredGrid* a_grid)
 336 {
 337   vec3_t C;
 338   a_grid->GetPoint(a_id_node,C.data());
 339   vtkIdType id_maxlen=-1;
 340   double maxlen=-1;
 341   foreach(vtkIdType neighbour,n2n[a_id_node])
 342   {
 343     vec3_t M;
 344     a_grid->GetPoint(neighbour,M.data());
 345     double len=(M-C).abs();
 346     if(maxlen<0 or len>maxlen)
 347     {
 348       maxlen=len;
 349       id_maxlen=neighbour;
 350     }
 351   }
 352   return(id_maxlen);
 353 }
 354
 355 bool Operation::SwapCells(vtkUnstructuredGrid* a_grid, stencil_t S)
 356 {
 357   bool swap = false;
 358   if (S.valid) {
 359     vec3_t x3[4], x3_0(0,0,0);
 360     vec2_t x[4];
 361     for (int k = 0; k < 4; ++k) {
 362       a_grid->GetPoints()->GetPoint(S.p[k], x3[k].data());
 363       x3_0 += x3[k];
 364     };
 365     vec3_t n1 = triNormal(x3[0], x3[1], x3[3]);
 366     vec3_t n2 = triNormal(x3[1], x3[2], x3[3]);
 367     n1.normalise();
 368     n2.normalise();
 369     if ( (n1*n2) > 0.8) {
 370       vec3_t n = n1 + n2;
 371       n.normalise();
 372       vec3_t ex = orthogonalVector(n);
 373       vec3_t ey = ex.cross(n);
 374       for (int k = 0; k < 4; ++k) {
 375         x[k] = vec2_t(x3[k]*ex, x3[k]*ey);
 376       };
 377       vec2_t r1, r2, r3, u1, u2, u3;
 378       r1 = 0.5*(x[0] + x[1]); u1 = turnLeft(x[1] - x[0]);
 379       r2 = 0.5*(x[1] + x[2]); u2 = turnLeft(x[2] - x[1]);
 380       r3 = 0.5*(x[1] + x[3]); u3 = turnLeft(x[3] - x[1]);
 381       double k, l;
 382       vec2_t xm1, xm2;
 383       bool ok = true;
 384       if (intersection(k, l, r1, u1, r3, u3)) {
 385         xm1 = r1 + k*u1;
 386         if (intersection(k, l, r2, u2, r3, u3)) {
 387           xm2 = r2 + k*u2;
 388         } else {
 389           ok = false;
 390         };
 391       } else {
 392         ok = false;
 393         swap = true;
 394       };
 395       if (ok) {
 396         if ((xm1 - x[2]).abs() < (xm1 - x[0]).abs()) {
 397           swap = true;
 398         };
 399         if ((xm2 - x[0]).abs() < (xm2 - x[2]).abs()) {
 400           swap = true;
 401         };
 402       };
 403     };
 404   };
 405   if (swap) {
 406     vtkIdType new_pts1[3], new_pts2[3];
 407     new_pts1[0] = S.p[1];
 408     new_pts1[1] = S.p[2];
 409     new_pts1[2] = S.p[0];
 410     new_pts2[0] = S.p[2];
 411     new_pts2[1] = S.p[3];
 412     new_pts2[2] = S.p[0];
 413     a_grid->ReplaceCell(S.id_cell1, 3, new_pts1);
 414     a_grid->ReplaceCell(S.id_cell2, 3, new_pts2);
 415   };
 416   return(swap);
 417 }
 418
 419 void Operation::quad2triangle(vtkUnstructuredGrid* src,vtkIdType quadcell)
 420 {
 421   vtkIdType type_cell = grid->GetCellType(quadcell);
 422   cout<<"It's a "<<type_cell<<endl;
 423   if(type_cell==VTK_QUAD)
 424   {
 425     cout_grid(cout,src,true,true,true,true);
 426     EG_VTKSP(vtkUnstructuredGrid, dst);
 427     //src grid info
 428     int N_points=src->GetNumberOfPoints();
 429     int N_cells=src->GetNumberOfCells();
 430     allocateGrid(dst,N_cells+1,N_points);
 431
 432     for (vtkIdType id_node = 0; id_node < src->GetNumberOfPoints(); ++id_node) {
 433       vec3_t x;
 434       src->GetPoints()->GetPoint(id_node, x.data());
 435       dst->GetPoints()->SetPoint(id_node, x.data());
 436       copyNodeData(src, id_node, dst, id_node);
 437     };
 438     for (vtkIdType id_cell = 0; id_cell < src->GetNumberOfCells(); ++id_cell) {
 439       vtkIdType N_pts, *pts;
 440       vtkIdType type_cell = src->GetCellType(id_cell);
 441       src->GetCellPoints(id_cell, N_pts, pts);
 442       vtkIdType id_new_cell;
 443       vtkIdType id_new_cell1;
 444       vtkIdType id_new_cell2;
 445       if(id_cell!=quadcell)
 446       {
 447         id_new_cell = dst->InsertNextCell(type_cell, N_pts, pts);
 448         copyCellData(src, id_cell, dst, id_new_cell);
 449       }
 450       else
 451       {
 452         vtkIdType triangle1[3];
 453         vtkIdType triangle2[3];
 454         triangle1[0]=pts[1];
 455         triangle1[1]=pts[3];
 456         triangle1[2]=pts[0];
 457         triangle2[0]=pts[3];
 458         triangle2[1]=pts[1];
 459         triangle2[2]=pts[2];
 460         id_new_cell1 = dst->InsertNextCell(VTK_TRIANGLE, 3, triangle1);
 461         copyCellData(src, id_cell, dst, id_new_cell1);
 462         id_new_cell2 = dst->InsertNextCell(VTK_TRIANGLE, 3, triangle2);
 463         copyCellData(src, id_cell, dst, id_new_cell2);
 464         stencil_t S;
 465         S.id_cell1=id_new_cell1;
 466         S.id_cell2=id_new_cell2;
 467         S.p[0]=pts[0];
 468         S.p[1]=pts[1];
 469         S.p[2]=pts[2];
 470         S.p[3]=pts[3];
 471         S.valid=true;
 472         SwapCells(dst,S);
 473       }
 474     };
 475     cout_grid(cout,dst,true,true,true,true);
 476     makeCopy(dst, src);
 477   }//end of if quad
 478 }
 479
 480 void Operation::quad2triangle(vtkUnstructuredGrid* src,vtkIdType quadcell,vtkIdType MovingPoint)
 481 {
 482   vtkIdType type_cell = grid->GetCellType(quadcell);
 483   cout<<"It's a "<<type_cell<<endl;
 484   if(type_cell==VTK_QUAD)
 485   {
 486     cout_grid(cout,src,true,true,true,true);
 487     EG_VTKSP(vtkUnstructuredGrid, dst);
 488     //src grid info
 489     int N_points=src->GetNumberOfPoints();
 490     int N_cells=src->GetNumberOfCells();
 491     allocateGrid(dst,N_cells+1,N_points);
 492
 493     for (vtkIdType id_node = 0; id_node < src->GetNumberOfPoints(); ++id_node) {
 494       vec3_t x;
 495       src->GetPoints()->GetPoint(id_node, x.data());
 496       dst->GetPoints()->SetPoint(id_node, x.data());
 497       copyNodeData(src, id_node, dst, id_node);
 498     };
 499     for (vtkIdType id_cell = 0; id_cell < src->GetNumberOfCells(); ++id_cell) {
 500       vtkIdType N_pts, *pts;
 501       src->GetCellPoints(id_cell, N_pts, pts);
 502       vtkIdType type_cell = src->GetCellType(id_cell);
 503       vtkIdType id_new_cell;
 504       vtkIdType id_new_cell1;
 505       vtkIdType id_new_cell2;
 506       if(id_cell!=quadcell)
 507       {
 508         id_new_cell = dst->InsertNextCell(type_cell, N_pts, pts);
 509         copyCellData(src, id_cell, dst, id_new_cell);
 510       }
 511       else
 512       {
 513         vtkIdType triangle1[3];
 514         vtkIdType triangle2[3];
 515         if(MovingPoint==pts[0] || MovingPoint==pts[2])
 516         {
 517           triangle1[0]=pts[1];
 518           triangle1[1]=pts[3];
 519           triangle1[2]=pts[0];
 520           triangle2[0]=pts[3];
 521           triangle2[1]=pts[1];
 522           triangle2[2]=pts[2];
 523         }
 524         else
 525         {
 526           triangle1[0]=pts[2];
 527           triangle1[1]=pts[0];
 528           triangle1[2]=pts[1];
 529           triangle2[0]=pts[0];
 530           triangle2[1]=pts[2];
 531           triangle2[2]=pts[3];
 532         }
 533         id_new_cell1 = dst->InsertNextCell(VTK_TRIANGLE, 3, triangle1);
 534         copyCellData(src, id_cell, dst, id_new_cell1);
 535         id_new_cell2 = dst->InsertNextCell(VTK_TRIANGLE, 3, triangle2);
 536         copyCellData(src, id_cell, dst, id_new_cell2);
 537       }
 538     };
 539     cout_grid(cout,dst,true,true,true,true);
 540     makeCopy(dst, src);
 541   }//end of if quad
 542 }
 543
 544 int Operation::NumberOfCommonPoints(vtkIdType node1, vtkIdType node2, bool& IsTetra)
 545 {
 546 //   QVector< QSet< int > >     n2n
 547   QSet <int> node1_neighbours=n2n[node1];
 548   QSet <int> node2_neighbours=n2n[node2];
 549   QSet <int> intersection=node1_neighbours.intersect(node2_neighbours);
 550   int N=intersection.size();
 551   IsTetra=false;
 552   if(N==2)
 553   {
 554     QSet<int>::const_iterator p1=intersection.begin();
 555     QSet<int>::const_iterator p2=p1+1;
 556     vtkIdType intersection1=_nodes[*p1];
 557     vtkIdType intersection2=_nodes[*p2];
 558     if(n2n[intersection1].contains(intersection2))//if there's an edge between intersection1 and intersection2
 559     {
 560       //check if (node1,intersection1,intersection2) and (node2,intersection1,intersection2) are defined as cells!
 561   //     QVector< QSet< int > >         n2c
 562       QSet< int > S1=n2c[intersection1];
 563       QSet< int > S2=n2c[intersection2];
 564       QSet< int > Si=S1.intersect(S2);
 565       int counter=0;
 566       foreach(vtkIdType C,Si){
 567         vtkIdType N_pts, *pts;
 568         grid->GetCellPoints(C, N_pts, pts);
 569         for(int i=0;i<N_pts;i++)
 570         {
 571           if(pts[i]==node1 || pts[i]==node2) counter++;
 572         }
 573       }
 574       if(counter>=2) IsTetra=true;
 575     }
 576   }
 577   return(N);
 578 }
 579
 580 // vtkIdType Operation::FindSnapPoint(vtkUnstructuredGrid *src, vtkIdType DeadNode)
 581 // {
 582 //   return(0);
 583 // }
 584
 585 bool Operation::DeletePoint(vtkUnstructuredGrid *src, vtkIdType DeadNode)
 586 {
 587   EG_VTKSP(vtkUnstructuredGrid, dst);
 588
 589     //src grid info
 590   int N_points=src->GetNumberOfPoints();
 591   int N_cells=src->GetNumberOfCells();
 592   int N_newpoints=-1;
 593   int N_newcells=0;
 594
 595   if(DeadNode<0 || DeadNode>=N_points)
 596   {
 597     cout<<"Warning: Point out of range: DeadNode="<<DeadNode<<" N_points="<<N_points<<endl;
 598     return(false);
 599   }
 600
 601   QSet <vtkIdType> DeadCells;
 602   QSet <vtkIdType> MutatedCells;
 603   QSet <vtkIdType> MutilatedCells;
 604
 605   vtkIdType SnapPoint=-1;
 606   //Find closest point to DeadNode
 607 //   vtkIdType SnapPoint = getClosestNode(DeadNode,src);//DeadNode moves to SnapPoint
 608
 609   foreach(vtkIdType PSP, n2n[DeadNode])
 610   {
 611     bool IsValidSnapPoint=true;
 612
 613     cout<<"====>PSP="<<PSP<<endl;
 614     bool IsTetra=true;
 615     if(NumberOfCommonPoints(DeadNode,PSP,IsTetra)>2)//common point check
 616     {
 617       cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
 618       IsValidSnapPoint=false;
 619     }
 620     if(IsTetra)//tetra check
 621     {
 622       cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
 623       IsValidSnapPoint=false;
 624     }
 625
 626     //count number of points and cells to remove + analyse cell transformations
 627     N_newpoints=-1;
 628     N_newcells=0;
 629     DeadCells.clear();
 630     MutatedCells.clear();
 631     MutilatedCells.clear();
 632     foreach(vtkIdType C, n2c[DeadNode])//loop through potentially dead cells
 633     {
 634       cout<<"C="<<C<<endl;
 635       //get points around cell
 636       vtkIdType N_pts, *pts;
 637       src->GetCellPoints(C, N_pts, pts);
 638
 639       bool ContainsSnapPoint=false;
 640       bool invincible=false;
 641       for(int i=0;i<N_pts;i++)
 642       {
 643         cout<<"pts["<<i<<"]="<<pts[i]<<" and PSP="<<PSP<<endl;
 644         if(pts[i]==PSP) {ContainsSnapPoint=true;}
 645         if(pts[i]!=DeadNode && pts[i]!=PSP &&  n2c[pts[i]].size()<=1) invincible=true;
 646       }
 647       if(ContainsSnapPoint)
 648       {
 649         if(N_pts==3)//dead cell
 650         {
 651           if(invincible)//Check that empty lines aren't left behind when a cell is killed
 652           {
 653             cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
 654             IsValidSnapPoint=false;
 655           }
 656           else
 657           {
 658             DeadCells.insert(C);
 659             N_newcells-=1;
 660             cout<<"cell "<<C<<" has been pwned!"<<endl;
 661           }
 662         }
 663   /*      else if(N_pts==4)//mutilated cell
 664         {
 665           MutilatedCells.insert(C);
 666           cout<<"cell "<<C<<" has lost a limb!"<<endl;
 667         }*/
 668         else
 669         {
 670           cout<<"RED ALERT: Xenomorph detected!"<<endl;
 671           EG_BUG;
 672         }
 673       }
 674       else
 675       {
 676         vtkIdType src_N_pts, *src_pts;
 677         src->GetCellPoints(C, src_N_pts, src_pts);
 678
 679         if(src_N_pts!=3)
 680         {
 681           cout<<"RED ALERT: Xenomorph detected!"<<endl;
 682           EG_BUG;
 683         }
 684
 685         vtkIdType OldTriangle[3];
 686         vtkIdType NewTriangle[3];
 687
 688         for(int i=0;i<src_N_pts;i++)
 689         {
 690           OldTriangle[i]=src_pts[i];
 691           NewTriangle[i]=( (src_pts[i]==DeadNode) ? PSP : src_pts[i] );
 692         }
 693         vec3_t Old_N= triNormal(src, OldTriangle[0], OldTriangle[1], OldTriangle[2]);
 694         vec3_t New_N= triNormal(src, NewTriangle[0], NewTriangle[1], NewTriangle[2]);
 695         double OldArea=Old_N.abs();
 696         double NewArea=New_N.abs();
 697         double scal=Old_N*New_N;
 698         double cross=(Old_N.cross(New_N)).abs();//double-cross on Nar Shadaa B-)
 699
 700         cout<<"OldArea="<<OldArea<<endl;
 701         cout<<"NewArea="<<NewArea<<endl;
 702         cout<<"scal="<<scal<<endl;
 703         cout<<"cross="<<cross<<endl;
 704
 705         if(Old_N*New_N<Old_N*Old_N*1./100.)//area + inversion check
 706         {
 707           cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
 708           IsValidSnapPoint=false;
 709         }
 710
 711   /*      if(NewArea<OldArea*1./100.)
 712         {
 713           cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
 714           IsValidSnapPoint=false;
 715         }
 716
 717         if(abs(cross)>10e-4)
 718         {
 719           cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
 720           IsValidSnapPoint=false;
 721         }*/
 722
 723         //mutated cell
 724         MutatedCells.insert(C);
 725         cout<<"cell "<<C<<" has been infected!"<<endl;
 726       }
 727     }
 728
 729     if(N_cells+N_newcells<=0)//survivor check
 730     {
 731       cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
 732       IsValidSnapPoint=false;
 733     }
 734 /*    if(EmptyVolume(DeadNode,PSP))//simplified volume check
 735     {
 736       cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
 737       IsValidSnapPoint=false;
 738     }*/
 739     if(IsValidSnapPoint) {SnapPoint=PSP; break;}
 740   }//end of loop through potential SnapPoints
 741
 742   cout<<"===>SNAPPOINT="<<SnapPoint<<endl;
 743   if(SnapPoint<0) {cout<<"Sorry no possible SnapPoint found."<<endl; return(false);}
 744
 745   //allocate
 746   cout<<"N_points="<<N_points<<endl;
 747   cout<<"N_cells="<<N_cells<<endl;
 748   cout<<"N_newpoints="<<N_newpoints<<endl;
 749   cout<<"N_newcells="<<N_newcells<<endl;
 750   allocateGrid(dst,N_cells+N_newcells,N_points+N_newpoints);
 751
 752   //vector used to redefine the new point IDs
 753   QVector <vtkIdType> OffSet(N_points);
 754
 755   //copy undead points
 756   vtkIdType dst_id_node=0;
 757   for (vtkIdType src_id_node = 0; src_id_node < N_points; src_id_node++) {//loop through src points
 758     if(src_id_node!=DeadNode)//if the node isn't dead, copy it
 759     {
 760       vec3_t x;
 761       src->GetPoints()->GetPoint(src_id_node, x.data());
 762       dst->GetPoints()->SetPoint(dst_id_node, x.data());
 763       copyNodeData(src, src_id_node, dst, dst_id_node);
 764       OffSet[src_id_node]=src_id_node-dst_id_node;
 765       dst_id_node++;
 766     }
 767   };
 768
 769   cout<<"DeadCells="<<DeadCells<<endl;
 770   cout<<"MutatedCells="<<MutatedCells<<endl;
 771   cout<<"MutilatedCells="<<MutilatedCells<<endl;
 772
 773   //Copy undead cells
 774   for (vtkIdType id_cell = 0; id_cell < src->GetNumberOfCells(); ++id_cell) {//loop through src cells
 775     if(!DeadCells.contains(id_cell))//if the cell isn't dead
 776     {
 777       vtkIdType src_N_pts, *src_pts;
 778       vtkIdType dst_N_pts, *dst_pts;
 779       src->GetCellPoints(id_cell, src_N_pts, src_pts);
 780
 781       vtkIdType type_cell = src->GetCellType(id_cell);
 782       cout<<"-->id_cell="<<id_cell<<endl;
 783       for(int i=0;i<src_N_pts;i++) cout<<"src_pts["<<i<<"]="<<src_pts[i]<<endl;
 784 //       src->GetCellPoints(id_cell, dst_N_pts, dst_pts);
 785       dst_N_pts=src_N_pts;
 786       dst_pts=new vtkIdType[dst_N_pts];
 787       if(MutatedCells.contains(id_cell))//mutated cell
 788       {
 789         cout<<"processing mutated cell "<<id_cell<<endl;
 790         for(int i=0;i<src_N_pts;i++)
 791         {
 792           if(src_pts[i]==DeadNode) {
 793             cout<<"SnapPoint="<<SnapPoint<<endl;
 794             cout<<"OffSet[SnapPoint]="<<OffSet[SnapPoint]<<endl;
 795             cout<<"src_pts["<<i<<"]="<<src_pts[i]<<endl;
 796             dst_pts[i]=SnapPoint-OffSet[SnapPoint];
 797           }
 798           else dst_pts[i]=src_pts[i]-OffSet[src_pts[i]];
 799         }
 800         cout<<"--->dst_pts:"<<endl;
 801         for(int i=0;i<dst_N_pts;i++) cout<<"dst_pts["<<i<<"]="<<dst_pts[i]<<endl;
 802
 803       }
 804       else if(MutilatedCells.contains(id_cell))//mutilated cell
 805       {
 806         cout<<"processing mutilated cell "<<id_cell<<endl;
 807
 808         if(type_cell==VTK_QUAD) {
 809           type_cell=VTK_TRIANGLE;
 810           dst_N_pts-=1;
 811         }
 812         else {cout<<"FATAL ERROR: Unknown mutilated cell detected! It is not a quad! Potential xenomorph infestation!"<<endl;EG_BUG;}
 813         //merge points
 814         int j=0;
 815         for(int i=0;i<src_N_pts;i++)
 816         {
 817           if(src_pts[i]==SnapPoint) { dst_pts[j]=SnapPoint-OffSet[SnapPoint];j++; }//SnapPoint
 818           else if(src_pts[i]!=DeadNode) { dst_pts[j]=src_pts[i]-OffSet[src_pts[i]];j++; }//pre-snap/dead + post-snap/dead
 819           //do nothing in case of DeadNode
 820         }
 821       }
 822       else//normal cell
 823       {
 824         cout<<"processing normal cell "<<id_cell<<endl;
 825         for(int i=0;i<src_N_pts;i++)
 826         {
 827           dst_pts[i]=src_pts[i]-OffSet[src_pts[i]];
 828         }
 829       }
 830       //copy the cell
 831       vtkIdType id_new_cell = dst->InsertNextCell(type_cell, dst_N_pts, dst_pts);
 832       copyCellData(src, id_cell, dst, id_new_cell);
 833       cout<<"===Copying cell "<<id_cell<<" to "<<id_new_cell<<"==="<<endl;
 834       cout<<"src_pts:"<<endl;
 835       for(int i=0;i<src_N_pts;i++) cout<<"src_pts["<<i<<"]="<<src_pts[i]<<endl;
 836       cout<<"dst_pts:"<<endl;
 837       for(int i=0;i<dst_N_pts;i++) cout<<"dst_pts["<<i<<"]="<<dst_pts[i]<<endl;
 838       cout<<"OffSet="<<OffSet<<endl;
 839       cout<<"===Copying cell end==="<<endl;
 840       delete dst_pts;
 841     }
 842   };
 843 //   cout_grid(cout,dst,true,true,true,true);
 844   makeCopy(dst, src);
 845   return(true);
 846 }
 847
 848 bool Operation::EmptyVolume(vtkIdType DeadNode, vtkIdType PSP)
 849 {
 850   c2c[DeadNode];
 851   c2c[PSP];
 852   return(true);
 853 }
 854
 855 vec3_t Operation::GetCenter(vtkIdType cellId, double& R)
 856 {
 857   vtkIdType *pts, Npts;
 858   grid->GetCellPoints(cellId, Npts, pts);
 859
 860   vec3_t x(0,0,0);
 861   for (vtkIdType i = 0; i < Npts; ++i) {
 862     vec3_t xp;
 863     grid->GetPoints()->GetPoint(pts[i], xp.data());
 864     x += double(1)/Npts * xp;
 865   };
 866
 867   R = 1e99;
 868   for (vtkIdType i = 0; i < Npts; ++i) {
 869     vec3_t xp;
 870     grid->GetPoints()->GetPoint(pts[i], xp.data());
 871     R = min(R, 0.25*(xp-x).abs());
 872   };
 873
 874   return(x);
 875 }
 876
 877 bool Operation::getNeighbours(vtkIdType Boss, QVector <vtkIdType>& Peons, int BC)
 878 {
 879 //   QVector <vtkIdType> Peons;
 880
 881   QSet <int> S1=n2c[Boss];
 882   cout<<"S1="<<S1<<endl;
 883   foreach(vtkIdType PN,n2n[Boss])
 884   {
 885     cout<<"PN="<<PN<<endl;
 886     QSet <int> S2=n2c[PN];
 887     cout<<"S2="<<S2<<endl;
 888     QSet <int> Si=S2.intersect(S1);
 889     cout<<"PN="<<PN<<" Si="<<Si<<endl;
 890     if(Si.size()<2)//only one common cell
 891     {
 892       Peons.push_back(PN);
 893     }
 894     else
 895     {
 896       QSet <int> bc_set;
 897       foreach(vtkIdType C,Si)
 898       {
 899         EG_VTKDCC(vtkIntArray, cell_code, grid, "cell_code");
 900         int bc=cell_code->GetValue(C);
 901         cout<<"C="<<C<<" bc="<<bc<<endl;
 902         bc_set.insert(bc);
 903       }
 904       if(bc_set.size()>1)//2 different boundary codes
 905       {
 906         Peons.push_back(PN);
 907       }
 908     }
 909   }
 910   if(Peons.size()==2)
 911   {
 912 /*    Peon1=Peons[0];
 913     Peon2=Peons[1];*/
 914     return(true);
 915   }
 916   else
 917   {
 918     int N=n2n[Boss].size();
 919     QVector <vtkIdType> neighbours(N);
 920     qCopy(n2n[Boss].begin(), n2n[Boss].end(), neighbours.begin());
 921
 922     double alphamin_value;
 923     vtkIdType alphamin_i;
 924     vtkIdType alphamin_j;
 925     bool first=true;
 926
 927     for(int i=0;i<N;i++)
 928     {
 929       for(int j=i+1;j<N;j++)
 930       {
 931         double alpha=deviation(grid,neighbours[i],Boss,neighbours[j]);
 932 //         cout<<"alpha("<<neighbours[i]<<","<<Boss<<","<<neighbours[j]<<")="<<alpha<<endl;
 933         if(first) {
 934           alphamin_value=alpha;
 935           alphamin_i=i;
 936           alphamin_j=j;
 937           first=false;
 938         }
 939         else
 940         {
 941           if(alpha<alphamin_value)
 942           {
 943             alphamin_value=alpha;
 944             alphamin_i=i;
 945             alphamin_j=j;
 946           }
 947         }
 948       }
 949     }
 950 //     cout<<"alphamin_value="<<alphamin_value<<endl;
 951
 952     Peons.resize(2);
 953     Peons[0]=neighbours[alphamin_i];
 954     Peons[1]=neighbours[alphamin_j];
 955     return(true);
 956 /*    cout<<"FATAL ERROR: number of neighbours != 2"<<endl;
 957     EG_BUG;*/
 958   }
 959   return(false);
 960 }
 961
 962 int Operation::UpdateMeshDensity()
 963 {
 964   cout<<"===UpdateMeshDensity START==="<<endl;
 965
 966   getAllSurfaceCells(cells,grid);
 967   EG_VTKDCC(vtkIntArray, cell_code, grid, "cell_code");
 968   EG_VTKDCN(vtkDoubleArray, node_meshdensity, grid, "node_meshdensity");
 969   getNodesFromCells(cells, nodes, grid);
 970   setGrid(grid);
 971   setCells(cells);
 972
 973   if(DebugLevel>5) cout<<"cells.size()="<<cells.size()<<endl;
 974
 975   EG_VTKDCN(vtkDoubleArray, node_meshdensity_current, grid, "node_meshdensity_current");
 976   foreach(vtkIdType node,nodes)
 977   {
 978     double L=CurrentVertexAvgDist(node,n2n,grid);
 979     double D=1./L;
 980     node_meshdensity_current->SetValue(node, D);
 981   }
 982   cout<<"===UpdateMeshDensity END==="<<endl;
 983   return(0);
 984 }
 985
 986 // Special structure for marking vertices
 987 typedef struct _vtkMeshVertex
 988 {
 989   char      type;
 990   vtkIdList *edges; // connected edges (list of connected point ids)
 991 } vtkMeshVertex, *vtkMeshVertexPtr;
 992
 993 int Operation::UpdateNodeType()
 994 {
 995   if(DebugLevel>47) cout<<"this->FeatureAngle="<<this->FeatureAngle<<endl;
 996   if(DebugLevel>47) cout<<"this->EdgeAngle="<<this->EdgeAngle<<endl;
 997 //   cout<<"===UpdateNodeType START==="<<endl;
 998
 999   getAllSurfaceCells(cells,grid);
1000   if(DebugLevel>5) cout<<"cells.size()="<<cells.size()<<endl;
1001
1002   EG_VTKSP(vtkPolyData, pdata);
1003   //   addToPolyData(m_SelectedCells, pdata, grid);
1004   addToPolyData(cells, pdata, grid);
1005
1006   vtkPolyData* input=pdata;
1007
1008   vtkPolyData *source = 0;
1009
1010   vtkIdType numPts, numCells, i, numPolys;
1011   int j, k;
1012   vtkIdType npts = 0;
1013   vtkIdType *pts = 0;
1014   vtkIdType p1, p2;
1015   double x[3], y[3], deltaX[3], xNew[3], conv, maxDist, dist, factor;
1016   double x1[3], x2[3], x3[3], l1[3], l2[3];
1017   double CosFeatureAngle; //Cosine of angle between adjacent polys
1018   double CosEdgeAngle; // Cosine of angle between adjacent edges
1019   double closestPt[3], dist2, *w = NULL;
1020   int iterationNumber, abortExecute;
1021   vtkIdType numSimple=0, numBEdges=0, numFixed=0, numFEdges=0;
1022   vtkPolyData *inMesh, *Mesh;
1023   vtkPoints *inPts;
1024   vtkCellArray *inVerts, *inLines, *inPolys;
1025   vtkPoints *newPts;
1026   vtkMeshVertexPtr Verts;
1027   vtkCellLocator *cellLocator=NULL;
1028
1029     // Check input
1030     //
1031   numPts=input->GetNumberOfPoints();
1032   numCells=input->GetNumberOfCells();
1033   if (numPts < 1 || numCells < 1)
1034   {
1035     cout<<"No data to smooth!"<<endl;
1036     return 1;
1037   }
1038
1039   CosFeatureAngle =
1040     cos((double) vtkMath::DegreesToRadians() * this->FeatureAngle);
1041   CosEdgeAngle = cos((double) vtkMath::DegreesToRadians() * this->EdgeAngle);
1042
1043   if(DebugLevel>5) {
1044     cout<<"Smoothing " << numPts << " vertices, " << numCells
1045       << " cells with:\n"
1046       << "\tConvergence= " << this->Convergence << "\n"
1047       << "\tIterations= " << this->NumberOfIterations << "\n"
1048       << "\tRelaxation Factor= " << this->RelaxationFactor << "\n"
1049       << "\tEdge Angle= " << this->EdgeAngle << "\n"
1050       << "\tBoundary Smoothing " << (this->BoundarySmoothing ? "On\n" : "Off\n")
1051       << "\tFeature Edge Smoothing " << (this->FeatureEdgeSmoothing ? "On\n" : "Off\n")
1052       << "\tError Scalars " << (this->GenerateErrorScalars ? "On\n" : "Off\n")
1053       << "\tError Vectors " << (this->GenerateErrorVectors ? "On\n" : "Off\n")<<endl;
1054   }
1055     // Peform topological analysis. What we're gonna do is build a connectivity
1056     // array of connected vertices. The outcome will be one of three
1057     // classifications for a vertex: VTK_SIMPLE_VERTEX, VTK_FIXED_VERTEX. or
1058     // VTK_EDGE_VERTEX. Simple vertices are smoothed using all connected
1059     // vertices. FIXED vertices are never smoothed. Edge vertices are smoothed
1060     // using a subset of the attached vertices.
1061     //
1062   if(DebugLevel>5) cout<<"===>Analyze topology==="<<endl;
1063   if(DebugLevel>5) cout<<"Analyzing topology..."<<endl;
1064   if(DebugLevel>5) cout<<"0:numPts="<<numPts<<endl;
1065   Verts = new vtkMeshVertex[numPts];
1066   for (i=0; i<numPts; i++)
1067   {
1068     if(DebugLevel>5) cout<<"0:VTK_SIMPLE_VERTEX"<<endl;
1069     Verts[i].type = VTK_SIMPLE_VERTEX; //can smooth
1070     Verts[i].edges = NULL;
1071   }
1072
1073   inPts = input->GetPoints();
1074   conv = this->Convergence * input->GetLength();
1075
1076   if(DebugLevel>5) cout<<"==polygons and triangle strips=="<<endl;
1077     // now polygons and triangle strips-------------------------------
1078   inPolys=input->GetPolys();
1079   numPolys = inPolys->GetNumberOfCells();
1080
1081   if(DebugLevel>5) cout<<"numPolys="<<numPolys<<endl;
1082
1083   if ( numPolys > 0 )
1084   { //build cell structure
1085     vtkCellArray *polys;
1086     vtkIdType cellId;
1087     int numNei, nei, edge;
1088     vtkIdType numNeiPts;
1089     vtkIdType *neiPts;
1090     double normal[3], neiNormal[3];
1091     vtkIdList *neighbors;
1092
1093     neighbors = vtkIdList::New();
1094     neighbors->Allocate(VTK_CELL_SIZE);
1095
1096     inMesh = vtkPolyData::New();
1097     inMesh->SetPoints(inPts);
1098     inMesh->SetPolys(inPolys);
1099     Mesh = inMesh;
1100
1101     Mesh->BuildLinks(); //to do neighborhood searching
1102     polys = Mesh->GetPolys();
1103
1104     for (cellId=0, polys->InitTraversal(); polys->GetNextCell(npts,pts);
1105          cellId++)
1106     {
1107       if(DebugLevel>5) cout<<"->cellId="<<cellId<<endl;
1108       for (i=0; i < npts; i++)
1109       {
1110         if(DebugLevel>5) cout<<"-->i="<<i<<endl;
1111         p1 = pts[i];
1112         p2 = pts[(i+1)%npts];
1113
1114         if ( Verts[p1].edges == NULL )
1115         {
1116           Verts[p1].edges = vtkIdList::New();
1117           Verts[p1].edges->Allocate(16,6);
1118         }
1119         if ( Verts[p2].edges == NULL )
1120         {
1121           Verts[p2].edges = vtkIdList::New();
1122           Verts[p2].edges->Allocate(16,6);
1123         }
1124
1125         Mesh->GetCellEdgeNeighbors(cellId,p1,p2,neighbors);
1126         numNei = neighbors->GetNumberOfIds();
1127         if(DebugLevel>5) cout<<"-->numNei="<<numNei<<endl;
1128
1129         edge = VTK_SIMPLE_VERTEX;
1130         if ( numNei == 0 )
1131         {
1132           edge = VTK_BOUNDARY_EDGE_VERTEX;
1133         }
1134
1135         else if ( numNei >= 2 )
1136         {
1137             // check to make sure that this edge hasn't been marked already
1138           for (j=0; j < numNei; j++)
1139           {
1140             if ( neighbors->GetId(j) < cellId )
1141             {
1142               break;
1143             }
1144           }
1145           if ( j >= numNei )
1146           {
1147             edge = VTK_FEATURE_EDGE_VERTEX;
1148           }
1149         }
1150
1151         else if ( numNei == 1 && (nei=neighbors->GetId(0)) > cellId )
1152         {
1153           vtkPolygon::ComputeNormal(inPts,npts,pts,normal);
1154           Mesh->GetCellPoints(nei,numNeiPts,neiPts);
1155           vtkPolygon::ComputeNormal(inPts,numNeiPts,neiPts,neiNormal);
1156
1157           if ( this->FeatureEdgeSmoothing &&
1158                vtkMath::Dot(normal,neiNormal) <= CosFeatureAngle )
1159           {
1160             edge = VTK_FEATURE_EDGE_VERTEX;
1161           }
1162         }
1163         else // a visited edge; skip rest of analysis
1164         {
1165           continue;
1166         }
1167
1168         if ( edge && Verts[p1].type == VTK_SIMPLE_VERTEX )
1169         {
1170           Verts[p1].edges->Reset();
1171           Verts[p1].edges->InsertNextId(p2);
1172           Verts[p1].type = edge;
1173         }
1174         else if ( (edge && Verts[p1].type == VTK_BOUNDARY_EDGE_VERTEX) ||
1175                   (edge && Verts[p1].type == VTK_FEATURE_EDGE_VERTEX) ||
1176                   (!edge && Verts[p1].type == VTK_SIMPLE_VERTEX ) )
1177         {
1178           Verts[p1].edges->InsertNextId(p2);
1179           if ( Verts[p1].type && edge == VTK_BOUNDARY_EDGE_VERTEX )
1180           {
1181             Verts[p1].type = VTK_BOUNDARY_EDGE_VERTEX;
1182           }
1183         }
1184
1185         if ( edge && Verts[p2].type == VTK_SIMPLE_VERTEX )
1186         {
1187           Verts[p2].edges->Reset();
1188           Verts[p2].edges->InsertNextId(p1);
1189           Verts[p2].type = edge;
1190         }
1191         else if ( (edge && Verts[p2].type == VTK_BOUNDARY_EDGE_VERTEX ) ||
1192                   (edge && Verts[p2].type == VTK_FEATURE_EDGE_VERTEX) ||
1193                   (!edge && Verts[p2].type == VTK_SIMPLE_VERTEX ) )
1194         {
1195           Verts[p2].edges->InsertNextId(p1);
1196           if ( Verts[p2].type && edge == VTK_BOUNDARY_EDGE_VERTEX )
1197           {
1198             Verts[p2].type = VTK_BOUNDARY_EDGE_VERTEX;
1199           }
1200         }
1201       }
1202     }
1203
1204     inMesh->Delete();
1205
1206     neighbors->Delete();
1207   }//if strips or polys
1208
1209     //post-process edge vertices to make sure we can smooth them
1210   for (i=0; i<numPts; i++)
1211   {
1212     if ( Verts[i].type == VTK_SIMPLE_VERTEX )
1213     {
1214       numSimple++;
1215     }
1216
1217     else if ( Verts[i].type == VTK_FIXED_VERTEX )
1218     {
1219       numFixed++;
1220     }
1221
1222     else if ( Verts[i].type == VTK_FEATURE_EDGE_VERTEX ||
1223               Verts[i].type == VTK_BOUNDARY_EDGE_VERTEX )
1224     { //see how many edges; if two, what the angle is
1225
1226       if ( !this->BoundarySmoothing &&
1227            Verts[i].type == VTK_BOUNDARY_EDGE_VERTEX )
1228       {
1229         if(DebugLevel>5) cout<<"Verts[i].type = VTK_FIXED_VERTEX; 4"<<endl;
1230         Verts[i].type = VTK_FIXED_VERTEX;
1231         numBEdges++;
1232       }
1233
1234       else if ( (npts = Verts[i].edges->GetNumberOfIds()) != 2 )
1235       {
1236         if(DebugLevel>5) cout<<"Verts["<<i<<"].type = VTK_FIXED_VERTEX; 5"<<endl;
1237         Verts[i].type = VTK_FIXED_VERTEX;
1238         numFixed++;
1239       }
1240
1241       else //check angle between edges
1242       {
1243         inPts->GetPoint(Verts[i].edges->GetId(0),x1);
1244         inPts->GetPoint(i,x2);
1245         inPts->GetPoint(Verts[i].edges->GetId(1),x3);
1246
1247         for (k=0; k<3; k++)
1248         {
1249           l1[k] = x2[k] - x1[k];
1250           l2[k] = x3[k] - x2[k];
1251         }
1252         if ( vtkMath::Normalize(l1) >= 0.0 &&
1253              vtkMath::Normalize(l2) >= 0.0 &&
1254              vtkMath::Dot(l1,l2) < CosEdgeAngle)
1255         {
1256           if(DebugLevel>5) cout<<"Verts["<<i<<"].type = VTK_FIXED_VERTEX; 6"<<endl;
1257           Verts[i].type = VTK_FIXED_VERTEX;
1258           numFixed++;
1259         }
1260         else
1261         {
1262           if ( Verts[i].type == VTK_FEATURE_EDGE_VERTEX )
1263           {
1264             numFEdges++;
1265           }
1266           else
1267           {
1268             numBEdges++;
1269           }
1270         }
1271       }//if along edge
1272     }//if edge vertex
1273   }//for all points
1274
1275   if(DebugLevel>5) {
1276     cout<<"Found\n\t" << numSimple << " simple vertices\n\t"
1277       << numFEdges << " feature edge vertices\n\t"
1278       << numBEdges << " boundary edge vertices\n\t"
1279       << numFixed << " fixed vertices\n\t"<<endl;
1280     cout<<"1:numPts="<<numPts<<endl;
1281   }
1282
1283   for (i=0; i<numPts; i++)
1284   {
1285     if(DebugLevel>5) cout<<"Verts["<<i<<"].type="<<VertexType2Str(Verts[i].type)<<endl;
1286     if(Verts[i].edges != NULL && (npts = Verts[i].edges->GetNumberOfIds()) > 0)
1287     {
1288       for (j=0; j<npts; j++)
1289       {
1290         if(DebugLevel>5) cout<<"Verts["<<i<<"].edges->GetId("<<j<<")="<<Verts[i].edges->GetId(j)<<endl;
1291       }//for all connected points
1292     }
1293   }
1294
1295   //Copy node type info from Verts
1296   EG_VTKDCN(vtkCharArray, node_type, grid, "node_type");
1297   if(DebugLevel>5) cout<<"nodes.size()="<<nodes.size()<<endl;
1298   foreach(vtkIdType node,nodes)
1299   {
1300     if(DebugLevel>5) cout<<"Verts["<<node<<"].type="<<VertexType2Str(Verts[node].type)<<endl;
1301     node_type->SetValue(node,Verts[node].type);
1302   }
1303
1304   //free up connectivity storage
1305   for (int i=0; i<numPts; i++)
1306   {
1307     if ( Verts[i].edges != NULL )
1308     {
1309       Verts[i].edges->Delete();
1310       Verts[i].edges = NULL;
1311     }
1312   }
1313   delete [] Verts;
1314
1315   return(0);
1316 }
1317 //End of UpdateNodeType
1318
1319 // DEFINITIONS:
1320 // Normal cell: nothing has changed
1321 // Dead cell: the cell does not exist anymore
1322 // Mutated cell: the cell's form has changed
1323 // Mutilated cell: the cell has less points than before
1324
1325 vtkIdType Operation::FindSnapPoint(vtkUnstructuredGrid *src, vtkIdType DeadNode,QSet <vtkIdType> & DeadCells,QSet <vtkIdType> & MutatedCells,QSet <vtkIdType> & MutilatedCells, int& N_newpoints, int& N_newcells)
1326 {
1327   getAllSurfaceCells(cells,src);
1328   getNodesFromCells(cells, nodes, src);
1329   setGrid(src);
1330   setCells(cells);
1331
1332   UpdateNodeType();
1333
1334   EG_VTKDCN(vtkCharArray, node_type, src, "node_type");
1335   if(node_type->GetValue(DeadNode)==VTK_FIXED_VERTEX)
1336   {
1337     cout<<"Sorry, unable to remove fixed vertex."<<endl;
1338     return(-1);
1339   }
1340
1341     //src grid info
1342   int N_points=src->GetNumberOfPoints();
1343   int N_cells=src->GetNumberOfCells();
1344   N_newpoints=-1;
1345   N_newcells=0;
1346
1347   vtkIdType SnapPoint=-1;
1348
1349   foreach(vtkIdType PSP, n2n[DeadNode])
1350   {
1351     bool IsValidSnapPoint=true;
1352
1353     if(DebugLevel>10) cout<<"====>PSP="<<PSP<<endl;
1354     bool IsTetra=true;
1355     if(NumberOfCommonPoints(DeadNode,PSP,IsTetra)>2)//common point check
1356     {
1357       if(DebugLevel>10) cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
1358       IsValidSnapPoint=false;
1359     }
1360     if(IsTetra)//tetra check
1361     {
1362       if(DebugLevel>10) cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
1363       IsValidSnapPoint=false;
1364     }
1365
1366     //count number of points and cells to remove + analyse cell transformations
1367     N_newpoints=-1;
1368     N_newcells=0;
1369     DeadCells.clear();
1370     MutatedCells.clear();
1371     MutilatedCells.clear();
1372     foreach(vtkIdType C, n2c[DeadNode])//loop through potentially dead cells
1373     {
1374       //get points around cell
1375       vtkIdType N_pts, *pts;
1376       src->GetCellPoints(C, N_pts, pts);
1377
1378       bool ContainsSnapPoint=false;
1379       bool invincible=false;
1380       for(int i=0;i<N_pts;i++)
1381       {
1382         if(DebugLevel>10) cout<<"pts["<<i<<"]="<<pts[i]<<" and PSP="<<PSP<<endl;
1383         if(pts[i]==PSP) {ContainsSnapPoint=true;}
1384         if(pts[i]!=DeadNode && pts[i]!=PSP &&  n2c[pts[i]].size()<=1) invincible=true;
1385       }
1386       if(ContainsSnapPoint)
1387       {
1388         if(N_pts==3)//dead cell
1389         {
1390           if(invincible)//Check that empty lines aren't left behind when a cell is killed
1391           {
1392             if(DebugLevel>10) cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
1393             IsValidSnapPoint=false;
1394           }
1395           else
1396           {
1397             DeadCells.insert(C);
1398             N_newcells-=1;
1399             if(DebugLevel>10) cout<<"cell "<<C<<" has been pwned!"<<endl;
1400           }
1401         }
1402         else
1403         {
1404           cout<<"RED ALERT: Xenomorph detected!"<<endl;
1405           EG_BUG;
1406         }
1407       }
1408       else
1409       {
1410         vtkIdType src_N_pts, *src_pts;
1411         src->GetCellPoints(C, src_N_pts, src_pts);
1412
1413         if(src_N_pts!=3)
1414         {
1415           cout<<"RED ALERT: Xenomorph detected!"<<endl;
1416           EG_BUG;
1417         }
1418
1419         vtkIdType OldTriangle[3];
1420         vtkIdType NewTriangle[3];
1421
1422         for(int i=0;i<src_N_pts;i++)
1423         {
1424           OldTriangle[i]=src_pts[i];
1425           NewTriangle[i]=( (src_pts[i]==DeadNode) ? PSP : src_pts[i] );
1426         }
1427         vec3_t Old_N= triNormal(src, OldTriangle[0], OldTriangle[1], OldTriangle[2]);
1428         vec3_t New_N= triNormal(src, NewTriangle[0], NewTriangle[1], NewTriangle[2]);
1429         double OldArea=Old_N.abs();
1430         double NewArea=New_N.abs();
1431         double scal=Old_N*New_N;
1432         double cross=(Old_N.cross(New_N)).abs();//double-cross on Nar Shadaa B-)
1433
1434         if(DebugLevel>10) {
1435           cout<<"OldArea="<<OldArea<<endl;
1436           cout<<"NewArea="<<NewArea<<endl;
1437           cout<<"scal="<<scal<<endl;
1438           cout<<"cross="<<cross<<endl;
1439         }
1440
1441         if(Old_N*New_N<Old_N*Old_N*1./100.)//area + inversion check
1442         {
1443           if(DebugLevel>10) cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
1444           IsValidSnapPoint=false;
1445         }
1446
1447         //mutated cell
1448         MutatedCells.insert(C);
1449         if(DebugLevel>10) cout<<"cell "<<C<<" has been infected!"<<endl;
1450       }
1451     }
1452
1453     if(N_cells+N_newcells<=0)//survivor check
1454     {
1455       if(DebugLevel>10) cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
1456       IsValidSnapPoint=false;
1457     }
1458
1459     if(node_type->GetValue(DeadNode)==VTK_BOUNDARY_EDGE_VERTEX && node_type->GetValue(PSP)==VTK_SIMPLE_VERTEX)
1460     {
1461       if(DebugLevel>10) cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
1462       IsValidSnapPoint=false;
1463     }
1464
1465     if(node_type->GetValue(DeadNode)==VTK_BOUNDARY_EDGE_VERTEX)
1466     {
1467       int BC=0;
1468       QVector <vtkIdType> Peons;
1469       getNeighbours(DeadNode, Peons, BC);
1470       if(!Peons.contains(PSP))
1471       {
1472         if(DebugLevel>0) cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
1473         IsValidSnapPoint=false;
1474       }
1475     }
1476
1477     if(node_type->GetValue(DeadNode)==VTK_FEATURE_EDGE_VERTEX && node_type->GetValue(PSP)==VTK_SIMPLE_VERTEX)
1478     {
1479       if(DebugLevel>10) cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
1480       IsValidSnapPoint=false;
1481     }
1482
1483     if(node_type->GetValue(DeadNode)==VTK_FEATURE_EDGE_VERTEX)
1484     {
1485       int BC=0;
1486       QVector <vtkIdType> Peons;
1487       getNeighbours(DeadNode, Peons, BC);
1488       if(!Peons.contains(PSP))
1489       {
1490         if(DebugLevel>0) cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
1491         IsValidSnapPoint=false;
1492       }
1493     }
1494
1495     if(IsValidSnapPoint) {SnapPoint=PSP; break;}
1496   }//end of loop through potential SnapPoints
1497
1498   if(DebugLevel>10)
1499   {
1500     cout<<"AT FINDSNAPPOINT EXIT"<<endl;
1501     cout<<"N_points="<<N_points<<endl;
1502     cout<<"N_cells="<<N_cells<<endl;
1503     cout<<"N_newpoints="<<N_newpoints<<endl;
1504     cout<<"N_newcells="<<N_newcells<<endl;
1505   }
1506   return(SnapPoint);
1507 }
1508 //End of FindSnapPoint
1509
1510 bool Operation::DeletePoint_2(vtkUnstructuredGrid *src, vtkIdType DeadNode, int& N_newpoints, int& N_newcells)
1511 {
1512   getAllSurfaceCells(cells,src);
1513 //   getNodesFromCells(cells, nodes, src);
1514   setGrid(src);
1515   setCells(cells);
1516   UpdateNodeType();
1517
1518   //src grid info
1519   int N_points=src->GetNumberOfPoints();
1520   int N_cells=src->GetNumberOfCells();
1521   N_newpoints=-1;
1522   N_newcells=0;
1523
1524   if(DeadNode<0 || DeadNode>=N_points)
1525   {
1526     cout<<"Warning: Point out of range: DeadNode="<<DeadNode<<" N_points="<<N_points<<endl;
1527     return(false);
1528   }
1529
1530   QSet <vtkIdType> DeadCells;
1531   QSet <vtkIdType> MutatedCells;
1532   QSet <vtkIdType> MutilatedCells;
1533
1534   if(DebugLevel>10) {
1535     cout<<"BEFORE FINDSNAPPOINT"<<endl;
1536     cout<<"N_points="<<N_points<<endl;
1537     cout<<"N_cells="<<N_cells<<endl;
1538     cout<<"N_newpoints="<<N_newpoints<<endl;
1539     cout<<"N_newcells="<<N_newcells<<endl;
1540   }
1541   vtkIdType SnapPoint=FindSnapPoint(src,DeadNode,DeadCells,MutatedCells,MutilatedCells, N_newpoints, N_newcells);
1542
1543   if(DebugLevel>0) cout<<"===>DeadNode="<<DeadNode<<" moving to SNAPPOINT="<<SnapPoint<<" DebugLevel="<<DebugLevel<<endl;
1544   if(SnapPoint<0) {cout<<"Sorry no possible SnapPoint found."<<endl; return(false);}
1545
1546   //allocate
1547   if(DebugLevel>10) {
1548     cout<<"BEFORE ALLOCATION"<<endl;
1549     cout<<"N_points="<<N_points<<endl;
1550     cout<<"N_cells="<<N_cells<<endl;
1551     cout<<"N_newpoints="<<N_newpoints<<endl;
1552     cout<<"N_newcells="<<N_newcells<<endl;
1553   }
1554   N_points=src->GetNumberOfPoints();
1555   N_cells=src->GetNumberOfCells();
1556
1557   if(DebugLevel>47) {
1558     cout<<"N_points="<<N_points<<endl;
1559     cout<<"N_cells="<<N_cells<<endl;
1560     cout<<"N_newpoints="<<N_newpoints<<endl;
1561     cout<<"N_newcells="<<N_newcells<<endl;
1562   }
1563   EG_VTKSP(vtkUnstructuredGrid,dst);
1564   allocateGrid(dst,N_cells+N_newcells,N_points+N_newpoints);
1565
1566   //vector used to redefine the new point IDs
1567   QVector <vtkIdType> OffSet(N_points);
1568
1569   //copy undead points
1570   vtkIdType dst_id_node=0;
1571   for (vtkIdType src_id_node = 0; src_id_node < N_points; src_id_node++) {//loop through src points
1572     if(src_id_node!=DeadNode)//if the node isn't dead, copy it
1573     {
1574       vec3_t x;
1575       src->GetPoints()->GetPoint(src_id_node, x.data());
1576       dst->GetPoints()->SetPoint(dst_id_node, x.data());
1577       copyNodeData(src, src_id_node, dst, dst_id_node);
1578       OffSet[src_id_node]=src_id_node-dst_id_node;
1579       dst_id_node++;
1580     }
1581     else
1582     {
1583       if(DebugLevel>0) cout<<"src_id_node="<<src_id_node<<" dst_id_node="<<dst_id_node<<endl;
1584     }
1585   };
1586   if(DebugLevel>10) {
1587     cout<<"DeadCells="<<DeadCells<<endl;
1588     cout<<"MutatedCells="<<MutatedCells<<endl;
1589     cout<<"MutilatedCells="<<MutilatedCells<<endl;
1590   }
1591   //Copy undead cells
1592   for (vtkIdType id_cell = 0; id_cell < src->GetNumberOfCells(); ++id_cell) {//loop through src cells
1593     if(!DeadCells.contains(id_cell))//if the cell isn't dead
1594     {
1595       vtkIdType src_N_pts, *src_pts;
1596       vtkIdType dst_N_pts, *dst_pts;
1597       src->GetCellPoints(id_cell, src_N_pts, src_pts);
1598
1599       vtkIdType type_cell = src->GetCellType(id_cell);
1600       if(DebugLevel>10) cout<<"-->id_cell="<<id_cell<<endl;
1601       if(DebugLevel>10) for(int i=0;i<src_N_pts;i++) cout<<"src_pts["<<i<<"]="<<src_pts[i]<<endl;
1602 //       src->GetCellPoints(id_cell, dst_N_pts, dst_pts);
1603       dst_N_pts=src_N_pts;
1604       dst_pts=new vtkIdType[dst_N_pts];
1605       if(MutatedCells.contains(id_cell))//mutated cell
1606       {
1607         if(DebugLevel>10) cout<<"processing mutated cell "<<id_cell<<endl;
1608         for(int i=0;i<src_N_pts;i++)
1609         {
1610           if(src_pts[i]==DeadNode) {
1611             if(DebugLevel>10) {
1612               cout<<"SnapPoint="<<SnapPoint<<endl;
1613               cout<<"OffSet[SnapPoint]="<<OffSet[SnapPoint]<<endl;
1614               cout<<"src_pts["<<i<<"]="<<src_pts[i]<<endl;
1615             }
1616             dst_pts[i]=SnapPoint-OffSet[SnapPoint];
1617           }
1618           else dst_pts[i]=src_pts[i]-OffSet[src_pts[i]];
1619         }
1620         if(DebugLevel>10) cout<<"--->dst_pts:"<<endl;
1621         if(DebugLevel>10) for(int i=0;i<dst_N_pts;i++) cout<<"dst_pts["<<i<<"]="<<dst_pts[i]<<endl;
1622
1623       }
1624       else if(MutilatedCells.contains(id_cell))//mutilated cell
1625       {
1626         if(DebugLevel>10) cout<<"processing mutilated cell "<<id_cell<<endl;
1627
1628         if(type_cell==VTK_QUAD) {
1629           type_cell=VTK_TRIANGLE;
1630           dst_N_pts-=1;
1631         }
1632         else {cout<<"FATAL ERROR: Unknown mutilated cell detected! It is not a quad! Potential xenomorph infestation!"<<endl;EG_BUG;}
1633         //merge points
1634         int j=0;
1635         for(int i=0;i<src_N_pts;i++)
1636         {
1637           if(src_pts[i]==SnapPoint) { dst_pts[j]=SnapPoint-OffSet[SnapPoint];j++; }//SnapPoint
1638           else if(src_pts[i]!=DeadNode) { dst_pts[j]=src_pts[i]-OffSet[src_pts[i]];j++; }//pre-snap/dead + post-snap/dead
1639           //do nothing in case of DeadNode
1640         }
1641       }
1642       else//normal cell
1643       {
1644         if(DebugLevel>10) cout<<"processing normal cell "<<id_cell<<endl;
1645         for(int i=0;i<src_N_pts;i++)
1646         {
1647           dst_pts[i]=src_pts[i]-OffSet[src_pts[i]];
1648         }
1649       }
1650       //copy the cell
1651       vtkIdType id_new_cell = dst->InsertNextCell(type_cell, dst_N_pts, dst_pts);
1652       copyCellData(src, id_cell, dst, id_new_cell);
1653       if(DebugLevel>10) {
1654         cout<<"===Copying cell "<<id_cell<<" to "<<id_new_cell<<"==="<<endl;
1655         cout<<"src_pts:"<<endl;
1656         for(int i=0;i<src_N_pts;i++) cout<<"src_pts["<<i<<"]="<<src_pts[i]<<endl;
1657         cout<<"dst_pts:"<<endl;
1658         for(int i=0;i<dst_N_pts;i++) cout<<"dst_pts["<<i<<"]="<<dst_pts[i]<<endl;
1659         cout<<"OffSet="<<OffSet<<endl;
1660         cout<<"===Copying cell end==="<<endl;
1661       }
1662       delete dst_pts;
1663     }
1664   };
1665 //   cout_grid(cout,dst,true,true,true,true);
1666   makeCopy(dst, src);
1667   return(true);
1668 }
1669 //End of DeletePoint_2
1670
1671 void Operation::TxtSave(QString a_filename)
1672 {
1673   cout << a_filename.toAscii().data() << endl;
1674   ofstream file;
1675   file.open(a_filename.toAscii().data());
1676   cout_grid(file,grid,true,true,true,true);
1677   file.close();
1678 }
1679
1680 void Operation::DualSave(QString a_filename)
1681 {
1682   TxtSave(a_filename+".txt");
1683   GuiMainWindow::pointer()->QuickSave(a_filename+".vtu");
1684 }