createspecialmapping.cpp

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