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