[Bf-blender-cvs] [3af3a2c] temp_merge_gooseberry_hair: Alternative new solver for cloth using the Eigen CG solver instead of a custom built solver.
Lukas Tönne
noreply at git.blender.org
Mon Jan 19 20:48:15 CET 2015
Commit: 3af3a2c9afbcb9d089d5f5a3fe5e40614728985f
Author: Lukas Tönne
Date: Fri Sep 5 12:48:49 2014 +0200
Branches: temp_merge_gooseberry_hair
https://developer.blender.org/rB3af3a2c9afbcb9d089d5f5a3fe5e40614728985f
Alternative new solver for cloth using the Eigen CG solver instead of
a custom built solver.
The old cloth solver is broken unfortunately. Eigen is a designated
linear algebra library and very likely their implementation is a lot
better (can't compare until it's implemented though).
Only basic gravity is active atm, spring forces, external force fields,
damping and volumetric friction have to be added back by converting
the data into the Eigen format.
===================================================================
M source/blender/blenkernel/CMakeLists.txt
M source/blender/blenkernel/intern/implicit.c
A source/blender/blenkernel/intern/implicit_eigen.cpp
===================================================================
diff --git a/source/blender/blenkernel/CMakeLists.txt b/source/blender/blenkernel/CMakeLists.txt
index fb717b1..098b118 100644
--- a/source/blender/blenkernel/CMakeLists.txt
+++ b/source/blender/blenkernel/CMakeLists.txt
@@ -45,8 +45,9 @@ set(INC
../../../intern/mikktspace
../../../intern/raskter
../../../intern/smoke/extern
- ../../../extern/libmv
../../../intern/atomic
+ ../../../extern/libmv
+ ../../../extern/Eigen3
# XXX - BAD LEVEL CALL WM_api.h
../windowmanager
@@ -106,6 +107,7 @@ set(SRC
intern/image.c
intern/image_gen.c
intern/implicit.c
+ intern/implicit_eigen.cpp
intern/ipo.c
intern/key.c
intern/lamp.c
diff --git a/source/blender/blenkernel/intern/implicit.c b/source/blender/blenkernel/intern/implicit.c
index 227bd2a..4081a60 100644
--- a/source/blender/blenkernel/intern/implicit.c
+++ b/source/blender/blenkernel/intern/implicit.c
@@ -29,6 +29,9 @@
* \ingroup bke
*/
+//#define IMPLICIT_SOLVER_BLENDER
+
+#ifdef IMPLICIT_SOLVER_BLENDER
#include "MEM_guardedalloc.h"
@@ -2283,3 +2286,5 @@ void implicit_set_positions(ClothModifierData *clmd)
if (G.debug_value > 0)
printf("implicit_set_positions\n");
}
+
+#endif /* IMPLICIT_SOLVER_BLENDER */
diff --git a/source/blender/blenkernel/intern/implicit_eigen.cpp b/source/blender/blenkernel/intern/implicit_eigen.cpp
new file mode 100644
index 0000000..9bf31a3
--- /dev/null
+++ b/source/blender/blenkernel/intern/implicit_eigen.cpp
@@ -0,0 +1,999 @@
+/*
+ * ***** BEGIN GPL LICENSE BLOCK *****
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ *
+ * The Original Code is Copyright (C) Blender Foundation
+ * All rights reserved.
+ *
+ * The Original Code is: all of this file.
+ *
+ * Contributor(s): Lukas Toenne
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ */
+
+/** \file blender/blenkernel/intern/implicit_eigen.cpp
+ * \ingroup bke
+ */
+
+#define IMPLICIT_SOLVER_EIGEN
+
+#ifdef IMPLICIT_SOLVER_EIGEN
+
+#include <Eigen/Sparse>
+
+#include "MEM_guardedalloc.h"
+
+extern "C" {
+#include "DNA_scene_types.h"
+#include "DNA_object_types.h"
+#include "DNA_object_force.h"
+#include "DNA_meshdata_types.h"
+#include "DNA_texture_types.h"
+
+#include "BLI_math.h"
+#include "BLI_linklist.h"
+#include "BLI_utildefines.h"
+
+#include "BKE_cloth.h"
+#include "BKE_collision.h"
+#include "BKE_effect.h"
+#include "BKE_global.h"
+}
+
+
+/* ==== hash functions for debugging ==== */
+static unsigned int hash_int_2d(unsigned int kx, unsigned int ky)
+{
+#define rot(x,k) (((x)<<(k)) | ((x)>>(32-(k))))
+
+ unsigned int a, b, c;
+
+ a = b = c = 0xdeadbeef + (2 << 2) + 13;
+ a += kx;
+ b += ky;
+
+ c ^= b; c -= rot(b,14);
+ a ^= c; a -= rot(c,11);
+ b ^= a; b -= rot(a,25);
+ c ^= b; c -= rot(b,16);
+ a ^= c; a -= rot(c,4);
+ b ^= a; b -= rot(a,14);
+ c ^= b; c -= rot(b,24);
+
+ return c;
+
+#undef rot
+}
+
+static int hash_vertex(int type, int vertex)
+{
+ return hash_int_2d((unsigned int)type, (unsigned int)vertex);
+}
+
+static int hash_collpair(int type, CollPair *collpair)
+{
+ return hash_int_2d((unsigned int)type, hash_int_2d((unsigned int)collpair->face1, (unsigned int)collpair->face2));
+}
+/* ================ */
+
+
+typedef float Scalar;
+
+typedef Eigen::SparseMatrix<Scalar> lMatrix;
+
+typedef Eigen::VectorXf lVector;
+
+typedef Eigen::Triplet<Scalar> Triplet;
+typedef std::vector<Triplet> TripletList;
+
+typedef Eigen::ConjugateGradient<lMatrix, Eigen::Lower, Eigen::DiagonalPreconditioner<Scalar> > ConjugateGradient;
+using Eigen::ComputationInfo;
+
+
+BLI_INLINE void reserve_diagonal(lMatrix &m, int num)
+{
+ m.reserve(Eigen::VectorXi::Constant(m.cols(), num));
+}
+
+BLI_INLINE float *lVector_v3(lVector &v, int index)
+{
+ return v.data() + 3 * index;
+}
+
+BLI_INLINE const float *lVector_v3(const lVector &v, int index)
+{
+ return v.data() + 3 * index;
+}
+
+struct Implicit_Data {
+ Implicit_Data(int numverts)
+ {
+ resize(numverts);
+ }
+
+ void resize(int numverts)
+ {
+ this->numverts = numverts;
+ int tot = 3 * numverts;
+
+ M.resize(tot, tot);
+ dFdV.resize(tot, tot);
+ dFdX.resize(tot, tot);
+
+ X.resize(tot);
+ Xnew.resize(tot);
+ V.resize(tot);
+ Vnew.resize(tot);
+ F.resize(tot);
+
+ B.resize(tot);
+ A.resize(tot, tot);
+
+ dV.resize(tot);
+ z.resize(tot);
+ S.resize(tot, tot);
+ }
+
+ int numverts;
+
+ /* inputs */
+ lMatrix M; /* masses */
+ lMatrix dFdV, dFdX; /* force jacobians */
+
+ /* motion state data */
+ lVector X, Xnew; /* positions */
+ lVector V, Vnew; /* velocities */
+ lVector F; /* forces */
+
+ /* internal solver data */
+ lVector B; /* B for A*dV = B */
+ lMatrix A; /* A for A*dV = B */
+
+ lVector dV; /* velocity change (solution of A*dV = B) */
+ lVector z; /* target velocity in constrained directions */
+ lMatrix S; /* filtering matrix for constraints */
+};
+
+static bool simulate_implicit_euler(Implicit_Data *id, float dt)
+{
+ ConjugateGradient cg;
+ cg.setMaxIterations(100);
+ cg.setTolerance(0.01f);
+
+ id->A = id->M - dt * id->dFdV - dt*dt * id->dFdX;
+ cg.compute(id->A);
+
+ id->B = dt * id->F + dt*dt * id->dFdX * id->V;
+ id->dV = cg.solve(id->B);
+
+ id->Vnew = id->V + id->dV;
+
+ return cg.info() != Eigen::Success;
+}
+
+static void cloth_calc_force(ClothModifierData *clmd, lVector &F, lMatrix &dFdX, lMatrix &dFdV, const lVector &X, const lVector &V, const lMatrix &M, ListBase *effectors, float time)
+{
+ Cloth *cloth = clmd->clothObject;
+ unsigned int numverts = cloth->numverts;
+ float spring_air = clmd->sim_parms->Cvi * 0.01f; /* viscosity of air scaled in percent */
+ float gravity[3] = {0,0,0};
+// lVector diagonal(3*numverts);
+// diagonal.setZero
+ TripletList coeff_dFdX, coeff_dFdV;
+
+ /* set dFdX jacobi matrix to zero */
+ dFdX.setZero();
+ /* set dFdV jacobi matrix diagonal entries to -spring_air */
+ dFdV.setZero();
+// dFdV.setIdentity();
+// initdiag_bfmatrix(dFdV, tm2);
+
+ /* global acceleration (gravitation) */
+ if (clmd->scene->physics_settings.flag & PHYS_GLOBAL_GRAVITY) {
+ /* scale gravity force
+ * XXX this factor looks totally arbitrary ... what is this? lukas_t
+ */
+ mul_v3_v3fl(gravity, clmd->scene->physics_settings.gravity, 0.001f * clmd->sim_parms->effector_weights->global_gravity);
+ }
+
+ for (int i = 0; i < numverts; ++i) {
+ copy_v3_v3(lVector_v3(F, i), gravity);
+ }
+
+#if 0
+ /* Collect forces and derivatives: F, dFdX, dFdV */
+ Cloth *cloth = clmd->clothObject;
+ unsigned int i = 0;
+ float spring_air = clmd->sim_parms->Cvi * 0.01f; /* viscosity of air scaled in percent */
+ float gravity[3] = {0.0f, 0.0f, 0.0f};
+ float tm2[3][3] = {{0}};
+ MFace *mfaces = cloth->mfaces;
+ unsigned int numverts = cloth->numverts;
+ LinkNode *search;
+ lfVector *winvec;
+ EffectedPoint epoint;
+
+ tm2[0][0] = tm2[1][1] = tm2[2][2] = -spring_air;
+
+ /* global acceleration (gravitation) */
+ if (clmd->scene->physics_settings.flag & PHYS_GLOBAL_GRAVITY) {
+ copy_v3_v3(gravity, clmd->scene->physics_settings.gravity);
+ mul_fvector_S(gravity, gravity, 0.001f * clmd->sim_parms->effector_weights->global_gravity); /* scale gravity force */
+ }
+
+ /* set dFdX jacobi matrix to zero */
+ init_bfmatrix(dFdX, ZERO);
+ /* set dFdX jacobi matrix diagonal entries to -spring_air */
+ initdiag_bfmatrix(dFdV, tm2);
+
+ init_lfvector(lF, gravity, numverts);
+
+ hair_volume_forces(clmd, lF, lX, lV, numverts);
+
+ /* multiply lF with mass matrix
+ * force = mass * acceleration (in this case: gravity)
+ */
+ for (i = 0; i < numverts; i++) {
+ float temp[3];
+ copy_v3_v3(temp, lF[i]);
+ mul_fmatrix_fvector(lF[i], M[i].m, temp);
+ }
+
+ submul_lfvectorS(lF, lV, spring_air, numverts);
+
+ /* handle external forces like wind */
+ if (effectors) {
+ // 0 = force, 1 = normalized force
+ winvec = create_lfvector(cloth->numverts);
+
+ if (!winvec)
+ printf("winvec: out of memory in implicit.c\n");
+
+ // precalculate wind forces
+ for (i = 0; i < cloth->numverts; i++) {
+ pd_point_from_loc(clmd->scene, (float*)lX[i], (float*)lV[i], i, &epoint);
+ pdDoEffectors(effectors, NULL, clmd->sim_parms->effector_weights, &epoint, winvec[i], NULL);
+ }
+
+ for (i = 0; i < cloth->numfaces; i++) {
+ float trinormal[3] = {0, 0, 0}; // normalized triangle normal
+ float triunnormal[3] = {0, 0, 0}; // not-normalized-triangle normal
+ float tmp[3] = {0, 0, 0};
+ float factor = (mfaces[i].v4) ? 0.25 : 1.0 / 3.0;
+ factor *= 0.02f;
+
+ // calculate face normal
+ if (mfaces[i].v4)
+ CalcFloat4(lX[mfaces[i].v1], lX[mfaces[i].v2], lX[mfaces[i].v3], lX[mfaces[i].v4], triunnormal);
+ else
+ CalcFloat(lX[mfaces[i].v1], lX[mfaces[i].v2], lX[mfaces[i].v3], triunnormal);
+
+ normalize_v3_v3(trinormal, triunnormal);
+
+ // add wind from v1
+ copy_v3_v3(tmp, trinormal);
+ mul_v3_fl(tmp, calculateVertexWindForce(winvec[mfaces[i].v1], triunnormal));
+ VECADDS(lF[mfaces[i].v1], lF[mfaces[i].v1], tmp, factor);
+
+ // add wind from v2
+ copy_v3_v3(tmp, trinormal);
+ mul_v3_fl(tmp, calculateVertexWindForce(winvec[mfaces[i].v2], triunnormal));
+ VECADDS(lF[mfaces[i].v2], lF[mfaces[i].v2], tmp, factor);
+
+ // add wind from v3
+ copy_v3_v3(tmp, trinormal);
+ mul_v3_fl(tmp, calculateVertexWindForce(winvec[mfaces[i].v3], triunnormal));
+ VECADDS(lF[mfaces[i].v3], lF[mfaces[i].v3], tmp, factor);
+
+ // add wind from v4
+ if (mfaces[i].v4) {
+ copy_v3_v3(tmp, trinormal);
+ mul_v3_fl(tmp, calculateVertexWindForce(winvec[mfaces[i].v4], triunnormal));
+ VECADDS(lF[mfaces[i].v4], lF[mfaces[i].v4], tmp, factor);
+
@@ Diff output truncated at 10240 characters. @@
More information about the Bf-blender-cvs
mailing list