src/laplacesmoother.cpp

   1 //
   2 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   3 // +                                                                      +
   4 // + This file is part of enGrid.                                         +
   5 // +                                                                      +
   6 // + Copyright 2008-2010 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 #include "guimainwindow.h"
  28
  29 using namespace GeometryTools;
  30
  31 LaplaceSmoother::LaplaceSmoother() : SurfaceOperation()
  32 {
  33   DebugLevel = 0;
  34   setQuickSave(true);
  35   m_UseProjection = true;
  36 //   m_UseNormalCorrection = false;
  37   getSet("surface meshing", "under relaxation for smoothing", 0.5, m_UnderRelaxation);
  38   getSet("surface meshing", "correct curvature (experimental)", false, m_correctCurvature);
  39   m_NoCheck = false;
  40   m_ProjectionIterations = 20;
  41 }
  42
  43 bool LaplaceSmoother::setNewPosition(vtkIdType id_node, vec3_t x_new)
  44 {
  45   using namespace GeometryTools;
  46
  47   vec3_t x_old;
  48   m_Grid->GetPoint(id_node, x_old.data());
  49   m_Grid->GetPoints()->SetPoint(id_node, x_new.data());
  50
  51   bool move = true;
  52   if(m_NoCheck) return move;
  53
  54   vec3_t n(0,0,0);
  55   QVector<vec3_t> cell_normals(m_Part.n2cGSize(id_node));
  56   double A_max = 0;//area of the biggest neighbour cell of id_node
  57   for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
  58     double A = fabs(GeometryTools::cellVA(m_Grid, m_Part.n2cGG(id_node, i)));
  59     A_max = max(A, A_max);
  60     cell_normals[i] = GeometryTools::cellNormal(m_Grid, m_Part.n2cGG(id_node, i));
  61     cell_normals[i].normalise();
  62   }
  63   int N = 0;
  64   for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
  65     double A = fabs(GeometryTools::cellVA(m_Grid, m_Part.n2cGG(id_node, i)));
  66     if (A > 0.01*A_max) {
  67       n += cell_normals[i];
  68       ++N;
  69     }
  70   }
  71   if (N == 0) {
  72     EG_BUG;
  73     move = false;
  74   } else {
  75     n.normalise();
  76     double L_max = 0;// distance to the furthest neighbour node of id_node
  77     for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
  78       vec3_t xn;
  79       m_Grid->GetPoint(m_Part.n2nGG(id_node, i), xn.data());
  80       double L = (xn - x_old).abs();
  81       L_max = max(L, L_max);
  82     }
  83
  84     vec3_t x_summit;
  85     if(m_correctCurvature) {
  86       // better for mesher with interpolation
  87       x_summit = x_new + L_max*n;
  88     }
  89     else {
  90     // better for mesher without interpolation
  91       x_summit = x_old + L_max*n;
  92     }
  93
  94     for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
  95       vec3_t x[3];
  96       vtkIdType N_pts, *pts;
  97       m_Grid->GetCellPoints(m_Part.n2cGG(id_node, i), N_pts, pts);
  98       if (N_pts != 3) {
  99         EG_BUG;
 100       }
 101       for (int j = 0; j < N_pts; ++j) {
 102         m_Grid->GetPoint(pts[j], x[j].data());
 103       }
 104       if (GeometryTools::tetraVol(x[0], x[1], x[2], x_summit, false) <= 0) {
 105 //         saveGrid(m_Grid, "after_move1");
 106 //         qWarning()<<"Cannot move point "<<id_node<<" because negative tetraVol.";
 107         move = false;
 108 //         m_Grid->GetPoints()->SetPoint(id_node, x_old.data());
 109 //         saveGrid(m_Grid, "before_move1");
 110 //         EG_BUG;
 111         break;
 112       }
 113     }
 114   }
 115   if (move) {
 116     for (int i = 0; i < cell_normals.size(); ++i) {
 117       for (int j = 0; j < cell_normals.size(); ++j) {
 118         if (cell_normals[i]*cell_normals[j] < -1000*0.1) {//Why -1000*0.1?
 119           saveGrid(m_Grid, "after_move2");
 120           qWarning()<<"Cannot move point "<<id_node<<" because opposite cell_normals.";
 121           move = false;
 122           m_Grid->GetPoints()->SetPoint(id_node, x_old.data());
 123           saveGrid(m_Grid, "before_move2");
 124           EG_BUG;
 125           break;
 126         }
 127       }
 128     }
 129   }
 130
 131   if (!move) {
 132     // comment this out if you want points to always move
 133     m_Grid->GetPoints()->SetPoint(id_node, x_old.data());
 134 //     saveGrid(m_Grid, "before_move");
 135   }
 136   return move;
 137 }
 138
 139 bool LaplaceSmoother::moveNode(vtkIdType id_node, vec3_t &Dx)
 140 {
 141   vec3_t x_old;
 142   m_Grid->GetPoint(id_node, x_old.data());
 143   bool moved = false;
 144   for (int i_relaxation = 0; i_relaxation < 1; ++i_relaxation) {
 145     vec3_t x_new = x_old + Dx;
 146     if (m_UseProjection) {
 147       int i_nodes = m_Part.localNode(id_node);
 148       if (m_NodeToBc[i_nodes].size() == 1) {
 149         int bc = m_NodeToBc[i_nodes][0];
 150         x_new = GuiMainWindow::pointer()->getSurfProj(bc)->project(x_new, id_node);
 151       } else {
 152         for (int i_proj_iter = 0; i_proj_iter < m_ProjectionIterations; ++i_proj_iter) {
 153           foreach (int bc, m_NodeToBc[i_nodes]) {
 154             x_new = GuiMainWindow::pointer()->getSurfProj(bc)->project(x_new, id_node);
 155           }
 156         }
 157       }
 158     }
 159     if (setNewPosition(id_node, x_new)) {
 160       moved = true;
 161       Dx = x_new - x_old;
 162       break;
 163     }
 164     Dx *= 0.5;
 165   }
 166   return moved;
 167 }
 168
 169
 170 void LaplaceSmoother::operate()
 171 {
 172 //   qDebug()<<"LaplaceSmoother::operate() called";
 173
 174   QSet<int> bcs;
 175   GuiMainWindow::pointer()->getAllBoundaryCodes(bcs);
 176   if (m_UseProjection) {
 177     foreach (int bc, bcs) {
 178       GuiMainWindow::pointer()->getSurfProj(bc)->setForegroundGrid(m_Grid);
 179       GuiMainWindow::pointer()->getSurfProj(bc)->setCorrectCurvature(m_correctCurvature);
 180     }
 181   }
 182   UpdatePotentialSnapPoints(false, false);
 183   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
 184   EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type" );
 185   QVector<vtkIdType> smooth_node(m_Grid->GetNumberOfPoints(), false);
 186   {
 187     l2g_t nodes = m_Part.getNodes();
 188     foreach (vtkIdType id_node, nodes) {
 189       smooth_node[id_node] = true;
 190     }
 191   }
 192   setAllSurfaceCells();
 193   l2g_t  nodes = m_Part.getNodes();
 194   g2l_t _nodes = m_Part.getLocalNodes();
 195   m_NodeToBc.resize(nodes.size());
 196   for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
 197     QSet<int> bcs;
 198     for (int j = 0; j < m_Part.n2cLSize(i_nodes); ++j) {
 199       bcs.insert(cell_code->GetValue(m_Part.n2cLG(i_nodes, j)));
 200     }
 201     m_NodeToBc[i_nodes].resize(bcs.size());
 202     qCopy(bcs.begin(), bcs.end(), m_NodeToBc[i_nodes].begin());
 203   }
 204
 205   QVector<vec3_t> x_new(nodes.size());
 206
 207   for (int i_iter = 0; i_iter < m_NumberOfIterations; ++i_iter) {
 208
 209     m_Success = true;
 210
 211 //     QVector<vec3_t> node_normals;
 212 //     if (m_UseNormalCorrection) {
 213 //       node_normals.fill(vec3_t(0,0,0), m_Grid->GetNumberOfPoints());
 214 //       vec3_t n(0,0,0);
 215 //       for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 216 //         for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
 217 //           node_normals[id_node] += GeometryTools::cellNormal(m_Grid, m_Part.n2cGG(id_node, i));
 218 //         }
 219 //         node_normals[id_node].normalise();
 220 //       }
 221 //     }
 222
 223     for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
 224       vtkIdType id_node = nodes[i_nodes];
 225       if (smooth_node[id_node] && node_type->GetValue(id_node) != VTK_FIXED_VERTEX) {
 226         if (node_type->GetValue(id_node) != VTK_FIXED_VERTEX) {
 227           QVector<vtkIdType> snap_points = getPotentialSnapPoints(id_node);
 228 //           vec3_t n(0,0,0);
 229           if (snap_points.size() > 0) {
 230             vec3_t x_old;
 231             vec3_t x;
 232             x_new[i_nodes] = vec3_t(0,0,0);
 233             m_Grid->GetPoint(id_node, x_old.data());
 234             foreach (vtkIdType id_snap_node, snap_points) {
 235               m_Grid->GetPoint(id_snap_node, x.data());
 236               x_new[i_nodes] += x;
 237 //               n += node_normals[id_snap_node];
 238             }
 239 //             n.normalise();
 240             x_new[i_nodes] *= 1.0/snap_points.size();
 241
 242 //             if (m_UseNormalCorrection) {
 243 //               vec3_t dx = x_new[i_nodes] - x_old;
 244 //               dx = (dx*n)*n;
 245 //               x_new[i_nodes] -= dx;
 246 //             }
 247
 248             vec3_t Dx = x_new[i_nodes] - x_old;
 249             Dx *= m_UnderRelaxation;
 250             if (moveNode(id_node, Dx)) {
 251               x_new[i_nodes] = x_old + Dx;
 252             } else {
 253               x_new[i_nodes] = x_old;
 254               m_Success = false;
 255             }
 256           }
 257         }
 258       }
 259     }
 260     if (m_Success) {
 261       break;
 262     }
 263   }
 264 }