src/libengrid/surfacemeshsmoother.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 "surfacemeshsmoother.h"
  23
  24 SurfaceMeshSmoother::SurfaceMeshSmoother()
  25 {
  26   m_UseEstimatedPlane = false;
  27   m_Cad = NULL;
  28   m_Prec = 1e-3;
  29   m_NumSteps = 5;
  30   m_UseSimpleCentreScheme = false;
  31 }
  32
  33 void SurfaceMeshSmoother::prepareEstimatedPlane()
  34 {
  35   m_N0 = vec3_t(0, 0, 0);
  36   m_X0 = vec3_t(0, 0, 0);
  37   int count = 0;
  38   EG_FORALL_CELLS(id_cell, m_Grid) {
  39     if (isSurface(id_cell, m_Grid)) {
  40       m_N0 += GeometryTools::cellNormal(m_Grid, id_cell);
  41       m_X0 += cellCentre(m_Grid, id_cell);
  42       ++count;
  43     }
  44   }
  45   m_N0.normalise();
  46   m_X0 *= 1.0/count;
  47   m_UseEstimatedPlane = true;
  48 }
  49
  50 void SurfaceMeshSmoother::prepareCadInterface(CadInterface *cad)
  51 {
  52   m_Cad = cad;
  53   m_UseEstimatedPlane = false;
  54 }
  55
  56 vec2_t SurfaceMeshSmoother::transform32(vec3_t x)
  57 {
  58   x -= m_X0;
  59   vec3_t x2 = m_M32*x;
  60   return vec2_t(x2[0], x2[1]);
  61 }
  62
  63 vec3_t SurfaceMeshSmoother::transform23(vec2_t x)
  64 {
  65   vec3_t x3(x[0], x[1], 0);
  66   return m_X0 + m_M23*x3;
  67 }
  68
  69 vec3_t SurfaceMeshSmoother::smoothNode(vtkIdType id_node)
  70 {
  71   // establish 2D coordinate system
  72   vec3_t x_old;
  73   m_Grid->GetPoint(id_node, x_old.data());
  74   if (!m_UseEstimatedPlane) {
  75     m_X0 = m_Cad->snapNode(id_node, x_old);
  76     m_N0 = m_Cad->getLastNormal();
  77   }
  78   m_N0.normalise();
  79   m_M32[0] = GeometryTools::orthogonalVector(m_N0);
  80   m_M32[1] = m_N0.cross(m_M32[0]);
  81   m_M32[2] = m_N0;
  82   m_M23 = m_M32.inverse();
  83
  84   // only smooth simple vertices for now
  85   //EG_STOPDATE("2015-06-01");
  86   EG_VTKDCN(vtkCharArray_t, node_type, m_Grid, "node_type");
  87   if (node_type->GetValue(id_node) != EG_SIMPLE_VERTEX) {
  88     return x_old;
  89   }
  90
  91   // build limiting polygon
  92   m_Limit.clear();
  93   for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
  94     vtkIdType id_cell = m_Part.n2cGG(id_node, i);
  95     EG_GET_CELL(id_cell, m_Grid);
  96     for (int j = 0; j < num_pts; ++j) {
  97       vtkIdType id_node1 = pts[num_pts - 1];
  98       vtkIdType id_node2 = pts[j];
  99       if (j > 0) {
 100         id_node1 = pts[j - 1];
 101       }
 102       if (id_node1 != id_node && id_node2 != id_node) {
 103         vec3_t x1, x2;
 104         m_Grid->GetPoint(id_node1, x1.data());
 105         m_Grid->GetPoint(id_node2, x2.data());
 106         m_Limit << QPair<vec2_t,vec2_t>(transform32(x1), transform32(x2));
 107       }
 108     }
 109   }
 110
 111   vec2_t x(0, 0);
 112   if (m_UseSimpleCentreScheme) {
 113
 114     vec2_t v(0, 0);
 115     for (int i = 0; i < m_Limit.size(); ++i) {
 116       v += 0.5*(m_Limit[i].first + m_Limit[i].second);
 117     }
 118     vec2_t c = v;
 119     c *= 1.0/m_Limit.size();
 120     v.normalise();
 121     if (checkVector(v)) {
 122       double k1, k2;
 123       computeLimits(x, v, k1, k2);
 124       if (k1 > 0) {
 125         double w = 0.9;
 126         x = (w*k2 + (1-w)*k1)*v;
 127       } else {
 128         x = c;
 129       }
 130     }
 131
 132   } else {
 133
 134     // compute smallest segment length
 135     m_MinLength = EG_LARGE_REAL;
 136     for (int i = 0; i < m_Limit.size(); ++i) {
 137       vec2_t a = m_Limit[i].first;
 138       vec2_t b = m_Limit[i].second;
 139       m_MinLength = min(m_MinLength, (a - b).abs());
 140     }
 141
 142     int count = 0;
 143     vec2_t x1(0,0,0);
 144     vec2_t x2 = x1;
 145     do {
 146       x1 = x2;
 147       x2 = iteration(x1);
 148       ++count;
 149     } while ((x2 - x1).abs() > m_Prec*m_MinLength && count < 10);
 150     x = x2;
 151
 152   }
 153
 154   vec3_t x_smooth = transform23(x);
 155   if (!m_UseEstimatedPlane) {
 156     x_smooth = m_Cad->snapNode(id_node, x_smooth);
 157   }
 158   return x_smooth;
 159 }
 160
 161 vec2_t SurfaceMeshSmoother::iteration(vec2_t x)
 162 {
 163   vec2_t v = gradError(x);
 164   v.normalise();
 165   if (!checkVector(v)) {
 166     return x;
 167   }
 168   double k1, k2;
 169   computeLimits(x, v, k1, k2);
 170   double l1 = k1;
 171   double l2 = k2;
 172   while (l2 - l1 > m_Prec) {
 173     int i = 0;
 174     double dk = (l2 - l1)/m_NumSteps;
 175     while (i < m_NumSteps - 2) {
 176       double j1 = l1 + i*dk;
 177       double j2 = j1 + dk;
 178       double err1 = error(x + j1*v);
 179       double err2 = error(x + j2*v);
 180       if (err2 > err1) {
 181         l1 = j1;
 182         l2 = j2;
 183         break;
 184       }
 185       ++i;
 186     }
 187     l1 += i*dk;
 188     l2  = l1 + dk;
 189   }
 190   return x + 0.5*(l1 + l2)*v;
 191 }
 192
 193 double SurfaceMeshSmoother::error(vec2_t x)
 194 {
 195   if (m_Limit.size() == 0) {
 196     EG_BUG;
 197   }
 198
 199   // min and max area
 200   double A_min =  EG_LARGE_REAL;
 201   double A_max = -EG_LARGE_REAL;
 202   for (int i = 0; i < m_Limit.size(); ++i) {
 203     vec2_t a = m_Limit[i].first - x;
 204     vec2_t b = m_Limit[i].second - x;
 205     double A = a[0]*b[1] - a[1]*b[0];
 206     A_min = min(A_min, A);
 207     A_max = max(A_max, A);
 208   }
 209   double err = A_max - A_min;
 210   return err;
 211   /*
 212   //A0 /= m_Limit.count();
 213
 214   // compute error
 215   double err = 0;
 216   for (int i = 0; i < m_Limit.size(); ++i) {
 217     vec2_t a = m_Limit[i].first - x;
 218     vec2_t b = m_Limit[i].second - x;
 219     double A = a[0]*b[1] - a[1]*b[0];
 220     err += sqr(A - A0);
 221   }
 222
 223   return err;
 224   */
 225 }
 226
 227 vec2_t SurfaceMeshSmoother::gradError(vec2_t x)
 228 {
 229   double d = 1e-3*m_MinLength;
 230   vec2_t grad;
 231   grad[0] = (error(x + vec2_t(0.5*d, 0)) - error(x - vec2_t(0.5*d, 0)))/d;
 232   grad[1] = (error(x + vec2_t(0, 0.5*d)) - error(x - vec2_t(0, 0.5*d)))/d;
 233   return grad;
 234 }
 235
 236 void SurfaceMeshSmoother::computeLimits(vec2_t x, vec2_t v, double &k1, double &k2)
 237 {
 238   int count = 0;
 239   k1 = EG_LARGE_REAL;
 240   k2 = -k1;
 241   for (int i = 0; i < m_Limit.size(); ++i) {
 242     vec2_t a = m_Limit[i].first;
 243     vec2_t b = m_Limit[i].second;
 244     double k, l;
 245     GeometryTools::intersection(k, l, x, v, a, (b - a));
 246     if (l >= 0 && l <= 1) {
 247       k1 = min(k1, k);
 248       k2 = max(k2, k);
 249       ++count;
 250     }
 251   }
 252   if (count < 2) {
 253     //EG_BUG;
 254   }
 255 }
 256
 257 void SurfaceMeshSmoother::operate()
 258 {
 259   EG_BUG;
 260 }