src/guicreateboundarylayer.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 "guicreateboundarylayer.h"
  24 #include "guimainwindow.h"
  25 #include "seedsimpleprismaticlayer.h"
  26 #include "gridsmoother.h"
  27 #include "createvolumemesh.h"
  28 #include "swaptriangles.h"
  29 #include "deletetetras.h"
  30 #include "deletecells.h"
  31 #include "meshpartition.h"
  32
  33 GuiCreateBoundaryLayer::GuiCreateBoundaryLayer()
  34 {
  35   getSet("boundary layer", "maximal relative error", 0.01, err_max);
  36   getSet("boundary layer", "maximal number of smoothing iterations", 10, max_iter);
  37 }
  38
  39 void GuiCreateBoundaryLayer::before()
  40 {
  41   populateBoundaryCodes(ui.listWidgetBC);
  42   populateVolumes(ui.listWidgetVC);
  43 }
  44
  45 void GuiCreateBoundaryLayer::operate()
  46 {
  47   getSelectedItems(ui.listWidgetBC, boundary_codes);
  48   QString volume_name = getSelectedVolume(ui.listWidgetVC);
  49   VolumeDefinition V = GuiMainWindow::pointer()->getVol(volume_name);
  50   foreach (int bc, boundary_codes) {
  51     if (V.getSign(bc) == 0) {
  52       QString msg;
  53       msg.setNum(bc);
  54       msg = "Boundary code " + msg + " is not part of the volume '" + volume_name +"'.";
  55       EG_ERR_RETURN(msg);
  56     }
  57   }
  58
  59   EG_VTKSP(vtkUnstructuredGrid, rest_grid);
  60   {
  61     EG_VTKSP(vtkUnstructuredGrid, vol_grid);
  62     MeshPartition volume(volume_name);
  63     MeshPartition rest(grid);
  64     rest.setRemainder(volume);
  65     volume.setVolumeOrientation();
  66     volume.extractToVtkGrid(vol_grid);
  67     rest.extractToVtkGrid(rest_grid);
  68     makeCopy(vol_grid, grid);
  69   }
  70
  71   setAllCells();
  72   l2g_t  nodes = getPartNodes();
  73   l2g_t  cells = getPartCells();
  74   g2l_t _nodes = getPartLocalNodes();
  75   l2l_t  n2c   = getPartN2C();
  76   getSurfaceCells(boundary_codes, layer_cells, grid);
  77
  78   cout << "\n\ncreating boundary layer mesh)" << endl;
  79
  80   {
  81     EG_VTKDCN(vtkIntArray, node_status, grid, "node_status");
  82     EG_VTKDCN(vtkIntArray, node_layer,  grid, "node_layer");
  83     EG_VTKDCC(vtkIntArray, bc,          grid, "cell_code");
  84     foreach(vtkIdType id_node, nodes) {
  85       node_status->SetValue(id_node, 0);
  86       QSet<int> bcs;
  87       foreach (int i_neigh_cells, n2c[_nodes[id_node]]) {
  88         vtkIdType id_neigh_cell = cells[i_neigh_cells];
  89         if (isSurface(id_neigh_cell, grid)) {
  90           if (boundary_codes.contains(bc->GetValue(id_neigh_cell))) {
  91             bcs.insert(bc->GetValue(id_neigh_cell));
  92           }
  93         }
  94       }
  95       if (bcs.size() >= 2) {
  96         node_status->SetValue(id_node, 1);
  97       }
  98       node_layer->SetValue(id_node, -1);
  99     }
 100     foreach (vtkIdType id_cell, layer_cells) {
 101       vtkIdType N_pts, *pts;
 102       grid->GetCellPoints(id_cell, N_pts, pts);
 103       for (int i_pts = 0; i_pts < N_pts; ++i_pts) {
 104         node_layer->SetValue(pts[i_pts], 0);
 105       }
 106     }
 107   }
 108
 109   cout << "preparing prismatic layer" << endl;
 110
 111   GridSmoother smooth;
 112   smooth.setGrid(grid);
 113   smooth.setBoundaryCodes(boundary_codes);
 114   smooth.prismsOn();
 115   //smooth.setNumIterations(5);
 116
 117   SeedSimplePrismaticLayer seed_layer;
 118
 119   CreateVolumeMesh vol;
 120   vol.setGrid(grid);
 121   SwapTriangles swap;
 122   swap.setGrid(grid);
 123   swap.setBoundaryCodes(boundary_codes);
 124   DeleteTetras del;
 125   del.setGrid(grid);
 126
 127   seed_layer.setGrid(grid);
 128   del();
 129   vol();
 130   seed_layer.setAllCells();
 131   seed_layer.setLayerCells(layer_cells);
 132   seed_layer.setBoundaryCodes(boundary_codes);
 133   seed_layer();
 134   seed_layer.getLayerCells(layer_cells);
 135
 136   for (int j = 0; j < max_iter; ++j) {
 137     cout << "improving prismatic layer -> iteration " << j+1 << "/" << max_iter << endl;
 138     smooth.setAllCells();
 139     smooth();
 140     del.setAllCells();
 141     del();
 142     swap();
 143     vol.setTraceCells(layer_cells);
 144     vol();
 145     vol.getTraceCells(layer_cells);
 146     if (smooth.improvement() < err_max) break;
 147   }
 148   //smooth.setAllCells();
 149   //smooth();
 150
 151   /*
 152   del();
 153   vol();
 154   cout << "finalising prismatic layer" << endl;
 155   smooth.setAllCells();
 156   smooth();
 157   */
 158
 159   {
 160     EG_VTKDCC(vtkIntArray, cell_code, grid, "cell_code");
 161     for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
 162       if (isVolume(id_cell, grid)) {
 163         cell_code->SetValue(id_cell, V.getVC());
 164       }
 165     }
 166   }
 167
 168   {
 169     MeshPartition volume(grid, true);
 170     MeshPartition rest(rest_grid, true);
 171     volume.addPartition(rest);
 172   }
 173   resetOrientation(grid);
 174   createIndices(grid);
 175 }