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