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