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