src/libengrid/insertpoints.cpp

   1 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   2 // +                                                                      +
   3 // + This file is part of enGrid.                                         +
   4 // +                                                                      +
   5 // + Copyright 2008-2014 enGits GmbH                                      +
   6 // +                                                                      +
   7 // + enGrid is free software: you can redistribute it and/or modify       +
   8 // + it under the terms of the GNU General Public License as published by +
   9 // + the Free Software Foundation, either version 3 of the License, or    +
  10 // + (at your option) any later version.                                  +
  11 // +                                                                      +
  12 // + enGrid is distributed in the hope that it will be useful,            +
  13 // + but WITHOUT ANY WARRANTY; without even the implied warranty of       +
  14 // + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the        +
  15 // + GNU General Public License for more details.                         +
  16 // +                                                                      +
  17 // + You should have received a copy of the GNU General Public License    +
  18 // + along with enGrid. If not, see <http://www.gnu.org/licenses/>.       +
  19 // +                                                                      +
  20 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  21 #include "insertpoints.h"
  22
  23 #include "guimainwindow.h"
  24
  25 #include <vtkCharArray.h>
  26
  27 #include <QTime>
  28
  29 InsertPoints::InsertPoints() : SurfaceOperation()
  30 {
  31   setQuickSave(true);
  32   getSet("surface meshing", "point insertion threshold", 1, m_Threshold);
  33 }
  34
  35 void InsertPoints::operate()
  36 {
  37   int N1 = m_Grid->GetNumberOfPoints();
  38   insertPoints();
  39   int N2 = m_Grid->GetNumberOfPoints();
  40   m_NumInserted = N2 - N1;
  41 }
  42
  43 ///\todo Adapt this code for multiple volumes.
  44 int InsertPoints::insertPoints()
  45 {
  46   setAllSurfaceCells();
  47   l2g_t  cells = getPartCells();
  48   g2l_t _cells = getPartLocalCells();
  49
  50   updateNodeInfo();
  51
  52   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
  53   EG_VTKDCN(vtkDoubleArray, characteristic_length_desired, m_Grid, "node_meshdensity_desired");
  54
  55   int num_newpoints=0;
  56   int num_newcells=0;
  57
  58   QVector <int> marked_cells(cells.size(), 0);
  59   QVector <stencil_t> stencil_vector(cells.size());
  60
  61   // find potential edges for splitting
  62   QList<edge_t> edges;
  63   for (int i = 0; i < cells.size(); ++i) {
  64     vtkIdType id_cell = cells[i];
  65
  66     // if cell is selected and a triangle
  67     if (m_BoundaryCodes.contains(cell_code->GetValue(id_cell)) && (m_Grid->GetCellType(id_cell) == VTK_TRIANGLE)) {
  68       int j_split = -1;
  69       double L_max = 0;
  70       vtkIdType N_pts, *pts;
  71       m_Grid->GetCellPoints(id_cell, N_pts, pts);
  72
  73       //find best side to split (longest)
  74       for (int j = 0; j < 3; ++j) {
  75
  76         //check if neighbour cell on this side is also selected
  77         stencil_t S = getStencil(id_cell, j);
  78         bool selected_edge = true;
  79         for (int i_cell_neighbour=1;i_cell_neighbour<S.id_cell.size();i_cell_neighbour++) {
  80           vtkIdType id_cell_neighbour = S.id_cell[i_cell_neighbour];
  81           if (!m_BoundaryCodes.contains(cell_code->GetValue(id_cell_neighbour)) || S.type_cell[i_cell_neighbour] != VTK_TRIANGLE) {
  82             selected_edge=false;
  83           }
  84         } // end of loop through neighbour cells
  85
  86         if (selected_edge) {
  87           vtkIdType id_node1 = pts[j];
  88           vtkIdType id_node2 = pts[(j+1)%N_pts];
  89           double L  = distance(m_Grid, id_node1, id_node2);
  90           double L1 = characteristic_length_desired->GetValue(id_node1);
  91           double L2 = characteristic_length_desired->GetValue(id_node2);
  92           if (L > m_Threshold*min(L1,L2)) {
  93             if (L > L_max) {
  94               j_split = j;
  95               L_max = L;
  96             }
  97           }
  98         }
  99       } // end of loop through sides
 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     int i_cells1 = _cells[E.S.id_cell[0]];
 117     bool all_unmarked = true;
 118     for (int i = 0; i < E.S.id_cell.size(); i++) {
 119       int i_cells = _cells[E.S.id_cell[i]];
 120       if (marked_cells[i_cells]) {
 121         all_unmarked=false;
 122       }
 123     }
 124     if (all_unmarked) {
 125       stencil_vector[i_cells1] = E.S;
 126       marked_cells[i_cells1] = 1;
 127       for (int i = 1; i < E.S.id_cell.size(); i++) {
 128         int i_cells = _cells[E.S.id_cell[i]];
 129         marked_cells[i_cells] = 2;
 130       }
 131       ++num_newpoints;
 132       num_newcells += E.S.id_cell.size();
 133     }
 134   }
 135
 136   //initialize grid_tmp
 137   int l_N_points=m_Grid->GetNumberOfPoints();
 138   int l_N_cells=m_Grid->GetNumberOfCells();
 139   EG_VTKSP(vtkUnstructuredGrid,grid_tmp);
 140   allocateGrid(grid_tmp, l_N_cells + num_newcells, l_N_points + num_newpoints);
 141   makeCopyNoAlloc(m_Grid, grid_tmp);
 142
 143   //initialize new node counter
 144   vtkIdType id_new_node = l_N_points;
 145
 146   //actor
 147   for (int i = 0; i < cells.size(); i++) {
 148     if (marked_cells[i] == 1) {
 149       stencil_t S = stencil_vector[i];
 150
 151       //calculate midpoint
 152       vec3_t A,B;
 153       grid_tmp->GetPoint(S.p1,A.data());
 154       grid_tmp->GetPoint(S.p2,B.data());
 155       vec3_t M=0.5*(A+B);
 156
 157       //add point
 158       grid_tmp->GetPoints()->SetPoint(id_new_node, M.data());
 159       copyNodeData(grid_tmp,S.p1,grid_tmp,id_new_node); ///\todo maybe trouble
 160
 161       // inserted edge point = type of the edge on which it is inserted
 162       EG_VTKDCN(vtkCharArray, node_type, grid_tmp, "node_type");
 163       node_type->SetValue(id_new_node, getNewNodeType(S) );
 164
 165       // insert new cells
 166       //four new triangles
 167       int N = S.id_cell.size();
 168       QVector< QVector<vtkIdType> > pts_triangle(2*N,QVector<vtkIdType>(3));
 169
 170       for(int i_triangle = 0; i_triangle < N; i_triangle++) {
 171         vtkIdType *pts, N_pts;
 172         grid_tmp->GetCellPoints(S.id_cell[i_triangle], N_pts, pts);
 173
 174         bool direct;
 175         for(int i_pts = 0; i_pts < N_pts; i_pts++) {
 176           if (pts[i_pts] == S.p1 ) {
 177             if (pts[(i_pts+1)%N_pts] == S.p2 ) {
 178               direct = true;
 179             } else {
 180               direct = false;
 181             }
 182           }
 183         }
 184         if (direct) {
 185           pts_triangle[i_triangle][0] = S.p1;
 186           pts_triangle[i_triangle][1] = id_new_node;
 187           pts_triangle[i_triangle][2] = S.id_node[i_triangle];
 188
 189           pts_triangle[i_triangle+N][0] = id_new_node;
 190           pts_triangle[i_triangle+N][1] = S.p2;
 191           pts_triangle[i_triangle+N][2] = S.id_node[i_triangle];
 192         } else {
 193           pts_triangle[i_triangle][0] = S.p2;
 194           pts_triangle[i_triangle][1] = id_new_node;
 195           pts_triangle[i_triangle][2] = S.id_node[i_triangle];
 196
 197           pts_triangle[i_triangle+N][0] = id_new_node;
 198           pts_triangle[i_triangle+N][1] = S.p1;
 199           pts_triangle[i_triangle+N][2] = S.id_node[i_triangle];
 200         }
 201
 202         grid_tmp->ReplaceCell(S.id_cell[i_triangle] , 3, pts_triangle[i_triangle].data());
 203         vtkIdType newCellId = grid_tmp->InsertNextCell(VTK_TRIANGLE,3,pts_triangle[i_triangle+N].data());
 204         copyCellData(grid_tmp,S.id_cell[i_triangle],grid_tmp,newCellId);
 205       }
 206
 207       //increment ID
 208       id_new_node++;
 209     }
 210   }
 211
 212   //update grid
 213   makeCopy(grid_tmp,m_Grid);
 214
 215   return(0);
 216 }
 217
 218 char InsertPoints::getNewNodeType(stencil_t S)
 219 {
 220   vtkIdType id_node1 = S.p1;
 221   vtkIdType id_node2 = S.p2;
 222
 223   EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type");
 224   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
 225   char type1 = node_type->GetValue(id_node1);
 226   char type2 = node_type->GetValue(id_node2);
 227   if (type1 == EG_SIMPLE_VERTEX || type2 == EG_SIMPLE_VERTEX) {
 228     return EG_SIMPLE_VERTEX;
 229   } else {
 230     QVector <vtkIdType> PSP1 = getPotentialSnapPoints(id_node1);
 231     QVector <vtkIdType> PSP2 = getPotentialSnapPoints(id_node2);
 232     bool complex = false;
 233     if (PSP1.contains(id_node2) || PSP2.contains(id_node1)) {
 234       complex = true;
 235     } else if (type1 >= EG_BOUNDARY_EDGE_VERTEX && type2 >= EG_BOUNDARY_EDGE_VERTEX) {
 236       if (!PSP1.contains(id_node2) && !PSP2.contains(id_node1) && S.id_cell.size() == 2) {
 237         int bc1 = cell_code->GetValue(S.id_cell[0]);
 238         int bc2 = cell_code->GetValue(S.id_cell[1]);
 239         if (bc1 != bc2) {
 240           return EG_BOUNDARY_EDGE_VERTEX;
 241         }
 242         vec3_t n1 = cellNormal(m_Grid, S.id_cell[0]);
 243         vec3_t n2 = cellNormal(m_Grid, S.id_cell[1]);
 244         if (GeometryTools::angle(n1, n2) > m_FeatureAngle) {
 245           return EG_FEATURE_EDGE_VERTEX;
 246         }
 247       }
 248     }
 249     if (complex) {
 250       if (S.id_cell.size() < 1) {
 251         return EG_BOUNDARY_EDGE_VERTEX;
 252       } else if (S.id_cell.size() == 1) {
 253         EG_ERR_RETURN("Invalid surface mesh. Check this with 'Tools -> Check surface integrity'.");
 254       } else {
 255         if (cell_code->GetValue(S.id_cell[0]) != cell_code->GetValue(S.id_cell[1])) {
 256           return EG_BOUNDARY_EDGE_VERTEX;
 257         } else {
 258           if (isFeatureNode(id_node1) || isFeatureNode(id_node2)) {
 259             return EG_FEATURE_EDGE_VERTEX;
 260           } else {
 261             return EG_SIMPLE_VERTEX;
 262           }
 263         }
 264       }
 265     } else {
 266       return EG_SIMPLE_VERTEX;
 267     }
 268   }
 269 }
 270