Commit: e57d3af847d6f618ab4e03d80b341c7eb9da0db6
Author: Lukas Tönne
Date:   Sun Oct 12 16:14:15 2014 +0200
Branches: gooseberry
https://developer.blender.org/rBe57d3af847d6f618ab4e03d80b341c7eb9da0db6

Ported the remaining implicit solver functions for Eigen.

Also added a couple of utility wrapper functions for Eigen types to make
interfacing with plain float arrays and blenlib math easier.

===================================================================

M       source/blender/physics/intern/BPH_mass_spring.cpp
M       source/blender/physics/intern/implicit.h
M       source/blender/physics/intern/implicit_blender.c
M       source/blender/physics/intern/implicit_eigen.cpp

===================================================================

diff --git a/source/blender/physics/intern/BPH_mass_spring.cpp 
b/source/blender/physics/intern/BPH_mass_spring.cpp
index ceb76af..875b7bf 100644
--- a/source/blender/physics/intern/BPH_mass_spring.cpp
+++ b/source/blender/physics/intern/BPH_mass_spring.cpp
@@ -620,7 +620,10 @@ static void cloth_calc_force(ClothModifierData *clmd, 
float UNUSED(frame), ListB
                /* scale gravity force */
                mul_v3_v3fl(gravity, clmd->scene->physics_settings.gravity, 
0.001f * clmd->sim_parms->effector_weights->global_gravity);
        }
-       BPH_mass_spring_force_gravity(data, gravity);
+       vert = cloth->verts;
+       for (i = 0; i < cloth->numverts; i++, vert++) {
+               BPH_mass_spring_force_gravity(data, i, vert->mass, gravity);
+       }
 #endif
 
        cloth_calc_volume_force(clmd);
