Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
v_vector.h
Go to the documentation of this file.
1 /*=========================================================================
2  Original version -- Copyright (c) 2006 Sandia Corporation.
3  All rights reserved.
4  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
5 
6  This software is distributed WITHOUT ANY WARRANTY; without even
7  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
8  PURPOSE. See the above copyright notice for more information.
9 
10 =========================================================================*/
11 
12 #ifndef VERDICT_VECTOR
13 #define VERDICT_VECTOR
14 
15 #include "moab/verdict.h"
16 #include <cmath>
17 #include <cassert>
18 
19 // computes the dot product of 3d vectors
20 inline double dot_product( double vec1[], double vec2[] )
21 {
22  double answer = vec1[0] * vec2[0] + vec1[1] * vec2[1] + vec1[2] * vec2[2];
23  return answer;
24 }
25 
26 // normalize a vector
27 inline void normalize( double vec[] )
28 {
29  double x = sqrt( vec[0] * vec[0] + vec[1] * vec[1] + vec[2] * vec[2] );
30 
31  vec[0] /= x;
32  vec[1] /= x;
33  vec[2] /= x;
34 }
35 
36 // computes the cross product
37 inline double* cross_product( double vec1[], double vec2[], double answer[] )
38 {
39  answer[0] = vec1[1] * vec2[2] - vec1[2] * vec2[1];
40  answer[1] = vec1[2] * vec2[0] - vec1[0] * vec2[2];
41  answer[2] = vec1[0] * vec2[1] - vec1[1] * vec2[0];
42  return answer;
43 }
44 
45 // computes the length of a vector
46 inline double length( double vec[] )
47 {
48  return sqrt( vec[0] * vec[0] + vec[1] * vec[1] + vec[2] * vec[2] );
49 }
50 
51 // computes the square length of a vector
52 inline double length_squared( double vec[] )
53 {
54  return ( vec[0] * vec[0] + vec[1] * vec[1] + vec[2] * vec[2] );
55 }
56 
57 // computes the interior angle between 2 vectors in degrees
58 inline double interior_angle( double vec1[], double vec2[] )
59 {
60  double cosAngle, angleRad;
61  double length1 = length( vec1 );
62  double length2 = length( vec2 );
63  assert( ( length1 > 0.0 && length2 > 0.0 ) );
64 
65  cosAngle = dot_product( vec1, vec2 ) / ( length1 * length2 );
66 
67  if( ( cosAngle > 1.0 ) && ( cosAngle < 1.0001 ) )
68  {
69  cosAngle = 1.0;
70  }
71  else if( cosAngle < -1.0 && cosAngle > -1.0001 )
72  {
73  cosAngle = -1.0;
74  }
75  else
76  {
77  assert( cosAngle < 1.0001 && cosAngle > -1.0001 );
78  }
79 
80  angleRad = acos( cosAngle );
81  return ( ( angleRad * 180. ) / VERDICT_PI );
82 }
83 
84 #endif