src/insertpoints.cpp

   1 //
   2 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   3 // +                                                                      +
   4 // + This file is part of enGrid.                                         +
   5 // +                                                                      +
   6 // + Copyright 2008,2009 Oliver Gloth                                     +
   7 // +                                                                      +
   8 // + enGrid is free software: you can redistribute it and/or modify       +
   9 // + it under the terms of the GNU General Public License as published by +
  10 // + the Free Software Foundation, either version 3 of the License, or    +
  11 // + (at your option) any later version.                                  +
  12 // +                                                                      +
  13 // + enGrid is distributed in the hope that it will be useful,            +
  14 // + but WITHOUT ANY WARRANTY; without even the implied warranty of       +
  15 // + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the        +
  16 // + GNU General Public License for more details.                         +
  17 // +                                                                      +
  18 // + You should have received a copy of the GNU General Public License    +
  19 // + along with enGrid. If not, see <http://www.gnu.org/licenses/>.       +
  20 // +                                                                      +
  21 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  22 //
  23 #include "insertpoints.h"
  24
  25 #include "guimainwindow.h"
  26
  27 #include <vtkCharArray.h>
  28
  29 #include <QTime>
  30
  31 InsertPoints::InsertPoints()
  32 : SurfaceOperation()
  33 {
  34   setQuickSave(true);
  35 }
  36
  37 void InsertPoints::operate()
  38 {
  39   int N1 = grid->GetNumberOfPoints();
  40   if(insert_FP) insert_FP_all();
  41   if(insert_EP) insert_EP_all();
  42   int N2 = grid->GetNumberOfPoints();
  43   m_NumInserted = N2 - N1;
  44 }
  45
  46 bool InsertPoints::insert_fieldpoint(vtkIdType id_cell)
  47 {
  48   double Fred1=1.0/sqrt(3);
  49   double Qmin=1.1;//1.189;
  50   double total=0;
  51   for(int i=0;i<3;i++)
  52   {
  53     vtkIdType cell = getPartC2C()[id_cell][i];
  54     if( cell != -1 ) total += Q_L(cell);
  55   }
  56   return ( Q_L(id_cell)>1.0/Fred1 && total>3*Qmin );
  57 }
  58
  59 bool InsertPoints::insert_edgepoint(vtkIdType id_node1, vtkIdType id_node2)
  60 {
  61   EG_BUG;
  62   bool result = distance(grid, id_node1, id_node2) > 0.5 * ( desiredEdgeLength(id_node1) + desiredEdgeLength(id_node2) );
  63   return ( result );
  64 }
  65
  66 bool InsertPoints::SplitSide(vtkIdType id_cell,int side)
  67 {
  68   vtkIdType N_pts,*pts;
  69   grid->GetCellPoints(id_cell,N_pts,pts);
  70   return( insert_edgepoint(pts[side],pts[(side+1)%N_pts]) );
  71 }
  72
  73 int InsertPoints::insert_FP_all()
  74 {
  75   //cout<<"===insert_FP_all START==="<<endl;
  76   QTime start = QTime::currentTime();
  77
  78   setAllSurfaceCells();
  79   UpdatePotentialSnapPoints(true);
  80
  81   QVector <vtkIdType> l_SelectedCells;
  82   getSurfaceCells(m_bcs, l_SelectedCells, grid);
  83
  84   QVector <bool> l_marked_cells(l_SelectedCells.size());
  85
  86   int l_N_newpoints=0;
  87   int l_N_newcells=0;
  88
  89   //counter
  90   for(int i_cell=0;i_cell<l_SelectedCells.size();i_cell++)
  91   {
  92     vtkIdType id_cell = l_SelectedCells[i_cell];
  93     if( insert_fieldpoint(id_cell) )
  94     {
  95       l_marked_cells[i_cell] = true;
  96       l_N_newcells += 2;
  97       l_N_newpoints += 1;
  98     }
  99   }
 100
 101   //initialize grid_tmp
 102   int l_N_points = grid->GetNumberOfPoints();
 103   int l_N_cells = grid->GetNumberOfCells();
 104   EG_VTKSP(vtkUnstructuredGrid,grid_tmp);
 105   allocateGrid(grid_tmp,l_N_cells+l_N_newcells,l_N_points+l_N_newpoints);
 106   makeCopyNoAlloc(grid, grid_tmp);
 107
 108   //initialize new node counter
 109   vtkIdType l_newNodeId = l_N_points;
 110
 111   //actor
 112   for(int i_cell=0;i_cell<l_SelectedCells.size();i_cell++)
 113   {
 114     vtkIdType id_cell = l_SelectedCells[i_cell];
 115     if( l_marked_cells[i_cell] )
 116     {
 117       vtkIdType N_pts, *pts;
 118       grid->GetCellPoints(id_cell, N_pts, pts);
 119       vec3_t C(0,0,0);
 120       for(int i=0;i<N_pts;i++)
 121       {
 122         vec3_t corner;
 123         grid->GetPoints()->GetPoint(pts[i], corner.data());
 124         C+=corner;
 125       }
 126       C=(1/(double)N_pts)*C;
 127
 128       //C=project(C);
 129       grid_tmp->GetPoints()->SetPoint(l_newNodeId,C.data());
 130       copyNodeData(grid_tmp,pts[0],grid_tmp,l_newNodeId);
 131       EG_VTKDCN(vtkCharArray, node_type, grid_tmp, "node_type");
 132       node_type->SetValue(l_newNodeId, VTK_SIMPLE_VERTEX);
 133
 134       for(int i=0;i<N_pts;i++)
 135       {
 136         vtkIdType pts_triangle[3];
 137         pts_triangle[0]=pts[i];
 138         pts_triangle[1]=pts[(i+1)%N_pts];
 139         pts_triangle[2]=l_newNodeId;
 140         if(i==0)
 141         {
 142           grid_tmp->ReplaceCell(id_cell , 3, pts_triangle);
 143         }
 144         else
 145         {
 146           vtkIdType newCellId = grid_tmp->InsertNextCell(VTK_TRIANGLE,3,pts_triangle);
 147           copyCellData(grid_tmp,id_cell,grid_tmp,newCellId);
 148         }
 149       }
 150       l_newNodeId++;
 151     }
 152   }
 153
 154   //update grid
 155   makeCopy(grid_tmp,grid);
 156
 157   //cout << start.msecsTo(QTime::currentTime()) << " milliseconds elapsed" << endl;
 158   //cout<<"===insert_FP_all END==="<<endl;
 159   return(0);
 160 }
 161
 162 int InsertPoints::insert_EP_all()
 163 {
 164   l2g_t  cells = getPartCells();
 165   g2l_t _cells = getPartLocalCells();
 166
 167   //cout<<"===insert_EP_all START==="<<endl;
 168   QTime start = QTime::currentTime();
 169
 170   setAllSurfaceCells();
 171   UpdatePotentialSnapPoints(true);
 172
 173   EG_VTKDCC(vtkIntArray, cell_code, grid, "cell_code");
 174
 175   int num_newpoints=0;
 176   int num_newcells=0;
 177
 178   QVector <int> marked_cells(cells.size(), 0);
 179   QVector <stencil_t> stencil_vector(cells.size());
 180
 181   //counter
 182   for (int i = 0; i < cells.size(); ++i) {
 183     vtkIdType id_cell = cells[i];
 184     if (m_bcs.contains(cell_code->GetValue(id_cell)) && (grid->GetCellType(id_cell) == VTK_TRIANGLE)) {//if selected and triangle cell
 185       int j_split = -1;
 186       double L_max = 0;
 187       vtkIdType N_pts, *pts;
 188       grid->GetCellPoints(id_cell, N_pts, pts);
 189       for (int j = 0; j < 3; ++j) {
 190         vtkIdType id_node1 = pts[j];
 191         vtkIdType id_node2 = pts[(j+1)%N_pts];
 192         double L = distance(grid, id_node1, id_node2);
 193         if (L > max(desiredEdgeLength(id_node1), desiredEdgeLength(id_node2))) {
 194           if (L > L_max) {
 195             j_split = j;
 196             L_max = L;
 197           }
 198         }
 199       }
 200       if (j_split != -1) {
 201         stencil_t S = getStencil(id_cell, j_split);
 202         if (S.twocells && (S.neighbour_type == VTK_TRIANGLE)) {
 203           if (!marked_cells[i]  && !marked_cells[_cells[S.id_cell2]]) {
 204             stencil_vector[i] = S;
 205             marked_cells[i] = 1;
 206             marked_cells[_cells[S.id_cell2]] = 2;
 207             num_newpoints++;
 208             num_newcells += 2;
 209           }
 210         } else if (!S.twocells) {
 211           if (!marked_cells[i]) {
 212             stencil_vector[i] = S;
 213             marked_cells[i] = 1;
 214             num_newpoints++;
 215             num_newcells+=1;
 216           }
 217         }
 218       } //end of loop through sides
 219     } //end of if selected and triangle cell
 220   } //end of counter loop
 221
 222   //initialize grid_tmp
 223   int l_N_points=grid->GetNumberOfPoints();
 224   int l_N_cells=grid->GetNumberOfCells();
 225   EG_VTKSP(vtkUnstructuredGrid,grid_tmp);
 226   allocateGrid(grid_tmp, l_N_cells + num_newcells, l_N_points + num_newpoints);
 227   makeCopyNoAlloc(grid, grid_tmp);
 228
 229   //initialize new node counter
 230   vtkIdType l_newNodeId = l_N_points;
 231
 232   //actor
 233   for (int i = 0; i < cells.size(); i++) {
 234     if (marked_cells[i] == 1) {
 235       stencil_t S = stencil_vector[i];
 236
 237       //calculate midpoint
 238       vec3_t A,B;
 239       grid_tmp->GetPoint(S.p[1],A.data());
 240       grid_tmp->GetPoint(S.p[3],B.data());
 241       vec3_t M=0.5*(A+B);
 242
 243       //project point
 244       //M=project(M);
 245       //add point
 246       grid_tmp->GetPoints()->SetPoint(l_newNodeId, M.data());
 247       copyNodeData(grid_tmp,S.p[1],grid_tmp,l_newNodeId);
 248
 249       // inserted edge point = type of the edge on which it is inserted
 250       EG_VTKDCN(vtkCharArray, node_type, grid_tmp, "node_type");
 251       node_type->SetValue(l_newNodeId, getNewNodeType(S) );
 252
 253       if(S.twocells && S.neighbour_type==VTK_TRIANGLE){//2 triangles
 254         //four new triangles
 255         vtkIdType pts_triangle[4][3];
 256         for(int i=0;i<4;i++)
 257         {
 258           pts_triangle[i][0]=S.p[i];
 259           pts_triangle[i][1]=S.p[(i+1)%4];
 260           pts_triangle[i][2]=l_newNodeId;
 261         }
 262
 263         grid_tmp->ReplaceCell(S.id_cell1 , 3, pts_triangle[0]);
 264         grid_tmp->ReplaceCell(S.id_cell2 , 3, pts_triangle[1]);
 265
 266         vtkIdType newCellId;
 267
 268         newCellId = grid_tmp->InsertNextCell(VTK_TRIANGLE,3,pts_triangle[2]);
 269         copyCellData(grid_tmp,S.id_cell2,grid_tmp,newCellId);
 270
 271         newCellId = grid_tmp->InsertNextCell(VTK_TRIANGLE,3,pts_triangle[3]);
 272         copyCellData(grid_tmp,S.id_cell1,grid_tmp,newCellId);
 273       }
 274       else if(!S.twocells) {//1 triangle
 275         //two new triangles
 276         vtkIdType pts_triangle[2][3];
 277         pts_triangle[0][0]=S.p[0];
 278         pts_triangle[0][1]=S.p[1];
 279         pts_triangle[0][2]=l_newNodeId;
 280         pts_triangle[1][0]=S.p[3];
 281         pts_triangle[1][1]=S.p[0];
 282         pts_triangle[1][2]=l_newNodeId;
 283
 284         grid_tmp->ReplaceCell(S.id_cell1 , 3, pts_triangle[0]);
 285
 286         vtkIdType newCellId;
 287         newCellId = grid_tmp->InsertNextCell(VTK_TRIANGLE,3,pts_triangle[1]);
 288         copyCellData(grid_tmp,S.id_cell1,grid_tmp,newCellId);
 289       }
 290       else {
 291         cout<<"I DON'T KNOW HOW TO SPLIT THIS CELL!!!"<<endl;
 292         EG_BUG;
 293       }
 294
 295       //increment ID
 296       l_newNodeId++;
 297     }
 298   }
 299
 300   //update grid
 301   makeCopy(grid_tmp,grid);
 302
 303   //cout << start.msecsTo(QTime::currentTime()) << " milliseconds elapsed" << endl;
 304   //cout<<"===insert_EP_all END==="<<endl;
 305   return(0);
 306 }
 307
 308 char InsertPoints::getNewNodeType(stencil_t S)
 309 {
 310   vtkIdType id_node1 = S.p[1];
 311   vtkIdType id_node2 = S.p[3];
 312
 313   EG_VTKDCN(vtkCharArray, node_type, grid, "node_type");
 314   if( node_type->GetValue(id_node1)==VTK_SIMPLE_VERTEX || node_type->GetValue(id_node2)==VTK_SIMPLE_VERTEX ) {
 315     return VTK_SIMPLE_VERTEX;
 316   }
 317   else {
 318     QVector <vtkIdType> PSP = getPotentialSnapPoints(id_node1);
 319     if( PSP.contains(id_node2) ) {
 320       EG_VTKDCC(vtkIntArray, cell_code, grid, "cell_code");
 321       if( cell_code->GetValue(S.id_cell1) != cell_code->GetValue(S.id_cell2) ) {
 322         return VTK_BOUNDARY_EDGE_VERTEX;
 323       }
 324       else {
 325         return VTK_FEATURE_EDGE_VERTEX;
 326       }
 327     }
 328     else {
 329       return VTK_SIMPLE_VERTEX;
 330     }
 331   }
 332 }
 333
 334        //============================================
 335       ///@@@  TODO: Update node info (densities+type) Still necessary?
 336 // EG_VTKDCN(vtkIntArray, node_specified_density, grid_tmp, "node_specified_density");//density index from table
 337 // EG_VTKDCN(vtkDoubleArray, node_meshdensity_desired, grid_tmp, "node_meshdensity_desired");//what we want
 338 // EG_VTKDCN(vtkDoubleArray, node_meshdensity_current, grid_tmp, "node_meshdensity_current");//what we have
 339 // EG_VTKDCN(vtkCharArray, node_type, grid_tmp, "node_type");//node type
 340       //============================================
 341
 342 //       //part 1
 343 //       node_type->SetValue(l_newNodeId,VTK_SIMPLE_VERTEX);
 344 //
 345 //       //part 2
 346 //       double total_dist=0;
 347 //       double avg_dist=0;
 348 //       for(int i=0;i<N_pts;i++)
 349 //       {
 350 //         double dist=(corner[i]-C).abs();
 351 //         total_dist+=dist;
 352 //         node_meshdensity_current->SetValue(pts[i],NewCurrentMeshDensity(pts[i],dist));
 353 //       }
 354 //       avg_dist=total_dist/(double)N_pts;
 355 //       node_meshdensity_current->SetValue(l_newNodeId,1./avg_dist);
 356 //
 357 //       //part 3
 358 //       VertexMeshDensity nodeVMD;
 359 //       nodeVMD.type=node_type->GetValue(l_newNodeId);
 360 //       nodeVMD.density=0;
 361 //       nodeVMD.CurrentNode=l_newNodeId;
 362 //       EG_VTKDCC(vtkIntArray, cell_code, grid, "cell_code");
 363 //       nodeVMD.BCmap[cell_code->GetValue(id_cell)]=2;
 364 //
 365 //       int idx=VMDvector.indexOf(nodeVMD);
 366 //       node_specified_density->SetValue(l_newNodeId, idx);
 367 //
 368 //       //part 4
 369 //       if(idx!=-1)//specified
 370 //       {
 371 //         node_meshdensity_desired->SetValue(l_newNodeId, VMDvector[idx].density);
 372 //       }
 373 //       else//unspecified
 374 //       {
 375 //         double D=DesiredMeshDensity(l_newNodeId);
 376 //         node_meshdensity_desired->SetValue(l_newNodeId, D);
 377 //       }
 378
 379       //============================================