src/surfacesmoother.cpp

   1 #include "surfacesmoother.h"
   2
   3 #include <QString>
   4 #include <QTextStream>
   5 #include <vtkCharArray.h>
   6
   7 #include "smoothingutilities.h"
   8
   9 #include "swaptriangles.h"
  10 #include "laplacesmoother.h"
  11 #include "guimainwindow.h"
  12
  13 #include <iostream>
  14 using namespace std;
  15
  16 SurfaceSmoother::SurfaceSmoother()
  17 {
  18    DebugLevel=0;
  19 }
  20
  21 int SurfaceSmoother::Process()
  22 {
  23   int i_iter=0;
  24   for(i_iter=0;i_iter<NumberOfIterations;i_iter++)//TODO:Optimize this loop
  25   {
  26     cout<<"===ITERATION NB "<<i_iter<<"/"<<NumberOfIterations<<"==="<<endl;
  27
  28     m_total_N_newpoints=0;
  29     m_total_N_newcells=0;
  30
  31     getAllSurfaceCells(m_AllCells,m_grid);
  32     getSurfaceCells(m_bcs, m_SelectedCells, m_grid);
  33     cout<<"m_AllCells.size()="<<m_AllCells.size()<<endl;
  34
  35     EG_VTKDCC(vtkIntArray, cell_code, m_grid, "cell_code");
  36
  37     m_SelectedNodes.clear();
  38     getSurfaceNodes(m_bcs,m_SelectedNodes,m_grid);
  39     getNodesFromCells(m_AllCells, nodes, m_grid);
  40     setGrid(m_grid);
  41     setCells(m_AllCells);
  42
  43     cout<<"m_AllCells.size()="<<m_AllCells.size()<<endl;
  44
  45     //Phase D: edit points
  46     cout<<"===Phase D==="<<endl;
  47     N_inserted_FP=0;
  48     N_inserted_EP=0;
  49     N_removed_FP=0;
  50     N_removed_EP=0;
  51
  52     //Method 1
  53 //     FullEdit();
  54
  55     //Method 2
  56 /*    if(insert_FP) insert_FP_all();
  57     if(insert_EP) insert_EP_all();
  58     if(remove_FP) remove_FP_all();
  59     if(remove_EP) remove_EP_all();*/
  60
  61     //Method 3
  62     if(insert_FP) {
  63       UpdateDesiredMeshDensity();
  64       insert_FP_all();
  65       DualSave("/data1/home/mtaverne/Geometries/simulations/SurfaceTests/insert_FP-post-insert");
  66       if(DoSwap) SwapFunction();
  67       DualSave("/data1/home/mtaverne/Geometries/simulations/SurfaceTests/insert_FP-post-swap-1");
  68       if(DoLaplaceSmoothing) SmoothFunction();
  69       DualSave("/data1/home/mtaverne/Geometries/simulations/SurfaceTests/insert_FP-post-laplace");
  70       if(DoSwap) SwapFunction();
  71       DualSave("/data1/home/mtaverne/Geometries/simulations/SurfaceTests/insert_FP-post-swap-2");
  72     }
  73     DualSave("/data1/home/mtaverne/Geometries/simulations/SurfaceTests/post-insert_FP");
  74
  75     if(insert_EP) {
  76       UpdateDesiredMeshDensity();
  77       insert_EP_all();
  78       DualSave("/data1/home/mtaverne/Geometries/simulations/SurfaceTests/insert_EP-post-insert");
  79       if(DoSwap) SwapFunction();
  80       DualSave("/data1/home/mtaverne/Geometries/simulations/SurfaceTests/insert_EP-post-swap");
  81       if(DoLaplaceSmoothing) SmoothFunction();
  82       DualSave("/data1/home/mtaverne/Geometries/simulations/SurfaceTests/insert_EP-post-laplace");
  83     }
  84     DualSave("/data1/home/mtaverne/Geometries/simulations/SurfaceTests/post-insert_EP");
  85
  86     if(remove_FP) {
  87       UpdateDesiredMeshDensity();
  88       remove_FP_all_3();
  89       if(DoSwap) SwapFunction();
  90       if(DoLaplaceSmoothing) SmoothFunction();
  91     }
  92     DualSave("/data1/home/mtaverne/Geometries/simulations/SurfaceTests/post-remove_FP");
  93
  94     if(remove_EP) {
  95       UpdateDesiredMeshDensity();
  96       remove_EP_all_3();
  97       if(DoSwap) SwapFunction();
  98       if(DoLaplaceSmoothing) SmoothFunction();
  99     }
 100     DualSave("/data1/home/mtaverne/Geometries/simulations/SurfaceTests/post-remove_EP");
 101
 102 /*    if(DoSwap) SwapFunction();
 103     if(DoLaplaceSmoothing) SmoothFunction();*/
 104
 105     cout<<"===Summary==="<<endl;
 106     cout<<"N_inserted_FP="<<N_inserted_FP<<endl;
 107     cout<<"N_inserted_EP="<<N_inserted_EP<<endl;
 108     cout<<"N_removed_FP="<<N_removed_FP<<endl;
 109     cout<<"N_removed_EP="<<N_removed_EP<<endl;
 110
 111     cout<<"N_points="<<N_points<<endl;
 112     cout<<"N_cells="<<N_cells<<endl;
 113     cout<<"m_total_N_newpoints="<<m_total_N_newpoints<<endl;
 114     cout<<"m_total_N_newcells="<<m_total_N_newcells<<endl;
 115     cout<<"============"<<endl;
 116
 117 //     if(m_total_N_newpoints==0 && m_total_N_newcells==0) break;
 118     if(N_inserted_FP==0 && N_inserted_EP==0 && N_removed_FP==0 && N_removed_EP==0) break;
 119   }
 120
 121   cout<<"i_iter/NumberOfIterations="<<i_iter<<"/"<<NumberOfIterations<<endl;
 122   UpdateDesiredMeshDensity();
 123   UpdateMeshDensity();
 124   if(i_iter<NumberOfIterations) cout<<"WARNING: Exited before finishing all iterations."<<endl;
 125   return 1;
 126 }
 127 //end of process
 128
 129 int SurfaceSmoother::UpdateDesiredMeshDensity()
 130 {
 131     //Phase B : define desired mesh density
 132   cout<<"=== UpdateDesiredMeshDensity ==="<<endl;
 133
 134   getAllSurfaceCells(m_AllCells,m_grid);
 135   getSurfaceCells(m_bcs, m_SelectedCells, m_grid);
 136   cout<<"m_AllCells.size()="<<m_AllCells.size()<<endl;
 137
 138   EG_VTKDCC(vtkIntArray, cell_code, m_grid, "cell_code");
 139
 140   m_SelectedNodes.clear();
 141   getSurfaceNodes(m_bcs,m_SelectedNodes,m_grid);
 142   getNodesFromCells(m_AllCells, nodes, m_grid);
 143   getNodesFromCells(m_AllCells, m_AllNodes, m_grid);
 144
 145   setGrid(m_grid);
 146   setCells(m_AllCells);
 147
 148   cout<<"m_AllCells.size()="<<m_AllCells.size()<<endl;
 149
 150   UpdateNodeType_all();
 151   EG_VTKDCN(vtkCharArray, node_type, m_grid, "node_type");
 152   EG_VTKDCN(vtkDoubleArray, node_meshdensity, m_grid, "node_meshdensity");
 153   EG_VTKDCN(vtkIntArray, node_specified_density, m_grid, "node_specified_density");
 154
 155 /*  //Phase A : Calculate current mesh density
 156   cout<<"===Phase A==="<<endl;
 157
 158   foreach(vtkIdType node,m_SelectedNodes)
 159   {
 160     VertexMeshDensity nodeVMD = getVMD(node,node_type->GetValue(node));
 161     int idx=VMDvector.indexOf(nodeVMD);
 162     if(DebugLevel>3) cout<<"idx="<<idx<<endl;
 163     if(idx!=-1)//specified
 164     {
 165       node_meshdensity->SetValue(node, VMDvector[idx].density);
 166     }
 167     else//unspecified
 168     {
 169       double L=CurrentVertexAvgDist(node,n2n,m_grid);
 170       double D=1./L;
 171       node_meshdensity->SetValue(node, D);
 172     }
 173   }*/
 174
 175   double diff=Convergence_meshdensity+1;
 176   if(DebugLevel>3) cout<<"before loop: diff="<<diff<<endl;
 177   bool first=true;
 178   int iter=0;
 179   do {
 180     if(DebugLevel>2) cout<<"--->diff="<<diff<<endl;
 181     first=true;
 182     foreach(vtkIdType node,m_AllNodes)
 183     {
 184       if(DebugLevel>2) cout<<"======>"<<endl;
 185       VertexMeshDensity nodeVMD = getVMD(node,node_type->GetValue(node));
 186       int idx=VMDvector.indexOf(nodeVMD);
 187       node_specified_density->SetValue(node, idx);
 188       if(DebugLevel>2) cout<<"------>idx="<<idx<<endl;
 189       if(idx!=-1)//specified
 190       {
 191         node_meshdensity->SetValue(node, VMDvector[idx].density);
 192       }
 193       else//unspecified
 194       {
 195         double D=DesiredMeshDensity(node,n2n,m_grid);
 196         if(first) {
 197           if(DebugLevel>2) {
 198             cout<<"------>FIRST:"<<endl;
 199             cout<<"------>D="<<D<<endl;
 200             cout<<"------>node_meshdensity->GetValue("<<node<<")="<<node_meshdensity->GetValue(node)<<endl;
 201             cout<<"------>D-node_meshdensity->GetValue("<<node<<")="<<D-node_meshdensity->GetValue(node)<<endl;
 202             cout<<"------>diff=abs(D-node_meshdensity->GetValue("<<node<<"))="<<abs(D-node_meshdensity->GetValue(node))<<endl;
 203           }
 204           diff=abs(D-node_meshdensity->GetValue(node));
 205           first=false;
 206         }
 207         else {
 208           if(DebugLevel>2) {
 209             cout<<"------>NOT FIRST:"<<endl;
 210             cout<<"------>D="<<D<<endl;
 211             cout<<"------>node_meshdensity->GetValue("<<node<<")="<<node_meshdensity->GetValue(node)<<endl;
 212             cout<<"------>D-node_meshdensity->GetValue("<<node<<")="<<D-node_meshdensity->GetValue(node)<<endl;
 213             cout<<"------>diff=abs(D-node_meshdensity->GetValue("<<node<<"))="<<abs(D-node_meshdensity->GetValue(node))<<endl;
 214             cout<<"------>diff="<<diff<<endl;
 215             cout<<"------>max(abs(D-node_meshdensity->GetValue("<<node<<")),diff)="<<max(abs(D-node_meshdensity->GetValue(node)),diff)<<endl;
 216           }
 217           diff=max(abs(D-node_meshdensity->GetValue(node)),diff);
 218         }
 219         node_meshdensity->SetValue(node, D);
 220       }
 221       if(DebugLevel>2) cout<<"======>"<<endl;
 222     }
 223     iter++;
 224   } while(diff>Convergence_meshdensity && !first && iter<maxiter_density);// if first=true, it means no new mesh density has been defined (all densities specified)
 225   cout<<"iter="<<iter<<endl;
 226   if(iter>=maxiter_density) cout<<"WARNING: Desired convergence factor has not been reached!"<<endl;
 227   return(0);
 228 }
 229
 230 int SurfaceSmoother::SwapFunction()
 231 {
 232   //Phase E : Delaunay swap
 233   QSet<int> bcs_complement=complementary_bcs(m_bcs,m_grid,cells);
 234   cout<<"m_bcs="<<m_bcs<<endl;
 235   cout<<"bcs_complement="<<bcs_complement<<endl;
 236
 237   SwapTriangles swap;
 238   swap.setGrid(m_grid);
 239   swap.setBoundaryCodes(bcs_complement);
 240   swap();
 241   return(0);
 242 }
 243
 244 int SurfaceSmoother::SmoothFunction()
 245 {
 246   cout<<"=== SmoothFunction START ==="<<endl;
 247   //Phase F : translate points to smooth grid
 248   //4 possibilities
 249   //vtk smooth 1
 250   //vtk smooth 2
 251   //laplacian smoothing with projection
 252   //Roland smoothing with projection
 253
 254   //laplacian smoothing with projection
 255   LaplaceSmoother Lap;
 256   Lap.SetInput(m_bcs,m_grid);
 257   Lap.SetNumberOfIterations(N_SmoothIterations);
 258   Lap();
 259   cout<<"=== SmoothFunction END ==="<<endl;
 260   return(0);
 261 }
 262
 263 VertexMeshDensity SurfaceSmoother::getVMD(vtkIdType node, char VertexType)
 264 {
 265   VertexMeshDensity VMD;
 266   VMD.type=VertexType;
 267   VMD.density=0;
 268   VMD.CurrentNode=node;
 269   EG_VTKDCC(vtkIntArray, cell_code, m_grid, "cell_code");
 270 /*  createNodeMapping(nodes, _nodes, m_grid);
 271   createNodeToCell(m_AllCells, nodes, _nodes, n2c, m_grid);*/
 272
 273   QSet <int> bc;
 274   foreach(vtkIdType C, n2c[node])
 275   {
 276     bc.insert(cell_code->GetValue(C));
 277     VMD.BCmap[cell_code->GetValue(C)]=2;
 278   }
 279   VMD.BClist.resize(bc.size());
 280   qCopy(bc.begin(),bc.end(),VMD.BClist.begin());
 281   qSort(VMD.BClist.begin(),VMD.BClist.end());
 282   return(VMD);
 283 }
 284
 285 int SurfaceSmoother::insert_FP_counter()
 286 {
 287   cout<<"===insert_FP_counter() START==="<<endl;
 288   foreach(vtkIdType id_cell, m_SelectedCells)
 289   {
 290     if( !marked_cells[id_cell] && insert_fieldpoint(id_cell) )
 291     {
 292       if(DebugLevel>0) cout<<"inserting a field point "<<id_cell<<endl;
 293       N_inserted_FP++;
 294       marked_cells[id_cell]=true;
 295       N_newcells+=2;
 296       N_newpoints+=1;
 297     }
 298   }
 299   cout<<"===insert_FP_counter() END==="<<endl;
 300   return(0);
 301 }
 302
 303 int SurfaceSmoother::insert_EP_counter()
 304 {
 305   cout<<"===insert_EP_counter() START==="<<endl;
 306
 307   //Phase C: Prepare edge_map
 308   cout<<"===Phase C==="<<endl;
 309   edge_map.clear();
 310   vtkIdType edgeId=1;
 311   foreach(vtkIdType node1,m_SelectedNodes)
 312   {
 313 //       cout<<"node1="<<node1<<endl;
 314     foreach(vtkIdType node2,n2n[node1])
 315     {
 316       if(edge_map[OrderedPair(node1,node2)]==0) { //this edge hasn't been numbered yet
 317         edge_map[OrderedPair(node1,node2)]=edgeId;edgeId++;
 318       }
 319     }
 320   }
 321   cout<<"edge_map.size()="<<edge_map.size()<<endl;
 322
 323   StencilVector.clear();
 324   QMapIterator< pair<vtkIdType,vtkIdType>, vtkIdType> edge_map_iter(edge_map);
 325       //rewind the iterator
 326   edge_map_iter.toFront ();
 327       //start loop
 328   while (edge_map_iter.hasNext()) {
 329     edge_map_iter.next();
 330     vtkIdType node1=edge_map_iter.key().first;
 331     vtkIdType node2=edge_map_iter.key().second;
 332     if(DebugLevel>10) cout << "--->(" << node1 << "," << node2 << ")" << ": " << edge_map_iter.value() << endl;
 333     QSet <int> stencil_cells_set;
 334     QVector <int> stencil_cells_vector;
 335     stencil_cells_set=n2c[node1];
 336     stencil_cells_set.intersect(n2c[node2]);
 337     if(DebugLevel>10) cout<<"stencil_cells_set="<<stencil_cells_set<<endl;
 338
 339     stencil_cells_vector.resize(stencil_cells_set.size());
 340     qCopy(stencil_cells_set.begin(),stencil_cells_set.end(),stencil_cells_vector.begin());
 341     if(DebugLevel>10) cout<<"stencil_cells_vector="<<stencil_cells_vector<<endl;
 342
 343     vtkIdType id_cell=stencil_cells_vector[0];
 344     int SideToSplit = getSide(id_cell,m_grid,node1,node2);
 345     if(DebugLevel>10) cout<<"SideToSplit="<<SideToSplit<<endl;
 346     if(DebugLevel>10) cout<<"c2c[id_cell][SideToSplit]="<<c2c[id_cell][SideToSplit]<<endl;
 347     if(DebugLevel>10) for(int i=0;i<3;i++) cout<<"c2c[id_cell]["<<i<<"]="<<c2c[id_cell][i]<<endl;
 348     stencil_t S=getStencil(id_cell,SideToSplit);
 349
 350     bool stencil_marked=false;
 351     foreach(vtkIdType C,stencil_cells_vector)
 352     {
 353       if(marked_cells[C]) stencil_marked=true;
 354     }
 355     if(DebugLevel>10) cout<<"stencil_marked="<<stencil_marked<<endl;
 356     if(DebugLevel>10) cout<<"insert_edgepoint(node1,node2)="<<insert_edgepoint(node1,node2)<<endl;
 357
 358     if( !stencil_marked && insert_edgepoint(node1,node2) )
 359     {
 360       if(DebugLevel>1) cout<<"inserting an edge point "<< "(" << node1 << "," << node2 << ")" << ": " << edge_map_iter.value() << endl;
 361       N_inserted_EP++;
 362       foreach(vtkIdType C,stencil_cells_vector) marked_cells[C]=true;
 363       StencilVector.push_back(S);
 364
 365       if(stencil_cells_vector.size()==2)//2 cells around the edge
 366       {
 367         N_newcells+=2;
 368         N_newpoints+=1;
 369       }
 370       else//1 cell around the edge
 371       {
 372         N_newcells+=1;
 373         N_newpoints+=1;
 374       }
 375     }
 376     if(DebugLevel>10) cout <<"--->end of edge processing"<<endl;
 377   }
 378   cout<<"===insert_EP_counter() END==="<<endl;
 379   return(0);
 380 }
 381
 382 int SurfaceSmoother::remove_FP_counter()
 383 {
 384   cout<<"===remove_FP_counter() START==="<<endl;
 385   cout<<"marked_cells="<<marked_cells<<endl;
 386 //   cout<<"hitlist="<<hitlist<<endl;
 387   cout<<"hitlist.size()="<<hitlist.size()<<endl;
 388   cout<<"N_newcells="<<N_newcells<<endl;
 389   cout<<"N_newpoints="<<N_newpoints<<endl;
 390   cout<<"N_removed_FP="<<N_removed_FP<<endl;
 391
 392   UpdateNodeType_all();
 393   EG_VTKDCN(vtkCharArray, node_type, m_grid, "node_type");
 394   foreach(vtkIdType node,m_SelectedNodes)
 395   {
 396     if(node_type->GetValue(node)==VTK_SIMPLE_VERTEX)
 397     {
 398       bool marked=false;
 399       foreach(vtkIdType C,n2c[node])
 400       {
 401         if(marked_cells[C]) marked=true;
 402       }
 403
 404       QSet <vtkIdType> DeadCells;
 405       QSet <vtkIdType> MutatedCells;
 406       QSet <vtkIdType> MutilatedCells;
 407       if( !marked && remove_fieldpoint(node) && FindSnapPoint(m_grid,node,DeadCells,MutatedCells,MutilatedCells, N_newpoints, N_newcells)!=-1)
 408       {
 409         if(DebugLevel>1) cout<<"removing field point "<<node<<endl;
 410         N_removed_FP++;
 411         hitlist[node]=1;
 412         foreach(vtkIdType C,n2c[node]) marked_cells[C]=true;
 413         N_newcells-=2;
 414         N_newpoints-=1;
 415       }
 416     }
 417   }
 418   cout<<"===remove_FP_counter() END==="<<endl;
 419   return(0);
 420 }
 421
 422 int SurfaceSmoother::remove_EP_counter()
 423 {
 424   cout<<"===remove_EP_counter() START==="<<endl;
 425   UpdateNodeType_all();
 426   EG_VTKDCN(vtkCharArray, node_type, m_grid, "node_type");
 427   foreach(vtkIdType node,m_SelectedNodes)
 428   {
 429     if(node_type->GetValue(node)==VTK_BOUNDARY_EDGE_VERTEX)
 430     {
 431       bool marked=false;
 432       foreach(vtkIdType C,n2c[node])
 433       {
 434         if(marked_cells[C]) marked=true;
 435       }
 436       QSet <vtkIdType> DeadCells;
 437       QSet <vtkIdType> MutatedCells;
 438       QSet <vtkIdType> MutilatedCells;
 439       if( !marked && remove_edgepoint(node) && FindSnapPoint(m_grid,node,DeadCells,MutatedCells,MutilatedCells, N_newpoints, N_newcells)!=-1)
 440       {
 441         cout<<"removing edge point "<<node<<endl;
 442         N_removed_EP++;
 443         hitlist[node]=2;
 444         foreach(vtkIdType C,n2c[node]) marked_cells[C]=true;
 445         if(n2n[node].size()==4)//4 cells around the edge
 446         {
 447           N_newcells-=2;
 448           N_newpoints-=1;
 449         }
 450         else//2 cells around the edge
 451         {
 452           N_newcells-=1;
 453           N_newpoints-=1;
 454         }
 455       }
 456     }
 457   }
 458   cout<<"===remove_EP_counter() END==="<<endl;
 459   return(0);
 460 }
 461
 462 int SurfaceSmoother::insert_FP_actor(vtkUnstructuredGrid* grid_tmp)
 463 {
 464   cout<<"===insert_FP_actor START==="<<endl;
 465
 466   EG_VTKDCC(vtkIntArray, cell_code_tmp, grid_tmp, "cell_code");
 467   foreach(vtkIdType id_cell, m_SelectedCells)
 468   {
 469 /*    if(marked_cells[id_cell]) cout<<"--->marked_cells["<<id_cell<<"]=TRUE"<<endl;
 470     else cout<<"--->marked_cells["<<id_cell<<"]=FALSE"<<endl;*/
 471
 472     if( !marked_cells[id_cell] && insert_fieldpoint(id_cell) )
 473     {
 474       if(DebugLevel>0) cout<<"inserting a field point "<<id_cell<<endl;
 475       vtkIdType newBC=cell_code_tmp->GetValue(id_cell);
 476       if(DebugLevel>42) cout<<"id_cell="<<id_cell<<" newBC="<<newBC<<endl;
 477
 478       vtkIdType N_pts, *pts;
 479       m_grid->GetCellPoints(id_cell, N_pts, pts);
 480       vec3_t C(0,0,0);
 481
 482       int N_neighbours=N_pts;
 483       if(DebugLevel>42) cout<<"N_neighbours="<<N_neighbours<<endl;
 484       vec3_t corner[4];
 485       vtkIdType pts_triangle[4][3];
 486       for(int i=0;i<N_neighbours;i++)
 487       {
 488         m_grid->GetPoints()->GetPoint(pts[i], corner[i].data());
 489         C+=corner[i];
 490       }
 491       C=(1/(double)N_neighbours)*C;
 492       addPoint(grid_tmp,m_newNodeId,C.data());
 493       vtkIdType intmidpoint=m_newNodeId;
 494       m_newNodeId++;
 495
 496       for(int i=0;i<N_neighbours;i++)
 497       {
 498         pts_triangle[i][0]=pts[i];
 499         pts_triangle[i][1]=pts[(i+1)%N_neighbours];
 500         pts_triangle[i][2]=intmidpoint;
 501         if(i==0)
 502         {
 503           grid_tmp->ReplaceCell(id_cell , 3, pts_triangle[0]);
 504           cell_code_tmp->SetValue(id_cell, newBC);
 505         }
 506         else
 507         {
 508           vtkIdType newCellId = grid_tmp->InsertNextCell(VTK_TRIANGLE,3,pts_triangle[i]);
 509           cell_code_tmp->SetValue(newCellId, newBC);
 510         }
 511       }
 512
 513     }
 514   }
 515   cout<<"===insert_FP_actor END==="<<endl;
 516   return(0);
 517 }
 518
 519 int SurfaceSmoother::insert_EP_actor(vtkUnstructuredGrid* grid_tmp)
 520 {
 521   cout<<"===insert_EP_actor START==="<<endl;
 522
 523   EG_VTKDCC(vtkIntArray, cell_code_tmp, grid_tmp, "cell_code");
 524   foreach(stencil_t S,StencilVector)
 525   {
 526     if(DebugLevel>10) cout<<"S="<<S<<endl;
 527     vec3_t A,B;
 528     grid_tmp->GetPoint(S.p[1],A.data());
 529     grid_tmp->GetPoint(S.p[3],B.data());
 530     vec3_t M=0.5*(A+B);
 531     addPoint(grid_tmp,m_newNodeId,M.data());
 532     if(DebugLevel>0) cout<<"NEW EDGE POINT: "<<m_newNodeId<<endl;
 533
 534     vtkIdType pts_triangle[4][3];
 535
 536     if(S.valid){//there is a neighbour cell
 537       if(DebugLevel>10) cout<<"marked_cells["<<S.id_cell1<<"]=true;"<<endl;
 538       if(DebugLevel>10) cout<<"marked_cells["<<S.id_cell2<<"]=true;"<<endl;
 539       marked_cells[S.id_cell1]=true;
 540       marked_cells[S.id_cell2]=true;
 541
 542       for(int i=0;i<4;i++)
 543       {
 544         pts_triangle[i][0]=S.p[i];
 545         pts_triangle[i][1]=S.p[(i+1)%4];
 546         pts_triangle[i][2]=m_newNodeId;
 547       }
 548
 549       int bc1=cell_code_tmp->GetValue(S.id_cell1);
 550       int bc2=cell_code_tmp->GetValue(S.id_cell2);
 551
 552       grid_tmp->ReplaceCell(S.id_cell1 , 3, pts_triangle[0]);
 553       cell_code_tmp->SetValue(S.id_cell1, bc1);
 554
 555       grid_tmp->ReplaceCell(S.id_cell2 , 3, pts_triangle[1]);
 556       cell_code_tmp->SetValue(S.id_cell2, bc2);
 557
 558       vtkIdType newCellId;
 559       newCellId = grid_tmp->InsertNextCell(VTK_TRIANGLE,3,pts_triangle[2]);
 560       cell_code_tmp->SetValue(newCellId, bc2);
 561       newCellId = grid_tmp->InsertNextCell(VTK_TRIANGLE,3,pts_triangle[3]);
 562       cell_code_tmp->SetValue(newCellId, bc1);
 563     }
 564     else{//there is no neighbour cell
 565       if(DebugLevel>10) cout<<"marked_cells["<<S.id_cell1<<"]=true;"<<endl;
 566       marked_cells[S.id_cell1]=true;
 567
 568       pts_triangle[0][0]=S.p[0];
 569       pts_triangle[0][1]=S.p[1];
 570       pts_triangle[0][2]=m_newNodeId;
 571       pts_triangle[3][0]=S.p[3];
 572       pts_triangle[3][1]=S.p[0];
 573       pts_triangle[3][2]=m_newNodeId;
 574
 575       int bc1=cell_code_tmp->GetValue(S.id_cell1);
 576
 577       grid_tmp->ReplaceCell(S.id_cell1 , 3, pts_triangle[0]);
 578       cell_code_tmp->SetValue(S.id_cell1, bc1);
 579
 580       vtkIdType newCellId;
 581       newCellId = grid_tmp->InsertNextCell(VTK_TRIANGLE,3,pts_triangle[3]);
 582       cell_code_tmp->SetValue(newCellId, bc1);
 583     }
 584
 585     m_newNodeId++;
 586   }
 587   cout<<"===insert_EP_actor END==="<<endl;
 588   return(0);
 589 }
 590
 591 int SurfaceSmoother::remove_FP_actor(vtkUnstructuredGrid* grid_tmp)
 592 {
 593   cout<<"===remove_FP_actor START==="<<endl;
 594   abort();
 595
 596   foreach(vtkIdType node,m_SelectedNodes)
 597   {
 598     if(hitlist[node]==1)
 599     {
 600
 601     }
 602     bool marked=false;
 603     foreach(vtkIdType C,n2c[node])
 604     {
 605       if(marked_cells[C]) marked=true;
 606     }
 607     if( !marked && remove_fieldpoint(node) )
 608     {
 609       if(DebugLevel>1) cout<<"removing field point "<<node<<endl;
 610       foreach(vtkIdType C,n2c[node]) marked_cells[C]=true;
 611       //TODO: Special copy function, leaving out nodes to remove
 612     }
 613   }
 614   cout<<"===remove_FP_actor END==="<<endl;
 615   return(0);
 616 }
 617
 618 int SurfaceSmoother::remove_EP_actor(vtkUnstructuredGrid* grid_tmp)
 619 {
 620   cout<<"===remove_EP_actor START==="<<endl;
 621   abort();
 622
 623   foreach(vtkIdType node,m_SelectedNodes)
 624   {
 625     bool marked=false;
 626     foreach(vtkIdType C,n2c[node])
 627     {
 628       if(marked_cells[C]) marked=true;
 629     }
 630     if( !marked && remove_edgepoint(node) )
 631     {
 632       if(DebugLevel>1) cout<<"removing edge point "<<node<<endl;
 633       foreach(vtkIdType C,n2c[node]) marked_cells[C]=true;
 634       if(n2n[node].size()==4)//4 cells around the edge
 635       {
 636
 637       }
 638       else//2 cells around the edge
 639       {
 640
 641       }
 642     }
 643   }
 644   cout<<"===remove_EP_actor END==="<<endl;
 645   return(0);
 646 }
 647
 648 int SurfaceSmoother::insert_FP_all()
 649 {
 650   cout<<"===insert_FP_all START==="<<endl;
 651
 652   getAllSurfaceCells(m_AllCells,m_grid);
 653   getSurfaceCells(m_bcs, m_SelectedCells, m_grid);
 654   EG_VTKDCC(vtkIntArray, cell_code, m_grid, "cell_code");
 655   EG_VTKDCN(vtkDoubleArray, node_meshdensity, m_grid, "node_meshdensity");
 656   m_SelectedNodes.clear();
 657   getSurfaceNodes(m_bcs,m_SelectedNodes,m_grid);
 658   getNodesFromCells(m_AllCells, nodes, m_grid);
 659   setGrid(m_grid);
 660   setCells(m_AllCells);
 661   cout<<"m_AllCells.size()="<<m_AllCells.size()<<endl;
 662
 663   N_inserted_FP=0;
 664
 665   N_points=m_grid->GetNumberOfPoints();
 666   N_cells=m_grid->GetNumberOfCells();
 667   N_newpoints=0;
 668   N_newcells=0;
 669
 670   marked_cells.clear();
 671   marked_nodes.clear();
 672
 673   insert_FP_counter();
 674
 675     //unmark cells (TODO: optimize)
 676   marked_cells.clear();
 677     //init grid_tmp
 678   N_points=m_grid->GetNumberOfPoints();
 679   N_cells=m_grid->GetNumberOfCells();
 680   cout<<"N_points="<<N_points<<endl;
 681   cout<<"N_cells="<<N_cells<<endl;
 682   cout<<"N_newpoints="<<N_newpoints<<endl;
 683   cout<<"N_newcells="<<N_newcells<<endl;
 684   EG_VTKSP(vtkUnstructuredGrid,grid_tmp);
 685   allocateGrid(grid_tmp,N_cells+N_newcells,N_points+N_newpoints);
 686   m_total_N_newpoints+=N_newpoints; m_total_N_newcells+=N_newcells;
 687
 688   makeCopyNoAlloc(m_grid, grid_tmp);
 689     //initialize new node counter
 690   m_newNodeId=N_points;
 691
 692   insert_FP_actor(grid_tmp);
 693
 694   makeCopy(grid_tmp,m_grid);
 695   cout<<"===insert_FP_all END==="<<endl;
 696   return(0);
 697 }
 698
 699 int SurfaceSmoother::insert_EP_all()
 700 {
 701   cout<<"===insert_EP_all START==="<<endl;
 702
 703   getAllSurfaceCells(m_AllCells,m_grid);
 704   getSurfaceCells(m_bcs, m_SelectedCells, m_grid);
 705   EG_VTKDCC(vtkIntArray, cell_code, m_grid, "cell_code");
 706   EG_VTKDCN(vtkDoubleArray, node_meshdensity, m_grid, "node_meshdensity");
 707   m_SelectedNodes.clear();
 708   getSurfaceNodes(m_bcs,m_SelectedNodes,m_grid);
 709   getNodesFromCells(m_AllCells, nodes, m_grid);
 710   setGrid(m_grid);
 711   setCells(m_AllCells);
 712   cout<<"m_AllCells.size()="<<m_AllCells.size()<<endl;
 713
 714   N_inserted_EP=0;
 715
 716   N_points=m_grid->GetNumberOfPoints();
 717   N_cells=m_grid->GetNumberOfCells();
 718   N_newpoints=0;
 719   N_newcells=0;
 720
 721   marked_cells.clear();
 722   marked_nodes.clear();
 723
 724   insert_EP_counter();
 725
 726     //unmark cells (TODO: optimize)
 727   marked_cells.clear();
 728     //init grid_tmp
 729   N_points=m_grid->GetNumberOfPoints();
 730   N_cells=m_grid->GetNumberOfCells();
 731   cout<<"N_points="<<N_points<<endl;
 732   cout<<"N_cells="<<N_cells<<endl;
 733   cout<<"N_newpoints="<<N_newpoints<<endl;
 734   cout<<"N_newcells="<<N_newcells<<endl;
 735   EG_VTKSP(vtkUnstructuredGrid,grid_tmp);
 736   allocateGrid(grid_tmp,N_cells+N_newcells,N_points+N_newpoints);
 737   m_total_N_newpoints+=N_newpoints; m_total_N_newcells+=N_newcells;
 738
 739   makeCopyNoAlloc(m_grid, grid_tmp);
 740     //initialize new node counter
 741   m_newNodeId=N_points;
 742
 743   insert_EP_actor(grid_tmp);
 744
 745   makeCopy(grid_tmp,m_grid);
 746
 747   cout<<"===insert_EP_all END==="<<endl;
 748   return(0);
 749 }
 750
 751 int SurfaceSmoother::remove_FP_all()
 752 {
 753   cout<<"===remove_FP_all START==="<<endl;
 754
 755   getAllSurfaceCells(m_AllCells,m_grid);
 756   getSurfaceCells(m_bcs, m_SelectedCells, m_grid);
 757   EG_VTKDCC(vtkIntArray, cell_code, m_grid, "cell_code");
 758   EG_VTKDCN(vtkDoubleArray, node_meshdensity, m_grid, "node_meshdensity");
 759   m_SelectedNodes.clear();
 760   getSurfaceNodes(m_bcs,m_SelectedNodes,m_grid);
 761   getNodesFromCells(m_AllCells, nodes, m_grid);
 762   setGrid(m_grid);
 763   setCells(m_AllCells);
 764   cout<<"m_AllCells.size()="<<m_AllCells.size()<<endl;
 765
 766   N_removed_FP=0;
 767
 768   N_points=m_grid->GetNumberOfPoints();
 769   N_cells=m_grid->GetNumberOfCells();
 770   N_newpoints=0;
 771   N_newcells=0;
 772
 773   hitlist.resize(N_points);
 774   offset.resize(N_points);
 775
 776   marked_cells.clear();
 777   marked_nodes.clear();
 778
 779   remove_FP_counter();
 780
 781     //unmark cells (TODO: optimize)
 782   marked_cells.clear();
 783     //init grid_tmp
 784   N_points=m_grid->GetNumberOfPoints();
 785   N_cells=m_grid->GetNumberOfCells();
 786   cout<<"N_points="<<N_points<<endl;
 787   cout<<"N_cells="<<N_cells<<endl;
 788   cout<<"N_newpoints="<<N_newpoints<<endl;
 789   cout<<"N_newcells="<<N_newcells<<endl;
 790   EG_VTKSP(vtkUnstructuredGrid,grid_tmp);
 791   allocateGrid(grid_tmp,N_cells+N_newcells,N_points+N_newpoints);
 792   m_total_N_newpoints+=N_newpoints; m_total_N_newcells+=N_newcells;
 793
 794   makeCopyNoAlloc(m_grid, grid_tmp);
 795     //initialize new node counter
 796   m_newNodeId=N_points;
 797
 798   remove_FP_actor(grid_tmp);
 799
 800   makeCopy(grid_tmp,m_grid);
 801
 802   cout<<"===remove_FP_all END==="<<endl;
 803   return(0);
 804 }
 805
 806 int SurfaceSmoother::remove_EP_all()
 807 {
 808   cout<<"===remove_EP_all START==="<<endl;
 809
 810   getAllSurfaceCells(m_AllCells,m_grid);
 811   getSurfaceCells(m_bcs, m_SelectedCells, m_grid);
 812   EG_VTKDCC(vtkIntArray, cell_code, m_grid, "cell_code");
 813   EG_VTKDCN(vtkDoubleArray, node_meshdensity, m_grid, "node_meshdensity");
 814   m_SelectedNodes.clear();
 815   getSurfaceNodes(m_bcs,m_SelectedNodes,m_grid);
 816   getNodesFromCells(m_AllCells, nodes, m_grid);
 817   setGrid(m_grid);
 818   setCells(m_AllCells);
 819   cout<<"m_AllCells.size()="<<m_AllCells.size()<<endl;
 820
 821   N_removed_EP=0;
 822
 823   N_points=m_grid->GetNumberOfPoints();
 824   N_cells=m_grid->GetNumberOfCells();
 825   N_newpoints=0;
 826   N_newcells=0;
 827
 828   hitlist.resize(N_points);
 829   offset.resize(N_points);
 830
 831   marked_cells.clear();
 832   marked_nodes.clear();
 833
 834   remove_EP_counter();
 835
 836     //unmark cells (TODO: optimize)
 837   marked_cells.clear();
 838     //init grid_tmp
 839   N_points=m_grid->GetNumberOfPoints();
 840   N_cells=m_grid->GetNumberOfCells();
 841   cout<<"N_points="<<N_points<<endl;
 842   cout<<"N_cells="<<N_cells<<endl;
 843   cout<<"N_newpoints="<<N_newpoints<<endl;
 844   cout<<"N_newcells="<<N_newcells<<endl;
 845   EG_VTKSP(vtkUnstructuredGrid,grid_tmp);
 846   allocateGrid(grid_tmp,N_cells+N_newcells,N_points+N_newpoints);
 847   m_total_N_newpoints+=N_newpoints; m_total_N_newcells+=N_newcells;
 848
 849   makeCopyNoAlloc(m_grid, grid_tmp);
 850     //initialize new node counter
 851   m_newNodeId=N_points;
 852
 853   remove_EP_actor(grid_tmp);
 854
 855   makeCopy(grid_tmp,m_grid);
 856
 857   cout<<"===remove_EP_all END==="<<endl;
 858   return(0);
 859 }
 860
 861 int SurfaceSmoother::FullEdit()
 862 {
 863   cout<<"===FullEdit START==="<<endl;
 864
 865   N_inserted_FP=0;
 866   N_inserted_EP=0;
 867   N_removed_FP=0;
 868   N_removed_EP=0;
 869
 870   N_points=m_grid->GetNumberOfPoints();
 871   N_cells=m_grid->GetNumberOfCells();
 872   N_newpoints=0;
 873   N_newcells=0;
 874
 875   hitlist.resize(N_points);
 876   offset.resize(N_points);
 877
 878   marked_cells.clear();
 879   marked_nodes.clear();
 880
 881   if(insert_FP) insert_FP_counter();
 882   if(insert_EP) insert_EP_counter();
 883   if(remove_FP) remove_FP_counter();
 884   if(remove_EP) remove_EP_counter();
 885
 886   cout<<"================="<<endl;
 887   cout<<"hitlist.size()="<<hitlist.size()<<endl;
 888 //   cout<<"hitlist="<<hitlist<<endl;
 889   cout<<"================="<<endl;
 890
 891     //unmark cells (TODO: optimize)
 892   marked_cells.clear();
 893     //init grid_tmp
 894   N_points=m_grid->GetNumberOfPoints();
 895   N_cells=m_grid->GetNumberOfCells();
 896   cout<<"N_points="<<N_points<<endl;
 897   cout<<"N_cells="<<N_cells<<endl;
 898   cout<<"N_newpoints="<<N_newpoints<<endl;
 899   cout<<"N_newcells="<<N_newcells<<endl;
 900   EG_VTKSP(vtkUnstructuredGrid,grid_tmp);
 901   allocateGrid(grid_tmp,N_cells+N_newcells,N_points+N_newpoints);
 902   m_total_N_newpoints+=N_newpoints; m_total_N_newcells+=N_newcells;
 903
 904   makeCopyNoAlloc(m_grid, grid_tmp);//TODO: This will not work if the size of the grid is reduced!
 905     //initialize new node counter
 906   m_newNodeId=N_points;
 907
 908   if(insert_FP) insert_FP_actor(grid_tmp);
 909   if(insert_EP) insert_EP_actor(grid_tmp);
 910
 911   cout<<"================="<<endl;
 912   cout<<"hitlist.size()="<<hitlist.size()<<endl;
 913 //   cout<<"hitlist="<<hitlist<<endl;
 914   cout<<"================="<<endl;
 915   if(remove_FP) remove_FP_actor(grid_tmp);
 916   if(remove_EP) remove_EP_actor(grid_tmp);
 917
 918   makeCopy(grid_tmp,m_grid);
 919
 920   cout<<"===FullEdit END==="<<endl;
 921   return(0);
 922 }
 923
 924 int SurfaceSmoother::remove_EP_all_2()
 925 {
 926   cout<<"===remove_EP_all_2 START==="<<endl;
 927   getAllSurfaceCells(m_AllCells,m_grid);
 928   getSurfaceCells(m_bcs, m_SelectedCells, m_grid);
 929   EG_VTKDCC(vtkIntArray, cell_code, m_grid, "cell_code");
 930   EG_VTKDCN(vtkDoubleArray, node_meshdensity, m_grid, "node_meshdensity");
 931   m_SelectedNodes.clear();
 932   getSurfaceNodes(m_bcs,m_SelectedNodes,m_grid);
 933   getNodesFromCells(m_AllCells, nodes, m_grid);
 934   setGrid(m_grid);
 935   setCells(m_AllCells);
 936   cout<<"m_AllCells.size()="<<m_AllCells.size()<<endl;
 937
 938   N_removed_EP=0;
 939
 940   N_points=m_grid->GetNumberOfPoints();
 941   N_cells=m_grid->GetNumberOfCells();
 942   N_newpoints=0;
 943   N_newcells=0;
 944
 945   hitlist.clear();
 946   offset.clear();
 947   hitlist.resize(N_points);
 948   offset.resize(N_points);
 949
 950   marked_cells.clear();
 951   marked_nodes.clear();
 952
 953   remove_EP_counter();
 954   cout<<"================="<<endl;
 955   cout<<"hitlist.size()="<<hitlist.size()<<endl;
 956 //   cout<<"hitlist="<<hitlist<<endl;
 957   cout<<"================="<<endl;
 958
 959   int kills=0;
 960   int contracts=0;
 961   for(int i=0;i<hitlist.size();i++)
 962   {
 963     if(hitlist[i]==2){
 964       contracts++;
 965       if(DebugLevel>47) cout<<"Deleting point "<<i<<" currently known as "<<i-kills<<endl;
 966
 967       QString num1;num1.setNum(i);
 968       QString num2;num2.setNum(i-kills);
 969 //       GuiMainWindow::pointer()->QuickSave("pre-deleting_"+num1+"_"+num2+".vtu");
 970
 971       bool DelResult=DeletePoint_2(m_grid,i-kills,N_newpoints,N_newcells);
 972       m_total_N_newpoints+=N_newpoints; m_total_N_newcells+=N_newcells;
 973
 974       if(DelResult)
 975       {
 976         kills++;
 977         if(DebugLevel>47) cout<<"Kill successful"<<endl;
 978       }
 979       else
 980       {
 981         if(DebugLevel>47) cout<<"Kill failed"<<endl;
 982         N_removed_EP--;
 983       }
 984
 985 //       GuiMainWindow::pointer()->QuickSave("post-deleting_"+num1+"_"+num2+".vtu");
 986
 987     }
 988   }
 989   cout<<"Killed: "<<kills<<"/"<<contracts<<endl;
 990   if(kills!=contracts) {cout<<"MISSION FAILED"<<endl;EG_BUG;}
 991   cout<<"===remove_EP_all_2 END==="<<endl;
 992   return(0);
 993 }
 994
 995
 996 //count all to remove, then remove them one by one
 997 int SurfaceSmoother::remove_FP_all_2()
 998 {
 999   cout<<"===remove_FP_all_2 START==="<<endl;
1000 /*  cout<<"+++++++"<<endl;
1001   cout_grid(cout,m_grid,true,true,true,true);
1002   cout<<"+++++++"<<endl;*/
1003
1004   getAllSurfaceCells(m_AllCells,m_grid);
1005   getSurfaceCells(m_bcs, m_SelectedCells, m_grid);
1006   EG_VTKDCC(vtkIntArray, cell_code, m_grid, "cell_code");
1007   EG_VTKDCN(vtkDoubleArray, node_meshdensity, m_grid, "node_meshdensity");
1008   m_SelectedNodes.clear();
1009   getSurfaceNodes(m_bcs,m_SelectedNodes,m_grid);
1010   getNodesFromCells(m_AllCells, nodes, m_grid);
1011   setGrid(m_grid);
1012   setCells(m_AllCells);
1013   cout<<"m_AllCells.size()="<<m_AllCells.size()<<endl;
1014
1015   N_removed_FP=0;
1016
1017   N_points=m_grid->GetNumberOfPoints();
1018   N_cells=m_grid->GetNumberOfCells();
1019   N_newpoints=0;
1020   N_newcells=0;
1021
1022   hitlist.clear();
1023   offset.clear();
1024   hitlist.resize(N_points);
1025   offset.resize(N_points);
1026
1027   marked_cells.clear();
1028   marked_nodes.clear();
1029
1030 //   DualSave("pre-counter");
1031   remove_FP_counter();
1032 //   DualSave("post-counter");
1033
1034 //   cout_grid(cout,m_grid);
1035   cout<<"================="<<endl;
1036   cout<<"hitlist.size()="<<hitlist.size()<<endl;
1037 //   cout<<"hitlist="<<hitlist<<endl;
1038   cout<<"================="<<endl;
1039
1040   int kills=0;
1041   int contracts=0;
1042   for(int i=0;i<hitlist.size();i++)
1043   {
1044     if(hitlist[i]==1){
1045       contracts++;
1046       if(DebugLevel>47) cout<<"Deleting point "<<i<<" currently known as "<<i-kills<<endl;
1047
1048       QString num1;num1.setNum(i);
1049       QString num2;num2.setNum(i-kills);
1050 //       GuiMainWindow::pointer()->QuickSave("pre-deleting_"+num1+"_"+num2+".vtu");
1051
1052       bool DelResult=DeletePoint_2(m_grid,i-kills,N_newpoints,N_newcells);
1053       m_total_N_newpoints+=N_newpoints; m_total_N_newcells+=N_newcells;
1054
1055       if(DelResult)
1056       {
1057         kills++;
1058         if(DebugLevel>47) cout<<"Kill successful"<<endl;
1059       }
1060       else
1061       {
1062         if(DebugLevel>47) cout<<"Kill failed"<<endl;
1063         N_removed_FP--;
1064       }
1065
1066 //       GuiMainWindow::pointer()->QuickSave("post-deleting_"+num1+"_"+num2+".vtu");
1067
1068     }
1069   }
1070   cout<<"Killed: "<<kills<<"/"<<contracts<<endl;
1071   if(kills!=contracts) {cout<<"MISSION FAILED"<<endl;EG_BUG;}
1072   cout<<"===remove_FP_all_2 END==="<<endl;
1073   return(0);
1074 }
1075
1076 //count all to remove, then remove them all at once
1077 int SurfaceSmoother::remove_FP_all_3()
1078 {
1079   cout<<"===remove_FP_all_3 START==="<<endl;
1080
1081   getAllSurfaceCells(m_AllCells,m_grid);
1082   getSurfaceCells(m_bcs, m_SelectedCells, m_grid);
1083   EG_VTKDCC(vtkIntArray, cell_code, m_grid, "cell_code");
1084   EG_VTKDCN(vtkDoubleArray, node_meshdensity, m_grid, "node_meshdensity");
1085   m_SelectedNodes.clear();
1086   getSurfaceNodes(m_bcs,m_SelectedNodes,m_grid);
1087   getNodesFromCells(m_AllCells, nodes, m_grid);
1088   setGrid(m_grid);
1089   setCells(m_AllCells);
1090   cout<<"m_AllCells.size()="<<m_AllCells.size()<<endl;
1091
1092   N_removed_FP=0;
1093
1094   N_points=m_grid->GetNumberOfPoints();
1095   N_cells=m_grid->GetNumberOfCells();
1096   N_newpoints=0;
1097   N_newcells=0;
1098
1099   hitlist.clear();
1100   offset.clear();
1101   hitlist.resize(N_points);
1102   offset.resize(N_points);
1103
1104   marked_cells.clear();
1105   marked_nodes.clear();
1106
1107   remove_FP_counter();
1108   cout<<"================="<<endl;
1109   cout<<"hitlist.size()="<<hitlist.size()<<endl;
1110   cout<<"================="<<endl;
1111
1112   QSet <vtkIdType> DeadNodes;
1113   for(vtkIdType i=0;i<hitlist.size();i++)
1114   {
1115     if(hitlist[i]==1) DeadNodes.insert(i);
1116   }
1117   int N_newpoints=0;
1118   int N_newcells=0;
1119   DeleteSetOfPoints(m_grid, DeadNodes, N_newpoints, N_newcells);
1120   cout<<"N_newpoints="<<N_newpoints<<endl;
1121   cout<<"N_newcells="<<N_newcells<<endl;
1122
1123   int kills=-N_newpoints;
1124   int contracts=DeadNodes.size();
1125   cout<<"Killed: "<<kills<<"/"<<contracts<<endl;
1126   if(kills!=contracts) {cout<<"MISSION FAILED"<<endl;EG_BUG;}
1127   cout<<"===remove_FP_all_3 END==="<<endl;
1128   return(0);
1129 }
1130
1131 //count all to remove, then remove them all at once
1132 int SurfaceSmoother::remove_EP_all_3()
1133 {
1134   cout<<"===remove_EP_all_3 START==="<<endl;
1135
1136   getAllSurfaceCells(m_AllCells,m_grid);
1137   getSurfaceCells(m_bcs, m_SelectedCells, m_grid);
1138   EG_VTKDCC(vtkIntArray, cell_code, m_grid, "cell_code");
1139   EG_VTKDCN(vtkDoubleArray, node_meshdensity, m_grid, "node_meshdensity");
1140   m_SelectedNodes.clear();
1141   getSurfaceNodes(m_bcs,m_SelectedNodes,m_grid);
1142   getNodesFromCells(m_AllCells, nodes, m_grid);
1143   setGrid(m_grid);
1144   setCells(m_AllCells);
1145   cout<<"m_AllCells.size()="<<m_AllCells.size()<<endl;
1146
1147   N_removed_EP=0;
1148
1149   N_points=m_grid->GetNumberOfPoints();
1150   N_cells=m_grid->GetNumberOfCells();
1151   N_newpoints=0;
1152   N_newcells=0;
1153
1154   hitlist.clear();
1155   offset.clear();
1156   hitlist.resize(N_points);
1157   offset.resize(N_points);
1158
1159   marked_cells.clear();
1160   marked_nodes.clear();
1161
1162   remove_EP_counter();
1163   cout<<"================="<<endl;
1164   cout<<"hitlist.size()="<<hitlist.size()<<endl;
1165   cout<<"================="<<endl;
1166
1167   QSet <vtkIdType> DeadNodes;
1168   for(vtkIdType i=0;i<hitlist.size();i++)
1169   {
1170     if(hitlist[i]==1) DeadNodes.insert(i);
1171   }
1172   int N_newpoints=0;
1173   int N_newcells=0;
1174   DeleteSetOfPoints(m_grid, DeadNodes, N_newpoints, N_newcells);
1175   cout<<"N_newpoints="<<N_newpoints<<endl;
1176   cout<<"N_newcells="<<N_newcells<<endl;
1177
1178   int kills=-N_newpoints;
1179   int contracts=DeadNodes.size();
1180   cout<<"Killed: "<<kills<<"/"<<contracts<<endl;
1181   if(kills!=contracts) {cout<<"MISSION FAILED"<<endl;EG_BUG;}
1182   cout<<"===remove_EP_all_3 END==="<<endl;
1183   return(0);
1184 }