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       vec3_t x0 = proj->project(x_new, id_node, true, m_NodeNormal[id_node]);
  96       if (convex) {
  97         x = x0 - L*m_NodeNormal[id_node];
  98       } else {
  99         x = x0 + L*m_NodeNormal[id_node];
 100       }
 101       vec3_t n = proj->lastProjNormal();
 102       if (!proj->lastProjFailed()) {
 103         double d = 2*L/tan(0.5*m_FeatureAngle);
 104         static const int num_steps = 36;
 105         double D_alpha = 2*M_PI/num_steps;
 106         vec3_t v;
 107
 108         v = GeometryTools::orthogonalVector(m_NodeNormal[id_node]);
 109         int num_miss = 0;
 110         int num_hit = 0;
 111         vec3_t x_corner(0,0,0);
 112         for (int i = 0; i < num_steps; ++i) {
 113           v = GeometryTools::rotate(v, m_NodeNormal[id_node], D_alpha);
 114           vec3_t xp = proj->project(x, id_node, true, v, true);
 115           if (proj->lastProjFailed()) {
 116             ++num_miss;
 117           } else {
 118             double l = (x - xp).abs();
 119             if (l < d) {
 120               ++num_hit;
 121               x_corner += xp;
 122             } else {
 123               ++num_miss;
 124             }
 125           }
 126         }
 127         if (num_miss == 0 && num_hit > 0) {
 128           x_corner *= 1.0/num_hit;
 129           x_new = proj->project(x_corner, id_node, true, m_NodeNormal[id_node]);
 130         }
 131       }
 132
 133     } else {
 134
 135       // "magic" vector to displace node for re-projection
 136       vec3_t x0  = x_new;
 137       vec3_t x1  = proj->project(x_new, id_node, m_CorrectCurvature);
 138       vec3_t n   = proj->lastProjNormal();
 139       double l   = cl->GetValue(id_node);
 140       double eps = 0.01*l;
 141       vec3_t mv  = l*m_NodeNormal[id_node];
 142       mv -= (n*mv)*n;
 143
 144       if (checkVector(mv)) {
 145         double L1 = 0;
 146         double L2 = 1;//m_FeatureMagic*cl->GetValue(id_node);
 147         bool flipped = false;
 148         double amp = 0.1;
 149         int i = 0;
 150         int hits = 0;
 151         while (i < 30 && amp < 10) {
 152           x_new = x1 + 0.5*amp*(L1 + L2)*mv;// + 2*eps*n;
 153           //x_new = proj->findClosest(x_new, id_node, n);
 154           x_new = proj->project(x_new, id_node, m_CorrectCurvature, n);
 155           double displacement = fabs((x_new - x1)*n);
 156           if (displacement > eps || proj->lastProjFailed()) {
 157             L2 = 0.5*(L1 + L2);
 158             ++i;
 159           } else {
 160
 161             // if there is no significant displacement after the first iteration
 162             // the node is probably in a smooth region of the surface
 163             // ==> stop here
 164             //
 165             if (i == 0) {
 166               if (!flipped) {
 167                 mv *= -1;
 168                 flipped = true;
 169               } else {
 170                 amp *= 1.5;
 171                 mv *= -1;
 172                 flipped = false;
 173               }
 174             } else {
 175               L1 = 0.5*(L1 + L2);
 176               ++hits;
 177               ++i;
 178             }
 179           }
 180         }
 181         if (hits > 0) {
 182           x_new = x1 + L1*m_FeatureMagic*amp*mv;
 183           x_new = proj->project(x_new, id_node, m_CorrectCurvature, n);
 184           if (proj->lastProjFailed()) {
 185             cout << "bad!" << endl;
 186           }
 187         } else {
 188           x_new = x0;
 189         }
 190       }
 191     }
 192   }
 193 }
 194
 195 bool LaplaceSmoother::moveNode(vtkIdType id_node, vec3_t &Dx)
 196 {
 197   if (!checkVector(Dx)) {
 198     return false;
 199   }
 200   EG_VTKDCN(vtkDoubleArray, cl, m_Grid, "node_meshdensity_desired");
 201   vec3_t x_old;
 202   m_Grid->GetPoint(id_node, x_old.data());
 203   bool moved = false;
 204
 205   for (int i_relaxation = 0; i_relaxation < 10; ++i_relaxation) {
 206     vec3_t x_new = x_old + Dx;
 207     if (m_UseProjection) {
 208       int i_nodes = m_Part.localNode(id_node);
 209       if (m_NodeToBc[i_nodes].size() == 1 || m_UniformSnapPoints) {
 210         int bc = m_NodeToBc[i_nodes][0];
 211         x_new = GuiMainWindow::pointer()->getSurfProj(bc)->project(x_new, id_node, m_CorrectCurvature);
 212         featureCorrection(id_node, GuiMainWindow::pointer()->getSurfProj(bc), x_new);
 213       } else {
 214         for (int i_proj_iter = 0; i_proj_iter < m_ProjectionIterations; ++i_proj_iter) {
 215           foreach (int bc, m_NodeToBc[i_nodes]) {
 216             x_new = GuiMainWindow::pointer()->getSurfProj(bc)->project(x_new, id_node, m_CorrectCurvature);
 217           }
 218         }
 219
 220         for (int i_proj_iter = 0; i_proj_iter < m_ProjectionIterations; ++i_proj_iter) {
 221           if (m_CorrectCurvature) {
 222             foreach (int bc, m_NodeToBc[i_nodes]) {
 223               x_new = GuiMainWindow::pointer()->getSurfProj(bc)->correctCurvature(GuiMainWindow::pointer()->getSurfProj(bc)->lastProjTriangle(), x_new);
 224             }
 225           }
 226         }
 227
 228       }
 229     }
 230
 231     // compute the minimal length of any edge adjacent to this node
 232     // .. This will be used to limit the node movement.
 233     // .. Hopefully jammed topologies can be avoided this way.
 234     //
 235     EG_VTKDCN(vtkDoubleArray, cl, m_Grid, "node_meshdensity_desired");
 236     vec3_t x_old;
 237     m_Grid->GetPoint(id_node, x_old.data());
 238     double L_min = cl->GetValue(id_node);
 239     for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
 240       vtkIdType id_neigh = m_Part.n2nGG(id_node, i);
 241       vec3_t x_neigh;
 242       m_Grid->GetPoint(id_neigh, x_neigh.data());
 243       L_min = min(L_min, (x_old - x_neigh).abs());
 244     }
 245
 246     // limit node displacement
 247     vec3_t dx = x_new - x_old;
 248     if (dx.abs() > m_Limit*L_min) {
 249       x_new -= dx;
 250       dx.normalise();
 251       x_new += m_Limit*L_min*dx;
 252     }
 253
 254     if (setNewPosition(id_node, x_new)) {
 255       m_Grid->GetPoints()->SetPoint(id_node, x_new.data());
 256       moved = true;
 257       Dx = x_new - x_old;
 258       break;
 259     }
 260     Dx *= 0.5;
 261   }
 262   return moved;
 263 }
 264
 265 void LaplaceSmoother::fixNodes(const QVector<bool> &fixnodes)
 266 {
 267   if (fixnodes.size() != m_Grid->GetNumberOfPoints()) {
 268     EG_BUG;
 269   }
 270   m_Fixed = fixnodes;
 271 }
 272
 273 void LaplaceSmoother::operate()
 274 {
 275   if (m_BCodeFeatureDefinition) {
 276     m_FeatureMagic = 0.0;
 277     m_NoCheck = false;
 278   } else {
 279     m_NoCheck = true;
 280   }
 281   if (m_Fixed.size() != m_Grid->GetNumberOfPoints()) {
 282     m_Fixed.fill(false, m_Grid->GetNumberOfPoints());
 283   }
 284   QVector<int> bcs;
 285   GuiMainWindow::pointer()->getAllBoundaryCodes(bcs);
 286   if (m_UseProjection) {
 287     foreach (int bc, bcs) {
 288       GuiMainWindow::pointer()->getSurfProj(bc)->setForegroundGrid(m_Grid);
 289     }
 290   }
 291   updateNodeInfo();
 292   EG_VTKDCC(vtkIntArray,    cell_code, m_Grid, "cell_code");
 293   EG_VTKDCN(vtkCharArray,   node_type, m_Grid, "node_type" );
 294   EG_VTKDCN(vtkDoubleArray, cl,        m_Grid, "node_meshdensity_desired");
 295   QVector<vtkIdType> smooth_node(m_Grid->GetNumberOfPoints(), false);
 296   {
 297     l2g_t nodes = m_Part.getNodes();
 298     foreach (vtkIdType id_node, nodes) {
 299       smooth_node[id_node] = true;
 300     }
 301   }
 302   setAllSurfaceCells();
 303   l2g_t  nodes = m_Part.getNodes();
 304   m_NodeToBc.resize(nodes.size());
 305   for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
 306     QSet<int> bcs;
 307     for (int j = 0; j < m_Part.n2cLSize(i_nodes); ++j) {
 308       bcs.insert(cell_code->GetValue(m_Part.n2cLG(i_nodes, j)));
 309     }
 310     m_NodeToBc[i_nodes].resize(bcs.size());
 311     qCopy(bcs.begin(), bcs.end(), m_NodeToBc[i_nodes].begin());
 312   }
 313
 314   QVector<vec3_t> x_new(nodes.size());
 315
 316   QVector<bool> blocked(nodes.size(), false);
 317   for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
 318     for (int i = 0; i < m_Part.n2cLSize(i_nodes); ++i) {
 319       vtkIdType type = m_Grid->GetCellType(m_Part.n2cLG(i_nodes, i));
 320       if (!m_AllowedCellTypes.contains(type)) {
 321         blocked[i_nodes] = true;
 322         break;
 323       }
 324     }
 325   }
 326
 327   for (int i_iter = 0; i_iter < m_NumberOfIterations; ++i_iter) {
 328     m_Success = true;
 329     computeNormals();
 330
 331     for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
 332       vtkIdType id_node = nodes[i_nodes];
 333       if (!m_Fixed[id_node] && !blocked[i_nodes]) {
 334         if (smooth_node[id_node] && node_type->GetValue(id_node) != EG_FIXED_VERTEX) {
 335           if (node_type->GetValue(id_node) != EG_FIXED_VERTEX) {
 336             QVector<vtkIdType> snap_points = getPotentialSnapPoints(id_node);
 337             vec3_t n(0,0,0);
 338             vec3_t x_old;
 339             m_Grid->GetPoint(id_node, x_old.data());
 340
 341             if (snap_points.size() > 0) {
 342               vec3_t x;
 343               x_new[i_nodes] = vec3_t(0,0,0);
 344               double w_tot = 0;
 345               double L_min = 1e99;
 346               foreach (vtkIdType id_snap_node, snap_points) {
 347                 m_Grid->GetPoint(id_snap_node, x.data());
 348                 double w = 1.0;
 349                 w_tot += w;
 350                 x_new[i_nodes] += w*x;
 351                 n += m_NodeNormal[id_snap_node];
 352                 double L = (x - x_old).abs();
 353                 L_min = min(L, L_min);
 354               }
 355               n.normalise();
 356               x_new[i_nodes] *= 1.0/w_tot;
 357             } else {
 358               x_new[i_nodes] = x_old;
 359             }
 360
 361             if (m_UseNormalCorrection) {
 362               vec3_t dx = x_new[i_nodes] - x_old;
 363               double scal = dx*m_NodeNormal[id_node];
 364               x_new[i_nodes] += scal*m_NodeNormal[id_node];
 365             }
 366             vec3_t Dx = x_new[i_nodes] - x_old;
 367             //Dx *= m_UnderRelaxation;
 368             if (moveNode(id_node, Dx)) {
 369               x_new[i_nodes] = x_old + m_UnderRelaxation*Dx;
 370             } else {
 371               x_new[i_nodes] = x_old;
 372               //m_Success = false;
 373             }
 374             m_Grid->GetPoints()->SetPoint(id_node, x_old.data());
 375           }
 376         }
 377       }
 378     }
 379     for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
 380       vtkIdType id_node = nodes[i_nodes];
 381       m_Grid->GetPoints()->SetPoint(id_node, x_new[id_node].data());
 382     }
 383
 384     if (m_Success) {
 385       break;
 386     }
 387   }
 388 }