src/surfaceprojection.h

   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 #ifndef SURFACEPROJECTION_H
  24 #define SURFACEPROJECTION_H
  25
  26 class SurfaceProjection;
  27
  28 #include "egvtkobject.h"
  29 #include "octree.h"
  30 #include "guimainwindow.h"
  31
  32 class SurfaceProjection : public EgVtkObject
  33 {
  34
  35 private: // data-types
  36
  37   struct Triangle
  38   {
  39     vtkIdType id_a, id_b, id_c;
  40     vec3_t a, b, c;
  41     vec3_t g1, g2, g3;
  42     mat3_t G, GI;
  43     double A;
  44     double smallest_length;
  45   };
  46
  47   struct Edge
  48   {
  49     vec3_t a, b;
  50     vec3_t v;
  51     vec3_t na, nb;
  52     double La, Lb;
  53   };
  54
  55 private: // attributes
  56
  57   vtkUnstructuredGrid*   m_BGrid;
  58   QVector<vtkIdType>     m_ProjTriangles;
  59   vtkUnstructuredGrid*   m_FGrid;
  60   QVector<double>        m_EdgeLength;
  61   QVector<vtkIdType>     m_Cells;
  62   QVector<vtkIdType>     m_Nodes;
  63   QVector<vec3_t>        m_NodeNormals;
  64   QVector<Triangle>      m_Triangles;
  65   QVector<QVector<int> > m_N2N;
  66   Octree                 m_OTGrid;
  67   QVector<double>        m_G;
  68   QVector<bool>          m_GSet;
  69   double                 m_Relax;
  70   double                 m_DistExp;
  71   double                 m_DistWeight;
  72   double                 m_DirWeight;
  73   double                 m_DirExp;
  74   double                 m_WeightOffset;
  75   double                 m_MinOTLength;
  76   int                    m_MaxOTCells;
  77   double                 m_Length;
  78   int                    m_MaxIter;
  79   double                 m_ConvLimit;
  80   double                 m_RadiusFactor;
  81   bool                   m_UseLevelSet;
  82   int                    m_NumDirect;
  83   int                    m_NumFull;
  84
  85 private: // methods
  86
  87   template <class C>
  88   void setBackgroundGrid_setupGrid(vtkUnstructuredGrid* grid, const C& cells);
  89   void setBackgroundGrid_initOctree();
  90   void setBackgroundGrid_refineFromNodes();
  91   void setBackgroundGrid_refineFromEdges();
  92   void setBackgroundGrid_refineFromFaces();
  93   void setBackgroundGrid_initLevelSet();
  94   void setBackgroundGrid_computeLevelSet();
  95
  96   vec3_t calcGradG(vec3_t x);
  97   double calcG(vec3_t x);
  98
  99   vec3_t projectWithLevelSet(vec3_t x);
 100
 101   bool   projectOnTriangle(vec3_t xp, vec3_t &xi, vec3_t &ri, double &d, const Triangle& T);
 102   vec3_t projectWithGeometry(vec3_t x, vtkIdType id_node);
 103   vec3_t correctCurvature(int i_tri, vec3_t r);
 104
 105 public: // methods
 106
 107   SurfaceProjection();
 108
 109   template <class C>
 110   void setBackgroundGrid(vtkUnstructuredGrid* grid, const C& cells);
 111
 112   void setForegroundGrid(vtkUnstructuredGrid* grid) {m_FGrid = grid; }
 113
 114   vec3_t project(vec3_t x, vtkIdType id_node = -1);
 115   int getNumOctreeCells() { return m_OTGrid.getNumCells(); }
 116   void writeOctree(QString file_name);
 117   bool usesLevelSet() { return m_UseLevelSet; };
 118
 119   int getNumDirectProjections() { return m_NumDirect; }
 120   int getNumFullSearches() { return m_NumFull; }
 121
 122 };
 123
 124 template <class C>
 125 void SurfaceProjection::setBackgroundGrid(vtkUnstructuredGrid* grid, const C& cells)
 126 {
 127   setBackgroundGrid_setupGrid(grid, cells);
 128   if (m_UseLevelSet) {
 129     setBackgroundGrid_initOctree();
 130     setBackgroundGrid_refineFromNodes();
 131     setBackgroundGrid_refineFromEdges();
 132     setBackgroundGrid_refineFromFaces();
 133     setBackgroundGrid_computeLevelSet();
 134   }
 135 }
 136
 137 template <class C>
 138 void SurfaceProjection::setBackgroundGrid_setupGrid(vtkUnstructuredGrid* grid, const C& cells)
 139 {
 140   QVector<vtkIdType> nodes;
 141   getNodesFromCells(cells, nodes, grid);
 142   allocateGrid(m_BGrid, cells.size(), nodes.size());
 143   QVector<vtkIdType> _nodes(grid->GetNumberOfPoints());
 144   vtkIdType id_new_node = 0;
 145   foreach (vtkIdType id_node, nodes) {
 146     vec3_t x;
 147     grid->GetPoints()->GetPoint(id_node, x.data());
 148     m_BGrid->GetPoints()->SetPoint(id_new_node, x.data());
 149     _nodes[id_node] = id_new_node;
 150     ++id_new_node;
 151   }
 152   foreach (vtkIdType id_cell, cells) {
 153     vtkIdType N_pts, *pts;
 154     vtkIdType type_cell = grid->GetCellType(id_cell);
 155     grid->GetCellPoints(id_cell, N_pts, pts);
 156     vtkIdType new_pts[N_pts];
 157     for (int i = 0; i < N_pts; ++i) {
 158       new_pts[i] = _nodes[pts[i]];
 159     }
 160     m_BGrid->InsertNextCell(type_cell, N_pts, new_pts);
 161   }
 162   getAllCells(m_Cells, m_BGrid);
 163   getNodesFromCells(m_Cells, m_Nodes, m_BGrid);
 164   QVector<int> m_LNodes(m_Nodes.size());
 165   for (int i = 0; i < m_LNodes.size(); ++i) {
 166     m_LNodes[i] = i;
 167   }
 168   createNodeToNode(m_Cells, m_Nodes, m_LNodes, m_N2N, m_BGrid);
 169   m_EdgeLength.fill(1e99, m_BGrid->GetNumberOfPoints());
 170   foreach (vtkIdType id_node, m_Nodes) {
 171     vec3_t x;
 172     m_BGrid->GetPoints()->GetPoint(id_node, x.data());
 173     foreach (vtkIdType id_neigh, m_N2N[id_node]) {
 174       vec3_t xn;
 175       m_BGrid->GetPoints()->GetPoint(id_neigh, xn.data());
 176       m_EdgeLength[id_node] = min(m_EdgeLength[id_node], (x-xn).abs());
 177     }
 178     if (m_N2N[id_node].size() < 2) {
 179       EG_BUG;
 180     }
 181   }
 182   // create m_Triangles
 183   m_Triangles.resize(m_BGrid->GetNumberOfCells());
 184   for (vtkIdType id_cell = 0; id_cell < m_BGrid->GetNumberOfCells(); ++id_cell) {
 185     vtkIdType Npts, *pts;
 186     m_BGrid->GetCellPoints(id_cell, Npts, pts);
 187     if (Npts == 3) {
 188       m_BGrid->GetPoints()->GetPoint(pts[0], m_Triangles[id_cell].a.data());
 189       m_BGrid->GetPoints()->GetPoint(pts[1], m_Triangles[id_cell].b.data());
 190       m_BGrid->GetPoints()->GetPoint(pts[2], m_Triangles[id_cell].c.data());
 191       m_Triangles[id_cell].id_a = pts[0];
 192       m_Triangles[id_cell].id_b = pts[1];
 193       m_Triangles[id_cell].id_c = pts[2];
 194       m_Triangles[id_cell].g1 = m_Triangles[id_cell].b - m_Triangles[id_cell].a;
 195       m_Triangles[id_cell].g2 = m_Triangles[id_cell].c - m_Triangles[id_cell].a;
 196       m_Triangles[id_cell].g3 = m_Triangles[id_cell].g1.cross(m_Triangles[id_cell].g2);
 197       m_Triangles[id_cell].A  = 0.5*m_Triangles[id_cell].g3.abs();
 198       m_Triangles[id_cell].g3.normalise();
 199       m_Triangles[id_cell].G.column(0, m_Triangles[id_cell].g1);
 200       m_Triangles[id_cell].G.column(1, m_Triangles[id_cell].g2);
 201       m_Triangles[id_cell].G.column(2, m_Triangles[id_cell].g3);
 202       m_Triangles[id_cell].GI = m_Triangles[id_cell].G.inverse();
 203       m_Triangles[id_cell].smallest_length = (m_Triangles[id_cell].b - m_Triangles[id_cell].a).abs();
 204       m_Triangles[id_cell].smallest_length = min(m_Triangles[id_cell].smallest_length, (m_Triangles[id_cell].c - m_Triangles[id_cell].b).abs());
 205       m_Triangles[id_cell].smallest_length = min(m_Triangles[id_cell].smallest_length, (m_Triangles[id_cell].a - m_Triangles[id_cell].c).abs());
 206     } else {
 207       EG_ERR_RETURN("only triangles allowed at the moment");
 208     }
 209   }
 210
 211   // compute node normals
 212   m_NodeNormals.resize(m_BGrid->GetNumberOfPoints());
 213   for (vtkIdType id_node = 0; id_node < m_BGrid->GetNumberOfPoints(); ++id_node) {
 214     m_NodeNormals[id_node] = vec3_t(0,0,0);
 215   }
 216   foreach (Triangle T, m_Triangles) {
 217     m_NodeNormals[T.id_a] += T.A*T.g3;
 218     m_NodeNormals[T.id_b] += T.A*T.g3;
 219     m_NodeNormals[T.id_c] += T.A*T.g3;
 220   }
 221   for (vtkIdType id_node = 0; id_node < m_BGrid->GetNumberOfPoints(); ++id_node) {
 222     m_NodeNormals[id_node].normalise();
 223   }
 224
 225   // compute maximum angle per node
 226   QVector<double> min_cos(m_BGrid->GetNumberOfPoints(), 1.0);
 227   foreach (Triangle T, m_Triangles) {
 228     double cosa = T.g3*m_NodeNormals[T.id_a];
 229     double cosb = T.g3*m_NodeNormals[T.id_b];
 230     double cosc = T.g3*m_NodeNormals[T.id_c];
 231     min_cos[T.id_a] = min(cosa, min_cos[T.id_a]);
 232     min_cos[T.id_b] = min(cosb, min_cos[T.id_b]);
 233     min_cos[T.id_c] = min(cosc, min_cos[T.id_c]);
 234   }
 235   for (vtkIdType id_node = 0; id_node < m_BGrid->GetNumberOfPoints(); ++id_node) {
 236     double s = sqrt(1.0 - sqr(min(1 - 1e-20, min_cos[id_node])));
 237     m_EdgeLength[id_node] *= m_RadiusFactor*min_cos[id_node]/s;
 238   }
 239
 240 }
 241
 242 #endif // SURFACEPROJECTION_H