src/libengrid/vtkEgEliminateShortEdges.cxx

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