VectorAnalysis.h 6.98 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42
//=============================================================================
//
//  Copyright (c) Kitware, Inc.
//  All rights reserved.
//  See LICENSE.txt for details.
//
//  This software is distributed WITHOUT ANY WARRANTY; without even
//  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
//  PURPOSE.  See the above copyright notice for more information.
//
//  Copyright 2015 Sandia Corporation.
//  Copyright 2015 UT-Battelle, LLC.
//  Copyright 2015 Los Alamos National Security.
//
//  Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
//  the U.S. Government retains certain rights in this software.
//  Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
//  Laboratory (LANL), the U.S. Government retains certain rights in
//  this software.
//
//=============================================================================
#ifndef vtk_m_VectorAnalysis_h
#define vtk_m_VectorAnalysis_h

// This header file defines math functions that deal with linear albegra funcitons

#include <vtkm/Math.h>
#include <vtkm/Types.h>
#include <vtkm/TypeTraits.h>
#include <vtkm/VecTraits.h>

namespace vtkm {


// ----------------------------------------------------------------------------
/// \brief Returns the linear interpolation of two values based on weight
///
/// \c Lerp interpolates return the linerar interpolation of v0 and v1 based on w. v0
/// and v1 are scalars or vectors of same length. w can either be a scalar or a
/// vector of the same length as x and y. If w is outside [0,1] then lerp
/// extrapolates. If w=0 => v0 is returned if w=1 => v1 is returned.
///
43
template<typename ValueType, typename WeightType>
44
VTKM_EXEC_CONT
45
ValueType Lerp(const ValueType &value0,
46
               const ValueType &value1,
47
               const WeightType &weight)
48
{
49
  return static_cast<ValueType>((WeightType(1)-weight)*value0+weight*value1);
50
}
51
template<typename ValueType, vtkm::IdComponent N, typename WeightType>
52
VTKM_EXEC_CONT
53 54 55
vtkm::Vec<ValueType,N> Lerp(const vtkm::Vec<ValueType,N> &value0,
                            const vtkm::Vec<ValueType,N> &value1,
                            const WeightType &weight)
56
{
57
  return (WeightType(1)-weight)*value0+weight*value1;
58 59
}
template<typename ValueType, vtkm::IdComponent N>
60
VTKM_EXEC_CONT
61 62 63 64
vtkm::Vec<ValueType,N> Lerp(const vtkm::Vec<ValueType,N> &value0,
                            const vtkm::Vec<ValueType,N> &value1,
                            const vtkm::Vec<ValueType,N> &weight)
{
65 66
  static const vtkm::Vec<ValueType,N> One(ValueType(1));
  return (One-weight)*value0+weight*value1;
67 68 69 70 71 72 73 74 75 76
}

// ----------------------------------------------------------------------------
/// \brief Returns the square of the magnitude of a vector.
///
/// It is usually much faster to compute the square of the magnitude than the
/// square, so you should use this function in place of Magnitude or RMagnitude
/// when possible.
///
template<typename T>
77
VTKM_EXEC_CONT
78 79 80 81 82 83 84 85 86
typename vtkm::VecTraits<T>::ComponentType
MagnitudeSquared(const T &x)
{
  return vtkm::dot(x,x);
}

// ----------------------------------------------------------------------------
namespace detail {
template<typename T>
87
VTKM_EXEC_CONT
88 89 90 91 92 93
T MagnitudeTemplate(T x, vtkm::TypeTraitsScalarTag)
{
  return vtkm::Abs(x);
}

template<typename T>
94
VTKM_EXEC_CONT
95 96 97 98 99 100 101 102
vtkm::FloatDefault MagnitudeTemplate(const T &x, vtkm::TypeTraitsVectorTag)
{
  return vtkm::Sqrt(vtkm::MagnitudeSquared(x));
}

template<vtkm::IdComponent Size>
VTKM_EXEC_CONT
vtkm::Float64 MagnitudeTemplate(const vtkm::Vec<vtkm::Float64, Size> &x)
103 104 105
{
  return vtkm::Sqrt(vtkm::MagnitudeSquared(x));
}
106

107 108 109 110 111 112 113 114 115 116 117
} // namespace detail

/// \brief Returns the magnitude of a vector.
///
/// It is usually much faster to compute MagnitudeSquared, so that should be
/// substituted when possible (unless you are just going to take the square
/// root, which would be besides the point). On some hardware it is also faster
/// to find the reciprocal magnitude, so RMagnitude should be used if you
/// actually plan to divide by the magnitude.
///
template<typename T>
118
VTKM_EXEC_CONT
119
vtkm::FloatDefault Magnitude(const T &x)
120 121 122 123
{
  return detail::MagnitudeTemplate(
        x, typename vtkm::TypeTraits<T>::DimensionalityTag());
}
124 125 126 127 128 129
template<vtkm::IdComponent Size>
VTKM_EXEC_CONT
vtkm::Float64 Magnitude(const vtkm::Vec<vtkm::Float64, Size> &x)
{
  return detail::MagnitudeTemplate(x);
}
130 131 132 133

// ----------------------------------------------------------------------------
namespace detail {
template<typename T>
134
VTKM_EXEC_CONT
135 136 137 138 139 140
T RMagnitudeTemplate(T x, vtkm::TypeTraitsScalarTag)
{
  return T(1)/vtkm::Abs(x);
}

template<typename T>
141
VTKM_EXEC_CONT
142 143 144 145 146 147 148 149 150 151 152 153 154
typename vtkm::VecTraits<T>::ComponentType
RMagnitudeTemplate(const T &x, vtkm::TypeTraitsVectorTag)
{
  return vtkm::RSqrt(vtkm::MagnitudeSquared(x));
}
} // namespace detail

/// \brief Returns the reciprocal magnitude of a vector.
///
/// On some hardware RMagnitude is faster than Magnitude, but neither is
/// as fast as MagnitudeSquared.
///
template<typename T>
155
VTKM_EXEC_CONT
156 157 158 159 160 161 162 163 164 165
typename vtkm::VecTraits<T>::ComponentType
RMagnitude(const T &x)
{
  return detail::RMagnitudeTemplate(
        x, typename vtkm::TypeTraits<T>::DimensionalityTag());
}

// ----------------------------------------------------------------------------
namespace detail {
template<typename T>
166
VTKM_EXEC_CONT
167 168 169 170 171 172
T NormalTemplate(T x, vtkm::TypeTraitsScalarTag)
{
  return vtkm::CopySign(T(1), x);
}

template<typename T>
173
VTKM_EXEC_CONT
174 175 176 177 178 179 180 181 182 183 184
T NormalTemplate(const T &x, vtkm::TypeTraitsVectorTag)
{
  return vtkm::RMagnitude(x)*x;
}
} // namespace detail

/// \brief Returns a normalized version of the given vector.
///
/// The resulting vector points in the same direction but has unit length.
///
template<typename T>
185
VTKM_EXEC_CONT
186 187 188 189 190 191 192 193 194 195 196 197
T Normal(const T &x)
{
  return detail::NormalTemplate(
        x, typename vtkm::TypeTraits<T>::DimensionalityTag());
}

// ----------------------------------------------------------------------------
/// \brief Changes a vector to be normal.
///
/// The given vector is scaled to be unit length.
///
template<typename T>
198
VTKM_EXEC_CONT
199 200 201 202 203 204 205 206 207
void Normalize(T &x)
{
  x = vtkm::Normal(x);
}

// ----------------------------------------------------------------------------
/// \brief Find the cross product of two vectors.
///
template<typename T>
208
VTKM_EXEC_CONT
209 210 211 212 213 214 215 216 217 218 219 220 221 222 223
vtkm::Vec<T,3> Cross(const vtkm::Vec<T,3> &x, const vtkm::Vec<T,3> &y)
{
  return vtkm::Vec<T,3>(x[1]*y[2] - x[2]*y[1],
                        x[2]*y[0] - x[0]*y[2],
                        x[0]*y[1] - x[1]*y[0]);
}

//-----------------------------------------------------------------------------
/// \brief Find the normal of a triangle.
///
/// Given three coordinates in space, which, unless degenerate, uniquely define
/// a triangle and the plane the triangle is on, returns a vector perpendicular
/// to that triangle/plane.
///
template<typename T>
224
VTKM_EXEC_CONT
225 226 227 228 229 230 231 232 233 234 235
vtkm::Vec<T,3> TriangleNormal(const vtkm::Vec<T,3> &a,
                              const vtkm::Vec<T,3> &b,
                              const vtkm::Vec<T,3> &c)
{
  return vtkm::Cross(b-a, c-a);
}


} // namespace vtkm

#endif //vtk_m_VectorAnalysis_h