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