src/libengrid/laplacesmoother.cpp

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