diff --git a/source/blender/physics/intern/implicit.h 
b/source/blender/physics/intern/implicit.h
index 4659bb0..287c064 100644
--- a/source/blender/physics/intern/implicit.h
+++ b/source/blender/physics/intern/implicit.h
@@ -59,6 +59,7 @@ extern "C" {
 //#define IMPLICIT_ENABLE_EIGEN_DEBUG
 
 struct Implicit_Data;
+struct ImplicitSolverInput;
 struct SimDebugData;
 
 typedef struct ImplicitSolverResult {
@@ -118,8 +119,6 @@ void BPH_mass_spring_set_velocity(struct Implicit_Data 
*data, int index, const f
 void BPH_mass_spring_get_motion_state(struct Implicit_Data *data, int index, 
float x[3], float v[3]);
 void BPH_mass_spring_set_vertex_mass(struct Implicit_Data *data, int index, 
float mass);
 
-int BPH_mass_spring_add_block(struct Implicit_Data *data, int v1, int v2);
-
 void BPH_mass_spring_clear_constraints(struct Implicit_Data *data);
 void BPH_mass_spring_add_constraint_ndof0(struct Implicit_Data *data, int 
index, const float dV[3]);
 void BPH_mass_spring_add_constraint_ndof1(struct Implicit_Data *data, int 
index, const float c1[3], const float c2[3], const float dV[3]);
@@ -131,9 +130,9 @@ void BPH_mass_spring_apply_result(struct Implicit_Data 
*data);
 /* Clear the force vector at the beginning of the time step */
 void BPH_mass_spring_clear_forces(struct Implicit_Data *data);
 /* Fictitious forces introduced by moving coordinate systems */
-void BPH_mass_spring_force_reference_frame(struct Implicit_Data *data, int 
index, const float acceleration[3], const float omega[3], const float 
domega_dt[3]);
+void BPH_mass_spring_force_reference_frame(struct Implicit_Data *data, int 
index, const float acceleration[3], const float omega[3], const float 
domega_dt[3], float mass);
 /* Simple uniform gravity force */
-void BPH_mass_spring_force_gravity(struct Implicit_Data *data, const float 
g[3]);
+void BPH_mass_spring_force_gravity(struct Implicit_Data *data, int index, 
float mass, const float g[3]);
 /* Global drag force (velocity damping) */
 void BPH_mass_spring_force_drag(struct Implicit_Data *data, float drag);
 /* Custom external force */
diff --git a/source/blender/physics/intern/implicit_blender.c 
b/source/blender/physics/intern/implicit_blender.c
index 741aeb0..7421a7a 100644
--- a/source/blender/physics/intern/implicit_blender.c
+++ b/source/blender/physics/intern/implicit_blender.c
@@ -925,7 +925,8 @@ int BPH_mass_spring_solver_numvert(Implicit_Data *id)
 
 void BPH_mass_spring_solver_debug_data(Implicit_Data *id, struct SimDebugData 
*debug_data)
 {
-       id->debug_data = debug_data;
+       if (id)
+               id->debug_data = debug_data;
 }
 
 /* ==== Transformation from/to root reference frames ==== */
@@ -1386,7 +1387,7 @@ void BPH_mass_spring_set_vertex_mass(Implicit_Data *data, 
int index, float mass)
        mul_m3_fl(data->M[index].m, mass);
 }
 
-int BPH_mass_spring_add_block(Implicit_Data *data, int v1, int v2)
+static int BPH_mass_spring_add_block(Implicit_Data *data, int v1, int v2)
 {
        int s = data->M[0].vcount + data->num_blocks; /* index from array start 
*/
        BLI_assert(s < data->M[0].vcount + data->M[0].scount);
@@ -1465,7 +1466,7 @@ void BPH_mass_spring_clear_forces(Implicit_Data *data)
        data->num_blocks = 0;
 }
 
-void BPH_mass_spring_force_reference_frame(Implicit_Data *data, int index, 
const float acceleration[3], const float omega[3], const float domega_dt[3])
+void BPH_mass_spring_force_reference_frame(Implicit_Data *data, int index, 
const float acceleration[3], const float omega[3], const float domega_dt[3], 
float mass)
 {
 #ifdef CLOTH_ROOT_FRAME
        float acc[3], w[3], dwdt[3];
@@ -1487,7 +1488,7 @@ void BPH_mass_spring_force_reference_frame(Implicit_Data 
*data, int index, const
        sub_v3_v3(f, coriolis);
        sub_v3_v3(f, centrifugal);
        
-       mul_m3_v3(data->M[index].m, f); /* F = m * a */
+       mul_v3_fl(f, mass); /* F = m * a */
        
        cross_v3_identity(deuler, dwdt);
        cross_v3_identity(dcoriolis, w);
@@ -1497,11 +1498,11 @@ void 
BPH_mass_spring_force_reference_frame(Implicit_Data *data, int index, const
        
        add_m3_m3m3(dfdx, deuler, dcentrifugal);
        negate_m3(dfdx);
-       mul_m3_m3m3(dfdx, data->M[index].m, dfdx);
+       mul_m3_fl(dfdx, mass);
        
        copy_m3_m3(dfdv, dcoriolis);
        negate_m3(dfdv);
-       mul_m3_m3m3(dfdv, data->M[index].m, dfdv);
+       mul_m3_fl(dfdv, mass);
        
        add_v3_v3(data->F[index], f);
        add_m3_m3m3(data->dFdX[index].m, data->dFdX[index].m, dfdx);
@@ -1515,19 +1516,14 @@ void 
BPH_mass_spring_force_reference_frame(Implicit_Data *data, int index, const
 #endif
 }
 
-void BPH_mass_spring_force_gravity(Implicit_Data *data, const float g[3])
+void BPH_mass_spring_force_gravity(Implicit_Data *data, int index, float mass, 
const float g[3])
 {
-       int i, numverts = data->M[0].vcount;
-       /* multiply F with mass matrix
-        * force = mass * acceleration (in this case: gravity)
-        */
-       for (i = 0; i < numverts; i++) {
-               float f[3];
-               world_to_root_v3(data, i, f, g);
-               mul_m3_v3(data->M[i].m, f);
-               
-               add_v3_v3(data->F[i], f);
-       }
+       /* force = mass * acceleration (in this case: gravity) */
+       float f[3];
+       world_to_root_v3(data, index, f, g);
+       mul_v3_fl(f, mass);
+       
+       add_v3_v3(data->F[index], f);
 }
 
 void BPH_mass_spring_force_drag(Implicit_Data *data, float drag)
diff --git a/source/blender/physics/intern/implicit_eigen.cpp 
b/source/blender/physics/intern/implicit_eigen.cpp
index 402ffcb..a6148b6 100644
--- a/source/blender/physics/intern/implicit_eigen.cpp
+++ b/source/blender/physics/intern/implicit_eigen.cpp
@@ -111,7 +111,7 @@ public:
                        coeffRef(k) = v[k];
        }
        
-       fVector &operator = (const ctype &v)
+       fVector& operator = (const ctype &v)
        {
                for (int k = 0; k < 3; ++k)
                        coeffRef(k) = v[k];
@@ -142,7 +142,7 @@ public:
                                coeffRef(l, k) = v[k][l];
        }
        
-       fMatrix &operator = (const ctype &v)
+       fMatrix& operator = (const ctype &v)
        {
                for (int k = 0; k < 3; ++k)
                        for (int l = 0; l < 3; ++l)
@@ -156,27 +156,63 @@ public:
        }
 };
 
-typedef Eigen::VectorXf lVector;
+/* Extension of dense Eigen vectors,
+ * providing 3-float block access for blenlib math functions
+ */
+class lVector : public Eigen::VectorXf {
+public:
+       typedef Eigen::VectorXf base_t;
+       
+       lVector()
+       {
+       }
+       
+       template <typename T>
+       lVector& operator = (T rhs)
+       {
+               base_t::operator=(rhs);
+               return *this;
+       }
+       
+       float* v3(int vertex)
+       {
+               return &coeffRef(3 * vertex);
+       }
+       
+       const float* v3(int vertex) const
+       {
+               return &coeffRef(3 * vertex);
+       }
+};
 
 typedef Eigen::Triplet<Scalar> Triplet;
 typedef std::vector<Triplet> TripletList;
 
 typedef Eigen::SparseMatrix<Scalar> lMatrix;
 
+/* Constructor type that provides more convenient handling of Eigen triplets
+ * for efficient construction of sparse 3x3 block matrices.
+ * This should be used for building lMatrix instead of writing to such lMatrix 
directly (which is very inefficient).
+ * After all elements have been defined using the set() method, the actual 
matrix can be filled using construct().
+ */
 struct lMatrixCtor {
-       lMatrixCtor(int numverts) :
-           m_numverts(numverts)
+       lMatrixCtor()
+       {
+       }
+       
+       void reset()
+       {
+               m_trips.clear();
+       }
+       
+       void reserve(int 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)
+       void add(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)
@@ -184,15 +220,22 @@ struct lMatrixCtor {
                                m_trips.push_back(Triplet(i + k, j + l, 
m.coeff(l, k)));
        }
        
-       inline lMatrix construct() const
+       void sub(int i, int j, const fMatrix &m)
+       {
+               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 void construct(lMatrix &m)
        {
-               lMatrix m(m_numverts, m_numverts);
                m.setFromTriplets(m_trips.begin(), m_trips.end());
-               return m;
+               m_trips.clear();
        }
        
 private:
-       const int m_numverts;
        TripletList m_trips;
 };
 
@@ -247,6 +290,7 @@ BLI_INLINE const float *lVector_v3(const lVector &v, int 
vertex)
        return v.data() + 3 * vertex;
 }
 
+#if 0
 BLI_INLINE void triplets_m3(TripletList &tlist, float m[3][3], int i, int j)
 {
        i *= 3;
@@ -289,6 +333,7 @@ BLI_INLINE void lMatrix_sub_triplets(lMatrix &r, const 
TripletList &tlist)
        t.setFromTriplets(tlist.begin(), tlist.end());
        r -= t;
 }
+#endif
 
 BLI_INLINE void outerproduct(float r[3][3], const float a[3], const float b[3])
 {
@@ -297,11 +342,50 @@ BLI_INLINE void outerproduct(float r[3][3], const float 
a[3], const float b[3])
        mul_v3_v3fl(r[2], a, b[2]);
 }
 
+BLI_INLINE void cross_m3_v3m3(float r[3][3], const float v[3], float m[3][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]);
+}
+
+BLI_INLINE void cross_v3_identity(float r[3][3], const float v[3])
+{
+       r[0][0] = 0.0f;         r[1][0] = v[2];         r[2][0] = -v[1];
+       r[0][1] = -v[2];        r[1][1] = 0.0f;         r[2][1] = v[0];
+       r[0][2] = v[1];         r[1][2] = -v[0];        r[2][2] = 0.0f;
+}
+
+BLI_INLINE void madd_m3_m3fl(float r[3][3], float m[3][3], float f)
+{
+       r[0][0] += m[0][0] * f;
+       r[0][1] += m[0][1] * f;
+       r[0][2] += m[0][2] * f;
+       r[1][0] += m[1][0] * f;
+       r[1][1] += m[1][1] * f;
+       r[1][2] += m[1][2] * f;
+       r[2][0] += m[2][0] * f;
+       r[2][1] += m[2][1] * f;
+       r[2][2] += m[2][2] * f;
+}
+
+BLI_INLINE void madd_m3_m3m3fl(float r[3][3], float a[3][3], float b[3][3], 
float f)
+{
+       r[0][0] = a[0][0] + b[0][0] * f;
+       r[0][1] = a[0][1] + b[0][1] * f;


@@ Diff output truncated at 10240 characters. @@

_______________________________________________
Bf-blender-cvs mailing list
[email protected]
http://lists.blender.org/mailman/listinfo/bf-blender-cvs

Reply via email to