[Bf-blender-cvs] [854600a] hair_immediate_fixes: Alternative new solver for cloth using the Eigen CG solver instead of a custom built solver.
Lukas Tönne
noreply at git.blender.org
Fri Sep 5 12:50:45 CEST 2014
Commit: 854600aeed5d53736b8e009ae08c8872875aaa25
Author: Lukas Tönne
Date: Fri Sep 5 12:48:49 2014 +0200
Branches: hair_immediate_fixes
https://developer.blender.org/rB854600aeed5d53736b8e009ae08c8872875aaa25
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 0caf7d1..460333c 100644
--- a/source/blender/blenkernel/CMakeLists.txt
+++ b/source/blender/blenkernel/CMakeLists.txt
@@ -43,8 +43,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
@@ -103,6 +104,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 ae85579..14a26ec 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"
@@ -2477,3 +2480,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