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