src/libengrid/gridsmoother.cpp

   1 //
   2 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   3 // +                                                                      +
   4 // + This file is part of enGrid.                                         +
   5 // +                                                                      +
   6 // + Copyright 2008-2012 enGits GmbH                                     +
   7 // +                                                                      +
   8 // + enGrid is free software: you can redistribute it and/or modify       +
   9 // + it under the terms of the GNU General Public License as published by +
  10 // + the Free Software Foundation, either version 3 of the License, or    +
  11 // + (at your option) any later version.                                  +
  12 // +                                                                      +
  13 // + enGrid is distributed in the hope that it will be useful,            +
  14 // + but WITHOUT ANY WARRANTY; without even the implied warranty of       +
  15 // + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the        +
  16 // + GNU General Public License for more details.                         +
  17 // +                                                                      +
  18 // + You should have received a copy of the GNU General Public License    +
  19 // + along with enGrid. If not, see <http://www.gnu.org/licenses/>.       +
  20 // +                                                                      +
  21 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  22 //
  23 #include "gridsmoother.h"
  24 #include "guimainwindow.h"
  25 #include "elements.h"
  26 #include "optimisenormalvector.h"
  27 #include "pointfinder.h"
  28
  29 #include <QTime>
  30
  31 GridSmoother::GridSmoother()
  32 {
  33   m_NumIterations          = 5;
  34   m_NumRelaxations         = 5;
  35   m_NumBoundaryCorrections = 50;
  36   m_DesiredStretching      = 1.2;
  37   m_FirstCall = true;
  38
  39   getSet("boundary layer", "number of smoothing sub-iterations",       5,     m_NumIterations);
  40   getSet("boundary layer", "use strict prism checking",                false, m_StrictPrismChecking);
  41   getSet("boundary layer", "number of normal vector relax iterations", 10,    m_NumNormalRelaxations);
  42   getSet("boundary layer", "number of layer height relax iterations",  3,     m_NumHeightRelaxations);
  43   getSet("boundary layer", "radar angle",                              45,    m_RadarAngle);
  44   getSet("boundary layer", "maximal layer height in gaps",             0.2,   m_MaxHeightInGaps);
  45
  46   //m_CritAngle = GeometryTools::deg2rad(m_CritAngle);
  47 }
  48
  49 void GridSmoother::markNodes()
  50 {
  51   m_NodeMarked.fill(false,m_Grid->GetNumberOfPoints());
  52   QVector<bool> new_mark(m_Grid->GetNumberOfPoints());
  53   for (int i_iterations = 0; i_iterations < 1; ++i_iterations) {
  54     qCopy(m_NodeMarked.begin(),m_NodeMarked.end(),new_mark.begin());
  55     for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
  56       bool mark_cell = false;
  57       vtkIdType type_cell, N_pts, *pts;
  58       type_cell = m_Grid->GetCellType(id_cell);
  59       m_Grid->GetCellPoints(id_cell, N_pts, pts);
  60       if (type_cell == VTK_WEDGE) {
  61         mark_cell = true;
  62       } else {
  63         for (int i_pts = 0; i_pts < N_pts; ++i_pts) {
  64           if (m_NodeMarked[pts[i_pts]]) {
  65             mark_cell = true;
  66           }
  67         }
  68       }
  69       if (mark_cell) {
  70         for (int i_pts = 0; i_pts < N_pts; ++i_pts) {
  71           new_mark[pts[i_pts]] = true;
  72         }
  73       }
  74     }
  75     qCopy(new_mark.begin(),new_mark.end(),m_NodeMarked.begin());
  76   }
  77   QSet<int> free_bcs = m_BoundaryCodes + m_LayerAdjacentBoundaryCodes;
  78   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
  79   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
  80     if (isSurface(id_cell, m_Grid)) {
  81       if (!free_bcs.contains(cell_code->GetValue(id_cell))) {
  82         vtkIdType N_pts, *pts;
  83         m_Grid->GetCellPoints(id_cell, N_pts, pts);
  84         for (int i_pts = 0; i_pts < N_pts; ++i_pts) {
  85           m_NodeMarked[pts[i_pts]] = false;
  86         }
  87       }
  88     }
  89   }
  90   m_NumMarkedNodes = 0;
  91   QVector<vtkIdType> nodes = m_Part.getNodes();
  92   foreach (vtkIdType id_node, nodes) {
  93     if (id_node < 0) EG_BUG;
  94     if (id_node > m_Grid->GetNumberOfPoints()) EG_BUG;
  95     if (m_NodeMarked[id_node]) {
  96       ++m_NumMarkedNodes;
  97     }
  98   }
  99 }
 100
 101 bool GridSmoother::setNewPosition(vtkIdType id_node, vec3_t x_new)
 102 {
 103   using namespace GeometryTools;
 104
 105   vec3_t x_old;
 106   m_Grid->GetPoint(id_node, x_old.data());
 107   m_Grid->GetPoints()->SetPoint(id_node, x_new.data());
 108   bool move = true;
 109
 110   if (move) {
 111     Elements E;
 112
 113     l2g_t cells = m_Part.getCells();
 114     l2l_t n2c = m_Part.getN2C();
 115
 116     foreach (int i_cells, n2c[id_node]) {
 117       vtkIdType id_cell = cells[i_cells];
 118       vtkIdType type_cell = m_Grid->GetCellType(id_cell);
 119
 120       if (type_cell == VTK_TRIANGLE) {
 121         vtkIdType N_pts, *pts;
 122         vec3_t x[3];
 123         m_Grid->GetCellPoints(id_cell, N_pts, pts);
 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         vtkIdType N_pts, *pts;
 149         vec3_t x[4];
 150         m_Grid->GetCellPoints(id_cell, N_pts, pts);
 151         for (int i = 0; i < 4; ++i) {
 152           m_Grid->GetPoint(pts[i], x[i].data());
 153         }
 154         double L_max = 0;
 155         for (int i = 0; i < 4; ++i) {
 156           for (int j = 0; j < 4; ++j) {
 157             if (i != j) {
 158               L_max = max(L_max, (x[i]-x[j]).abs());
 159             }
 160           }
 161         }
 162         if (GeometryTools::cellVA(m_Grid, id_cell) < 1e-3*L_max*L_max*L_max) {
 163           move = false;
 164         }
 165       }
 166
 167       if (type_cell == VTK_WEDGE && m_StrictPrismChecking) {
 168         vtkIdType N_pts, *pts;
 169         vec3_t xtet[4];
 170         m_Grid->GetCellPoints(id_cell, N_pts, pts);
 171         bool ok = true;
 172         for (int i = 0; i < 4; ++i) {     // variation
 173           ok = true;
 174           for (int j = 0; j < 3; ++j) {   // tetrahedron
 175             for (int k = 0; k < 4; ++k) { // node
 176               m_Grid->GetPoint(pts[E.priTet(i,j,k)], xtet[k].data());
 177             }
 178             double V = GeometryTools::tetraVol(xtet[0], xtet[1], xtet[2], xtet[3]);
 179             if (V <= 0) {
 180               ok = false;
 181             }
 182           }
 183           if (ok) {
 184             break;
 185           }
 186         }
 187         if (!ok) {
 188           move = false;
 189         }
 190       }
 191     }
 192   }
 193   if (!move) {
 194     m_Grid->GetPoints()->SetPoint(id_node, x_old.data());
 195   }
 196   return move;
 197 }
 198
 199 void GridSmoother::correctDx(int i_nodes, vec3_t &Dx)
 200 {
 201   l2g_t nodes = m_Part.getNodes();
 202   l2l_t n2c = m_Part.getN2C();
 203   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
 204   vec3_t x_old;
 205   m_Grid->GetPoint(nodes[i_nodes],x_old.data());
 206   for (int i_boundary_correction = 0; i_boundary_correction < m_NumBoundaryCorrections; ++i_boundary_correction) {
 207     foreach (vtkIdType id_cell, n2c[i_nodes]) {
 208       if (isSurface(id_cell, m_Grid)) {
 209         int bc = cell_code->GetValue(id_cell);
 210         vec3_t x_new = x_old + Dx;
 211         x_new = GuiMainWindow::pointer()->getSurfProj(bc)->projectRestricted(x_new, nodes[i_nodes]);
 212         Dx = x_new - x_old;
 213       } else {
 214         if (m_Grid->GetCellType(id_cell) == VTK_WEDGE) {
 215           vtkIdType N_pts, *pts;
 216           m_Grid->GetCellPoints(id_cell, N_pts, pts);
 217           vtkIdType id_surf_node = -1;
 218           if (pts[3] == nodes[i_nodes]) id_surf_node = pts[0];
 219           if (pts[4] == nodes[i_nodes]) id_surf_node = pts[1];
 220           if (pts[5] == nodes[i_nodes]) id_surf_node = pts[2];
 221           if (id_surf_node != -1) {
 222             vec3_t x0,x1,x2;
 223             m_Grid->GetPoint(pts[0],x0.data());
 224             m_Grid->GetPoint(pts[1],x1.data());
 225             m_Grid->GetPoint(pts[2],x2.data());
 226             vec3_t a = x1-x0;
 227             vec3_t b = x2-x0;
 228             vec3_t c = b-a;
 229             double L = (a.abs()+b.abs()+c.abs())/3.0;
 230             vec3_t n = b.cross(a);
 231             n.normalise();
 232             vec3_t x_old;
 233             m_Grid->GetPoint(nodes[i_nodes],x_old.data());
 234             vec3_t x_new = x_old + Dx - x0;
 235             if ( (n*x_new) <= 0 ) {
 236               x_new -= (x_new*n)*n;
 237               x_new += 1e-4*L*n;
 238               x_new += x0;
 239               Dx = x_new - x_old;
 240             }
 241           }
 242         }
 243       }
 244     }
 245   }
 246 }
 247
 248 bool GridSmoother::moveNode(int i_nodes, vec3_t &Dx)
 249 {
 250   m_CollisionDetected = false;
 251   l2g_t nodes = m_Part.getNodes();
 252   vtkIdType id_node = nodes[i_nodes];
 253   if (id_node == 27029) {
 254     cout << "break-point" << endl;
 255   }
 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   }
 266   return moved;
 267 }
 268
 269 void GridSmoother::computeNormals()
 270 {
 271   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
 272   m_NodeNormal.fill(vec3_t(0,0,0), m_Grid->GetNumberOfPoints());
 273   QVector<int> num_bcs(m_Grid->GetNumberOfPoints());
 274   QVector<OptimiseNormalVector> n_opt(m_Grid->GetNumberOfPoints());
 275   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 276     QSet<int> bcs;
 277     for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
 278       vtkIdType id_cell = m_Part.n2cGG(id_node, i);
 279       if (isSurface(id_cell, m_Grid)) {
 280         int bc = cell_code->GetValue(id_cell);
 281         if (m_BoundaryCodes.contains(bc)) {
 282           bcs.insert(bc);
 283         }
 284       }
 285     }
 286     num_bcs[id_node] = bcs.size();
 287     QVector<vec3_t> normal(num_bcs[id_node], vec3_t(0,0,0));
 288     QMap<int,int> bcmap;
 289     int i_bc = 0;
 290     foreach (int bc, bcs) {
 291       bcmap[bc] = i_bc;
 292       ++i_bc;
 293     }
 294     for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
 295       vtkIdType id_cell = m_Part.n2cGG(id_node, i);
 296       if (isSurface(id_cell, m_Grid)) {
 297         int bc = cell_code->GetValue(id_cell);
 298         vtkIdType N_pts, *pts;
 299         m_Grid->GetCellPoints(id_cell, N_pts, pts);
 300         vec3_t a, b, c;
 301         for (int j = 0; j < N_pts; ++j) {
 302           if (pts[j] == id_node) {
 303             m_Grid->GetPoint(pts[j], a.data());
 304             if (j > 0) {
 305               m_Grid->GetPoint(pts[j-1], b.data());
 306             } else {
 307               m_Grid->GetPoint(pts[N_pts-1], b.data());
 308             }
 309             if (j < N_pts - 1) {
 310               m_Grid->GetPoint(pts[j+1], c.data());
 311             } else {
 312               m_Grid->GetPoint(pts[0], c.data());
 313             }
 314           }
 315         }
 316         vec3_t u = b - a;
 317         vec3_t v = c - a;
 318         double alpha = GeometryTools::angle(u, v);
 319         vec3_t n = u.cross(v);
 320         n.normalise();
 321         if (m_BoundaryCodes.contains(bc)) {
 322           normal[bcmap[bc]] += alpha*n;
 323           n_opt[id_node].addFace(n);
 324         } else {
 325           n_opt[id_node].addConstraint(n);
 326         }
 327       }
 328     }
 329     for (int i = 0; i < num_bcs[id_node]; ++i) {
 330       normal[i].normalise();
 331     }
 332     if (num_bcs[id_node] > 0) {
 333       if (num_bcs[id_node] > 1) {
 334         if (num_bcs[id_node] == 3) {
 335           for (int i = 0; i < num_bcs[id_node]; ++i) {
 336             for (int j = i + 1; j < num_bcs[id_node]; ++j) {
 337               vec3_t n = normal[i] + normal[j];
 338               n.normalise();
 339               m_NodeNormal[id_node] += n;
 340             }
 341           }
 342         } else {
 343           for (int i = 0; i < num_bcs[id_node]; ++i) {
 344             m_NodeNormal[id_node] += normal[i];
 345           }
 346         }
 347       } else {
 348         m_NodeNormal[id_node] = normal[0];
 349       }
 350       m_NodeNormal[id_node].normalise();
 351       m_NodeNormal[id_node] = n_opt[id_node](m_NodeNormal[id_node]);
 352       m_NodeNormal[id_node].normalise();
 353     }
 354   }
 355
 356   relaxNormalVectors();
 357
 358   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 359     if (m_NodeNormal[id_node].abs() < 0.1) {
 360       vec3_t n(0,0,0);
 361       int N = 0;
 362       for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
 363         vtkIdType id_cell = m_Part.n2cGG(id_node, i);
 364         if (isSurface(id_cell, m_Grid)) {
 365           n += GeometryTools::cellNormal(m_Grid, id_cell);
 366           ++N;
 367         }
 368       }
 369       if (N) {
 370         n.normalise();
 371         m_NodeNormal[id_node] = n;
 372       }
 373     }
 374     if (num_bcs[id_node] > 1) {
 375       m_NodeNormal[id_node] = n_opt[id_node](m_NodeNormal[id_node]);
 376     }
 377   }
 378
 379 }
 380
 381 void GridSmoother::relaxNormalVectors()
 382 {
 383   m_SurfNode.fill(false, m_Grid->GetNumberOfPoints());
 384   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 385     for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
 386       if (isSurface(m_Part.n2cGG(id_node,i), m_Grid)) {
 387         m_SurfNode[id_node] = true;
 388         break;
 389       }
 390     }
 391   }
 392   EG_VTKDCC(vtkIntArray, bc, m_Grid, "cell_code");
 393   QVector<QSet<int> > n2bc(m_Grid->GetNumberOfPoints());
 394   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 395     if (isSurface(id_cell, m_Grid)) {
 396       vtkIdType N_pts, *pts;
 397       m_Grid->GetCellPoints(id_cell, N_pts, pts);
 398       for (int i = 0; i < N_pts; ++i) {
 399         if (m_SurfNode[pts[i]]) {
 400           n2bc[pts[i]].insert(bc->GetValue(id_cell));
 401         }
 402       }
 403     }
 404   }
 405   QVector<int> num_bcs(m_Grid->GetNumberOfPoints(),0);
 406   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 407     if (m_SurfNode[id_node]) {
 408       QList<vtkIdType> snap_points;
 409       foreach (int bc, n2bc[id_node]) {
 410         if (m_BoundaryCodes.contains(bc)) {
 411           ++num_bcs[id_node];
 412         }
 413       }
 414     }
 415   }
 416   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 417     if (num_bcs[id_node] == 0) {
 418       m_SurfNode[id_node] = false;
 419       n2bc[id_node].clear();
 420     }
 421   }
 422   for (int iter = 0; iter < m_NumNormalRelaxations; ++iter) {
 423     QVector<vec3_t> n_new(m_NodeNormal.size(), vec3_t(0,0,0));
 424     for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 425       if (m_SurfNode[id_node] ) {
 426         QList<vtkIdType> snap_points;
 427         for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
 428           vtkIdType id_neigh = m_Part.n2nGG(id_node,i);
 429           if (num_bcs[id_node] <= num_bcs[id_neigh]) {
 430             if (!m_SurfNode[id_neigh]) {
 431               cout << id_node << ',' << m_Part.n2nGG(id_node,i) << ',' << num_bcs[id_node] << ',' << num_bcs[id_neigh] << endl;
 432               EG_BUG;
 433             }
 434             snap_points.append(id_neigh);
 435           }
 436         }
 437         if (snap_points.size() > 0) {
 438           n_new[id_node] = vec3_t(0,0,0);
 439           foreach (vtkIdType id_snap, snap_points) {
 440             n_new[id_node] += m_NodeNormal[id_snap];
 441           }
 442           n_new[id_node].normalise();
 443         } else {
 444           n_new[id_node] = m_NodeNormal[id_node];
 445         }
 446         if (n_new[id_node].abs() < 0.1) {
 447           cout << id_node << ',' << n_new[id_node] << ',' << num_bcs[id_node] << ',' << n2bc[id_node] << endl;
 448           EG_BUG;
 449         }
 450       }
 451     }
 452     m_NodeNormal = n_new;
 453     correctNormalVectors();
 454     /*
 455     QString num;
 456     num.setNum(iter);
 457     num = "normals_" + num;
 458     writeDebugFile(num);
 459     */
 460   }
 461 }
 462
 463 void GridSmoother::getRules()
 464 {
 465   QString rules_text = GuiMainWindow::pointer()->getXmlSection("engrid/blayer/rules");
 466   QStringList rules = rules_text.split(";", QString::SkipEmptyParts);
 467   foreach (QString rule, rules) {
 468     rule = rule.trimmed();
 469     QStringList parts = rule.split("=");
 470     if (parts.count() > 1) {
 471       rule_t rule;
 472       QString left = parts[0].trimmed();
 473       rule.h = parts[1].trimmed().toDouble();
 474       QStringList rows = left.split("<OR>");
 475       foreach (QString row, rows) {
 476         QStringList cols = row.split("<AND>");
 477         foreach (QString col, cols) {
 478           rule.bcs.insert(col.toInt());
 479         }
 480       }
 481     }
 482   }
 483 }
 484
 485 void GridSmoother::computeDesiredHeights()
 486 {
 487   // first pass (intial height)
 488   EG_VTKDCN(vtkDoubleArray, cl, m_Grid, "node_meshdensity_desired" );
 489   m_Height.fill(0, m_Grid->GetNumberOfPoints());
 490   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 491     if (m_SurfNode[id_node]) {
 492       m_Height[id_node] = cl->GetValue(id_node);
 493       // if undefined: compute height from surrounding edges
 494       if (m_Height[id_node] < 1e-99) {
 495         m_Height[id_node] = 0;
 496         int N = 0;
 497         vec3_t x;
 498         m_Grid->GetPoint(id_node, x.data());
 499         for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
 500           vtkIdType id_neigh_node = m_Part.n2nGG(id_node, i);
 501           if (m_SurfNode[id_neigh_node]) {
 502             ++N;
 503             vec3_t xn;
 504             m_Grid->GetPoint(id_neigh_node, xn.data());
 505             m_Height[id_node] += (x - xn).abs();
 506           }
 507         }
 508         if (N == 0) {
 509           EG_BUG;
 510         }
 511         m_Height[id_node] /= N;
 512       }
 513     }
 514   }
 515   QVector<double> Dh_max = m_Height;
 516   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 517     Dh_max[id_node] *= m_FarRatio;
 518   }
 519
 520   getRules();
 521   // second pass (correct with absolute height if required)
 522   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 523     if (m_SurfNode[id_node]) {
 524       m_Height[id_node] = m_Blending*m_AbsoluteHeight + (1.0-m_Blending)*m_RelativeHeight*m_Height[id_node];
 525       double h_rule = 1e99;
 526       foreach (rule_t rule, m_Rules) {
 527         bool apply_rule = true;
 528         foreach (int bc, rule.bcs) {
 529           if (!m_Part.hasBC(id_node, bc)) {
 530             apply_rule = false;
 531             break;
 532           }
 533         }
 534         if (apply_rule) {
 535           h_rule = min(h_rule, rule.h);
 536         }
 537       }
 538       if (h_rule < 1e98) {
 539         //m_Height[id_node] = h_rule;
 540       }
 541     }
 542   }
 543
 544
 545   QVector<double> h(m_Grid->GetNumberOfPoints(), 0);
 546   QVector<double> h_opt(m_Grid->GetNumberOfPoints(), 0);
 547   int num_layers = 0;
 548   double err_sq_max = 1e99;
 549   m_NumLayers = 0;
 550   do {
 551     ++num_layers;
 552     double err_sq = 0;
 553     double err_max_iter = 0;
 554     for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 555       if (m_SurfNode[id_node]) {
 556         double Dh = m_Height[id_node];
 557         h[id_node] = m_Height[id_node];
 558         double total_ratio = Dh_max[id_node]/Dh;
 559         double stretch = pow(total_ratio, 1.0/num_layers);
 560         for (int i = 1; i < num_layers; ++i) {
 561           Dh *= stretch;
 562           h[id_node] += Dh;
 563         }
 564         double err;
 565         if (fabs(stretch) > fabs(m_DesiredStretching)) {
 566           err = fabs(stretch - m_DesiredStretching)/m_DesiredStretching;
 567         } else {
 568           err = fabs(stretch - m_DesiredStretching)/stretch;
 569         }
 570         err_max_iter = max(err_max_iter, err);
 571         err_sq += err*err;
 572       }
 573     }
 574     if (err_sq < err_sq_max) {
 575       m_NumLayers = num_layers;
 576       h_opt = h;
 577       err_sq_max = err_sq;
 578     }
 579   } while (num_layers < 200);
 580   m_Height = h_opt;
 581 }
 582
 583 void GridSmoother::computeHeights()
 584 {
 585   EG_VTKDCC(vtkIntArray, bc, m_Grid, "cell_code");
 586
 587   {
 588     QString blayer_txt = GuiMainWindow::pointer()->getXmlSection("blayer/global");
 589     QTextStream s(&blayer_txt);
 590     if (!s.atEnd()) s >> m_AbsoluteHeight;
 591     if (!s.atEnd()) s >> m_RelativeHeight;
 592     if (!s.atEnd()) s >> m_Blending;
 593     if (!s.atEnd()) s >> m_DesiredStretching;
 594     if (!s.atEnd()) s >> m_NumLayers;
 595     if (!s.atEnd()) s >> m_FarRatio;
 596   }
 597   computeDesiredHeights();
 598   {
 599     QString blayer_txt = "";
 600     QTextStream s(&blayer_txt);
 601     s << m_AbsoluteHeight << " ";
 602     s << m_RelativeHeight << " ";
 603     s << m_Blending << " ";
 604     s << m_DesiredStretching << " ";
 605     s << m_NumLayers << " ";
 606     s << m_FarRatio << " ";
 607     GuiMainWindow::pointer()->setXmlSection("blayer/global", blayer_txt);
 608   }
 609
 610   // third pass (gaps)
 611   QList<vtkIdType> search_nodes;
 612   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 613     bool append_node = m_SurfNode[id_node];
 614     if (!append_node) {
 615       for (int i = 0; i < m_Part.n2bcGSize(id_node); ++i) {
 616         vtkIdType id_cell = m_Part.n2cGG(id_node, i);
 617         if (isSurface(id_cell, m_Grid)) {
 618           if (!m_LayerAdjacentBoundaryCodes.contains(bc->GetValue(id_cell))) {
 619             append_node = true;
 620             break;
 621           }
 622         }
 623       }
 624     }
 625     if (append_node) {
 626       search_nodes.append(id_node);
 627     }
 628   }
 629
 630   QVector<vec3_t> points(search_nodes.size());
 631   for (int i = 0; i < search_nodes.size(); ++i) {
 632     m_Grid->GetPoint(search_nodes[i], points[i].data());
 633   }
 634   PointFinder pfind;
 635   pfind.setPoints(points);
 636
 637   for (vtkIdType id_node1 = 0; id_node1 < m_Grid->GetNumberOfPoints(); ++id_node1) {
 638     if (m_SurfNode[id_node1]) {
 639       const vec3_t& n1 = m_NodeNormal[id_node1];
 640       vec3_t x1;
 641       m_Grid->GetPoint(id_node1, x1.data());
 642       QVector<int> close_points;
 643       pfind.getClosePoints(x1, close_points, 2*m_Height[id_node1]/m_MaxHeightInGaps);
 644       foreach (int i, close_points) {
 645         vtkIdType id_node2 = search_nodes[i];
 646         if (id_node1 != id_node2) {
 647           vec3_t x2;
 648           m_Grid->GetPoint(id_node2, x2.data());
 649           vec3_t Dx = x2 - x1;
 650           double a = Dx*n1;
 651           if (a > 0) {
 652             double b = Dx.abs();
 653             double alpha = 180.0/M_PI*acos(a/b);
 654             if (alpha < m_RadarAngle) {
 655               m_Height[id_node1] = min(m_Height[id_node1], m_MaxHeightInGaps*a);
 656             }
 657           }
 658         }
 659       }
 660     }
 661   }
 662
 663   // fourth pass (smoothing)
 664   for (int iter = 0; iter < m_NumHeightRelaxations; ++iter) {
 665     QVector<double> h_new = m_Height;
 666     for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 667       if (m_SurfNode[id_node]) {
 668         int N = 0;
 669         h_new[id_node] = 0;
 670         for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
 671           vtkIdType id_neigh = m_Part.n2nGG(id_node,i);
 672           if (m_SurfNode[id_neigh]) {
 673             h_new[id_node] += m_Height[id_neigh];
 674             ++N;
 675           }
 676         }
 677         if (N == 0) {
 678           EG_BUG;
 679         }
 680         h_new[id_node] /= N;
 681       }
 682     }
 683     m_Height = h_new;
 684   }
 685 }
 686
 687 void GridSmoother::correctNormalVectors()
 688 {
 689   EG_VTKDCC(vtkIntArray, bc, m_Grid, "cell_code");
 690   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 691     if (m_NodeNormal[id_node].abs() > 0.1) {
 692       for (int iter = 0; iter < m_NumBoundaryCorrections; ++iter) {
 693         for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
 694           vtkIdType id_cell = m_Part.n2cGG(id_node,i);
 695           if (isSurface(id_cell, m_Grid)) {
 696             if (!m_BoundaryCodes.contains(bc->GetValue(id_cell))) {
 697               vec3_t v = m_NodeNormal[id_node];
 698               vec3_t n = GeometryTools::cellNormal(m_Grid, id_cell);
 699               n.normalise();
 700               v -= (n*m_NodeNormal[id_node])*n;
 701               v.normalise();
 702               m_NodeNormal[id_node] = v;
 703             }
 704           }
 705         }
 706       }
 707     }
 708   }
 709 }
 710
 711 void GridSmoother::computeFeet()
 712 {
 713   m_IdFoot.fill(-1, m_Grid->GetNumberOfPoints());
 714   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 715     if (m_Grid->GetCellType(id_cell) == VTK_WEDGE) {
 716       vtkIdType N_pts, *pts;
 717       m_Grid->GetCellPoints(id_cell, N_pts, pts);
 718       m_IdFoot[pts[3]] = pts[0];
 719       m_IdFoot[pts[4]] = pts[1];
 720       m_IdFoot[pts[5]] = pts[2];
 721       m_NodeMarked[pts[0]] = false;
 722       m_NodeMarked[pts[1]] = false;
 723       m_NodeMarked[pts[2]] = false;
 724     }
 725   }
 726 }
 727
 728 void GridSmoother::simpleNodeMovement(int i_nodes)
 729 {
 730   vtkIdType id_node = m_Part.globalNode(i_nodes);
 731   vtkIdType id_foot = m_IdFoot[id_node];
 732   vec3_t x_new(0,0,0), x_node;
 733   m_Grid->GetPoint(id_node, x_node.data());
 734
 735   vec3_t x_surf(0,0,0);
 736   if (m_SurfNode[id_node]) {
 737     int N = 0;
 738     for (int i = 0; i < m_Part.n2nLSize(i_nodes); ++i) {
 739       vtkIdType id_neigh = m_Part.n2nLG(i_nodes,i);
 740       if (m_SurfNode[id_neigh]) {
 741         vec3_t x_neigh;
 742         m_Grid->GetPoint(id_neigh, x_neigh.data());
 743         x_surf += x_neigh;
 744         ++N;
 745       }
 746     }
 747     if (N == 0) {
 748       EG_BUG;
 749     } else {
 750       x_surf *= 1.0/N;
 751     }
 752   }
 753
 754   if (id_foot != -1) {
 755     vec3_t x_foot, x_node;
 756     m_Grid->GetPoint(id_foot, x_foot.data());
 757     //double H = m_Blending*m_AbsoluteHeight + (1.0-m_Blending)*m_RelativeHeight*m_Height[id_foot];
 758     double H = m_Height[id_foot];
 759     x_new = x_foot + H*m_NodeNormal[id_foot];
 760   } else {
 761     if (m_SurfNode[id_node]) {
 762       x_new = x_surf;
 763     } else {
 764       if (m_Part.n2nLSize(i_nodes) == 0) {
 765         EG_BUG;
 766       }
 767       for (int i = 0; i < m_Part.n2nLSize(i_nodes); ++i) {
 768         vtkIdType id_neigh = m_Part.n2nLG(i_nodes,i);
 769         vec3_t x_neigh;
 770         m_Grid->GetPoint(id_neigh, x_neigh.data());
 771         x_new += x_neigh;
 772       }
 773       x_new *= 1.0/m_Part.n2nLSize(i_nodes);
 774     }
 775   }
 776   vec3_t Dx = x_new - x_node;
 777   correctDx(i_nodes, Dx);
 778   moveNode(i_nodes, Dx);
 779 }
 780
 781 void GridSmoother::operate()
 782 {
 783   if (m_FirstCall) {
 784     markNodes();
 785     computeNormals();
 786     computeFeet();
 787     computeHeights();
 788     m_FirstCall = false;
 789   }
 790   l2g_t nodes = m_Part.getNodes();
 791   for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
 792     if (m_NodeMarked[nodes[i_nodes]]) {
 793       simpleNodeMovement(i_nodes);
 794     }
 795   }
 796 }
 797
 798 void GridSmoother::writeDebugFile(QString file_name)
 799 {
 800   QVector<vtkIdType> bcells;
 801   QSet<int> bcs = getAllBoundaryCodes(m_Grid);
 802   getSurfaceCells(bcs, bcells, m_Grid);
 803   MeshPartition bpart(m_Grid);
 804   bpart.setCells(bcells);
 805   file_name = GuiMainWindow::pointer()->getCwd() + "/" + file_name + ".vtk";
 806   QFile file(file_name);
 807   file.open(QIODevice::WriteOnly);
 808   QTextStream f(&file);
 809   f << "# vtk DataFile Version 2.0\n";
 810   f << "m_NodeNormal\n";
 811   f << "ASCII\n";
 812   f << "DATASET UNSTRUCTURED_GRID\n";
 813   f << "POINTS " << bpart.getNumberOfNodes() << " float\n";
 814   for (int i = 0; i < bpart.getNumberOfNodes(); ++i) {
 815     vec3_t x;
 816     vtkIdType id_node = bpart.globalNode(i);
 817     m_Grid->GetPoint(id_node, x.data());
 818     f << x[0] << " " << x[1] << " " << x[2] << "\n";
 819   }
 820   f << "CELLS " << bpart.getNumberOfCells() << " ";
 821   int N = 0;
 822   for (int i = 0; i < bpart.getNumberOfCells(); ++ i) {
 823     vtkIdType id_cell = bpart.globalCell(i);
 824     vtkIdType N_pts, *pts;
 825     m_Grid->GetCellPoints(id_cell, N_pts, pts);
 826     N += 1 + N_pts;
 827   }
 828   f << N << "\n";
 829   for (int i = 0; i < bpart.getNumberOfCells(); ++ i) {
 830     vtkIdType id_cell = bpart.globalCell(i);
 831     vtkIdType N_pts, *pts;
 832     m_Grid->GetCellPoints(id_cell, N_pts, pts);
 833     QVector<int> nds(N_pts);
 834     f << N_pts;
 835     for (int j = 0; j < N_pts; ++j) {
 836       f << " " << bpart.localNode(pts[j]);
 837     }
 838     f << "\n";
 839   }
 840
 841   f << "CELL_TYPES " << bpart.getNumberOfCells() << "\n";
 842   for (int i = 0; i < bpart.getNumberOfCells(); ++ i) {
 843     vtkIdType id_cell = bpart.globalCell(i);
 844     vtkIdType type_cell = m_Grid->GetCellType(id_cell);
 845     f << type_cell << "\n";
 846   }
 847   f << "POINT_DATA " << bpart.getNumberOfNodes() << "\n";
 848   f << "VECTORS N float\n";
 849
 850   for (int i = 0; i < bpart.getNumberOfNodes(); ++i) {
 851     vtkIdType id_node = bpart.globalNode(i);
 852     f << m_NodeNormal[id_node][0] << " " << m_NodeNormal[id_node][1] << " " << m_NodeNormal[id_node][2] << "\n";
 853   }
 854 }
 855
 856