src/libengrid/createboundarylayershell.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
  22 #include "createboundarylayershell.h"
  23
  24 #include "createvolumemesh.h"
  25 #include "swaptriangles.h"
  26 #include "deletetetras.h"
  27 #include "deletecells.h"
  28 #include "meshpartition.h"
  29 #include "deletevolumegrid.h"
  30 #include "laplacesmoother.h"
  31 #include "surfacemeshsmoother.h"
  32 #include "deletecells.h"
  33 #include "stitchholes.h"
  34
  35 CreateBoundaryLayerShell::CreateBoundaryLayerShell()
  36 {
  37   m_RestGrid      = vtkSmartPointer<vtkUnstructuredGrid>::New();
  38   m_OriginalGrid  = vtkSmartPointer<vtkUnstructuredGrid>::New();
  39   m_PrismaticGrid = vtkSmartPointer<vtkUnstructuredGrid>::New();
  40   m_Success = false;
  41 }
  42
  43 void CreateBoundaryLayerShell::prepare()
  44 {
  45   m_Part.trackGrid(m_Grid);
  46
  47   DeleteVolumeGrid delete_volume;
  48   delete_volume.setGrid(m_Grid);
  49   delete_volume.setAllCells();
  50   delete_volume();
  51
  52   readSettings();
  53   setAllCells();
  54   getSurfaceCells(m_BoundaryLayerCodes, layer_cells, m_Grid);
  55
  56   // fill m_LayerAdjacentBoundaryCodes
  57   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
  58   foreach (vtkIdType id_cell, layer_cells) {
  59     for (int i = 0; i < m_Part.c2cGSize(id_cell); ++i) {
  60       vtkIdType id_neigh = m_Part.c2cGG(id_cell, i);
  61       int bc = cell_code->GetValue(id_neigh);
  62       if (!m_BoundaryLayerCodes.contains(bc)) {
  63         m_LayerAdjacentBoundaryCodes.insert(bc);
  64       }
  65     }
  66   }
  67
  68   // compute normals and origins of adjacent planes
  69   m_LayerAdjacentNormals.clear();
  70   m_LayerAdjacentOrigins.clear();
  71   foreach (int bc, m_LayerAdjacentBoundaryCodes) {
  72     double L = EG_LARGE_REAL;
  73     vec3_t n0(0, 0, 0);
  74     vec3_t x0(0, 0, 0);
  75     double total_area = 0;
  76     for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
  77       if (isSurface(id_cell, m_Grid) && cell_code->GetValue(id_cell) == bc) {
  78         vec3_t n = cellNormal(m_Grid, id_cell);
  79         double A = n.abs();
  80         total_area += A;
  81         n0 += n;
  82         x0 += A*cellCentre(m_Grid, id_cell);
  83         L = min(L, sqrt(4*A/sqrt(3.0)));
  84       }
  85     }
  86     n0.normalise();
  87     x0 *= 1.0/total_area;
  88     m_LayerAdjacentNormals[bc] = n0;
  89     m_LayerAdjacentOrigins[bc] = x0;
  90   }
  91
  92   computeBoundaryLayerVectors();
  93   makeCopy(m_Grid, m_OriginalGrid);
  94 }
  95
  96 QList<vtkIdType> CreateBoundaryLayerShell::correctAdjacentBC(int bc, int num_levels)
  97 {
  98   cout << "correcting boundary \"" << qPrintable(GuiMainWindow::pointer()->getBC(bc).getName()) << "\"" << endl;
  99
 100   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
 101   EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type");
 102   double scal_min = -1;
 103   int    count    = 0;
 104
 105   SurfaceMeshSmoother smooth;
 106   smooth.setGrid(m_Grid);
 107   smooth.useSimpleCentreScheme();
 108
 109   // mark num_levels levels of nodes next to the boundary layer
 110   QVector<bool> marked(m_Grid->GetNumberOfPoints(), false);
 111   QList<vtkIdType> marked_nodes;
 112   EG_FORALL_NODES(id_node, m_Grid) {
 113     if (m_BoundaryLayerNode[id_node]) {
 114       marked[id_node] = true;
 115       marked_nodes << id_node;
 116     }
 117   }
 118   for (int level = 0; level < num_levels; ++level) {
 119     QList<vtkIdType> new_marked_nodes;
 120     foreach (vtkIdType id_node, marked_nodes) {
 121       for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
 122         vtkIdType id_neigh = m_Part.n2nGG(id_node, i);
 123         if (!marked[id_neigh]) {
 124           marked[id_neigh] = true;
 125           new_marked_nodes << id_neigh;
 126         }
 127       }
 128     }
 129     marked_nodes = new_marked_nodes;
 130   }
 131
 132   QList<vtkIdType> cells;
 133   EG_FORALL_CELLS(id_cell, m_Grid) {
 134     if (isSurface(id_cell, m_Grid)) {
 135       if (cell_code->GetValue(id_cell) == bc) {
 136         cells << id_cell;
 137       }
 138     }
 139   }
 140
 141   CadInterface *cad = GuiMainWindow::pointer()->getCadInterface(bc);
 142   smooth.setCells(cells);
 143   smooth.prepareCadInterface(cad);
 144
 145   QList<vtkIdType> bad_nodes;
 146   int new_num_bad = m_Grid->GetNumberOfPoints();
 147   int old_num_bad = 0;
 148   updateNodeInfo();
 149   do {
 150     old_num_bad = new_num_bad;
 151     scal_min = 1;
 152     bad_nodes.clear();
 153     for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 154       if (!m_BoundaryLayerNode[id_node]){
 155         bool node_has_bc = false;
 156         bool all_bcs_adjacent = true;
 157         for (int i = 0; i < m_Part.n2bcGSize(id_node); ++i) {
 158           int bcn = m_Part.n2bcG(id_node, i);
 159           if (bcn == bc) {
 160             node_has_bc = true;
 161           }
 162           if (!m_LayerAdjacentBoundaryCodes.contains(bcn)) {
 163             all_bcs_adjacent = false;
 164           }
 165         }
 166         if (node_has_bc && all_bcs_adjacent && marked[id_node]) {
 167           if (m_Part.n2bcGSize(id_node) == 1) {
 168             vec3_t xs = smooth.smoothNode(id_node);
 169             xs = cad->snapNode(id_node, xs);
 170             if (!checkVector(xs)) {
 171               EG_ERR_RETURN("error while correcting adjacent boundaries");
 172             }
 173             m_Grid->GetPoints()->SetPoint(id_node, xs.data());
 174           } else if (m_Part.n2bcGSize(id_node) == 2 && node_type->GetValue(id_node) != EG_FIXED_VERTEX) {
 175             int count = 0;
 176             vec3_t xs(0, 0, 0);
 177             for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
 178               vtkIdType id_neigh = m_Part.n2nGG(id_node, i);
 179               if (m_Part.n2bcGSize(id_neigh) > 1) {
 180                 vec3_t x;
 181                 m_Grid->GetPoint(id_neigh, x.data());
 182                 xs += x;
 183                 ++count;
 184               }
 185             }
 186             if (count < 2) {
 187               vec3_t x;
 188               m_Grid->GetPoint(id_node, x.data());
 189               xs += x;
 190               ++count;
 191             }
 192             xs *= 1.0/count;
 193             xs = cad->snapToEdge(xs);
 194             m_Grid->GetPoints()->SetPoint(id_node, xs.data());
 195           }
 196
 197           bool node_bad = false;
 198           for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
 199             vtkIdType id_cell = m_Part.n2cGG(id_node, i);
 200             if (isSurface(id_cell, m_Grid)) {
 201               if (cell_code->GetValue(id_cell) == bc) {
 202                 CadInterface *cad = GuiMainWindow::pointer()->getCadInterface(cell_code->GetValue(id_cell));
 203                 cad->snap(cellCentre(m_Grid, id_cell));
 204                 vec3_t n = cellNormal(m_Grid, id_cell);
 205                 n.normalise();
 206                 double scal = n*cad->getLastNormal();
 207                 scal_min = min(scal_min, scal);
 208                 if (scal < 0.5 && !node_bad) {
 209                   bad_nodes << id_node;
 210                   node_bad = true;
 211                 }
 212               }
 213             }
 214           }
 215         }
 216       }
 217     }
 218     new_num_bad = bad_nodes.size();
 219     //reduceSurface();
 220     //cout << "  " << bad_nodes.size() << " node defects" << endl;
 221     ++count;
 222   } while (scal_min < 0.5 && count < 1000);
 223   if (scal_min < 0.5) {
 224     //writeGrid(grid, "adjacent_bc_failure");
 225     //return false;
 226     //EG_ERR_RETURN("failed to correct adjacent surfaces");
 227   }
 228   return bad_nodes;
 229 }
 230
 231 void CreateBoundaryLayerShell::createLayerNodes(vtkIdType id_node)
 232 {
 233   vec3_t x1;
 234   m_Grid->GetPoint(id_node, x1.data());
 235   m_PrismaticGrid->GetPoints()->SetPoint(m_ShellNodeMap[id_node], x1.data());
 236   vec3_t x2 = x1 + m_BoundaryLayerVectors[id_node];
 237
 238   double H  = m_BoundaryLayerVectors[id_node].abs();
 239   double h  = H*(1.0 - m_StretchingRatio)/(1.0 - pow(m_StretchingRatio, m_NumLayers));
 240   vec3_t dx = (1.0/H)*m_BoundaryLayerVectors[id_node];
 241   vec3_t x  = x1;
 242   m_PrismaticGrid->GetPoints()->SetPoint(m_ShellNodeMap[id_node], x1.data());
 243   for (int i = 1; i < m_NumLayers; ++i) {
 244     x += h*dx;
 245     h *= m_StretchingRatio;
 246     m_PrismaticGrid->GetPoints()->SetPoint(i*m_ShellPart.getNumberOfNodes() + m_ShellNodeMap[id_node], x.data());
 247   }
 248   m_PrismaticGrid->GetPoints()->SetPoint(m_NumLayers*m_ShellPart.getNumberOfNodes() + m_ShellNodeMap[id_node], x2.data());
 249   m_Grid->GetPoints()->SetPoint(id_node, x2.data());
 250 }
 251
 252 void CreateBoundaryLayerShell::createPrismaticGrid()
 253 {
 254   QVector<vtkIdType> original_triangles, shell_triangles;
 255   getSurfaceCells(m_BoundaryLayerCodes, original_triangles, m_OriginalGrid);
 256   getSurfaceCells(m_BoundaryLayerCodes, shell_triangles, m_Grid);
 257   {
 258     MeshPartition part(m_Grid);
 259     part.setCells(shell_triangles);
 260     allocateGrid(m_PrismaticGrid, (m_NumLayers + 1)*part.getNumberOfCells(), (m_NumLayers + 1)*part.getNumberOfNodes());
 261   }
 262
 263   m_ShellNodeMap.fill(-1, m_Grid->GetNumberOfPoints());
 264   m_ShellPart.setGrid(m_Grid);
 265   m_ShellPart.setCells(shell_triangles);
 266   for (int i = 0; i < m_ShellPart.getNumberOfNodes(); ++i) {
 267     m_ShellNodeMap[m_ShellPart.globalNode(i)] = i;
 268   }
 269
 270   QVector<QSet<int> > n2bc(m_PrismaticGrid->GetNumberOfPoints());
 271
 272   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 273     if (m_ShellNodeMap[id_node] != -1) {
 274       createLayerNodes(id_node);
 275       for (int i = 0; i < m_Part.n2bcGSize(id_node); ++i) {
 276         n2bc[m_ShellNodeMap[id_node]].insert(m_Part.n2bcG(id_node, i));
 277         //n2bc[m_ShellNodeMap[id_node] + m_ShellPart.getNumberOfNodes()].insert(m_Part.n2bcG(id_node, i));
 278       }
 279     }
 280   }
 281
 282   QList<QVector<vtkIdType> > adjacent_edges;
 283
 284   // create prismatic cells and prepare adjacent quad faces
 285   //
 286   foreach (vtkIdType id_cell, shell_triangles) {
 287     vtkIdType num_pts, *pts;
 288     m_Grid->GetCellPoints(id_cell, num_pts, pts);
 289     vtkIdType tri_pts[3], pri_pts[6];
 290     for (int i_pts = 0; i_pts < 3; ++i_pts) {
 291       if (m_ShellNodeMap[pts[i_pts]] < 0) {
 292         EG_BUG;
 293       }
 294       if (m_ShellNodeMap[pts[i_pts]] >= m_ShellPart.getNumberOfNodes()) {
 295         EG_BUG;
 296       }
 297       QVector<vtkIdType> edge(4);
 298       edge[1] = m_ShellNodeMap[pts[i_pts]];
 299       edge[2] = m_ShellNodeMap[pts[0]];
 300       if (i_pts < 2) {
 301         edge[2] = m_ShellNodeMap[pts[i_pts+1]];
 302       }
 303       QSet<int> edge_codes = m_LayerAdjacentBoundaryCodes;
 304       edge_codes.intersect(n2bc[edge[1]]);
 305       edge_codes.intersect(n2bc[edge[2]]);
 306       if (edge_codes.size() == 1) {
 307         edge[0] = *edge_codes.begin();
 308         edge[3] = id_cell;
 309         adjacent_edges.append(edge);
 310       }
 311       tri_pts[i_pts] = m_ShellNodeMap[pts[i_pts]];
 312     }
 313     vtkIdType id_tri = m_PrismaticGrid->InsertNextCell(VTK_TRIANGLE, 3, tri_pts);
 314     copyCellData(m_Grid, id_cell, m_PrismaticGrid, id_tri);
 315     for (int i_layer = 0; i_layer < m_NumLayers; ++i_layer) {
 316       for (int i_pts = 0; i_pts < 3; ++i_pts) {
 317         pri_pts[i_pts]     = m_ShellNodeMap[pts[i_pts]] + i_layer*m_ShellPart.getNumberOfNodes();
 318         pri_pts[i_pts + 3] = m_ShellNodeMap[pts[i_pts]] + (i_layer + 1)*m_ShellPart.getNumberOfNodes();
 319       }
 320       vtkIdType id_pri = m_PrismaticGrid->InsertNextCell(VTK_WEDGE, 6, pri_pts);
 321     }
 322   }
 323
 324   // create quads on adjacent boundary faces
 325   //
 326   EG_VTKSP(vtkUnstructuredGrid, noquad_grid);
 327   makeCopy(m_PrismaticGrid, noquad_grid);
 328   allocateGrid(m_PrismaticGrid, m_PrismaticGrid->GetNumberOfCells() + m_NumLayers*adjacent_edges.size(), m_PrismaticGrid->GetNumberOfPoints());
 329   makeCopyNoAlloc(noquad_grid, m_PrismaticGrid);
 330
 331   EG_VTKDCC(vtkIntArray, cell_code, m_PrismaticGrid, "cell_code");
 332
 333   foreach (QVector<vtkIdType> edge, adjacent_edges) {
 334     vtkIdType qua_pts[4];
 335     for (int i_layer = 0; i_layer < m_NumLayers; ++i_layer) {
 336       qua_pts[0] = edge[2] + i_layer*m_ShellPart.getNumberOfNodes();
 337       qua_pts[1] = edge[1] + i_layer*m_ShellPart.getNumberOfNodes();
 338       qua_pts[2] = edge[1] + (i_layer + 1)*m_ShellPart.getNumberOfNodes();
 339       qua_pts[3] = edge[2] + (i_layer + 1)*m_ShellPart.getNumberOfNodes();
 340       vtkIdType id_qua = m_PrismaticGrid->InsertNextCell(VTK_QUAD, 4, qua_pts);
 341       copyCellData(m_Grid, edge[3], m_PrismaticGrid, id_qua);
 342       cell_code->SetValue(id_qua, edge[0]);
 343     }
 344   }
 345 }
 346
 347 void CreateBoundaryLayerShell::reduceSurface()
 348 {
 349   RemovePoints remove_points;
 350   MeshPartition part;
 351   part.setGrid(m_Grid);
 352   part.setAllCells();
 353   remove_points.setMeshPartition(part);
 354   remove_points.setBoundaryCodes(m_LayerAdjacentBoundaryCodes);
 355   remove_points.setUpdatePSPOn();
 356   remove_points.setThreshold(3);
 357   QVector<bool> fix(m_Grid->GetNumberOfPoints(), false);
 358   EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type");
 359   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 360     // reset all node types to EG_SIMPLE_VERTEX to trigger a recalculation of the node types
 361     node_type->SetValue(id_node, EG_SIMPLE_VERTEX);
 362     for (int i = 0; i < part.n2cGSize(id_node); ++i) {
 363       if (m_Grid->GetCellType(part.n2cGG(id_node, i)) == VTK_QUAD) {
 364         fix[id_node] = true;
 365       }
 366     }
 367   }
 368   remove_points.fixNodes(fix);
 369   remove_points();
 370 }
 371
 372 void CreateBoundaryLayerShell::smoothSurface()
 373 {
 374   LaplaceSmoother smooth;
 375   MeshPartition part;
 376   part.setGrid(m_Grid);
 377   part.setAllCells();
 378   smooth.setMeshPartition(part);
 379   smooth.setNumberOfIterations(2);
 380   smooth.setBoundaryCodes(m_LayerAdjacentBoundaryCodes);
 381   QVector<bool> fix(m_Grid->GetNumberOfPoints(), false);
 382   EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type");
 383   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 384     // reset all node types to EG_SIMPLE_VERTEX to trigger a recalculation of the node types
 385     node_type->SetValue(id_node, EG_SIMPLE_VERTEX);
 386     for (int i = 0; i < part.n2cGSize(id_node); ++i) {
 387       if (m_Grid->GetCellType(part.n2cGG(id_node, i)) == VTK_QUAD) {
 388         fix[id_node] = true;
 389       }
 390     }
 391   }
 392   /*
 393   for (int layer = 0; layer < 3; ++layer) {
 394     QVector<bool> tmp = fix;
 395     for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 396       if (!tmp[id_node]) {
 397         for (int i = 0; i < part.n2nGSize(id_node); ++i) {
 398           fix[part.n2nGG(id_node, i)] = false;
 399         }
 400       }
 401     }
 402   }
 403   */
 404   smooth.fixNodes(fix);
 405   smooth();
 406 }
 407
 408
 409 void CreateBoundaryLayerShell::operate()
 410 {
 411   prepare();
 412   writeBoundaryLayerVectors("blayer");
 413   createPrismaticGrid();
 414   m_Success = true;
 415   m_Part.trackGrid(m_Grid);
 416
 417   foreach (int bc, m_LayerAdjacentBoundaryCodes) {
 418     QList<vtkIdType> bad_nodes = correctAdjacentBC(bc);
 419     if (bad_nodes.size() > 0) {
 420       bool fixable = true;
 421
 422       QVector<bool> is_bad_cell(m_Grid->GetNumberOfCells(), false);
 423       foreach (vtkIdType id_node, bad_nodes) {
 424         if (m_Part.n2bcGSize(id_node) != 1) {
 425           cout << "node " << id_node << " cannot be fixed." << endl;
 426           fixable = false;
 427           break;
 428         }
 429         for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
 430           is_bad_cell[m_Part.n2cGG(id_node, i)] = true;
 431         }
 432       }
 433       if (fixable) {
 434         QList<vtkIdType> bad_cells;
 435         EG_FORALL_CELLS (id_cell, m_Grid) {
 436           if (is_bad_cell[id_cell]) {
 437             bad_cells << id_cell;
 438           }
 439         }
 440
 441         DeleteCells del;
 442         del.setGrid(m_Grid);
 443         del.setCellsToDelete(bad_cells);
 444         del();
 445         StitchHoles stitch(bc);
 446         stitch.setGrid(m_Grid);
 447         stitch();
 448       } else {
 449         m_Success = false;
 450         cout << "adjacent patch cannot be corrected!" << endl;
 451       }
 452
 453     }
 454   }
 455   SwapTriangles swap;
 456   swap.setGrid(m_Grid);
 457   QSet<int> swap_codes = getAllBoundaryCodes(m_Grid);
 458   swap_codes -= m_LayerAdjacentBoundaryCodes;
 459   swap.setBoundaryCodes(swap_codes);
 460   swap.setVerboseOff();
 461
 462   for (int iter = 0; iter < 5; ++iter) {
 463     cout << "correcting adjacent boundaries\n" << "  iteration: " << iter + 1 << endl;
 464     swap();
 465     smoothSurface();
 466     reduceSurface();
 467     swap();
 468   }
 469
 470 }
 471