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   bool result = distance(grid, id_node1, id_node2) > 0.5 * ( desiredEdgeLength(id_node1) + desiredEdgeLength(id_node2) );
  62   return ( result );
  63 }
  64
  65 bool InsertPoints::SplitSide(vtkIdType id_cell,int side)
  66 {
  67   vtkIdType N_pts,*pts;
  68   grid->GetCellPoints(id_cell,N_pts,pts);
  69   return( insert_edgepoint(pts[side],pts[(side+1)%N_pts]) );
  70 }
  71
  72 int InsertPoints::insert_FP_all()
  73 {
  74   //cout<<"===insert_FP_all START==="<<endl;
  75   QTime start = QTime::currentTime();
  76
  77   setAllSurfaceCells();
  78   UpdatePotentialSnapPoints(true);
  79
  80   QVector <vtkIdType> l_SelectedCells;
  81   getSurfaceCells(m_bcs, l_SelectedCells, grid);
  82
  83   QVector <bool> l_marked_cells(l_SelectedCells.size());
  84
  85   int l_N_newpoints=0;
  86   int l_N_newcells=0;
  87
  88   //counter
  89   for(int i_cell=0;i_cell<l_SelectedCells.size();i_cell++)
  90   {
  91     vtkIdType id_cell = l_SelectedCells[i_cell];
  92     if( insert_fieldpoint(id_cell) )
  93     {
  94       l_marked_cells[i_cell] = true;
  95       l_N_newcells += 2;
  96       l_N_newpoints += 1;
  97     }
  98   }
  99
 100   //initialize grid_tmp
 101   int l_N_points = grid->GetNumberOfPoints();
 102   int l_N_cells = grid->GetNumberOfCells();
 103   EG_VTKSP(vtkUnstructuredGrid,grid_tmp);
 104   allocateGrid(grid_tmp,l_N_cells+l_N_newcells,l_N_points+l_N_newpoints);
 105   makeCopyNoAlloc(grid, grid_tmp);
 106
 107   //initialize new node counter
 108   vtkIdType l_newNodeId = l_N_points;
 109
 110   //actor
 111   for(int i_cell=0;i_cell<l_SelectedCells.size();i_cell++)
 112   {
 113     vtkIdType id_cell = l_SelectedCells[i_cell];
 114     if( l_marked_cells[i_cell] )
 115     {
 116       vtkIdType N_pts, *pts;
 117       grid->GetCellPoints(id_cell, N_pts, pts);
 118       vec3_t C(0,0,0);
 119       for(int i=0;i<N_pts;i++)
 120       {
 121         vec3_t corner;
 122         grid->GetPoints()->GetPoint(pts[i], corner.data());
 123         C+=corner;
 124       }
 125       C=(1/(double)N_pts)*C;
 126
 127       //C=project(C);
 128       grid_tmp->GetPoints()->SetPoint(l_newNodeId,C.data());
 129       copyNodeData(grid_tmp,pts[0],grid_tmp,l_newNodeId);
 130       EG_VTKDCN(vtkCharArray, node_type, grid_tmp, "node_type");
 131       node_type->SetValue(l_newNodeId, VTK_SIMPLE_VERTEX);
 132
 133       for(int i=0;i<N_pts;i++)
 134       {
 135         vtkIdType pts_triangle[3];
 136         pts_triangle[0]=pts[i];
 137         pts_triangle[1]=pts[(i+1)%N_pts];
 138         pts_triangle[2]=l_newNodeId;
 139         if(i==0)
 140         {
 141           grid_tmp->ReplaceCell(id_cell , 3, pts_triangle);
 142         }
 143         else
 144         {
 145           vtkIdType newCellId = grid_tmp->InsertNextCell(VTK_TRIANGLE,3,pts_triangle);
 146           copyCellData(grid_tmp,id_cell,grid_tmp,newCellId);
 147         }
 148       }
 149       l_newNodeId++;
 150     }
 151   }
 152
 153   //update grid
 154   makeCopy(grid_tmp,grid);
 155
 156   //cout << start.msecsTo(QTime::currentTime()) << " milliseconds elapsed" << endl;
 157   //cout<<"===insert_FP_all END==="<<endl;
 158   return(0);
 159 }
 160
 161 int InsertPoints::insert_EP_all()
 162 {
 163   l2g_t  cells = getPartCells();
 164   g2l_t _cells = getPartLocalCells();
 165
 166   //cout<<"===insert_EP_all START==="<<endl;
 167   QTime start = QTime::currentTime();
 168
 169   setAllSurfaceCells();
 170   UpdatePotentialSnapPoints(true);
 171
 172   EG_VTKDCC(vtkIntArray, cell_code, grid, "cell_code");
 173
 174   int num_newpoints=0;
 175   int num_newcells=0;
 176
 177   QVector <int> marked_cells(cells.size());
 178   QVector <stencil_t> stencil_vector(cells.size());
 179
 180   //counter
 181   for (int i = 0; i < cells.size(); ++i) {
 182     vtkIdType id_cell = cells[i];
 183     if (m_bcs.contains(cell_code->GetValue(id_cell)) && (grid->GetCellType(id_cell) == VTK_TRIANGLE)) {//if selected and triangle cell
 184       int j_split = -1;
 185       double L_max = 0;
 186       vtkIdType N_pts, *pts;
 187       grid->GetCellPoints(id_cell, N_pts, pts);
 188       for (int j = 0; j < 3; ++j) {
 189         vtkIdType id_node1 = pts[j];
 190         vtkIdType id_node2 = pts[(j+1)%N_pts];
 191         double L = distance(grid, id_node1, id_node2);
 192         if (L > max(desiredEdgeLength(id_node1), desiredEdgeLength(id_node2))) {
 193           if (L > L_max) {
 194             j_split = j;
 195             L_max = L;
 196           }
 197         }
 198       }
 199       if (j_split != -1) {
 200         stencil_t S = getStencil(id_cell, j_split);
 201         if (S.twocells && (S.neighbour_type == VTK_TRIANGLE)) {
 202           if (!marked_cells[i]  && !marked_cells[_cells[S.id_cell2]]) {
 203             stencil_vector[i] = S;
 204             marked_cells[i] = 1;
 205             marked_cells[_cells[S.id_cell2]] = 2;
 206             num_newpoints++;
 207             num_newcells += 2;
 208           }
 209         } else if (!S.twocells) {
 210           if (!marked_cells[i]) {
 211             stencil_vector[i] = S;
 212             marked_cells[i] = 1;
 213             num_newpoints++;
 214             num_newcells+=1;
 215           }
 216         }
 217       } //end of loop through sides
 218     } //end of if selected and triangle cell
 219   } //end of counter loop
 220
 221   //initialize grid_tmp
 222   int l_N_points=grid->GetNumberOfPoints();
 223   int l_N_cells=grid->GetNumberOfCells();
 224   EG_VTKSP(vtkUnstructuredGrid,grid_tmp);
 225   allocateGrid(grid_tmp, l_N_cells + num_newcells, l_N_points + num_newpoints);
 226   makeCopyNoAlloc(grid, grid_tmp);
 227
 228   //initialize new node counter
 229   vtkIdType l_newNodeId = l_N_points;
 230
 231   //actor
 232   for (int i = 0; i < cells.size(); i++) {
 233     if (marked_cells[i] == 1) {
 234       stencil_t S = stencil_vector[i];
 235
 236       //calculate midpoint
 237       vec3_t A,B;
 238       grid_tmp->GetPoint(S.p[1],A.data());
 239       grid_tmp->GetPoint(S.p[3],B.data());
 240       vec3_t M=0.5*(A+B);
 241
 242       //project point
 243       //M=project(M);
 244       //add point
 245       grid_tmp->GetPoints()->SetPoint(l_newNodeId, M.data());
 246       copyNodeData(grid_tmp,S.p[1],grid_tmp,l_newNodeId);
 247
 248       // inserted edge point = type of the edge on which it is inserted
 249       EG_VTKDCN(vtkCharArray, node_type, grid_tmp, "node_type");
 250       node_type->SetValue(l_newNodeId, getNewNodeType(S) );
 251
 252       if(S.twocells && S.neighbour_type==VTK_TRIANGLE){//2 triangles
 253         //four new triangles
 254         vtkIdType pts_triangle[4][3];
 255         for(int i=0;i<4;i++)
 256         {
 257           pts_triangle[i][0]=S.p[i];
 258           pts_triangle[i][1]=S.p[(i+1)%4];
 259           pts_triangle[i][2]=l_newNodeId;
 260         }
 261
 262         grid_tmp->ReplaceCell(S.id_cell1 , 3, pts_triangle[0]);
 263         grid_tmp->ReplaceCell(S.id_cell2 , 3, pts_triangle[1]);
 264
 265         vtkIdType newCellId;
 266
 267         newCellId = grid_tmp->InsertNextCell(VTK_TRIANGLE,3,pts_triangle[2]);
 268         copyCellData(grid_tmp,S.id_cell2,grid_tmp,newCellId);
 269
 270         newCellId = grid_tmp->InsertNextCell(VTK_TRIANGLE,3,pts_triangle[3]);
 271         copyCellData(grid_tmp,S.id_cell1,grid_tmp,newCellId);
 272       }
 273       else if(!S.twocells) {//1 triangle
 274         //two new triangles
 275         vtkIdType pts_triangle[2][3];
 276         pts_triangle[0][0]=S.p[0];
 277         pts_triangle[0][1]=S.p[1];
 278         pts_triangle[0][2]=l_newNodeId;
 279         pts_triangle[1][0]=S.p[3];
 280         pts_triangle[1][1]=S.p[0];
 281         pts_triangle[1][2]=l_newNodeId;
 282
 283         grid_tmp->ReplaceCell(S.id_cell1 , 3, pts_triangle[0]);
 284
 285         vtkIdType newCellId;
 286         newCellId = grid_tmp->InsertNextCell(VTK_TRIANGLE,3,pts_triangle[1]);
 287         copyCellData(grid_tmp,S.id_cell1,grid_tmp,newCellId);
 288       }
 289       else {
 290         cout<<"I DON'T KNOW HOW TO SPLIT THIS CELL!!!"<<endl;
 291         EG_BUG;
 292       }
 293
 294       //increment ID
 295       l_newNodeId++;
 296     }
 297   }
 298
 299   //update grid
 300   makeCopy(grid_tmp,grid);
 301
 302   //cout << start.msecsTo(QTime::currentTime()) << " milliseconds elapsed" << endl;
 303   //cout<<"===insert_EP_all END==="<<endl;
 304   return(0);
 305 }
 306
 307 char InsertPoints::getNewNodeType(stencil_t S)
 308 {
 309   vtkIdType id_node1 = S.p[1];
 310   vtkIdType id_node2 = S.p[3];
 311
 312   EG_VTKDCN(vtkCharArray, node_type, grid, "node_type");
 313   if( node_type->GetValue(id_node1)==VTK_SIMPLE_VERTEX || node_type->GetValue(id_node2)==VTK_SIMPLE_VERTEX ) {
 314     return VTK_SIMPLE_VERTEX;
 315   }
 316   else {
 317     QVector <vtkIdType> PSP = getPotentialSnapPoints(id_node1);
 318     if( PSP.contains(id_node2) ) {
 319       EG_VTKDCC(vtkIntArray, cell_code, grid, "cell_code");
 320       if( cell_code->GetValue(S.id_cell1) != cell_code->GetValue(S.id_cell2) ) {
 321         return VTK_BOUNDARY_EDGE_VERTEX;
 322       }
 323       else {
 324         return VTK_FEATURE_EDGE_VERTEX;
 325       }
 326     }
 327     else {
 328       return VTK_SIMPLE_VERTEX;
 329     }
 330   }
 331 }
 332
 333        //============================================
 334       ///@@@  TODO: Update node info (densities+type) Still necessary?
 335 // EG_VTKDCN(vtkIntArray, node_specified_density, grid_tmp, "node_specified_density");//density index from table
 336 // EG_VTKDCN(vtkDoubleArray, node_meshdensity_desired, grid_tmp, "node_meshdensity_desired");//what we want
 337 // EG_VTKDCN(vtkDoubleArray, node_meshdensity_current, grid_tmp, "node_meshdensity_current");//what we have
 338 // EG_VTKDCN(vtkCharArray, node_type, grid_tmp, "node_type");//node type
 339       //============================================
 340
 341 //       //part 1
 342 //       node_type->SetValue(l_newNodeId,VTK_SIMPLE_VERTEX);
 343 //
 344 //       //part 2
 345 //       double total_dist=0;
 346 //       double avg_dist=0;
 347 //       for(int i=0;i<N_pts;i++)
 348 //       {
 349 //         double dist=(corner[i]-C).abs();
 350 //         total_dist+=dist;
 351 //         node_meshdensity_current->SetValue(pts[i],NewCurrentMeshDensity(pts[i],dist));
 352 //       }
 353 //       avg_dist=total_dist/(double)N_pts;
 354 //       node_meshdensity_current->SetValue(l_newNodeId,1./avg_dist);
 355 //
 356 //       //part 3
 357 //       VertexMeshDensity nodeVMD;
 358 //       nodeVMD.type=node_type->GetValue(l_newNodeId);
 359 //       nodeVMD.density=0;
 360 //       nodeVMD.CurrentNode=l_newNodeId;
 361 //       EG_VTKDCC(vtkIntArray, cell_code, grid, "cell_code");
 362 //       nodeVMD.BCmap[cell_code->GetValue(id_cell)]=2;
 363 //
 364 //       int idx=VMDvector.indexOf(nodeVMD);
 365 //       node_specified_density->SetValue(l_newNodeId, idx);
 366 //
 367 //       //part 4
 368 //       if(idx!=-1)//specified
 369 //       {
 370 //         node_meshdensity_desired->SetValue(l_newNodeId, VMDvector[idx].density);
 371 //       }
 372 //       else//unspecified
 373 //       {
 374 //         double D=DesiredMeshDensity(l_newNodeId);
 375 //         node_meshdensity_desired->SetValue(l_newNodeId, D);
 376 //       }
 377
 378       //============================================