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