src/libengrid/triangle.cpp

   1 //
   2 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   3 // +                                                                      +
   4 // + This file is part of enGrid.                                         +
   5 // +                                                                      +
   6 // + Copyright 2008-2013 enGits GmbH                                      +
   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 #include "triangle.h"
  24 #include "geometrytools.h"
  25 #include "engrid.h"
  26 #include "utilities.h"
  27 #include "egvtkobject.h"
  28
  29 Triangle::Triangle() : EgVtkObject()
  30 {
  31   setDefaults();
  32   setupTriangle();
  33 }
  34
  35 Triangle::Triangle(vec3_t a, vec3_t b, vec3_t c) : EgVtkObject()
  36 {
  37   setDefaults();
  38   m_Xa = a;
  39   m_Xb = b;
  40   m_Xc = c;
  41   setupTriangle();
  42 }
  43
  44 Triangle::Triangle(vtkUnstructuredGrid* grid, vtkIdType id_a, vtkIdType id_b, vtkIdType id_c) : EgVtkObject()
  45 {
  46   setDefaults();
  47   m_IdA = id_a;
  48   m_IdB = id_b;
  49   m_IdC = id_c;
  50   grid->GetPoints()->GetPoint(id_a, m_Xa.data());
  51   grid->GetPoints()->GetPoint(id_b, m_Xb.data());
  52   grid->GetPoints()->GetPoint(id_c, m_Xc.data());
  53   setupTriangle();
  54 }
  55
  56 Triangle::Triangle(vtkUnstructuredGrid* grid, vtkIdType id_cell)  : EgVtkObject()
  57 {
  58   setDefaults();
  59   vtkIdType Npts, *pts;
  60   grid->GetCellPoints(id_cell, Npts, pts);
  61   if (Npts == 3) {
  62     m_IdA = pts[0];
  63     m_IdB = pts[1];
  64     m_IdC = pts[2];
  65     grid->GetPoints()->GetPoint(m_IdA, m_Xa.data());
  66     grid->GetPoints()->GetPoint(m_IdB, m_Xb.data());
  67     grid->GetPoints()->GetPoint(m_IdC, m_Xc.data());
  68     setupTriangle();
  69   } else {
  70     EG_ERR_RETURN("only triangles allowed at the moment");
  71   }
  72 }
  73
  74 void Triangle::setDefaults()
  75 {
  76   m_IdA = 0;
  77   m_IdB = 0;
  78   m_IdC = 0;
  79   m_Xa = vec3_t(0, 0, 0);
  80   m_Xb = vec3_t(0, 1, 0);
  81   m_Xc = vec3_t(0, 0, 1);
  82   m_NormalA = vec3_t(0, 0, 0);
  83   m_NormalB = vec3_t(0, 0, 0);
  84   m_NormalC = vec3_t(0, 0, 0);
  85 }
  86
  87 void Triangle::setupTriangle()
  88 {
  89   m_HasNeighbour.fill(false, 6);
  90
  91   m_G1 = m_Xb - m_Xa;
  92   m_G2 = m_Xc - m_Xa;
  93   m_G3 = m_G1.cross(m_G2);
  94
  95   if(m_G3.abs2() <= 0) {
  96     m_Valid = false;
  97   } else {
  98     m_Valid = true;
  99   }
 100
 101   if(!checkVector(m_G3)) {
 102     qWarning() << "m_G1=" << m_G1;
 103     qWarning() << "m_G2=" << m_G2;
 104     qWarning() << "m_G3=" << m_G3;
 105     EG_BUG;
 106   }
 107
 108   m_A  = 0.5 * m_G3.abs();
 109   if (m_Valid) {
 110     m_G3.normalise();
 111   }
 112
 113   if(!checkVector(m_G3)) {
 114     qWarning() << "m_G1="<<m_G1;
 115     qWarning() << "m_G2="<<m_G2;
 116     qWarning() << "m_G3="<<m_G3;
 117     qWarning() << "m_Xa="<<m_Xa;
 118     qWarning() << "m_Xb="<<m_Xb;
 119     qWarning() << "m_Xc="<<m_Xc;
 120     this->saveTriangle("crash");
 121     EG_BUG;
 122   }
 123
 124   m_G.column(0, m_G1);
 125   m_G.column(1, m_G2);
 126   m_G.column(2, m_G3);
 127   m_GI = m_G.inverse();
 128
 129   m_SmallestLength = (m_Xb - m_Xa).abs();
 130   m_SmallestLength = min(m_SmallestLength, (m_Xc - m_Xb).abs());
 131   m_SmallestLength = min(m_SmallestLength, (m_Xa - m_Xc).abs());
 132
 133   // compute minimal height
 134   double ha = (GeometryTools::projectPointOnEdge(m_Xa, m_Xb, (m_Xc - m_Xb)) - m_Xa).abs();
 135   double hb = (GeometryTools::projectPointOnEdge(m_Xb, m_Xa, (m_Xc - m_Xa)) - m_Xb).abs();
 136   double hc = (GeometryTools::projectPointOnEdge(m_Xc, m_Xa, (m_Xb - m_Xa)) - m_Xc).abs();
 137   m_SmallestHeight = min(ha, min(hb, hc));
 138
 139 }
 140
 141 vec3_t Triangle::local3DToGlobal3D(vec3_t r)
 142 {
 143   return m_Xa + m_G*r;
 144 }
 145
 146 vec3_t Triangle::global3DToLocal3D(vec3_t x)
 147 {
 148   vec3_t tmp = x - m_Xa;
 149   return m_GI*tmp;
 150 }
 151
 152 vec3_t Triangle::local2DToGlobal3D(vec2_t r)
 153 {
 154   return local3DToGlobal3D(vec3_t(r[0], r[1], 0));
 155 }
 156
 157 vec2_t Triangle::global3DToLocal2D(vec3_t x)
 158 {
 159   vec3_t r = global3DToLocal3D(x);
 160   return vec2_t(r[0], r[1]);
 161 }
 162
 163 bool Triangle::snapOntoTriangle(vec3_t xp, vec3_t &xi, vec3_t &ri, double &d, int& side, bool restrict_to_triangle)
 164 {
 165   side = -1;
 166   double scal = (xp - this->m_Xa) * this->m_G3;
 167   vec3_t x1, x2;
 168   if (scal > 0) {
 169     x1 = xp + this->m_G3;
 170     x2 = xp - scal * this->m_G3 - this->m_G3;
 171   } else {
 172     x1 = xp - this->m_G3;
 173     x2 = xp - scal * this->m_G3 + this->m_G3;
 174   }
 175   // (xi,ri) gets set to the intersection of the line with the plane here!
 176   bool intersects_face = GeometryTools::intersectEdgeAndTriangle(this->m_Xa, this->m_Xb, this->m_Xc, x1, x2, xi, ri);
 177   vec3_t xi_free = xi;
 178   if (intersects_face) {
 179     vec3_t dx = xp - this->m_Xa;
 180     d = fabs(dx * this->m_G3);
 181   } else {
 182     double kab = GeometryTools::intersection(this->m_Xa, this->m_Xb - this->m_Xa, xp, this->m_Xb - this->m_Xa);
 183     double kac = GeometryTools::intersection(this->m_Xa, this->m_Xc - this->m_Xa, xp, this->m_Xc - this->m_Xa);
 184     double kbc = GeometryTools::intersection(this->m_Xb, this->m_Xc - this->m_Xb, xp, this->m_Xc - this->m_Xb);
 185
 186     double dab = (this->m_Xa + kab * (this->m_Xb - this->m_Xa) - xp).abs2();
 187     double dac = (this->m_Xa + kac * (this->m_Xc - this->m_Xa) - xp).abs2();
 188     double dbc = (this->m_Xb + kbc * (this->m_Xc - this->m_Xb) - xp).abs2();
 189     double da = (this->m_Xa - xp).abs2();
 190     double db = (this->m_Xb - xp).abs2();
 191     double dc = (this->m_Xc - xp).abs2();
 192
 193     bool set = false;
 194     d = 1e99;//max(max(max(max(max(dab,dac),dbc),da),db),dc);
 195
 196     if (dab < d) {
 197       if ((kab >= 0) && (kab <= 1)) {
 198         xi = this->m_Xa + kab * (this->m_Xb - this->m_Xa);
 199         ri = vec3_t(kab, 0, 0);
 200         d = dab;
 201         set = true;
 202         side = 0;
 203       }
 204     }
 205     if (dbc < d) {
 206       if ((kbc >= 0) && (kbc <= 1)) {
 207         xi = this->m_Xb + kbc * (this->m_Xc - this->m_Xb);
 208         ri = vec3_t(1 - kbc, kbc, 0);
 209         d = dbc;
 210         set = true;
 211         side = 1;
 212       }
 213     }
 214     if (dac < d) {
 215       if ((kac >= 0) && (kac <= 1)) {
 216         xi = this->m_Xa + kac * (this->m_Xc - this->m_Xa);
 217         ri = vec3_t(0, kac, 0);
 218         d = dac;
 219         set = true;
 220         side = 2;
 221       }
 222     }
 223     if (da < d) {
 224       xi = this->m_Xa;
 225       ri = vec3_t(0, 0);
 226       d = da;
 227       set = true;
 228       side = 3;
 229     }
 230     if (db < d) {
 231       xi = this->m_Xb;
 232       ri = vec3_t(1, 0);
 233       d = db;
 234       set = true;
 235       side = 4;
 236     }
 237     if (dc < d) {
 238       xi = this->m_Xc;
 239       ri = vec3_t(0, 1);
 240       d = dc;
 241       set = true;
 242       side = 5;
 243     }
 244     if (!set) {
 245       EG_BUG;
 246     }
 247     d = sqrt(d);
 248   }
 249   if (!intersects_face && !restrict_to_triangle) {
 250     xi = xi_free;
 251   }
 252   if (xi[0] > 1e98) { // should never happen
 253     EG_BUG;
 254   }
 255   return intersects_face;
 256 }
 257
 258 void Triangle::saveTriangle(QString filename)
 259 {
 260   int N_cells = 1;
 261   int N_points = 3;
 262
 263   EG_VTKSP(vtkUnstructuredGrid, triangle_grid);
 264   allocateGrid(triangle_grid , N_cells, N_points);
 265
 266   vtkIdType node_count = 0;
 267   int cell_count = 0;
 268
 269   vtkIdType pts[3];
 270   triangle_grid->GetPoints()->SetPoint(node_count, m_Xa.data()); pts[0]=node_count; node_count++;
 271   triangle_grid->GetPoints()->SetPoint(node_count, m_Xb.data()); pts[1]=node_count; node_count++;
 272   triangle_grid->GetPoints()->SetPoint(node_count, m_Xc.data()); pts[2]=node_count; node_count++;
 273
 274   triangle_grid->InsertNextCell(VTK_TRIANGLE,3,pts);cell_count++;
 275
 276   saveGrid(triangle_grid, filename+"_triangle_grid");
 277 }
 278
 279 void Triangle::setNormals(vec3_t na, vec3_t nb, vec3_t nc)
 280 {
 281   m_NormalA = na;
 282   m_NormalB = nb;
 283   m_NormalC = nc;
 284   // compute normal vectors in local coordinate system
 285   m_RNormalA = global3DToLocal3D(a() + nA());
 286   m_RNormalB = global3DToLocal3D(a() + nB());
 287   m_RNormalC = global3DToLocal3D(a() + nC());
 288 }