src/libengrid/insertpoints.cpp

   1 //
   2 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   3 // +                                                                      +
   4 // + This file is part of enGrid.                                         +
   5 // +                                                                      +
   6 // + Copyright 2008-2010 enGits GmbH                                     +
   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
  54   UpdatePotentialSnapPoints(true);
  55
  56   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
  57   EG_VTKDCN(vtkDoubleArray, characteristic_length_desired, m_Grid, "node_meshdensity_desired");
  58
  59   int num_newpoints=0;
  60   int num_newcells=0;
  61
  62   QVector <int> marked_cells(cells.size(), 0);
  63   QVector <stencil_t> stencil_vector(cells.size());
  64
  65   // find potential edges for splitting
  66   QList<edge_t> edges;
  67   for (int i = 0; i < cells.size(); ++i) {
  68     vtkIdType id_cell = cells[i];
  69
  70     // if cell is selected and a triangle
  71     if (m_BoundaryCodes.contains(cell_code->GetValue(id_cell)) && (m_Grid->GetCellType(id_cell) == VTK_TRIANGLE)) {
  72       int j_split = -1;
  73       double L_max = 0;
  74       vtkIdType N_pts, *pts;
  75       m_Grid->GetCellPoints(id_cell, N_pts, pts);
  76       //find best side to split (longest)
  77       for (int j = 0; j < 3; ++j) {
  78         //check if neighbour cell on this side is also selected
  79         stencil_t S = getStencil(id_cell, j);
  80         bool selected_edge = true;
  81         for(int i_cell_neighbour=1;i_cell_neighbour<S.id_cell.size();i_cell_neighbour++) {
  82           vtkIdType id_cell_neighbour = S.id_cell[i_cell_neighbour];
  83           if( !m_BoundaryCodes.contains(cell_code->GetValue(id_cell_neighbour)) || S.type_cell[i_cell_neighbour] != VTK_TRIANGLE) selected_edge=false;
  84         }// end of loop through neighbour cells
  85         if(selected_edge) {
  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       }// end of loop through sides
  99       if (j_split != -1) {
 100         stencil_t S = getStencil(id_cell, j_split);
 101         edge_t E;
 102         E.S = S;
 103         E.L1 = characteristic_length_desired->GetValue(S.p1);
 104         E.L2 = characteristic_length_desired->GetValue(S.p2);
 105         E.L12 = distance(m_Grid, S.p1, S.p2);
 106         edges.push_back(E);
 107       }
 108     }
 109   }
 110
 111   qSort(edges);
 112
 113   //counter
 114   foreach (edge_t E, edges) {
 115     int i_cells1 = _cells[E.S.id_cell[0]];
 116     bool all_unmarked = true;
 117     for(int i=0; i<E.S.id_cell.size();i++) {
 118       int i_cells = _cells[E.S.id_cell[i]];
 119       if (marked_cells[i_cells]) all_unmarked=false;
 120     }
 121     if(all_unmarked) {
 122       stencil_vector[i_cells1] = E.S;
 123       marked_cells[i_cells1] = 1;
 124       for(int i=1; i<E.S.id_cell.size();i++) {
 125         int i_cells = _cells[E.S.id_cell[i]];
 126         marked_cells[i_cells] = 2;
 127       }
 128       ++num_newpoints;
 129       num_newcells += E.S.id_cell.size();
 130     }
 131   }
 132
 133   //initialize grid_tmp
 134   int l_N_points=m_Grid->GetNumberOfPoints();
 135   int l_N_cells=m_Grid->GetNumberOfCells();
 136   EG_VTKSP(vtkUnstructuredGrid,grid_tmp);
 137   allocateGrid(grid_tmp, l_N_cells + num_newcells, l_N_points + num_newpoints);
 138   makeCopyNoAlloc(m_Grid, grid_tmp);
 139
 140   //initialize new node counter
 141   vtkIdType id_new_node = l_N_points;
 142
 143   //actor
 144   for (int i = 0; i < cells.size(); i++) {
 145     if (marked_cells[i] == 1) {
 146       stencil_t S = stencil_vector[i];
 147
 148       //calculate midpoint
 149       vec3_t A,B;
 150       grid_tmp->GetPoint(S.p1,A.data());
 151       grid_tmp->GetPoint(S.p2,B.data());
 152       vec3_t M=0.5*(A+B);
 153
 154       //add point
 155       grid_tmp->GetPoints()->SetPoint(id_new_node, M.data());
 156       copyNodeData(grid_tmp,S.p1,grid_tmp,id_new_node); ///\todo maybe trouble
 157
 158       // inserted edge point = type of the edge on which it is inserted
 159       EG_VTKDCN(vtkCharArray, node_type, grid_tmp, "node_type");
 160       node_type->SetValue(id_new_node, getNewNodeType(S) );
 161
 162       // insert new cells
 163       //four new triangles
 164       int N = S.id_cell.size();
 165       vtkIdType pts_triangle[2*N][3];
 166
 167       for(int i_triangle=0; i_triangle<N; i_triangle++) {
 168         vtkIdType *pts, N_pts;
 169         grid_tmp->GetCellPoints(S.id_cell[i_triangle], N_pts, pts);
 170
 171         bool direct;
 172         for(int i_pts = 0; i_pts<N_pts; i_pts++) {
 173           if( pts[i_pts] == S.p1 ) {
 174             if( pts[(i_pts+1)%N_pts] == S.p2 ) direct = true;
 175             else direct = false;
 176           }
 177         }
 178         if(direct) {
 179           pts_triangle[i_triangle][0] = S.p1;
 180           pts_triangle[i_triangle][1] = id_new_node;
 181           pts_triangle[i_triangle][2] = S.id_node[i_triangle];
 182
 183           pts_triangle[i_triangle+N][0] = id_new_node;
 184           pts_triangle[i_triangle+N][1] = S.p2;
 185           pts_triangle[i_triangle+N][2] = S.id_node[i_triangle];
 186         }
 187         else {
 188           pts_triangle[i_triangle][0] = S.p2;
 189           pts_triangle[i_triangle][1] = id_new_node;
 190           pts_triangle[i_triangle][2] = S.id_node[i_triangle];
 191
 192           pts_triangle[i_triangle+N][0] = id_new_node;
 193           pts_triangle[i_triangle+N][1] = S.p1;
 194           pts_triangle[i_triangle+N][2] = S.id_node[i_triangle];
 195         }
 196
 197         grid_tmp->ReplaceCell(S.id_cell[i_triangle] , 3, pts_triangle[i_triangle]);
 198         vtkIdType newCellId = grid_tmp->InsertNextCell(VTK_TRIANGLE,3,pts_triangle[i_triangle+N]);
 199         copyCellData(grid_tmp,S.id_cell[i_triangle],grid_tmp,newCellId);
 200       }
 201
 202       //increment ID
 203       id_new_node++;
 204     }
 205   }
 206
 207   //update grid
 208   makeCopy(grid_tmp,m_Grid);
 209
 210   return(0);
 211 }
 212
 213 char InsertPoints::getNewNodeType(stencil_t S)
 214 {
 215 //   cout<<"S="<<S<<endl;
 216
 217   vtkIdType id_node1 = S.p1;
 218   vtkIdType id_node2 = S.p2;
 219
 220   if (S.id_cell.size() != 2) {
 221     EG_ERR_RETURN("The surface mesh is not water-tight");
 222   }
 223
 224   EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type");
 225   if( node_type->GetValue(id_node1)==VTK_SIMPLE_VERTEX || node_type->GetValue(id_node2)==VTK_SIMPLE_VERTEX ) {
 226     return VTK_SIMPLE_VERTEX;
 227   } else {
 228     QVector <vtkIdType> PSP = getPotentialSnapPoints(id_node1);
 229     if( PSP.contains(id_node2) ) {
 230       EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
 231       if(S.id_cell.size()<1) {
 232         return VTK_BOUNDARY_EDGE_VERTEX;
 233       } else {
 234         if( cell_code->GetValue(S.id_cell[0]) != cell_code->GetValue(S.id_cell[1]) ) {
 235           return VTK_BOUNDARY_EDGE_VERTEX;
 236         } else {
 237           return VTK_FEATURE_EDGE_VERTEX;
 238         }
 239       }
 240     } else {
 241       return VTK_SIMPLE_VERTEX;
 242     }
 243   }
 244 }
 245