createspecialmapping.cpp

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