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 "geometrytools.h"
  23
  24 vtkStandardNewMacro(vtkEgEliminateShortEdges);
  25
  26 vtkEgEliminateShortEdges::vtkEgEliminateShortEdges()
  27 {
  28   max_ratio = 1000.0;
  29   max_length = 1e99;
  30 };
  31
  32 void vtkEgEliminateShortEdges::CheckEdges()
  33 {
  34   new_node_idx.fill(-1, m_Input->GetNumberOfPoints());
  35   delete_cell.fill(false, m_Input->GetNumberOfCells());
  36   QVector<bool> marked(m_Input->GetNumberOfPoints(), false);
  37   for (vtkIdType cellId = 0; cellId < m_Input->GetNumberOfCells(); ++cellId) {
  38     int cellType = m_Input->GetCellType(cellId);
  39     if ((cellType == VTK_TRIANGLE) || (cellType == VTK_QUAD)) {
  40       vtkIdType *pts, Npts;
  41       m_Input->GetCellPoints(cellId, Npts, pts);
  42       double L_av = 0;
  43       vector<vec3_t> x(Npts);
  44       int N_del = 0;
  45       double Lmin = 1e99;
  46       int imin = 0;
  47       for (int i = 0; i < Npts; ++i) {
  48         vtkIdType p1 = pts[i];
  49         vtkIdType p2;
  50         if (i < Npts-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 (Npts == 4) {
  72           A = GeometryTools::quadArea(x[0],x[1],x[2],x[3]);
  73         };
  74         L_av /= Npts;
  75         for (int i = 0; i < Npts; ++i) {
  76           vtkIdType p1 = pts[i];
  77           vtkIdType p2;
  78           if (i < Npts-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       vtkIdType *pts, Npts;
 114       m_Input->GetCellPoints(cellId, Npts, pts);
 115       for (int i = 0; i < Npts; ++i) {
 116         vtkIdType p1 = pts[i];
 117         vtkIdType p2;
 118         if (i < Npts-1) p2 = pts[i+1];
 119         else            p2 = pts[0];
 120         if (p1 == p2) {
 121           EG_ERR_RETURN("bug encountered");
 122         };
 123         if (new_node_idx[p1] == p2) {
 124           delete_cell[cellId] = true;
 125         };
 126         if (new_node_idx[p2] == p1) {
 127           delete_cell[cellId] = true;
 128         };
 129       };
 130     };
 131   };
 132 };
 133
 134 void vtkEgEliminateShortEdges::CopyPoints()
 135 {
 136   node_mapping.fill(-1, m_Input->GetNumberOfPoints());
 137   vtkIdType newPointId = 0;
 138   for (vtkIdType pointId = 0; pointId < m_Input->GetNumberOfPoints(); ++pointId) {
 139     if(new_node_idx[pointId] < 0) {
 140       node_mapping[pointId] = newPointId;
 141       vec3_t x;
 142       m_Input->GetPoints()->GetPoint(pointId,x.data());
 143       m_Output->GetPoints()->SetPoint(newPointId,x.data());
 144       ++newPointId;
 145     };
 146   };
 147   for (vtkIdType pointId = 0; pointId < m_Input->GetNumberOfPoints(); ++pointId) {
 148     if(new_node_idx[pointId] >= 0) {
 149       if (node_mapping[new_node_idx[pointId]] < 0) {
 150         EG_ERR_RETURN("bug encountered");
 151       };
 152       node_mapping[pointId] = node_mapping[new_node_idx[pointId]];
 153     };
 154   };
 155 };
 156
 157 void vtkEgEliminateShortEdges::CopyCells()
 158 {
 159   for (vtkIdType cellId = 0; cellId < m_Input->GetNumberOfCells(); ++cellId) {
 160     if(!delete_cell[cellId]) {
 161       vtkIdType *old_pts;
 162       vtkIdType  Npts;
 163       m_Input->GetCellPoints(cellId, Npts, old_pts);
 164       vtkIdType *new_pts = new vtkIdType[Npts];
 165       for (int i = 0; i < Npts; ++i) {
 166         new_pts[i] = node_mapping[old_pts[i]];
 167       };
 168       vtkIdType newCellId = m_Output->InsertNextCell(m_Input->GetCellType(cellId), Npts, new_pts);
 169       copyCellData(m_Input, cellId, m_Output, newCellId);
 170       delete [] new_pts;
 171     };
 172   };
 173 };
 174
 175 void vtkEgEliminateShortEdges::ExecuteEg()
 176 {
 177   N_eliminated = 0;
 178   N_new_points = 0;
 179   N_new_cells  = 0;
 180   CheckEdges();
 181   CheckCells();
 182   int N = 0;
 183   for (vtkIdType pointId = 0; pointId < m_Input->GetNumberOfPoints(); ++pointId) {
 184     if(new_node_idx[pointId] < 0) {
 185       ++N_new_points;
 186     };
 187   };
 188   for (vtkIdType cellId = 0; cellId < m_Input->GetNumberOfCells(); ++cellId) {
 189     if(!delete_cell[cellId]) {
 190       ++N_new_cells;
 191     } else {
 192       ++N;
 193     };
 194   };
 195   allocateGrid(m_Output, N_new_cells, N_new_points);
 196   CopyPoints();
 197   CopyCells();
 198   UpdateCellIndex(m_Output);
 199   N_eliminated = N;
 200 };