Commit: 4e38d9df9ed851b4f034be69fb0dcd6d6576219d
Author: Lukas Tönne
Date: Sun Oct 12 12:45:07 2014 +0200
Branches: gooseberry
https://developer.blender.org/rB4e38d9df9ed851b4f034be69fb0dcd6d6576219d
Updating Eigen implicit dynamics solver implementation to adhere to the
new mass-spring solver API.
===================================================================
M source/blender/physics/intern/implicit.h
M source/blender/physics/intern/implicit_eigen.cpp
===================================================================
diff --git a/source/blender/physics/intern/implicit.h
b/source/blender/physics/intern/implicit.h
index ad714de..4659bb0 100644
--- a/source/blender/physics/intern/implicit.h
+++ b/source/blender/physics/intern/implicit.h
@@ -42,8 +42,8 @@
extern "C" {
#endif
-//#define IMPLICIT_SOLVER_EIGEN
-#define IMPLICIT_SOLVER_BLENDER
+#define IMPLICIT_SOLVER_EIGEN
+//#define IMPLICIT_SOLVER_BLENDER
#define CLOTH_ROOT_FRAME /* enable use of root frame coordinate transform */
diff --git a/source/blender/physics/intern/implicit_eigen.cpp
b/source/blender/physics/intern/implicit_eigen.cpp
index d391907..402ffcb 100644
--- a/source/blender/physics/intern/implicit_eigen.cpp
+++ b/source/blender/physics/intern/implicit_eigen.cpp
@@ -36,6 +36,12 @@
//#define USE_EIGEN_CORE
#define USE_EIGEN_CONSTRAINED_CG
+#ifdef __GNUC__
+# pragma GCC diagnostic push
+/* XXX suppress verbose warnings in eigen */
+# pragma GCC diagnostic ignored "-Wlogical-op"
+#endif
+
#ifndef IMPLICIT_ENABLE_EIGEN_DEBUG
#ifdef NDEBUG
#define IMPLICIT_NDEBUG
@@ -58,6 +64,10 @@
#endif
#endif
+#ifdef __GNUC__
+# pragma GCC diagnostic pop
+#endif
+
#include "MEM_guardedalloc.h"
extern "C" {
@@ -81,13 +91,111 @@ extern "C" {
typedef float Scalar;
-typedef Eigen::SparseMatrix<Scalar> lMatrix;
+static float I[3][3] = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
+static float ZERO[3][3] = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
+
+/* slightly extended Eigen vector class
+ * with conversion to/from plain C float array
+ */
+class fVector : public Eigen::Vector3f {
+public:
+ typedef float *ctype;
+
+ fVector()
+ {
+ }
+
+ fVector(const ctype &v)
+ {
+ for (int k = 0; k < 3; ++k)
+ coeffRef(k) = v[k];
+ }
+
+ fVector &operator = (const ctype &v)
+ {
+ for (int k = 0; k < 3; ++k)
+ coeffRef(k) = v[k];
+ return *this;
+ }
+
+ operator ctype()
+ {
+ return data();
+ }
+};
+
+/* slightly extended Eigen matrix class
+ * with conversion to/from plain C float array
+ */
+class fMatrix : public Eigen::Matrix3f {
+public:
+ typedef float (*ctype)[3];
+
+ fMatrix()
+ {
+ }
+
+ fMatrix(const ctype &v)
+ {
+ for (int k = 0; k < 3; ++k)
+ for (int l = 0; l < 3; ++l)
+ coeffRef(l, k) = v[k][l];
+ }
+
+ fMatrix &operator = (const ctype &v)
+ {
+ for (int k = 0; k < 3; ++k)
+ for (int l = 0; l < 3; ++l)
+ coeffRef(l, k) = v[k][l];
+ return *this;
+ }
+
+ operator ctype()
+ {
+ return (ctype)data();
+ }
+};
typedef Eigen::VectorXf lVector;
typedef Eigen::Triplet<Scalar> Triplet;
typedef std::vector<Triplet> TripletList;
+typedef Eigen::SparseMatrix<Scalar> lMatrix;
+
+struct lMatrixCtor {
+ lMatrixCtor(int numverts) :
+ m_numverts(numverts)
+ {
+ /* reserve for diagonal entries */
+ m_trips.reserve(numverts * 9);
+ }
+
+ int numverts() const { return m_numverts; }
+
+ void set(int i, int j, const fMatrix &m)
+ {
+ BLI_assert(i >= 0 && i < m_numverts);
+ BLI_assert(j >= 0 && j < m_numverts);
+ i *= 3;
+ j *= 3;
+ for (int k = 0; k < 3; ++k)
+ for (int l = 0; l < 3; ++l)
+ m_trips.push_back(Triplet(i + k, j + l,
m.coeff(l, k)));
+ }
+
+ inline lMatrix construct() const
+ {
+ lMatrix m(m_numverts, m_numverts);
+ m.setFromTriplets(m_trips.begin(), m_trips.end());
+ return m;
+ }
+
+private:
+ const int m_numverts;
+ TripletList m_trips;
+};
+
#ifdef USE_EIGEN_CORE
typedef Eigen::ConjugateGradient<lMatrix, Eigen::Lower,
Eigen::DiagonalPreconditioner<Scalar> > ConjugateGradient;
#endif
@@ -124,9 +232,6 @@ static void print_lmatrix(const lMatrix &m)
}
}
-static float I[3][3] = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
-static float ZERO[3][3] = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
-
BLI_INLINE void lMatrix_reserve_elems(lMatrix &m, int num)
{
m.reserve(Eigen::VectorXi::Constant(m.cols(), num));
@@ -185,40 +290,6 @@ BLI_INLINE void lMatrix_sub_triplets(lMatrix &r, const
TripletList &tlist)
r -= t;
}
-#if 0
-BLI_INLINE void lMatrix_copy_m3(lMatrix &r, float m[3][3], int i, int j)
-{
- i *= 3;
- j *= 3;
- for (int l = 0; l < 3; ++l) {
- for (int k = 0; k < 3; ++k) {
- r.coeffRef(i + k, j + l) = m[k][l];
- }
- }
-}
-
-BLI_INLINE void lMatrix_add_m3(lMatrix &r, float m[3][3], int i, int j)
-{
- lMatrix tmp(r.cols(), r.cols());
- lMatrix_copy_m3(tmp, m, i, j);
- r += tmp;
-}
-
-BLI_INLINE void lMatrix_sub_m3(lMatrix &r, float m[3][3], int i, int j)
-{
- lMatrix tmp(r.cols(), r.cols());
- lMatrix_copy_m3(tmp, m, i, j);
- r -= tmp;
-}
-
-BLI_INLINE void lMatrix_madd_m3(lMatrix &r, float m[3][3], float s, int i, int
j)
-{
- lMatrix tmp(r.cols(), r.cols());
- lMatrix_copy_m3(tmp, m, i, j);
- r += s * tmp;
-}
-#endif
-
BLI_INLINE void outerproduct(float r[3][3], const float a[3], const float b[3])
{
mul_v3_v3fl(r[0], a, b[0]);
@@ -226,21 +297,11 @@ BLI_INLINE void outerproduct(float r[3][3], const float
a[3], const float b[3])
mul_v3_v3fl(r[2], a, b[2]);
}
-struct RootTransform {
- float loc[3];
- float rot[3][3];
-
- float vel[3];
- float omega[3];
-
- float acc[3];
- float domega_dt[3];
-};
-
struct Implicit_Data {
- typedef std::vector<RootTransform> RootTransforms;
+ typedef std::vector<fMatrix> fMatrixVector;
- Implicit_Data(int numverts)
+ Implicit_Data(int numverts) :
+ M(numverts)
{
resize(numverts);
}
@@ -250,11 +311,10 @@ struct Implicit_Data {
this->numverts = numverts;
int tot = 3 * numverts;
- M.resize(tot, tot);
dFdV.resize(tot, tot);
dFdX.resize(tot, tot);
- root.resize(numverts);
+ tfm.resize(numverts, I);
X.resize(tot);
Xnew.resize(tot);
@@ -273,11 +333,11 @@ struct Implicit_Data {
int numverts;
/* inputs */
- lMatrix M; /* masses */
+ lMatrixCtor M; /* masses */
lVector F; /* forces */
lMatrix dFdV, dFdX; /* force jacobians */
- RootTransforms root; /* root transforms */
+ fMatrixVector tfm; /* local coordinate transform */
/* motion state data */
lVector X, Xnew; /* positions */
@@ -290,182 +350,39 @@ struct Implicit_Data {
lVector dV; /* velocity change
(solution of A*dV = B) */
lVector z; /* target velocity in
constrained directions */
lMatrix S; /* filtering matrix for
constraints */
+
+ struct SimDebugData *debug_data;
};
-/* ==== Transformation of Moving Reference Frame ====
- * x_world, v_world, f_world, a_world, dfdx_world, dfdv_world : state
variables in world space
- * x_root, v_root, f_root, a_root, dfdx_root, dfdv_root : state
variables in root space
- *
- * x0 : translation of the root frame (hair root location)
- * v0 : linear velocity of the root frame
- * a0 : acceleration of the root frame
- * R : rotation matrix of the root frame
- * w : angular velocity of the root frame
- * dwdt : angular acceleration of the root frame
- */
-
-/* x_root = R^T * x_world */
-BLI_INLINE void loc_world_to_root(float r[3], const float v[3], const
RootTransform &root)
-{
- sub_v3_v3v3(r, v, root.loc);
- mul_transposed_m3_v3((float (*)[3])root.rot, r);
-}
+/* ==== Transformation from/to root reference frames ==== */
-/* x_world = R * x_root */
-BLI_INLINE void loc_root_to_world(float r[3], const float v[3], const
RootTransform &root)
+BLI_INLINE void world_to_root_v3(Implicit_Data *data, int index, float r[3],
const float v[3])
{
copy_v3_v3(r, v);
- mul_m3_v3((float (*)[3])root.rot, r);
- add_v3_v3(r, root.loc);
-}
-
-/* v_root = cross(w, x_root) + R^T*(v_world - v0) */
-BLI_INLINE void vel_world_to_root(float r[3], const float x_root[3], const
float v[3], const RootTransform &root)
-{
- float angvel[3];
- cross_v3_v3v3(angvel, root.omega, x_root);
-
- sub_v3_v3v3(r, v, root.vel);
- mul_transposed_m3_v3((float (*)[3])root.rot, r);
- add_v3_v3(r, angvel);
-}
-
-/* v_world = R*(v_root - cross(w, x_root)) + v0 */
-BLI_INLINE void vel_root_to_world(float r[3], const float x_root[3], const
float v[3], const RootTransform &root)
-{
- float angvel[3];
- cross_v3_v3v3(angvel, root.omega, x_root);
-
- sub_v3_v3v3(r, v, angvel);
- mul_m3_v3((float (*)[3])root.rot, r);
- add_v3_v3(r, root.vel);
-}
-
-/* a_root = -cross(dwdt, x_root) - 2*cross(w, v_root) - cross(w, cross(w,
x_root)) + R^T*(a_world - a0) */
-BLI_INLINE void force_world_to_root(float r[3], const float x_root[3], const
float v_root[3], const float force[3], float mass, const RootTransform &root)
-{
- float euler[3], coriolis[3], centrifugal[3], rotvel[3];
-
- cross_v3_v3v3(euler, root.domega_dt, x_root);
- cross_v3_v3v3(coriolis, root.omega, v_root);
- mul_v3_fl(coriolis, 2.0f);
- cross_v3_v3v3(rotvel, root.omega, x_root);
- cross_v3_v3v3(centrifugal, root.omega, rotvel);
-
- madd_v3_v3v3fl(r, force, root.acc, mass);
- mul_transposed_m3_v3((float (*)[3])root.rot, r);
- madd_v3_v3fl(r, euler, mass);
- madd_v3_v3fl(r, coriolis, mass);
- madd_v3_v3fl(r, centrifugal, mass);
-}
-
-/* a_world = R*[ a_root + cross(dwdt, x_root) + 2*cross(w, v_root) + cross(w,
cross(w, x_root)) ] + a0 */
-BLI_INLINE void force_root_to_world(float r[3], const float x_root[3], const
float v_root[3], const float force[3], float mass, const RootTransform &root)
-{
- float euler[3], coriolis[3], centrifugal[3], rotvel[3];
-
- cross_v3_v3v3(euler, root.domega_dt, x_root);
- cross_v3_v3v3(coriolis, root.omega, v_root);
- mul_v3_fl(coriolis, 2.0f);
- cross_v3_v3v3(rotvel, root.omega, x_root);
- cross_v3_v3v3(centrifugal, root.omega, rotvel);
-
- madd_v3_v3v3fl(r, force, euler, mass);
- madd_v3_v3fl(r, coriolis, mass);
- madd_v3_v3fl(r, centrifugal, mass);
- mul_m3_v3((float (*)[3])root.rot, r);
- madd_v3_v3fl(r, root.acc, mass);
-}
-
-BLI_INLINE void acc_world_to_root(float r[3], const float x_root[3], const
float v_root[3], const float acc[3], const RootTransform &root)
-{
- force_world_to_root(r, x_root, v_root, acc, 1.0f, root);
-}
-
-BLI_INLINE void acc_root_to_world(float r[3], const float x_root[3], const
float v_root[3], const float acc[3], const RootTransform &root)
-{
- force_root_to_world(r, x_root, v_root, acc, 1.0f, root);
+ mul_transposed_m3_v3(data->tfm[index], r);
}
-BLI_INLINE void cross_m3_v3m3(float r[3][3], const float v[3], float m[3][3])
+BLI_INLINE void root_to_world_v3(Implicit_Data *data, int index, float r[3],
const float v[3])
{
- cross_v3_v3v3(r[0], v, m[0]);
- cross_v3_v3v3(r[1], v, m[1]);
- cross_v3_v3v3(r[2], v, m[2]);
+ mul_v3_m
@@ Diff output truncated at 10240 characters. @@
_______________________________________________
Bf-blender-cvs mailing list
[email protected]
http://lists.blender.org/mailman/listinfo/bf-blender-cvs
