src/gridsmoother.cpp

   1 //
   2 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   3 // +                                                                      +
   4 // + This file is part of enGrid.                                         +
   5 // +                                                                      +
   6 // + Copyright 2008-2010 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
  26 #include "elements.h"
  27
  28 #include <QTime>
  29
  30 GridSmoother::GridSmoother()
  31 {
  32   m_NumIterations          = 5;
  33   m_NumRelaxations         = 5;
  34   m_NumBoundaryCorrections = 50;
  35
  36   getSet("boundary layer", "number of smoothing sub-iterations",       5,     m_NumIterations);
  37   getSet("boundary layer", "prism layer clearance",                    0.2,   m_LayerClearance);
  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",                              30,    m_RadarAngle);
  42   getSet("boundary layer", "maximal layer height in gaps",             0.2,   m_MaxHeightInGaps);
  43
  44   //m_CritAngle = GeometryTools::deg2rad(m_CritAngle);
  45 }
  46
  47 void GridSmoother::markNodes()
  48 {
  49   m_NodeMarked.fill(false,m_Grid->GetNumberOfPoints());
  50   QVector<bool> new_mark(m_Grid->GetNumberOfPoints());
  51   for (int i_iterations = 0; i_iterations < 4; ++i_iterations) {
  52     qCopy(m_NodeMarked.begin(),m_NodeMarked.end(),new_mark.begin());
  53     for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
  54       bool mark_cell = false;
  55       vtkIdType type_cell, N_pts, *pts;
  56       type_cell = m_Grid->GetCellType(id_cell);
  57       m_Grid->GetCellPoints(id_cell, N_pts, pts);
  58       if (type_cell == VTK_WEDGE) {
  59         mark_cell = true;
  60       } else {
  61         for (int i_pts = 0; i_pts < N_pts; ++i_pts) {
  62           if (m_NodeMarked[pts[i_pts]]) {
  63             mark_cell = true;
  64           }
  65         }
  66       }
  67       if (mark_cell) {
  68         for (int i_pts = 0; i_pts < N_pts; ++i_pts) {
  69           new_mark[pts[i_pts]] = true;
  70         }
  71       }
  72     }
  73     qCopy(new_mark.begin(),new_mark.end(),m_NodeMarked.begin());
  74   }
  75   m_NumMarkedNodes = 0;
  76   QVector<vtkIdType> nodes = m_Part.getNodes();
  77   foreach (vtkIdType id_node, nodes) {
  78     if (id_node < 0) EG_BUG;
  79     if (id_node > m_Grid->GetNumberOfPoints()) EG_BUG;
  80     if (m_NodeMarked[id_node]) {
  81       ++m_NumMarkedNodes;
  82     }
  83   }
  84 }
  85
  86 bool GridSmoother::noCollision(vtkIdType id_node)
  87 {
  88   bool cleared = true;
  89   vtkIdType id_foot = m_IdFoot[id_node];
  90   if (id_foot != -1) {
  91     QSet<vtkIdType> front_neigh;
  92     QSet<vtkIdType> own_prisms;
  93     for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
  94       vtkIdType id_neigh = m_Part.n2nGG(id_node, i);
  95       if (m_IdFoot[id_neigh] != -1) {
  96         front_neigh.insert(id_neigh);
  97       } else {
  98         QSet<vtkIdType> neigh_prisms;
  99         for (int j = 0; j < m_Part.n2cGSize(id_neigh); ++j) {
 100           vtkIdType id_cell = m_Part.n2cGG(id_neigh, j);
 101           if (m_Grid->GetCellType(id_cell) == VTK_WEDGE) {
 102             neigh_prisms.insert(id_cell);
 103           }
 104         }
 105         if (neigh_prisms.size() == 0) {
 106           for (int j = 0; j < m_Part.n2nGSize(id_neigh); ++j) {
 107             vtkIdType id_next_neigh = m_Part.n2nGG(id_neigh, j);
 108             if (m_IdFoot[id_next_neigh] != -1) {
 109               front_neigh.insert(id_next_neigh);
 110             }
 111           }
 112         }
 113       }
 114     }
 115     for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
 116       vtkIdType id_cell = m_Part.n2cGG(id_node, i);
 117       if (m_Grid->GetCellType(id_cell) == VTK_WEDGE) {
 118         own_prisms.insert(id_cell);
 119       }
 120     }
 121     foreach (vtkIdType id_neigh, front_neigh) {
 122       QSet<vtkIdType> neigh_prisms;
 123       for (int i = 0; i < m_Part.n2cGSize(id_neigh); ++i) {
 124         vtkIdType id_cell = m_Part.n2cGG(id_neigh, i);
 125         if (m_Grid->GetCellType(id_cell) == VTK_WEDGE) {
 126           neigh_prisms.insert(id_cell);
 127         }
 128       }
 129       QSet<vtkIdType> all_prisms = neigh_prisms + own_prisms;
 130       if (all_prisms.size() == neigh_prisms.size() + own_prisms.size()) { // here is a collision candidate!
 131         vec3_t n1(0,0,0);
 132         foreach (vtkIdType id_cell, own_prisms) {
 133           vtkIdType N_pts, *pts;
 134           m_Grid->GetCellPoints(id_cell, N_pts, pts);
 135           vec3_t a, b, c;
 136           m_Grid->GetPoint(pts[0], a.data());
 137           m_Grid->GetPoint(pts[2], b.data());
 138           m_Grid->GetPoint(pts[1], c.data());
 139           vec3_t u = b - a;
 140           vec3_t v = c - a;
 141           vec3_t n = u.cross(v);
 142           n.normalise();
 143           double alpha = GeometryTools::angle(u, v);
 144           n1 += alpha*n;
 145         }
 146         n1.normalise();
 147         vec3_t n2(0,0,0);
 148         foreach (vtkIdType id_cell, neigh_prisms) {
 149           vtkIdType N_pts, *pts;
 150           m_Grid->GetCellPoints(id_cell, N_pts, pts);
 151           vec3_t a, b, c;
 152           m_Grid->GetPoint(pts[0], a.data());
 153           m_Grid->GetPoint(pts[2], b.data());
 154           m_Grid->GetPoint(pts[1], c.data());
 155           vec3_t u = b - a;
 156           vec3_t v = c - a;
 157           vec3_t n = u.cross(v);
 158           n.normalise();
 159           double alpha = GeometryTools::angle(u, v);
 160           n2 += alpha*n;
 161         }
 162         n2.normalise();
 163         if (n1*n2 < 0) {
 164           vec3_t x1, x2, xf1, xf2;
 165           m_Grid->GetPoint(id_node, x1.data());
 166           m_Grid->GetPoint(id_neigh, x2.data());
 167           m_Grid->GetPoint(m_IdFoot[id_node], xf1.data());
 168           m_Grid->GetPoint(m_IdFoot[id_neigh], xf2.data());
 169           double l1 = (x1-xf1)*n1;
 170           double l2 = (x2-xf2)*n2;
 171           if (n1*(x2-x1) < m_LayerClearance*l1) {
 172             cleared = false;
 173           }
 174           if (n2*(x1-x2) < m_LayerClearance*l2) {
 175             cleared = false;
 176           }
 177         }
 178       }
 179     }
 180   }
 181   if (!cleared) {
 182     m_CollisionDetected = true;
 183   }
 184   return cleared;
 185 }
 186
 187 bool GridSmoother::setNewPosition(vtkIdType id_node, vec3_t x_new)
 188 {
 189   using namespace GeometryTools;
 190
 191   vec3_t x_old;
 192   m_Grid->GetPoint(id_node, x_old.data());
 193   m_Grid->GetPoints()->SetPoint(id_node, x_new.data());
 194   bool move = noCollision(id_node);
 195
 196   if (move) {
 197     Elements E;
 198
 199     l2g_t cells = m_Part.getCells();
 200     l2l_t n2c = m_Part.getN2C();
 201
 202     foreach (int i_cells, n2c[id_node]) {
 203       vtkIdType id_cell = cells[i_cells];
 204       vtkIdType type_cell = m_Grid->GetCellType(id_cell);
 205       if (type_cell == VTK_TETRA) {
 206         if (GeometryTools::cellVA(m_Grid, id_cell) < 0) {
 207           move = false;
 208         }
 209       }
 210
 211       if (type_cell == VTK_WEDGE && m_StrictPrismChecking) {
 212         vtkIdType N_pts, *pts;
 213         vec3_t xtet[4];
 214         m_Grid->GetCellPoints(id_cell, N_pts, pts);
 215         bool ok = true;
 216         for (int i = 0; i < 4; ++i) {     // variation
 217           ok = true;
 218           for (int j = 0; j < 3; ++j) {   // tetrahedron
 219             for (int k = 0; k < 4; ++k) { // node
 220               m_Grid->GetPoint(pts[E.priTet(i,j,k)], xtet[k].data());
 221             }
 222             double V = GeometryTools::tetraVol(xtet[0], xtet[1], xtet[2], xtet[3]);
 223             if (V <= 0) {
 224               ok = false;
 225             }
 226           }
 227           if (ok) {
 228             break;
 229           }
 230         }
 231         if (!ok) {
 232           move = false;
 233         }
 234       }
 235     }
 236   }
 237   if (!move) {
 238     m_Grid->GetPoints()->SetPoint(id_node, x_old.data());
 239   }
 240   return move;
 241 }
 242
 243 void GridSmoother::correctDx(int i_nodes, vec3_t &Dx)
 244 {
 245   l2g_t nodes = m_Part.getNodes();
 246   l2l_t n2c = m_Part.getN2C();
 247   for (int i_boundary_correction = 0; i_boundary_correction < m_NumBoundaryCorrections; ++i_boundary_correction) {
 248     foreach (vtkIdType id_cell, n2c[i_nodes]) {
 249       if (isSurface(id_cell, m_Grid)) {
 250         double A = GeometryTools::cellVA(m_Grid, id_cell);
 251         if (A > 1e-20) {
 252           vec3_t n = GeometryTools::cellNormal(m_Grid, id_cell);
 253           n.normalise();
 254           Dx -= (n*Dx)*n;
 255         }
 256       } else {
 257         if (m_Grid->GetCellType(id_cell) == VTK_WEDGE) {
 258           vtkIdType N_pts, *pts;
 259           m_Grid->GetCellPoints(id_cell, N_pts, pts);
 260           vtkIdType id_surf_node = -1;
 261           if (pts[3] == nodes[i_nodes]) id_surf_node = pts[0];
 262           if (pts[4] == nodes[i_nodes]) id_surf_node = pts[1];
 263           if (pts[5] == nodes[i_nodes]) id_surf_node = pts[2];
 264           if (id_surf_node != -1) {
 265             vec3_t x0,x1,x2;
 266             m_Grid->GetPoint(pts[0],x0.data());
 267             m_Grid->GetPoint(pts[1],x1.data());
 268             m_Grid->GetPoint(pts[2],x2.data());
 269             vec3_t a = x1-x0;
 270             vec3_t b = x2-x0;
 271             vec3_t c = b-a;
 272             double L = (a.abs()+b.abs()+c.abs())/3.0;
 273             vec3_t n = b.cross(a);
 274             n.normalise();
 275             vec3_t x_old;
 276             m_Grid->GetPoint(nodes[i_nodes],x_old.data());
 277             vec3_t x_new = x_old + Dx - x0;
 278             if ( (n*x_new) <= 0 ) {
 279               x_new -= (x_new*n)*n;
 280               x_new += 1e-4*L*n;
 281               x_new += x0;
 282               Dx = x_new - x_old;
 283             }
 284           }
 285         }
 286       }
 287     }
 288   }
 289 }
 290
 291 bool GridSmoother::moveNode(int i_nodes, vec3_t &Dx)
 292 {
 293   m_CollisionDetected = false;
 294   l2g_t nodes = m_Part.getNodes();
 295   vtkIdType id_node = nodes[i_nodes];
 296   vec3_t x_old;
 297   m_Grid->GetPoint(id_node, x_old.data());
 298   bool moved = false;
 299   for (int i_relaxation = 0; i_relaxation < m_NumRelaxations; ++i_relaxation) {
 300     if (setNewPosition(id_node, x_old + Dx)) {
 301       moved = true;
 302       break;
 303     }
 304     Dx *= 0.5;
 305   }
 306   return moved;
 307 }
 308
 309 void GridSmoother::computeNormals()
 310 {
 311   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
 312   m_NodeNormal.fill(vec3_t(0,0,0), m_Grid->GetNumberOfPoints());
 313   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 314     QSet<int> bcs;
 315     for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
 316       vtkIdType id_cell = m_Part.n2cGG(id_node, i);
 317       if (isSurface(id_cell, m_Grid)) {
 318         int bc = cell_code->GetValue(id_cell);
 319         if (m_BoundaryCodes.contains(bc)) {
 320           bcs.insert(bc);
 321         }
 322       }
 323     }
 324     int num_bcs = bcs.size();
 325     QVector<vec3_t> normal(num_bcs, vec3_t(0,0,0));
 326     QMap<int,int> bcmap;
 327     int i_bc = 0;
 328     foreach (int bc, bcs) {
 329       bcmap[bc] = i_bc;
 330       ++i_bc;
 331     }
 332     for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
 333       vtkIdType id_cell = m_Part.n2cGG(id_node, i);
 334       if (isSurface(id_cell, m_Grid)) {
 335         int bc = cell_code->GetValue(id_cell);
 336         if (m_BoundaryCodes.contains(bc)) {
 337           vtkIdType N_pts, *pts;
 338           m_Grid->GetCellPoints(id_cell, N_pts, pts);
 339           vec3_t a, b, c;
 340           for (int j = 0; j < N_pts; ++j) {
 341             if (pts[j] == id_node) {
 342               m_Grid->GetPoint(pts[j], a.data());
 343               if (j > 0) {
 344                 m_Grid->GetPoint(pts[j-1], b.data());
 345               } else {
 346                 m_Grid->GetPoint(pts[N_pts-1], b.data());
 347               }
 348               if (j < N_pts - 1) {
 349                 m_Grid->GetPoint(pts[j+1], c.data());
 350               } else {
 351                 m_Grid->GetPoint(pts[0], c.data());
 352               }
 353             }
 354           }
 355           vec3_t u = b - a;
 356           vec3_t v = c - a;
 357           double alpha = GeometryTools::angle(u, v);
 358           vec3_t n = u.cross(v);
 359           n.normalise();
 360           normal[bcmap[bc]] += alpha*n;
 361         }
 362       }
 363     }
 364     for (int i = 0; i < num_bcs; ++i) {
 365       normal[i].normalise();
 366     }
 367     if (num_bcs > 0) {
 368       if (num_bcs > 1) {
 369         if (num_bcs == 3) {
 370           for (int i = 0; i < num_bcs; ++i) {
 371             for (int j = i + 1; j < num_bcs; ++j) {
 372               vec3_t n = normal[i] + normal[j];
 373               n.normalise();
 374               m_NodeNormal[id_node] += n;
 375             }
 376           }
 377         } else {
 378           for (int i = 0; i < num_bcs; ++i) {
 379             m_NodeNormal[id_node] += normal[i];
 380           }
 381         }
 382       } else {
 383         m_NodeNormal[id_node] = normal[0];
 384       }
 385       m_NodeNormal[id_node].normalise();
 386     }
 387   }
 388   /*
 389   m_IsSharpNode.fill(false, m_Grid->GetNumberOfPoints());
 390   for (int i = 0; i < m_Part.getNumberOfCells(); ++i) {
 391     vtkIdType id_cell1 = m_Part.globalCell(i);
 392     if (isSurface(id_cell1, m_Grid)) {
 393       if (m_BoundaryCodes.contains(cell_code->GetValue(id_cell1))) {
 394         vec3_t n1 = cellNormal(m_Grid, id_cell1);
 395         n1.normalise();
 396         vtkIdType N_pts, *pts;
 397         m_Grid->GetCellPoints(id_cell1, N_pts, pts);
 398         for (int j = 0; j < m_Part.c2cLSize(i); ++j) {
 399           vtkIdType id_cell2 = m_Part.c2cLG(i, j);
 400           if (m_BoundaryCodes.contains(cell_code->GetValue(id_cell2))) {
 401             vec3_t n2 = cellNormal(m_Grid, id_cell2);
 402             n2.normalise();
 403             if (GeometryTools::angle(n1, n2) > m_CritAngle) {
 404               vtkIdType id_node1 = pts[j];
 405               vtkIdType id_node2 = pts[0];
 406               if (j < m_Part.c2cLSize(i) - 1) {
 407                 id_node2 = pts[j+1];
 408               }
 409               vec3_t x_node1, x_node2;
 410               m_Grid->GetPoint(id_node1, x_node1.data());
 411               m_Grid->GetPoint(id_node2, x_node2.data());
 412               vec3_t x_corner = 0.5*(x_node1 + x_node2);
 413               vec3_t x_cell1 = cellCentre(m_Grid, id_cell1);
 414               vec3_t x_cell2 = cellCentre(m_Grid, id_cell2);
 415               vec3_t v_cell1 = x_cell2 - x_cell1;
 416               vec3_t v_cell2 = x_cell1 - x_cell2;
 417               if (n1*v_cell1 > 0 || n2*v_cell2 > 0) {
 418                 m_IsSharpNode[id_node1] = true;
 419                 m_IsSharpNode[id_node2] = true;
 420               }
 421             }
 422           }
 423         }
 424       }
 425     }
 426   }
 427   */
 428   /*
 429   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 430     if (m_IsSharpNode[id_node]) {
 431       QVector<vec3_t> normals;
 432       for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
 433         vtkIdType id_cell = m_Part.n2cGG(id_node, i);
 434         if (isSurface(id_cell, m_Grid)) {
 435           if (m_BoundaryCodes.contains(cell_code->GetValue(id_cell))) {
 436             vec3_t n = cellNormal(m_Grid, id_cell);
 437             n.normalise();
 438             normals.push_back(n);
 439           }
 440         }
 441       }
 442       int N = 0;
 443       foreach (vec3_t n1, normals) {
 444         foreach (vec3_t n2, normals) {
 445           if (GeometryTools::angle(n1, n2) > m_CritAngle) {
 446             ++N;
 447           }
 448         }
 449       }
 450     }
 451   }
 452   */
 453
 454   relaxNormalVectors();
 455
 456 }
 457
 458 void GridSmoother::relaxNormalVectors()
 459 {
 460   m_SurfNode.fill(false, m_Grid->GetNumberOfPoints());
 461   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 462     for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
 463       if (isSurface(m_Part.n2cGG(id_node,i), m_Grid)) {
 464         m_SurfNode[id_node] = true;
 465         break;
 466       }
 467     }
 468   }
 469   EG_VTKDCC(vtkIntArray, bc, m_Grid, "cell_code");
 470   QVector<QSet<int> > n2bc(m_Grid->GetNumberOfPoints());
 471   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 472     if (isSurface(id_cell, m_Grid)) {
 473       vtkIdType N_pts, *pts;
 474       m_Grid->GetCellPoints(id_cell, N_pts, pts);
 475       for (int i = 0; i < N_pts; ++i) {
 476         if (m_SurfNode[pts[i]]) {
 477           n2bc[pts[i]].insert(bc->GetValue(id_cell));
 478         }
 479       }
 480     }
 481   }
 482   QVector<int> num_bcs(m_Grid->GetNumberOfPoints(),0);
 483   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 484     if (m_SurfNode[id_node]) {
 485       QList<vtkIdType> snap_points;
 486       foreach (int bc, n2bc[id_node]) {
 487         if (m_BoundaryCodes.contains(bc)) {
 488           ++num_bcs[id_node];
 489         }
 490       }
 491     }
 492   }
 493   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 494     if (num_bcs[id_node] == 0) {
 495       m_SurfNode[id_node] = false;
 496       n2bc[id_node].clear();
 497     }
 498   }
 499   for (int iter = 0; iter < m_NumNormalRelaxations; ++iter) {
 500     QVector<vec3_t> n_new(m_NodeNormal.size(), vec3_t(0,0,0));
 501     for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 502       if (m_SurfNode[id_node] ) {
 503         QList<vtkIdType> snap_points;
 504         for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
 505           vtkIdType id_neigh = m_Part.n2nGG(id_node,i);
 506           if (num_bcs[id_node] <= num_bcs[id_neigh]) {
 507             if (!m_SurfNode[id_neigh]) {
 508               cout << id_node << ',' << m_Part.n2nGG(id_node,i) << ',' << num_bcs[id_node] << ',' << num_bcs[id_neigh] << endl;
 509               EG_BUG;
 510             }
 511             snap_points.append(id_neigh);
 512           }
 513         }
 514         if (snap_points.size() > 0) {
 515           n_new[id_node] = vec3_t(0,0,0);
 516           foreach (vtkIdType id_snap, snap_points) {
 517             n_new[id_node] += m_NodeNormal[id_snap];
 518           }
 519           n_new[id_node].normalise();
 520         } else {
 521           n_new[id_node] = m_NodeNormal[id_node];
 522         }
 523         if (n_new[id_node].abs() < 0.1) {
 524           cout << id_node << ',' << n_new[id_node] << ',' << num_bcs[id_node] << ',' << n2bc[id_node] << endl;
 525           EG_BUG;
 526         }
 527       }
 528     }
 529     m_NodeNormal = n_new;
 530     correctNormalVectors();
 531     QString num;
 532     num.setNum(iter);
 533     num = "normals_" + num;
 534     writeDebugFile(num);
 535   }
 536
 537   //writeDebugFile("normals");
 538   //EG_BUG;
 539 }
 540
 541 void GridSmoother::computeHeights()
 542 {
 543   // first pass (intial heigh)
 544   EG_VTKDCN(vtkDoubleArray, cl, m_Grid, "node_meshdensity_desired" );
 545   m_Height.fill(0, m_Grid->GetNumberOfPoints());
 546   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 547     if (m_SurfNode[id_node]) {
 548       m_Height[id_node] = cl->GetValue(id_node);
 549       // if undefined: compute height from surrounding edges
 550       if (m_Height[id_node] < 1e-99) {
 551         m_Height[id_node] = 0;
 552         int N = 0;
 553         vec3_t x;
 554         m_Grid->GetPoint(id_node, x.data());
 555         for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
 556           vtkIdType id_neigh_node = m_Part.n2nGG(id_node, i);
 557           if (m_SurfNode[id_neigh_node]) {
 558             ++N;
 559             vec3_t xn;
 560             m_Grid->GetPoint(id_neigh_node, xn.data());
 561             m_Height[id_node] += (x - xn).abs();
 562           }
 563         }
 564         if (N == 0) {
 565           EG_BUG;
 566         }
 567         m_Height[id_node] /= N;
 568       }
 569     }
 570   }
 571
 572   // second pass (correct with absolute height if required)
 573   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 574     if (m_SurfNode[id_node]) {
 575       m_Height[id_node] = m_Blending*m_AbsoluteHeight + (1.0-m_Blending)*m_RelativeHeight*m_Height[id_node];
 576     }
 577   }
 578
 579   // third pass (gaps)
 580   QList<vtkIdType> surf_nodes;
 581   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 582     if (m_SurfNode[id_node]) {
 583       surf_nodes.append(id_node);
 584     }
 585   }
 586   foreach (vtkIdType id_node1, surf_nodes) {
 587     foreach (vtkIdType id_node2, surf_nodes) {
 588       if (id_node1 != id_node2) {
 589         const vec3_t& n1 = m_NodeNormal[id_node1];
 590         const vec3_t& n2 = m_NodeNormal[id_node2];
 591         vec3_t x1, x2;
 592         m_Grid->GetPoint(id_node1, x1.data());
 593         m_Grid->GetPoint(id_node2, x2.data());
 594         vec3_t Dx = x2 - x1;
 595         double a = Dx*n1;
 596         if (a > 0) {
 597           double b = Dx.abs();
 598           double alpha = 180.0/M_PI*acos(a/b);
 599           if (alpha < m_RadarAngle) {
 600             m_Height[id_node1] = min(m_Height[id_node1], m_MaxHeightInGaps*a);
 601           }
 602         }
 603       }
 604     }
 605   }
 606
 607   // fourth pass (smoothing)
 608   for (int iter = 0; iter < m_NumHeightRelaxations; ++iter) {
 609     QVector<double> h_new = m_Height;
 610     for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 611       if (m_SurfNode[id_node]) {
 612         int N = 0;
 613         h_new[id_node] = 0;
 614         for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
 615           vtkIdType id_neigh = m_Part.n2nGG(id_node,i);
 616           if (m_SurfNode[id_neigh]) {
 617             h_new[id_node] += m_Height[id_neigh];
 618             ++N;
 619           }
 620         }
 621         if (N == 0) {
 622           EG_BUG;
 623         }
 624         h_new[id_node] /= N;
 625       }
 626     }
 627     m_Height = h_new;
 628   }
 629 }
 630
 631 void GridSmoother::correctNormalVectors()
 632 {
 633   EG_VTKDCC(vtkIntArray, bc, m_Grid, "cell_code");
 634   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 635     if (m_NodeNormal[id_node].abs() > 0.1) {
 636       for (int iter = 0; iter < m_NumBoundaryCorrections; ++iter) {
 637         for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
 638           vtkIdType id_cell = m_Part.n2cGG(id_node,i);
 639           if (isSurface(id_cell, m_Grid)) {
 640             if (!m_BoundaryCodes.contains(bc->GetValue(id_cell))) {
 641               vec3_t v = m_NodeNormal[id_node];
 642               vec3_t n = GeometryTools::cellNormal(m_Grid, id_cell);
 643               n.normalise();
 644               v -= (n*m_NodeNormal[id_node])*n;
 645               v.normalise();
 646               m_NodeNormal[id_node] = v;
 647             }
 648           }
 649         }
 650       }
 651     }
 652   }
 653 }
 654
 655 void GridSmoother::computeFeet()
 656 {
 657   m_IdFoot.fill(-1, m_Grid->GetNumberOfPoints());
 658   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 659     if (m_Grid->GetCellType(id_cell) == VTK_WEDGE) {
 660       vtkIdType N_pts, *pts;
 661       m_Grid->GetCellPoints(id_cell, N_pts, pts);
 662       m_IdFoot[pts[3]] = pts[0];
 663       m_IdFoot[pts[4]] = pts[1];
 664       m_IdFoot[pts[5]] = pts[2];
 665       m_NodeMarked[pts[0]] = false;
 666       m_NodeMarked[pts[1]] = false;
 667       m_NodeMarked[pts[2]] = false;
 668     }
 669   }
 670 }
 671
 672 void GridSmoother::simpleNodeMovement(int i_nodes)
 673 {
 674   vtkIdType id_node = m_Part.globalNode(i_nodes);
 675   vtkIdType id_foot = m_IdFoot[id_node];
 676   vec3_t x_new(0,0,0), x_node;
 677   m_Grid->GetPoint(id_node, x_node.data());
 678
 679   vec3_t x_surf(0,0,0);
 680   if (m_SurfNode[id_node]) {
 681     int N = 0;
 682     for (int i = 0; i < m_Part.n2nLSize(i_nodes); ++i) {
 683       vtkIdType id_neigh = m_Part.n2nLG(i_nodes,i);
 684       if (m_SurfNode[id_neigh]) {
 685         vec3_t x_neigh;
 686         m_Grid->GetPoint(id_neigh, x_neigh.data());
 687         x_surf += x_neigh;
 688         ++N;
 689       }
 690     }
 691     if (N == 0) {
 692       EG_BUG;
 693     } else {
 694       x_surf *= 1.0/N;
 695     }
 696   }
 697
 698   if (id_foot != -1) {
 699     vec3_t x_foot, x_node;
 700     m_Grid->GetPoint(id_foot, x_foot.data());
 701     //double H = m_Blending*m_AbsoluteHeight + (1.0-m_Blending)*m_RelativeHeight*m_Height[id_foot];
 702     double H = m_Height[id_foot];
 703     x_new = x_foot + H*m_NodeNormal[id_foot];
 704   } else {
 705     if (m_SurfNode[id_node]) {
 706       x_new = x_surf;
 707     } else {
 708       if (m_Part.n2nLSize(i_nodes) == 0) {
 709         EG_BUG;
 710       }
 711       for (int i = 0; i < m_Part.n2nLSize(i_nodes); ++i) {
 712         vtkIdType id_neigh = m_Part.n2nLG(i_nodes,i);
 713         vec3_t x_neigh;
 714         m_Grid->GetPoint(id_neigh, x_neigh.data());
 715         x_new += x_neigh;
 716       }
 717       x_new *= 1.0/m_Part.n2nLSize(i_nodes);
 718     }
 719   }
 720   vec3_t Dx = x_new - x_node;
 721   correctDx(i_nodes, Dx);
 722   moveNode(i_nodes, Dx);
 723 }
 724
 725 void GridSmoother::operate()
 726 {
 727   markNodes();
 728   computeNormals();
 729   computeFeet();
 730   computeHeights();
 731   l2g_t nodes = m_Part.getNodes();
 732   for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
 733     if (m_NodeMarked[nodes[i_nodes]]) {
 734       simpleNodeMovement(i_nodes);
 735     }
 736   }
 737 }
 738
 739 void GridSmoother::writeDebugFile(QString file_name)
 740 {
 741   QVector<vtkIdType> bcells;
 742   getSurfaceCells(m_BoundaryCodes, bcells, m_Grid);
 743   MeshPartition bpart(m_Grid);
 744   bpart.setCells(bcells);
 745   /*
 746   QVector<vtkIdType> foot2field(bpart.getNumberOfNodes(), -1);
 747   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 748     vtkIdType id_foot = m_IdFoot[id_node];
 749     if (id_foot != -1) {
 750       if (bpart.localNode(id_foot) == -1) {
 751         EG_BUG;
 752       }
 753       foot2field[bpart.localNode(id_foot)] = id_node;
 754     }
 755   }
 756   */
 757   file_name = GuiMainWindow::pointer()->getCwd() + "/" + file_name + ".vtk";
 758   QFile file(file_name);
 759   file.open(QIODevice::WriteOnly);
 760   QTextStream f(&file);
 761   f << "# vtk DataFile Version 2.0\n";
 762   f << "m_NodeNormal\n";
 763   f << "ASCII\n";
 764   f << "DATASET UNSTRUCTURED_GRID\n";
 765   f << "POINTS " << bpart.getNumberOfNodes() << " float\n";
 766   for (int i = 0; i < bpart.getNumberOfNodes(); ++i) {
 767     vec3_t x;
 768     vtkIdType id_node = bpart.globalNode(i);
 769     m_Grid->GetPoint(id_node, x.data());
 770     f << x[0] << " " << x[1] << " " << x[2] << "\n";
 771   }
 772   f << "CELLS " << bpart.getNumberOfCells() << " " << 4*bpart.getNumberOfCells() << "\n";
 773   for (int i = 0; i < bpart.getNumberOfCells(); ++ i) {
 774     vtkIdType id_cell = bpart.globalCell(i);
 775     vtkIdType N_pts, *pts;
 776     if (m_Grid->GetCellType(id_cell) != VTK_TRIANGLE) {
 777       EG_BUG;
 778     }
 779     m_Grid->GetCellPoints(id_cell, N_pts, pts);
 780     QVector<int> nds(3);
 781     for (int j = 0; j < 3; ++j) {
 782       nds[j] = bpart.localNode(pts[j]);
 783     }
 784     f << "3 " << nds[0] << " " << nds[1] << " " << nds[2] << "\n";
 785   }
 786
 787   f << "CELL_TYPES " << bpart.getNumberOfCells() << "\n";
 788   for (int i = 0; i < bpart.getNumberOfCells(); ++ i) {
 789     f << VTK_TRIANGLE << "\n";
 790   }
 791   f << "POINT_DATA " << bpart.getNumberOfNodes() << "\n";
 792   f << "VECTORS N float\n";
 793
 794   for (int i = 0; i < bpart.getNumberOfNodes(); ++i) {
 795     vtkIdType id_node = bpart.globalNode(i);
 796     f << m_NodeNormal[id_node][0] << " " << m_NodeNormal[id_node][1] << " " << m_NodeNormal[id_node][2] << "\n";
 797   }
 798 }
 799
 800