src/libengrid/guidivideboundarylayer.cpp

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