src/libengrid/gridsmoother.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 "gridsmoother.h"
  22 #include "engrid.h"
  23 #include "guimainwindow.h"
  24 #include "elements.h"
  25 #include "optimisenormalvector.h"
  26 #include "pointfinder.h"
  27
  28 #include <QTime>
  29
  30 GridSmoother::GridSmoother()
  31 {
  32   m_NumIterations          = 5;
  33   m_NumRelaxations         = 5;
  34   m_NumBoundaryCorrections = 50;
  35   m_DesiredStretching      = 1.2;
  36   m_FirstCall = true;
  37
  38   getSet("boundary layer", "number of smoothing sub-iterations",       5,     m_NumIterations);
  39   getSet("boundary layer", "use strict prism checking",                false, m_StrictPrismChecking);
  40   getSet("boundary layer", "number of normal vector relax iterations", 10,    m_NumNormalRelaxations);
  41   getSet("boundary layer", "number of layer height relax iterations",  3,     m_NumHeightRelaxations);
  42   getSet("boundary layer", "radar angle",                              45,    m_RadarAngle);
  43   getSet("boundary layer", "maximal layer height in gaps",             0.2,   m_MaxHeightInGaps);
  44   getSet("boundary layer", "relative face size (lower limit)",         0.5,   m_FaceSizeLowerLimit);
  45   getSet("boundary layer", "relative face size (upper limit)",         2.0,   m_FaceSizeUpperLimit);
  46   getSet("boundary layer", "angle between top and bottom face",        45.0,  m_FaceAngleLimit);
  47
  48   m_FaceAngleLimit = deg2rad(m_FaceAngleLimit);
  49
  50   //m_CritAngle = GeometryTools::deg2rad(m_CritAngle);
  51 }
  52
  53 void GridSmoother::markNodes()
  54 {
  55   m_NodeMarked.fill(false,m_Grid->GetNumberOfPoints());
  56   QVector<bool> new_mark(m_Grid->GetNumberOfPoints());
  57   for (int i_iterations = 0; i_iterations < 1; ++i_iterations) {
  58     qCopy(m_NodeMarked.begin(),m_NodeMarked.end(),new_mark.begin());
  59     for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
  60       bool mark_cell = false;
  61       EG_GET_CELL(id_cell, m_Grid);
  62       if (type_cell == VTK_WEDGE) {
  63         mark_cell = true;
  64       } else {
  65         for (int i_pts = 0; i_pts < num_pts; ++i_pts) {
  66           if (m_NodeMarked[pts[i_pts]]) {
  67             mark_cell = true;
  68           }
  69         }
  70       }
  71       if (mark_cell) {
  72         for (int i_pts = 0; i_pts < num_pts; ++i_pts) {
  73           new_mark[pts[i_pts]] = true;
  74         }
  75       }
  76     }
  77     qCopy(new_mark.begin(),new_mark.end(),m_NodeMarked.begin());
  78   }
  79   QSet<int> free_bcs = m_BoundaryCodes + m_LayerAdjacentBoundaryCodes;
  80   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
  81   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
  82     if (isSurface(id_cell, m_Grid)) {
  83       if (!free_bcs.contains(cell_code->GetValue(id_cell))) {
  84         EG_GET_CELL(id_cell, m_Grid);
  85         for (int i_pts = 0; i_pts < num_pts; ++i_pts) {
  86           m_NodeMarked[pts[i_pts]] = false;
  87         }
  88       }
  89     }
  90   }
  91   m_NumMarkedNodes = 0;
  92   QVector<vtkIdType> nodes = m_Part.getNodes();
  93   foreach (vtkIdType id_node, nodes) {
  94     if (id_node < 0) EG_BUG;
  95     if (id_node > m_Grid->GetNumberOfPoints()) EG_BUG;
  96     if (m_NodeMarked[id_node]) {
  97       ++m_NumMarkedNodes;
  98     }
  99   }
 100 }
 101
 102 bool GridSmoother::setNewPosition(vtkIdType id_node, vec3_t x_new)
 103 {
 104   using namespace GeometryTools;
 105
 106   vec3_t x_old;
 107   m_Grid->GetPoint(id_node, x_old.data());
 108   m_Grid->GetPoints()->SetPoint(id_node, x_new.data());
 109   bool move = true;
 110
 111   if (move) {
 112     Elements E;
 113
 114     l2g_t cells = m_Part.getCells();
 115     l2l_t n2c = m_Part.getN2C();
 116
 117     foreach (int i_cells, n2c[id_node]) {
 118       vtkIdType id_cell = cells[i_cells];
 119       vtkIdType type_cell = m_Grid->GetCellType(id_cell);
 120
 121       if (type_cell == VTK_TRIANGLE) {
 122         vec3_t x[3];
 123         EG_GET_CELL(id_cell, m_Grid);
 124         for (int i = 0; i < 3; ++i) {
 125           m_Grid->GetPoint(pts[i], x[i].data());
 126         }
 127         double L_max = 0;
 128         for (int i = 0; i < 3; ++i) {
 129           for (int j = 0; j < 3; ++j) {
 130             if (i != j) {
 131               L_max = max(L_max, (x[i]-x[j]).abs());
 132             }
 133           }
 134         }
 135         double A = GeometryTools::cellVA(m_Grid, id_cell);
 136         if (A < 1e-3*L_max*L_max) {
 137           move = false;
 138         } else if (m_NodeNormal[id_node].abs() > 0.1) {
 139           vec3_t n = GeometryTools::triNormal(x[0], x[1], x[2]);
 140           n.normalise();
 141           if (n*m_NodeNormal[id_node] < 0.1) {
 142             move = false;
 143           }
 144         }
 145       }
 146
 147       if (type_cell == VTK_TETRA) {
 148         vec3_t x[4];
 149         EG_GET_CELL(id_cell, m_Grid);
 150         for (int i = 0; i < 4; ++i) {
 151           m_Grid->GetPoint(pts[i], x[i].data());
 152         }
 153         double L_max = 0;
 154         for (int i = 0; i < 4; ++i) {
 155           for (int j = 0; j < 4; ++j) {
 156             if (i != j) {
 157               L_max = max(L_max, (x[i]-x[j]).abs());
 158             }
 159           }
 160         }
 161         if (GeometryTools::cellVA(m_Grid, id_cell) < 1e-3*L_max*L_max*L_max) {
 162           move = false;
 163         }
 164       }
 165
 166       if (type_cell == VTK_WEDGE && m_StrictPrismChecking) {
 167         vec3_t xtet[4];
 168         EG_GET_CELL(id_cell, m_Grid);
 169         bool ok = true;
 170         for (int i = 0; i < 4; ++i) {     // variation
 171           ok = true;
 172           for (int j = 0; j < 3; ++j) {   // tetrahedron
 173             for (int k = 0; k < 4; ++k) { // node
 174               m_Grid->GetPoint(pts[E.priTet(i,j,k)], xtet[k].data());
 175             }
 176             double V = GeometryTools::tetraVol(xtet[0], xtet[1], xtet[2], xtet[3]);
 177             if (V <= 0) {
 178               ok = false;
 179             }
 180           }
 181           if (ok) {
 182             break;
 183           }
 184         }
 185         if (!ok) {
 186           move = false;
 187         }
 188       }
 189     }
 190   }
 191   if (!move) {
 192     m_Grid->GetPoints()->SetPoint(id_node, x_old.data());
 193   }
 194   return move;
 195 }
 196
 197 void GridSmoother::correctDx(int i_nodes, vec3_t &Dx)
 198 {
 199   l2g_t nodes = m_Part.getNodes();
 200   l2l_t n2c = m_Part.getN2C();
 201   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
 202   vec3_t x_old;
 203   m_Grid->GetPoint(nodes[i_nodes],x_old.data());
 204   for (int i_boundary_correction = 0; i_boundary_correction < m_NumBoundaryCorrections; ++i_boundary_correction) {
 205     foreach (vtkIdType id_cell, n2c[i_nodes]) {
 206       if (isSurface(id_cell, m_Grid)) {
 207         int bc = cell_code->GetValue(id_cell);
 208         vec3_t x_new = x_old + Dx;
 209         x_new = GuiMainWindow::pointer()->getCadInterface(bc)->snapNode(nodes[i_nodes], x_new, false);
 210         Dx = x_new - x_old;
 211       } else {
 212         if (m_Grid->GetCellType(id_cell) == VTK_WEDGE) {
 213           EG_GET_CELL(id_cell, m_Grid);
 214           vtkIdType id_surf_node = -1;
 215           if (pts[3] == nodes[i_nodes]) id_surf_node = pts[0];
 216           if (pts[4] == nodes[i_nodes]) id_surf_node = pts[1];
 217           if (pts[5] == nodes[i_nodes]) id_surf_node = pts[2];
 218           if (id_surf_node != -1) {
 219             vec3_t x0,x1,x2;
 220             m_Grid->GetPoint(pts[0],x0.data());
 221             m_Grid->GetPoint(pts[1],x1.data());
 222             m_Grid->GetPoint(pts[2],x2.data());
 223             vec3_t a = x1-x0;
 224             vec3_t b = x2-x0;
 225             vec3_t c = b-a;
 226             double L = (a.abs()+b.abs()+c.abs())/3.0;
 227             vec3_t n = b.cross(a);
 228             n.normalise();
 229             vec3_t x_old;
 230             m_Grid->GetPoint(nodes[i_nodes],x_old.data());
 231             vec3_t x_new = x_old + Dx - x0;
 232             if ( (n*x_new) <= 0 ) {
 233               x_new -= (x_new*n)*n;
 234               x_new += 1e-4*L*n;
 235               x_new += x0;
 236               Dx = x_new - x_old;
 237             }
 238           }
 239         }
 240       }
 241     }
 242   }
 243 }
 244
 245 bool GridSmoother::moveNode(int i_nodes, vec3_t &Dx)
 246 {
 247   m_CollisionDetected = false;
 248   l2g_t nodes = m_Part.getNodes();
 249   vtkIdType id_node = nodes[i_nodes];
 250   vec3_t x_old;
 251   m_Grid->GetPoint(id_node, x_old.data());
 252   bool moved = false;
 253   for (int i_relaxation = 0; i_relaxation < m_NumRelaxations; ++i_relaxation) {
 254     if (setNewPosition(id_node, x_old + Dx)) {
 255       moved = true;
 256       break;
 257     }
 258     Dx *= 0.5;
 259     correctDx(i_nodes, Dx);
 260   }
 261   return moved;
 262 }
 263
 264 void GridSmoother::computeNormals()
 265 {
 266   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
 267   m_NodeNormal.fill(vec3_t(0,0,0), m_Grid->GetNumberOfPoints());
 268   QVector<int> num_bcs(m_Grid->GetNumberOfPoints());
 269   QVector<OptimiseNormalVector> n_opt(m_Grid->GetNumberOfPoints());
 270   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 271     QSet<int> bcs;
 272     for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
 273       vtkIdType id_cell = m_Part.n2cGG(id_node, i);
 274       if (isSurface(id_cell, m_Grid)) {
 275         int bc = cell_code->GetValue(id_cell);
 276         if (m_BoundaryCodes.contains(bc)) {
 277           bcs.insert(bc);
 278         }
 279       }
 280     }
 281     num_bcs[id_node] = bcs.size();
 282     QVector<vec3_t> normal(num_bcs[id_node], vec3_t(0,0,0));
 283     QMap<int,int> bcmap;
 284     int i_bc = 0;
 285     foreach (int bc, bcs) {
 286       bcmap[bc] = i_bc;
 287       ++i_bc;
 288     }
 289     for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
 290       vtkIdType id_cell = m_Part.n2cGG(id_node, i);
 291       if (isSurface(id_cell, m_Grid)) {
 292         int bc = cell_code->GetValue(id_cell);
 293         EG_GET_CELL(id_cell, m_Grid);
 294         vec3_t a, b, c;
 295         for (int j = 0; j < num_pts; ++j) {
 296           if (pts[j] == id_node) {
 297             m_Grid->GetPoint(pts[j], a.data());
 298             if (j > 0) {
 299               m_Grid->GetPoint(pts[j-1], b.data());
 300             } else {
 301               m_Grid->GetPoint(pts[num_pts-1], b.data());
 302             }
 303             if (j < num_pts - 1) {
 304               m_Grid->GetPoint(pts[j+1], c.data());
 305             } else {
 306               m_Grid->GetPoint(pts[0], c.data());
 307             }
 308           }
 309         }
 310         vec3_t u = b - a;
 311         vec3_t v = c - a;
 312         double alpha = GeometryTools::angle(u, v);
 313         vec3_t n = u.cross(v);
 314         n.normalise();
 315         if (m_BoundaryCodes.contains(bc)) {
 316           normal[bcmap[bc]] += alpha*n;
 317           n_opt[id_node].addFace(n);
 318         } else {
 319           n_opt[id_node].addConstraint(n);
 320         }
 321       }
 322     }
 323     for (int i = 0; i < num_bcs[id_node]; ++i) {
 324       normal[i].normalise();
 325     }
 326     if (num_bcs[id_node] > 0) {
 327       if (num_bcs[id_node] > 1) {
 328         if (num_bcs[id_node] == 3) {
 329           for (int i = 0; i < num_bcs[id_node]; ++i) {
 330             for (int j = i + 1; j < num_bcs[id_node]; ++j) {
 331               vec3_t n = normal[i] + normal[j];
 332               n.normalise();
 333               m_NodeNormal[id_node] += n;
 334             }
 335           }
 336         } else {
 337           for (int i = 0; i < num_bcs[id_node]; ++i) {
 338             m_NodeNormal[id_node] += normal[i];
 339           }
 340         }
 341       } else {
 342         m_NodeNormal[id_node] = normal[0];
 343       }
 344       m_NodeNormal[id_node].normalise();
 345       m_NodeNormal[id_node] = n_opt[id_node](m_NodeNormal[id_node]);
 346       m_NodeNormal[id_node].normalise();
 347     }
 348   }
 349
 350   m_GeoNormal = m_NodeNormal;
 351   relaxNormalVectors();
 352
 353   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 354     if (m_NodeNormal[id_node].abs() < 0.1) {
 355       vec3_t n(0,0,0);
 356       int N = 0;
 357       for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
 358         vtkIdType id_cell = m_Part.n2cGG(id_node, i);
 359         if (isSurface(id_cell, m_Grid)) {
 360           n += GeometryTools::cellNormal(m_Grid, id_cell);
 361           ++N;
 362         }
 363       }
 364       if (N) {
 365         n.normalise();
 366         m_NodeNormal[id_node] = n;
 367       }
 368       if (num_bcs[id_node] > 1) {
 369         m_NodeNormal[id_node] = n_opt[id_node](m_NodeNormal[id_node]);
 370       }
 371     }
 372   }
 373 }
 374
 375 void GridSmoother::relaxNormalVectors()
 376 {
 377   m_SurfNode.fill(false, m_Grid->GetNumberOfPoints());
 378   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 379     for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
 380       if (isSurface(m_Part.n2cGG(id_node,i), m_Grid)) {
 381         m_SurfNode[id_node] = true;
 382         break;
 383       }
 384     }
 385   }
 386   EG_VTKDCC(vtkIntArray, bc, m_Grid, "cell_code");
 387   QVector<QSet<int> > n2bc(m_Grid->GetNumberOfPoints());
 388   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 389     if (isSurface(id_cell, m_Grid)) {
 390       EG_GET_CELL(id_cell, m_Grid);
 391       for (int i = 0; i < num_pts; ++i) {
 392         if (m_SurfNode[pts[i]] && m_BoundaryCodes.contains(bc->GetValue(id_cell))) {
 393           n2bc[pts[i]].insert(bc->GetValue(id_cell));
 394         }
 395       }
 396     }
 397   }
 398   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 399     if (n2bc[id_node].size() == 0) {
 400       m_SurfNode[id_node] = false;
 401     }
 402   }
 403   for (int iter = 0; iter < m_NumNormalRelaxations; ++iter) {
 404     QVector<vec3_t> n_new(m_NodeNormal.size(), vec3_t(0,0,0));
 405     for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 406       if (m_SurfNode[id_node] ) {
 407         QList<vtkIdType> snap_points;
 408         for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
 409           vtkIdType id_neigh = m_Part.n2nGG(id_node,i);
 410           bool append = true;
 411           foreach (int bc, n2bc[id_node]) {
 412             if (!n2bc[id_neigh].contains(bc)) {
 413               append = false;
 414               break;
 415             }
 416           }
 417           if (append) {
 418             if (!m_SurfNode[id_neigh]) {
 419               EG_BUG;
 420             }
 421             snap_points.append(id_neigh);
 422           }
 423         }
 424         if (snap_points.size() > 0) {
 425           n_new[id_node] = vec3_t(0,0,0);
 426           foreach (vtkIdType id_snap, snap_points) {
 427             n_new[id_node] += m_NodeNormal[id_snap];
 428           }
 429           n_new[id_node].normalise();
 430         } else {
 431           n_new[id_node] = m_NodeNormal[id_node];
 432         }
 433         if (n_new[id_node].abs() < 0.1) {
 434           EG_BUG;
 435         }
 436       }
 437     }
 438     m_NodeNormal = n_new;
 439     correctNormalVectors();
 440   }
 441 }
 442
 443 void GridSmoother::getRules()
 444 {
 445   QString rules_text = GuiMainWindow::pointer()->getXmlSection("engrid/blayer/rules");
 446   QStringList rules = rules_text.split(";", QString::SkipEmptyParts);
 447   foreach (QString rule, rules) {
 448     rule = rule.trimmed();
 449     QStringList parts = rule.split("=");
 450     if (parts.count() > 1) {
 451       rule_t rule;
 452       QString left = parts[0].trimmed();
 453       rule.h = parts[1].trimmed().toDouble();
 454       QStringList rows = left.split("<OR>");
 455       foreach (QString row, rows) {
 456         QStringList cols = row.split("<AND>");
 457         foreach (QString col, cols) {
 458           rule.bcs.insert(col.toInt());
 459         }
 460       }
 461     }
 462   }
 463 }
 464
 465 void GridSmoother::computeDesiredHeights()
 466 {
 467   // first pass (intial height)
 468   EG_VTKDCN(vtkDoubleArray, cl, m_Grid, "node_meshdensity_desired" );
 469   m_Height.fill(0, m_Grid->GetNumberOfPoints());
 470   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 471     if (m_SurfNode[id_node]) {
 472       m_Height[id_node] = cl->GetValue(id_node);
 473       // if undefined: compute height from surrounding edges
 474       if (m_Height[id_node] < 1e-99) {
 475         m_Height[id_node] = 0;
 476         int N = 0;
 477         vec3_t x;
 478         m_Grid->GetPoint(id_node, x.data());
 479         for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
 480           vtkIdType id_neigh_node = m_Part.n2nGG(id_node, i);
 481           if (m_SurfNode[id_neigh_node]) {
 482             ++N;
 483             vec3_t xn;
 484             m_Grid->GetPoint(id_neigh_node, xn.data());
 485             m_Height[id_node] += (x - xn).abs();
 486           }
 487         }
 488         if (N == 0) {
 489           EG_BUG;
 490         }
 491         m_Height[id_node] /= N;
 492       }
 493     }
 494   }
 495   QVector<double> Dh_max = m_Height;
 496   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 497     Dh_max[id_node] *= m_FarRatio;
 498   }
 499
 500   getRules();
 501   // second pass (correct with absolute height if required)
 502   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 503     if (m_SurfNode[id_node]) {
 504       m_Height[id_node] = m_Blending*m_AbsoluteHeight + (1.0-m_Blending)*m_RelativeHeight*m_Height[id_node];
 505       double h_rule = 1e99;
 506       foreach (rule_t rule, m_Rules) {
 507         bool apply_rule = true;
 508         foreach (int bc, rule.bcs) {
 509           if (!m_Part.hasBC(id_node, bc)) {
 510             apply_rule = false;
 511             break;
 512           }
 513         }
 514         if (apply_rule) {
 515           h_rule = min(h_rule, rule.h);
 516         }
 517       }
 518       if (h_rule < 1e98) {
 519         //m_Height[id_node] = h_rule;
 520       }
 521     }
 522   }
 523
 524
 525   QVector<double> h(m_Grid->GetNumberOfPoints(), 0);
 526   QVector<double> h_opt(m_Grid->GetNumberOfPoints(), 0);
 527   int num_layers = 0;
 528   double err_sq_max = 1e99;
 529   m_NumLayers = 0;
 530   do {
 531     ++num_layers;
 532     double err_sq = 0;
 533     double err_max_iter = 0;
 534     for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 535       if (m_SurfNode[id_node]) {
 536         double Dh = m_Height[id_node];
 537         h[id_node] = m_Height[id_node];
 538         double total_ratio = Dh_max[id_node]/Dh;
 539         double stretch = pow(total_ratio, 1.0/num_layers);
 540         for (int i = 1; i < num_layers; ++i) {
 541           Dh *= stretch;
 542           h[id_node] += Dh;
 543         }
 544         double err;
 545         if (fabs(stretch) > fabs(m_DesiredStretching)) {
 546           err = fabs(stretch - m_DesiredStretching)/m_DesiredStretching;
 547         } else {
 548           err = fabs(stretch - m_DesiredStretching)/stretch;
 549         }
 550         err_max_iter = max(err_max_iter, err);
 551         err_sq += err*err;
 552       }
 553     }
 554     if (err_sq < err_sq_max) {
 555       m_NumLayers = num_layers;
 556       h_opt = h;
 557       err_sq_max = err_sq;
 558     }
 559   } while (num_layers < 200);
 560   m_Height = h_opt;
 561
 562   // correct with angle between face normal and propagation direction (node normals)
 563   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 564     if (m_SurfNode[id_node]) {
 565       int N = 0;
 566       double scale = 0;
 567       for (int j = 0; j < m_Part.n2cGSize(id_node); ++j) {
 568         vtkIdType id_cell = m_Part.n2cGG(id_node, j);
 569         if (isSurface(id_cell, m_Grid)) {
 570           scale += m_NodeNormal[id_node]*cellNormal(m_Grid, id_cell).normalise();
 571         }
 572       }
 573       if (N > 0) {
 574         scale /= N;
 575         m_Height[id_node] /= scale;
 576       }
 577     }
 578   }
 579 }
 580
 581 bool GridSmoother::faceFine(vtkIdType id_face, double scale)
 582 {
 583   EG_GET_CELL(id_face, m_Grid);
 584   if (type_cell != VTK_TRIANGLE) {
 585     EG_BUG;
 586   }
 587   QVector<vec3_t> x1(num_pts);
 588   for (vtkIdType i = 0; i < num_pts; ++i) {
 589     m_Grid->GetPoint(pts[i], x1[i].data());
 590   }
 591   vec3_t n1 = triNormal(x1[0], x1[1], x1[2]);
 592   QVector<vec3_t> x2(num_pts);
 593   for (vtkIdType i = 0; i < num_pts; ++i) {
 594     x2[i] = x1[i] + scale*m_Height[pts[i]]*m_NodeNormal[pts[i]];
 595   }
 596   vec3_t n2 = triNormal(x2[0], x2[1], x2[2]);
 597   double A1 = n1.abs();
 598   double A2 = n2.abs();
 599   if (A2/A1 >= m_FaceSizeLowerLimit && A2/A1 <= m_FaceSizeUpperLimit){
 600     if (angle(n1, n2) < m_FaceAngleLimit) {
 601       return true;
 602     }
 603   }
 604   return false;
 605 }
 606
 607 void GridSmoother::computeHeights()
 608 {
 609   EG_VTKDCC(vtkIntArray, bc, m_Grid, "cell_code");
 610
 611   {
 612     QString blayer_txt = GuiMainWindow::pointer()->getXmlSection("blayer/global");
 613     QTextStream s(&blayer_txt);
 614     if (!s.atEnd()) s >> m_AbsoluteHeight;
 615     if (!s.atEnd()) s >> m_RelativeHeight;
 616     if (!s.atEnd()) s >> m_Blending;
 617     if (!s.atEnd()) s >> m_DesiredStretching;
 618     if (!s.atEnd()) s >> m_FarRatio;
 619     if (!s.atEnd()) s >> m_NumLayers;
 620     if (!s.atEnd()) s >> m_NumHeightRelaxations;
 621     if (!s.atEnd()) s >> m_NumNormalRelaxations;
 622   }
 623   computeDesiredHeights();
 624   {
 625     QString blayer_txt = "";
 626     QTextStream s(&blayer_txt);
 627     s << m_AbsoluteHeight << " ";
 628     s << m_RelativeHeight << " ";
 629     s << m_Blending << " ";
 630     s << m_DesiredStretching << " ";
 631     s << m_FarRatio << " ";
 632     s << m_NumLayers << " ";
 633     s << m_NumHeightRelaxations << " ";
 634     s << m_NumNormalRelaxations << " ";
 635     GuiMainWindow::pointer()->setXmlSection("blayer/global", blayer_txt);
 636   }
 637
 638   // gaps
 639   {
 640
 641     // prepare search list of potential "crash" partner points
 642     QList<vtkIdType> search_nodes;
 643     for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 644       bool append_node = m_SurfNode[id_node];
 645       if (!append_node) {
 646         for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
 647           vtkIdType id_cell = m_Part.n2cGG(id_node, i);
 648           if (isSurface(id_cell, m_Grid)) {
 649             if (!m_LayerAdjacentBoundaryCodes.contains(bc->GetValue(id_cell))) {
 650               append_node = true;
 651               break;
 652             }
 653           }
 654         }
 655       }
 656       if (append_node) {
 657         search_nodes.append(id_node);
 658       }
 659     }
 660
 661     // find close points
 662     QVector<vec3_t> points(search_nodes.size());
 663     for (int i = 0; i < search_nodes.size(); ++i) {
 664       m_Grid->GetPoint(search_nodes[i], points[i].data());
 665     }
 666     PointFinder pfind;
 667     pfind.setPoints(points);
 668
 669     // check for potential collisions
 670     for (vtkIdType id_node1 = 0; id_node1 < m_Grid->GetNumberOfPoints(); ++id_node1) {
 671       if (m_SurfNode[id_node1]) {
 672         vec3_t x1;
 673         m_Grid->GetPoint(id_node1, x1.data());
 674         QVector<int> close_points;
 675         pfind.getClosePoints(x1, close_points, 20*m_Height[id_node1]/m_MaxHeightInGaps);
 676         foreach (int i, close_points) {
 677
 678           // maybe check for topological neighbours and exclude them from the search ...
 679
 680           vtkIdType id_node2 = search_nodes[i];
 681           if (id_node1 != id_node2) {
 682             const vec3_t& n1 = m_NodeNormal[id_node1];
 683             vec3_t x1, x2;
 684             m_Grid->GetPoint(id_node1, x1.data());
 685             m_Grid->GetPoint(id_node2, x2.data());
 686             vec3_t Dx = x2 - x1;
 687             double a = Dx*n1;
 688             if (a > 0) {
 689               double b = Dx.abs();
 690               double alpha = 180.0/M_PI*acos(a/b); /// @todo This is very slow; look at alternatives!
 691               if (alpha < m_RadarAngle) {
 692                 m_Height[id_node1] = min(m_Height[id_node1], m_MaxHeightInGaps*a);
 693               }
 694             }
 695           }
 696         }
 697       }
 698     }
 699   }
 700
 701   // limit face size difference (neighbour collisions)
 702   {
 703     bool done;
 704     do {
 705       done = true;
 706       for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 707         if (isSurface(id_cell, m_Grid)) {
 708           EG_GET_CELL(id_cell, m_Grid);
 709           bool check_face = true;
 710           for (vtkIdType i = 0; i < num_pts; ++i) {
 711             if (!m_SurfNode[pts[i]]) {
 712               check_face = false;
 713               break;
 714             }
 715           }
 716           if (check_face) {
 717             if (!faceFine(id_cell, 1)) {
 718               done = false;
 719               double scale1 = 0;
 720               double scale2 = 1;
 721               while (scale2 - scale1 > 1e-3) {
 722                 if (faceFine(id_cell, 0.5*(scale1 + scale2))) {
 723                   scale1 = 0.5*(scale1 + scale2);
 724                 } else {
 725                   scale2 = 0.5*(scale1 + scale2);
 726                 }
 727               }
 728               for (vtkIdType i = 0; i < num_pts; ++i) {
 729                 m_Height[pts[i]] *= 0.5*(scale1 + scale2);
 730               }
 731             }
 732           }
 733         }
 734       }
 735     } while (!done);
 736   }
 737
 738   // smoothing
 739   QVector<int> num_bcs(m_Grid->GetNumberOfPoints(),0);
 740   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 741     if (m_SurfNode[id_node]) {
 742       for (int i = 0; i < m_Part.n2bcGSize(id_node); ++i) {
 743         int bc = m_Part.n2bcG(id_node, i);
 744         if (m_BoundaryCodes.contains(bc)) {
 745           ++num_bcs[id_node];
 746         }
 747       }
 748     }
 749   }
 750   for (int iter = 0; iter < m_NumHeightRelaxations; ++iter) {
 751     QVector<double> h_new = m_Height;
 752     for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 753       if (m_SurfNode[id_node]) {
 754         QList<vtkIdType> snap_points;
 755         for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
 756           vtkIdType id_neigh = m_Part.n2nGG(id_node,i);
 757           if (num_bcs[id_node] <= num_bcs[id_neigh]) {
 758             if (!m_SurfNode[id_neigh]) {
 759               EG_BUG;
 760             }
 761             snap_points.append(id_neigh);
 762           }
 763         }
 764         if (snap_points.size() > 0) {
 765           h_new[id_node] = 0;
 766           if (snap_points.size() == 2) {
 767             snap_points.append(id_node);
 768           }
 769           foreach (vtkIdType id_snap, snap_points) {
 770             h_new[id_node] += m_Height[id_snap];
 771           }
 772           h_new[id_node] /= snap_points.size();
 773           //h_new[id_node] = min(h_new[id_node], m_Height[id_node]);
 774         } else {
 775           h_new[id_node] = m_Height[id_node];
 776         }
 777         //h_new[id_node] = min(h_new[id_node], height_limit[id_node]);
 778       }
 779     }
 780     m_Height = h_new;
 781   }
 782
 783
 784 }
 785
 786 void GridSmoother::correctNormalVectors()
 787 {
 788   EG_VTKDCC(vtkIntArray, bc, m_Grid, "cell_code");
 789   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 790     if (m_NodeNormal[id_node].abs() > 0.1) {
 791       for (int iter = 0; iter < m_NumBoundaryCorrections; ++iter) {
 792         for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
 793           vtkIdType id_cell = m_Part.n2cGG(id_node,i);
 794           if (isSurface(id_cell, m_Grid)) {
 795             if (!m_BoundaryCodes.contains(bc->GetValue(id_cell))) {
 796               vec3_t v = m_NodeNormal[id_node];
 797               vec3_t n = GeometryTools::cellNormal(m_Grid, id_cell);
 798               n.normalise();
 799               v -= (n*m_NodeNormal[id_node])*n;
 800               v.normalise();
 801               m_NodeNormal[id_node] = v;
 802             }
 803           }
 804         }
 805       }
 806     }
 807   }
 808 }
 809
 810 void GridSmoother::computeFeet()
 811 {
 812   m_IdFoot.fill(-1, m_Grid->GetNumberOfPoints());
 813   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 814     if (m_Grid->GetCellType(id_cell) == VTK_WEDGE) {
 815       EG_GET_CELL(id_cell, m_Grid);
 816       m_IdFoot[pts[3]] = pts[0];
 817       m_IdFoot[pts[4]] = pts[1];
 818       m_IdFoot[pts[5]] = pts[2];
 819       m_NodeMarked[pts[0]] = false;
 820       m_NodeMarked[pts[1]] = false;
 821       m_NodeMarked[pts[2]] = false;
 822     }
 823   }
 824 }
 825
 826 void GridSmoother::simpleNodeMovement(int i_nodes)
 827 {
 828   vtkIdType id_node = m_Part.globalNode(i_nodes);
 829   vtkIdType id_foot = m_IdFoot[id_node];
 830   vec3_t x_new(0,0,0), x_node;
 831   m_Grid->GetPoint(id_node, x_node.data());
 832
 833   vec3_t x_surf(0,0,0);
 834   if (m_SurfNode[id_node]) {
 835     int N = 0;
 836     for (int i = 0; i < m_Part.n2nLSize(i_nodes); ++i) {
 837       vtkIdType id_neigh = m_Part.n2nLG(i_nodes,i);
 838       if (m_SurfNode[id_neigh]) {
 839         vec3_t x_neigh;
 840         m_Grid->GetPoint(id_neigh, x_neigh.data());
 841         x_surf += x_neigh;
 842         ++N;
 843       }
 844     }
 845     if (N == 0) {
 846       EG_BUG;
 847     } else {
 848       x_surf *= 1.0/N;
 849     }
 850   }
 851
 852   if (id_foot != -1) {
 853     vec3_t x_foot, x_node;
 854     m_Grid->GetPoint(id_foot, x_foot.data());
 855     //double H = m_Blending*m_AbsoluteHeight + (1.0-m_Blending)*m_RelativeHeight*m_Height[id_foot];
 856     double H = m_Height[id_foot];
 857     x_new = x_foot + H*m_NodeNormal[id_foot];
 858   } else {
 859     if (m_SurfNode[id_node]) {
 860       x_new = x_surf;
 861     } else {
 862       if (m_Part.n2nLSize(i_nodes) == 0) {
 863         EG_BUG;
 864       }
 865       for (int i = 0; i < m_Part.n2nLSize(i_nodes); ++i) {
 866         vtkIdType id_neigh = m_Part.n2nLG(i_nodes,i);
 867         vec3_t x_neigh;
 868         m_Grid->GetPoint(id_neigh, x_neigh.data());
 869         x_new += x_neigh;
 870       }
 871       x_new *= 1.0/m_Part.n2nLSize(i_nodes);
 872     }
 873   }
 874   vec3_t Dx = x_new - x_node;
 875   correctDx(i_nodes, Dx);
 876   moveNode(i_nodes, Dx);
 877 }
 878
 879 void GridSmoother::operate()
 880 {
 881   if (m_FirstCall) {
 882     markNodes();
 883     computeNormals();
 884     computeFeet();
 885     computeHeights();
 886     m_FirstCall = false;
 887   }
 888   l2g_t nodes = m_Part.getNodes();
 889   for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
 890     if (m_NodeMarked[nodes[i_nodes]]) {
 891       simpleNodeMovement(i_nodes);
 892     }
 893   }
 894 }
 895
 896 void GridSmoother::writeDebugFile(QString file_name)
 897 {
 898   QVector<vtkIdType> bcells;
 899   QSet<int> bcs = getAllBoundaryCodes(m_Grid);
 900   getSurfaceCells(bcs, bcells, m_Grid);
 901   MeshPartition bpart(m_Grid);
 902   bpart.setCells(bcells);
 903   file_name = GuiMainWindow::pointer()->getCwd() + "/" + file_name + ".vtk";
 904   QFile file(file_name);
 905   file.open(QIODevice::WriteOnly);
 906   QTextStream f(&file);
 907   f << "# vtk DataFile Version 2.0\n";
 908   f << "m_NodeNormal\n";
 909   f << "ASCII\n";
 910   f << "DATASET UNSTRUCTURED_GRID\n";
 911   f << "POINTS " << bpart.getNumberOfNodes() << " float\n";
 912   for (int i = 0; i < bpart.getNumberOfNodes(); ++i) {
 913     vec3_t x;
 914     vtkIdType id_node = bpart.globalNode(i);
 915     m_Grid->GetPoint(id_node, x.data());
 916     f << x[0] << " " << x[1] << " " << x[2] << "\n";
 917   }
 918   f << "CELLS " << bpart.getNumberOfCells() << " ";
 919   int N = 0;
 920   for (int i = 0; i < bpart.getNumberOfCells(); ++ i) {
 921     vtkIdType id_cell = bpart.globalCell(i);
 922     EG_GET_CELL(id_cell, m_Grid);
 923     N += 1 + num_pts;
 924   }
 925   f << N << "\n";
 926   for (int i = 0; i < bpart.getNumberOfCells(); ++ i) {
 927     vtkIdType id_cell = bpart.globalCell(i);
 928     EG_GET_CELL(id_cell, m_Grid);
 929     f << num_pts;
 930     for (int j = 0; j < num_pts; ++j) {
 931       f << " " << bpart.localNode(pts[j]);
 932     }
 933     f << "\n";
 934   }
 935
 936   f << "CELL_TYPES " << bpart.getNumberOfCells() << "\n";
 937   for (int i = 0; i < bpart.getNumberOfCells(); ++ i) {
 938     vtkIdType id_cell = bpart.globalCell(i);
 939     vtkIdType type_cell = m_Grid->GetCellType(id_cell);
 940     f << type_cell << "\n";
 941   }
 942   f << "POINT_DATA " << bpart.getNumberOfNodes() << "\n";
 943
 944   f << "VECTORS N float\n";
 945   for (int i = 0; i < bpart.getNumberOfNodes(); ++i) {
 946     vtkIdType id_node = bpart.globalNode(i);
 947     f << m_NodeNormal[id_node][0] << " " << m_NodeNormal[id_node][1] << " " << m_NodeNormal[id_node][2] << "\n";
 948   }
 949
 950   f << "SCALARS h float\n";
 951   f << "LOOKUP_TABLE default\n";
 952   for (int i = 0; i < bpart.getNumberOfNodes(); ++i) {
 953     vtkIdType id_node = bpart.globalNode(i);
 954     f << m_Height[id_node] << "\n";
 955   }
 956 }
 957
 958