src/libengrid/triangularcadinterface.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 "triangularcadinterface.h"
  22 #include "guimainwindow.h"
  23
  24 vtkIdType TriangularCadInterface::m_LastPindex = 0;
  25
  26 TriangularCadInterface::TriangularCadInterface()
  27 {
  28   m_BGrid = vtkUnstructuredGrid::New();
  29   this->setGrid(m_BGrid);
  30   m_CritDistance = 0.1;
  31   setName("triangulated geometry CAD interface");
  32 }
  33
  34 void TriangularCadInterface::computeSurfaceCurvature()
  35 {
  36   m_Radius.fill(1e99, m_BGrid->GetNumberOfPoints());
  37   for (vtkIdType id_node = 0; id_node < m_BGrid->GetNumberOfPoints(); ++id_node) {
  38     vec3_t x1;
  39     m_BGrid->GetPoint(id_node, x1.data());
  40     for (int i = 0; i < m_BPart.n2nGSize(id_node); ++i) {
  41       vtkIdType id_neigh = m_BPart.n2nGG(id_node, i);
  42       double scal_prod = max(-1.0, min(1.0, m_NodeNormals[id_node]*m_NodeNormals[id_neigh]));
  43       double alpha = max(1e-3, acos(scal_prod));
  44       if (alpha > 1e-3) {
  45         vec3_t x2;
  46         m_BGrid->GetPoint(id_neigh, x2.data());
  47         double a = (x1 - x2).abs();
  48         m_Radius[id_node] = min(m_Radius[id_node], a/alpha);
  49       }
  50     }
  51   }
  52
  53   // compute weighted (distance) average of radii
  54   QVector<double> R_new(m_BGrid->GetNumberOfPoints(), 1e99);
  55   for (vtkIdType id_node = 0; id_node < m_BGrid->GetNumberOfPoints(); ++id_node) {
  56     vec3_t x1;
  57     m_BGrid->GetPoint(id_node, x1.data());
  58     int N = m_BPart.n2nGSize(id_node);
  59     QVector<double> L(N);
  60     QVector<double> R(N, -1);
  61     double Lmax = 0;
  62     bool average = false;
  63     for (int i = 0; i < N; ++i) {
  64       vtkIdType id_neigh = m_BPart.n2nGG(id_node, i);
  65       vec3_t x2;
  66       m_BGrid->GetPoint(id_neigh, x2.data());
  67       L[i] = (x2 - x1).abs();
  68       if (m_Radius[id_neigh] < 1e90 && L[i] > 0) {
  69         R[i] = m_Radius[id_neigh];
  70         Lmax = max(Lmax, L[i]);
  71         average = true;
  72       }
  73     }
  74     if (average) {
  75       R_new[id_node] = 0;
  76       double total_weight = 0;
  77       for (int i = 0; i < N; ++i) {
  78         if (R[i] > 0) {
  79           R_new[id_node] += Lmax/L[i] * R[i];
  80           total_weight += Lmax/L[i];
  81           //R_new[id_node] += R[i];
  82           //total_weight += 1.0;
  83         }
  84       }
  85       R_new[id_node] /= total_weight;
  86     }
  87   }
  88   m_Radius = R_new;
  89 }
  90
  91 void TriangularCadInterface::updateBackgroundGridInfo()
  92 {
  93   getAllCells(m_Cells, m_BGrid);
  94   getNodesFromCells(m_Cells, m_Nodes, m_BGrid);
  95   setBoundaryCodes(GuiMainWindow::pointer()->getAllBoundaryCodes());
  96   setAllCells();
  97   readVMD();
  98   updatePotentialSnapPoints();
  99   l2l_t  c2c   = getPartC2C();
 100   g2l_t _cells = getPartLocalCells();
 101   QVector<int> m_LNodes(m_Nodes.size());
 102   for (int i = 0; i < m_LNodes.size(); ++i) {
 103     m_LNodes[i] = i;
 104   }
 105   createNodeToNode(m_Cells, m_Nodes, m_LNodes, m_N2N, m_BGrid);
 106
 107   // create m_Triangles
 108   m_Triangles.resize(m_BGrid->GetNumberOfCells());
 109   for (vtkIdType id_cell = 0; id_cell < m_BGrid->GetNumberOfCells(); ++id_cell) {
 110     m_Triangles[id_cell] = Triangle(m_BGrid, id_cell);
 111     for (int i = 0; i < 3; i++) {
 112       int i_cell = _cells[id_cell];
 113       if (c2c[i_cell][i] < 0) {
 114         m_Triangles[id_cell].setNeighbourFalse(i);
 115       } else {
 116         m_Triangles[id_cell].setNeighbourTrue(i);
 117       }
 118     }
 119   }
 120
 121   // compute node normals
 122   m_NodeNormals.resize(m_BGrid->GetNumberOfPoints());
 123   for (vtkIdType id_node = 0; id_node < m_BGrid->GetNumberOfPoints(); ++id_node) {
 124     m_NodeNormals[id_node] = vec3_t(0, 0, 0);
 125   }
 126
 127   foreach(Triangle T, m_Triangles) {
 128     double angle_a = GeometryTools::angle(m_BGrid, T.idC(), T.idA(), T.idB());
 129     double angle_b = GeometryTools::angle(m_BGrid, T.idA(), T.idB(), T.idC());
 130     double angle_c = GeometryTools::angle(m_BGrid, T.idB(), T.idC(), T.idA());
 131     if (isnan(angle_a) || isinf(angle_a)) EG_BUG;
 132     if (isnan(angle_b) || isinf(angle_b)) EG_BUG;
 133     if (isnan(angle_c) || isinf(angle_c)) EG_BUG;
 134     if (!checkVector(T.g3())) {
 135       qWarning() << "T.g3 = " << T.g3();
 136       EG_BUG;
 137     }
 138     m_NodeNormals[T.idA()] += angle_a * T.g3();
 139     m_NodeNormals[T.idB()] += angle_b * T.g3();
 140     m_NodeNormals[T.idC()] += angle_c * T.g3();
 141     if (!checkVector(m_NodeNormals[T.idA()])) EG_BUG;
 142     if (!checkVector(m_NodeNormals[T.idB()])) EG_BUG;
 143     if (!checkVector(m_NodeNormals[T.idC()])) EG_BUG;
 144   }
 145   for (vtkIdType id_node = 0; id_node < m_BGrid->GetNumberOfPoints(); ++id_node) {
 146     m_NodeNormals[id_node].normalise();
 147   }
 148   for (int i = 0; i < m_Triangles.size(); ++i) {
 149     Triangle &T = m_Triangles[i];
 150     T.setNormals(m_NodeNormals[T.idA()], m_NodeNormals[T.idB()], m_NodeNormals[T.idC()]);
 151   }
 152
 153   m_BPart.setGrid(m_BGrid);
 154   m_BPart.setAllCells();
 155   computeSurfaceCurvature();
 156 }
 157
 158 void TriangularCadInterface::searchNewTriangle(vec3_t xp, vtkIdType &id_tri, vec3_t &x_proj, vec3_t &r_proj, bool &on_triangle)
 159 {
 160   x_proj = vec3_t(1e99, 1e99, 1e99);
 161   r_proj = vec3_t(0, 0, 0);
 162   vec3_t xi, ri;
 163   double d_min      = 1e99;
 164   bool   x_proj_set = false;
 165   on_triangle = false;
 166   QVector<vtkIdType> candidate_faces;
 167   m_FaceFinder.getCloseFaces(xp, candidate_faces);
 168   candidate_faces.clear();
 169   if (candidate_faces.size() == 0) {
 170     // backup method -- really inefficient!
 171     candidate_faces.resize(m_Triangles.size());
 172     for (vtkIdType id_triangle = 0; id_triangle < m_Triangles.size(); ++id_triangle) {
 173       candidate_faces[id_triangle] = id_triangle;
 174     }
 175   }
 176   foreach (vtkIdType id_triangle, candidate_faces) {
 177     Triangle T = m_Triangles[id_triangle];
 178     double d;
 179     int side;
 180     bool intersects = T.snapOntoTriangle(xp, xi, ri, d, side, true);
 181     if (d >= 1e99) {
 182       EG_BUG;
 183     }
 184     if (d < d_min && (intersects || !on_triangle)) {
 185       x_proj = xi;
 186       x_proj_set = true;
 187       r_proj = ri;
 188       d_min = d;
 189       id_tri = id_triangle;
 190       on_triangle = intersects;
 191     }
 192   }
 193   if (!x_proj_set) {
 194     EG_BUG;
 195   }
 196 }
 197
 198 vec3_t TriangularCadInterface::correctCurvature(vtkIdType proj_triangle, vec3_t x)
 199 {
 200   vec3_t x_corr = x;
 201   if (proj_triangle != -1) {
 202     Triangle &T = m_Triangles[proj_triangle];
 203     vec3_t rx = T.global3DToLocal3D(x);
 204     double w1 = 1.0 - rx[0] - rx[1];
 205     double w2 = rx[0];
 206     double w3 = rx[1];
 207     double k1 = GeometryTools::intersection(x, T.g3(), T.a(), T.nA());
 208     double k2 = GeometryTools::intersection(x, T.g3(), T.b(), T.nB());
 209     double k3 = GeometryTools::intersection(x, T.g3(), T.c(), T.nC());
 210     double S = 0.5;
 211     double k = w1*k1 + w2*k2 + w3*k3;
 212     k -= rx[2];
 213     x_corr = x + S*k*T.g3();
 214     if (!checkVector(x_corr)) {
 215       x_corr = x;
 216     }
 217   }
 218   return x_corr;
 219 }
 220
 221 vec3_t TriangularCadInterface::correctCurvature(vec3_t x)
 222 {
 223   vec3_t x_proj(1e99, 1e99, 1e99);
 224   vec3_t r_proj(0, 0, 0);
 225   vtkIdType proj_triangle;
 226   bool on_triangle = false;
 227   searchNewTriangle(x, proj_triangle, x_proj, r_proj, on_triangle);
 228   if (proj_triangle >= m_Triangles.size()) {
 229     EG_BUG;
 230   }
 231   return correctCurvature(proj_triangle, x)  ;
 232 }
 233
 234 vec3_t TriangularCadInterface::snap(vec3_t x, bool correct_curvature)
 235 {
 236   vec3_t x_proj(1e99, 1e99, 1e99);
 237   vec3_t r_proj(0, 0, 0);
 238
 239   bool on_triangle = false;
 240
 241   vtkIdType proj_triangle;
 242   searchNewTriangle(x, proj_triangle, x_proj, r_proj, on_triangle);
 243   if (proj_triangle >= m_Triangles.size()) {
 244     EG_BUG;
 245   }
 246   Triangle T = m_Triangles[proj_triangle];
 247   vec3_t xi, ri;
 248   double d;
 249   int side;
 250   T.snapOntoTriangle(x, xi, ri, d, side, true);
 251   x_proj = xi;
 252   if (x_proj[0] > 1e98) { // should never happen
 253     EG_BUG;
 254   }
 255   r_proj = ri;
 256   if (!checkVector(x_proj)) {
 257     x_proj = x;
 258     m_Failed = true;
 259   } else {
 260     vec2_t r = T.global3DToLocal2D(x_proj);
 261     double Ra   = m_Radius[T.idA()];
 262     double Rb   = m_Radius[T.idB()];
 263     double Rc   = m_Radius[T.idC()];
 264     double R    = min(Ra, min(Rb, Rc));
 265     double Rmax = max(Ra, max(Rb, Rc));
 266     if (Rmax < 1e90) {
 267        R = Ra + r[0]*(Rb - Ra) + r[1]*(Rc - Ra);
 268     }
 269     m_LastRadius = R;
 270     m_LastNormal = T.g3();
 271     if (correct_curvature) {
 272       vec3_t x_corr = correctCurvature(proj_triangle, x_proj);
 273       if (checkVector(x_corr)) {
 274         x_proj = x_corr;
 275       }
 276     }
 277     m_Failed = false;
 278   }
 279   return x_proj;
 280 }
 281
 282 vtkIdType TriangularCadInterface::getProjTriangle(vtkIdType id_node)
 283 {
 284   EG_VTKDCN(vtkLongArray_t, pi, m_FGrid, "node_pindex");
 285   vtkIdType pindex = pi->GetValue(id_node);
 286   vtkIdType proj_triangle = -1;
 287   if (pindex < 0) {
 288     pindex = m_LastPindex;
 289     pi->SetValue(id_node, pindex);
 290     ++m_LastPindex;
 291     m_Pindex[pindex] = -1;
 292   } else {
 293     proj_triangle = m_Pindex[pindex];
 294   }
 295   return proj_triangle;
 296 }
 297
 298 void TriangularCadInterface::setProjTriangle(vtkIdType id_node, vtkIdType proj_triangle)
 299 {
 300   EG_VTKDCN(vtkLongArray_t, pi, m_FGrid, "node_pindex");
 301   vtkIdType pindex = pi->GetValue(id_node);
 302   if (pindex < 0) {
 303     pindex = m_LastPindex;
 304     pi->SetValue(id_node, pindex);
 305     ++m_LastPindex;
 306   }
 307   m_Pindex[pindex] = proj_triangle;
 308 }
 309
 310 vec3_t TriangularCadInterface::snapNode(vtkIdType id_node, vec3_t x, bool correct_curvature)
 311 {
 312   //return snap(x, correct_curvature);
 313
 314   if (!checkVector(x)) {
 315     EG_BUG;
 316   }
 317
 318   if (id_node < 0) {
 319     EG_BUG;
 320   }
 321   if (id_node >= m_FGrid->GetNumberOfPoints()) {
 322     EG_BUG;
 323   }
 324   vec3_t x_proj(1e99, 1e99, 1e99);
 325   vec3_t r_proj(0, 0, 0);
 326   bool x_proj_set = false;
 327
 328   bool on_triangle = false;
 329
 330   vtkIdType proj_triangle = -1;
 331   if (id_node != -1) {
 332     //proj_triangle = getProjTriangle(id_node);
 333   }
 334
 335   if (proj_triangle == -1) {
 336     searchNewTriangle(x, proj_triangle, x_proj, r_proj, on_triangle);
 337     if (id_node != -1) {
 338       setProjTriangle(id_node, proj_triangle);
 339     }
 340   }
 341   if (proj_triangle >= m_Triangles.size()) {
 342     EG_BUG;
 343   }
 344   Triangle T = m_Triangles[proj_triangle];
 345   vec3_t xi, ri;
 346   double d;
 347   int side;
 348   bool intersects = T.snapOntoTriangle(x, xi, ri, d, side, true);
 349   if (!intersects || (d > m_CritDistance*T.smallestLength())) {
 350     searchNewTriangle(x, proj_triangle, x_proj, r_proj, on_triangle);
 351     T = m_Triangles[proj_triangle];
 352     T.snapOntoTriangle(x, xi, ri, d, side, true);
 353     if (id_node != -1) {
 354       setProjTriangle(id_node, proj_triangle);
 355     }
 356   }
 357   x_proj = xi;
 358   x_proj_set = true;
 359   if (x_proj[0] > 1e98) { // should never happen
 360     EG_BUG;
 361   }
 362   r_proj = ri;
 363   on_triangle = intersects;
 364   if (!checkVector(x_proj)) {
 365     x_proj = x;
 366     m_Failed = true;
 367   } else {
 368     vec2_t r = T.global3DToLocal2D(x_proj);
 369     double Ra   = m_Radius[T.idA()];
 370     double Rb   = m_Radius[T.idB()];
 371     double Rc   = m_Radius[T.idC()];
 372     double R    = min(Ra, min(Rb, Rc));
 373     double Rmax = max(Ra, max(Rb, Rc));
 374     if (Rmax < 1e90) {
 375        R = Ra + r[0]*(Rb - Ra) + r[1]*(Rc - Ra);
 376     }
 377     m_LastRadius = R;
 378     m_LastNormal = T.g3();
 379     if (correct_curvature) {
 380       vec3_t x_corr = correctCurvature(proj_triangle, x_proj);
 381       if (checkVector(x_corr)) {
 382         x_proj = x_corr;
 383       }
 384     }
 385     m_Failed = false;
 386   }
 387   return x_proj;
 388 }
 389