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