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 "globalnodegraphinterface.h"
  30
  31 using namespace GeometryTools;
  32
  33 LaplaceSmoother::LaplaceSmoother() : SurfaceOperation()
  34 {
  35   DebugLevel = 0;
  36   setQuickSave(true);
  37   m_UseProjection = true;
  38 //   m_UseNormalCorrection = false;
  39   getSet("surface meshing", "under relaxation for smoothing", 0.5, m_UnderRelaxation);
  40   getSet("surface meshing", "feature magic", 0.0, m_FeatureMagic);
  41   getSet("surface meshing", "smoothing limiter", 1.0, m_Limit);
  42   getSet("surface meshing", "use uniform smoothing", false, m_UniformSnapPoints);
  43   m_Limit = min(1.0, max(0.0, m_Limit));
  44   m_NoCheck = false;
  45   m_ProjectionIterations = 50;
  46   m_AllowedCellTypes.clear();
  47   m_AllowedCellTypes.insert(VTK_TRIANGLE);
  48 }
  49
  50 bool LaplaceSmoother::setNewPosition(vtkIdType id_node, vec3_t x_new)
  51 {
  52   using namespace GeometryTools;
  53
  54   vec3_t x_old;
  55   m_Grid->GetPoint(id_node, x_old.data());
  56   m_Grid->GetPoints()->SetPoint(id_node, x_new.data());
  57
  58   bool move = true;
  59   if(m_NoCheck) return move;
  60
  61   // compute the extrusion vector to compute the tetrahedrons for volume checking
  62   // start with an average of all adjacent cell normals and count the number of
  63   // adjacent boundary codes (for one boundary code an alternative vector can be
  64   // computed with the help of a SurfaceProjection
  65   //
  66   vec3_t n(0,0,0);
  67   QVector<vec3_t> cell_normals(m_Part.n2cGSize(id_node));
  68   QSet<int> bcs;
  69   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
  70   double A_max = 0; //area of the biggest neighbour cell of id_node
  71   for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
  72     double A = fabs(GeometryTools::cellVA(m_Grid, m_Part.n2cGG(id_node, i)));
  73     A_max = max(A, A_max);
  74     cell_normals[i] = GeometryTools::cellNormal(m_Grid, m_Part.n2cGG(id_node, i));
  75     cell_normals[i].normalise();
  76     bcs.insert(cell_code->GetValue(m_Part.n2cGG(id_node, i)));
  77   }
  78   int N = 0;
  79   for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
  80     double A = fabs(GeometryTools::cellVA(m_Grid, m_Part.n2cGG(id_node, i)));
  81     if (A > 0.01*A_max) {
  82       n += cell_normals[i];
  83       ++N;
  84     }
  85   }
  86   /*
  87   SurfaceProjection* proj = NULL;
  88   if (bcs.size() == 1) {
  89     proj = GuiMainWindow::pointer()->getSurfProj(*bcs.begin());
  90     if (proj) {
  91       proj->project(x_new, id_node, false, n);
  92       n = proj->lastProjNormal();
  93     }
  94   }
  95   */
  96   if (N == 0) {
  97     //EG_BUG;
  98     move = false;
  99   } else {
 100     n.normalise();
 101     double L_max = 0;// distance to the furthest neighbour node of id_node
 102     for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
 103       vec3_t xn;
 104       m_Grid->GetPoint(m_Part.n2nGG(id_node, i), xn.data());
 105       double L = (xn - x_old).abs();
 106       L_max = max(L, L_max);
 107     }
 108
 109     vec3_t x_summit;
 110     if(m_CorrectCurvature) {
 111       // better for mesher with interpolation
 112       x_summit = x_new + L_max*n;
 113     }
 114     else {
 115     // better for mesher without interpolation
 116       x_summit = x_old + L_max*n;
 117     }
 118
 119     for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
 120       vec3_t x[3];
 121       vtkIdType N_pts, *pts;
 122       m_Grid->GetCellPoints(m_Part.n2cGG(id_node, i), N_pts, pts);
 123       if (N_pts != 3) {
 124         //EG_BUG;
 125         move = false;
 126         break;
 127       }
 128       for (int j = 0; j < N_pts; ++j) {
 129         m_Grid->GetPoint(pts[j], x[j].data());
 130       }
 131       if (GeometryTools::tetraVol(x[0], x[1], x[2], x_summit, false) <= 0) {
 132         move = false;
 133         break;
 134       }
 135     }
 136   }
 137   /*
 138   if (move) {
 139     for (int i = 0; i < cell_normals.size(); ++i) {
 140       for (int j = 0; j < cell_normals.size(); ++j) {
 141         if (cell_normals[i]*cell_normals[j] < -100.0) {//Why -100.0?
 142           saveGrid(m_Grid, "after_move2");
 143           qWarning()<<"Cannot move point "<<id_node<<" because opposite cell_normals.";
 144           move = false;
 145           m_Grid->GetPoints()->SetPoint(id_node, x_old.data());
 146           saveGrid(m_Grid, "before_move2");
 147           EG_BUG;
 148           break;
 149         }
 150       }
 151     }
 152   }
 153   */
 154   if (!move) {
 155     m_Grid->GetPoints()->SetPoint(id_node, x_old.data());
 156   }
 157   return move;
 158 }
 159
 160 void LaplaceSmoother::featureCorrection(vtkIdType id_node, SurfaceProjection* proj, vec3_t &x_new)
 161 {
 162   if (m_FeatureMagic > 0) {
 163     EG_VTKDCN(vtkDoubleArray, cl, m_Grid, "node_meshdensity_desired");
 164     EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type");
 165     //if (node_type->GetValue(id_node) == VTK_FEATURE_EDGE_VERTEX) {
 166     {
 167       // "magic" vector to displace node for re-projection
 168       vec3_t magic_vector = m_NodeNormal[id_node];
 169
 170       vec3_t x0 = x_new;
 171       vec3_t x1 = proj->project(x_new, id_node, m_CorrectCurvature);
 172       vec3_t n = proj->lastProjNormal();
 173
 174       // check if mesh normal and projection normal are aligned
 175       // .. try to displace node slightly in order to get a proper projection normal
 176       //
 177       if (fabs(n*magic_vector) > 0.99) {
 178         x_new += 0.1*cl->GetValue(id_node)*magic_vector;
 179         x1 = proj->project(x_new, id_node, m_CorrectCurvature);
 180         n = proj->lastProjNormal();
 181       }
 182
 183       if (fabs(n*magic_vector) <= 0.99) {
 184
 185         // start the procedure if the vectors are not aligned
 186         //
 187         double L1 = 0;
 188         double L2 = m_FeatureMagic*cl->GetValue(id_node);
 189         for (int i = 0; i < 30; ++i) {
 190           x_new = x1 - 0.5*(L1 + L2)*magic_vector;
 191           x_new = proj->project(x_new, id_node, m_CorrectCurvature, n);
 192           double displacement = fabs((x_new - x1)*n);
 193           if (displacement > 0.01*cl->GetValue(id_node)) {
 194             L2 = 0.5*(L1 + L2);
 195           } else {
 196
 197             // if there is no significant displacement after the first iteration
 198             // the node is probably in a smooth region of the surface
 199             // ==> stop here
 200             //
 201             if (i == 0) {
 202               x_new = x0;
 203               break;
 204             }
 205
 206             L1 = 0.5*(L1 + L2);
 207           }
 208         }
 209         x_new = x1 - L1*magic_vector;
 210         x_new = proj->project(x_new, id_node, m_CorrectCurvature, n);
 211       } else {
 212
 213         // If they are still aligned it is an awkward situation.
 214         // .. The node might be in a corner already
 215         // .. skip iteration in this case
 216         //
 217         x_new = x0;
 218
 219       }
 220     }
 221   }
 222 }
 223
 224 bool LaplaceSmoother::moveNode(vtkIdType id_node, vec3_t &Dx)
 225 {
 226   if (!checkVector(Dx)) {
 227     return false;
 228   }
 229   EG_VTKDCN(vtkDoubleArray, cl, m_Grid, "node_meshdensity_desired");
 230   vec3_t x_old;
 231   m_Grid->GetPoint(id_node, x_old.data());
 232   bool moved = false;
 233   for (int i_relaxation = 0; i_relaxation < 1; ++i_relaxation) {
 234     vec3_t x_new = x_old + Dx;
 235     if (m_UseProjection) {
 236       int i_nodes = m_Part.localNode(id_node);
 237       if (m_NodeToBc[i_nodes].size() == 1 || m_UniformSnapPoints) {
 238         int bc = m_NodeToBc[i_nodes][0];
 239         x_new = GuiMainWindow::pointer()->getSurfProj(bc)->project(x_new, id_node, m_CorrectCurvature);
 240         featureCorrection(id_node, GuiMainWindow::pointer()->getSurfProj(bc), x_new);
 241       } else {
 242         for (int i_proj_iter = 0; i_proj_iter < m_ProjectionIterations; ++i_proj_iter) {
 243           foreach (int bc, m_NodeToBc[i_nodes]) {
 244             x_new = GuiMainWindow::pointer()->getSurfProj(bc)->project(x_new, id_node, m_CorrectCurvature);
 245           }
 246         }
 247
 248         for (int i_proj_iter = 0; i_proj_iter < m_ProjectionIterations; ++i_proj_iter) {
 249           if (m_CorrectCurvature) {
 250             foreach (int bc, m_NodeToBc[i_nodes]) {
 251               x_new = GuiMainWindow::pointer()->getSurfProj(bc)->correctCurvature(GuiMainWindow::pointer()->getSurfProj(bc)->lastProjTriangle(), x_new);
 252             }
 253           }
 254         }
 255
 256       }
 257     }
 258
 259     // compute the minimal length of any edge adjacent to this node
 260     // .. This will be used to limit the node movement.
 261     // .. Hopefully jammed topologies can be avoided this way.
 262     //
 263     EG_VTKDCN(vtkDoubleArray, cl, m_Grid, "node_meshdensity_desired");
 264     vec3_t x_old;
 265     m_Grid->GetPoint(id_node, x_old.data());
 266     double L_min = cl->GetValue(id_node);
 267     for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
 268       vtkIdType id_neigh = m_Part.n2nGG(id_node, i);
 269       vec3_t x_neigh;
 270       m_Grid->GetPoint(id_neigh, x_neigh.data());
 271       L_min = min(L_min, (x_old - x_neigh).abs());
 272     }
 273
 274     // limit node displacement
 275     vec3_t dx = x_new - x_old;
 276     if (dx.abs() > m_Limit*L_min) {
 277       x_new -= dx;
 278       dx.normalise();
 279       x_new += m_Limit*L_min*dx;
 280     }
 281
 282     if (setNewPosition(id_node, x_new)) {
 283       //if (m_NodeMovementCheck.moveNode(id_node, x_new)) {
 284       moved = true;
 285       Dx = x_new - x_old;
 286       break;
 287     }
 288     Dx *= 0.5;
 289   }
 290   return moved;
 291 }
 292
 293 void LaplaceSmoother::fixNodes(const QVector<bool> &fixnodes)
 294 {
 295   if (fixnodes.size() != m_Grid->GetNumberOfPoints()) {
 296     EG_BUG;
 297   }
 298   m_Fixed = fixnodes;
 299 }
 300
 301 void LaplaceSmoother::operate()
 302 {
 303   if (m_Fixed.size() != m_Grid->GetNumberOfPoints()) {
 304     m_Fixed.fill(false, m_Grid->GetNumberOfPoints());
 305   }
 306   QVector<int> bcs;
 307   GuiMainWindow::pointer()->getAllBoundaryCodes(bcs);
 308   if (m_UseProjection) {
 309     foreach (int bc, bcs) {
 310       GuiMainWindow::pointer()->getSurfProj(bc)->setForegroundGrid(m_Grid);
 311     }
 312   }
 313   UpdatePotentialSnapPoints(false, false);
 314   EG_VTKDCC(vtkIntArray,    cell_code, m_Grid, "cell_code");
 315   EG_VTKDCN(vtkCharArray,   node_type, m_Grid, "node_type" );
 316   EG_VTKDCN(vtkDoubleArray, cl,        m_Grid, "node_meshdensity_desired");
 317   QVector<vtkIdType> smooth_node(m_Grid->GetNumberOfPoints(), false);
 318   {
 319     l2g_t nodes = m_Part.getNodes();
 320     foreach (vtkIdType id_node, nodes) {
 321       smooth_node[id_node] = true;
 322     }
 323   }
 324   setAllSurfaceCells();
 325   l2g_t  nodes = m_Part.getNodes();
 326   m_NodeToBc.resize(nodes.size());
 327   for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
 328     QSet<int> bcs;
 329     for (int j = 0; j < m_Part.n2cLSize(i_nodes); ++j) {
 330       bcs.insert(cell_code->GetValue(m_Part.n2cLG(i_nodes, j)));
 331     }
 332     m_NodeToBc[i_nodes].resize(bcs.size());
 333     qCopy(bcs.begin(), bcs.end(), m_NodeToBc[i_nodes].begin());
 334   }
 335
 336   QVector<vec3_t> x_new(nodes.size());
 337
 338   QVector<bool> blocked(nodes.size(), false);
 339   for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
 340     for (int i = 0; i < m_Part.n2cLSize(i_nodes); ++i) {
 341       vtkIdType type = m_Grid->GetCellType(m_Part.n2cLG(i_nodes, i));
 342       if (!m_AllowedCellTypes.contains(type)) {
 343         blocked[i_nodes] = true;
 344         break;
 345       }
 346     }
 347   }
 348
 349   for (int i_iter = 0; i_iter < m_NumberOfIterations; ++i_iter) {
 350     m_Success = true;
 351     computeNormals();
 352     for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
 353       vtkIdType id_node = nodes[i_nodes];
 354       if (!m_Fixed[id_node] && !blocked[i_nodes]) {
 355         if (smooth_node[id_node] && node_type->GetValue(id_node) != VTK_FIXED_VERTEX) {
 356           if (node_type->GetValue(id_node) != VTK_FIXED_VERTEX) {
 357             QVector<vtkIdType> snap_points = getPotentialSnapPoints(id_node);
 358             vec3_t n(0,0,0);
 359             if (snap_points.size() > 0) {
 360               vec3_t x_old;
 361               vec3_t x;
 362               x_new[i_nodes] = vec3_t(0,0,0);
 363               m_Grid->GetPoint(id_node, x_old.data());
 364               double w_tot = 0;
 365               double L_min = 1e99;
 366               foreach (vtkIdType id_snap_node, snap_points) {
 367                 m_Grid->GetPoint(id_snap_node, x.data());
 368                 double w = 1.0;
 369                 w_tot += w;
 370                 x_new[i_nodes] += w*x;
 371                 n += m_NodeNormal[id_snap_node];
 372                 double L = (x - x_old).abs();
 373                 L_min = min(L, L_min);
 374               }
 375               n.normalise();
 376               x_new[i_nodes] *= 1.0/w_tot;
 377
 378               if (m_UseNormalCorrection) {
 379                 vec3_t dx = x_new[i_nodes] - x_old;
 380                 double scal = dx*m_NodeNormal[id_node];
 381                 x_new[i_nodes] += scal*m_NodeNormal[id_node];
 382               }
 383               vec3_t Dx = x_new[i_nodes] - x_old;
 384               Dx *= m_UnderRelaxation;
 385               if (moveNode(id_node, Dx)) {
 386                 x_new[i_nodes] = x_old + Dx;
 387               } else {
 388                 x_new[i_nodes] = x_old;
 389                 m_Success = false;
 390               }
 391             }
 392           }
 393         }
 394       }
 395     }
 396     if (m_Success) {
 397       break;
 398     }
 399   }
 400 }