src/triangle.cpp

   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 #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   setDefaults();
  31   setupTriangle();
  32 }
  33
  34 Triangle::Triangle(vec3_t a_a, vec3_t a_b, vec3_t a_c) : EgVtkObject() {
  35   setDefaults();
  36   m_a = a_a;
  37   m_b = a_b;
  38   m_c = a_c;
  39   setupTriangle();
  40 }
  41
  42 Triangle::Triangle(vtkUnstructuredGrid* a_grid, vtkIdType a_id_a, vtkIdType a_id_b, vtkIdType a_id_c) : EgVtkObject() {
  43   setDefaults();
  44   m_id_a = a_id_a;
  45   m_id_b = a_id_b;
  46   m_id_c = a_id_c;
  47   a_grid->GetPoints()->GetPoint(m_id_a, m_a.data());
  48   a_grid->GetPoints()->GetPoint(m_id_b, m_b.data());
  49   a_grid->GetPoints()->GetPoint(m_id_c, m_c.data());
  50   setupTriangle();
  51 }
  52
  53 Triangle::Triangle(vtkUnstructuredGrid* a_grid, vtkIdType a_id_cell)  : EgVtkObject() {
  54   setDefaults();
  55   vtkIdType Npts, *pts;
  56   a_grid->GetCellPoints(a_id_cell, Npts, pts);
  57   if (Npts == 3) {
  58     m_id_a = pts[0];
  59     m_id_b = pts[1];
  60     m_id_c = pts[2];
  61     a_grid->GetPoints()->GetPoint(m_id_a, m_a.data());
  62     a_grid->GetPoints()->GetPoint(m_id_b, m_b.data());
  63     a_grid->GetPoints()->GetPoint(m_id_c, m_c.data());
  64     setupTriangle();
  65   } else {
  66     EG_ERR_RETURN("only triangles allowed at the moment");
  67   }
  68 }
  69
  70 void Triangle::setDefaults()
  71 {
  72   m_id_a = 0;
  73   m_id_b = 0;
  74   m_id_c = 0;
  75   m_a = vec3_t(0, 0, 0);
  76   m_b = vec3_t(0, 1, 0);
  77   m_c = vec3_t(0, 0, 1);
  78   m_Normal_a = vec3_t(0, 0, 0);
  79   m_Normal_b = vec3_t(0, 0, 0);
  80   m_Normal_c = vec3_t(0, 0, 0);
  81 }
  82
  83 void Triangle::setupTriangle() {
  84   m_has_neighbour.resize(6);
  85   m_has_neighbour[0] = false;
  86   m_has_neighbour[1] = false;
  87   m_has_neighbour[2] = false;
  88   m_has_neighbour[3] = false;
  89   m_has_neighbour[4] = false;
  90   m_has_neighbour[5] = false;
  91
  92   m_g1 = m_b - m_a;
  93   m_g2 = m_c - m_a;
  94   m_g3 = m_g1.cross(m_g2);
  95
  96   if(m_g3.abs2()<=0) {
  97     m_Valid = false;
  98   }
  99   else {
 100     m_Valid = true;
 101   }
 102
 103   if(!checkVector(m_g3)) {
 104     qWarning()<<"m_g1="<<m_g1;
 105     qWarning()<<"m_g2="<<m_g2;
 106     qWarning()<<"m_g3="<<m_g3;
 107     vec3_t foo = vec3_t(0,0,0);
 108     vec3_t bar = vec3_t(0,0,0);
 109     qWarning()<<"foo.cross(bar)="<<foo.cross(bar);
 110     EG_BUG;
 111   }
 112
 113   m_A  = 0.5 * m_g3.abs();
 114   if(m_Valid) m_g3.normalise();
 115
 116   if(!checkVector(m_g3)) {
 117     qWarning()<<"m_g1="<<m_g1;
 118     qWarning()<<"m_g2="<<m_g2;
 119     qWarning()<<"m_g3="<<m_g3;
 120     qWarning()<<"m_a="<<m_a;
 121     qWarning()<<"m_b="<<m_b;
 122     qWarning()<<"m_c="<<m_c;
 123     this->saveTriangle("crash");
 124     EG_BUG;
 125   }
 126
 127   m_G.column(0, m_g1);
 128   m_G.column(1, m_g2);
 129   m_G.column(2, m_g3);
 130   m_GI = m_G.inverse();
 131
 132   m_smallest_length = (m_b - m_a).abs();
 133   m_smallest_length = min(m_smallest_length, (m_c - m_b).abs());
 134   m_smallest_length = min(m_smallest_length, (m_a - m_c).abs());
 135 }
 136
 137 vec3_t Triangle::local3DToGlobal3D(vec3_t l_M) {
 138   return m_a + m_G*l_M;
 139 }
 140
 141 vec3_t Triangle::global3DToLocal3D(vec3_t g_M) {
 142   vec3_t tmp = g_M - m_a;
 143   return m_GI*tmp;
 144 }
 145
 146 vec3_t Triangle::local2DToGlobal3D(vec2_t l_M) {
 147   return local3DToGlobal3D(vec3_t(l_M[0], l_M[1], 0));
 148 }
 149
 150 vec2_t Triangle::global3DToLocal2D(vec3_t g_M) {
 151   vec3_t l_M = global3DToLocal3D(g_M);
 152   return vec2_t(l_M[0], l_M[1]);
 153 }
 154
 155 bool Triangle::projectOnTriangle(vec3_t xp, vec3_t &xi, vec3_t &ri, double &d, int& side, bool restrict_to_triangle) {
 156   side = -1;
 157   double scal = (xp - this->m_a) * this->m_g3;
 158   vec3_t x1, x2;
 159   if (scal > 0) {
 160     x1 = xp + this->m_g3;
 161     x2 = xp - scal * this->m_g3 - this->m_g3;
 162   } else {
 163     x1 = xp - this->m_g3;
 164     x2 = xp - scal * this->m_g3 + this->m_g3;
 165   }
 166   // (xi,ri) gets set to the intersection of the line with the plane here!
 167   bool intersects_face = GeometryTools::intersectEdgeAndTriangle(this->m_a, this->m_b, this->m_c, x1, x2, xi, ri);
 168   if (intersects_face || !restrict_to_triangle) {
 169     vec3_t dx = xp - this->m_a;
 170     d = fabs(dx * this->m_g3);
 171   } else {
 172     double kab = GeometryTools::intersection(this->m_a, this->m_b - this->m_a, xp, this->m_b - this->m_a);
 173     double kac = GeometryTools::intersection(this->m_a, this->m_c - this->m_a, xp, this->m_c - this->m_a);
 174     double kbc = GeometryTools::intersection(this->m_b, this->m_c - this->m_b, xp, this->m_c - this->m_b);
 175
 176     double dab = (this->m_a + kab * (this->m_b - this->m_a) - xp).abs();
 177     double dac = (this->m_a + kac * (this->m_c - this->m_a) - xp).abs();
 178     double dbc = (this->m_b + kbc * (this->m_c - this->m_b) - xp).abs();
 179     double da = (this->m_a - xp).abs();
 180     double db = (this->m_b - xp).abs();
 181     double dc = (this->m_c - xp).abs();
 182
 183     bool set = false;
 184     d = 1e99;//max(max(max(max(max(dab,dac),dbc),da),db),dc);
 185
 186     if (dab < d) {
 187       if ((kab >= 0) && (kab <= 1)) {
 188         xi = this->m_a + kab * (this->m_b - this->m_a);
 189         ri = vec3_t(kab, 0, 0);
 190         d = dab;
 191         set = true;
 192         side = 0;
 193       }
 194     }
 195     if (dbc < d) {
 196       if ((kbc >= 0) && (kbc <= 1)) {
 197         xi = this->m_b + kbc * (this->m_c - this->m_b);
 198         ri = vec3_t(1 - kbc, kbc, 0);
 199         d = dbc;
 200         set = true;
 201         side = 1;
 202       }
 203     }
 204     if (dac < d) {
 205       if ((kac >= 0) && (kac <= 1)) {
 206         xi = this->m_a + kac * (this->m_c - this->m_a);
 207         ri = vec3_t(0, kac, 0);
 208         d = dac;
 209         set = true;
 210         side = 2;
 211       }
 212     }
 213     if (da < d) {
 214       xi = this->m_a;
 215       ri = vec3_t(0, 0);
 216       d = da;
 217       set = true;
 218       side = 3;
 219     }
 220     if (db < d) {
 221       xi = this->m_b;
 222       ri = vec3_t(1, 0);
 223       d = db;
 224       set = true;
 225       side = 4;
 226     }
 227     if (dc < d) {
 228       xi = this->m_c;
 229       ri = vec3_t(0, 1);
 230       d = dc;
 231       set = true;
 232       side = 5;
 233     }
 234     if (!set) {
 235       EG_BUG;
 236     }
 237   }
 238   if (xi[0] > 1e98) { // should never happen
 239     EG_BUG;
 240   }
 241   /*  if (not( 0<=ri[0] && ri[0]<=1 && 0<=ri[1] && ri[1]<=1 && ri[2]==0 )) {
 242       qWarning()<<"ri="<<ri;
 243       EG_BUG;
 244     }*/
 245   return intersects_face;
 246 }
 247
 248 void Triangle::saveTriangle(QString filename)
 249 {
 250   int N_cells = 1;
 251   int N_points = 3;
 252
 253   EG_VTKSP(vtkUnstructuredGrid, triangle_grid);
 254   allocateGrid(triangle_grid , N_cells, N_points);
 255
 256   vtkIdType node_count = 0;
 257   int cell_count = 0;
 258
 259   vtkIdType pts[3];
 260   triangle_grid->GetPoints()->SetPoint(node_count, m_a.data()); pts[0]=node_count; node_count++;
 261   triangle_grid->GetPoints()->SetPoint(node_count, m_b.data()); pts[1]=node_count; node_count++;
 262   triangle_grid->GetPoints()->SetPoint(node_count, m_c.data()); pts[2]=node_count; node_count++;
 263
 264   triangle_grid->InsertNextCell(VTK_TRIANGLE,3,pts);cell_count++;
 265
 266   saveGrid(triangle_grid, filename+"_triangle_grid");
 267 }