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