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 #include "cgaltricadinterface.h"
  35 #include "correctsurfaceorientation.h"
  36
  37
  38 CreateBoundaryLayerShell::DeleteBadNodes::DeleteBadNodes(QList<vtkIdType> bad_nodes)
  39 {
  40   m_BadNodes = bad_nodes;
  41   m_PerformGeometricChecks = false;
  42 }
  43
  44 bool CreateBoundaryLayerShell::DeleteBadNodes::checkEdge(vtkIdType id_node1, vtkIdType id_node2)
  45 {
  46   EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type");
  47   if (m_BadNodes.contains(id_node1) && !m_BadNodes.contains(id_node2)) {
  48     return true;
  49   }
  50   return false;
  51 }
  52
  53 CreateBoundaryLayerShell::CreateBoundaryLayerShell()
  54 {
  55   m_RestGrid      = vtkSmartPointer<vtkUnstructuredGrid>::New();
  56   m_OriginalGrid  = vtkSmartPointer<vtkUnstructuredGrid>::New();
  57   m_PrismaticGrid = vtkSmartPointer<vtkUnstructuredGrid>::New();
  58   m_Success = false;
  59 }
  60
  61 void CreateBoundaryLayerShell::prepare()
  62 {
  63   m_Part.trackGrid(m_Grid);
  64
  65   m_OriginalNodeNormals.resize(m_Grid->GetNumberOfPoints());
  66   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
  67     m_OriginalNodeNormals[id_node] = m_Part.globalNormal(id_node);
  68   }
  69
  70   DeleteVolumeGrid delete_volume;
  71   delete_volume.setGrid(m_Grid);
  72   delete_volume.setAllCells();
  73   delete_volume();
  74
  75   readSettings();
  76   setAllCells();
  77   getSurfaceCells(m_BoundaryLayerCodes, layer_cells, m_Grid);
  78
  79   // fill m_LayerAdjacentBoundaryCodes
  80   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
  81   foreach (vtkIdType id_cell, layer_cells) {
  82     for (int i = 0; i < m_Part.c2cGSize(id_cell); ++i) {
  83       vtkIdType id_neigh = m_Part.c2cGG(id_cell, i);
  84       int bc = cell_code->GetValue(id_neigh);
  85       if (!m_BoundaryLayerCodes.contains(bc)) {
  86         m_LayerAdjacentBoundaryCodes.insert(bc);
  87       }
  88     }
  89   }
  90
  91   // compute normals and origins of adjacent planes
  92   m_LayerAdjacentNormals.clear();
  93   m_LayerAdjacentOrigins.clear();
  94   foreach (int bc, m_LayerAdjacentBoundaryCodes) {
  95     double L = EG_LARGE_REAL;
  96     vec3_t n0(0, 0, 0);
  97     vec3_t x0(0, 0, 0);
  98     double total_area = 0;
  99     for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 100       if (isSurface(id_cell, m_Grid) && cell_code->GetValue(id_cell) == bc) {
 101         vec3_t n = cellNormal(m_Grid, id_cell);
 102         double A = n.abs();
 103         total_area += A;
 104         n0 += n;
 105         x0 += A*cellCentre(m_Grid, id_cell);
 106         L = min(L, sqrt(4*A/sqrt(3.0)));
 107       }
 108     }
 109     n0.normalise();
 110     x0 *= 1.0/total_area;
 111     m_LayerAdjacentNormals[bc] = n0;
 112     m_LayerAdjacentOrigins[bc] = x0;
 113   }
 114
 115   computeBoundaryLayerVectors();
 116   makeCopy(m_Grid, m_OriginalGrid);
 117 }
 118
 119 QList<vtkIdType> CreateBoundaryLayerShell::findBadNodes(int bc)
 120 {
 121   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
 122   EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type");
 123   CgalTriCadInterface cad(m_ShellGrid);
 124   QList<vtkIdType> bad_nodes;
 125   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 126     if (!m_BoundaryLayerNode[id_node]){
 127       bool node_has_bc = false;
 128       bool all_bcs_adjacent = true;
 129       for (int i = 0; i < m_Part.n2bcGSize(id_node); ++i) {
 130         int bcn = m_Part.n2bcG(id_node, i);
 131         if (bcn == bc) {
 132           node_has_bc = true;
 133         }
 134         if (!m_LayerAdjacentBoundaryCodes.contains(bcn)) {
 135           all_bcs_adjacent = false;
 136         }
 137       }
 138       if (node_has_bc && all_bcs_adjacent && node_type->GetValue(id_node) != EG_FIXED_VERTEX) {
 139         vec3_t x1;
 140         m_Grid->GetPoint(id_node, x1.data());
 141         vec3_t x2 = cad.snap(x1);
 142         if ((x2-x1)*cad.getLastNormal() <= 0) {
 143           bad_nodes << id_node;
 144         }
 145       }
 146     }
 147   }
 148   return bad_nodes;
 149 }
 150
 151 QList<vtkIdType> CreateBoundaryLayerShell::correctAdjacentBC(int bc, int num_levels)
 152 {
 153   cout << "correcting boundary \"" << qPrintable(GuiMainWindow::pointer()->getBC(bc).getName()) << "\" with " << num_levels << " levels" << endl;
 154
 155   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
 156   EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type");
 157   double scal_min = -1;
 158   int    count    = 0;
 159
 160   SurfaceMeshSmoother smooth;
 161   smooth.setGrid(m_Grid);
 162   smooth.useSimpleCentreScheme();
 163
 164   // mark num_levels levels of nodes next to the boundary layer
 165   QVector<bool> marked(m_Grid->GetNumberOfPoints(), false);
 166   QList<vtkIdType> marked_nodes;
 167   EG_FORALL_NODES(id_node, m_Grid) {
 168     if (m_BoundaryLayerNode[id_node]) {
 169       marked[id_node] = true;
 170       marked_nodes << id_node;
 171     }
 172   }
 173   for (int level = 0; level < num_levels; ++level) {
 174     QList<vtkIdType> new_marked_nodes;
 175     foreach (vtkIdType id_node, marked_nodes) {
 176       for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
 177         vtkIdType id_neigh = m_Part.n2nGG(id_node, i);
 178         if (!marked[id_neigh]) {
 179           marked[id_neigh] = true;
 180           new_marked_nodes << id_neigh;
 181         }
 182       }
 183     }
 184     marked_nodes = new_marked_nodes;
 185   }
 186
 187   QList<vtkIdType> cells;
 188   EG_FORALL_CELLS(id_cell, m_Grid) {
 189     if (isSurface(id_cell, m_Grid)) {
 190       if (cell_code->GetValue(id_cell) == bc) {
 191         cells << id_cell;
 192       }
 193     }
 194   }
 195
 196   CadInterface *cad = GuiMainWindow::pointer()->getCadInterface(bc);
 197   smooth.setCells(cells);
 198   smooth.prepareCadInterface(cad);
 199
 200   QList<vtkIdType> bad_nodes;
 201   int new_num_bad = m_Grid->GetNumberOfPoints();
 202   int old_num_bad = 0;
 203   updateNodeInfo();
 204   do {
 205     old_num_bad = new_num_bad;
 206     scal_min = 1;
 207     bad_nodes.clear();
 208     for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 209       if (!m_BoundaryLayerNode[id_node]){
 210         bool node_has_bc = false;
 211         bool all_bcs_adjacent = true;
 212         for (int i = 0; i < m_Part.n2bcGSize(id_node); ++i) {
 213           int bcn = m_Part.n2bcG(id_node, i);
 214           if (bcn == bc) {
 215             node_has_bc = true;
 216           }
 217           if (!m_LayerAdjacentBoundaryCodes.contains(bcn)) {
 218             all_bcs_adjacent = false;
 219           }
 220         }
 221         if (node_has_bc && all_bcs_adjacent && marked[id_node] && node_type->GetValue(id_node) != EG_FIXED_VERTEX) {
 222
 223           // one boundary condition and EG_SIMPLE_VERTEX
 224           //
 225           if (m_Part.n2bcGSize(id_node) == 1) {
 226             vec3_t xs(0, 0, 0); // = smooth.smoothNode(id_node);
 227             int count = 0;
 228             for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
 229               vtkIdType id_cell = m_Part.n2cGG(id_node, i);
 230               if (cell_code->GetValue(id_cell) == bc) {
 231                 vec3_t x = cellCentre(m_Grid, id_cell);
 232                 xs += x;
 233                 ++count;
 234               }
 235             }
 236             if (count > 0) {
 237               xs *= 1.0/count;
 238               if (node_type->GetValue(id_node) == EG_SIMPLE_VERTEX) {
 239                 vec3_t n = m_Part.globalNormal(id_node);
 240                 if (m_OriginalNodeNormals[id_node]*n < 0) {
 241                   n = m_OriginalNodeNormals[id_node];
 242                 }
 243                 xs = cad->snapWithNormal(xs, m_OriginalNodeNormals[id_node]);
 244               } else {
 245                 xs = cad->snapToEdge(xs);
 246               }
 247               if (!checkVector(xs)) {
 248                 EG_ERR_RETURN("error while correcting adjacent boundaries");
 249               }
 250               m_Grid->GetPoints()->SetPoint(id_node, xs.data());
 251             }
 252           }
 253
 254           // one boundary condition and EG_FEATURE_EDGE_VERTEX
 255           //
 256           else if (m_Part.n2bcGSize(id_node) == 1 && node_type->GetValue(id_node) == EG_FEATURE_EDGE_VERTEX) {
 257             int count = 0;
 258             vec3_t xs(0, 0, 0);
 259             for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
 260               vtkIdType id_neigh = m_Part.n2nGG(id_node, i);
 261               if (node_type->GetValue(id_neigh) == EG_FEATURE_EDGE_VERTEX) {
 262                 vec3_t x;
 263                 m_Grid->GetPoint(id_neigh, x.data());
 264                 xs += x;
 265                 ++count;
 266               }
 267             }
 268             if (count < 2) {
 269               vec3_t x;
 270               m_Grid->GetPoint(id_node, x.data());
 271               xs += x;
 272               ++count;
 273             }
 274             xs *= 1.0/count;
 275             xs = cad->snapToEdge(xs);
 276             m_Grid->GetPoints()->SetPoint(id_node, xs.data());
 277           }
 278
 279           // two boundary conditions and no EG_FIXED_VERTEX
 280           //
 281           else if (m_Part.n2bcGSize(id_node) == 2 && node_type->GetValue(id_node) != EG_FIXED_VERTEX) {
 282             int count = 0;
 283             vec3_t xs(0, 0, 0);
 284             for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
 285               vtkIdType id_neigh = m_Part.n2nGG(id_node, i);
 286               if (m_Part.n2bcGSize(id_neigh) > 1) {
 287                 vec3_t x;
 288                 m_Grid->GetPoint(id_neigh, x.data());
 289                 xs += x;
 290                 ++count;
 291               }
 292             }
 293             if (count < 2) {
 294               vec3_t x;
 295               m_Grid->GetPoint(id_node, x.data());
 296               xs += x;
 297               ++count;
 298             }
 299             xs *= 1.0/count;
 300             xs = cad->snapToEdge(xs);
 301             m_Grid->GetPoints()->SetPoint(id_node, xs.data());
 302           }
 303
 304           bool node_bad = false;
 305           for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
 306             vtkIdType id_cell = m_Part.n2cGG(id_node, i);
 307             if (isSurface(id_cell, m_Grid)) {
 308               if (cell_code->GetValue(id_cell) == bc) {
 309                 CadInterface *cad = GuiMainWindow::pointer()->getCadInterface(cell_code->GetValue(id_cell));
 310                 cad->snap(cellCentre(m_Grid, id_cell));
 311                 vec3_t n = cellNormal(m_Grid, id_cell);
 312                 n.normalise();
 313                 double scal = n*cad->getLastNormal();
 314                 scal_min = min(scal_min, scal);
 315                 if (scal < 0.5 && !node_bad) {
 316                   bad_nodes << id_node;
 317                   node_bad = true;
 318                 }
 319               }
 320             }
 321           }
 322         }
 323       }
 324     }
 325     new_num_bad = bad_nodes.size();
 326     ++count;
 327   } while (scal_min < 0.5 && count < 20);
 328   cout << "  " << bad_nodes.size() << " node defects" << endl;
 329   return bad_nodes;
 330 }
 331
 332 void CreateBoundaryLayerShell::createLayerNodes(vtkIdType id_node)
 333 {
 334   vec3_t x1;
 335   m_Grid->GetPoint(id_node, x1.data());
 336   m_PrismaticGrid->GetPoints()->SetPoint(m_ShellNodeMap[id_node], x1.data());
 337   vec3_t x2 = x1 + m_BoundaryLayerVectors[id_node];
 338
 339   double H  = m_BoundaryLayerVectors[id_node].abs();
 340   double h  = H*(1.0 - m_StretchingRatio)/(1.0 - pow(m_StretchingRatio, m_NumLayers));
 341   vec3_t dx = (1.0/H)*m_BoundaryLayerVectors[id_node];
 342   vec3_t x  = x1;
 343   m_PrismaticGrid->GetPoints()->SetPoint(m_ShellNodeMap[id_node], x1.data());
 344   for (int i = 1; i < m_NumLayers; ++i) {
 345     x += h*dx;
 346     h *= m_StretchingRatio;
 347     m_PrismaticGrid->GetPoints()->SetPoint(i*m_ShellPart.getNumberOfNodes() + m_ShellNodeMap[id_node], x.data());
 348   }
 349   m_PrismaticGrid->GetPoints()->SetPoint(m_NumLayers*m_ShellPart.getNumberOfNodes() + m_ShellNodeMap[id_node], x2.data());
 350   m_Grid->GetPoints()->SetPoint(id_node, x2.data());
 351 }
 352
 353 void CreateBoundaryLayerShell::createPrismaticGrid()
 354 {
 355   QVector<vtkIdType> original_triangles, shell_triangles;
 356   getSurfaceCells(m_BoundaryLayerCodes, original_triangles, m_OriginalGrid);
 357   getSurfaceCells(m_BoundaryLayerCodes, shell_triangles, m_Grid);
 358   {
 359     MeshPartition part(m_Grid);
 360     part.setCells(shell_triangles);
 361     allocateGrid(m_PrismaticGrid, (m_NumLayers + 1)*part.getNumberOfCells(), (m_NumLayers + 1)*part.getNumberOfNodes());
 362   }
 363
 364   m_ShellNodeMap.fill(-1, m_Grid->GetNumberOfPoints());
 365   m_ShellPart.setGrid(m_Grid);
 366   m_ShellPart.setCells(shell_triangles);
 367   for (int i = 0; i < m_ShellPart.getNumberOfNodes(); ++i) {
 368     m_ShellNodeMap[m_ShellPart.globalNode(i)] = i;
 369   }
 370
 371   QVector<QSet<int> > n2bc(m_PrismaticGrid->GetNumberOfPoints());
 372
 373   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 374     if (m_ShellNodeMap[id_node] != -1) {
 375       createLayerNodes(id_node);
 376       for (int i = 0; i < m_Part.n2bcGSize(id_node); ++i) {
 377         n2bc[m_ShellNodeMap[id_node]].insert(m_Part.n2bcG(id_node, i));
 378         //n2bc[m_ShellNodeMap[id_node] + m_ShellPart.getNumberOfNodes()].insert(m_Part.n2bcG(id_node, i));
 379       }
 380     }
 381   }
 382
 383   QList<QVector<vtkIdType> > adjacent_edges;
 384
 385   // create prismatic cells and prepare adjacent quad faces
 386   //
 387   foreach (vtkIdType id_cell, shell_triangles) {
 388     vtkIdType num_pts, *pts;
 389     m_Grid->GetCellPoints(id_cell, num_pts, pts);
 390     vtkIdType tri_pts[3], pri_pts[6];
 391     for (int i_pts = 0; i_pts < 3; ++i_pts) {
 392       if (m_ShellNodeMap[pts[i_pts]] < 0) {
 393         EG_BUG;
 394       }
 395       if (m_ShellNodeMap[pts[i_pts]] >= m_ShellPart.getNumberOfNodes()) {
 396         EG_BUG;
 397       }
 398       QVector<vtkIdType> edge(4);
 399       edge[1] = m_ShellNodeMap[pts[i_pts]];
 400       edge[2] = m_ShellNodeMap[pts[0]];
 401       if (i_pts < 2) {
 402         edge[2] = m_ShellNodeMap[pts[i_pts+1]];
 403       }
 404       QSet<int> edge_codes = m_LayerAdjacentBoundaryCodes;
 405       edge_codes.intersect(n2bc[edge[1]]);
 406       edge_codes.intersect(n2bc[edge[2]]);
 407       if (edge_codes.size() == 1) {
 408         edge[0] = *edge_codes.begin();
 409         edge[3] = id_cell;
 410         adjacent_edges.append(edge);
 411       }
 412       tri_pts[i_pts] = m_ShellNodeMap[pts[i_pts]];
 413     }
 414     vtkIdType id_tri = m_PrismaticGrid->InsertNextCell(VTK_TRIANGLE, 3, tri_pts);
 415     copyCellData(m_Grid, id_cell, m_PrismaticGrid, id_tri);
 416     for (int i_layer = 0; i_layer < m_NumLayers; ++i_layer) {
 417       for (int i_pts = 0; i_pts < 3; ++i_pts) {
 418         pri_pts[i_pts]     = m_ShellNodeMap[pts[i_pts]] + i_layer*m_ShellPart.getNumberOfNodes();
 419         pri_pts[i_pts + 3] = m_ShellNodeMap[pts[i_pts]] + (i_layer + 1)*m_ShellPart.getNumberOfNodes();
 420       }
 421       vtkIdType id_pri = m_PrismaticGrid->InsertNextCell(VTK_WEDGE, 6, pri_pts);
 422     }
 423   }
 424
 425   // create quads on adjacent boundary faces
 426   //
 427   EG_VTKSP(vtkUnstructuredGrid, noquad_grid);
 428   makeCopy(m_PrismaticGrid, noquad_grid);
 429   allocateGrid(m_PrismaticGrid, m_PrismaticGrid->GetNumberOfCells() + m_NumLayers*adjacent_edges.size(), m_PrismaticGrid->GetNumberOfPoints());
 430   makeCopyNoAlloc(noquad_grid, m_PrismaticGrid);
 431
 432   EG_VTKDCC(vtkIntArray, cell_code, m_PrismaticGrid, "cell_code");
 433
 434   foreach (QVector<vtkIdType> edge, adjacent_edges) {
 435     vtkIdType qua_pts[4];
 436     for (int i_layer = 0; i_layer < m_NumLayers; ++i_layer) {
 437       qua_pts[0] = edge[2] + i_layer*m_ShellPart.getNumberOfNodes();
 438       qua_pts[1] = edge[1] + i_layer*m_ShellPart.getNumberOfNodes();
 439       qua_pts[2] = edge[1] + (i_layer + 1)*m_ShellPart.getNumberOfNodes();
 440       qua_pts[3] = edge[2] + (i_layer + 1)*m_ShellPart.getNumberOfNodes();
 441       vtkIdType id_qua = m_PrismaticGrid->InsertNextCell(VTK_QUAD, 4, qua_pts);
 442       copyCellData(m_Grid, edge[3], m_PrismaticGrid, id_qua);
 443       cell_code->SetValue(id_qua, edge[0]);
 444     }
 445   }
 446 }
 447
 448 void CreateBoundaryLayerShell::reduceSurface()
 449 {
 450   RemovePoints remove_points;
 451   MeshPartition part;
 452   part.setGrid(m_Grid);
 453   part.setAllCells();
 454   remove_points.setMeshPartition(part);
 455   remove_points.setBoundaryCodes(m_LayerAdjacentBoundaryCodes);
 456   remove_points.setUpdatePSPOn();
 457   remove_points.setThreshold(3);
 458   QVector<bool> fix(m_Grid->GetNumberOfPoints(), false);
 459   EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type");
 460   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 461     // reset all node types to EG_SIMPLE_VERTEX to trigger a recalculation of the node types
 462     node_type->SetValue(id_node, EG_SIMPLE_VERTEX);
 463     for (int i = 0; i < part.n2cGSize(id_node); ++i) {
 464       if (m_Grid->GetCellType(part.n2cGG(id_node, i)) == VTK_QUAD) {
 465         fix[id_node] = true;
 466       }
 467     }
 468   }
 469   remove_points.fixNodes(fix);
 470   remove_points();
 471 }
 472
 473 void CreateBoundaryLayerShell::smoothSurface()
 474 {
 475   LaplaceSmoother smooth;
 476   MeshPartition part;
 477   part.setGrid(m_Grid);
 478   part.setAllCells();
 479   smooth.setMeshPartition(part);
 480   smooth.setNumberOfIterations(2);
 481   smooth.setBoundaryCodes(m_LayerAdjacentBoundaryCodes);
 482   QVector<bool> fix(m_Grid->GetNumberOfPoints(), false);
 483   EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type");
 484   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 485     // reset all node types to EG_SIMPLE_VERTEX to trigger a recalculation of the node types
 486     node_type->SetValue(id_node, EG_SIMPLE_VERTEX);
 487     for (int i = 0; i < part.n2cGSize(id_node); ++i) {
 488       if (m_Grid->GetCellType(part.n2cGG(id_node, i)) == VTK_QUAD) {
 489         fix[id_node] = true;
 490       }
 491     }
 492   }
 493   smooth.fixNodes(fix);
 494   smooth();
 495 }
 496
 497 void CreateBoundaryLayerShell::operate()
 498 {
 499   prepare();
 500   createPrismaticGrid();
 501   m_Success = true;
 502   m_Part.trackGrid(m_Grid);
 503
 504   foreach (int bc, m_LayerAdjacentBoundaryCodes) {
 505     QList<vtkIdType> bad_nodes;
 506     int levels = 1;
 507     QVector<vec3_t> x_save(m_Grid->GetNumberOfPoints());
 508     do {
 509       ++levels;
 510       for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 511         m_Grid->GetPoint(id_node, x_save[id_node].data());
 512       }
 513       QList<vtkIdType> last_bad = bad_nodes;
 514       bad_nodes = correctAdjacentBC(bc, levels);
 515       if (bad_nodes.size() >= last_bad.size() && last_bad.size() > 0) {
 516         bad_nodes = last_bad;
 517         for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 518           m_Grid->GetPoints()->SetPoint(id_node, x_save[id_node].data());
 519         }
 520         break;
 521       }
 522     } while (bad_nodes.size() > 0);
 523   }
 524
 525   QList<vtkIdType> bad_cells;
 526   foreach (int bc, m_LayerAdjacentBoundaryCodes) {
 527     QList<vtkIdType> bad_nodes = findBadNodes(bc);
 528     if (bad_nodes.size() > 0) {
 529       bool fixable = true;
 530       QVector<bool> is_bad_cell(m_Grid->GetNumberOfCells(), false);
 531       foreach (vtkIdType id_node, bad_nodes) {
 532         if (m_Part.n2bcGSize(id_node) != 1) {
 533           cout << "node " << id_node << " cannot be fixed." << endl;
 534           fixable = false;
 535           break;
 536         }
 537         for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
 538           is_bad_cell[m_Part.n2cGG(id_node, i)] = true;
 539         }
 540       }
 541       if (fixable) {
 542         EG_FORALL_CELLS (id_cell, m_Grid) {
 543           if (is_bad_cell[id_cell]) {
 544             bad_cells << id_cell;
 545           }
 546         }
 547       } else {
 548         m_Success = false;
 549         cout << "adjacent patch cannot be corrected!" << endl;
 550         return;
 551       }
 552     }
 553   }
 554
 555   cout << "deleting "<< bad_cells.size() << "bad cells" << endl;
 556   DeleteCells del;
 557   del.setGrid(m_Grid);
 558   del.setCellsToDelete(bad_cells);
 559   del();
 560
 561   QMap<int,int> orgdir, curdir, voldir;
 562   EG_VTKDCC(vtkIntArray, cell_orgdir, m_Grid, "cell_orgdir");
 563   EG_VTKDCC(vtkIntArray, cell_curdir, m_Grid, "cell_curdir");
 564   EG_VTKDCC(vtkIntArray, cell_voldir, m_Grid, "cell_voldir");
 565   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
 566   foreach (int bc, m_LayerAdjacentBoundaryCodes) {
 567     for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 568       if (cell_code->GetValue(id_cell) == bc) {
 569         orgdir[bc] = cell_orgdir->GetValue(id_cell);
 570         curdir[bc] = cell_curdir->GetValue(id_cell);
 571         voldir[bc] = cell_voldir->GetValue(id_cell);
 572         break;
 573       }
 574     }
 575     StitchHoles stitch(bc);
 576     stitch.setGrid(m_Grid);
 577     stitch();
 578   }
 579   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 580     int bc = cell_code->GetValue(id_cell);
 581     if (m_LayerAdjacentBoundaryCodes.contains(bc)) {
 582       cell_orgdir->SetValue(id_cell, orgdir[bc]);
 583       cell_curdir->SetValue(id_cell, curdir[bc]);
 584       cell_voldir->SetValue(id_cell, voldir[bc]);
 585     }
 586   }
 587
 588   // check surface orientation
 589   {
 590     CorrectSurfaceOrientation surf_check;
 591     surf_check.setGrid(m_Grid);
 592     surf_check.setStart(0);
 593     for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 594       if (!m_LayerAdjacentBoundaryCodes.contains(cell_code->GetValue(id_cell))) {
 595         surf_check.setStart(id_cell);
 596         break;
 597       }
 598     }
 599     surf_check();
 600   }
 601
 602   SwapTriangles swap;
 603   swap.setGrid(m_Grid);
 604   QSet<int> swap_codes = getAllBoundaryCodes(m_Grid);
 605   swap_codes -= m_LayerAdjacentBoundaryCodes;
 606   swap.setBoundaryCodes(swap_codes);
 607   swap.setVerboseOff();
 608
 609   for (int iter = 0; iter < 5; ++iter) {
 610     cout << "correcting adjacent boundaries\n" << "  iteration: " << iter + 1 << endl;
 611     swap();
 612     smoothSurface();
 613     reduceSurface();
 614     swap();
 615   }
 616 }
 617