HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
SIM_DerVector3.h
Go to the documentation of this file.
1 /*
2  * PROPRIETARY INFORMATION. This software is proprietary to
3  * Side Effects Software Inc., and is not to be reproduced,
4  * transmitted, or disclosed in any way without written permission.
5  *
6  */
7 
8 #ifndef __SIM_DerVector3_h__
9 #define __SIM_DerVector3_h__
10 
11 #include "SIM_API.h"
12 #include <UT/UT_Matrix3.h>
13 #include <UT/UT_Vector3.h>
14 
15 class SIM_DerScalar;
16 
17 /// This class defines a 3D vector and its partial derivative w.r.t. another
18 /// 3D vector. It uses automatic differentiation to maintain the dependency
19 /// upon the derivative vector as arithmetic operations are performed.
20 /// The derivative of a vector-valued function is, of course, a Jacobian
21 /// matrix.
22 ///
23 /// By performing a sequence of arithmetic operations on this
24 /// vector class after initializing its derivative appropriately, you can
25 /// easily keep track of the effect of those operations on the derivative.
26 /// Independent variables can be included in an equation using the
27 /// conventional UT_Vector3 and fpreal types, and dependent variables can
28 /// use the SIM_DerVector3 and SIM_DerScalar types.
29 ///
30 /// It is inspired by Eitan Grinspun's class for the same purpose,
31 /// described in his 2003 SCA paper on Discrete Shells.
33 {
34 public:
36  /// Initialize to a constant vector, with no derivative.
37  explicit SIM_DerVector3(const UT_Vector3 &v) : myV(v), myD(0.f)
38  { }
39  /// Initialize to a vector with a derivative. This is particularly
40  /// useful for initializing the variables themselves, where D=I.
42  const UT_Matrix3 &D) : myV(v), myD(D)
43  { }
44 
45  // Default copy constructor is fine.
46  //SIM_DerVector3(const SIM_DerVector3 &rhs);
47 
48  /// The vector v.
49  const UT_Vector3&v() const
50  { return myV; }
51 
52  /// Derivative matrix, dv/dx.
53  /// The entries of the matrix are laid out like a typical Jacobian:
54  ///
55  /// [ dv1/dx1 dv1/dx2 dv1/dx3 ]
56  /// [ dv2/dx1 dv2/dx2 dv2/dx3 ]
57  /// [ dv3/dx1 dv3/dx2 dv3/dx3 ]
58  ///
59  /// [ dv1/dx ]
60  /// = [ dv2/dx ]
61  /// [ dv3/dx ]
62  ///
63  /// = [ dv/dx1 dv/dx2 dv/dx3 ]
64  const UT_Matrix3&D() const
65  { return myD; }
66 
67  // Default assignment operator is fine.
68  // SIM_DerVector3 operator=(const SIM_DerVector3 &rhs);
69 
71  {
72  return SIM_DerVector3(-v(), -D());
73  }
75  {
76  // d(v1+v2)/dx = dv1/dx + dv2/dx
77  return SIM_DerVector3(v() + rhs.v(), D() + rhs.D());
78  }
80  {
81  return SIM_DerVector3(v() + rhs, D());
82  }
84  {
85  // d(v1-v2)/dx = dv1/dx - dv2/dx
86  return SIM_DerVector3(v() - rhs.v(), D() - rhs.D());
87  }
89  {
90  return SIM_DerVector3(v() - rhs, D());
91  }
92  SIM_DerVector3 operator*(const SIM_DerScalar &rhs) const;
94  {
95  // d(v*s)/dx = s*dv/dx + v * (ds/dx)^T
96  return SIM_DerVector3(v() * scalar, D() * scalar);
97  }
99  { return operator=((*this) + rhs); }
101  { return operator=((*this) + rhs); }
103  { return operator=((*this) - rhs); }
105  { return operator=((*this) - rhs); }
107  { return operator=((*this) * rhs); }
109  { return operator=((*this) * rhs); }
110 
111  SIM_DerScalar dot(const SIM_DerVector3 &rhs) const;
112  SIM_DerScalar dot(const UT_Vector3 &rhs) const;
113  SIM_DerVector3 cross(const SIM_DerVector3 &rhs) const;
114  SIM_DerVector3 cross(const UT_Vector3 &rhs) const;
115  SIM_DerScalar length2() const;
116  SIM_DerScalar length() const;
117  SIM_DerVector3 normalize() const;
118 
119 
120  // Matrix corresponding to a vector cross-product.
121  // a x b = S(a) * b
122  // The matrix is skew-symmetric.
123  static UT_Matrix3 S(const UT_Vector3 &v)
124  {
125  return UT_Matrix3( 0, -v.z(), v.y(),
126  v.z(), 0, -v.x(),
127  -v.y(), v.x(), 0);
128  }
129 
130 private:
131  UT_Vector3 myV;
132  UT_Matrix3 myD;
133 };
134 
135 #include "SIM_DerScalar.h"
136 
137 inline
138 SIM_DerVector3 operator+(const UT_Vector3 &lhs, const SIM_DerVector3 &rhs);
139 inline
140 SIM_DerVector3 operator-(const UT_Vector3 &lhs, const SIM_DerVector3 &rhs);
141 inline
143 inline
145 inline
147 inline
149 inline
150 SIM_DerScalar dot(const SIM_DerVector3 &lhs, const SIM_DerVector3 &rhs);
151 inline
152 SIM_DerScalar dot(const SIM_DerVector3 &lhs, const UT_Vector3 &rhs);
153 inline
154 SIM_DerScalar dot(const UT_Vector3 &lhs, const SIM_DerVector3 &rhs);
155 inline
156 SIM_DerVector3 cross(const SIM_DerVector3 &lhs, const SIM_DerVector3 &rhs);
157 inline
158 SIM_DerVector3 cross(const SIM_DerVector3 &lhs, const UT_Vector3 &rhs);
159 inline
160 SIM_DerVector3 cross(const UT_Vector3 &lhs, const SIM_DerVector3 &rhs);
161 
162 
163 
164 inline SIM_DerVector3
166 {
167  // d(v*s)/dx = s*dv/dx + v * (ds/dx)^T
168  UT_Matrix3 newD(D());
169  newD *= rhs.v();
170  newD.outerproductUpdate(1, v(), rhs.D());
171 
172  return SIM_DerVector3(v() * rhs.v(), newD);
173 }
174 
175 inline SIM_DerScalar
177 {
178  // d(v1.v2)/dx = v1^T * dv2/dx + v2^T * dv1/dx
179  // Note the rowvector*matrix multiplication.
180  return SIM_DerScalar(::dot(v(), rhs.v()),
181  ::rowVecMult(rhs.v(), D()) +
182  ::rowVecMult(v(), rhs.D()));
183 }
184 
185 inline SIM_DerScalar
187 {
188  return SIM_DerScalar(::dot(v(), rhs),
189  ::rowVecMult(rhs, D()));
190 }
191 
192 // d(v1 x v2)/dx = dv1/dx x v2 + v1 x dv2/dx
193 // = -v2 x dv1/dx + v1 x dv2/dx
194 // = S(-v2) dv1/dx + S(v1) dv2/dx
195 inline SIM_DerVector3
197 {
198  return SIM_DerVector3(::cross(v(), rhs.v()),
199  S(-rhs.v()) * D() + S(v()) * rhs.D());
200 }
201 
202 // d(v1 x v2)/dx = dv1/dx x v2 + v1 x dv2/dx
203 // = -v2 x dv1/dx + v1 x dv2/dx
204 // = S(-v2) dv1/dx + S(v1) dv2/dx
205 inline SIM_DerVector3
207 {
208  return SIM_DerVector3(::cross(v(), rhs), S(-rhs) * D());
209 }
210 
211 // d|v|^2/dx = d|v.v|/dx
212 // = 2 * v * dv/dx
213 inline SIM_DerScalar
215 {
216  return SIM_DerScalar(v().length2(),
217  2 * ::rowVecMult(v(), D()));
218 }
219 
220 // d|v|/dx = d((v.v)^.5)/dx
221 // = .5 / |v| * d(v.v)/dx
222 // = v / |v| * dv/dx
223 
224 // Because it includes a square root, there is a discontinuity at the
225 // origin. Like the square root, I approximate using a zero derivative at
226 // the origin.
227 inline SIM_DerScalar
229 {
230  const fpreal tol = 1e-5;
231  const fpreal len = v().length();
232  if( len < tol )
233  return SIM_DerScalar(len);
234  else
235  return SIM_DerScalar(len, ::rowVecMult(v() / len, D()));
236 }
237 
238 inline SIM_DerVector3
240 {
241  // TODO: we can make this more efficient... can't we?
242  return (*this)/length();
243 }
244 
245 
246 
247 
248 
249 inline SIM_DerVector3
250 operator+(const UT_Vector3 &lhs, const SIM_DerVector3 &rhs)
251 {
252  return rhs + lhs;
253 }
254 
255 inline SIM_DerVector3
256 operator-(const UT_Vector3 &lhs, const SIM_DerVector3 &rhs)
257 {
258  return SIM_DerVector3(lhs - rhs.v(), -rhs.D());
259 }
260 
261 inline SIM_DerVector3
263 {
264  return v * s;
265 }
266 
267 inline SIM_DerVector3
269 {
270  return v * s;
271 }
272 
273 inline SIM_DerVector3
275 {
276  return v * s.inverse();
277 }
278 
279 inline SIM_DerVector3
281 {
282  return v * (1/s);
283 }
284 
285 inline SIM_DerScalar
286 dot(const SIM_DerVector3 &lhs, const SIM_DerVector3 &rhs)
287 {
288  return lhs.dot(rhs);
289 }
290 
291 inline SIM_DerScalar
292 dot(const SIM_DerVector3 &lhs, const UT_Vector3 &rhs)
293 {
294  return lhs.dot(rhs);
295 }
296 
297 inline SIM_DerScalar
298 dot(const UT_Vector3 &lhs, const SIM_DerVector3 &rhs)
299 {
300  return rhs.dot(lhs);
301 }
302 
303 inline SIM_DerVector3
304 cross(const SIM_DerVector3 &lhs, const SIM_DerVector3 &rhs)
305 {
306  return lhs.cross(rhs);
307 }
308 
309 inline SIM_DerVector3
310 cross(const SIM_DerVector3 &lhs, const UT_Vector3 &rhs)
311 {
312  return lhs.cross(rhs);
313 }
314 
315 inline SIM_DerVector3
316 cross(const UT_Vector3 &lhs, const SIM_DerVector3 &rhs)
317 {
318  // d(v1 x v2)/dx = dv1/dx x v2 + v1 x dv2/dx
319  // = -v2 x dv1/dx + v1 x dv2/dx
320  // = S(-v2) dv1/dx + S(v1) dv2/dx
321  return SIM_DerVector3(::cross(lhs, rhs.v()),
322  SIM_DerVector3::S(lhs) * rhs.D());
323 }
324 #endif
UT_Vector3T< T > rowVecMult(const UT_Vector3T< T > &v, const UT_Matrix3T< S > &m)
Definition: UT_Matrix3.h:1295
Mat3< typename promote< S, T >::type > operator*(S scalar, const Mat3< T > &m)
Multiply each element of the given matrix by scalar and return the result.
Definition: Mat3.h:609
SIM_DerVector3 & operator-=(const UT_Vector3 &rhs)
SIM_DerVector3 operator+(const SIM_DerVector3 &rhs) const
SIM_DerScalar dot(const SIM_DerVector3 &lhs, const SIM_DerVector3 &rhs)
SIM_DerVector3 & operator+=(const UT_Vector3 &rhs)
SIM_DerVector3 & operator*=(const SIM_DerScalar &rhs)
const GLdouble * v
Definition: glcorearb.h:836
Mat3< typename promote< T0, T1 >::type > operator+(const Mat3< T0 > &m0, const Mat3< T1 > &m1)
Add corresponding elements of m0 and m1 and return the result.
Definition: Mat3.h:625
static UT_Matrix3 S(const UT_Vector3 &v)
SIM_DerScalar dot(const SIM_DerVector3 &rhs) const
SYS_FORCE_INLINE T & x(void)
Definition: UT_Vector3.h:498
SIM_DerScalar inverse() const
SIM_DerVector3(const UT_Vector3 &v)
Initialize to a constant vector, with no derivative.
SIM_DerVector3 normalize() const
SYS_FORCE_INLINE T & z(void)
Definition: UT_Vector3.h:502
SIM_DerVector3 & operator+=(const SIM_DerVector3 &rhs)
SYS_FORCE_INLINE void outerproductUpdate(T b, const UT_Vector3F &v1, const UT_Vector3F &v2)
Definition: UT_Matrix3.h:429
GLfloat f
Definition: glcorearb.h:1925
SIM_DerVector3 operator-() const
UT_Matrix3T< float > UT_Matrix3
SIM_DerScalar length2() const
fpreal v() const
Definition: SIM_DerScalar.h:46
SIM_DerScalar length() const
Mat3< typename promote< T0, T1 >::type > operator-(const Mat3< T0 > &m0, const Mat3< T1 > &m1)
Subtract corresponding elements of m0 and m1 and return the result.
Definition: Mat3.h:635
SIM_DerVector3 operator*(fpreal scalar) const
SIM_DerVector3(const UT_Vector3 &v, const UT_Matrix3 &D)
SIM_DerVector3 & operator*=(const fpreal rhs)
SIM_DerVector3 cross(const SIM_DerVector3 &rhs) const
SIM_DerVector3 & operator-=(const SIM_DerVector3 &rhs)
const UT_Vector3 & D() const
Definition: SIM_DerScalar.h:49
GridType::Ptr normalize(const GridType &grid, bool threaded, InterruptT *interrupt)
Normalize the vectors of the given vector-valued grid.
SIM_DerVector3 operator*(const SIM_DerScalar &s, const SIM_DerVector3 &v)
SYS_FORCE_INLINE T & y(void)
Definition: UT_Vector3.h:500
double fpreal
Definition: SYS_Types.h:270
const UT_Matrix3 & D() const
const UT_Vector3 & v() const
The vector v.
SIM_DerVector3 operator*(const SIM_DerScalar &rhs) const
#define SIM_API
Definition: SIM_API.h:10
SIM_DerVector3 operator+(const UT_Vector3 &rhs) const
SYS_FORCE_INLINE Storage::MathFloat length() const
SIM_DerVector3 operator-(const SIM_DerVector3 &rhs) const
SIM_DerVector3 operator/(const SIM_DerVector3 &v, const SIM_DerScalar &s)
SIM_DerVector3 cross(const SIM_DerVector3 &lhs, const SIM_DerVector3 &rhs)
SIM_DerVector3 operator-(const UT_Vector3 &rhs) const
GLuint GLsizei GLsizei * length
Definition: glcorearb.h:794