createspecialmapping.cpp

   1 #include "createspecialmapping.h"
   2
   3 #include <QString>
   4 #include <QTextStream>
   5 #include <vtkCharArray.h>
   6
   7 #include "smoothingutilities.h"
   8
   9 #include "swaptriangles.h"
  10 #include "laplacesmoother.h"
  11 #include "guimainwindow.h"
  12
  13 #include <iostream>
  14 using namespace std;
  15
  16 CreateSpecialMapping::CreateSpecialMapping()
  17 {
  18   DebugLevel=1;
  19 }
  20
  21 int CreateSpecialMapping::Process()
  22 {
  23   int i_iter=0;
  24   for(i_iter=0;i_iter<NumberOfIterations;i_iter++)//TODO:Optimize this loop
  25   {
  26     cout<<"===ITERATION NB "<<i_iter<<"/"<<NumberOfIterations<<"==="<<endl;
  27
  28     m_total_N_newpoints=0;
  29     m_total_N_newcells=0;
  30
  31     getAllSurfaceCells(m_AllCells,m_grid);
  32     getSurfaceCells(m_bcs, m_SelectedCells, m_grid);
  33     cout<<"m_AllCells.size()="<<m_AllCells.size()<<endl;
  34
  35     EG_VTKDCC(vtkIntArray, cell_code, m_grid, "cell_code");
  36     EG_VTKDCN(vtkDoubleArray, node_meshdensity, m_grid, "node_meshdensity");
  37
  38     m_SelectedNodes.clear();
  39     getSurfaceNodes(m_bcs,m_SelectedNodes,m_grid);
  40     getNodesFromCells(m_AllCells, nodes, m_grid);
  41     setGrid(m_grid);
  42     setCells(m_AllCells);
  43
  44     cout<<"m_AllCells.size()="<<m_AllCells.size()<<endl;
  45
  46     UpdateNodeType();
  47     EG_VTKDCN(vtkCharArray, node_type, m_grid, "node_type");
  48
  49     //Phase A : Calculate current mesh density
  50     cout<<"===Phase A==="<<endl;
  51
  52     foreach(vtkIdType node,m_SelectedNodes)
  53     {
  54       VertexMeshDensity nodeVMD = getVMD(node,node_type->GetValue(node));
  55       int idx=VMDvector.indexOf(nodeVMD);
  56       if(DebugLevel>3) cout<<"idx="<<idx<<endl;
  57       if(idx!=-1)//specified
  58       {
  59         node_meshdensity->SetValue(node, VMDvector[idx].density);
  60       }
  61       else//unspecified
  62       {
  63         double L=CurrentVertexAvgDist(node,n2n,m_grid);
  64         double D=1./L;
  65         node_meshdensity->SetValue(node, D);
  66       }
  67     }
  68
  69     //Phase B : define desired mesh density
  70     cout<<"===Phase B==="<<endl;
  71     double diff=Convergence_meshdensity+1;
  72     if(DebugLevel>3) cout<<"before loop: diff="<<diff<<endl;
  73     bool first=true;
  74     int iter=0;
  75     int maxiter=100;
  76     do {
  77       if(DebugLevel>2) cout<<"--->diff="<<diff<<endl;
  78       first=true;
  79       foreach(vtkIdType node,m_SelectedNodes)
  80       {
  81         if(DebugLevel>2) cout<<"======>"<<endl;
  82         VertexMeshDensity nodeVMD = getVMD(node,node_type->GetValue(node));
  83         int idx=VMDvector.indexOf(nodeVMD);
  84         if(DebugLevel>2) cout<<"------>idx="<<idx<<endl;
  85         if(idx!=-1)//specified
  86         {
  87           node_meshdensity->SetValue(node, VMDvector[idx].density);
  88         }
  89         else//unspecified
  90         {
  91           double D=DesiredMeshDensity(node,n2n,m_grid);
  92           if(first) {
  93             if(DebugLevel>2) {
  94               cout<<"------>FIRST:"<<endl;
  95               cout<<"------>D="<<D<<endl;
  96               cout<<"------>node_meshdensity->GetValue("<<node<<")="<<node_meshdensity->GetValue(node)<<endl;
  97               cout<<"------>D-node_meshdensity->GetValue("<<node<<")="<<D-node_meshdensity->GetValue(node)<<endl;
  98               cout<<"------>diff=abs(D-node_meshdensity->GetValue("<<node<<"))="<<abs(D-node_meshdensity->GetValue(node))<<endl;
  99             }
 100             diff=abs(D-node_meshdensity->GetValue(node));
 101             first=false;
 102           }
 103           else {
 104             if(DebugLevel>2) {
 105               cout<<"------>NOT FIRST:"<<endl;
 106               cout<<"------>D="<<D<<endl;
 107               cout<<"------>node_meshdensity->GetValue("<<node<<")="<<node_meshdensity->GetValue(node)<<endl;
 108               cout<<"------>D-node_meshdensity->GetValue("<<node<<")="<<D-node_meshdensity->GetValue(node)<<endl;
 109               cout<<"------>diff=abs(D-node_meshdensity->GetValue("<<node<<"))="<<abs(D-node_meshdensity->GetValue(node))<<endl;
 110               cout<<"------>diff="<<diff<<endl;
 111               cout<<"------>max(abs(D-node_meshdensity->GetValue("<<node<<")),diff)="<<max(abs(D-node_meshdensity->GetValue(node)),diff)<<endl;
 112             }
 113             diff=max(abs(D-node_meshdensity->GetValue(node)),diff);
 114           }
 115           node_meshdensity->SetValue(node, D);
 116         }
 117         if(DebugLevel>2) cout<<"======>"<<endl;
 118       }
 119       iter++;
 120     } while(diff>Convergence_meshdensity && !first && iter<maxiter);// if first=true, it means no new mesh density has been defined (all densities specified)
 121     cout<<"iter="<<iter<<endl;
 122     if(iter>=maxiter) cout<<"WARNING: Desired convergence factor has not been reached!"<<endl;
 123
 124     //Phase C: Prepare edge_map
 125     cout<<"===Phase C==="<<endl;
 126     edge_map.clear();
 127     vtkIdType edgeId=1;
 128     foreach(vtkIdType node1,m_SelectedNodes)
 129     {
 130 //       cout<<"node1="<<node1<<endl;
 131       foreach(vtkIdType node2,n2n[node1])
 132       {
 133         if(edge_map[OrderedPair(node1,node2)]==0) { //this edge hasn't been numbered yet
 134           edge_map[OrderedPair(node1,node2)]=edgeId;edgeId++;
 135         }
 136       }
 137     }
 138     cout<<"edge_map.size()="<<edge_map.size()<<endl;
 139
 140     //Phase D: edit points
 141     cout<<"===Phase D==="<<endl;
 142     N_inserted_FP=0;
 143     N_inserted_EP=0;
 144     N_removed_FP=0;
 145     N_removed_EP=0;
 146
 147     //Method 1
 148 //     FullEdit();
 149
 150     //Method 2
 151 /*    if(insert_FP) insert_FP_all();
 152     if(insert_EP) insert_EP_all();
 153     if(remove_FP) remove_FP_all();
 154     if(remove_EP) remove_EP_all();*/
 155
 156     ofstream file1,file2,file3,file4;
 157
 158     //Method 3
 159     if(insert_FP) {
 160       insert_FP_all();
 161 /*      file1.open ("file1.txt");
 162       cout_grid(file1,m_grid,true,true,true,true);
 163       file1.close();*/
 164     }
 165
 166     if(insert_EP) {
 167       insert_EP_all();
 168 /*      file2.open ("file2.txt");
 169       cout_grid(file2,m_grid,true,true,true,true);
 170       file2.close();*/
 171     }
 172
 173     if(remove_FP) {
 174       remove_FP_all_2();
 175 /*      file3.open ("file3.txt");
 176       cout_grid(file3,m_grid,true,true,true,true);
 177       file3.close();*/
 178     }
 179
 180     if(remove_EP) {
 181       remove_EP_all_2();
 182 /*      file4.open ("file4.txt");
 183       cout_grid(file4,m_grid,true,true,true,true);
 184       file4.close();*/
 185     }
 186
 187     //Phase E : Delaunay swap
 188     if(DoSwap) {
 189       QSet<int> bcs_complement=complementary_bcs(m_bcs,m_grid,cells);
 190       cout<<"m_bcs="<<m_bcs<<endl;
 191       cout<<"bcs_complement="<<bcs_complement<<endl;
 192
 193       SwapTriangles swap;
 194       swap.setGrid(m_grid);
 195       swap.setBoundaryCodes(bcs_complement);
 196       swap();
 197     }
 198
 199       //Phase F : translate points to smooth grid
 200       //4 possibilities
 201       //vtk smooth 1
 202       //vtk smooth 2
 203       //laplacian smoothing with projection
 204       //Roland smoothing with projection
 205
 206       //laplacian smoothing with projection
 207     if(DoLaplaceSmoothing) {
 208       LaplaceSmoother Lap;
 209       Lap.SetInput(m_bcs,m_grid);
 210       Lap.SetNumberOfIterations(N_SmoothIterations);
 211       Lap();
 212     }
 213
 214     cout<<"===Summary==="<<endl;
 215     cout<<"N_inserted_FP="<<N_inserted_FP<<endl;
 216     cout<<"N_inserted_EP="<<N_inserted_EP<<endl;
 217     cout<<"N_removed_FP="<<N_removed_FP<<endl;
 218     cout<<"N_removed_EP="<<N_removed_EP<<endl;
 219
 220     cout<<"N_points="<<N_points<<endl;
 221     cout<<"N_cells="<<N_cells<<endl;
 222     cout<<"m_total_N_newpoints="<<m_total_N_newpoints<<endl;
 223     cout<<"m_total_N_newcells="<<m_total_N_newcells<<endl;
 224     cout<<"============"<<endl;
 225
 226     if(m_total_N_newpoints==0 && m_total_N_newcells==0) break;
 227
 228   }
 229
 230   cout<<"i_iter/NumberOfIterations="<<i_iter<<"/"<<NumberOfIterations<<endl;
 231   UpdateMeshDensity();
 232   return 1;
 233 }
 234 //end of process
 235
 236 VertexMeshDensity CreateSpecialMapping::getVMD(vtkIdType node, char VertexType)
 237 {
 238   VertexMeshDensity VMD;
 239   VMD.type=VertexType;
 240   VMD.density=0;
 241   VMD.CurrentNode=node;
 242   EG_VTKDCC(vtkIntArray, cell_code, m_grid, "cell_code");
 243 /*  createNodeMapping(nodes, _nodes, m_grid);
 244   createNodeToCell(m_AllCells, nodes, _nodes, n2c, m_grid);*/
 245
 246   QSet <int> bc;
 247   foreach(vtkIdType C, n2c[node])
 248   {
 249     bc.insert(cell_code->GetValue(C));
 250   }
 251   VMD.BClist.resize(bc.size());
 252   qCopy(bc.begin(),bc.end(),VMD.BClist.begin());
 253   qSort(VMD.BClist.begin(),VMD.BClist.end());
 254   return(VMD);
 255 }
 256
 257 int CreateSpecialMapping::insert_FP_counter()
 258 {
 259   cout<<"===insert_FP_counter() START==="<<endl;
 260   foreach(vtkIdType id_cell, m_SelectedCells)
 261   {
 262     if( !marked_cells[id_cell] && insert_fieldpoint(id_cell) )
 263     {
 264       cout<<"inserting a field point "<<id_cell<<endl;
 265       N_inserted_FP++;
 266       marked_cells[id_cell]=true;
 267       N_newcells+=2;
 268       N_newpoints+=1;
 269     }
 270   }
 271   cout<<"===insert_FP_counter() END==="<<endl;
 272   return(0);
 273 }
 274
 275 int CreateSpecialMapping::insert_EP_counter()
 276 {
 277   cout<<"===insert_EP_counter() START==="<<endl;
 278   StencilVector.clear();
 279   QMapIterator< pair<vtkIdType,vtkIdType>, vtkIdType> edge_map_iter(edge_map);
 280       //rewind the iterator
 281   edge_map_iter.toFront ();
 282       //start loop
 283   while (edge_map_iter.hasNext()) {
 284     edge_map_iter.next();
 285     vtkIdType node1=edge_map_iter.key().first;
 286     vtkIdType node2=edge_map_iter.key().second;
 287     if(DebugLevel>10) cout << "--->(" << node1 << "," << node2 << ")" << ": " << edge_map_iter.value() << endl;
 288     QSet <int> stencil_cells_set;
 289     QVector <int> stencil_cells_vector;
 290     stencil_cells_set=n2c[node1];
 291     stencil_cells_set.intersect(n2c[node2]);
 292     if(DebugLevel>10) cout<<"stencil_cells_set="<<stencil_cells_set<<endl;
 293
 294     stencil_cells_vector.resize(stencil_cells_set.size());
 295     qCopy(stencil_cells_set.begin(),stencil_cells_set.end(),stencil_cells_vector.begin());
 296     if(DebugLevel>10) cout<<"stencil_cells_vector="<<stencil_cells_vector<<endl;
 297
 298     vtkIdType id_cell=stencil_cells_vector[0];
 299     int SideToSplit = getSide(id_cell,m_grid,node1,node2);
 300     if(DebugLevel>10) cout<<"SideToSplit="<<SideToSplit<<endl;
 301     if(DebugLevel>10) cout<<"c2c[id_cell][SideToSplit]="<<c2c[id_cell][SideToSplit]<<endl;
 302     if(DebugLevel>10) for(int i=0;i<3;i++) cout<<"c2c[id_cell]["<<i<<"]="<<c2c[id_cell][i]<<endl;
 303     stencil_t S=getStencil(id_cell,SideToSplit);
 304
 305     bool stencil_marked=false;
 306     foreach(vtkIdType C,stencil_cells_vector)
 307     {
 308       if(marked_cells[C]) stencil_marked=true;
 309     }
 310     if(DebugLevel>10) cout<<"stencil_marked="<<stencil_marked<<endl;
 311     if(DebugLevel>10) cout<<"insert_edgepoint(node1,node2)="<<insert_edgepoint(node1,node2)<<endl;
 312
 313     if( !stencil_marked && insert_edgepoint(node1,node2) )
 314     {
 315       if(DebugLevel>1) cout<<"inserting an edge point "<< "(" << node1 << "," << node2 << ")" << ": " << edge_map_iter.value() << endl;
 316       N_inserted_EP++;
 317       foreach(vtkIdType C,stencil_cells_vector) marked_cells[C]=true;
 318       StencilVector.push_back(S);
 319
 320       if(stencil_cells_vector.size()==2)//2 cells around the edge
 321       {
 322         N_newcells+=2;
 323         N_newpoints+=1;
 324       }
 325       else//1 cell around the edge
 326       {
 327         N_newcells+=1;
 328         N_newpoints+=1;
 329       }
 330     }
 331     if(DebugLevel>10) cout <<"--->end of edge processing"<<endl;
 332   }
 333   cout<<"===insert_EP_counter() END==="<<endl;
 334   return(0);
 335 }
 336
 337 int CreateSpecialMapping::remove_FP_counter()
 338 {
 339   cout<<"===remove_FP_counter() START==="<<endl;
 340   UpdateNodeType();
 341   EG_VTKDCN(vtkCharArray, node_type, m_grid, "node_type");
 342   foreach(vtkIdType node,m_SelectedNodes)
 343   {
 344     if(node_type->GetValue(node)==VTK_SIMPLE_VERTEX)
 345     {
 346       bool marked=false;
 347       foreach(vtkIdType C,n2c[node])
 348       {
 349         if(marked_cells[C]) marked=true;
 350       }
 351
 352       QSet <vtkIdType> DeadCells;
 353       QSet <vtkIdType> MutatedCells;
 354       QSet <vtkIdType> MutilatedCells;
 355       if( !marked && remove_fieldpoint(node) && FindSnapPoint(m_grid,node,DeadCells,MutatedCells,MutilatedCells)!=-1)
 356       {
 357         if(DebugLevel>1) cout<<"removing field point "<<node<<endl;
 358         N_removed_FP++;
 359         hitlist[node]=1;
 360         foreach(vtkIdType C,n2c[node]) marked_cells[C]=true;
 361         N_newcells-=2;
 362         N_newpoints-=1;
 363       }
 364     }
 365   }
 366   cout<<"===remove_FP_counter() END==="<<endl;
 367   return(0);
 368 }
 369
 370 int CreateSpecialMapping::remove_EP_counter()
 371 {
 372   cout<<"===remove_EP_counter() START==="<<endl;
 373   UpdateNodeType();
 374   EG_VTKDCN(vtkCharArray, node_type, m_grid, "node_type");
 375   foreach(vtkIdType node,m_SelectedNodes)
 376   {
 377     if(node_type->GetValue(node)==VTK_BOUNDARY_EDGE_VERTEX)
 378     {
 379       bool marked=false;
 380       foreach(vtkIdType C,n2c[node])
 381       {
 382         if(marked_cells[C]) marked=true;
 383       }
 384       QSet <vtkIdType> DeadCells;
 385       QSet <vtkIdType> MutatedCells;
 386       QSet <vtkIdType> MutilatedCells;
 387       if( !marked && remove_edgepoint(node) && FindSnapPoint(m_grid,node,DeadCells,MutatedCells,MutilatedCells)!=-1)
 388       {
 389         cout<<"removing edge point "<<node<<endl;
 390         N_removed_EP++;
 391         hitlist[node]=2;
 392         foreach(vtkIdType C,n2c[node]) marked_cells[C]=true;
 393         if(n2n[node].size()==4)//4 cells around the edge
 394         {
 395           N_newcells-=2;
 396           N_newpoints-=1;
 397         }
 398         else//2 cells around the edge
 399         {
 400           N_newcells-=1;
 401           N_newpoints-=1;
 402         }
 403       }
 404     }
 405   }
 406   cout<<"===remove_EP_counter() END==="<<endl;
 407   return(0);
 408 }
 409
 410 int CreateSpecialMapping::insert_FP_actor(vtkUnstructuredGrid* grid_tmp)
 411 {
 412   cout<<"===insert_FP_actor START==="<<endl;
 413
 414   EG_VTKDCC(vtkIntArray, cell_code_tmp, grid_tmp, "cell_code");
 415   foreach(vtkIdType id_cell, m_SelectedCells)
 416   {
 417 /*    if(marked_cells[id_cell]) cout<<"--->marked_cells["<<id_cell<<"]=TRUE"<<endl;
 418     else cout<<"--->marked_cells["<<id_cell<<"]=FALSE"<<endl;*/
 419
 420     if( !marked_cells[id_cell] && insert_fieldpoint(id_cell) )
 421     {
 422       cout<<"inserting a field point "<<id_cell<<endl;
 423       vtkIdType newBC=cell_code_tmp->GetValue(id_cell);
 424       cout<<"id_cell="<<id_cell<<" newBC="<<newBC<<endl;
 425
 426       vtkIdType N_pts, *pts;
 427       m_grid->GetCellPoints(id_cell, N_pts, pts);
 428       vec3_t C(0,0,0);
 429
 430       int N_neighbours=N_pts;
 431       cout<<"N_neighbours="<<N_neighbours<<endl;
 432       vec3_t corner[4];
 433       vtkIdType pts_triangle[4][3];
 434       for(int i=0;i<N_neighbours;i++)
 435       {
 436         m_grid->GetPoints()->GetPoint(pts[i], corner[i].data());
 437         C+=corner[i];
 438       }
 439       C=(1/(double)N_neighbours)*C;
 440       addPoint(grid_tmp,m_newNodeId,C.data());
 441       vtkIdType intmidpoint=m_newNodeId;
 442       m_newNodeId++;
 443
 444       for(int i=0;i<N_neighbours;i++)
 445       {
 446         pts_triangle[i][0]=pts[i];
 447         pts_triangle[i][1]=pts[(i+1)%N_neighbours];
 448         pts_triangle[i][2]=intmidpoint;
 449         if(i==0)
 450         {
 451           grid_tmp->ReplaceCell(id_cell , 3, pts_triangle[0]);
 452           cell_code_tmp->SetValue(id_cell, newBC);
 453         }
 454         else
 455         {
 456           vtkIdType newCellId = grid_tmp->InsertNextCell(VTK_TRIANGLE,3,pts_triangle[i]);
 457           cell_code_tmp->SetValue(newCellId, newBC);
 458         }
 459       }
 460
 461     }
 462   }
 463   cout<<"===insert_FP_actor END==="<<endl;
 464   return(0);
 465 }
 466
 467 int CreateSpecialMapping::insert_EP_actor(vtkUnstructuredGrid* grid_tmp)
 468 {
 469   cout<<"===insert_EP_actor START==="<<endl;
 470
 471   EG_VTKDCC(vtkIntArray, cell_code_tmp, grid_tmp, "cell_code");
 472   foreach(stencil_t S,StencilVector)
 473   {
 474     if(DebugLevel>10) cout<<"S="<<S<<endl;
 475     vec3_t A,B;
 476     grid_tmp->GetPoint(S.p[1],A.data());
 477     grid_tmp->GetPoint(S.p[3],B.data());
 478     vec3_t M=0.5*(A+B);
 479     addPoint(grid_tmp,m_newNodeId,M.data());
 480
 481     vtkIdType pts_triangle[4][3];
 482
 483     if(S.valid){//there is a neighbour cell
 484       if(DebugLevel>10) cout<<"marked_cells["<<S.id_cell1<<"]=true;"<<endl;
 485       if(DebugLevel>10) cout<<"marked_cells["<<S.id_cell2<<"]=true;"<<endl;
 486       marked_cells[S.id_cell1]=true;
 487       marked_cells[S.id_cell2]=true;
 488
 489       for(int i=0;i<4;i++)
 490       {
 491         pts_triangle[i][0]=S.p[i];
 492         pts_triangle[i][1]=S.p[(i+1)%4];
 493         pts_triangle[i][2]=m_newNodeId;
 494       }
 495
 496       int bc1=cell_code_tmp->GetValue(S.id_cell1);
 497       int bc2=cell_code_tmp->GetValue(S.id_cell2);
 498
 499       grid_tmp->ReplaceCell(S.id_cell1 , 3, pts_triangle[0]);
 500       cell_code_tmp->SetValue(S.id_cell1, bc1);
 501
 502       grid_tmp->ReplaceCell(S.id_cell2 , 3, pts_triangle[1]);
 503       cell_code_tmp->SetValue(S.id_cell2, bc2);
 504
 505       vtkIdType newCellId;
 506       newCellId = grid_tmp->InsertNextCell(VTK_TRIANGLE,3,pts_triangle[2]);
 507       cell_code_tmp->SetValue(newCellId, bc2);
 508       newCellId = grid_tmp->InsertNextCell(VTK_TRIANGLE,3,pts_triangle[3]);
 509       cell_code_tmp->SetValue(newCellId, bc1);
 510     }
 511     else{//there is no neighbour cell
 512       if(DebugLevel>10) cout<<"marked_cells["<<S.id_cell1<<"]=true;"<<endl;
 513       marked_cells[S.id_cell1]=true;
 514
 515       pts_triangle[0][0]=S.p[0];
 516       pts_triangle[0][1]=S.p[1];
 517       pts_triangle[0][2]=m_newNodeId;
 518       pts_triangle[3][0]=S.p[3];
 519       pts_triangle[3][1]=S.p[0];
 520       pts_triangle[3][2]=m_newNodeId;
 521
 522       int bc1=cell_code_tmp->GetValue(S.id_cell1);
 523
 524       grid_tmp->ReplaceCell(S.id_cell1 , 3, pts_triangle[0]);
 525       cell_code_tmp->SetValue(S.id_cell1, bc1);
 526
 527       vtkIdType newCellId;
 528       newCellId = grid_tmp->InsertNextCell(VTK_TRIANGLE,3,pts_triangle[3]);
 529       cell_code_tmp->SetValue(newCellId, bc1);
 530     }
 531
 532     m_newNodeId++;
 533   }
 534   cout<<"===insert_EP_actor END==="<<endl;
 535   return(0);
 536 }
 537
 538 int CreateSpecialMapping::remove_FP_actor(vtkUnstructuredGrid* grid_tmp)
 539 {
 540   cout<<"===remove_FP_actor START==="<<endl;
 541
 542   foreach(vtkIdType node,m_SelectedNodes)
 543   {
 544     if(hitlist[node]==1)
 545     {
 546
 547     }
 548     bool marked=false;
 549     foreach(vtkIdType C,n2c[node])
 550     {
 551       if(marked_cells[C]) marked=true;
 552     }
 553     if( !marked && remove_fieldpoint(node) )
 554     {
 555       if(DebugLevel>1) cout<<"removing field point "<<node<<endl;
 556       foreach(vtkIdType C,n2c[node]) marked_cells[C]=true;
 557       //TODO: Special copy function, leaving out nodes to remove
 558     }
 559   }
 560   cout<<"===remove_FP_actor END==="<<endl;
 561   return(0);
 562 }
 563
 564 int CreateSpecialMapping::remove_EP_actor(vtkUnstructuredGrid* grid_tmp)
 565 {
 566   cout<<"===remove_EP_actor START==="<<endl;
 567
 568   foreach(vtkIdType node,m_SelectedNodes)
 569   {
 570     bool marked=false;
 571     foreach(vtkIdType C,n2c[node])
 572     {
 573       if(marked_cells[C]) marked=true;
 574     }
 575     if( !marked && remove_edgepoint(node) )
 576     {
 577       cout<<"removing edge point "<<node<<endl;
 578       foreach(vtkIdType C,n2c[node]) marked_cells[C]=true;
 579       if(n2n[node].size()==4)//4 cells around the edge
 580       {
 581
 582       }
 583       else//2 cells around the edge
 584       {
 585
 586       }
 587     }
 588   }
 589   cout<<"===remove_EP_actor END==="<<endl;
 590   return(0);
 591 }
 592
 593 int CreateSpecialMapping::insert_FP_all()
 594 {
 595   cout<<"===insert_FP_all START==="<<endl;
 596
 597   getAllSurfaceCells(m_AllCells,m_grid);
 598   getSurfaceCells(m_bcs, m_SelectedCells, m_grid);
 599   EG_VTKDCC(vtkIntArray, cell_code, m_grid, "cell_code");
 600   EG_VTKDCN(vtkDoubleArray, node_meshdensity, m_grid, "node_meshdensity");
 601   m_SelectedNodes.clear();
 602   getSurfaceNodes(m_bcs,m_SelectedNodes,m_grid);
 603   getNodesFromCells(m_AllCells, nodes, m_grid);
 604   setGrid(m_grid);
 605   setCells(m_AllCells);
 606   cout<<"m_AllCells.size()="<<m_AllCells.size()<<endl;
 607
 608   N_inserted_FP=0;
 609
 610   N_points=m_grid->GetNumberOfPoints();
 611   N_cells=m_grid->GetNumberOfCells();
 612   N_newpoints=0;
 613   N_newcells=0;
 614
 615   marked_cells.clear();
 616   marked_nodes.clear();
 617
 618   insert_FP_counter();
 619
 620     //unmark cells (TODO: optimize)
 621   marked_cells.clear();
 622     //init grid_tmp
 623   N_points=m_grid->GetNumberOfPoints();
 624   N_cells=m_grid->GetNumberOfCells();
 625   cout<<"N_points="<<N_points<<endl;
 626   cout<<"N_cells="<<N_cells<<endl;
 627   cout<<"N_newpoints="<<N_newpoints<<endl;
 628   cout<<"N_newcells="<<N_newcells<<endl;
 629   EG_VTKSP(vtkUnstructuredGrid,grid_tmp);
 630   allocateGrid(grid_tmp,N_cells+N_newcells,N_points+N_newpoints);
 631   m_total_N_newpoints+=N_newpoints; m_total_N_newcells+=N_newcells;
 632
 633   makeCopyNoAlloc(m_grid, grid_tmp);
 634     //initialize new node counter
 635   m_newNodeId=N_points;
 636
 637   insert_FP_actor(grid_tmp);
 638
 639   makeCopy(grid_tmp,m_grid);
 640   cout<<"===insert_FP_all END==="<<endl;
 641   return(0);
 642 }
 643
 644 int CreateSpecialMapping::insert_EP_all()
 645 {
 646   cout<<"===insert_EP_all START==="<<endl;
 647
 648   getAllSurfaceCells(m_AllCells,m_grid);
 649   getSurfaceCells(m_bcs, m_SelectedCells, m_grid);
 650   EG_VTKDCC(vtkIntArray, cell_code, m_grid, "cell_code");
 651   EG_VTKDCN(vtkDoubleArray, node_meshdensity, m_grid, "node_meshdensity");
 652   m_SelectedNodes.clear();
 653   getSurfaceNodes(m_bcs,m_SelectedNodes,m_grid);
 654   getNodesFromCells(m_AllCells, nodes, m_grid);
 655   setGrid(m_grid);
 656   setCells(m_AllCells);
 657   cout<<"m_AllCells.size()="<<m_AllCells.size()<<endl;
 658
 659   N_inserted_EP=0;
 660
 661   N_points=m_grid->GetNumberOfPoints();
 662   N_cells=m_grid->GetNumberOfCells();
 663   N_newpoints=0;
 664   N_newcells=0;
 665
 666   marked_cells.clear();
 667   marked_nodes.clear();
 668
 669   insert_EP_counter();
 670
 671     //unmark cells (TODO: optimize)
 672   marked_cells.clear();
 673     //init grid_tmp
 674   N_points=m_grid->GetNumberOfPoints();
 675   N_cells=m_grid->GetNumberOfCells();
 676   cout<<"N_points="<<N_points<<endl;
 677   cout<<"N_cells="<<N_cells<<endl;
 678   cout<<"N_newpoints="<<N_newpoints<<endl;
 679   cout<<"N_newcells="<<N_newcells<<endl;
 680   EG_VTKSP(vtkUnstructuredGrid,grid_tmp);
 681   allocateGrid(grid_tmp,N_cells+N_newcells,N_points+N_newpoints);
 682   m_total_N_newpoints+=N_newpoints; m_total_N_newcells+=N_newcells;
 683
 684   makeCopyNoAlloc(m_grid, grid_tmp);
 685     //initialize new node counter
 686   m_newNodeId=N_points;
 687
 688   insert_EP_actor(grid_tmp);
 689
 690   makeCopy(grid_tmp,m_grid);
 691
 692   cout<<"===insert_EP_all END==="<<endl;
 693   return(0);
 694 }
 695
 696 int CreateSpecialMapping::remove_FP_all()
 697 {
 698   cout<<"===remove_FP_all START==="<<endl;
 699
 700   getAllSurfaceCells(m_AllCells,m_grid);
 701   getSurfaceCells(m_bcs, m_SelectedCells, m_grid);
 702   EG_VTKDCC(vtkIntArray, cell_code, m_grid, "cell_code");
 703   EG_VTKDCN(vtkDoubleArray, node_meshdensity, m_grid, "node_meshdensity");
 704   m_SelectedNodes.clear();
 705   getSurfaceNodes(m_bcs,m_SelectedNodes,m_grid);
 706   getNodesFromCells(m_AllCells, nodes, m_grid);
 707   setGrid(m_grid);
 708   setCells(m_AllCells);
 709   cout<<"m_AllCells.size()="<<m_AllCells.size()<<endl;
 710
 711   N_removed_FP=0;
 712
 713   N_points=m_grid->GetNumberOfPoints();
 714   N_cells=m_grid->GetNumberOfCells();
 715   N_newpoints=0;
 716   N_newcells=0;
 717
 718   hitlist.resize(N_points);
 719   offset.resize(N_points);
 720
 721   marked_cells.clear();
 722   marked_nodes.clear();
 723
 724   remove_FP_counter();
 725
 726     //unmark cells (TODO: optimize)
 727   marked_cells.clear();
 728     //init grid_tmp
 729   N_points=m_grid->GetNumberOfPoints();
 730   N_cells=m_grid->GetNumberOfCells();
 731   cout<<"N_points="<<N_points<<endl;
 732   cout<<"N_cells="<<N_cells<<endl;
 733   cout<<"N_newpoints="<<N_newpoints<<endl;
 734   cout<<"N_newcells="<<N_newcells<<endl;
 735   EG_VTKSP(vtkUnstructuredGrid,grid_tmp);
 736   allocateGrid(grid_tmp,N_cells+N_newcells,N_points+N_newpoints);
 737   m_total_N_newpoints+=N_newpoints; m_total_N_newcells+=N_newcells;
 738
 739   makeCopyNoAlloc(m_grid, grid_tmp);
 740     //initialize new node counter
 741   m_newNodeId=N_points;
 742
 743   remove_FP_actor(grid_tmp);
 744
 745   makeCopy(grid_tmp,m_grid);
 746
 747   cout<<"===remove_FP_all END==="<<endl;
 748   return(0);
 749 }
 750
 751 int CreateSpecialMapping::remove_EP_all()
 752 {
 753   cout<<"===remove_EP_all START==="<<endl;
 754
 755   getAllSurfaceCells(m_AllCells,m_grid);
 756   getSurfaceCells(m_bcs, m_SelectedCells, m_grid);
 757   EG_VTKDCC(vtkIntArray, cell_code, m_grid, "cell_code");
 758   EG_VTKDCN(vtkDoubleArray, node_meshdensity, m_grid, "node_meshdensity");
 759   m_SelectedNodes.clear();
 760   getSurfaceNodes(m_bcs,m_SelectedNodes,m_grid);
 761   getNodesFromCells(m_AllCells, nodes, m_grid);
 762   setGrid(m_grid);
 763   setCells(m_AllCells);
 764   cout<<"m_AllCells.size()="<<m_AllCells.size()<<endl;
 765
 766   N_removed_EP=0;
 767
 768   N_points=m_grid->GetNumberOfPoints();
 769   N_cells=m_grid->GetNumberOfCells();
 770   N_newpoints=0;
 771   N_newcells=0;
 772
 773   hitlist.resize(N_points);
 774   offset.resize(N_points);
 775
 776   marked_cells.clear();
 777   marked_nodes.clear();
 778
 779   remove_EP_counter();
 780
 781     //unmark cells (TODO: optimize)
 782   marked_cells.clear();
 783     //init grid_tmp
 784   N_points=m_grid->GetNumberOfPoints();
 785   N_cells=m_grid->GetNumberOfCells();
 786   cout<<"N_points="<<N_points<<endl;
 787   cout<<"N_cells="<<N_cells<<endl;
 788   cout<<"N_newpoints="<<N_newpoints<<endl;
 789   cout<<"N_newcells="<<N_newcells<<endl;
 790   EG_VTKSP(vtkUnstructuredGrid,grid_tmp);
 791   allocateGrid(grid_tmp,N_cells+N_newcells,N_points+N_newpoints);
 792   m_total_N_newpoints+=N_newpoints; m_total_N_newcells+=N_newcells;
 793
 794   makeCopyNoAlloc(m_grid, grid_tmp);
 795     //initialize new node counter
 796   m_newNodeId=N_points;
 797
 798   remove_EP_actor(grid_tmp);
 799
 800   makeCopy(grid_tmp,m_grid);
 801
 802   cout<<"===remove_EP_all END==="<<endl;
 803   return(0);
 804 }
 805
 806 int CreateSpecialMapping::FullEdit()
 807 {
 808   cout<<"===FullEdit START==="<<endl;
 809
 810   N_inserted_FP=0;
 811   N_inserted_EP=0;
 812   N_removed_FP=0;
 813   N_removed_EP=0;
 814
 815   N_points=m_grid->GetNumberOfPoints();
 816   N_cells=m_grid->GetNumberOfCells();
 817   N_newpoints=0;
 818   N_newcells=0;
 819
 820   hitlist.resize(N_points);
 821   offset.resize(N_points);
 822
 823   marked_cells.clear();
 824   marked_nodes.clear();
 825
 826   if(insert_FP) insert_FP_counter();
 827   if(insert_EP) insert_EP_counter();
 828   if(remove_FP) remove_FP_counter();
 829   if(remove_EP) remove_EP_counter();
 830
 831   cout<<"================="<<endl;
 832   cout<<"hitlist="<<hitlist<<endl;
 833   cout<<"================="<<endl;
 834
 835     //unmark cells (TODO: optimize)
 836   marked_cells.clear();
 837     //init grid_tmp
 838   N_points=m_grid->GetNumberOfPoints();
 839   N_cells=m_grid->GetNumberOfCells();
 840   cout<<"N_points="<<N_points<<endl;
 841   cout<<"N_cells="<<N_cells<<endl;
 842   cout<<"N_newpoints="<<N_newpoints<<endl;
 843   cout<<"N_newcells="<<N_newcells<<endl;
 844   EG_VTKSP(vtkUnstructuredGrid,grid_tmp);
 845   allocateGrid(grid_tmp,N_cells+N_newcells,N_points+N_newpoints);
 846   m_total_N_newpoints+=N_newpoints; m_total_N_newcells+=N_newcells;
 847
 848   makeCopyNoAlloc(m_grid, grid_tmp);//TODO: This will not work if the size of the grid is reduced!
 849     //initialize new node counter
 850   m_newNodeId=N_points;
 851
 852   if(insert_FP) insert_FP_actor(grid_tmp);
 853   if(insert_EP) insert_EP_actor(grid_tmp);
 854
 855   cout<<"================="<<endl;
 856   cout<<"hitlist="<<hitlist<<endl;
 857   cout<<"================="<<endl;
 858   if(remove_FP) remove_FP_actor(grid_tmp);
 859   if(remove_EP) remove_EP_actor(grid_tmp);
 860
 861   makeCopy(grid_tmp,m_grid);
 862
 863   cout<<"===FullEdit END==="<<endl;
 864   return(0);
 865 }
 866
 867 // vtkIdType CreateSpecialMapping::FindSnapPoint(vtkUnstructuredGrid *src, vtkIdType DeadNode)
 868 vtkIdType CreateSpecialMapping::FindSnapPoint(vtkUnstructuredGrid *src, vtkIdType DeadNode,QSet <vtkIdType> & DeadCells,QSet <vtkIdType> & MutatedCells,QSet <vtkIdType> & MutilatedCells)
 869 {
 870   getAllSurfaceCells(m_AllCells,src);
 871   getSurfaceCells(m_bcs, m_SelectedCells, src);
 872   m_SelectedNodes.clear();
 873   getSurfaceNodes(m_bcs,m_SelectedNodes,src);
 874   getNodesFromCells(m_AllCells, nodes, src);
 875   setGrid(src);
 876   setCells(m_AllCells);
 877
 878   UpdateNodeType();
 879
 880   EG_VTKDCN(vtkCharArray, node_type, src, "node_type");
 881   if(node_type->GetValue(DeadNode)==VTK_FIXED_VERTEX)
 882   {
 883     cout<<"Sorry, unable to remove fixed vertex."<<endl;
 884     return(-1);
 885   }
 886
 887     //src grid info
 888   N_points=src->GetNumberOfPoints();
 889   N_cells=src->GetNumberOfCells();
 890   N_newpoints=-1;
 891   N_newcells=0;
 892
 893   vtkIdType SnapPoint=-1;
 894
 895   foreach(vtkIdType PSP, n2n[DeadNode])
 896   {
 897     bool IsValidSnapPoint=true;
 898
 899     if(DebugLevel>10) cout<<"====>PSP="<<PSP<<endl;
 900     bool IsTetra=true;
 901     if(NumberOfCommonPoints(DeadNode,PSP,IsTetra)>2)//common point check
 902     {
 903       if(DebugLevel>10) cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
 904       IsValidSnapPoint=false;
 905     }
 906     if(IsTetra)//tetra check
 907     {
 908       if(DebugLevel>10) cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
 909       IsValidSnapPoint=false;
 910     }
 911
 912     //count number of points and cells to remove + analyse cell transformations
 913     N_newpoints=-1;
 914     N_newcells=0;
 915     DeadCells.clear();
 916     MutatedCells.clear();
 917     MutilatedCells.clear();
 918     foreach(vtkIdType C, n2c[DeadNode])//loop through potentially dead cells
 919     {
 920       //get points around cell
 921       vtkIdType N_pts, *pts;
 922       src->GetCellPoints(C, N_pts, pts);
 923
 924       bool ContainsSnapPoint=false;
 925       bool invincible=false;
 926       for(int i=0;i<N_pts;i++)
 927       {
 928         if(DebugLevel>10) cout<<"pts["<<i<<"]="<<pts[i]<<" and PSP="<<PSP<<endl;
 929         if(pts[i]==PSP) {ContainsSnapPoint=true;}
 930         if(pts[i]!=DeadNode && pts[i]!=PSP &&  n2c[pts[i]].size()<=1) invincible=true;
 931       }
 932       if(ContainsSnapPoint)
 933       {
 934         if(N_pts==3)//dead cell
 935         {
 936           if(invincible)//Check that empty lines aren't left behind when a cell is killed
 937           {
 938             if(DebugLevel>10) cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
 939             IsValidSnapPoint=false;
 940           }
 941           else
 942           {
 943             DeadCells.insert(C);
 944             N_newcells-=1;
 945             if(DebugLevel>10) cout<<"cell "<<C<<" has been pwned!"<<endl;
 946           }
 947         }
 948         else
 949         {
 950           cout<<"RED ALERT: Xenomorph detected!"<<endl;
 951           EG_BUG;
 952         }
 953       }
 954       else
 955       {
 956         vtkIdType src_N_pts, *src_pts;
 957         src->GetCellPoints(C, src_N_pts, src_pts);
 958
 959         if(src_N_pts!=3)
 960         {
 961           cout<<"RED ALERT: Xenomorph detected!"<<endl;
 962           EG_BUG;
 963         }
 964
 965         vtkIdType OldTriangle[3];
 966         vtkIdType NewTriangle[3];
 967
 968         for(int i=0;i<src_N_pts;i++)
 969         {
 970           OldTriangle[i]=src_pts[i];
 971           NewTriangle[i]=( (src_pts[i]==DeadNode) ? PSP : src_pts[i] );
 972         }
 973         vec3_t Old_N= triNormal(src, OldTriangle[0], OldTriangle[1], OldTriangle[2]);
 974         vec3_t New_N= triNormal(src, NewTriangle[0], NewTriangle[1], NewTriangle[2]);
 975         double OldArea=Old_N.abs();
 976         double NewArea=New_N.abs();
 977         double scal=Old_N*New_N;
 978         double cross=(Old_N.cross(New_N)).abs();//double-cross on Nar Shadaa B-)
 979
 980         if(DebugLevel>10) {
 981           cout<<"OldArea="<<OldArea<<endl;
 982           cout<<"NewArea="<<NewArea<<endl;
 983           cout<<"scal="<<scal<<endl;
 984           cout<<"cross="<<cross<<endl;
 985         }
 986
 987         if(Old_N*New_N<Old_N*Old_N*1./100.)//area + inversion check
 988         {
 989           if(DebugLevel>10) cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
 990           IsValidSnapPoint=false;
 991         }
 992
 993         //mutated cell
 994         MutatedCells.insert(C);
 995         if(DebugLevel>10) cout<<"cell "<<C<<" has been infected!"<<endl;
 996       }
 997     }
 998
 999     if(N_cells+N_newcells<=0)//survivor check
1000     {
1001       if(DebugLevel>10) cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
1002       IsValidSnapPoint=false;
1003     }
1004     if(node_type->GetValue(DeadNode)==VTK_BOUNDARY_EDGE_VERTEX && node_type->GetValue(PSP)==VTK_SIMPLE_VERTEX)
1005     {
1006       if(DebugLevel>10) cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
1007       IsValidSnapPoint=false;
1008     }
1009
1010     if(IsValidSnapPoint) {SnapPoint=PSP; break;}
1011   }//end of loop through potential SnapPoints
1012
1013   if(DebugLevel>10)
1014   {
1015     cout<<"AT FINDSNAPPOINT EXIT"<<endl;
1016     cout<<"N_points="<<N_points<<endl;
1017     cout<<"N_cells="<<N_cells<<endl;
1018     cout<<"N_newpoints="<<N_newpoints<<endl;
1019     cout<<"N_newcells="<<N_newcells<<endl;
1020   }
1021   return(SnapPoint);
1022 }
1023
1024 bool CreateSpecialMapping::DeletePoint_2(vtkUnstructuredGrid *src, vtkIdType DeadNode)
1025 {
1026   getAllSurfaceCells(m_AllCells,src);
1027   getSurfaceCells(m_bcs, m_SelectedCells, src);
1028   m_SelectedNodes.clear();
1029   getSurfaceNodes(m_bcs,m_SelectedNodes,src);
1030   getNodesFromCells(m_AllCells, nodes, src);
1031   setGrid(src);
1032   setCells(m_AllCells);
1033
1034   //src grid info
1035   N_points=src->GetNumberOfPoints();
1036   N_cells=src->GetNumberOfCells();
1037   N_newpoints=-1;
1038   N_newcells=0;
1039
1040   if(DeadNode<0 || DeadNode>=N_points)
1041   {
1042     cout<<"Warning: Point out of range: DeadNode="<<DeadNode<<" N_points="<<N_points<<endl;
1043     return(false);
1044   }
1045
1046   QSet <vtkIdType> DeadCells;
1047   QSet <vtkIdType> MutatedCells;
1048   QSet <vtkIdType> MutilatedCells;
1049
1050   if(DebugLevel>10) {
1051     cout<<"BEFORE FINDSNAPPOINT"<<endl;
1052     cout<<"N_points="<<N_points<<endl;
1053     cout<<"N_cells="<<N_cells<<endl;
1054     cout<<"N_newpoints="<<N_newpoints<<endl;
1055     cout<<"N_newcells="<<N_newcells<<endl;
1056   }
1057   vtkIdType SnapPoint=FindSnapPoint(src,DeadNode,DeadCells,MutatedCells,MutilatedCells);
1058
1059   if(DebugLevel>0) cout<<"===>DeadNode="<<DeadNode<<" moving to SNAPPOINT="<<SnapPoint<<" DebugLevel="<<DebugLevel<<endl;
1060   if(SnapPoint<0) {cout<<"Sorry no possible SnapPoint found."<<endl; return(false);}
1061
1062   //allocate
1063   if(DebugLevel>10) {
1064     cout<<"BEFORE ALLOCATION"<<endl;
1065     cout<<"N_points="<<N_points<<endl;
1066     cout<<"N_cells="<<N_cells<<endl;
1067     cout<<"N_newpoints="<<N_newpoints<<endl;
1068     cout<<"N_newcells="<<N_newcells<<endl;
1069   }
1070   N_points=src->GetNumberOfPoints();
1071   N_cells=src->GetNumberOfCells();
1072   cout<<"N_points="<<N_points<<endl;
1073   cout<<"N_cells="<<N_cells<<endl;
1074   cout<<"N_newpoints="<<N_newpoints<<endl;
1075   cout<<"N_newcells="<<N_newcells<<endl;
1076   EG_VTKSP(vtkUnstructuredGrid,dst);
1077   allocateGrid(dst,N_cells+N_newcells,N_points+N_newpoints);
1078   m_total_N_newpoints+=N_newpoints; m_total_N_newcells+=N_newcells;
1079
1080   //vector used to redefine the new point IDs
1081   QVector <vtkIdType> OffSet(N_points);
1082
1083   //copy undead points
1084   vtkIdType dst_id_node=0;
1085   for (vtkIdType src_id_node = 0; src_id_node < N_points; src_id_node++) {//loop through src points
1086     if(src_id_node!=DeadNode)//if the node isn't dead, copy it
1087     {
1088       vec3_t x;
1089       src->GetPoints()->GetPoint(src_id_node, x.data());
1090       dst->GetPoints()->SetPoint(dst_id_node, x.data());
1091       copyNodeData(src, src_id_node, dst, dst_id_node);
1092       OffSet[src_id_node]=src_id_node-dst_id_node;
1093       dst_id_node++;
1094     }
1095     else
1096     {
1097       if(DebugLevel>0) cout<<"src_id_node="<<src_id_node<<" dst_id_node="<<dst_id_node<<endl;
1098     }
1099   };
1100   if(DebugLevel>10) {
1101     cout<<"DeadCells="<<DeadCells<<endl;
1102     cout<<"MutatedCells="<<MutatedCells<<endl;
1103     cout<<"MutilatedCells="<<MutilatedCells<<endl;
1104   }
1105   //Copy undead cells
1106   for (vtkIdType id_cell = 0; id_cell < src->GetNumberOfCells(); ++id_cell) {//loop through src cells
1107     if(!DeadCells.contains(id_cell))//if the cell isn't dead
1108     {
1109       vtkIdType src_N_pts, *src_pts;
1110       vtkIdType dst_N_pts, *dst_pts;
1111       src->GetCellPoints(id_cell, src_N_pts, src_pts);
1112
1113       vtkIdType type_cell = src->GetCellType(id_cell);
1114       if(DebugLevel>10) cout<<"-->id_cell="<<id_cell<<endl;
1115       if(DebugLevel>10) for(int i=0;i<src_N_pts;i++) cout<<"src_pts["<<i<<"]="<<src_pts[i]<<endl;
1116 //       src->GetCellPoints(id_cell, dst_N_pts, dst_pts);
1117       dst_N_pts=src_N_pts;
1118       dst_pts=new vtkIdType[dst_N_pts];
1119       if(MutatedCells.contains(id_cell))//mutated cell
1120       {
1121         if(DebugLevel>10) cout<<"processing mutated cell "<<id_cell<<endl;
1122         for(int i=0;i<src_N_pts;i++)
1123         {
1124           if(src_pts[i]==DeadNode) {
1125             if(DebugLevel>10) {
1126               cout<<"SnapPoint="<<SnapPoint<<endl;
1127               cout<<"OffSet[SnapPoint]="<<OffSet[SnapPoint]<<endl;
1128               cout<<"src_pts["<<i<<"]="<<src_pts[i]<<endl;
1129             }
1130             dst_pts[i]=SnapPoint-OffSet[SnapPoint];
1131           }
1132           else dst_pts[i]=src_pts[i]-OffSet[src_pts[i]];
1133         }
1134         if(DebugLevel>10) cout<<"--->dst_pts:"<<endl;
1135         if(DebugLevel>10) for(int i=0;i<dst_N_pts;i++) cout<<"dst_pts["<<i<<"]="<<dst_pts[i]<<endl;
1136
1137       }
1138       else if(MutilatedCells.contains(id_cell))//mutilated cell
1139       {
1140         if(DebugLevel>10) cout<<"processing mutilated cell "<<id_cell<<endl;
1141
1142         if(type_cell==VTK_QUAD) {
1143           type_cell=VTK_TRIANGLE;
1144           dst_N_pts-=1;
1145         }
1146         else {cout<<"FATAL ERROR: Unknown mutilated cell detected! It is not a quad! Potential xenomorph infestation!"<<endl;EG_BUG;}
1147         //merge points
1148         int j=0;
1149         for(int i=0;i<src_N_pts;i++)
1150         {
1151           if(src_pts[i]==SnapPoint) { dst_pts[j]=SnapPoint-OffSet[SnapPoint];j++; }//SnapPoint
1152           else if(src_pts[i]!=DeadNode) { dst_pts[j]=src_pts[i]-OffSet[src_pts[i]];j++; }//pre-snap/dead + post-snap/dead
1153           //do nothing in case of DeadNode
1154         }
1155       }
1156       else//normal cell
1157       {
1158         if(DebugLevel>10) cout<<"processing normal cell "<<id_cell<<endl;
1159         for(int i=0;i<src_N_pts;i++)
1160         {
1161           dst_pts[i]=src_pts[i]-OffSet[src_pts[i]];
1162         }
1163       }
1164       //copy the cell
1165       vtkIdType id_new_cell = dst->InsertNextCell(type_cell, dst_N_pts, dst_pts);
1166       copyCellData(src, id_cell, dst, id_new_cell);
1167       if(DebugLevel>10) {
1168         cout<<"===Copying cell "<<id_cell<<" to "<<id_new_cell<<"==="<<endl;
1169         cout<<"src_pts:"<<endl;
1170         for(int i=0;i<src_N_pts;i++) cout<<"src_pts["<<i<<"]="<<src_pts[i]<<endl;
1171         cout<<"dst_pts:"<<endl;
1172         for(int i=0;i<dst_N_pts;i++) cout<<"dst_pts["<<i<<"]="<<dst_pts[i]<<endl;
1173         cout<<"OffSet="<<OffSet<<endl;
1174         cout<<"===Copying cell end==="<<endl;
1175       }
1176       delete dst_pts;
1177     }
1178   };
1179 //   cout_grid(cout,dst,true,true,true,true);
1180   makeCopy(dst, src);
1181   return(true);
1182 }
1183
1184 int CreateSpecialMapping::remove_EP_all_2()
1185 {
1186   cout<<"===remove_EP_all_2 START==="<<endl;
1187   getAllSurfaceCells(m_AllCells,m_grid);
1188   getSurfaceCells(m_bcs, m_SelectedCells, m_grid);
1189   EG_VTKDCC(vtkIntArray, cell_code, m_grid, "cell_code");
1190   EG_VTKDCN(vtkDoubleArray, node_meshdensity, m_grid, "node_meshdensity");
1191   m_SelectedNodes.clear();
1192   getSurfaceNodes(m_bcs,m_SelectedNodes,m_grid);
1193   getNodesFromCells(m_AllCells, nodes, m_grid);
1194   setGrid(m_grid);
1195   setCells(m_AllCells);
1196   cout<<"m_AllCells.size()="<<m_AllCells.size()<<endl;
1197
1198   N_removed_EP=0;
1199
1200   N_points=m_grid->GetNumberOfPoints();
1201   N_cells=m_grid->GetNumberOfCells();
1202   N_newpoints=0;
1203   N_newcells=0;
1204
1205   hitlist.clear();
1206   offset.clear();
1207   hitlist.resize(N_points);
1208   offset.resize(N_points);
1209
1210   marked_cells.clear();
1211   marked_nodes.clear();
1212
1213   remove_EP_counter();
1214   cout<<"================="<<endl;
1215   cout<<"hitlist="<<hitlist<<endl;
1216   cout<<"================="<<endl;
1217
1218   int kills=0;
1219   int contracts=0;
1220   for(int i=0;i<hitlist.size();i++)
1221   {
1222     if(hitlist[i]==2){
1223       contracts++;
1224       cout<<"Deleting point "<<i<<" currently known as "<<i-kills<<endl;
1225
1226       QString num1;num1.setNum(i);
1227       QString num2;num2.setNum(i-kills);
1228       GuiMainWindow::pointer()->QuickSave("pre-deleting_"+num1+"_"+num2+".vtu");
1229
1230       if(DeletePoint_2(m_grid,i-kills))
1231       {
1232         kills++;
1233         cout<<"Kill successful"<<endl;
1234       }
1235       else
1236       {
1237         cout<<"Kill failed"<<endl;
1238         N_removed_EP--;
1239       }
1240
1241       GuiMainWindow::pointer()->QuickSave("post-deleting_"+num1+"_"+num2+".vtu");
1242
1243     }
1244   }
1245   cout<<"Killed: "<<kills<<"/"<<contracts<<endl;
1246   if(kills!=contracts) {cout<<"MISSION FAILED"<<endl;EG_BUG;}
1247   cout<<"===remove_EP_all_2 END==="<<endl;
1248   return(0);
1249 }
1250
1251 int CreateSpecialMapping::remove_FP_all_2()
1252 {
1253   cout<<"===remove_FP_all_2 START==="<<endl;
1254   cout<<"+++++++"<<endl;
1255   cout_grid(cout,m_grid,true,true,true,true);
1256   cout<<"+++++++"<<endl;
1257
1258   getAllSurfaceCells(m_AllCells,m_grid);
1259   getSurfaceCells(m_bcs, m_SelectedCells, m_grid);
1260   EG_VTKDCC(vtkIntArray, cell_code, m_grid, "cell_code");
1261   EG_VTKDCN(vtkDoubleArray, node_meshdensity, m_grid, "node_meshdensity");
1262   m_SelectedNodes.clear();
1263   getSurfaceNodes(m_bcs,m_SelectedNodes,m_grid);
1264   getNodesFromCells(m_AllCells, nodes, m_grid);
1265   setGrid(m_grid);
1266   setCells(m_AllCells);
1267   cout<<"m_AllCells.size()="<<m_AllCells.size()<<endl;
1268
1269   N_removed_FP=0;
1270
1271   N_points=m_grid->GetNumberOfPoints();
1272   N_cells=m_grid->GetNumberOfCells();
1273   N_newpoints=0;
1274   N_newcells=0;
1275
1276   hitlist.clear();
1277   offset.clear();
1278   hitlist.resize(N_points);
1279   offset.resize(N_points);
1280
1281   marked_cells.clear();
1282   marked_nodes.clear();
1283
1284   remove_FP_counter();
1285   cout_grid(cout,m_grid);
1286   cout<<"================="<<endl;
1287   cout<<"hitlist="<<hitlist<<endl;
1288   cout<<"================="<<endl;
1289
1290   int kills=0;
1291   int contracts=0;
1292   for(int i=0;i<hitlist.size();i++)
1293   {
1294     if(hitlist[i]==1){
1295       contracts++;
1296       cout<<"Deleting point "<<i<<" currently known as "<<i-kills<<endl;
1297
1298       QString num1;num1.setNum(i);
1299       QString num2;num2.setNum(i-kills);
1300       GuiMainWindow::pointer()->QuickSave("pre-deleting_"+num1+"_"+num2+".vtu");
1301
1302       if(DeletePoint_2(m_grid,i-kills))
1303       {
1304         kills++;
1305         cout<<"Kill successful"<<endl;
1306       }
1307       else
1308       {
1309         cout<<"Kill failed"<<endl;
1310         N_removed_FP--;
1311       }
1312
1313       GuiMainWindow::pointer()->QuickSave("post-deleting_"+num1+"_"+num2+".vtu");
1314
1315     }
1316   }
1317   cout<<"Killed: "<<kills<<"/"<<contracts<<endl;
1318   if(kills!=contracts) {cout<<"MISSION FAILED"<<endl;EG_BUG;}
1319   cout<<"===remove_FP_all_2 END==="<<endl;
1320   return(0);
1321 }