src/geometrytools.h

   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 #ifndef geometrytools_H
  24 #define geometrytools_H
  25
  26 #include "math/mathvector.h"
  27 #include "math/smallsquarematrix.h"
  28
  29 #include <vtkUnstructuredGrid.h>
  30
  31 namespace GeometryTools
  32 {
  33
  34
  35 double rad2deg( double rad ); ///< Converts radians to degrees
  36 double deg2rad( double deg ); ///< Converts degrees to radians
  37
  38 void rotate(vec3_t g1, vec3_t g2, vec3_t g3, vec3_t &b, double theta);
  39
  40 /** Rotates vector v around vector axis by an angle theta */
  41 vec3_t rotate(vec3_t v, vec3_t axis, double theta);
  42
  43 /** Returns a normalized vector orthogonal to v */
  44 vec3_t orthogonalVector(vec3_t v);
  45
  46 /** Calculates the intersection between a line and a plane.
  47  * @param x_straight A point of the line.
  48  * @param v_straight Direction of the line.
  49  * @param x_plane A point of the plane
  50  * @param n_plane Normal of the plane
  51  */
  52 double intersection(vec3_t x_straight, vec3_t v_straight,
  53                     vec3_t x_plane, vec3_t n_plane);
  54
  55 /** Calculates the intersection between a line and a plane.
  56  * @param x_straight A point of the line.
  57  * @param v_straight Direction of the line.
  58  * @param x_plane A point of the plane
  59  * @param u_plane A vector of the plane
  60  * @param v_plane Another vector of the plane not colinear with u_plane
  61  */
  62 double intersection(vec3_t x_straight, vec3_t v_straight,
  63                     vec3_t x_plane, vec3_t u_plane, vec3_t v_plane);
  64
  65 /**
  66  * Calculates the intersection of a segment [x1,x2] with a triangle (a,b,c)
  67  * Note: (xi,ri) will always be set to the intersection of the line (x1,x2) with the plane (a,b,c) even if the segment does not intersect the triangle!
  68  * @param a Input: Triangle point 1
  69  * @param b Input: Triangle point 2
  70  * @param c Input: Triangle point 3
  71  * @param x1 Input: Segment point 1
  72  * @param x2 Input: Segment point 2
  73  * @param xi Output: 3D Global coordinates of the intersection point
  74  * @param ri Output: 3D local triangle coordinates of the intersection point
  75  * @param tol Input: Relative tolerance: There can only be an intersection if:
  76  * 0-tol*(x1-x2).abs()<=k<=1+tol*(x1-x2).abs()
  77  * 0-tol<=ri[0]<=1+tol
  78  * 0-tol<=ri[1]<=1+tol
  79  * @return true if an intersection point was found, else false.
  80 */
  81 bool intersectEdgeAndTriangle(const vec3_t& a, const vec3_t& b, const vec3_t& c,
  82                               const vec3_t& x1, const vec3_t& x2, vec3_t& xi, vec3_t& ri, double tol = 1e-4);
  83
  84 bool isInsideTriangle(vec2_t t_M, double tol = 1e-4);
  85
  86 /** Calculates the intersection point M between the lines (r1,u1) and (r2,u2).
  87  * @param k1 Returned by reference. Verifies M = r1+k1*u1
  88  * @param k2 Returned by reference. Verifies M = r2+k2*u2
  89  * @param r1 point of line 1
  90  * @param r2 point of line 2
  91  * @param u1 direction of line 1
  92  * @param u2 direction of line 2
  93  * @return true if an intersection has been found, else false.
  94  */
  95 bool intersection (double &k1, double &k2, vec2_t r1, vec2_t u1, vec2_t r2, vec2_t u2);
  96
  97 void sliceTriangle(const vector<vec3_t> &Tin, vec3_t x, vec3_t n, vector<vector<vec3_t> > &Tout);
  98
  99 /** Returns the volume of a tetrahedron.
 100  * V= v1*(v2^v3) with vi=xi-x0
 101  * If neg = false and V<0, it will return V=-1e99, else it returns V.
 102  * @param x0 point 0 of the tetrahedron
 103  * @param x1 point 1 of the tetrahedron
 104  * @param x2 point 2 of the tetrahedron
 105  * @param x3 point 3 of the tetrahedron
 106  * @param neg If neg = false and V<0, it will return V=-1e99, else it returns V.
 107  * @return volume of the tetrahedron
 108  */
 109 double tetraVol(const vec3_t& x0, const vec3_t& x1, const vec3_t& x2, const vec3_t& x3, bool neg = false);
 110
 111 double pyraVol(vec3_t x1, vec3_t x2, vec3_t x3, vec3_t x4, vec3_t x5, bool neg = false);
 112
 113 double prismVol(vec3_t x1, vec3_t x2, vec3_t x3, vec3_t x4, vec3_t x5, vec3_t x6, bool neg = false);
 114
 115 double hexaVol(vec3_t x1, vec3_t x2, vec3_t x3, vec3_t x4, vec3_t x5, vec3_t x6, vec3_t x7, vec3_t x8,  bool neg = false);
 116
 117 double triArea(vec3_t x1, vec3_t x2, vec3_t x3);
 118
 119 double quadArea(vec3_t x1, vec3_t x2, vec3_t x3, vec3_t x4);
 120
 121 vec3_t triNormal(vec3_t x0, vec3_t x1, vec3_t x2);///< Returns the normal of the surface defined by x0,x1,x2
 122
 123 vec3_t quadNormal(vec3_t x0, vec3_t x1, vec3_t x2, vec3_t x3);///< Returns the normal of the surface defined by x0,x1,x2,x3
 124
 125 vec3_t triNormal(vtkUnstructuredGrid *grid, vtkIdType p1, vtkIdType p2, vtkIdType p3);
 126
 127 vec3_t quadNormal(vtkUnstructuredGrid *grid, vtkIdType p1, vtkIdType p2, vtkIdType p3, vtkIdType p4);
 128
 129 vec3_t cellNormal(vtkUnstructuredGrid *grid, vtkIdType i);
 130
 131 /// Returns the area or volume of a cell.
 132 double cellVA(vtkUnstructuredGrid *grid, vtkIdType cellId, bool neg = false);
 133
 134 inline vec2_t turnRight(const vec2_t &v)
 135 {
 136   vec2_t u;
 137   u[0] =  v[1];
 138   u[1] = -v[0];
 139   return u;
 140 };
 141
 142 inline vec2_t turnLeft(const vec2_t &v)
 143 {
 144   vec2_t u;
 145   u[0] = -v[1];
 146   u[1] =  v[0];
 147   return u;
 148 };
 149
 150 //polygon must be numbered clockwise
 151 inline bool IsConvex(vec3_t a,vec3_t b,vec3_t c,vec3_t d)
 152 {
 153   vec3_t u[4];
 154   u[0]=b-a;
 155   u[1]=c-b;
 156   u[2]=d-c;
 157   u[3]=a-d;
 158
 159   for(int i=0;i<4;i++) {
 160     vec3_t n=u[i].cross(u[(i+1)%4]);
 161     if(n[2]>0) return(false);
 162   }
 163   return(true);
 164 };
 165
 166 inline bool IsConvex(vec2_t a_2D,vec2_t b_2D,vec2_t c_2D,vec2_t d_2D)
 167 {
 168   vec3_t a_3D(a_2D[0],a_2D[1]);
 169   vec3_t b_3D(b_2D[0],b_2D[1]);
 170   vec3_t c_3D(c_2D[0],c_2D[1]);
 171   vec3_t d_3D(d_2D[0],d_2D[1]);
 172   return(IsConvex(a_3D,b_3D,c_3D,d_3D));
 173 };
 174
 175 /// return the angle with relation to another 3-vector
 176 double angle(const vec3_t & u, const vec3_t & v);
 177
 178 /** return the deviation p1->p2->p3 (angle(p2-p1,p3-p2)) */
 179 double deviation(vtkUnstructuredGrid *grid, vtkIdType p1, vtkIdType p2, vtkIdType p3);
 180
 181 /** return the angle p1,p2,p3 (angle(p1-p2,p3-p2)) */
 182 double angle(vtkUnstructuredGrid *grid, vtkIdType p1, vtkIdType p2, vtkIdType p3);
 183
 184 /** return the cosine of the angle between the normals of cell1 and cell2 */
 185 double CosAngle(vtkUnstructuredGrid *grid, vtkIdType cell1, vtkIdType cell2);
 186
 187 /** Returns the center of mass of cellId + passes the minimal and maximal center to corner distances by reference */
 188 vec3_t getCenter(vtkUnstructuredGrid *grid, vtkIdType cellId, double& Rmin, double& Rmax);
 189
 190 /** Returns the distance between id_node1 and id_node2 */
 191 double distance(vtkUnstructuredGrid *grid, vtkIdType id_node1, vtkIdType id_node2);
 192
 193 /** Returns the distance squared between id_node1 and id_node2 */
 194 double distance2(vtkUnstructuredGrid *grid, vtkIdType id_node1, vtkIdType id_node2);
 195
 196 /** area of the circumscribed circle of the triangle */
 197 double areaOfCircumscribedCircle(vtkUnstructuredGrid *grid, vtkIdType id_cell);
 198
 199 vec3_t getBarycentricCoordinates(double x, double y);
 200
 201 vec3_t intersectionOnPlane(vec3_t v, vec3_t A, vec3_t nA, vec3_t B, vec3_t nB);
 202
 203 /** Projects vector V onto plane (O,i,j)
 204  * @param V The vector to project
 205  * @param i A vector of the plane
 206  * @param j A vector of the plane
 207  * @return Returns a 2D vector (x,y) so that V = x*i +y*j
 208  */
 209 vec2_t projectVectorOnPlane(vec3_t V,vec3_t i,vec3_t j);
 210
 211 vec3_t projectPointOnEdge(const vec3_t& M,const vec3_t& A, const vec3_t& u);
 212
 213 vec3_t projectPointOnPlane(const vec3_t& M, const vec3_t& A, const vec3_t& N);
 214
 215 };
 216
 217 #endif