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=1;
  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     EG_VTKDCN(vtkDoubleArray, node_meshdensity, m_grid, "node_meshdensity");
  37
  38     m_SelectedNodes.clear();
  39     getSurfaceNodes(m_bcs,m_SelectedNodes,m_grid);
  40     getNodesFromCells(m_AllCells, nodes, m_grid);
  41     setGrid(m_grid);
  42     setCells(m_AllCells);
  43
  44     cout<<"m_AllCells.size()="<<m_AllCells.size()<<endl;
  45
  46     UpdateNodeType();
  47     EG_VTKDCN(vtkCharArray, node_type, m_grid, "node_type");
  48
  49     //Phase A : Calculate current mesh density
  50     cout<<"===Phase A==="<<endl;
  51
  52     foreach(vtkIdType node,m_SelectedNodes)
  53     {
  54       VertexMeshDensity nodeVMD = getVMD(node,node_type->GetValue(node));
  55       int idx=VMDvector.indexOf(nodeVMD);
  56       if(DebugLevel>3) cout<<"idx="<<idx<<endl;
  57       if(idx!=-1)//specified
  58       {
  59         node_meshdensity->SetValue(node, VMDvector[idx].density);
  60       }
  61       else//unspecified
  62       {
  63         double L=CurrentVertexAvgDist(node,n2n,m_grid);
  64         double D=1./L;
  65         node_meshdensity->SetValue(node, D);
  66       }
  67     }
  68
  69     //Phase B : define desired mesh density
  70     cout<<"===Phase B==="<<endl;
  71     double diff=Convergence_meshdensity+1;
  72     if(DebugLevel>3) cout<<"before loop: diff="<<diff<<endl;
  73     bool first=true;
  74     int iter=0;
  75     int maxiter=100;
  76     do {
  77       if(DebugLevel>2) cout<<"--->diff="<<diff<<endl;
  78       first=true;
  79       foreach(vtkIdType node,m_SelectedNodes)
  80       {
  81         if(DebugLevel>2) cout<<"======>"<<endl;
  82         VertexMeshDensity nodeVMD = getVMD(node,node_type->GetValue(node));
  83         int idx=VMDvector.indexOf(nodeVMD);
  84         if(DebugLevel>2) cout<<"------>idx="<<idx<<endl;
  85         if(idx!=-1)//specified
  86         {
  87           node_meshdensity->SetValue(node, VMDvector[idx].density);
  88         }
  89         else//unspecified
  90         {
  91           double D=DesiredMeshDensity(node,n2n,m_grid);
  92           if(first) {
  93             if(DebugLevel>2) {
  94               cout<<"------>FIRST:"<<endl;
  95               cout<<"------>D="<<D<<endl;
  96               cout<<"------>node_meshdensity->GetValue("<<node<<")="<<node_meshdensity->GetValue(node)<<endl;
  97               cout<<"------>D-node_meshdensity->GetValue("<<node<<")="<<D-node_meshdensity->GetValue(node)<<endl;
  98               cout<<"------>diff=abs(D-node_meshdensity->GetValue("<<node<<"))="<<abs(D-node_meshdensity->GetValue(node))<<endl;
  99             }
 100             diff=abs(D-node_meshdensity->GetValue(node));
 101             first=false;
 102           }
 103           else {
 104             if(DebugLevel>2) {
 105               cout<<"------>NOT FIRST:"<<endl;
 106               cout<<"------>D="<<D<<endl;
 107               cout<<"------>node_meshdensity->GetValue("<<node<<")="<<node_meshdensity->GetValue(node)<<endl;
 108               cout<<"------>D-node_meshdensity->GetValue("<<node<<")="<<D-node_meshdensity->GetValue(node)<<endl;
 109               cout<<"------>diff=abs(D-node_meshdensity->GetValue("<<node<<"))="<<abs(D-node_meshdensity->GetValue(node))<<endl;
 110               cout<<"------>diff="<<diff<<endl;
 111               cout<<"------>max(abs(D-node_meshdensity->GetValue("<<node<<")),diff)="<<max(abs(D-node_meshdensity->GetValue(node)),diff)<<endl;
 112             }
 113             diff=max(abs(D-node_meshdensity->GetValue(node)),diff);
 114           }
 115           node_meshdensity->SetValue(node, D);
 116         }
 117         if(DebugLevel>2) cout<<"======>"<<endl;
 118       }
 119       iter++;
 120     } while(diff>Convergence_meshdensity && !first && iter<maxiter);// if first=true, it means no new mesh density has been defined (all densities specified)
 121     cout<<"iter="<<iter<<endl;
 122     if(iter>=maxiter) cout<<"WARNING: Desired convergence factor has not been reached!"<<endl;
 123
 124     //Phase C: Prepare edge_map
 125     cout<<"===Phase C==="<<endl;
 126     edge_map.clear();
 127     vtkIdType edgeId=1;
 128     foreach(vtkIdType node1,m_SelectedNodes)
 129     {
 130 //       cout<<"node1="<<node1<<endl;
 131       foreach(vtkIdType node2,n2n[node1])
 132       {
 133         if(edge_map[OrderedPair(node1,node2)]==0) { //this edge hasn't been numbered yet
 134           edge_map[OrderedPair(node1,node2)]=edgeId;edgeId++;
 135         }
 136       }
 137     }
 138     cout<<"edge_map.size()="<<edge_map.size()<<endl;
 139
 140     //Phase D: edit points
 141     cout<<"===Phase D==="<<endl;
 142     N_inserted_FP=0;
 143     N_inserted_EP=0;
 144     N_removed_FP=0;
 145     N_removed_EP=0;
 146
 147     //Method 1
 148 //     FullEdit();
 149
 150     //Method 2
 151 /*    if(insert_FP) insert_FP_all();
 152     if(insert_EP) insert_EP_all();
 153     if(remove_FP) remove_FP_all();
 154     if(remove_EP) remove_EP_all();*/
 155
 156     //Method 3
 157     if(insert_FP) {
 158       insert_FP_all();
 159       DualSave("file1");
 160     }
 161
 162     if(insert_EP) {
 163       insert_EP_all();
 164       DualSave("file2");
 165     }
 166
 167     if(remove_FP) {
 168       remove_FP_all_2();
 169       DualSave("file3");
 170     }
 171
 172     if(remove_EP) {
 173       remove_EP_all_2();
 174       DualSave("file4");
 175     }
 176
 177     //Phase E : Delaunay swap
 178     if(DoSwap) {
 179       QSet<int> bcs_complement=complementary_bcs(m_bcs,m_grid,cells);
 180       cout<<"m_bcs="<<m_bcs<<endl;
 181       cout<<"bcs_complement="<<bcs_complement<<endl;
 182
 183       SwapTriangles swap;
 184       swap.setGrid(m_grid);
 185       swap.setBoundaryCodes(bcs_complement);
 186       swap();
 187     }
 188
 189       //Phase F : translate points to smooth grid
 190       //4 possibilities
 191       //vtk smooth 1
 192       //vtk smooth 2
 193       //laplacian smoothing with projection
 194       //Roland smoothing with projection
 195
 196       //laplacian smoothing with projection
 197     if(DoLaplaceSmoothing) {
 198       LaplaceSmoother Lap;
 199       Lap.SetInput(m_bcs,m_grid);
 200       Lap.SetNumberOfIterations(N_SmoothIterations);
 201       Lap();
 202     }
 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, N_newpoints, N_newcells)!=-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, N_newpoints, N_newcells)!=-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
 853   cout<<"===FullEdit END==="<<endl;
 854   return(0);
 855 }
 856
 857 int CreateSpecialMapping::remove_EP_all_2()
 858 {
 859   cout<<"===remove_EP_all_2 START==="<<endl;
 860   getAllSurfaceCells(m_AllCells,m_grid);
 861   getSurfaceCells(m_bcs, m_SelectedCells, m_grid);
 862   EG_VTKDCC(vtkIntArray, cell_code, m_grid, "cell_code");
 863   EG_VTKDCN(vtkDoubleArray, node_meshdensity, m_grid, "node_meshdensity");
 864   m_SelectedNodes.clear();
 865   getSurfaceNodes(m_bcs,m_SelectedNodes,m_grid);
 866   getNodesFromCells(m_AllCells, nodes, m_grid);
 867   setGrid(m_grid);
 868   setCells(m_AllCells);
 869   cout<<"m_AllCells.size()="<<m_AllCells.size()<<endl;
 870
 871   N_removed_EP=0;
 872
 873   N_points=m_grid->GetNumberOfPoints();
 874   N_cells=m_grid->GetNumberOfCells();
 875   N_newpoints=0;
 876   N_newcells=0;
 877
 878   hitlist.clear();
 879   offset.clear();
 880   hitlist.resize(N_points);
 881   offset.resize(N_points);
 882
 883   marked_cells.clear();
 884   marked_nodes.clear();
 885
 886   remove_EP_counter();
 887   cout<<"================="<<endl;
 888   cout<<"hitlist="<<hitlist<<endl;
 889   cout<<"================="<<endl;
 890
 891   int kills=0;
 892   int contracts=0;
 893   for(int i=0;i<hitlist.size();i++)
 894   {
 895     if(hitlist[i]==2){
 896       contracts++;
 897       cout<<"Deleting point "<<i<<" currently known as "<<i-kills<<endl;
 898
 899       QString num1;num1.setNum(i);
 900       QString num2;num2.setNum(i-kills);
 901 //       GuiMainWindow::pointer()->QuickSave("pre-deleting_"+num1+"_"+num2+".vtu");
 902
 903       bool DelResult=DeletePoint_2(m_grid,i-kills,N_newpoints,N_newcells);
 904       m_total_N_newpoints+=N_newpoints; m_total_N_newcells+=N_newcells;
 905
 906       if(DelResult)
 907       {
 908         kills++;
 909         cout<<"Kill successful"<<endl;
 910       }
 911       else
 912       {
 913         cout<<"Kill failed"<<endl;
 914         N_removed_EP--;
 915       }
 916
 917 //       GuiMainWindow::pointer()->QuickSave("post-deleting_"+num1+"_"+num2+".vtu");
 918
 919     }
 920   }
 921   cout<<"Killed: "<<kills<<"/"<<contracts<<endl;
 922   if(kills!=contracts) {cout<<"MISSION FAILED"<<endl;EG_BUG;}
 923   cout<<"===remove_EP_all_2 END==="<<endl;
 924   return(0);
 925 }
 926
 927 int CreateSpecialMapping::remove_FP_all_2()
 928 {
 929   cout<<"===remove_FP_all_2 START==="<<endl;
 930   cout<<"+++++++"<<endl;
 931   cout_grid(cout,m_grid,true,true,true,true);
 932   cout<<"+++++++"<<endl;
 933
 934   getAllSurfaceCells(m_AllCells,m_grid);
 935   getSurfaceCells(m_bcs, m_SelectedCells, m_grid);
 936   EG_VTKDCC(vtkIntArray, cell_code, m_grid, "cell_code");
 937   EG_VTKDCN(vtkDoubleArray, node_meshdensity, m_grid, "node_meshdensity");
 938   m_SelectedNodes.clear();
 939   getSurfaceNodes(m_bcs,m_SelectedNodes,m_grid);
 940   getNodesFromCells(m_AllCells, nodes, m_grid);
 941   setGrid(m_grid);
 942   setCells(m_AllCells);
 943   cout<<"m_AllCells.size()="<<m_AllCells.size()<<endl;
 944
 945   N_removed_FP=0;
 946
 947   N_points=m_grid->GetNumberOfPoints();
 948   N_cells=m_grid->GetNumberOfCells();
 949   N_newpoints=0;
 950   N_newcells=0;
 951
 952   hitlist.clear();
 953   offset.clear();
 954   hitlist.resize(N_points);
 955   offset.resize(N_points);
 956
 957   marked_cells.clear();
 958   marked_nodes.clear();
 959
 960   DualSave("pre-counter");
 961   remove_FP_counter();
 962   DualSave("post-counter");
 963
 964   cout_grid(cout,m_grid);
 965   cout<<"================="<<endl;
 966   cout<<"hitlist="<<hitlist<<endl;
 967   cout<<"================="<<endl;
 968
 969   int kills=0;
 970   int contracts=0;
 971   for(int i=0;i<hitlist.size();i++)
 972   {
 973     if(hitlist[i]==1){
 974       contracts++;
 975       cout<<"Deleting point "<<i<<" currently known as "<<i-kills<<endl;
 976
 977       QString num1;num1.setNum(i);
 978       QString num2;num2.setNum(i-kills);
 979 //       GuiMainWindow::pointer()->QuickSave("pre-deleting_"+num1+"_"+num2+".vtu");
 980
 981       bool DelResult=DeletePoint_2(m_grid,i-kills,N_newpoints,N_newcells);
 982       m_total_N_newpoints+=N_newpoints; m_total_N_newcells+=N_newcells;
 983
 984       if(DelResult)
 985       {
 986         kills++;
 987         cout<<"Kill successful"<<endl;
 988       }
 989       else
 990       {
 991         cout<<"Kill failed"<<endl;
 992         N_removed_FP--;
 993       }
 994
 995 //       GuiMainWindow::pointer()->QuickSave("post-deleting_"+num1+"_"+num2+".vtu");
 996
 997     }
 998   }
 999   cout<<"Killed: "<<kills<<"/"<<contracts<<endl;
1000   if(kills!=contracts) {cout<<"MISSION FAILED"<<endl;EG_BUG;}
1001   cout<<"===remove_FP_all_2 END==="<<endl;
1002   return(0);
1003 }