src/libengrid/facefinder.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 "facefinder.h"
  22 #include "engrid.h"
  23 #include "triangle.h"
  24
  25 FaceFinder::FaceFinder()
  26 {
  27   m_MinSize = 1.0;
  28   m_MaxFaces = 10;
  29   m_CollectedPointsIncrement = 10000;
  30   m_CollectedPoints.resize(10000);
  31   getSet("surface meshing", "use improved face search", false, m_UseImprovedFaceSearch);
  32 }
  33
  34 void FaceFinder::setGrid(vtkUnstructuredGrid *grid)
  35 {
  36   m_Grid = grid;
  37   m_Triangles.resize(m_Grid->GetNumberOfCells());
  38   m_Centres.resize(m_Grid->GetNumberOfCells());
  39   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
  40     if (m_Grid->GetCellType(id_cell) != VTK_TRIANGLE) {
  41       EG_BUG;
  42     }
  43     m_Triangles[id_cell] = Triangle(m_Grid, id_cell);
  44     m_Centres[id_cell] = cellCentre(m_Grid, id_cell);
  45   }
  46   {
  47     double bounds[6];
  48     m_Grid->GetBounds(bounds);
  49     vec3_t x1(bounds[0], bounds[2], bounds[4]);
  50     vec3_t x2(bounds[1], bounds[3], bounds[5]);
  51     vec3_t xc = 0.5*(x1 + x2);
  52     vec3_t Dx = x2 - xc;
  53     double dx_max = max(Dx[0], max(Dx[1], Dx[2]));
  54     Dx = vec3_t(dx_max, dx_max, dx_max);
  55     x1 = xc - 2*Dx;
  56     x2 = xc + 2*Dx;
  57     m_Octree.setBounds(x1, x2);
  58     if (m_UseImprovedFaceSearch) {
  59       m_MinSize = 0.005*Dx[0];
  60     } else {
  61       m_MinSize = 0.0001*Dx[0];
  62     }
  63   }
  64   m_Faces.resize(1);
  65   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
  66     m_Faces[0].append(id_cell);
  67   }
  68   if (m_UseImprovedFaceSearch) {
  69     m_CritLengthOctreeCell.fill(EG_LARGE_REAL, 1);
  70     foreach (vtkIdType id_cell, m_Faces[0]) {
  71       m_CritLengthOctreeCell[0] = min(m_CritLengthOctreeCell[0], calcCritLength(id_cell));
  72     }
  73     calcCritLengthForAllNodes();
  74   } else {
  75     m_CritLengthOctreeCell.fill(0, 1);
  76     foreach (vtkIdType id_cell, m_Faces[0]) {
  77       m_CritLengthOctreeCell[0] = max(m_CritLengthOctreeCell[0], calcCritLength(id_cell));
  78     }
  79   }
  80
  81   int N;
  82   do {
  83     N = refine();
  84   } while (N > 0);
  85
  86   int max_num_faces = 0;
  87   double ave_num_faces = 0;
  88   int num_buckets = 0;
  89   for (int i_cell = 0; i_cell < m_Octree.getNumCells(); ++i_cell) {
  90     if (!m_Octree.hasChildren(i_cell)) {
  91       int cell = i_cell;
  92       while (m_Faces[cell].size() == 0) {
  93         cell = m_Octree.getParent(cell);
  94         if (cell == -1) {
  95           EG_BUG;
  96         }
  97       }
  98       max_num_faces = max(m_Faces[cell].size(), max_num_faces);
  99       ave_num_faces += double(m_Faces[cell].size());
 100       ++num_buckets;
 101     }
 102   }
 103
 104   cout << "refining octree for FaceFinder ..." << endl;
 105   cout << "  minimal cell size: " << m_MinSize << endl;
 106   cout << "  maximal number of faces: " << max_num_faces << endl;
 107   cout << "  average number of faces: " << ave_num_faces/num_buckets << endl;
 108   cout << "  number of buckets: " << num_buckets << endl;
 109   //cout << "  largest critical length: " << m_CritLength[0] << endl;
 110
 111   //EG_VTKSP(vtkUnstructuredGrid, otg);
 112   //m_Octree.toVtkGrid(otg);
 113   //writeGrid(otg, "facefinder");
 114 }
 115
 116 double FaceFinder::calcCritLength(vtkIdType id_cell)
 117 {
 118   if (m_UseImprovedFaceSearch) {
 119     return max(m_MinSize, 10*m_Triangles[id_cell].smallestLength());
 120   }
 121   QVector<vec3_t> x;
 122   EG_GET_CELL(id_cell, m_Grid);
 123   x.resize(num_pts + 1);
 124   for (int i = 0; i < num_pts; ++i) {
 125     m_Grid->GetPoint(pts[i], x[i].data());
 126   }
 127   x[num_pts] = x[0];
 128   double L = 0;
 129   for (int i = 0; i < num_pts; ++i) {
 130     L = max(L, (x[i]-x[i+1]).abs());
 131   }
 132   return L;
 133 }
 134
 135 void FaceFinder::calcCritLengthForAllNodes()
 136 {
 137   m_CritLengthNode.fill(EG_LARGE_REAL, m_Grid->GetNumberOfPoints());
 138   for (vtkIdType id_face = 0; id_face < m_Grid->GetNumberOfCells(); ++id_face) {
 139     double l_crit = calcCritLength(id_face);
 140     EG_GET_CELL(id_face, m_Grid);
 141     for (int i = 0; i < num_pts; ++i) {
 142       m_CritLengthNode[pts[i]] = min(m_CritLengthNode[pts[i]], l_crit);
 143     }
 144   }
 145 }
 146
 147 void FaceFinder::getPointsOfFace(vtkIdType id_face)
 148 {
 149   m_NumCollectedPoints = 3;
 150   Triangle T = m_Triangles[id_face];
 151   m_CollectedPoints[0] = T.a();
 152   m_CollectedPoints[1] = T.b();
 153   m_CollectedPoints[2] = T.c();
 154   QVector<vec3_t> x_tri(4);
 155   QVector<vtkIdType> id_tri(4);
 156   x_tri[0] = T.a();
 157   x_tri[1] = T.b();
 158   x_tri[2] = T.c();
 159   x_tri[3] = T.a();
 160   id_tri[0] = T.idA();
 161   id_tri[1] = T.idB();
 162   id_tri[2] = T.idC();
 163   id_tri[3] = T.idA();
 164   for (int i = 0; i < 3; ++i) {
 165     vec3_t x1 = x_tri[i];
 166     vec3_t x2 = x_tri[i+1];
 167     double l  = (x2 - x1).abs();
 168     double h = 0.5*(m_CritLengthNode[id_tri[i]] + m_CritLengthNode[id_tri[i+1]]);
 169     int N = int(l/h) + 1;
 170     vec3_t dx = (1.0/N)*(x2 - x1);
 171     vec3_t x = x1 + dx;
 172     for (int j = 1; j < N; ++j) {
 173       if (m_NumCollectedPoints >= m_CollectedPoints.size()) {
 174         m_CollectedPoints.resize(m_CollectedPoints.size() + m_CollectedPointsIncrement);
 175       }
 176       m_CollectedPoints[m_NumCollectedPoints] = x;
 177       ++m_NumCollectedPoints;
 178       x += dx;
 179     }
 180   }
 181 }
 182
 183 int FaceFinder::refine()
 184 {
 185   int N = m_Octree.getNumCells();
 186   m_Octree.setSmoothTransitionOff();
 187   int N_new = 0;
 188   m_Octree.resetRefineMarks();
 189   for (int cell = 0; cell < m_Octree.getNumCells(); ++cell) {
 190     double dx = m_Octree.getDx(cell);
 191     if (!m_Octree.hasChildren(cell)  && m_Faces[cell].size() > m_MaxFaces  && dx > 2*m_MinSize && dx > 2*m_CritLengthOctreeCell[cell]) {
 192       m_Octree.markToRefine(cell);
 193       ++N_new;
 194     }
 195   }
 196   int N_mrk = 0;
 197   for (int cell = 0; cell < m_Octree.getNumCells(); ++cell) {
 198     if (m_Octree.markedForRefine(cell)) {
 199       ++N_mrk;
 200     }
 201   }
 202   N_new *= 8;
 203   m_Octree.refineAll();
 204   if (m_Octree.getNumCells() != N + N_new) {
 205     EG_BUG;
 206   }
 207   m_Faces.insert(N, N_new, QList<vtkIdType>());
 208   if (m_UseImprovedFaceSearch) {
 209     m_CritLengthOctreeCell.insert(N, N_new, 1e99);
 210   } else {
 211     m_CritLengthOctreeCell.insert(N, N_new, 0);
 212   }
 213   int N_min = m_Grid->GetNumberOfCells();
 214   int N_max = 0;
 215   int N_ave = 0;
 216   for (int cell = N; cell < m_Octree.getNumCells(); ++cell) {
 217     m_Timer << "  " << m_Octree.getNumCells() << " octree cells" << Timer::endl;
 218     if (m_Octree.getNumCells() != N + N_new) {
 219       EG_BUG;
 220     }
 221     int parent = m_Octree.getParent(cell);
 222     if (parent < 0) {
 223       EG_BUG;
 224     }
 225     foreach (vtkIdType id_cell, m_Faces[parent]) {
 226       if (m_UseImprovedFaceSearch) {
 227         getPointsOfFace(id_cell);
 228         for (int i = 0; i < m_NumCollectedPoints; ++i) {
 229           if (m_Octree.isInsideCell(cell, m_CollectedPoints[i], 2.0)) {
 230             m_Faces[cell].append(id_cell);
 231             m_CritLengthOctreeCell[cell] = min(m_CritLengthOctreeCell[cell], calcCritLength(id_cell));
 232             break;
 233           }
 234         }
 235       } else {
 236         for (int node = 0; node < 8; ++node) {
 237           vec3_t x = m_Octree.getNodePosition(cell, node);
 238           double d;
 239           d = (x - m_Centres[id_cell]).abs();
 240           if (d < m_Octree.getDx(cell)) {
 241             m_Faces[cell].append(id_cell);
 242             m_CritLengthOctreeCell[cell] = max(m_CritLengthOctreeCell[cell], calcCritLength(id_cell));
 243             break;
 244           }
 245         }
 246       }
 247     }
 248     N_min = min(N_min, m_Faces[cell].size());
 249     N_max = max(N_max, m_Faces[cell].size());
 250     N_ave += m_Faces[cell].size();
 251   }
 252   int N_cells = 0;
 253   int N_faces = 0;
 254   for (int cell = 0; cell < m_Octree.getNumCells(); ++cell) {
 255     if (!m_Octree.hasChildren(cell)) {
 256       ++N_cells;
 257       N_faces += m_Faces[cell].size();
 258     }
 259   }
 260   if (N_faces == 0) {
 261     EG_BUG;
 262   }
 263   return N_new;
 264 }
 265
 266 void FaceFinder::getCloseFaces(vec3_t x, QVector<vtkIdType> &faces)
 267 {
 268   if (m_Octree.isInsideBounds(x)) {
 269     int cell = m_Octree.findCell(x);
 270     if (cell < 0) {
 271       EG_BUG;
 272     }
 273     if (m_Octree.hasChildren(cell)) {
 274       EG_BUG;
 275     }
 276     while (m_Faces[cell].size() == 0 && m_Octree.getParent(cell) >= 0) {
 277       cell = m_Octree.getParent(cell);
 278     }
 279     faces.resize(m_Faces[cell].size());
 280     qCopy(m_Faces[cell].begin(), m_Faces[cell].end(), faces.begin());
 281   } else {
 282     faces.clear();
 283   }
 284 }
 285
 286 vtkIdType FaceFinder::getClosestFace(vec3_t x, double &L_min)
 287 {
 288   QVector<vtkIdType> faces;
 289   getCloseFaces(x, faces);
 290   vtkIdType id_close = -1;
 291   foreach (vtkIdType id_face, faces) {
 292     vec3_t xi, ri;
 293     int side;
 294     double L;
 295     m_Triangles[id_face].snapOntoTriangle(x, xi, ri, L, side, true);
 296     if (L < L_min) {
 297       L_min = L;
 298       id_close = id_face;
 299     }
 300   }
 301   return id_close;
 302 }
 303