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