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
  26 #include "geometrytools.h"
  27 using namespace GeometryTools;
  28
  29 #include <QApplication>
  30
  31 QSet<Operation*> Operation::garbage_operations;
  32
  33 void Operation::collectGarbage()
  34 {
  35   QSet<Operation*> delete_operations;
  36   foreach (Operation *op, garbage_operations) {
  37     if (!op->getThread().isRunning()) {
  38       delete_operations.insert(op);
  39       cout << "deleting Operation " << op << endl;
  40       delete op;
  41     };
  42   };
  43   foreach (Operation *op, delete_operations) {
  44     garbage_operations.remove(op);
  45   };
  46 };
  47
  48 Operation::Operation()
  49 {
  50   grid = NULL;
  51   gui = false;
  52   err = NULL;
  53   autoset = true;
  54 };
  55
  56 Operation::~Operation()
  57 {
  58   if (err) {
  59     err->display();
  60     delete err;
  61   };
  62 };
  63
  64 void Operation::del()
  65 {
  66   garbage_operations.insert(this);
  67 };
  68
  69 void OperationThread::run()
  70 {
  71   try {
  72     GuiMainWindow::lock();
  73     GuiMainWindow::pointer()->setBusy();
  74     op->operate();
  75   } catch (Error err) {
  76     op->err = new Error();
  77     *(op->err) = err;
  78   };
  79   GuiMainWindow::unlock();
  80   GuiMainWindow::pointer()->setIdle();
  81 };
  82
  83 void Operation::operator()()
  84 {
  85   if (gui) {
  86     if (GuiMainWindow::tryLock()) {
  87       checkGrid();
  88       thread.setOperation(this);
  89       GuiMainWindow::unlock();
  90       thread.start(QThread::LowPriority);
  91     } else {
  92       QMessageBox::warning(NULL, "not permitted", "Operation is not permitted while background process is running!");
  93     };
  94   } else {
  95     checkGrid();
  96     operate();
  97   };
  98   /*
  99   checkGrid();
 100   operate();
 101   */
 102 };
 103
 104 void Operation::setAllCells()
 105 {
 106   QVector<vtkIdType> all_cells;
 107   getAllCells(all_cells, grid);
 108   setCells(all_cells);
 109 };
 110
 111 void Operation::setAllVolumeCells()
 112 {
 113   QVector<vtkIdType> cells;
 114   getAllVolumeCells(cells, grid);
 115   setCells(cells);
 116 };
 117
 118 void Operation::setAllSurfaceCells()
 119 {
 120   QVector<vtkIdType> cells;
 121   getAllSurfaceCells(cells, grid);
 122   setCells(cells);
 123 };
 124
 125 void Operation::initMapping()
 126 {
 127   nodes_map.resize(nodes.size());
 128   for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
 129     nodes_map[i_nodes] = nodes[i_nodes];
 130   };
 131   cells_map.resize(cells.size());
 132   for (int i_cells = 0; i_cells < cells.size(); ++i_cells) {
 133     cells_map[i_cells] = cells[i_cells];
 134   };
 135 };
 136
 137 void Operation::checkGrid()
 138 {
 139   if (grid == NULL) {
 140     grid = GuiMainWindow::pointer()->getGrid();
 141   };
 142   if ((cells.size() == 0) && autoset) {
 143     setAllCells();
 144   };
 145 };
 146
 147 void Operation::updateActors()
 148 {
 149   mainWindow()->updateActors();
 150 };
 151
 152 GuiMainWindow* Operation::mainWindow()
 153 {
 154   return GuiMainWindow::pointer();
 155 };
 156
 157 void Operation::populateBoundaryCodes(QListWidget *lw)
 158 {
 159   QSet<int> bcs;
 160   mainWindow()->getAllBoundaryCodes(bcs);
 161   foreach(int bc, bcs) {
 162     QListWidgetItem *lwi = new QListWidgetItem(lw);
 163     lwi->setCheckState(Qt::Unchecked);
 164     QString text = "";
 165     QTextStream ts(&text);
 166     ts << bc;
 167     lwi->setText(text);
 168     lwi->setFlags(Qt::ItemIsUserCheckable | Qt::ItemIsEnabled);
 169   };
 170 };
 171
 172 stencil_t Operation::getStencil(vtkIdType id_cell1, int j1)
 173 {
 174   stencil_t S;
 175   S.valid = true;
 176   S.id_cell1 = id_cell1;
 177   if (c2c[_cells[id_cell1]][j1] != -1) {
 178     S.id_cell2 = cells[c2c[_cells[id_cell1]][j1]];
 179     if (grid->GetCellType(S.id_cell2) != VTK_TRIANGLE) {
 180       EG_BUG;
 181     };
 182     vtkIdType N1, N2, *pts1, *pts2;
 183     grid->GetCellPoints(S.id_cell1, N1, pts1);
 184     grid->GetCellPoints(S.id_cell2, N2, pts2);
 185     if      (j1 == 0) { S.p[0] = pts1[2]; S.p[1] = pts1[0]; S.p[3] = pts1[1]; }
 186     else if (j1 == 1) { S.p[0] = pts1[0]; S.p[1] = pts1[1]; S.p[3] = pts1[2]; }
 187     else if (j1 == 2) { S.p[0] = pts1[1]; S.p[1] = pts1[2]; S.p[3] = pts1[0]; };
 188     bool p2 = false;
 189     if (c2c[_cells[S.id_cell2]][0] != -1) {
 190       if (cells[c2c[_cells[S.id_cell2]][0]] == S.id_cell1) {
 191         S.p[2] = pts2[2];
 192         p2 = true;
 193       };
 194     };
 195     if (c2c[_cells[S.id_cell2]][1] != -1) {
 196       if (cells[c2c[_cells[S.id_cell2]][1]] == S.id_cell1) {
 197         S.p[2] = pts2[0];
 198         p2 = true;
 199       };
 200     };
 201     if (c2c[_cells[S.id_cell2]][2] != -1) {
 202       if (cells[c2c[_cells[S.id_cell2]][2]] == S.id_cell1) {
 203         S.p[2] = pts2[1];
 204         p2 = true;
 205       };
 206     };
 207     if (!p2) {
 208       EG_BUG;
 209     };
 210   } else {
 211     S.valid = false;
 212     S.id_cell2 = -1;
 213     vtkIdType N1, *pts1;
 214     grid->GetCellPoints(S.id_cell1, N1, pts1);
 215     if      (j1 == 0) { S.p[0] = pts1[2]; S.p[1] = pts1[0]; S.p[3] = pts1[1]; }
 216     else if (j1 == 1) { S.p[0] = pts1[0]; S.p[1] = pts1[1]; S.p[3] = pts1[2]; }
 217     else if (j1 == 2) { S.p[0] = pts1[1]; S.p[1] = pts1[2]; S.p[3] = pts1[0]; };
 218   };
 219   return S;
 220 };
 221
 222 ostream& operator<<(ostream &out, stencil_t S)
 223 {
 224   out<<"S.id_cell1="<<S.id_cell1<<" ";
 225   out<<"S.id_cell2="<<S.id_cell2<<" ";
 226   out<<"S.valid="<<S.valid<<" ";
 227   out<<"[";
 228   for(int i=0;i<4;i++){
 229     out<<S.p[i];
 230     if(i!=3) out<<",";
 231   }
 232   out<<"]";
 233   return(out);
 234 }
 235
 236 //////////////////////////////////////////////
 237 double CurrentVertexAvgDist(vtkIdType a_vertex,QVector< QSet< int > > &n2n,vtkUnstructuredGrid *a_grid)
 238 {
 239   double total_dist=0;
 240   double avg_dist=0;
 241   int N=n2n[a_vertex].size();
 242   vec3_t C;
 243   a_grid->GetPoint(a_vertex, C.data());
 244   foreach(int i,n2n[a_vertex])
 245   {
 246     vec3_t M;
 247     a_grid->GetPoint(i, M.data());
 248     total_dist+=(M-C).abs();
 249   }
 250   avg_dist=total_dist/(double)N;
 251   return(avg_dist);
 252 }
 253
 254 double CurrentMeshDensity(vtkIdType a_vertex,QVector< QSet< int > > &n2n,vtkUnstructuredGrid *a_grid)
 255 {
 256   double total_dist=0;
 257   double avg_dist=0;
 258   int N=n2n[a_vertex].size();
 259   vec3_t C;
 260   a_grid->GetPoint(a_vertex, C.data());
 261   foreach(int i,n2n[a_vertex])
 262   {
 263     vec3_t M;
 264     a_grid->GetPoint(i, M.data());
 265     total_dist+=(M-C).abs();
 266   }
 267   avg_dist=total_dist/(double)N;
 268   double avg_density=1./avg_dist;
 269   return(avg_density);
 270 }
 271
 272 double DesiredVertexAvgDist(vtkIdType a_vertex,QVector< QSet< int > > &n2n,vtkUnstructuredGrid *a_grid)
 273 {
 274   double total_dist=0;
 275   double avg_dist=0;
 276   EG_VTKDCN(vtkDoubleArray, node_meshdensity, a_grid, "node_meshdensity");
 277   int N=n2n[a_vertex].size();
 278   foreach(int i,n2n[a_vertex])
 279   {
 280     total_dist+=1./node_meshdensity->GetValue(i);
 281   }
 282   avg_dist=total_dist/(double)N;
 283   return(avg_dist);
 284 }
 285
 286 double DesiredMeshDensity(vtkIdType a_vertex,QVector< QSet< int > > &n2n,vtkUnstructuredGrid *a_grid)
 287 {
 288   double total_density=0;
 289   double avg_density=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_density+=node_meshdensity->GetValue(i);
 295   }
 296   avg_density=total_density/(double)N;
 297   return(avg_density);
 298 }
 299
 300 ///////////////////////////////////////////
 301 vtkIdType Operation::getClosestNode(vtkIdType a_id_node,vtkUnstructuredGrid* a_grid)
 302 {
 303   vec3_t C;
 304   a_grid->GetPoint(a_id_node,C.data());
 305   vtkIdType id_minlen=-1;
 306   double minlen=-1;
 307   foreach(vtkIdType neighbour,n2n[a_id_node])
 308   {
 309     vec3_t M;
 310     a_grid->GetPoint(neighbour,M.data());
 311     double len=(M-C).abs();
 312     if(minlen<0 or len<minlen)
 313     {
 314       minlen=len;
 315       id_minlen=neighbour;
 316     }
 317   }
 318   return(id_minlen);
 319 }
 320
 321 vtkIdType Operation::getFarthestNode(vtkIdType a_id_node,vtkUnstructuredGrid* a_grid)
 322 {
 323   vec3_t C;
 324   a_grid->GetPoint(a_id_node,C.data());
 325   vtkIdType id_maxlen=-1;
 326   double maxlen=-1;
 327   foreach(vtkIdType neighbour,n2n[a_id_node])
 328   {
 329     vec3_t M;
 330     a_grid->GetPoint(neighbour,M.data());
 331     double len=(M-C).abs();
 332     if(maxlen<0 or len>maxlen)
 333     {
 334       maxlen=len;
 335       id_maxlen=neighbour;
 336     }
 337   }
 338   return(id_maxlen);
 339 }
 340
 341 bool Operation::SwapCells(vtkUnstructuredGrid* a_grid, stencil_t S)
 342 {
 343   bool swap = false;
 344   if (S.valid) {
 345     vec3_t x3[4], x3_0(0,0,0);
 346     vec2_t x[4];
 347     for (int k = 0; k < 4; ++k) {
 348       a_grid->GetPoints()->GetPoint(S.p[k], x3[k].data());
 349       x3_0 += x3[k];
 350     };
 351     vec3_t n1 = triNormal(x3[0], x3[1], x3[3]);
 352     vec3_t n2 = triNormal(x3[1], x3[2], x3[3]);
 353     n1.normalise();
 354     n2.normalise();
 355     if ( (n1*n2) > 0.8) {
 356       vec3_t n = n1 + n2;
 357       n.normalise();
 358       vec3_t ex = orthogonalVector(n);
 359       vec3_t ey = ex.cross(n);
 360       for (int k = 0; k < 4; ++k) {
 361         x[k] = vec2_t(x3[k]*ex, x3[k]*ey);
 362       };
 363       vec2_t r1, r2, r3, u1, u2, u3;
 364       r1 = 0.5*(x[0] + x[1]); u1 = turnLeft(x[1] - x[0]);
 365       r2 = 0.5*(x[1] + x[2]); u2 = turnLeft(x[2] - x[1]);
 366       r3 = 0.5*(x[1] + x[3]); u3 = turnLeft(x[3] - x[1]);
 367       double k, l;
 368       vec2_t xm1, xm2;
 369       bool ok = true;
 370       if (intersection(k, l, r1, u1, r3, u3)) {
 371         xm1 = r1 + k*u1;
 372         if (intersection(k, l, r2, u2, r3, u3)) {
 373           xm2 = r2 + k*u2;
 374         } else {
 375           ok = false;
 376         };
 377       } else {
 378         ok = false;
 379         swap = true;
 380       };
 381       if (ok) {
 382         if ((xm1 - x[2]).abs() < (xm1 - x[0]).abs()) {
 383           swap = true;
 384         };
 385         if ((xm2 - x[0]).abs() < (xm2 - x[2]).abs()) {
 386           swap = true;
 387         };
 388       };
 389     };
 390   };
 391   if (swap) {
 392     vtkIdType new_pts1[3], new_pts2[3];
 393     new_pts1[0] = S.p[1];
 394     new_pts1[1] = S.p[2];
 395     new_pts1[2] = S.p[0];
 396     new_pts2[0] = S.p[2];
 397     new_pts2[1] = S.p[3];
 398     new_pts2[2] = S.p[0];
 399     a_grid->ReplaceCell(S.id_cell1, 3, new_pts1);
 400     a_grid->ReplaceCell(S.id_cell2, 3, new_pts2);
 401   };
 402   return(swap);
 403 }
 404
 405 void Operation::quad2triangle(vtkUnstructuredGrid* src,vtkIdType quadcell)
 406 {
 407   vtkIdType type_cell = grid->GetCellType(quadcell);
 408   cout<<"It's a "<<type_cell<<endl;
 409   if(type_cell==VTK_QUAD)
 410   {
 411     cout_grid(cout,src,true,true,true,true);
 412     EG_VTKSP(vtkUnstructuredGrid, dst);
 413     //src grid info
 414     int N_points=src->GetNumberOfPoints();
 415     int N_cells=src->GetNumberOfCells();
 416     allocateGrid(dst,N_cells+1,N_points);
 417
 418     for (vtkIdType id_node = 0; id_node < src->GetNumberOfPoints(); ++id_node) {
 419       vec3_t x;
 420       src->GetPoints()->GetPoint(id_node, x.data());
 421       dst->GetPoints()->SetPoint(id_node, x.data());
 422       copyNodeData(src, id_node, dst, id_node);
 423     };
 424     for (vtkIdType id_cell = 0; id_cell < src->GetNumberOfCells(); ++id_cell) {
 425       vtkIdType N_pts, *pts;
 426       vtkIdType type_cell = src->GetCellType(id_cell);
 427       src->GetCellPoints(id_cell, N_pts, pts);
 428       vtkIdType id_new_cell;
 429       vtkIdType id_new_cell1;
 430       vtkIdType id_new_cell2;
 431       if(id_cell!=quadcell)
 432       {
 433         id_new_cell = dst->InsertNextCell(type_cell, N_pts, pts);
 434         copyCellData(src, id_cell, dst, id_new_cell);
 435       }
 436       else
 437       {
 438         vtkIdType triangle1[3];
 439         vtkIdType triangle2[3];
 440         triangle1[0]=pts[1];
 441         triangle1[1]=pts[3];
 442         triangle1[2]=pts[0];
 443         triangle2[0]=pts[3];
 444         triangle2[1]=pts[1];
 445         triangle2[2]=pts[2];
 446         id_new_cell1 = dst->InsertNextCell(VTK_TRIANGLE, 3, triangle1);
 447         copyCellData(src, id_cell, dst, id_new_cell1);
 448         id_new_cell2 = dst->InsertNextCell(VTK_TRIANGLE, 3, triangle2);
 449         copyCellData(src, id_cell, dst, id_new_cell2);
 450         stencil_t S;
 451         S.id_cell1=id_new_cell1;
 452         S.id_cell2=id_new_cell2;
 453         S.p[0]=pts[0];
 454         S.p[1]=pts[1];
 455         S.p[2]=pts[2];
 456         S.p[3]=pts[3];
 457         S.valid=true;
 458         SwapCells(dst,S);
 459       }
 460     };
 461     cout_grid(cout,dst,true,true,true,true);
 462     makeCopy(dst, src);
 463   }//end of if quad
 464 }
 465
 466 void Operation::quad2triangle(vtkUnstructuredGrid* src,vtkIdType quadcell,vtkIdType MovingPoint)
 467 {
 468   vtkIdType type_cell = grid->GetCellType(quadcell);
 469   cout<<"It's a "<<type_cell<<endl;
 470   if(type_cell==VTK_QUAD)
 471   {
 472     cout_grid(cout,src,true,true,true,true);
 473     EG_VTKSP(vtkUnstructuredGrid, dst);
 474     //src grid info
 475     int N_points=src->GetNumberOfPoints();
 476     int N_cells=src->GetNumberOfCells();
 477     allocateGrid(dst,N_cells+1,N_points);
 478
 479     for (vtkIdType id_node = 0; id_node < src->GetNumberOfPoints(); ++id_node) {
 480       vec3_t x;
 481       src->GetPoints()->GetPoint(id_node, x.data());
 482       dst->GetPoints()->SetPoint(id_node, x.data());
 483       copyNodeData(src, id_node, dst, id_node);
 484     };
 485     for (vtkIdType id_cell = 0; id_cell < src->GetNumberOfCells(); ++id_cell) {
 486       vtkIdType N_pts, *pts;
 487       src->GetCellPoints(id_cell, N_pts, pts);
 488       vtkIdType type_cell = src->GetCellType(id_cell);
 489       vtkIdType id_new_cell;
 490       vtkIdType id_new_cell1;
 491       vtkIdType id_new_cell2;
 492       if(id_cell!=quadcell)
 493       {
 494         id_new_cell = dst->InsertNextCell(type_cell, N_pts, pts);
 495         copyCellData(src, id_cell, dst, id_new_cell);
 496       }
 497       else
 498       {
 499         vtkIdType triangle1[3];
 500         vtkIdType triangle2[3];
 501         if(MovingPoint==pts[0] || MovingPoint==pts[2])
 502         {
 503           triangle1[0]=pts[1];
 504           triangle1[1]=pts[3];
 505           triangle1[2]=pts[0];
 506           triangle2[0]=pts[3];
 507           triangle2[1]=pts[1];
 508           triangle2[2]=pts[2];
 509         }
 510         else
 511         {
 512           triangle1[0]=pts[2];
 513           triangle1[1]=pts[0];
 514           triangle1[2]=pts[1];
 515           triangle2[0]=pts[0];
 516           triangle2[1]=pts[2];
 517           triangle2[2]=pts[3];
 518         }
 519         id_new_cell1 = dst->InsertNextCell(VTK_TRIANGLE, 3, triangle1);
 520         copyCellData(src, id_cell, dst, id_new_cell1);
 521         id_new_cell2 = dst->InsertNextCell(VTK_TRIANGLE, 3, triangle2);
 522         copyCellData(src, id_cell, dst, id_new_cell2);
 523       }
 524     };
 525     cout_grid(cout,dst,true,true,true,true);
 526     makeCopy(dst, src);
 527   }//end of if quad
 528 }
 529
 530 int Operation::NumberOfCommonPoints(vtkIdType node1, vtkIdType node2, bool& IsTetra)
 531 {
 532 //   QVector< QSet< int > >     n2n
 533   QSet <int> node1_neighbours=n2n[node1];
 534   QSet <int> node2_neighbours=n2n[node2];
 535   QSet <int> intersection=node1_neighbours.intersect(node2_neighbours);
 536   int N=intersection.size();
 537   IsTetra=false;
 538   if(N==2)
 539   {
 540     QSet<int>::const_iterator p1=intersection.begin();
 541     QSet<int>::const_iterator p2=p1+1;
 542     vtkIdType intersection1=_nodes[*p1];
 543     vtkIdType intersection2=_nodes[*p2];
 544     if(n2n[intersection1].contains(intersection2))//if there's an edge between intersection1 and intersection2
 545     {
 546       //check if (node1,intersection1,intersection2) and (node2,intersection1,intersection2) are defined as cells!
 547   //     QVector< QSet< int > >         n2c
 548       QSet< int > S1=n2c[intersection1];
 549       QSet< int > S2=n2c[intersection2];
 550       QSet< int > Si=S1.intersect(S2);
 551       int counter=0;
 552       foreach(vtkIdType C,Si){
 553         vtkIdType N_pts, *pts;
 554         grid->GetCellPoints(C, N_pts, pts);
 555         for(int i=0;i<N_pts;i++)
 556         {
 557           if(pts[i]==node1 || pts[i]==node2) counter++;
 558         }
 559       }
 560       if(counter>=2) IsTetra=true;
 561     }
 562   }
 563   return(N);
 564 }
 565
 566 // vtkIdType Operation::FindSnapPoint(vtkUnstructuredGrid *src, vtkIdType DeadNode)
 567 // {
 568 //   return(0);
 569 // }
 570
 571 bool Operation::DeletePoint(vtkUnstructuredGrid *src, vtkIdType DeadNode)
 572 {
 573   EG_VTKSP(vtkUnstructuredGrid, dst);
 574
 575     //src grid info
 576   int N_points=src->GetNumberOfPoints();
 577   int N_cells=src->GetNumberOfCells();
 578   int N_newpoints=-1;
 579   int N_newcells=0;
 580
 581   if(DeadNode<0 || DeadNode>=N_points)
 582   {
 583     cout<<"Warning: Point out of range: DeadNode="<<DeadNode<<" N_points="<<N_points<<endl;
 584     return(false);
 585   }
 586
 587   QSet <vtkIdType> DeadCells;
 588   QSet <vtkIdType> MutatedCells;
 589   QSet <vtkIdType> MutilatedCells;
 590
 591   vtkIdType SnapPoint=-1;
 592   //Find closest point to DeadNode
 593 //   vtkIdType SnapPoint = getClosestNode(DeadNode,src);//DeadNode moves to SnapPoint
 594
 595   foreach(vtkIdType PSP, n2n[DeadNode])
 596   {
 597     bool IsValidSnapPoint=true;
 598
 599     cout<<"====>PSP="<<PSP<<endl;
 600     bool IsTetra=true;
 601     if(NumberOfCommonPoints(DeadNode,PSP,IsTetra)>2)//common point check
 602     {
 603       cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
 604       IsValidSnapPoint=false;
 605     }
 606     if(IsTetra)//tetra check
 607     {
 608       cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
 609       IsValidSnapPoint=false;
 610     }
 611
 612     //count number of points and cells to remove + analyse cell transformations
 613     N_newpoints=-1;
 614     N_newcells=0;
 615     DeadCells.clear();
 616     MutatedCells.clear();
 617     MutilatedCells.clear();
 618     foreach(vtkIdType C, n2c[DeadNode])//loop through potentially dead cells
 619     {
 620       cout<<"C="<<C<<endl;
 621       //get points around cell
 622       vtkIdType N_pts, *pts;
 623       src->GetCellPoints(C, N_pts, pts);
 624
 625       bool ContainsSnapPoint=false;
 626       bool invincible=false;
 627       for(int i=0;i<N_pts;i++)
 628       {
 629         cout<<"pts["<<i<<"]="<<pts[i]<<" and PSP="<<PSP<<endl;
 630         if(pts[i]==PSP) {ContainsSnapPoint=true;}
 631         if(pts[i]!=DeadNode && pts[i]!=PSP &&  n2c[pts[i]].size()<=1) invincible=true;
 632       }
 633       if(ContainsSnapPoint)
 634       {
 635         if(N_pts==3)//dead cell
 636         {
 637           if(invincible)//Check that empty lines aren't left behind when a cell is killed
 638           {
 639             cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
 640             IsValidSnapPoint=false;
 641           }
 642           else
 643           {
 644             DeadCells.insert(C);
 645             N_newcells-=1;
 646             cout<<"cell "<<C<<" has been pwned!"<<endl;
 647           }
 648         }
 649   /*      else if(N_pts==4)//mutilated cell
 650         {
 651           MutilatedCells.insert(C);
 652           cout<<"cell "<<C<<" has lost a limb!"<<endl;
 653         }*/
 654         else
 655         {
 656           cout<<"RED ALERT: Xenomorph detected!"<<endl;
 657           EG_BUG;
 658         }
 659       }
 660       else
 661       {
 662         vtkIdType src_N_pts, *src_pts;
 663         src->GetCellPoints(C, src_N_pts, src_pts);
 664
 665         if(src_N_pts!=3)
 666         {
 667           cout<<"RED ALERT: Xenomorph detected!"<<endl;
 668           EG_BUG;
 669         }
 670
 671         vtkIdType OldTriangle[3];
 672         vtkIdType NewTriangle[3];
 673
 674         for(int i=0;i<src_N_pts;i++)
 675         {
 676           OldTriangle[i]=src_pts[i];
 677           NewTriangle[i]=( (src_pts[i]==DeadNode) ? PSP : src_pts[i] );
 678         }
 679         vec3_t Old_N= triNormal(src, OldTriangle[0], OldTriangle[1], OldTriangle[2]);
 680         vec3_t New_N= triNormal(src, NewTriangle[0], NewTriangle[1], NewTriangle[2]);
 681         double OldArea=Old_N.abs();
 682         double NewArea=New_N.abs();
 683         double scal=Old_N*New_N;
 684         double cross=(Old_N.cross(New_N)).abs();//double-cross on Nar Shadaa B-)
 685
 686         cout<<"OldArea="<<OldArea<<endl;
 687         cout<<"NewArea="<<NewArea<<endl;
 688         cout<<"scal="<<scal<<endl;
 689         cout<<"cross="<<cross<<endl;
 690
 691         if(Old_N*New_N<Old_N*Old_N*1./100.)//area + inversion check
 692         {
 693           cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
 694           IsValidSnapPoint=false;
 695         }
 696
 697   /*      if(NewArea<OldArea*1./100.)
 698         {
 699           cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
 700           IsValidSnapPoint=false;
 701         }
 702
 703         if(abs(cross)>10e-4)
 704         {
 705           cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
 706           IsValidSnapPoint=false;
 707         }*/
 708
 709         //mutated cell
 710         MutatedCells.insert(C);
 711         cout<<"cell "<<C<<" has been infected!"<<endl;
 712       }
 713     }
 714
 715     if(N_cells+N_newcells<=0)//survivor check
 716     {
 717       cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
 718       IsValidSnapPoint=false;
 719     }
 720 /*    if(EmptyVolume(DeadNode,PSP))//simplified volume check
 721     {
 722       cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
 723       IsValidSnapPoint=false;
 724     }*/
 725     if(IsValidSnapPoint) {SnapPoint=PSP; break;}
 726   }//end of loop through potential SnapPoints
 727
 728   cout<<"===>SNAPPOINT="<<SnapPoint<<endl;
 729   if(SnapPoint<0) {cout<<"Sorry no possible SnapPoint found."<<endl; return(false);}
 730
 731   //allocate
 732   cout<<"N_points="<<N_points<<endl;
 733   cout<<"N_cells="<<N_cells<<endl;
 734   cout<<"N_newpoints="<<N_newpoints<<endl;
 735   cout<<"N_newcells="<<N_newcells<<endl;
 736   allocateGrid(dst,N_cells+N_newcells,N_points+N_newpoints);
 737
 738   //vector used to redefine the new point IDs
 739   QVector <vtkIdType> OffSet(N_points);
 740
 741   //copy undead points
 742   vtkIdType dst_id_node=0;
 743   for (vtkIdType src_id_node = 0; src_id_node < N_points; src_id_node++) {//loop through src points
 744     if(src_id_node!=DeadNode)//if the node isn't dead, copy it
 745     {
 746       vec3_t x;
 747       src->GetPoints()->GetPoint(src_id_node, x.data());
 748       dst->GetPoints()->SetPoint(dst_id_node, x.data());
 749       copyNodeData(src, src_id_node, dst, dst_id_node);
 750       OffSet[src_id_node]=src_id_node-dst_id_node;
 751       dst_id_node++;
 752     }
 753   };
 754
 755   cout<<"DeadCells="<<DeadCells<<endl;
 756   cout<<"MutatedCells="<<MutatedCells<<endl;
 757   cout<<"MutilatedCells="<<MutilatedCells<<endl;
 758
 759   //Copy undead cells
 760   for (vtkIdType id_cell = 0; id_cell < src->GetNumberOfCells(); ++id_cell) {//loop through src cells
 761     if(!DeadCells.contains(id_cell))//if the cell isn't dead
 762     {
 763       vtkIdType src_N_pts, *src_pts;
 764       vtkIdType dst_N_pts, *dst_pts;
 765       src->GetCellPoints(id_cell, src_N_pts, src_pts);
 766
 767       vtkIdType type_cell = src->GetCellType(id_cell);
 768       cout<<"-->id_cell="<<id_cell<<endl;
 769       for(int i=0;i<src_N_pts;i++) cout<<"src_pts["<<i<<"]="<<src_pts[i]<<endl;
 770 //       src->GetCellPoints(id_cell, dst_N_pts, dst_pts);
 771       dst_N_pts=src_N_pts;
 772       dst_pts=new vtkIdType[dst_N_pts];
 773       if(MutatedCells.contains(id_cell))//mutated cell
 774       {
 775         cout<<"processing mutated cell "<<id_cell<<endl;
 776         for(int i=0;i<src_N_pts;i++)
 777         {
 778           if(src_pts[i]==DeadNode) {
 779             cout<<"SnapPoint="<<SnapPoint<<endl;
 780             cout<<"OffSet[SnapPoint]="<<OffSet[SnapPoint]<<endl;
 781             cout<<"src_pts["<<i<<"]="<<src_pts[i]<<endl;
 782             dst_pts[i]=SnapPoint-OffSet[SnapPoint];
 783           }
 784           else dst_pts[i]=src_pts[i]-OffSet[src_pts[i]];
 785         }
 786         cout<<"--->dst_pts:"<<endl;
 787         for(int i=0;i<dst_N_pts;i++) cout<<"dst_pts["<<i<<"]="<<dst_pts[i]<<endl;
 788
 789       }
 790       else if(MutilatedCells.contains(id_cell))//mutilated cell
 791       {
 792         cout<<"processing mutilated cell "<<id_cell<<endl;
 793
 794         if(type_cell==VTK_QUAD) {
 795           type_cell=VTK_TRIANGLE;
 796           dst_N_pts-=1;
 797         }
 798         else {cout<<"FATAL ERROR: Unknown mutilated cell detected! It is not a quad! Potential xenomorph infestation!"<<endl;EG_BUG;}
 799         //merge points
 800         int j=0;
 801         for(int i=0;i<src_N_pts;i++)
 802         {
 803           if(src_pts[i]==SnapPoint) { dst_pts[j]=SnapPoint-OffSet[SnapPoint];j++; }//SnapPoint
 804           else if(src_pts[i]!=DeadNode) { dst_pts[j]=src_pts[i]-OffSet[src_pts[i]];j++; }//pre-snap/dead + post-snap/dead
 805           //do nothing in case of DeadNode
 806         }
 807       }
 808       else//normal cell
 809       {
 810         cout<<"processing normal cell "<<id_cell<<endl;
 811         for(int i=0;i<src_N_pts;i++)
 812         {
 813           dst_pts[i]=src_pts[i]-OffSet[src_pts[i]];
 814         }
 815       }
 816       //copy the cell
 817       vtkIdType id_new_cell = dst->InsertNextCell(type_cell, dst_N_pts, dst_pts);
 818       copyCellData(src, id_cell, dst, id_new_cell);
 819       cout<<"===Copying cell "<<id_cell<<" to "<<id_new_cell<<"==="<<endl;
 820       cout<<"src_pts:"<<endl;
 821       for(int i=0;i<src_N_pts;i++) cout<<"src_pts["<<i<<"]="<<src_pts[i]<<endl;
 822       cout<<"dst_pts:"<<endl;
 823       for(int i=0;i<dst_N_pts;i++) cout<<"dst_pts["<<i<<"]="<<dst_pts[i]<<endl;
 824       cout<<"OffSet="<<OffSet<<endl;
 825       cout<<"===Copying cell end==="<<endl;
 826       delete dst_pts;
 827     }
 828   };
 829 //   cout_grid(cout,dst,true,true,true,true);
 830   makeCopy(dst, src);
 831   return(true);
 832 }
 833
 834 bool Operation::EmptyVolume(vtkIdType DeadNode, vtkIdType PSP)
 835 {
 836   c2c[DeadNode];
 837   c2c[PSP];
 838   return(true);
 839 }
 840
 841 vec3_t Operation::GetCenter(vtkIdType cellId, double& R)
 842 {
 843   vtkIdType *pts, Npts;
 844   grid->GetCellPoints(cellId, Npts, pts);
 845
 846   vec3_t x(0,0,0);
 847   for (vtkIdType i = 0; i < Npts; ++i) {
 848     vec3_t xp;
 849     grid->GetPoints()->GetPoint(pts[i], xp.data());
 850     x += double(1)/Npts * xp;
 851   };
 852
 853   R = 1e99;
 854   for (vtkIdType i = 0; i < Npts; ++i) {
 855     vec3_t xp;
 856     grid->GetPoints()->GetPoint(pts[i], xp.data());
 857     R = min(R, 0.25*(xp-x).abs());
 858   };
 859
 860   return(x);
 861 }
 862
 863 bool Operation::getNeighbours(vtkIdType Boss, vtkIdType& Peon1, vtkIdType& Peon2, int BC)
 864 {
 865   QVector <vtkIdType> Peons;
 866
 867   QSet <int> S1=n2c[Boss];
 868   cout<<"S1="<<S1<<endl;
 869   foreach(vtkIdType PN,n2n[Boss])
 870   {
 871     cout<<"PN="<<PN<<endl;
 872     QSet <int> S2=n2c[PN];
 873     cout<<"S2="<<S2<<endl;
 874     QSet <int> Si=S2.intersect(S1);
 875     cout<<"PN="<<PN<<" Si="<<Si<<endl;
 876     if(Si.size()<2)//only one common cell
 877     {
 878       Peons.push_back(PN);
 879     }
 880     else
 881     {
 882       QSet <int> bc_set;
 883       foreach(vtkIdType C,Si)
 884       {
 885         EG_VTKDCC(vtkIntArray, cell_code, grid, "cell_code");
 886         int bc=cell_code->GetValue(C);
 887         cout<<"C="<<C<<" bc="<<bc<<endl;
 888         bc_set.insert(bc);
 889       }
 890       if(bc_set.size()>1)//2 different boundary codes
 891       {
 892         Peons.push_back(PN);
 893       }
 894     }
 895   }
 896   if(Peons.size()==2)
 897   {
 898     Peon1=Peons[0];
 899     Peon2=Peons[1];
 900     return(true);
 901   }
 902   else
 903   {
 904     cout<<"FATAL ERROR: number of neighbours != 2"<<endl;
 905     EG_BUG;
 906   }
 907   return(false);
 908 }
 909
 910 int Operation::UpdateMeshDensity()
 911 {
 912   cout<<"===UpdateMeshDensity START==="<<endl;
 913
 914   getAllSurfaceCells(cells,grid);
 915   EG_VTKDCC(vtkIntArray, cell_code, grid, "cell_code");
 916   EG_VTKDCN(vtkDoubleArray, node_meshdensity, grid, "node_meshdensity");
 917   getNodesFromCells(cells, nodes, grid);
 918   setGrid(grid);
 919   setCells(cells);
 920
 921   if(DebugLevel>5) cout<<"cells.size()="<<cells.size()<<endl;
 922
 923   EG_VTKDCN(vtkDoubleArray, node_meshdensity_current, grid, "node_meshdensity_current");
 924   foreach(vtkIdType node,nodes)
 925   {
 926     double L=CurrentVertexAvgDist(node,n2n,grid);
 927     double D=1./L;
 928     node_meshdensity_current->SetValue(node, D);
 929   }
 930   cout<<"===UpdateMeshDensity END==="<<endl;
 931   return(0);
 932 }