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