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
  28 Triangle::Triangle() {
  29   id_a = 0;
  30   id_b = 0;
  31   id_c = 0;
  32   a = vec3_t(0, 0, 0);
  33   b = vec3_t(0, 0, 0);
  34   c = vec3_t(0, 0, 0);
  35   setupTriangle();
  36 }
  37
  38 Triangle::Triangle(vec3_t a_a, vec3_t a_b, vec3_t a_c) {
  39   a = a_a;
  40   b = a_b;
  41   c = a_c;
  42   setupTriangle();
  43 }
  44
  45 Triangle::Triangle(vtkUnstructuredGrid* a_grid, vtkIdType a_id_a, vtkIdType a_id_b, vtkIdType a_id_c) {
  46   id_a = a_id_a;
  47   id_b = a_id_b;
  48   id_c = a_id_c;
  49   a_grid->GetPoints()->GetPoint(id_a, a.data());
  50   a_grid->GetPoints()->GetPoint(id_b, b.data());
  51   a_grid->GetPoints()->GetPoint(id_c, c.data());
  52   setupTriangle();
  53 }
  54
  55 Triangle::Triangle(vtkUnstructuredGrid* a_grid, vtkIdType a_id_cell) {
  56   vtkIdType Npts, *pts;
  57   a_grid->GetCellPoints(a_id_cell, Npts, pts);
  58   if (Npts == 3) {
  59     id_a = pts[0];
  60     id_b = pts[1];
  61     id_c = pts[2];
  62     a_grid->GetPoints()->GetPoint(id_a, a.data());
  63     a_grid->GetPoints()->GetPoint(id_b, b.data());
  64     a_grid->GetPoints()->GetPoint(id_c, c.data());
  65     setupTriangle();
  66   } else {
  67     EG_ERR_RETURN("only triangles allowed at the moment");
  68   }
  69 }
  70
  71 void Triangle::setupTriangle() {
  72   m_has_neighbour.resize(6);
  73   m_has_neighbour[0] = false;
  74   m_has_neighbour[1] = false;
  75   m_has_neighbour[2] = false;
  76   m_has_neighbour[3] = false;
  77   m_has_neighbour[4] = false;
  78   m_has_neighbour[5] = false;
  79
  80   g1 = b - a;
  81   g2 = c - a;
  82   g3 = g1.cross(g2);
  83
  84   A  = 0.5 * g3.abs();
  85   g3.normalise();
  86
  87   G.column(0, g1);
  88   G.column(1, g2);
  89   G.column(2, g3);
  90   GI = G.inverse();
  91
  92   smallest_length = (b - a).abs();
  93   smallest_length = min(smallest_length, (c - b).abs());
  94   smallest_length = min(smallest_length, (a - c).abs());
  95 }
  96
  97 vec3_t Triangle::local3DToGlobal3D(vec3_t l_M) {
  98   return a + G*l_M;
  99 }
 100
 101 vec3_t Triangle::global3DToLocal3D(vec3_t g_M) {
 102   vec3_t tmp = g_M - a;
 103   return GI*tmp;
 104 }
 105
 106 vec3_t Triangle::local2DToGlobal3D(vec2_t l_M) {
 107   return local3DToGlobal3D(vec3_t(l_M[0], l_M[1], 0));
 108 }
 109
 110 vec2_t Triangle::global3DToLocal2D(vec3_t g_M) {
 111   vec3_t l_M = global3DToLocal3D(g_M);
 112   return vec2_t(l_M[0], l_M[1]);
 113 }
 114
 115 bool Triangle::projectOnTriangle(vec3_t xp, vec3_t &xi, vec3_t &ri, double &d, int& side, bool restrict_to_triangle) {
 116   side = -1;
 117   double scal = (xp - this->a) * this->g3;
 118   vec3_t x1, x2;
 119   if (scal > 0) {
 120     x1 = xp + this->g3;
 121     x2 = xp - scal * this->g3 - this->g3;
 122   } else {
 123     x1 = xp - this->g3;
 124     x2 = xp - scal * this->g3 + this->g3;
 125   }
 126   // (xi,ri) gets set to the intersection of the line with the plane here!
 127   bool intersects_face = GeometryTools::intersectEdgeAndTriangle(this->a, this->b, this->c, x1, x2, xi, ri);
 128   if (intersects_face || !restrict_to_triangle) {
 129     vec3_t dx = xp - this->a;
 130     d = fabs(dx * this->g3);
 131   } else {
 132     double kab = GeometryTools::intersection(this->a, this->b - this->a, xp, this->b - this->a);
 133     double kac = GeometryTools::intersection(this->a, this->c - this->a, xp, this->c - this->a);
 134     double kbc = GeometryTools::intersection(this->b, this->c - this->b, xp, this->c - this->b);
 135
 136     double dab = (this->a + kab * (this->b - this->a) - xp).abs();
 137     double dac = (this->a + kac * (this->c - this->a) - xp).abs();
 138     double dbc = (this->b + kbc * (this->c - this->b) - xp).abs();
 139     double da = (this->a - xp).abs();
 140     double db = (this->b - xp).abs();
 141     double dc = (this->c - xp).abs();
 142
 143     bool set = false;
 144     d = 1e99;//max(max(max(max(max(dab,dac),dbc),da),db),dc);
 145
 146     if (dab < d) {
 147       if ((kab >= 0) && (kab <= 1)) {
 148         xi = this->a + kab * (this->b - this->a);
 149         ri = vec3_t(kab, 0, 0);
 150         d = dab;
 151         set = true;
 152         side = 0;
 153       }
 154     }
 155     if (dbc < d) {
 156       if ((kbc >= 0) && (kbc <= 1)) {
 157         xi = this->b + kbc * (this->c - this->b);
 158         ri = vec3_t(1 - kbc, kbc, 0);
 159         d = dbc;
 160         set = true;
 161         side = 1;
 162       }
 163     }
 164     if (dac < d) {
 165       if ((kac >= 0) && (kac <= 1)) {
 166         xi = this->a + kac * (this->c - this->a);
 167         ri = vec3_t(0, kac, 0);
 168         d = dac;
 169         set = true;
 170         side = 2;
 171       }
 172     }
 173     if (da < d) {
 174       xi = this->a;
 175       ri = vec3_t(0, 0);
 176       d = da;
 177       set = true;
 178       side = 3;
 179     }
 180     if (db < d) {
 181       xi = this->b;
 182       ri = vec3_t(1, 0);
 183       d = db;
 184       set = true;
 185       side = 4;
 186     }
 187     if (dc < d) {
 188       xi = this->c;
 189       ri = vec3_t(0, 1);
 190       d = dc;
 191       set = true;
 192       side = 5;
 193     }
 194     if (!set) {
 195       EG_BUG;
 196     }
 197   }
 198   if (xi[0] > 1e98) { // should never happen
 199     EG_BUG;
 200   }
 201   /*  if (not( 0<=ri[0] && ri[0]<=1 && 0<=ri[1] && ri[1]<=1 && ri[2]==0 )) {
 202       qWarning()<<"ri="<<ri;
 203       EG_BUG;
 204     }*/
 205   return intersects_face;
 206 }