src/libengrid/vtkEgEliminateShortEdges.cxx

   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 "vtkEgEliminateShortEdges.h"
  22 #include "engrid.h"
  23 #include "geometrytools.h"
  24
  25 vtkStandardNewMacro(vtkEgEliminateShortEdges);
  26
  27 vtkEgEliminateShortEdges::vtkEgEliminateShortEdges()
  28 {
  29   max_ratio = 1000.0;
  30   max_length = 1e99;
  31 };
  32
  33 void vtkEgEliminateShortEdges::CheckEdges()
  34 {
  35   new_node_idx.fill(-1, m_Input->GetNumberOfPoints());
  36   delete_cell.fill(false, m_Input->GetNumberOfCells());
  37   QVector<bool> marked(m_Input->GetNumberOfPoints(), false);
  38   for (vtkIdType cellId = 0; cellId < m_Input->GetNumberOfCells(); ++cellId) {
  39     int cellType = m_Input->GetCellType(cellId);
  40     if ((cellType == VTK_TRIANGLE) || (cellType == VTK_QUAD)) {
  41       EG_GET_CELL(cellId, m_Input);
  42       double L_av = 0;
  43       vector<vec3_t> x(num_pts);
  44       int N_del = 0;
  45       double Lmin = 1e99;
  46       int imin = 0;
  47       for (int i = 0; i < num_pts; ++i) {
  48         vtkIdType p1 = pts[i];
  49         vtkIdType p2;
  50         if (i < num_pts-1) p2 = pts[i+1];
  51         else            p2 = pts[0];
  52         if (p1 == p2) {
  53           EG_ERR_RETURN("bug encountered");
  54         };
  55         vec3_t x1, x2;
  56         m_Input->GetPoints()->GetPoint(p1, x1.data());
  57         m_Input->GetPoints()->GetPoint(p2, x2.data());
  58         double L = (x2-x1).abs();
  59         L_av += L;
  60         if (L < Lmin) {
  61           Lmin = L;
  62           imin = i;
  63         };
  64         x[i] = x1;
  65         if (marked[pts[i]]) {
  66           ++N_del;
  67         };
  68       };
  69       if (N_del == 0) {
  70         double A = GeometryTools::triArea(x[0],x[1],x[2]);
  71         if (num_pts == 4) {
  72           A = GeometryTools::quadArea(x[0],x[1],x[2],x[3]);
  73         };
  74         L_av /= num_pts;
  75         for (int i = 0; i < num_pts; ++i) {
  76           vtkIdType p1 = pts[i];
  77           vtkIdType p2;
  78           if (i < num_pts-1) p2 = pts[i+1];
  79           else            p2 = pts[0];
  80           vec3_t x1, x2;
  81           m_Input->GetPoints()->GetPoint(p1, x1.data());
  82           m_Input->GetPoints()->GetPoint(p2, x2.data());
  83           double L = (x2-x1).abs();
  84           bool delete_edge = false;
  85           if (L == 0) {
  86             delete_edge = true;
  87           } else if ((L_av*L_av/A > max_ratio) && (i == imin)) {
  88             if (L < max_length) {
  89               delete_edge = true;
  90             };
  91           };
  92           if (delete_edge) {
  93             new_node_idx[p1] = p2;
  94             marked[p1] = true;
  95             marked[p2] = true;
  96             if (cellType != VTK_TRIANGLE) {
  97               EG_ERR_RETURN("The present configuration cannot be handled yet.");
  98             };
  99             delete_cell[cellId] = true;
 100             break;
 101           };
 102         };
 103       };
 104     };
 105   };
 106 };
 107
 108 void vtkEgEliminateShortEdges::CheckCells()
 109 {
 110   for (vtkIdType cellId = 0; cellId < m_Input->GetNumberOfCells(); ++cellId) {
 111     int cellType = m_Input->GetCellType(cellId);
 112     if (cellType == VTK_TRIANGLE) {
 113       EG_GET_CELL(cellId, m_Input);
 114       for (int i = 0; i < num_pts; ++i) {
 115         vtkIdType p1 = pts[i];
 116         vtkIdType p2;
 117         if (i < num_pts-1) p2 = pts[i+1];
 118         else            p2 = pts[0];
 119         if (p1 == p2) {
 120           EG_ERR_RETURN("bug encountered");
 121         };
 122         if (new_node_idx[p1] == p2) {
 123           delete_cell[cellId] = true;
 124         };
 125         if (new_node_idx[p2] == p1) {
 126           delete_cell[cellId] = true;
 127         };
 128       };
 129     };
 130   };
 131 };
 132
 133 void vtkEgEliminateShortEdges::CopyPoints()
 134 {
 135   node_mapping.fill(-1, m_Input->GetNumberOfPoints());
 136   vtkIdType newPointId = 0;
 137   for (vtkIdType pointId = 0; pointId < m_Input->GetNumberOfPoints(); ++pointId) {
 138     if(new_node_idx[pointId] < 0) {
 139       node_mapping[pointId] = newPointId;
 140       vec3_t x;
 141       m_Input->GetPoints()->GetPoint(pointId,x.data());
 142       m_Output->GetPoints()->SetPoint(newPointId,x.data());
 143       ++newPointId;
 144     };
 145   };
 146   for (vtkIdType pointId = 0; pointId < m_Input->GetNumberOfPoints(); ++pointId) {
 147     if(new_node_idx[pointId] >= 0) {
 148       if (node_mapping[new_node_idx[pointId]] < 0) {
 149         EG_ERR_RETURN("bug encountered");
 150       };
 151       node_mapping[pointId] = node_mapping[new_node_idx[pointId]];
 152     };
 153   };
 154 };
 155
 156 void vtkEgEliminateShortEdges::CopyCells()
 157 {
 158   for (vtkIdType cellId = 0; cellId < m_Input->GetNumberOfCells(); ++cellId) {
 159     if(!delete_cell[cellId]) {
 160       EG_GET_CELL(cellId, m_Input);
 161       std::vector<vtkIdType> new_pts(num_pts);
 162       for (int i = 0; i < num_pts; ++i) {
 163         new_pts[i] = node_mapping[pts[i]];
 164       };
 165       auto newPtIds = idListFromVector(new_pts);
 166       vtkIdType newCellId = m_Output->InsertNextCell(m_Input->GetCellType(cellId), newPtIds);
 167       copyCellData(m_Input, cellId, m_Output, newCellId);
 168     };
 169   };
 170 };
 171
 172 void vtkEgEliminateShortEdges::ExecuteEg()
 173 {
 174   N_eliminated = 0;
 175   N_new_points = 0;
 176   N_new_cells  = 0;
 177   CheckEdges();
 178   CheckCells();
 179   int N = 0;
 180   for (vtkIdType pointId = 0; pointId < m_Input->GetNumberOfPoints(); ++pointId) {
 181     if(new_node_idx[pointId] < 0) {
 182       ++N_new_points;
 183     };
 184   };
 185   for (vtkIdType cellId = 0; cellId < m_Input->GetNumberOfCells(); ++cellId) {
 186     if(!delete_cell[cellId]) {
 187       ++N_new_cells;
 188     } else {
 189       ++N;
 190     };
 191   };
 192   allocateGrid(m_Output, N_new_cells, N_new_points);
 193   CopyPoints();
 194   CopyCells();
 195   updateCellIndex(m_Output);
 196   N_eliminated = N;
 197 };