src/libengrid/laplacesmoother.cpp

   1 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   2 // +                                                                      +
   3 // + This file is part of enGrid.                                         +
   4 // +                                                                      +
   5 // + Copyright 2008-2014 enGits GmbH                                      +
   6 // +                                                                      +
   7 // + enGrid is free software: you can redistribute it and/or modify       +
   8 // + it under the terms of the GNU General Public License as published by +
   9 // + the Free Software Foundation, either version 3 of the License, or    +
  10 // + (at your option) any later version.                                  +
  11 // +                                                                      +
  12 // + enGrid is distributed in the hope that it will be useful,            +
  13 // + but WITHOUT ANY WARRANTY; without even the implied warranty of       +
  14 // + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the        +
  15 // + GNU General Public License for more details.                         +
  16 // +                                                                      +
  17 // + You should have received a copy of the GNU General Public License    +
  18 // + along with enGrid. If not, see <http://www.gnu.org/licenses/>.       +
  19 // +                                                                      +
  20 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  21 #include "laplacesmoother.h"
  22 #include <vtkCellLocator.h>
  23 #include <vtkCharArray.h>
  24 #include <vtkGenericCell.h>
  25
  26 #include "guimainwindow.h"
  27 #include "localnodegraphinterface.h"
  28 #include "checkerboardgraphiterator.h"
  29 #include "snaptofeatures.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", 1.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   //m_StrictFeatureSnap = false;
  49 }
  50
  51 bool LaplaceSmoother::setNewPosition(vtkIdType id_node, vec3_t x_new)
  52 {
  53   using namespace GeometryTools;
  54
  55   vec3_t x_old;
  56   m_Grid->GetPoint(id_node, x_old.data());
  57   EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type");
  58   bool move = true;
  59   if (m_NoCheck) {
  60     m_Grid->GetPoints()->SetPoint(id_node, x_new.data());
  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(vtkCharArray, node_type, m_Grid, "node_type");
  86   char type = node_type->GetValue(id_node);
  87   if (type == EG_FEATURE_EDGE_VERTEX || type == EG_BOUNDARY_EDGE_VERTEX) {
  88     x_new = cad_interface->snapToEdge(x_new);
  89   }
  90 }
  91
  92 bool LaplaceSmoother::moveNode(vtkIdType id_node, vec3_t &Dx)
  93 {
  94   if (!checkVector(Dx)) {
  95     return false;
  96   }
  97   EG_VTKDCN(vtkDoubleArray, cl, m_Grid, "node_meshdensity_desired");
  98   vec3_t x_old;
  99   m_Grid->GetPoint(id_node, x_old.data());
 100   bool moved = false;
 101
 102   for (int i_relaxation = 0; i_relaxation < 10; ++i_relaxation) {
 103     vec3_t x_new = x_old + Dx;
 104     if (m_UseProjection) {
 105       int i_nodes = m_Part.localNode(id_node);
 106       int bc = m_NodeToBc[i_nodes][0];
 107       x_new = GuiMainWindow::pointer()->getCadInterface(bc)->snapNode(id_node, x_new, m_CorrectCurvature);
 108       featureCorrection(id_node, GuiMainWindow::pointer()->getCadInterface(bc), x_new);
 109     }
 110
 111     // compute the minimal length of any edge adjacent to this node
 112     // .. This will be used to limit the node movement.
 113     // .. Hopefully jammed topologies can be avoided this way.
 114     //
 115     EG_VTKDCN(vtkDoubleArray, cl, m_Grid, "node_meshdensity_desired");
 116     double L_min = cl->GetValue(id_node);
 117     for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
 118       vtkIdType id_neigh = m_Part.n2nGG(id_node, i);
 119       vec3_t x_neigh;
 120       m_Grid->GetPoint(id_neigh, x_neigh.data());
 121       L_min = min(L_min, (x_old - x_neigh).abs());
 122     }
 123
 124     if (setNewPosition(id_node, x_new)) {
 125       m_Grid->GetPoints()->SetPoint(id_node, x_new.data());
 126       moved = true;
 127       Dx = x_new - x_old;
 128       break;
 129     }
 130     Dx *= 0.5;
 131   }
 132   return moved;
 133 }
 134
 135 void LaplaceSmoother::fixNodes(const QVector<bool> &fixnodes)
 136 {
 137   if (fixnodes.size() != m_Grid->GetNumberOfPoints()) {
 138     EG_BUG;
 139   }
 140   m_Fixed = fixnodes;
 141 }
 142
 143 void LaplaceSmoother::operate()
 144 {
 145   m_NoCheck = false; // !!!!!!!!!!!!!!
 146   if (m_Fixed.size() != m_Grid->GetNumberOfPoints()) {
 147     m_Fixed.fill(false, m_Grid->GetNumberOfPoints());
 148   }
 149   QVector<int> bcs;
 150   GuiMainWindow::pointer()->getAllBoundaryCodes(bcs);
 151   if (m_UseProjection) {
 152     foreach (int bc, bcs) {
 153       GuiMainWindow::pointer()->getCadInterface(bc)->setForegroundGrid(m_Grid);
 154     }
 155   }
 156   updateNodeInfo();
 157   EG_VTKDCC(vtkIntArray,    cell_code, m_Grid, "cell_code");
 158   EG_VTKDCN(vtkCharArray,   node_type, m_Grid, "node_type" );
 159   EG_VTKDCN(vtkDoubleArray, cl,        m_Grid, "node_meshdensity_desired");
 160   QVector<vtkIdType> smooth_node(m_Grid->GetNumberOfPoints(), false);
 161   {
 162     l2g_t nodes = m_Part.getNodes();
 163     foreach (vtkIdType id_node, nodes) {
 164       smooth_node[id_node] = true;
 165     }
 166   }
 167   setAllSurfaceCells();
 168   l2g_t  nodes = m_Part.getNodes();
 169   m_NodeToBc.resize(nodes.size());
 170   for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
 171     QSet<int> bcs;
 172     for (int j = 0; j < m_Part.n2cLSize(i_nodes); ++j) {
 173       bcs.insert(cell_code->GetValue(m_Part.n2cLG(i_nodes, j)));
 174     }
 175     m_NodeToBc[i_nodes].resize(bcs.size());
 176     qCopy(bcs.begin(), bcs.end(), m_NodeToBc[i_nodes].begin());
 177   }
 178
 179   QVector<vec3_t> x_new(nodes.size());
 180
 181   QVector<bool> blocked(nodes.size(), false);
 182   for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
 183     if (node_type->GetValue(nodes[i_nodes]) == EG_FEATURE_CORNER_VERTEX) {
 184       blocked[i_nodes] = true;
 185     } else {
 186       for (int i = 0; i < m_Part.n2cLSize(i_nodes); ++i) {
 187         vtkIdType type = m_Grid->GetCellType(m_Part.n2cLG(i_nodes, i));
 188         if (!m_AllowedCellTypes.contains(type)) {
 189           blocked[i_nodes] = true;
 190           break;
 191         }
 192       }
 193     }
 194   }
 195
 196   for (int i_iter = 0; i_iter < m_NumberOfIterations; ++i_iter) {
 197     m_Success = true;
 198     computeNormals();
 199
 200     for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
 201       vtkIdType id_node = nodes[i_nodes];
 202       m_Grid->GetPoint(id_node, x_new[i_nodes].data());
 203       if (!m_Fixed[id_node] && !blocked[i_nodes]) {
 204         if (smooth_node[id_node] && node_type->GetValue(id_node) != EG_FIXED_VERTEX) {
 205           if (node_type->GetValue(id_node) != EG_FIXED_VERTEX) {
 206             QVector<vtkIdType> snap_points = getPotentialSnapPoints(id_node);
 207             vec3_t n(0,0,0);
 208             vec3_t x_old = x_new[i_nodes];
 209
 210             if (snap_points.size() > 0) {
 211               vec3_t x;
 212               x_new[i_nodes] = vec3_t(0,0,0);
 213               double w_tot = 0;
 214               double L_min = 1e99;
 215               foreach (vtkIdType id_snap_node, snap_points) {
 216                 m_Grid->GetPoint(id_snap_node, x.data());
 217                 double w = 1.0;
 218                 w_tot += w;
 219                 x_new[i_nodes] += w*x;
 220                 n += m_NodeNormal[id_snap_node];
 221                 double L = (x - x_old).abs();
 222                 L_min = min(L, L_min);
 223               }
 224               n.normalise();
 225               x_new[i_nodes] *= 1.0/w_tot;
 226             } else {
 227               x_new[i_nodes] = x_old;
 228             }
 229
 230             if (m_UseNormalCorrection) {
 231               vec3_t dx = x_new[i_nodes] - x_old;
 232               double scal = dx*m_NodeNormal[id_node];
 233               x_new[i_nodes] += scal*m_NodeNormal[id_node];
 234             }
 235             vec3_t Dx = x_new[i_nodes] - x_old;
 236             //Dx *= m_UnderRelaxation;
 237             if (moveNode(id_node, Dx)) {
 238               x_new[i_nodes] = x_old + m_UnderRelaxation*Dx;
 239             } else {
 240               x_new[i_nodes] = x_old;
 241               //m_Success = false;
 242             }
 243             m_Grid->GetPoints()->SetPoint(id_node, x_old.data());
 244           }
 245         }
 246       }
 247     }
 248     for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
 249       vtkIdType id_node = nodes[i_nodes];
 250       m_Grid->GetPoints()->SetPoint(id_node, x_new[i_nodes].data());
 251     }
 252 /*
 253     SnapToFeatures snap;
 254     snap();
 255 */
 256     if (m_Success) {
 257       break;
 258     }
 259   }
 260 }