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