src/gridsmoother.cpp

   1 //
   2 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   3 // +                                                                      +
   4 // + This file is part of enGrid.                                         +
   5 // +                                                                      +
   6 // + Copyright 2008,2009 Oliver Gloth                                     +
   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   N_iterations = 5;
  33   N_relaxations = 5;
  34   N_boundary_corrections = 20;
  35   N_search = 10;
  36   L_search = 0.05;
  37   smooth_prisms = true;
  38   dbg = false;
  39   F_old = 0;
  40   F_new = 0;
  41
  42   getSet("boundary layer", "relative height of boundary layer", 0.75, m_RelativeHeight);
  43   getSet("boundary layer", "under relaxation",                  1.00, m_UnderRelaxation);
  44   getSet("boundary layer", "maximal relative edge length",      1.50, m_MaxRelLength);
  45
  46   getSet("boundary layer", "number of smoothing sub-iterations", 5, N_iterations);
  47
  48   getSet("boundary layer", "use strict prism checking", true, m_StrictPrismChecking);
  49   getSet("boundary layer", "write debug file",          true, m_WriteDebugFile);
  50
  51   getErrSet("boundary layer", "layer height error",      200.0,  2.0, 2.0, 0.10, m_HeightError);
  52   getErrSet("boundary layer", "edge length error" ,      200.0,  0.1, 2.0, 0.10, m_EdgeLengthError);
  53   getErrSet("boundary layer", "edge direction error" ,     0.0,  0.0, 2.0, 0.00, m_EdgeDirectionError);
  54   getErrSet("boundary layer", "tetra error",              50.0,  0.1, 2.0, 0.05, m_TetraError);
  55   getErrSet("boundary layer", "parallel edges error",      1.0, 10.0, 2.0, 0.00, m_ParallelEdgesError);
  56   getErrSet("boundary layer", "parallel faces error",      1.0, 20.0, 2.0, 0.00, m_ParallelFacesError);
  57   getErrSet("boundary layer", "sharp nodes error",         0.0,  0.0, 2.0, 0.00, m_SharpNodesError);
  58   getErrSet("boundary layer", "sharp edges error",         0.0,  0.0, 2.0, 0.00, m_SharpEdgesError);
  59   getErrSet("boundary layer", "similar face area error",   0.0,  1.0, 2.0, 0.00, m_SimilarFaceAreaError);
  60   getErrSet("boundary layer", "feature line error",        0.0,  0.0, 2.0, 0.00, m_FeatureLineError);
  61
  62   m_CritAngle = GeometryTools::deg2rad(450.0);
  63
  64 }
  65
  66 void GridSmoother::computeNormals()
  67 {
  68   EG_VTKDCC(vtkIntArray, cell_code, grid, "cell_code");
  69   m_NodeNormal.fill(vec3_t(0,0,0), grid->GetNumberOfPoints());
  70   for (vtkIdType id_node = 0; id_node < grid->GetNumberOfPoints(); ++id_node) {
  71     QSet<int> bcs;
  72     for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
  73       vtkIdType id_cell = m_Part.n2cGG(id_node, i);
  74       if (isSurface(id_cell, grid)) {
  75         int bc = cell_code->GetValue(id_cell);
  76         if (m_BoundaryCodes.contains(bc)) {
  77           bcs.insert(bc);
  78         }
  79       }
  80     }
  81     int num_bcs = bcs.size();
  82     QVector<vec3_t> normal(num_bcs, vec3_t(0,0,0));
  83     QMap<int,int> bcmap;
  84     int i_bc = 0;
  85     foreach (int bc, bcs) {
  86       bcmap[bc] = i_bc;
  87       ++i_bc;
  88     }
  89     for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
  90       vtkIdType id_cell = m_Part.n2cGG(id_node, i);
  91       if (isSurface(id_cell, grid)) {
  92         int bc = cell_code->GetValue(id_cell);
  93         if (m_BoundaryCodes.contains(bc)) {
  94           vtkIdType N_pts, *pts;
  95           grid->GetCellPoints(id_cell, N_pts, pts);
  96           vec3_t a, b, c;
  97           for (int j = 0; j < N_pts; ++j) {
  98             if (pts[j] == id_node) {
  99               grid->GetPoint(pts[j], a.data());
 100               if (j > 0) {
 101                 grid->GetPoint(pts[j-1], b.data());
 102               } else {
 103                 grid->GetPoint(pts[N_pts-1], b.data());
 104               }
 105               if (j < N_pts - 1) {
 106                 grid->GetPoint(pts[j+1], c.data());
 107               } else {
 108                 grid->GetPoint(pts[0], c.data());
 109               }
 110             }
 111           }
 112           vec3_t u = b - a;
 113           vec3_t v = c - a;
 114           double alpha = GeometryTools::angle(u, v);
 115           vec3_t n = u.cross(v);
 116           n.normalise();
 117           normal[bcmap[bc]] += alpha*n;
 118         }
 119       }
 120     }
 121     for (int i = 0; i < num_bcs; ++i) {
 122       normal[i].normalise();
 123     }
 124     if (num_bcs > 0) {
 125       if (num_bcs > 1) {
 126         for (int i = 0; i < num_bcs; ++i) {
 127           for (int j = i + 1; j < num_bcs; ++j) {
 128             vec3_t n = normal[i] + normal[j];
 129             n.normalise();
 130             m_NodeNormal[id_node] += n;
 131           }
 132         }
 133       } else {
 134         m_NodeNormal[id_node] = normal[0];
 135       }
 136       m_NodeNormal[id_node].normalise();
 137     }
 138   }
 139 }
 140
 141 void GridSmoother::markNodes()
 142 {
 143   node_marked.fill(false,grid->GetNumberOfPoints());
 144   QVector<bool> new_mark(grid->GetNumberOfPoints());
 145   for (int i_iterations = 0; i_iterations < 4; ++i_iterations) {
 146     qCopy(node_marked.begin(),node_marked.end(),new_mark.begin());
 147     for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
 148       bool mark_cell = false;
 149       vtkIdType type_cell, N_pts, *pts;
 150       type_cell = grid->GetCellType(id_cell);
 151       grid->GetCellPoints(id_cell, N_pts, pts);
 152       if (type_cell == VTK_WEDGE) {
 153         mark_cell = true;
 154       } else {
 155         for (int i_pts = 0; i_pts < N_pts; ++i_pts) {
 156           if (node_marked[pts[i_pts]]) {
 157             mark_cell = true;
 158           }
 159         }
 160       }
 161       if (mark_cell) {
 162         for (int i_pts = 0; i_pts < N_pts; ++i_pts) {
 163           new_mark[pts[i_pts]] = true;
 164         }
 165       }
 166     }
 167     qCopy(new_mark.begin(),new_mark.end(),node_marked.begin());
 168   }
 169   N_marked_nodes = 0;
 170   QVector<vtkIdType> nodes = m_Part.getNodes();
 171   foreach (vtkIdType id_node, nodes) {
 172     if (id_node < 0) EG_BUG;
 173     if (id_node > grid->GetNumberOfPoints()) EG_BUG;
 174     if (node_marked[id_node]) {
 175       ++N_marked_nodes;
 176     }
 177   }
 178 }
 179
 180 bool GridSmoother::setNewPosition(vtkIdType id_node, vec3_t x_new)
 181 {
 182   using namespace GeometryTools;
 183
 184   vec3_t x_old;
 185   grid->GetPoint(id_node, x_old.data());
 186   grid->GetPoints()->SetPoint(id_node, x_new.data());
 187   bool move = true;
 188   Elements E;
 189   if (move) {
 190     l2g_t cells = getPartCells();
 191     l2l_t n2c   = getPartN2C();
 192     foreach (int i_cells, n2c[id_node]) {
 193       vtkIdType id_cell = cells[i_cells];
 194       vtkIdType type_cell = grid->GetCellType(id_cell);
 195       if (type_cell == VTK_TETRA) {
 196         if (GeometryTools::cellVA(grid, id_cell) < 0) {
 197           move = false;
 198           //if (dbg) cout << id_node << " : tetra negative" << endl;
 199         }
 200       }
 201       if (type_cell == VTK_WEDGE && m_StrictPrismChecking) {
 202         vtkIdType N_pts, *pts;
 203         vec3_t xtet[4];
 204         grid->GetCellPoints(id_cell, N_pts, pts);
 205         bool ok = true;
 206         for (int i = 0; i < 4; ++i) {     // variation
 207           ok = true;
 208           for (int j = 0; j < 3; ++j) {   // tetrahedron
 209             for (int k = 0; k < 4; ++k) { // node
 210               grid->GetPoint(pts[E.priTet(i,j,k)], xtet[k].data());
 211             }
 212             if (GeometryTools::tetraVol(xtet[0], xtet[1], xtet[2], xtet[3]) < 0) {
 213               ok = false;
 214               //if (dbg) cout << id_node << " : prism negative" << endl;
 215             }
 216           }
 217           if (ok) {
 218             break;
 219           }
 220         }
 221         if (!ok) {
 222           move = false;
 223         }
 224       }
 225     }
 226   }
 227   if (!move) {
 228     grid->GetPoints()->SetPoint(id_node, x_old.data());
 229   }
 230   return move;
 231 }
 232
 233 void GridSmoother::correctDx(int i_nodes, vec3_t &Dx)
 234 {
 235   l2g_t nodes = m_Part.getNodes();
 236   l2l_t n2c   = m_Part.getN2C();
 237
 238   for (int i_boundary_correction = 0; i_boundary_correction < N_boundary_corrections; ++i_boundary_correction) {
 239     foreach (vtkIdType id_cell, n2c[i_nodes]) {
 240       if (isSurface(id_cell, grid)) {
 241         double A = GeometryTools::cellVA(grid, id_cell);
 242         if (A > 1e-20) {
 243           vec3_t n = GeometryTools::cellNormal(grid, id_cell);
 244           n.normalise();
 245           Dx -= (n*Dx)*n;
 246         }
 247       } else {
 248         if (grid->GetCellType(id_cell) == VTK_WEDGE) {
 249           vtkIdType N_pts, *pts;
 250           grid->GetCellPoints(id_cell, N_pts, pts);
 251           vtkIdType id_surf_node = -1;
 252           if (pts[3] == nodes[i_nodes]) id_surf_node = pts[0];
 253           if (pts[4] == nodes[i_nodes]) id_surf_node = pts[1];
 254           if (pts[5] == nodes[i_nodes]) id_surf_node = pts[2];
 255           if (id_surf_node != -1) {
 256             vec3_t x0,x1,x2;
 257             grid->GetPoint(pts[0],x0.data());
 258             grid->GetPoint(pts[1],x1.data());
 259             grid->GetPoint(pts[2],x2.data());
 260             vec3_t a = x1-x0;
 261             vec3_t b = x2-x0;
 262             vec3_t c = b-a;
 263             double L = (a.abs()+b.abs()+c.abs())/3.0;
 264             vec3_t n = b.cross(a);
 265             n.normalise();
 266             vec3_t x_old;
 267             grid->GetPoint(nodes[i_nodes],x_old.data());
 268             vec3_t x_new = x_old + Dx - x0;
 269             if ( (n*x_new) <= 0 ) {
 270               x_new -= (x_new*n)*n;
 271               x_new += 1e-4*L*n;
 272               x_new += x0;
 273               Dx = x_new - x_old;
 274             }
 275           }
 276         }
 277       }
 278     }
 279   }
 280 }
 281
 282 bool GridSmoother::moveNode(int i_nodes, vec3_t &Dx)
 283 {
 284   l2g_t nodes = getPartNodes();
 285   vtkIdType id_node = nodes[i_nodes];
 286   vec3_t x_old;
 287   grid->GetPoint(id_node, x_old.data());
 288
 289   if (m_IdFoot[id_node] != -1) {
 290     vec3_t x_foot;
 291     grid->GetPoint(m_IdFoot[id_node], x_foot.data());
 292     Dx += x_old - x_foot;
 293     if (Dx.abs() > m_MaxRelLength*m_L[id_node]) {
 294       Dx.normalise();
 295       Dx *= m_MaxRelLength*m_L[id_node];
 296     }
 297     Dx -= x_old - x_foot;
 298   }
 299
 300   bool moved = false;
 301   for (int i_relaxation = 0; i_relaxation < N_relaxations; ++i_relaxation) {
 302     if (setNewPosition(id_node, x_old + Dx)) {
 303       moved = true;
 304       break;
 305     }
 306     Dx *= 0.1;
 307   }
 308   return moved;
 309 }
 310
 311 double GridSmoother::func(vec3_t x)
 312 {
 313   l2g_t nodes = getPartNodes();
 314   l2l_t n2c   = getPartN2C();
 315   l2g_t cells = getPartCells();
 316   l2l_t c2c   = getPartC2C();
 317
 318   vec3_t x_old;
 319   grid->GetPoint(nodes[i_nodes_opt], x_old.data());
 320   grid->GetPoints()->SetPoint(nodes[i_nodes_opt], x.data());
 321   double f       = 0;
 322   double f13     = 1.0/3.0;
 323   double f14     = 0.25;
 324
 325   vec3_t n_node(1,0,0);
 326   QList<vec3_t> n_pri;
 327
 328   EG_VTKDCN(vtkDoubleArray, cl, grid, "node_meshdensity_desired" );
 329
 330   int N_prisms = 0;
 331
 332   QSet<edge_t> edges;
 333
 334   foreach (int i_cells, n2c[i_nodes_opt]) {
 335     vtkIdType id_cell = cells[i_cells];
 336     if (isVolume(id_cell, grid)) {
 337       vtkIdType type_cell = grid->GetCellType(id_cell);
 338       vtkIdType N_pts, *pts;
 339       grid->GetCellPoints(id_cell, N_pts, pts);
 340       QVector<vec3_t> xn(N_pts);
 341       for (int i_pts = 0; i_pts < N_pts; ++i_pts) {
 342         grid->GetPoint(pts[i_pts],xn[i_pts].data());
 343       }
 344       if (type_cell == VTK_TETRA) {
 345         if (m_TetraError.active()) {
 346           double L = 0;
 347           L = max(L, (xn[0]-xn[1]).abs());
 348           L = max(L, (xn[0]-xn[2]).abs());
 349           L = max(L, (xn[0]-xn[3]).abs());
 350           L = max(L, (xn[1]-xn[2]).abs());
 351           L = max(L, (xn[1]-xn[3]).abs());
 352           L = max(L, (xn[2]-xn[3]).abs());
 353           double H = sqrt(2.0/3.0)*L;
 354
 355           vec3_t n012 = GeometryTools::triNormal(xn[0], xn[1], xn[2]);
 356           vec3_t n031 = GeometryTools::triNormal(xn[0], xn[3], xn[1]);
 357           vec3_t n023 = GeometryTools::triNormal(xn[0], xn[2], xn[3]);
 358           vec3_t n213 = GeometryTools::triNormal(xn[2], xn[1], xn[3]);
 359           n012.normalise();
 360           n031.normalise();
 361           n023.normalise();
 362           n213.normalise();
 363           double h = 1e99;
 364           h = min(h, fabs((xn[0]-xn[3])*n012));
 365           h = min(h, fabs((xn[0]-xn[2])*n031));
 366           h = min(h, fabs((xn[0]-xn[1])*n023));
 367           h = min(h, fabs((xn[0]-xn[3])*n213));
 368           h /= H;
 369           f += m_TetraError(h);
 370         }
 371       }
 372       if (type_cell == VTK_WEDGE) {
 373         ++N_prisms;
 374         vec3_t a  = xn[2]-xn[0];
 375         vec3_t b  = xn[1]-xn[0];
 376         vec3_t c  = xn[5]-xn[3];
 377         vec3_t d  = xn[4]-xn[3];
 378         vec3_t n_face[5];
 379         vec3_t x_face[5];
 380         n_face[0] = GeometryTools::triNormal(xn[0],xn[1],xn[2]);
 381         n_face[1] = GeometryTools::triNormal(xn[3],xn[5],xn[4]);
 382         n_face[2] = GeometryTools::quadNormal(xn[0],xn[3],xn[4],xn[1]);
 383         n_face[3] = GeometryTools::quadNormal(xn[1],xn[4],xn[5],xn[2]);
 384         n_face[4] = GeometryTools::quadNormal(xn[0],xn[2],xn[5],xn[3]);
 385         x_face[0] = f13*(xn[0]+xn[1]+xn[2]);
 386         x_face[1] = f13*(xn[3]+xn[4]+xn[5]);
 387         x_face[2] = f14*(xn[0]+xn[3]+xn[4]+xn[1]);
 388         x_face[3] = f14*(xn[1]+xn[4]+xn[5]+xn[2]);
 389         x_face[4] = f14*(xn[0]+xn[2]+xn[5]+xn[3]);
 390
 391         double A1 = 0.5*n_face[0].abs();
 392         double A2 = 0.5*n_face[1].abs();
 393         for (int i_face = 0; i_face < 5; ++i_face) {
 394           n_face[i_face].normalise();
 395         }
 396         m_IdFoot[pts[3]] = pts[0];
 397         m_IdFoot[pts[4]] = pts[1];
 398         m_IdFoot[pts[5]] = pts[2];
 399         vec3_t v0 = xn[0]-xn[3];
 400         vec3_t v1 = xn[1]-xn[4];
 401         vec3_t v2 = xn[2]-xn[5];
 402         double h0 = v0*n_face[0];
 403         double h1 = v1*n_face[0];
 404         double h2 = v2*n_face[0];
 405         if (m_HeightError.active() || m_EdgeLengthError.active() || m_EdgeDirectionError.active()) {
 406           double L = 0;
 407           int i_foot = -1;
 408           if (nodes[i_nodes_opt] == pts[3]) {
 409             i_foot = 0;
 410           }
 411           if (nodes[i_nodes_opt] == pts[4]) {
 412             i_foot = 1;
 413           }
 414           if (nodes[i_nodes_opt] == pts[5]) {
 415             i_foot = 2;
 416           }
 417           if (i_foot != -1) {
 418             L = m_RelativeHeight*cl->GetValue(pts[i_foot]);
 419             n_node = xn[i_foot + 3] - xn[i_foot];
 420             n_pri.append(n_face[1]);
 421           }
 422           double h = 0;
 423           double l = 0;
 424           vec3_t ve(0,0,0);
 425           if (i_foot != -1) {
 426             m_L[nodes[i_nodes_opt]] = L;
 427             if (i_foot == 0) {
 428               h = h0/L;
 429               l = -(v0*m_NodeNormal[pts[0]])/L;
 430               ve = v0;
 431             } else if (i_foot == 1) {
 432               h = h1/L;
 433               l = -(v1*m_NodeNormal[pts[1]])/L;
 434               ve = v1;
 435             } else if (i_foot == 2) {
 436               h = h2/L;
 437               l = -(v2*m_NodeNormal[pts[2]])/L;
 438               ve = v2;
 439             }
 440             ve *= -1;
 441             ve.normalise();
 442             f += m_HeightError(h);
 443             f += m_EdgeLengthError(l);
 444             f += m_EdgeDirectionError(angleX(m_NodeNormal[pts[i_foot]],ve));
 445           }
 446           if (m_ParallelEdgesError.active()) {
 447             if ((h0 > 0.01*L) && (h1 > 0.01*L) && (h2 > 0.01*L)) {
 448               v0.normalise();
 449               v1.normalise();
 450               v2.normalise();
 451               f += f13*m_ParallelEdgesError(angleX(v0,v1));
 452               f += f13*m_ParallelEdgesError(angleX(v0,v2));
 453               f += f13*m_ParallelEdgesError(angleX(v1,v2));
 454             }
 455           }
 456           if (m_ParallelFacesError.active()) {
 457             if ((h0 > 0.01*L) && (h1 > 0.01*L) && (h2 > 0.01*L)) {
 458               f += m_ParallelFacesError(angleX(-1*n_face[0],n_face[1]));
 459             }
 460           }
 461           if (m_SimilarFaceAreaError.active()) {
 462             if ((h0 > 0.01*L) && (h1 > 0.01*L) && (h2 > 0.01*L)) {
 463               f += m_SimilarFaceAreaError(fabs(1.0 - 2*(A1-A2)/(A1+A2)));
 464             }
 465           }
 466         }
 467
 468         if (m_SharpEdgesError.active()) {
 469           for (int j = 2; j <= 4; ++j) {
 470             vtkIdType id_ncell = c2c[id_cell][j];
 471             if (id_ncell != -1) {
 472               if (grid->GetCellType(id_ncell) == VTK_WEDGE) {
 473                 vtkIdType N_pts, *pts;
 474                 grid->GetCellPoints(id_ncell, N_pts, pts);
 475                 QVector<vec3_t> x(3);
 476                 for (int i_pts = 3; i_pts <= 5; ++i_pts) {
 477                   grid->GetPoint(pts[i_pts],x[i_pts-3].data());
 478                 }
 479                 vec3_t n = GeometryTools::triNormal(x[0], x[2], x[1]);
 480                 n.normalise();
 481                 f += f13*m_SharpEdgesError(0.5*(1 + n_face[1]*n));
 482               }
 483             }
 484           }
 485         }
 486       }
 487     }
 488   }
 489   grid->GetPoints()->SetPoint(nodes[i_nodes_opt], x_old.data());
 490   n_node.normalise();
 491
 492   //CHECK THAT THERE ARE NO BOUNDARY NODE NEIGHBOURS!!!!!!!!!!
 493
 494   if (m_SharpNodesError.active()) {
 495     if (n_pri.size() > 0) {
 496       bool apply_error = true;
 497       for (int i = 0; i < m_Part.n2nLSize(i_nodes_opt); ++i) {
 498         vtkIdType id_neigh_node = m_Part.n2nLG(i_nodes_opt, i);
 499         vtkIdType id_foot = m_IdFoot[nodes[i_nodes_opt]];
 500         int Nbcs = m_Node2BC[id_neigh_node].size();
 501         if (m_Node2BC[id_neigh_node].size() != 0 && id_neigh_node != m_IdFoot[nodes[i_nodes_opt]]) {
 502           apply_error = false;
 503           break;
 504         }
 505       }
 506       if (apply_error) {
 507         double e = 0;
 508         foreach (vec3_t n1, n_pri) {
 509           foreach (vec3_t n2, n_pri) {
 510             //e = max(e, m_SharpNodesError(angleX(n1, n2)));
 511             e = max(e, m_SharpNodesError(0.5*(1 + n1*n2)));
 512           }
 513         }
 514         f += e;
 515       }
 516     }
 517   }
 518
 519   // smooth feature lines
 520   vtkIdType id_node1 = nodes[i_nodes_opt];
 521   vtkIdType id_foot1 = m_IdFoot[id_node1];
 522   if (id_foot1 != -1 && m_FeatureLineError.active()) {
 523     if (m_Node2BC[id_foot1].size() == 2) {
 524       QVector<vtkIdType> ids(2, -1);
 525       int N = 0;
 526       for (int i = 0; i < m_Part.n2nGSize(id_node1); ++i) {
 527         vtkIdType id_node2 = m_Part.n2nGG(id_node1, i);
 528         vtkIdType id_foot2 = m_IdFoot[id_node2];
 529         if (id_foot2 != -1) {
 530           bool found = true;
 531           foreach (int bc, m_Node2BC[id_foot1]) {
 532             if (!m_Node2BC[id_foot2].contains(bc)) {
 533               found = false;
 534               break;
 535             }
 536           }
 537           if (found) {
 538             if (N < 2) {
 539               ids[N] = id_node2;
 540             }
 541             ++N;
 542           }
 543         }
 544       }
 545       if (N == 2) {
 546         vec3_t a, b;
 547         grid->GetPoint(ids[0], a.data());
 548         grid->GetPoint(ids[1], b.data());
 549         vec3_t u = x - a;
 550         vec3_t v = b - x;
 551         f += m_FeatureLineError(0.5*(1 + u*v));
 552       }
 553     }
 554   }
 555
 556   return f;
 557 }
 558
 559 void GridSmoother::resetStencil()
 560 {
 561   stencil.clear();
 562 }
 563
 564 void GridSmoother::addToStencil(double C, vec3_t x)
 565 {
 566   stencil_node_t sn;
 567   sn.x = x;
 568   sn.C = C;
 569   stencil.append(sn);
 570 }
 571
 572 void GridSmoother::operate()
 573 {
 574   markNodes();
 575   computeNormals();
 576
 577   EG_VTKDCC(vtkIntArray, bc,          grid, "cell_code");
 578   EG_VTKDCN(vtkIntArray, node_status, grid, "node_status");
 579   EG_VTKDCN(vtkIntArray, node_layer,  grid, "node_layer");
 580
 581   m_Node2BC.fill(QSet<int>(), grid->GetNumberOfPoints());
 582   for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
 583     if (isSurface(id_cell, grid)) {
 584       vtkIdType N_pts, *pts;
 585       grid->GetCellPoints(id_cell, N_pts, pts);
 586       for (int i_pts = 0; i_pts < N_pts; ++i_pts) {
 587         m_Node2BC[pts[i_pts]].insert(bc->GetValue(id_cell));
 588       }
 589     }
 590   }
 591
 592   m_IdFoot.fill(-1, grid->GetNumberOfPoints());
 593   m_L.fill(0, grid->GetNumberOfPoints());
 594
 595
 596   l2g_t  nodes = getPartNodes();
 597   g2l_t _nodes = getPartLocalNodes();
 598   l2g_t  cells = getPartCells();
 599   l2l_t  n2c   = getPartN2C();
 600   l2l_t  n2n   = getPartN2N();
 601
 602   QVector<bool> prism_node(nodes.size(),false);
 603   foreach (vtkIdType id_cell, cells) {
 604     if (grid->GetCellType(id_cell) == VTK_WEDGE) {
 605       vtkIdType N_pts, *pts;
 606       grid->GetCellPoints(id_cell, N_pts, pts);
 607       for (int i_pts = 0; i_pts < N_pts; ++i_pts) {
 608         if (_nodes[pts[i_pts]] != -1) {
 609           prism_node[_nodes[pts[i_pts]]] = true;
 610         }
 611       }
 612     }
 613   }
 614
 615   F_old = 0;
 616   F_max_old = 0;
 617   for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
 618     if (prism_node[i_nodes]) {
 619       vec3_t x;
 620       i_nodes_opt = i_nodes;
 621       grid->GetPoint(nodes[i_nodes], x.data());
 622       double f = func(x);
 623       F_old += f;
 624       F_max_old = max(F_max_old,f);
 625     }
 626   }
 627
 628   cout << "\nsmoothing volume mesh (" << N_marked_nodes << " nodes)" << endl;
 629   for (int i_iterations = 0; i_iterations < N_iterations; ++i_iterations) {
 630     cout << "iteration " << i_iterations+1 << "/" << N_iterations << endl;
 631     int N_blocked  = 0;
 632     int N_searched = 0;
 633     int N_illegal  = 0;
 634     int N1 = 0;
 635     int N2 = 0;
 636     int N3 = 0;
 637     QTime start = QTime::currentTime();
 638     for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
 639       dbg = false;
 640       bool smooth = node_marked[nodes[i_nodes]];
 641       if (smooth) {
 642         foreach (int bc, m_Node2BC[_nodes[i_nodes]]) {
 643           if (m_BoundaryCodes.contains(bc) || (m_BoundaryCodes.size() ==0)) {
 644             smooth = false;
 645           }
 646         }
 647         foreach (int i_cells, n2c[i_nodes]) {
 648           vtkIdType type_cell = grid->GetCellType(cells[i_cells]);
 649           if ((type_cell == VTK_WEDGE) && !smooth_prisms) {
 650             smooth = false;
 651           }
 652         }
 653       }
 654       if (smooth) {
 655         vtkIdType id_node = nodes[i_nodes];
 656         vec3_t xn;
 657         vec3_t x_old;
 658         grid->GetPoint(id_node, x_old.data());
 659         resetStencil();
 660         bool is_surf = m_Node2BC[_nodes[i_nodes]].size() > 0;
 661         int N = 0;
 662         foreach (int j_nodes, n2n[i_nodes]) {
 663           if (!is_surf || (m_Node2BC[_nodes[j_nodes]].size() > 0)) {
 664             vtkIdType id_neigh_node = nodes[j_nodes];
 665             grid->GetPoint(id_neigh_node, xn.data());
 666             addToStencil(1.0, xn);
 667             ++N;
 668           }
 669         }
 670         if (N == 0) {
 671           EG_BUG;
 672         }
 673         vec3_t x_new1 = vec3_t(0,0,0);
 674         sum_C = 0;
 675         L0 = 0;
 676         double L_min = 1e99;
 677         foreach (stencil_node_t sn, stencil) {
 678           sum_C += sn.C;
 679           x_new1 += sn.C*sn.x;
 680           double L = (sn.x - x_old).abs();
 681           L0 += sn.C*L;
 682           L_min = min(L_min, L);
 683         }
 684         L0 /= sum_C;
 685         x_new1 *= 1.0/sum_C;
 686         vec3_t Dx1 = x_new1 - x_old;
 687         setDeltas(1e-6*L0);
 688         //setDeltas(1e-6);
 689         i_nodes_opt = i_nodes;
 690         vec3_t Dx2(0,0,0);
 691         Dx2 = m_UnderRelaxation*optimise(x_old);
 692         vec3_t Dx3 = (-1e-4/func(x_old))*grad_f;
 693         correctDx(i_nodes, Dx1);
 694         correctDx(i_nodes, Dx2);
 695         correctDx(i_nodes, Dx3);
 696         vec3_t Dx = Dx1;
 697         double f  = func(x_old + Dx1);
 698         int _N1 = 1;
 699         int _N2 = 0;
 700         int _N3 = 0;
 701         if (f > func(x_old + Dx2)) {
 702           Dx = Dx2;
 703           f = func(x_old + Dx2);
 704           _N1 = 0;
 705           _N2 = 1;
 706           _N3 = 0;
 707         }
 708         if (f > func(x_old + Dx3)) {
 709           Dx = Dx3;
 710           _N1 = 0;
 711           _N2 = 0;
 712           _N3 = 1;
 713         }
 714
 715         if (!moveNode(i_nodes, Dx)) {
 716           // search for a better place
 717           vec3_t x_save = x_old;
 718           vec3_t ex(L_search*L0/N_search,0,0);
 719           vec3_t ey(0,L_search*L0/N_search,0);
 720           vec3_t ez(0,0,L_search*L0/N_search);
 721           vec3_t x_best = x_old;
 722           double f_min = func(x_old);
 723           bool found = false;
 724           bool illegal = false;
 725           if (!setNewPosition(id_node,x_old)) {
 726             illegal = true;
 727           }
 728           for (int i = -N_search; i <= N_search; ++i) {
 729             for (int j = -N_search; j <= N_search; ++j) {
 730               for (int k = -N_search; k <= N_search; ++k) {
 731                 if ((i != 0) || (j != 0) || (k != 0)) {
 732                   vec3_t Dx = double(i)*ex + double(j)*ey + double(k)*ez;
 733                   correctDx(i_nodes, Dx);
 734                   vec3_t x = x_old + x;
 735                   if (setNewPosition(id_node, x)) {
 736                     grid->GetPoints()->SetPoint(id_node, x_old.data());
 737                     double f = func(x);
 738                     if (f < f_min) {
 739                       f_min = f;
 740                       x_best = x;
 741                       found = true;
 742                       illegal = false;
 743                     }
 744                   }
 745                 }
 746               }
 747             }
 748           }
 749           if (found) {
 750             grid->GetPoints()->SetPoint(id_node, x_best.data());
 751             ++N_searched;
 752           } else {
 753             ++N_blocked;
 754             if (illegal) {
 755               ++N_illegal;
 756             }
 757           }
 758         } else {
 759           N1 += _N1;
 760           N2 += _N2;
 761           N3 += _N3;
 762         }
 763       }
 764     }
 765
 766     cout << N1 << " type 1 movements (simple)" << endl;
 767     cout << N2 << " type 2 movements (Newton)" << endl;
 768     cout << N3 << " type 3 movements (gradient)" << endl;
 769     cout << N_searched << " type X movements (search)" << endl;
 770     cout << N_blocked << " type 0 movements (failure)" << endl;
 771
 772     cout << start.secsTo(QTime::currentTime()) << " seconds elapsed" << endl;
 773     F_new = 0;
 774     F_max_new = 0;
 775
 776     m_HeightError.reset();
 777     m_EdgeLengthError.reset();
 778     m_TetraError.reset();
 779     m_SharpEdgesError.reset();
 780     m_SharpNodesError.reset();
 781     m_ParallelEdgesError.reset();
 782     m_ParallelFacesError.reset();
 783     m_SimilarFaceAreaError.reset();
 784     m_FeatureLineError.reset();
 785     for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
 786       if (prism_node[i_nodes]) {
 787         vec3_t x;
 788         i_nodes_opt = i_nodes;
 789         grid->GetPoint(nodes[i_nodes], x.data());
 790         double f = func(x);
 791         F_new += f;
 792         F_max_new = max(F_max_new,f);
 793       }
 794     }
 795     cout << "total error (old) = " << F_old << endl;
 796     cout << "total error (new) = " << F_new << endl;
 797     double f_old = max(1e-10,F_old);
 798     cout << "total improvement = " << 100*(1-F_new/f_old) << "%" << endl;
 799     printMaxErrors();
 800   }
 801   if (m_WriteDebugFile) {
 802     writeDebugFile("gridsmoother");
 803   }
 804   cout << "done" << endl;
 805 }
 806
 807 double GridSmoother::improvement()
 808 {
 809   double f_max_old = max(1e-10,F_max_old);
 810   double f_old = max(1e-10,F_old);
 811   double i2 = 1-F_new/f_old;
 812   return i2;
 813 }
 814
 815 void GridSmoother::printMaxErrors()
 816 {
 817   cout << "maximal height error = " << m_HeightError.maxError() << endl;
 818   cout << "maximal edge lenth error = " << m_EdgeLengthError.maxError() << endl;
 819   cout << "maximal tetra error = " << m_TetraError.maxError() << endl;
 820   cout << "maximal sharp nodes error = " << m_SharpNodesError.maxError() << endl;
 821   cout << "maximal sharp edges error = " << m_SharpEdgesError.maxError() << endl;
 822   cout << "maximal parallel edges error = " << m_ParallelEdgesError.maxError() << endl;
 823   cout << "maximal parallel faces error = " << m_ParallelFacesError.maxError() << endl;
 824   cout << "maximal face area error = " << m_SimilarFaceAreaError.maxError() << endl;
 825   cout << "maximal feature line error = " << m_FeatureLineError.maxError() << endl;
 826 }
 827
 828 void GridSmoother::writeDebugFile(QString file_name)
 829 {
 830   QVector<vtkIdType> bcells;
 831   getSurfaceCells(m_BoundaryCodes, bcells, grid);
 832   MeshPartition bpart(grid);
 833   bpart.setCells(bcells);
 834   QVector<vtkIdType> foot2field(bpart.getNumberOfNodes(), -1);
 835   for (vtkIdType id_node = 0; id_node < grid->GetNumberOfPoints(); ++id_node) {
 836     vtkIdType id_foot = m_IdFoot[id_node];
 837     if (id_foot != -1) {
 838       if (bpart.localNode(id_foot) == -1) {
 839         EG_BUG;
 840       }
 841       foot2field[bpart.localNode(id_foot)] = id_node;
 842     }
 843   }
 844   file_name = GuiMainWindow::pointer()->getCwd() + "/" + file_name + ".vtk";
 845   QFile file(file_name);
 846   file.open(QIODevice::WriteOnly);
 847   QTextStream f(&file);
 848   f << "# vtk DataFile Version 2.0\n";
 849   f << "m_NodeNormal\n";
 850   f << "ASCII\n";
 851   f << "DATASET UNSTRUCTURED_GRID\n";
 852   f << "POINTS " << 2*bpart.getNumberOfNodes() << " float\n";
 853   for (int i = 0; i < bpart.getNumberOfNodes(); ++i) {
 854     vec3_t x;
 855     vtkIdType id_node = bpart.globalNode(i);
 856     grid->GetPoint(id_node, x.data());
 857     f << x[0] << " " << x[1] << " " << x[2] << "\n";
 858   }
 859   for (int i = 0; i < bpart.getNumberOfNodes(); ++i) {
 860     vec3_t x;
 861     vtkIdType id_node = foot2field[i];
 862     grid->GetPoint(id_node, x.data());
 863     f << x[0] << " " << x[1] << " " << x[2] << "\n";
 864   }
 865   f << "CELLS " << 2*bpart.getNumberOfCells() << " " << 8*bpart.getNumberOfCells() << "\n";
 866   for (int i = 0; i < bpart.getNumberOfCells(); ++ i) {
 867     vtkIdType id_cell = bpart.globalCell(i);
 868     vtkIdType N_pts, *pts;
 869     if (grid->GetCellType(id_cell) != VTK_TRIANGLE) {
 870       EG_BUG;
 871     }
 872     grid->GetCellPoints(id_cell, N_pts, pts);
 873     QVector<int> nds1(3), nds2(3);
 874     for (int j = 0; j < 3; ++j) {
 875       nds1[j] = bpart.localNode(pts[j]);
 876       nds2[j] = nds1[j] + bpart.getNumberOfNodes();
 877     }
 878     f << "3 " << nds1[0] << " " << nds1[1] << " " << nds1[2] << "\n";
 879     f << "3 " << nds2[0] << " " << nds2[1] << " " << nds2[2] << "\n";
 880   }
 881
 882   f << "CELL_TYPES " << 2*bpart.getNumberOfCells() << "\n";
 883   for (int i = 0; i < 2*bpart.getNumberOfCells(); ++ i) {
 884     f << VTK_TRIANGLE << "\n";
 885   }
 886   f << "POINT_DATA " << 2*bpart.getNumberOfNodes() << "\n";
 887   f << "VECTORS N float\n";
 888
 889   for (int i = 0; i < bpart.getNumberOfNodes(); ++i) {
 890     vtkIdType id_node = bpart.globalNode(i);
 891     f << m_NodeNormal[id_node][0] << " " << m_NodeNormal[id_node][1] << " " << m_NodeNormal[id_node][2] << "\n";
 892   }
 893   for (int i = 0; i < bpart.getNumberOfNodes(); ++i) {
 894     f << "0 0 0\n";
 895   }
 896 }
 897