src/libengrid/laplacesmoother.cpp

   1 //
   2 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   3 // +                                                                      +
   4 // + This file is part of enGrid.                                         +
   5 // +                                                                      +
   6 // + Copyright 2008-2012 enGits GmbH                                     +
   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 "laplacesmoother.h"
  24 #include <vtkCellLocator.h>
  25 #include <vtkCharArray.h>
  26 #include <vtkGenericCell.h>
  27
  28 #include "guimainwindow.h"
  29 #include "localnodegraphinterface.h"
  30 #include "checkerboardgraphiterator.h"
  31
  32 using namespace GeometryTools;
  33
  34 LaplaceSmoother::LaplaceSmoother() : SurfaceOperation()
  35 {
  36   DebugLevel = 0;
  37   setQuickSave(true);
  38   m_UseProjection = true;
  39   //   m_UseNormalCorrection = false;
  40   getSet("surface meshing", "under relaxation for smoothing", 0.5, m_UnderRelaxation);
  41   getSet("surface meshing", "feature magic", 0.0, m_FeatureMagic);
  42   getSet("surface meshing", "smoothing limiter", 1.0, m_Limit);
  43   getSet("surface meshing", "use uniform smoothing", false, m_UniformSnapPoints);
  44   m_Limit = min(1.0, max(0.0, m_Limit));
  45   m_NoCheck = false;
  46   m_ProjectionIterations = 50;
  47   m_AllowedCellTypes.clear();
  48   m_AllowedCellTypes.insert(VTK_TRIANGLE);
  49   //m_StrictFeatureSnap = false;
  50 }
  51
  52 bool LaplaceSmoother::setNewPosition(vtkIdType id_node, vec3_t x_new)
  53 {
  54   using namespace GeometryTools;
  55
  56   vec3_t x_old;
  57   m_Grid->GetPoint(id_node, x_old.data());
  58   EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type");
  59   bool move = true;
  60   if(m_NoCheck) {
  61     return move;
  62   }
  63   QVector<vec3_t> old_cell_normals(m_Part.n2cGSize(id_node));
  64   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
  65   for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
  66     old_cell_normals[i] = GeometryTools::cellNormal(m_Grid, m_Part.n2cGG(id_node, i));
  67   }
  68   m_Grid->GetPoints()->SetPoint(id_node, x_new.data());
  69
  70   for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
  71     vec3_t n = GeometryTools::cellNormal(m_Grid, m_Part.n2cGG(id_node, i));
  72     if (n*old_cell_normals[i] < 0.2*old_cell_normals[i].abs2()) {
  73       move = false;
  74       break;
  75     }
  76   }
  77   if (!move) {
  78     m_Grid->GetPoints()->SetPoint(id_node, x_old.data());
  79   }
  80   return move;
  81 }
  82
  83 void LaplaceSmoother::featureCorrection(vtkIdType id_node, SurfaceProjection* proj, vec3_t &x_new)
  84 {
  85   if (m_FeatureMagic > 0) {
  86
  87     EG_VTKDCN(vtkDoubleArray, cl, m_Grid, "node_meshdensity_desired");
  88     EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type");
  89     bool convex = isConvexNode(id_node);
  90     if (node_type->GetValue(id_node) == EG_FEATURE_CORNER_VERTEX) {
  91
  92       vec3_t x;
  93       double L = 0.1*cl->GetValue(id_node);
  94
  95       // do not use curvature correction here
  96       // .. a proper CAD model does not need it
  97       // .. a discrete model (e.g. triangulation) will create bulges on features
  98       vec3_t x0 = proj->project(x_new, id_node, false, m_NodeNormal[id_node]);
  99
 100       if (convex) {
 101         x = x0 - L*m_NodeNormal[id_node];
 102       } else {
 103         x = x0 + L*m_NodeNormal[id_node];
 104       }
 105       vec3_t n = proj->lastProjNormal();
 106       if (!proj->lastProjFailed()) {
 107         double d = 2*L/tan(0.5*m_FeatureAngle);
 108         static const int num_steps = 36;
 109         double D_alpha = 2*M_PI/num_steps;
 110         vec3_t v;
 111
 112         v = GeometryTools::orthogonalVector(m_NodeNormal[id_node]);
 113         int num_miss = 0;
 114         int num_hit = 0;
 115         vec3_t x_corner(0,0,0);
 116         for (int i = 0; i < num_steps; ++i) {
 117           v = GeometryTools::rotate(v, m_NodeNormal[id_node], D_alpha);
 118           vec3_t xp = proj->project(x, id_node, true, v, true);
 119           if (proj->lastProjFailed()) {
 120             ++num_miss;
 121           } else {
 122             double l = (x - xp).abs();
 123             if (l < d) {
 124               ++num_hit;
 125               x_corner += xp;
 126             } else {
 127               ++num_miss;
 128             }
 129           }
 130         }
 131         if (num_miss == 0 && num_hit > 0) {
 132           x_corner *= 1.0/num_hit;
 133           x_new = proj->project(x_corner, id_node, true, m_NodeNormal[id_node]);
 134         }
 135       }
 136
 137     } else {
 138
 139       // "magic" vector to displace node for re-projection
 140       vec3_t x0  = x_new;
 141       vec3_t x1  = proj->project(x_new, id_node, m_CorrectCurvature);
 142       vec3_t n   = proj->lastProjNormal();
 143       double l   = cl->GetValue(id_node);
 144       double eps = 0.01*l;
 145       vec3_t mv  = l*m_NodeNormal[id_node];
 146       mv -= (n*mv)*n;
 147
 148       if (checkVector(mv)) {
 149         double L1 = 0;
 150         double L2 = 1;//m_FeatureMagic*cl->GetValue(id_node);
 151         bool flipped = false;
 152         double amp = 0.1;
 153         int i = 0;
 154         int hits = 0;
 155         while (i < 30 && amp < 10) {
 156           x_new = x1 + 0.5*amp*(L1 + L2)*mv;// + 2*eps*n;
 157           //x_new = proj->findClosest(x_new, id_node, n);
 158           x_new = proj->project(x_new, id_node, m_CorrectCurvature, n);
 159           double displacement = fabs((x_new - x1)*n);
 160           if (displacement > eps || proj->lastProjFailed()) {
 161             L2 = 0.5*(L1 + L2);
 162             ++i;
 163           } else {
 164
 165             // if there is no significant displacement after the first iteration
 166             // the node is probably in a smooth region of the surface
 167             // ==> stop here
 168             //
 169             if (i == 0) {
 170               if (!flipped) {
 171                 mv *= -1;
 172                 flipped = true;
 173               } else {
 174                 amp *= 1.5;
 175                 mv *= -1;
 176                 flipped = false;
 177               }
 178             } else {
 179               L1 = 0.5*(L1 + L2);
 180               ++hits;
 181               ++i;
 182             }
 183           }
 184         }
 185         if (hits > 0) {
 186           x_new = x1 + L1*m_FeatureMagic*amp*mv;
 187           x_new = proj->project(x_new, id_node, m_CorrectCurvature, n);
 188           if (proj->lastProjFailed()) {
 189             cout << "bad!" << endl;
 190           }
 191         } else {
 192           x_new = x0;
 193         }
 194       }
 195     }
 196   }
 197 }
 198
 199 bool LaplaceSmoother::moveNode(vtkIdType id_node, vec3_t &Dx)
 200 {
 201   if (!checkVector(Dx)) {
 202     return false;
 203   }
 204   EG_VTKDCN(vtkDoubleArray, cl, m_Grid, "node_meshdensity_desired");
 205   vec3_t x_old;
 206   m_Grid->GetPoint(id_node, x_old.data());
 207   bool moved = false;
 208
 209   for (int i_relaxation = 0; i_relaxation < 10; ++i_relaxation) {
 210     vec3_t x_new = x_old + Dx;
 211     if (m_UseProjection) {
 212       int i_nodes = m_Part.localNode(id_node);
 213       if (m_NodeToBc[i_nodes].size() == 1 || m_UniformSnapPoints) {
 214         int bc = m_NodeToBc[i_nodes][0];
 215         x_new = GuiMainWindow::pointer()->getSurfProj(bc)->project(x_new, id_node, m_CorrectCurvature);
 216         featureCorrection(id_node, GuiMainWindow::pointer()->getSurfProj(bc), x_new);
 217       } else {
 218         for (int i_proj_iter = 0; i_proj_iter < m_ProjectionIterations; ++i_proj_iter) {
 219           foreach (int bc, m_NodeToBc[i_nodes]) {
 220             x_new = GuiMainWindow::pointer()->getSurfProj(bc)->project(x_new, id_node, m_CorrectCurvature);
 221           }
 222         }
 223
 224         for (int i_proj_iter = 0; i_proj_iter < m_ProjectionIterations; ++i_proj_iter) {
 225           if (m_CorrectCurvature) {
 226             foreach (int bc, m_NodeToBc[i_nodes]) {
 227               x_new = GuiMainWindow::pointer()->getSurfProj(bc)->correctCurvature(GuiMainWindow::pointer()->getSurfProj(bc)->lastProjTriangle(), x_new);
 228             }
 229           }
 230         }
 231
 232       }
 233     }
 234
 235     // compute the minimal length of any edge adjacent to this node
 236     // .. This will be used to limit the node movement.
 237     // .. Hopefully jammed topologies can be avoided this way.
 238     //
 239     EG_VTKDCN(vtkDoubleArray, cl, m_Grid, "node_meshdensity_desired");
 240     vec3_t x_old;
 241     m_Grid->GetPoint(id_node, x_old.data());
 242     double L_min = cl->GetValue(id_node);
 243     for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
 244       vtkIdType id_neigh = m_Part.n2nGG(id_node, i);
 245       vec3_t x_neigh;
 246       m_Grid->GetPoint(id_neigh, x_neigh.data());
 247       L_min = min(L_min, (x_old - x_neigh).abs());
 248     }
 249
 250     // limit node displacement
 251     vec3_t dx = x_new - x_old;
 252     if (dx.abs() > m_Limit*L_min) {
 253       x_new -= dx;
 254       dx.normalise();
 255       x_new += m_Limit*L_min*dx;
 256     }
 257
 258     if (setNewPosition(id_node, x_new)) {
 259       m_Grid->GetPoints()->SetPoint(id_node, x_new.data());
 260       moved = true;
 261       Dx = x_new - x_old;
 262       break;
 263     }
 264     Dx *= 0.5;
 265   }
 266   return moved;
 267 }
 268
 269 void LaplaceSmoother::fixNodes(const QVector<bool> &fixnodes)
 270 {
 271   if (fixnodes.size() != m_Grid->GetNumberOfPoints()) {
 272     EG_BUG;
 273   }
 274   m_Fixed = fixnodes;
 275 }
 276
 277 void LaplaceSmoother::operate()
 278 {
 279   if (m_BCodeFeatureDefinition) {
 280     m_FeatureMagic = 0.0;
 281     m_NoCheck = false;
 282   } else {
 283     m_NoCheck = true;
 284   }
 285   if (m_Fixed.size() != m_Grid->GetNumberOfPoints()) {
 286     m_Fixed.fill(false, m_Grid->GetNumberOfPoints());
 287   }
 288   QVector<int> bcs;
 289   GuiMainWindow::pointer()->getAllBoundaryCodes(bcs);
 290   if (m_UseProjection) {
 291     foreach (int bc, bcs) {
 292       GuiMainWindow::pointer()->getSurfProj(bc)->setForegroundGrid(m_Grid);
 293     }
 294   }
 295   updateNodeInfo();
 296   EG_VTKDCC(vtkIntArray,    cell_code, m_Grid, "cell_code");
 297   EG_VTKDCN(vtkCharArray,   node_type, m_Grid, "node_type" );
 298   EG_VTKDCN(vtkDoubleArray, cl,        m_Grid, "node_meshdensity_desired");
 299   QVector<vtkIdType> smooth_node(m_Grid->GetNumberOfPoints(), false);
 300   {
 301     l2g_t nodes = m_Part.getNodes();
 302     foreach (vtkIdType id_node, nodes) {
 303       smooth_node[id_node] = true;
 304     }
 305   }
 306   setAllSurfaceCells();
 307   l2g_t  nodes = m_Part.getNodes();
 308   m_NodeToBc.resize(nodes.size());
 309   for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
 310     QSet<int> bcs;
 311     for (int j = 0; j < m_Part.n2cLSize(i_nodes); ++j) {
 312       bcs.insert(cell_code->GetValue(m_Part.n2cLG(i_nodes, j)));
 313     }
 314     m_NodeToBc[i_nodes].resize(bcs.size());
 315     qCopy(bcs.begin(), bcs.end(), m_NodeToBc[i_nodes].begin());
 316   }
 317
 318   QVector<vec3_t> x_new(nodes.size());
 319
 320   QVector<bool> blocked(nodes.size(), false);
 321   for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
 322     for (int i = 0; i < m_Part.n2cLSize(i_nodes); ++i) {
 323       vtkIdType type = m_Grid->GetCellType(m_Part.n2cLG(i_nodes, i));
 324       if (!m_AllowedCellTypes.contains(type)) {
 325         blocked[i_nodes] = true;
 326         break;
 327       }
 328     }
 329   }
 330
 331   for (int i_iter = 0; i_iter < m_NumberOfIterations; ++i_iter) {
 332     m_Success = true;
 333     computeNormals();
 334
 335     for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
 336       vtkIdType id_node = nodes[i_nodes];
 337       if (!m_Fixed[id_node] && !blocked[i_nodes]) {
 338         if (smooth_node[id_node] && node_type->GetValue(id_node) != EG_FIXED_VERTEX) {
 339           if (node_type->GetValue(id_node) != EG_FIXED_VERTEX) {
 340             QVector<vtkIdType> snap_points = getPotentialSnapPoints(id_node);
 341             vec3_t n(0,0,0);
 342             vec3_t x_old;
 343             m_Grid->GetPoint(id_node, x_old.data());
 344
 345             if (snap_points.size() > 0) {
 346               vec3_t x;
 347               x_new[i_nodes] = vec3_t(0,0,0);
 348               double w_tot = 0;
 349               double L_min = 1e99;
 350               foreach (vtkIdType id_snap_node, snap_points) {
 351                 m_Grid->GetPoint(id_snap_node, x.data());
 352                 double w = 1.0;
 353                 w_tot += w;
 354                 x_new[i_nodes] += w*x;
 355                 n += m_NodeNormal[id_snap_node];
 356                 double L = (x - x_old).abs();
 357                 L_min = min(L, L_min);
 358               }
 359               n.normalise();
 360               x_new[i_nodes] *= 1.0/w_tot;
 361             } else {
 362               x_new[i_nodes] = x_old;
 363             }
 364
 365             if (m_UseNormalCorrection) {
 366               vec3_t dx = x_new[i_nodes] - x_old;
 367               double scal = dx*m_NodeNormal[id_node];
 368               x_new[i_nodes] += scal*m_NodeNormal[id_node];
 369             }
 370             vec3_t Dx = x_new[i_nodes] - x_old;
 371             //Dx *= m_UnderRelaxation;
 372             if (moveNode(id_node, Dx)) {
 373               x_new[i_nodes] = x_old + m_UnderRelaxation*Dx;
 374             } else {
 375               x_new[i_nodes] = x_old;
 376               //m_Success = false;
 377             }
 378             m_Grid->GetPoints()->SetPoint(id_node, x_old.data());
 379           }
 380         }
 381       }
 382     }
 383     for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
 384       vtkIdType id_node = nodes[i_nodes];
 385       m_Grid->GetPoints()->SetPoint(id_node, x_new[id_node].data());
 386     }
 387
 388     if (m_Success) {
 389       break;
 390     }
 391   }
 392 }