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