src/libengrid/createhexibmesh.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 "createhexibmesh.h"
  22 #include "guimainwindow.h"
  23 #include "updatedesiredmeshdensity.h"
  24 #include "deletevolumegrid.h"
  25
  26 CreateHexIbMesh::CreateHexIbMesh()
  27 {
  28   m_MinDim = 30;
  29   m_InsidePosition = vec3_t(0, 0, 0);
  30 }
  31
  32 double CreateHexIbMesh::meshSize(vtkIdType id_face)
  33 {
  34   EG_VTKDCN(vtkDoubleArray, cl, m_Grid, "node_meshdensity_desired");
  35   vtkIdType *pts, num_pts;
  36   m_Grid->GetCellPoints(id_face, num_pts, pts);
  37   double h = 0;
  38   for (int i = 0; i < 3; ++i) {
  39     h = max(h,cl->GetValue(pts[i]));
  40   }
  41   return h;
  42 }
  43
  44 QString CreateHexIbMesh::bigIntText(long long int N)
  45 {
  46   QString text, num;
  47   text.setNum(N%1000);
  48   N /= 1000;
  49   int L = 3;
  50   while (N > 0) {
  51     num.setNum(N%1000);
  52     text = num + "," + text.rightJustified(L, '0');
  53     L += 4;
  54     N /= 1000;
  55   }
  56   return text;
  57 }
  58
  59 double CreateHexIbMesh::meshSize(const QList<vtkIdType> &faces)
  60 {
  61   double mesh_size = m_MaxEdgeLength;
  62   foreach (vtkIdType id_face, faces) {
  63     mesh_size = min(meshSize(id_face), mesh_size);
  64   }
  65   return mesh_size;
  66 }
  67
  68 int CreateHexIbMesh::refine()
  69 {
  70   int N = 0;
  71   m_Octree.resetRefineMarks();
  72   int old_num_cells = m_Octree.getNumCells();
  73   for (int cell = 0; cell < m_Octree.getNumCells(); ++cell) {
  74     if (!m_Octree.hasChildren(cell)) {
  75       double H = m_Octree.getDx(cell);
  76       if (H > 2*m_MeshSize[cell]*m_MinDim && H > 2*m_MinSize) {
  77         m_Octree.markToRefine(cell);
  78         ++N;
  79       }
  80     }
  81   }
  82   m_Octree.refineAll();
  83   m_Faces.insert(old_num_cells, m_Octree.getNumCells() - old_num_cells, QList<vtkIdType>());
  84   m_MeshSize.insert(old_num_cells, m_Octree.getNumCells() - old_num_cells, m_MaxEdgeLength);
  85   for (int cell = old_num_cells; cell < m_Octree.getNumCells(); ++cell) {
  86     foreach (vtkIdType id_face, m_Faces[m_Octree.getParent(cell)]) {
  87       vtkIdType *pts, num_pts;
  88       m_Grid->GetCellPoints(id_face, num_pts, pts);
  89       if (num_pts != 3) {
  90         EG_BUG;
  91       }
  92       QVector<vec3_t> tri(3);
  93       for (int i = 0; i < 3; ++i) {
  94         m_Grid->GetPoint(pts[i], tri[i].data());
  95       }
  96
  97       // compute scale factor to make sure that we have the required minimal number of layers
  98       double H = max(m_Octree.getDx(cell), max(m_Octree.getDy(cell), m_Octree.getDz(cell)));
  99       double h = meshSize(id_face)*m_MinNumLayersWithRequiredResolution;
 100       double scale = 1.0;
 101       if (h < H) {
 102         scale += h/H;
 103       }
 104       if (m_Octree.triangleIntersectsCell(cell, tri, scale)) {
 105         m_Faces[cell].append(id_face);
 106       }
 107     }
 108     m_MeshSize[cell] = min(meshSize(m_Faces[cell]), m_ELSManager.minEdgeLength(m_Octree.getCellCentre(cell)));
 109   }
 110   int num_supercells = 0;
 111   for (int cell = 0; cell < m_Octree.getNumCells(); ++cell) {
 112     if (!m_Octree.hasChildren(cell)) {
 113       ++num_supercells;
 114     }
 115   }
 116   cout << qPrintable(bigIntText(num_supercells)) << " super cells" << endl;
 117   return N;
 118 }
 119
 120 void CreateHexIbMesh::updateMeshSize()
 121 {
 122   // find smallest cell size (on highest level)
 123   double H = 1e99;
 124   for (int cell = 0; cell < m_Octree.getNumCells(); ++cell) {
 125     if (!m_Octree.hasChildren(cell)) {
 126       double l = m_Octree.getLevel(cell);
 127       H = min(H, m_MeshSize[cell]*pow(2.0, l));
 128     }
 129   }
 130
 131   // set discrete cell size according to highest level
 132   for (int cell = 0; cell < m_Octree.getNumCells(); ++cell) {
 133     if (!m_Octree.hasChildren(cell)) {
 134       double l = m_Octree.getLevel(cell);
 135       m_MeshSize[cell] = H/pow(2.0, l);
 136     }
 137   }
 138 }
 139
 140 void CreateHexIbMesh::findInsideCells(MeshPartition &part, QList<vtkIdType> &inside_cells)
 141 {
 142   vtkUnstructuredGrid *grid = part.getGrid();
 143   QVector<bool> is_inside(grid->GetNumberOfCells(), false);
 144
 145   // recompute faces for cells with no children with a minimal scale factor
 146   for (int cell = 0; cell < m_Octree.getNumCells(); ++cell) {
 147     if (m_Octree.getLevel(cell) > 0 && !m_Octree.hasChildren(cell)) {
 148       m_Faces[cell].clear();
 149       foreach (vtkIdType id_face, m_Faces[m_Octree.getParent(cell)]) {
 150         vtkIdType *pts, num_pts;
 151         m_Grid->GetCellPoints(id_face, num_pts, pts);
 152         QVector<vec3_t> tri(3);
 153         for (int i = 0; i < 3; ++i) {
 154           m_Grid->GetPoint(pts[i], tri[i].data());
 155         }
 156         if (m_Octree.triangleIntersectsCell(cell, tri, 1.01)) {
 157           m_Faces[cell].append(id_face);
 158         }
 159       }
 160     }
 161   }
 162
 163   // find closest cell to m_InsidePosition and mark boundary cells
 164   vtkIdType id_start = -1;
 165   {
 166     double min_dist = 1e99;
 167     for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
 168       vec3_t x = cellCentre(grid, id_cell);
 169       double dist = (m_InsidePosition - x).abs();
 170       if (dist < min_dist) {
 171         min_dist = dist;
 172         id_start = id_cell;
 173       }
 174       int cell = m_Octree.findCell(x);
 175       if (m_Faces[cell].size() > 0) {
 176         is_inside[id_cell] = true;
 177       }
 178     }
 179   }
 180   if (id_start == -1) {
 181     EG_BUG;
 182   }
 183   QList<vtkIdType> front_cells;
 184   front_cells << id_start;
 185   is_inside[id_start] = true;
 186   while (front_cells.size() != 0) {
 187     QList<vtkIdType> new_front_cells;
 188     foreach (vtkIdType id_cell, front_cells) {
 189       for (int i = 0; i < part.c2cGSize(id_cell); ++i) {
 190         vtkIdType id_neigh = part.c2cGG(id_cell, i);
 191         if (id_neigh >= 0) {
 192           if (!is_inside[id_neigh]) {
 193             is_inside[id_neigh] = true;
 194             new_front_cells << id_neigh;
 195           }
 196         }
 197       }
 198     }
 199     front_cells = new_front_cells;
 200   }
 201   inside_cells.clear();
 202   for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
 203     if (is_inside[id_cell]) {
 204       inside_cells << id_cell;
 205     }
 206   }
 207 }
 208
 209 void CreateHexIbMesh::operate()
 210 {
 211   DeleteVolumeGrid del_vol;
 212   del_vol();
 213   m_Part.setAllCells();
 214   cout << "refining grid" << endl;
 215   updateNodeInfo();
 216   UpdateDesiredMeshDensity update_mesh_density;
 217   update_mesh_density.readSettings();
 218   update_mesh_density.setRelaxationOff();
 219   update_mesh_density();
 220
 221   QString buffer = GuiMainWindow::pointer()->getXmlSection("engrid/surface/settings").replace("\n", " ");
 222   if (!buffer.isEmpty()) {
 223     QTextStream in(&buffer, QIODevice::ReadOnly);
 224     in >> m_MaxEdgeLength;
 225     in >> m_MinEdgeLength;
 226     in >> m_GrowthFactor;
 227   } else {
 228     m_MaxEdgeLength = 1000.0;
 229     m_MinEdgeLength = 0.0;
 230     m_GrowthFactor = 1.5;
 231   }
 232   m_ELSManager.read();
 233
 234   m_Octree.setSmoothTransitionOn();
 235   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 236     if (m_Grid->GetCellType(id_cell) != VTK_TRIANGLE) {
 237       EG_BUG;
 238     }
 239   }
 240   {
 241     double bounds[6];
 242     m_Grid->GetBounds(bounds);
 243     vec3_t x1(bounds[0], bounds[2], bounds[4]);
 244     vec3_t x2(bounds[1], bounds[3], bounds[5]);
 245     vec3_t xc = 0.5*(x1 + x2);
 246     vec3_t Dx = x2 - xc;
 247     double dx_max = max(Dx[0], max(Dx[1], Dx[2]));
 248     Dx = vec3_t(dx_max, dx_max, dx_max);
 249     x1 = xc - 2*Dx;
 250     x2 = xc + 2*Dx;
 251     m_Octree.setBounds(x1, x2);
 252     m_MinSize = 0;
 253   }
 254   m_Faces.resize(1);
 255   m_MeshSize.resize(1);
 256   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 257     m_Faces[0].append(id_cell);
 258   }
 259   m_MeshSize[0] = meshSize(m_Faces[0]);
 260
 261   int N;
 262   do {
 263     N = refine();
 264   } while (N > 0);
 265   updateMeshSize();
 266
 267   cout << "deleting outside super cells" << endl;
 268   EG_VTKSP(vtkUnstructuredGrid, otgrid);
 269   m_Octree.toVtkGridPolyhedral(otgrid, true);
 270   MeshPartition otpart(otgrid);
 271   otpart.setAllCells();
 272   QList<vtkIdType> add_cells;
 273   findInsideCells(otpart, add_cells);
 274   otpart.setCells(add_cells);
 275   EG_VTKDCC(vtkDoubleArray, cl, otgrid, "cell_subres");
 276   double fvcells = 0;
 277   int N_min = 1000000;
 278   int N_max = 0;
 279   foreach (vtkIdType id_cell, add_cells) {
 280     vec3_t x = cellCentre(otgrid, id_cell);
 281     int cell = m_Octree.findCell(x);
 282     double h = m_MeshSize[cell];
 283     cl->SetValue(id_cell, h);
 284     double Ni = m_Octree.getDx(cell)/h;
 285     double Nj = m_Octree.getDy(cell)/h;
 286     double Nk = m_Octree.getDx(cell)/h;
 287     fvcells += Ni*Nj*Nk;
 288     N_min = min(N_min, min(int(Ni), min(int(Nj), int(Nk))));
 289     N_max = max(N_max, max(int(Ni), max(int(Nj), int(Nk))));
 290   }
 291   cout << add_cells.size() << " super cells" << endl;
 292   if (fvcells > 1e12) {
 293     cout << fvcells << " finite volume cells!!!!" << endl;
 294   } else {
 295     long long int fvcells_li = fvcells;
 296     cout << qPrintable(bigIntText(fvcells_li)) << " finite volume cells" << endl;
 297   }
 298   cout << "smallest dimension: " << N_min << endl;
 299   cout << "largest dimension: " << N_max << endl;
 300   m_Part.addPartition(otpart);
 301   m_Part.setAllCells();
 302 }