src/libengrid/guidivideboundarylayer.cpp

   1 //
   2 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   3 // +                                                                      +
   4 // + This file is part of enGrid.                                         +
   5 // +                                                                      +
   6 // + Copyright 2008-2013 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 "guidivideboundarylayer.h"
  24 #include "math/linsolve.h"
  25
  26 #include "volumedefinition.h"
  27 #include "guimainwindow.h"
  28
  29 void GuiDivideBoundaryLayer::before()
  30 {
  31   populateBoundaryCodes(m_Ui.listWidgetBC);
  32   populateVolumes(m_Ui.listWidgetVC);
  33
  34   getSet("boundary layer", "first critical angle",  180, m_CritAngle1);
  35   getSet("boundary layer", "second critical angle", 270, m_CritAngle2);
  36
  37   QString blayer_txt = GuiMainWindow::pointer()->getXmlSection("blayer/global");
  38   QTextStream s(&blayer_txt);
  39   double v;
  40   if (!s.atEnd()) {
  41     s >> v;
  42     QString num;
  43     num.setNum(v);
  44     m_Ui.lineEditAbsolute->setText(num);
  45   }
  46   if (!s.atEnd()) {
  47     s >> v; // relative height
  48     m_Ui.doubleSpinBoxHeight->setValue(v);
  49   }
  50   if (!s.atEnd()) {
  51     s >> v; // blending
  52     m_Ui.doubleSpinBoxBlending->setValue(v);
  53   }
  54   if (!s.atEnd()) {
  55     s >> v;
  56     m_Ui.doubleSpinBoxStretching->setValue(v);
  57   }
  58   if (!s.atEnd()) {
  59     s >> v;
  60     m_Ui.doubleSpinBoxFarRatio->setValue(v);
  61   }
  62   if (!s.atEnd()) {
  63     int v;
  64     s >> v;
  65     m_Ui.spinBoxLayers->setValue(v);
  66   }
  67
  68   m_Ui.doubleSpinBoxCritAngle1->setValue(m_CritAngle1);
  69   m_Ui.doubleSpinBoxCritAngle2->setValue(m_CritAngle2);
  70
  71   m_RestGrid = vtkUnstructuredGrid::New();
  72 }
  73
  74 bool GuiDivideBoundaryLayer::findBoundaryLayer()
  75 {
  76   EG_VTKDCC(vtkIntArray, bc, m_Grid, "cell_code");
  77
  78   l2g_t cells = getPartCells();
  79   l2l_t c2c   = getPartC2C();
  80
  81   m_BoundaryCodes.clear();
  82   m_Pairs.clear();
  83   m_InsertCell.fill(true,m_Grid->GetNumberOfCells());
  84   m_NumPrisms = 0;
  85   m_NumQuads = 0;
  86   for (int i_cells = 0; i_cells < cells.size(); ++i_cells) {
  87     if (m_Grid->GetCellType(cells[i_cells]) == VTK_WEDGE) {
  88       ++m_NumPrisms;
  89       vtkIdType N_pts, *pts;
  90       m_Grid->GetCellPoints(cells[i_cells],N_pts,pts);
  91       for (int j = 0; j < 3; ++j) {
  92         m_Pairs.insert(QPair<vtkIdType,vtkIdType>(pts[j],pts[j+3]));
  93       }
  94       for (int j = 2; j < 5; ++j) {
  95         if (c2c[i_cells][j] != -1) {
  96           vtkIdType type_ncell = m_Grid->GetCellType(cells[c2c[i_cells][j]]);
  97           if ((type_ncell != VTK_WEDGE) && (type_ncell != VTK_QUAD)) {
  98             finalise();
  99             EG_ERR_RETURN("unable to identify boundary layer");
 100             return(false);
 101           }
 102         } else {
 103           vec3_t x;
 104           m_Grid->GetPoint(pts[0],x.data());
 105           cout << x << endl;
 106           EG_BUG;
 107         }
 108       }
 109       {
 110         vtkIdType id_tri = m_Part.c2cLG(i_cells, 0);
 111         if (m_Grid->GetCellType(id_tri) != VTK_TRIANGLE) {
 112           EG_BUG;
 113         }
 114         m_BoundaryCodes.insert(bc->GetValue(id_tri));
 115       }
 116       for (int j = 0; j < 2; ++j) {
 117         if (c2c[i_cells][j] != -1) {
 118           vtkIdType type_ncell = m_Grid->GetCellType(cells[c2c[i_cells][j]]);
 119           if (type_ncell == VTK_WEDGE) {
 120             finalise();
 121             EG_ERR_RETURN("the boundary layer seems to have been split already");
 122             return(false);
 123           }
 124         } else {
 125           EG_BUG;
 126         }
 127       }
 128     }
 129     if (m_Grid->GetCellType(cells[i_cells]) == VTK_QUAD) {
 130       ++m_NumQuads;
 131     }
 132   }
 133   if (m_NumPrisms == 0) {
 134     finalise();
 135     EG_ERR_RETURN("unable to identify boundary layer");
 136     return(false);
 137   }
 138   m_IsBlayerNode.clear();
 139   m_IsBlayerNode.fill(false, m_Grid->GetNumberOfPoints());
 140
 141   QPair<vtkIdType,vtkIdType> P;
 142   foreach (P, m_Pairs) {
 143     m_IsBlayerNode[P.second] = true;
 144   }
 145
 146   computeMaxConvexAngles();
 147   return(true);
 148 }
 149
 150 void GuiDivideBoundaryLayer::computeY1()
 151 {
 152   /*
 153   double s1 = 0.01;
 154   double s2 = 10*m_DesiredStretching;
 155   while (fabs(s1-s2) > 1e-4) {
 156     double s = 0.5*(s1+s2);
 157     for (int i = 2; i < m_Y.size(); ++i) {
 158       m_Y[i] = m_Y[i-1] + s*(m_Y[i-1]-m_Y[i-2]);
 159     }
 160     if (m_Y.last() < 1) {
 161       s1 = s;
 162     } else {
 163       s2 = s;
 164     }
 165   }
 166   */
 167   double C1 = 0.1;
 168   double C2 = 2.0;
 169   while (C2 - C1 > 1e-6) {
 170     double s = m_DesiredStretching;
 171     for (int i = 2; i <= m_NumLayers; ++i) {
 172       m_Y[i] = m_Y[i-1] + s*(m_Y[i-1] - m_Y[i-2]);
 173       s *= 0.5*(C1 + C2);
 174     }
 175     if (m_Y[m_NumLayers] > 1) {
 176       C2 = 0.5*(C1 + C2);
 177     } else {
 178       C1 = 0.5*(C1 + C2);
 179     }
 180   }
 181   m_Y.last() = 1;
 182 }
 183
 184 void GuiDivideBoundaryLayer::computeY2()
 185 {
 186   double C1 = 1.0;
 187   double C2 = 100.0;
 188   double y_target = m_Y[m_NumLayers - 1];
 189   while (C2 - C1 > 1e-6) {
 190     double s = 0.5*(C1 + C2);
 191     for (int i = 2; i < m_NumLayers; ++i) {
 192       m_Y[i] = m_Y[i-1] + s*(m_Y[i-1] - m_Y[i-2]);
 193     }
 194     if (m_Y[m_NumLayers - 1] > y_target) {
 195       C2 = 0.5*(C1 + C2);
 196     } else {
 197       C1 = 0.5*(C1 + C2);
 198     }
 199   }
 200   //m_Y[m_NumLayers - 1] = y_target;
 201 }
 202
 203 void GuiDivideBoundaryLayer::createEdges(vtkUnstructuredGrid *new_grid)
 204 {
 205   EG_VTKDCN(vtkDoubleArray, cl, m_Grid, "node_meshdensity_desired" );
 206   m_Edges.fill(QVector<vtkIdType>(m_NumLayers+1), m_Pairs.size());
 207   m_Old2Edge.fill(-1, m_Grid->GetNumberOfPoints());
 208   int N = 0;
 209   vtkIdType id_new_node = m_Grid->GetNumberOfPoints();
 210   QPair<vtkIdType,vtkIdType> P;
 211   double max_step = 0;
 212   double ymax     = 0;
 213   double ymin     = 1e99;
 214   foreach (P, m_Pairs) {
 215     m_Edges[N][0] = P.first;
 216     m_Edges[N][m_NumLayers] = P.second;
 217     m_Old2Edge[P.first] = N;
 218     m_Old2Edge[P.second] = N;
 219
 220     vec3_t x1,x2;
 221     m_Grid->GetPoint(P.first, x1.data());
 222     m_Grid->GetPoint(P.second, x2.data());
 223     vec3_t n = x2-x1;
 224     double alpha  = GeometryTools::rad2deg(m_MaxConvexAngle[P.first]);
 225     double cl_loc = cl->GetValue(P.first);
 226     if (cl_loc < 1e-99) {
 227       cl_loc = m_Part.getAverageSurfaceEdgeLength(P.first);
 228     }
 229     double h_rel = m_RelativeHeight*cl_loc/n.abs();
 230     {
 231       m_Y.resize(m_NumLayers + 1);
 232       m_Y[0] = 0;
 233       m_Y[1] = m_Blending*m_AbsoluteHeight/n.abs() + (1-m_Blending)*h_rel;
 234       computeY1();
 235     }
 236     if (alpha > m_CritAngle1) {
 237       double blend = min(1.0, (alpha - m_CritAngle1)/(m_CritAngle2 - m_CritAngle1));
 238       double far_ratio = blend*m_FarRatio + (1-blend)*(1.0 - m_Y[m_NumLayers - 1])*n.abs()/cl_loc;
 239       m_Y.resize(m_NumLayers + 1);
 240       m_Y[0] = 0;
 241       m_Y[1] = m_Blending*m_AbsoluteHeight/n.abs() + (1-m_Blending)*h_rel;
 242       m_Y[m_NumLayers - 1] = 1.0 - far_ratio*cl_loc/n.abs();
 243       m_Y[m_NumLayers] = 1.0;
 244       computeY2();
 245     }
 246     ymin = min(ymin, m_Y[1]*n.abs());
 247     ymax = max(ymax, m_Y[1]*n.abs());
 248     for (int i = 1; i < m_NumLayers; ++i) {
 249       vec3_t x = x1 + m_Y[i]*n;
 250       max_step = max(max_step,(m_Y[i+1]-m_Y[i])/(m_Y[i]-m_Y[i-1]));
 251       new_grid->GetPoints()->SetPoint(id_new_node, x.data());
 252       m_Edges[N][i] = id_new_node;
 253       ++id_new_node;
 254     }
 255     ++N;
 256   }
 257   cout << LINE;
 258   cout << "maximal increment : " << max_step << endl;
 259   cout << "min(y) : " << ymin << endl;
 260   cout << "max(y) : " << ymax << endl;
 261   cout << LINE;
 262 }
 263
 264 void GuiDivideBoundaryLayer::computeMaxConvexAngles()
 265 {
 266   EG_VTKDCC(vtkIntArray, bc, m_Grid, "cell_code");
 267   m_MaxConvexAngle.fill(0, m_Grid->GetNumberOfPoints());
 268   EG_FORALL_CELLS (id_cell1, m_Grid) {
 269     if (isSurface(id_cell1, m_Grid)) {
 270       if (m_BoundaryCodes.contains(bc->GetValue(id_cell1))) {
 271         EG_GET_CELL (id_cell1, m_Grid);
 272         QVector<vec3_t> x(num_pts + 1);
 273         for (int i = 0; i < num_pts; ++i) {
 274           m_Grid->GetPoint(pts[i], x[i].data());
 275         }
 276         x.last() = x.first();
 277         vec3_t n1 = cellNormal(m_Grid, id_cell1);
 278         n1.normalise();
 279         vec3_t x1 = cellCentre(m_Grid, id_cell1);
 280         for (int i = 0; i < num_pts; ++i) {
 281           vec3_t xe = 0.5*(x[i] + x[i+1]);
 282           vtkIdType id_cell2 = m_Part.c2cGG(id_cell1, i);
 283           if (m_BoundaryCodes.contains(bc->GetValue(id_cell2))) {
 284             vec3_t n2 = cellNormal(m_Grid, id_cell2);
 285             n2.normalise();
 286             double alpha = angle(n1, n2);
 287             vec3_t v1 = x1 - xe;
 288             vec3_t v2 = cellCentre(m_Grid, id_cell2) - xe;
 289             vec3_t v = v1 + v2;
 290             v.normalise();
 291             double sp = v*(n1+n2);
 292             if (sp > 0) {
 293               alpha = M_PI + alpha;
 294             } else {
 295               alpha = M_PI - alpha;
 296             }
 297             vtkIdType p1 = pts[i];
 298             vtkIdType p2 = pts[0];
 299             if (i < num_pts - 1) {
 300               p2 = pts[i+1];
 301             }
 302             m_MaxConvexAngle[p1] = max(m_MaxConvexAngle[p1], alpha);
 303             m_MaxConvexAngle[p2] = max(m_MaxConvexAngle[p2], alpha);
 304           }
 305         }
 306       }
 307     }
 308   }
 309 }
 310
 311 void GuiDivideBoundaryLayer::operate()
 312 {
 313   // set m_Grid to selected volume
 314   getSelectedItems(m_Ui.listWidgetBC, m_BoundaryCodes); // fill m_BoundaryCodes with values from listWidgetBC
 315   QString volume_name = getSelectedVolume(m_Ui.listWidgetVC);
 316   VolumeDefinition V = GuiMainWindow::pointer()->getVol(volume_name);
 317   foreach (int bc, m_BoundaryCodes) {
 318     qDebug()<<"V.getSign("<<bc<<")="<<V.getSign(bc);
 319     if (V.getSign(bc) == 0) {
 320       QString msg;
 321       msg.setNum(bc);
 322       msg = "Boundary code " + msg + " is not part of the volume '" + volume_name +"'.";
 323       EG_ERR_RETURN(msg);
 324     }
 325   }
 326
 327   EG_VTKSP(vtkUnstructuredGrid, rest_grid);
 328   {
 329     EG_VTKSP(vtkUnstructuredGrid, vol_grid);
 330     MeshPartition volume(volume_name);
 331     MeshPartition rest(m_Grid);
 332     rest.setRemainder(volume);
 333     volume.setVolumeOrientation();
 334     volume.extractToVtkGrid(vol_grid);
 335     rest.extractToVtkGrid(rest_grid);
 336     makeCopy(vol_grid, m_Grid);
 337   }
 338   setAllCells();
 339
 340   m_NumLayers         = m_Ui.spinBoxLayers->value();
 341   m_RelativeHeight    = m_Ui.doubleSpinBoxHeight->value();
 342   m_AbsoluteHeight    = m_Ui.lineEditAbsolute->text().toDouble();
 343   m_Blending          = m_Ui.doubleSpinBoxBlending->value();
 344   m_DesiredStretching = m_Ui.doubleSpinBoxStretching->value();
 345   m_FarRatio          = m_Ui.doubleSpinBoxFarRatio->value();
 346   m_CritAngle1        = m_Ui.doubleSpinBoxCritAngle1->value();
 347   m_CritAngle2        = m_Ui.doubleSpinBoxCritAngle2->value();
 348   cout << "dividing boundary layer into " << m_NumLayers << " layers:" << endl;
 349   if(findBoundaryLayer()) {
 350     EG_VTKSP(vtkUnstructuredGrid,new_grid);
 351     allocateGrid(new_grid, m_Grid->GetNumberOfCells() + (m_NumPrisms + m_NumQuads)*(m_NumLayers-1), m_Grid->GetNumberOfPoints() + m_Pairs.size()*(m_NumLayers-1));
 352
 353
 354     EG_VTKDCC(vtkIntArray, old_orgdir, m_Grid, "cell_orgdir");
 355     EG_VTKDCC(vtkIntArray, old_voldir, m_Grid, "cell_voldir");
 356     EG_VTKDCC(vtkIntArray, old_curdir, m_Grid, "cell_curdir");
 357     EG_VTKDCC(vtkIntArray, new_orgdir, new_grid, "cell_orgdir");
 358     EG_VTKDCC(vtkIntArray, new_voldir, new_grid, "cell_voldir");
 359     EG_VTKDCC(vtkIntArray, new_curdir, new_grid, "cell_curdir");
 360
 361     int orgdir = -99;
 362     int curdir = -99;
 363     int voldir = -99;
 364
 365     // copy existing mesh without prisms and adjacent cells
 366     vtkIdType id_new_node = 0;
 367     for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 368       vec3_t x;
 369       m_Grid->GetPoint(id_node, x.data());
 370       new_grid->GetPoints()->SetPoint(id_new_node, x.data());
 371       copyNodeData(m_Grid, id_node, new_grid, id_new_node);
 372       ++id_new_node;
 373     }
 374     vtkIdType id_new_cell;
 375     for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 376       vtkIdType N_pts, *pts;
 377       m_Grid->GetCellPoints(id_cell, N_pts, pts);
 378       bool insert_cell = true;
 379       if (m_Grid->GetCellType(id_cell) == VTK_WEDGE) {
 380         insert_cell = false;
 381       }
 382       if (m_Grid->GetCellType(id_cell) == VTK_QUAD) {
 383         insert_cell = false;
 384         if (orgdir != -99 && old_orgdir->GetValue(id_cell) != orgdir) {
 385           EG_BUG;
 386         }
 387         if (voldir != -99 && old_voldir->GetValue(id_cell) != voldir) {
 388           EG_BUG;
 389         }
 390         if (curdir != -99 && old_curdir->GetValue(id_cell) != curdir) {
 391           EG_BUG;
 392         }
 393         orgdir = old_orgdir->GetValue(id_cell);
 394         voldir = old_voldir->GetValue(id_cell);
 395         curdir = old_curdir->GetValue(id_cell);
 396       }
 397       if (insert_cell) {
 398         id_new_cell = new_grid->InsertNextCell(m_Grid->GetCellType(id_cell), N_pts, pts);
 399         copyCellData(m_Grid, id_cell, new_grid, id_new_cell);
 400       }
 401     }
 402
 403     // create divided boundary layer
 404     createEdges(new_grid);
 405
 406     for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 407       if (m_Grid->GetCellType(id_cell) == VTK_WEDGE) {
 408         vtkIdType N_pts, *pts;
 409         m_Grid->GetCellPoints(id_cell, N_pts, pts);
 410         for (int i = 0; i < m_NumLayers; ++i) {
 411           vtkIdType p[6];
 412           p[0] = m_Edges[m_Old2Edge[pts[0]]][i];
 413           p[1] = m_Edges[m_Old2Edge[pts[1]]][i];
 414           p[2] = m_Edges[m_Old2Edge[pts[2]]][i];
 415           p[3] = m_Edges[m_Old2Edge[pts[0]]][i+1];
 416           p[4] = m_Edges[m_Old2Edge[pts[1]]][i+1];
 417           p[5] = m_Edges[m_Old2Edge[pts[2]]][i+1];
 418           id_new_cell = new_grid->InsertNextCell(VTK_WEDGE, 6, p);
 419           copyCellData(m_Grid, id_cell, new_grid, id_new_cell);
 420         }
 421       }
 422       if (m_Grid->GetCellType(id_cell) == VTK_QUAD) {
 423         vtkIdType N_pts, *pts;
 424         m_Grid->GetCellPoints(id_cell, N_pts, pts);
 425         if ((m_Old2Edge[pts[0]] != -1) && (m_Old2Edge[pts[1]] != -1) && (m_Old2Edge[pts[2]] != -1) && (m_Old2Edge[pts[3]] != -1)) {
 426           for (int i = 0; i < m_NumLayers; ++i) {
 427             vtkIdType p[4];
 428             p[0] = m_Edges[m_Old2Edge[pts[0]]][i];
 429             p[1] = m_Edges[m_Old2Edge[pts[1]]][i];
 430             p[2] = m_Edges[m_Old2Edge[pts[1]]][i+1];
 431             p[3] = m_Edges[m_Old2Edge[pts[0]]][i+1];
 432             id_new_cell = new_grid->InsertNextCell(VTK_QUAD, 4, p);
 433             copyCellData(m_Grid, id_cell, new_grid, id_new_cell);
 434             new_orgdir->SetValue(id_new_cell, orgdir);
 435             new_voldir->SetValue(id_new_cell, voldir);
 436             new_curdir->SetValue(id_new_cell, curdir);
 437           }
 438         }
 439       }
 440     }
 441
 442     makeCopy(new_grid, m_Grid);
 443
 444     ///////////////////////////////////////////////////////////////
 445     // set m_Grid to modified selected volume + unselected volumes
 446     {
 447       MeshPartition volume(m_Grid, true);
 448       MeshPartition rest(rest_grid, true);
 449       volume.addPartition(rest);
 450     }
 451     resetOrientation(m_Grid);
 452     createIndices(m_Grid);
 453     ///////////////////////////////////////////////////////////////
 454
 455   }
 456 }
 457
 458 void GuiDivideBoundaryLayer::finalise()
 459 {
 460   // set m_Grid to modified selected volume + unselected volumes
 461   {
 462     MeshPartition volume(m_Grid, true);
 463     MeshPartition rest(m_RestGrid, true);
 464     volume.addPartition(rest);
 465   }
 466   resetOrientation(m_Grid);
 467   createIndices(m_Grid);
 468   m_RestGrid->Delete();
 469 }
 470