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() : SurfaceOperation()
  32 {
  33   setQuickSave(true);
  34   getSet("surface meshing", "point insertion threshold", 1, m_Threshold);
  35 }
  36
  37 void InsertPoints::operate()
  38 {
  39   int N1 = m_Grid->GetNumberOfPoints();
  40   insertPoints();
  41   int N2 = m_Grid->GetNumberOfPoints();
  42   m_NumInserted = N2 - N1;
  43 }
  44
  45 ///\todo Adapt this code for multiple volumes.
  46 int InsertPoints::insertPoints()
  47 {
  48   QTime start = QTime::currentTime();
  49
  50   setAllSurfaceCells();
  51   l2g_t  cells = getPartCells();
  52   g2l_t _cells = getPartLocalCells();
  53   g2l_t _nodes = getPartLocalNodes();
  54   l2l_t  n2c   = getPartN2C();
  55   l2l_t  c2c = getPartC2C();
  56
  57   UpdatePotentialSnapPoints(true);
  58
  59   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
  60   EG_VTKDCN(vtkDoubleArray, characteristic_length_desired, m_Grid, "node_meshdensity_desired");
  61
  62   int num_newpoints=0;
  63   int num_newcells=0;
  64
  65   QVector <int> marked_cells(cells.size(), 0);
  66   QVector <stencil_t> stencil_vector(cells.size());
  67
  68   // find potential edges for splitting
  69   QList<edge_t> edges;
  70   for (int i = 0; i < cells.size(); ++i) {
  71     vtkIdType id_cell = cells[i];
  72
  73     // if cell is selected and a triangle
  74     if (m_BoundaryCodes.contains(cell_code->GetValue(id_cell)) && (m_Grid->GetCellType(id_cell) == VTK_TRIANGLE)) {
  75       int j_split = -1;
  76       double L_max = 0;
  77       vtkIdType N_pts, *pts;
  78       m_Grid->GetCellPoints(id_cell, N_pts, pts);
  79       //find best side to split (longest)
  80       for (int j = 0; j < 3; ++j) {
  81         //check if neighbour cell on this side is also selected
  82         int i_cell_neighbour = c2c[_cells[id_cell]][j];
  83         if(i_cell_neighbour!=-1) {
  84           vtkIdType id_cell_neighbour = cells[i_cell_neighbour];
  85           if(m_BoundaryCodes.contains(cell_code->GetValue(id_cell_neighbour))) {
  86             vtkIdType id_node1 = pts[j];
  87             vtkIdType id_node2 = pts[(j+1)%N_pts];
  88             double L  = distance(m_Grid, id_node1, id_node2);
  89             double L1 = characteristic_length_desired->GetValue(id_node1);
  90             double L2 = characteristic_length_desired->GetValue(id_node2);
  91             if (L > m_Threshold*min(L1,L2)) {
  92               if (L > L_max) {
  93                 j_split = j;
  94                 L_max = L;
  95               }
  96             }
  97           }
  98         }
  99       }
 100       if (j_split != -1) {
 101         stencil_t S = getStencil(id_cell, j_split);
 102         edge_t E;
 103         E.S = S;
 104         E.L1 = characteristic_length_desired->GetValue(S.p1);
 105         E.L2 = characteristic_length_desired->GetValue(S.p2);
 106         E.L12 = distance(m_Grid, S.p1, S.p2);
 107         edges.push_back(E);
 108       }
 109     }
 110   }
 111
 112   qSort(edges);
 113
 114   //counter
 115   foreach (edge_t E, edges) {
 116     if (E.S.id_cell.size()==2 && (E.S.type_cell[1] == VTK_TRIANGLE)) {
 117       int i_cells1 = _cells[E.S.id_cell[0]];
 118       int i_cells2 = _cells[E.S.id_cell[1]];
 119       if (!marked_cells[i_cells1]  && !marked_cells[i_cells2]) {
 120         stencil_vector[i_cells1] = E.S;
 121         marked_cells[i_cells1] = 1;
 122         marked_cells[i_cells2] = 2;
 123         ++num_newpoints;
 124         num_newcells += 2;
 125       }
 126     } else if (E.S.id_cell.size()!=2) {
 127       int i_cells1 = _cells[E.S.id_cell[0]];
 128       if (!marked_cells[i_cells1]) {
 129         stencil_vector[i_cells1] = E.S;
 130         marked_cells[i_cells1] = 1;
 131         ++num_newpoints;
 132         num_newcells += 1;
 133       }
 134     }
 135   }
 136
 137   //initialize grid_tmp
 138   int l_N_points=m_Grid->GetNumberOfPoints();
 139   int l_N_cells=m_Grid->GetNumberOfCells();
 140   EG_VTKSP(vtkUnstructuredGrid,grid_tmp);
 141   allocateGrid(grid_tmp, l_N_cells + num_newcells, l_N_points + num_newpoints);
 142   makeCopyNoAlloc(m_Grid, grid_tmp);
 143
 144   //initialize new node counter
 145   vtkIdType id_new_node = l_N_points;
 146
 147   //actor
 148   for (int i = 0; i < cells.size(); i++) {
 149     if (marked_cells[i] == 1) {
 150       stencil_t S = stencil_vector[i];
 151
 152       //calculate midpoint
 153       vec3_t A,B;
 154       grid_tmp->GetPoint(S.p1,A.data());
 155       grid_tmp->GetPoint(S.p2,B.data());
 156       vec3_t M=0.5*(A+B);
 157
 158       //add point
 159       grid_tmp->GetPoints()->SetPoint(id_new_node, M.data());
 160       copyNodeData(grid_tmp,S.p1,grid_tmp,id_new_node); ///\todo maybe trouble
 161
 162       // inserted edge point = type of the edge on which it is inserted
 163       EG_VTKDCN(vtkCharArray, node_type, grid_tmp, "node_type");
 164       node_type->SetValue(id_new_node, getNewNodeType(S) );
 165
 166       if(S.id_cell.size()==2 && S.type_cell[1]==VTK_TRIANGLE) { //2 triangles
 167         //four new triangles
 168         vtkIdType pts_triangle[4][3];
 169
 170         pts_triangle[0][0] = S.id_node[0];
 171         pts_triangle[0][1] = S.p1;
 172         pts_triangle[0][2] = id_new_node;
 173
 174         pts_triangle[1][0] = S.p1;
 175         pts_triangle[1][1] = S.id_node[1];
 176         pts_triangle[1][2] = id_new_node;
 177
 178         pts_triangle[2][0] = S.id_node[1];
 179         pts_triangle[2][1] = S.p2;
 180         pts_triangle[2][2] = id_new_node;
 181
 182         pts_triangle[3][0] = S.p2;
 183         pts_triangle[3][1] = S.id_node[0];
 184         pts_triangle[3][2] = id_new_node;
 185
 186         grid_tmp->ReplaceCell(S.id_cell[0] , 3, pts_triangle[0]);
 187         grid_tmp->ReplaceCell(S.id_cell[1] , 3, pts_triangle[1]);
 188
 189         vtkIdType newCellId;
 190
 191         newCellId = grid_tmp->InsertNextCell(VTK_TRIANGLE,3,pts_triangle[2]);
 192         copyCellData(grid_tmp,S.id_cell[1],grid_tmp,newCellId);
 193
 194         newCellId = grid_tmp->InsertNextCell(VTK_TRIANGLE,3,pts_triangle[3]);
 195         copyCellData(grid_tmp,S.id_cell[0],grid_tmp,newCellId);
 196       } else if(S.id_cell.size()!=2) { //1 triangle
 197         //two new triangles
 198         vtkIdType pts_triangle[2][3];
 199         pts_triangle[0][0] = S.id_node[0];
 200         pts_triangle[0][1] = S.p1;
 201         pts_triangle[0][2] = id_new_node;
 202         pts_triangle[1][0] = S.p2;
 203         pts_triangle[1][1] = S.id_node[0];
 204         pts_triangle[1][2] = id_new_node;
 205
 206         grid_tmp->ReplaceCell(S.id_cell[0] , 3, pts_triangle[0]);
 207
 208         vtkIdType newCellId;
 209         newCellId = grid_tmp->InsertNextCell(VTK_TRIANGLE,3,pts_triangle[1]);
 210         copyCellData(grid_tmp,S.id_cell[0],grid_tmp,newCellId);
 211       } else {
 212         cout<<"I DON'T KNOW HOW TO SPLIT THIS CELL!!!"<<endl;
 213         EG_BUG;
 214       }
 215
 216       //increment ID
 217       id_new_node++;
 218     }
 219   }
 220
 221   //update grid
 222   makeCopy(grid_tmp,m_Grid);
 223
 224   return(0);
 225 }
 226
 227 char InsertPoints::getNewNodeType(stencil_t S)
 228 {
 229   vtkIdType id_node1 = S.p1;
 230   vtkIdType id_node2 = S.p2;
 231
 232   EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type");
 233   if( node_type->GetValue(id_node1)==VTK_SIMPLE_VERTEX || node_type->GetValue(id_node2)==VTK_SIMPLE_VERTEX ) {
 234     return VTK_SIMPLE_VERTEX;
 235   }
 236   else {
 237     QVector <vtkIdType> PSP = getPotentialSnapPoints(id_node1);
 238     if( PSP.contains(id_node2) ) {
 239       EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
 240       if( cell_code->GetValue(S.id_cell[0]) != cell_code->GetValue(S.id_cell[1]) ) {
 241         return VTK_BOUNDARY_EDGE_VERTEX;
 242       } else {
 243         return VTK_FEATURE_EDGE_VERTEX;
 244       }
 245     } else {
 246       return VTK_SIMPLE_VERTEX;
 247     }
 248   }
 249 }
 250