src/libengrid/geometrysmoother.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 //
  22 #include "geometrysmoother.h"
  23 #include "geometrytools.h"
  24
  25 GeometrySmoother::GeometrySmoother()
  26 {
  27   m_CornerAngle      = GeometryTools::deg2rad(20.0);
  28   m_RelaxationFactor = 0.1;
  29   m_NumIterations    = 10;
  30   m_ConvexRelaxation = 0.0;
  31 }
  32
  33 void GeometrySmoother::operate()
  34 {
  35   // copy all node coordinates into a QVector
  36   QVector<vec3_t> x(m_Grid->GetNumberOfPoints());
  37   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
  38     m_Grid->GetPoint(id_node, x[id_node].data());
  39   }
  40
  41   // compute "snap" points
  42   int N = 0;
  43   QVector<QList<vtkIdType> > snap_points(m_Grid->GetNumberOfPoints());
  44   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
  45     QList<vtkIdType> boundary_neighbours;
  46     for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
  47       vtkIdType id_neigh = m_Part.n2nGG(id_node, i);
  48       QList<vtkIdType> edge_faces;
  49       m_Part.getEdgeFaces(id_node, id_neigh, edge_faces);
  50       if (edge_faces.size() > 2) {
  51         EG_BUG;
  52       }
  53       else if (edge_faces.size() < 1) {
  54         EG_BUG;
  55       }
  56       else if (edge_faces.size() == 1) {
  57         boundary_neighbours << id_neigh;
  58       }
  59       else {
  60         snap_points[id_node] << id_neigh;
  61       }
  62     }
  63     if (boundary_neighbours.size() > 2) {
  64       EG_BUG;
  65     }
  66     if (boundary_neighbours.size() == 2) {
  67       vec3_t u = x[id_node] - x[boundary_neighbours[0]];
  68       vec3_t v = x[boundary_neighbours[1]] - x[id_node];
  69       vec3_t n = m_Part.globalNormal(id_node);
  70       u -= (n*u)*n;
  71       v -= (n*v)*n;
  72       if (GeometryTools::angle(u, v) > m_CornerAngle) {
  73         snap_points[id_node].clear();
  74         ++N;
  75       }
  76       else {
  77         snap_points[id_node] = boundary_neighbours;
  78       }
  79     }
  80   }
  81
  82   // the smoothing process starts here
  83   for (int iter = 0; iter < m_NumIterations; ++iter) {
  84     QVector<vec3_t> x_new = x;
  85     for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
  86       if (snap_points[id_node].size()) {
  87         double cf_corr = 0;
  88         if (m_ConvexRelaxation > 1e-6) {
  89           vec3_t n = m_Part.globalNormal(id_node);
  90           double cf = 0; // convexity factor ...
  91           foreach (vtkIdType id_snap, snap_points[id_node]) {
  92             vec3_t v = x[id_node] - x[id_snap];
  93             v.normalise();
  94             cf += v*n;
  95           }
  96           cf /= snap_points[id_node].size();
  97           cf_corr = min(1.0, max(-1.0, m_ConvexRelaxation*cf));
  98         }
  99         double w = m_RelaxationFactor*(1.0 + cf_corr)/snap_points[id_node].size();
 100         if (snap_points[id_node].size() == 2) {
 101           w *= 0.3333;
 102         }
 103         foreach (vtkIdType id_snap, snap_points[id_node]) {
 104           x_new[id_node] += w*(x[id_snap] - x[id_node]);
 105         }
 106       }
 107     }
 108     x = x_new;
 109   }
 110
 111   // copy node coordinates back to the VTK grid
 112   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 113     m_Grid->GetPoints()->SetPoint(id_node, x[id_node].data());
 114   }
 115 }