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 "engrid.h"
  24 #include "guimainwindow.h"
  25
  26 #include <vtkSignedCharArray.h>
  27
  28 #include <QTime>
  29
  30 InsertPoints::InsertPoints() : SurfaceOperation()
  31 {
  32   setQuickSave(true);
  33   getSet("surface meshing", "point insertion threshold", 1, m_Threshold);
  34 }
  35
  36 void InsertPoints::operate()
  37 {
  38   int N1 = m_Grid->GetNumberOfPoints();
  39   insertPoints();
  40   int N2 = m_Grid->GetNumberOfPoints();
  41   m_NumInserted = N2 - N1;
  42 }
  43
  44 ///\todo Adapt this code for multiple volumes.
  45 int InsertPoints::insertPoints()
  46 {
  47   setAllSurfaceCells();
  48   l2g_t  cells = getPartCells();
  49   g2l_t _cells = getPartLocalCells();
  50
  51   updateNodeInfo();
  52
  53   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
  54   EG_VTKDCN(vtkDoubleArray, characteristic_length_desired, m_Grid, "node_meshdensity_desired");
  55
  56   int num_newpoints=0;
  57   int num_newcells=0;
  58
  59   QVector <int> marked_cells(cells.size(), 0);
  60   QVector <stencil_t> stencil_vector(cells.size());
  61
  62   // find potential edges for splitting
  63   QList<edge_t> edges;
  64   for (int i = 0; i < cells.size(); ++i) {
  65     vtkIdType id_cell = cells[i];
  66
  67     // if cell is selected and a triangle
  68     if (m_BoundaryCodes.contains(cell_code->GetValue(id_cell)) && (m_Grid->GetCellType(id_cell) == VTK_TRIANGLE)) {
  69       int j_split = -1;
  70       double L_max = 0;
  71       EG_GET_CELL(id_cell, m_Grid);
  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)%num_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_t, 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         EG_GET_CELL(S.id_cell[i_triangle], grid_tmp);
 172
 173         bool direct;
 174         for(int i_pts = 0; i_pts < num_pts; i_pts++) {
 175           if (pts[i_pts] == S.p1 ) {
 176             if (pts[(i_pts+1)%num_pts] == S.p2 ) {
 177               direct = true;
 178             } else {
 179               direct = false;
 180             }
 181           }
 182         }
 183         if (direct) {
 184           pts_triangle[i_triangle][0] = S.p1;
 185           pts_triangle[i_triangle][1] = id_new_node;
 186           pts_triangle[i_triangle][2] = S.id_node[i_triangle];
 187
 188           pts_triangle[i_triangle+N][0] = id_new_node;
 189           pts_triangle[i_triangle+N][1] = S.p2;
 190           pts_triangle[i_triangle+N][2] = S.id_node[i_triangle];
 191         } else {
 192           pts_triangle[i_triangle][0] = S.p2;
 193           pts_triangle[i_triangle][1] = id_new_node;
 194           pts_triangle[i_triangle][2] = S.id_node[i_triangle];
 195
 196           pts_triangle[i_triangle+N][0] = id_new_node;
 197           pts_triangle[i_triangle+N][1] = S.p1;
 198           pts_triangle[i_triangle+N][2] = S.id_node[i_triangle];
 199         }
 200
 201         grid_tmp->ReplaceCell(S.id_cell[i_triangle] , 3, pts_triangle[i_triangle].data());
 202         auto new_pts = idListFromVector(pts_triangle[i_triangle+N]);
 203         vtkIdType newCellId = grid_tmp->InsertNextCell(VTK_TRIANGLE, new_pts);
 204         //vtkIdType newCellId = grid_tmp->InsertNextCell(VTK_TRIANGLE,3,pts_triangle[i_triangle+N].data());
 205         copyCellData(grid_tmp,S.id_cell[i_triangle],grid_tmp,newCellId);
 206       }
 207
 208       //increment ID
 209       id_new_node++;
 210     }
 211   }
 212
 213   //update grid
 214   makeCopy(grid_tmp,m_Grid);
 215
 216   return(0);
 217 }
 218
 219 char InsertPoints::getNewNodeType(stencil_t S)
 220 {
 221   vtkIdType id_node1 = S.p1;
 222   vtkIdType id_node2 = S.p2;
 223
 224   EG_VTKDCN(vtkCharArray_t, node_type, m_Grid, "node_type");
 225   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
 226   char type1 = node_type->GetValue(id_node1);
 227   char type2 = node_type->GetValue(id_node2);
 228   if (type1 == EG_SIMPLE_VERTEX || type2 == EG_SIMPLE_VERTEX) {
 229     return EG_SIMPLE_VERTEX;
 230   } else {
 231     QVector <vtkIdType> PSP1 = getPotentialSnapPoints(id_node1);
 232     QVector <vtkIdType> PSP2 = getPotentialSnapPoints(id_node2);
 233     bool complex = false;
 234     if (PSP1.contains(id_node2) || PSP2.contains(id_node1)) {
 235       complex = true;
 236     } else if (type1 >= EG_BOUNDARY_EDGE_VERTEX && type2 >= EG_BOUNDARY_EDGE_VERTEX) {
 237       if (!PSP1.contains(id_node2) && !PSP2.contains(id_node1) && S.id_cell.size() == 2) {
 238         int bc1 = cell_code->GetValue(S.id_cell[0]);
 239         int bc2 = cell_code->GetValue(S.id_cell[1]);
 240         if (bc1 != bc2) {
 241           return EG_BOUNDARY_EDGE_VERTEX;
 242         }
 243         vec3_t n1 = cellNormal(m_Grid, S.id_cell[0]);
 244         vec3_t n2 = cellNormal(m_Grid, S.id_cell[1]);
 245         if (GeometryTools::angle(n1, n2) > m_FeatureAngle) {
 246           return EG_FEATURE_EDGE_VERTEX;
 247         }
 248       }
 249     }
 250     if (complex) {
 251       if (S.id_cell.size() < 1) {
 252         return EG_BOUNDARY_EDGE_VERTEX;
 253       } else if (S.id_cell.size() == 1) {
 254         EG_ERR_RETURN("Invalid surface mesh. Check this with 'Tools -> Check surface integrity'.");
 255       } else {
 256         if (cell_code->GetValue(S.id_cell[0]) != cell_code->GetValue(S.id_cell[1])) {
 257           return EG_BOUNDARY_EDGE_VERTEX;
 258         } else {
 259           if (isFeatureNode(id_node1) || isFeatureNode(id_node2)) {
 260             return EG_FEATURE_EDGE_VERTEX;
 261           } else {
 262             return EG_SIMPLE_VERTEX;
 263           }
 264         }
 265       }
 266     } else {
 267       return EG_SIMPLE_VERTEX;
 268     }
 269   }
 270 }
 271