[Bf-blender-cvs] [e303e23] strand_nodes: Implementation of a basic constraint solver for strand length and root velocity.
Lukas Tönne
noreply at git.blender.org
Sun Aug 7 10:39:52 CEST 2016
Commit: e303e23d56200c43fd4d06c6aece309cfd9a6fd2
Author: Lukas Tönne
Date: Sun Aug 7 10:35:48 2016 +0200
Branches: strand_nodes
https://developer.blender.org/rBe303e23d56200c43fd4d06c6aece309cfd9a6fd2
Implementation of a basic constraint solver for strand length and root velocity.
This uses the Lagrange multiplier method for correcting the "velocity" of
vertex displacements as described in
"Linear-Time Dynamics using Lagrange Multipliers" (Baraff, 1996)
Future additions would include stiffness forces for creating more natural
hair shapes, and handling of collisions (deflection).
===================================================================
M source/blender/blenkernel/BKE_editstrands.h
M source/blender/editors/hair/hair_edit.c
M source/blender/physics/intern/eigen_utils.h
M source/blender/physics/intern/strands.cpp
===================================================================
diff --git a/source/blender/blenkernel/BKE_editstrands.h b/source/blender/blenkernel/BKE_editstrands.h
index 82258a3..a028308 100644
--- a/source/blender/blenkernel/BKE_editstrands.h
+++ b/source/blender/blenkernel/BKE_editstrands.h
@@ -35,6 +35,7 @@
#include "BLI_utildefines.h"
#include "DNA_customdata_types.h"
+#include "DNA_defs.h"
#include "BKE_customdata.h"
#include "BKE_editmesh.h"
@@ -61,6 +62,8 @@ typedef struct BMEditStrands {
int totfibers;
int flag;
+
+ float (*locs)[3] DNA_DEPRECATED; /* XXX debugging, remove */
} BMEditStrands;
/* BMEditStrands->flag */
diff --git a/source/blender/editors/hair/hair_edit.c b/source/blender/editors/hair/hair_edit.c
index 47d2a58..967aa00 100644
--- a/source/blender/editors/hair/hair_edit.c
+++ b/source/blender/editors/hair/hair_edit.c
@@ -343,6 +343,8 @@ static bool hair_stroke_apply(bContext *C, wmOperator *op, PointerRNA *itemptr)
mul_mat3_m4_v3(tool_data.imat, tool_data.delta);
BMEditStrandsLocations locs = BKE_editstrands_get_locations(edit);
+ if (!edit->locs)
+ edit->locs = MEM_dupallocN(locs);
for (step = 0; step < totsteps; ++step) {
bool step_updated = hair_brush_step(&tool_data);
diff --git a/source/blender/physics/intern/eigen_utils.h b/source/blender/physics/intern/eigen_utils.h
index ebfed7c..cd5f431 100644
--- a/source/blender/physics/intern/eigen_utils.h
+++ b/source/blender/physics/intern/eigen_utils.h
@@ -38,6 +38,7 @@
# pragma GCC diagnostic ignored "-Wlogical-op"
#endif
+#include <Eigen/Dense>
#include <Eigen/Sparse>
#include <Eigen/SVD>
#include <Eigen/src/Core/util/DisableStupidWarnings.h>
diff --git a/source/blender/physics/intern/strands.cpp b/source/blender/physics/intern/strands.cpp
index 0fde614..e23bfaa 100644
--- a/source/blender/physics/intern/strands.cpp
+++ b/source/blender/physics/intern/strands.cpp
@@ -29,6 +29,8 @@
* \ingroup bke
*/
+#include "iostream"
+
extern "C" {
#include "MEM_guardedalloc.h"
@@ -420,10 +422,122 @@ static void strands_solve_inverse_kinematics(Object *ob, BMEditStrands *edit, fl
}
#endif
+#ifdef STRAND_CONSTRAINT_LAGRANGEMULT
+
+//#define DO_DEBUG
+
+/* Solve edge constraints and collisions for a single strand based on
+ * "Linear-Time Dynamics using Lagrange Multipliers" (Baraff, 1996)
+ */
+static void strand_solve(BMesh *UNUSED(bm), BMVert *root, float (*orig)[3], int numverts,
+ const Eigen::Vector3f &root_v)
+{
+ using Eigen::Vector3f;
+ using Eigen::Matrix3f;
+
+ BMIter iter;
+ BMVert *vert;
+ int k;
+
+ /* compute unconstrained velocities by 1st order differencing */
+ VectorX v(3 * numverts);
+// VectorX L(numverts);
+ BM_ITER_STRANDS_ELEM_INDEX(vert, &iter, root, BM_VERTS_OF_STRAND, k) {
+ sub_v3_v3v3(&v.coeffRef(3*k), vert->co, orig[k]);
+// L.coeffRef(k) = BM_elem_float_data_named_get(&bm->vdata, vert, CD_PROP_FLT, CD_HAIR_SEGMENT_LENGTH);
+#ifdef DO_DEBUG
+ BKE_sim_debug_data_add_line(orig[k], vert->co, 0,0,1, "hair solve", 3874, BLI_ghashutil_ptrhash(root), k);
+#endif
+ }
+
+ /* "Mass" matrix can be understood as resistance to editing changes.
+ * XXX For now just using identity, in future more interesting things could be done here.
+ */
+// MatrixX M = MatrixX::Identity(3 * numverts, 3 * numverts);
+ /* XXX we actually only need the inverse of M, here just skip a pointless solve step */
+ MatrixX M_inv = MatrixX::Identity(3 * numverts, 3 * numverts);
+
+ /* Constraint matrix */
+ int numjoints_root = 3; /* root velocity constraint */
+ int numjoints_edges = numverts - 1; /* distance constraints */
+ int numjoints = numjoints_edges + numjoints_root;
+ MatrixX J = MatrixX::Zero(numjoints, 3 * numverts);
+ /* root velocity constraint */
+ J.block<3,3>(0, 0) = Matrix3f::Identity();
+ /* distance constraints */
+ BMVert *vertprev = NULL;
+ BM_ITER_STRANDS_ELEM_INDEX(vert, &iter, root, BM_VERTS_OF_STRAND, k) {
+ if (k > 0) {
+// float target_length = L[k];
+// if (target_length > 0.0f) {
+ Vector3 xa(vertprev->co);
+ Vector3 xb(vert->co);
+// Vector3f j = (xb - xa) / target_length;
+ Vector3f j = (xb - xa);
+ j.normalize();
+#ifdef DO_DEBUG
+ BKE_sim_debug_data_add_vector(vert->co, j.data(), 0,1,0, "hair solve", 3274, BLI_ghashutil_ptrhash(root), k);
+#endif
+
+ int i = k + 2; /* constraint index */
+ J.block<1,3>(i, 3*(k-1)) = -j.transpose();
+ J.block<1,3>(i, 3*(k) ) = j.transpose();
+// }
+ }
+
+ vertprev = vert;
+ }
+
+ /* A = J * M^-1 * J^T */
+ MatrixX A = J * M_inv * J.transpose();
+ /* b = -(J * M^-1 * F_ext + c) = -(J * v0 + c) */
+ VectorX c = VectorX::Zero(numjoints);
+ c.block<3,1>(0, 0) = root_v;
+ VectorX b = -(J * v + c);
+
+ /* Lagrange multipliers are the solution to A * lambda = b */
+ VectorX lambda = A.ldlt().solve(b);
+
+ /* calculate velocity correction by constraint forces */
+ VectorX dv = M_inv * J.transpose() * lambda;
+ v += dv;
+
+#ifdef DO_DEBUG
+ std::cout << "J = " << std::endl << J << std::endl;
+ std::cout << "A = " << std::endl << A << std::endl;
+ std::cout << "v = " << std::endl << v << std::endl;
+ std::cout << "b = " << std::endl << b << std::endl;
+ std::cout << "lambda = " << std::endl << lambda << std::endl;
+ BM_ITER_STRANDS_ELEM_INDEX(vert, &iter, root, BM_VERTS_OF_STRAND, k) {
+ BKE_sim_debug_data_add_vector(vert->co, &dv.coeff(3*k), 1,0,1, "hair solve", 3833, BLI_ghashutil_ptrhash(root), k);
+ BKE_sim_debug_data_add_vector(orig[k], &v.coeff(3*k), 0,1,1, "hair solve", 3811, BLI_ghashutil_ptrhash(root), k);
+ }
+#endif
+
+ BM_ITER_STRANDS_ELEM_INDEX(vert, &iter, root, BM_VERTS_OF_STRAND, k) {
+ add_v3_v3v3(vert->co, orig[k], &v.coeff(3*k));
+ }
+}
+
static void strands_solve_lagrange_multipliers(Object *ob, BMEditStrands *edit, float (*orig)[3])
{
+ using Eigen::Vector3f;
+ BMesh *bm = edit->base.bm;
+
+ BMIter iter;
+ BMVert *root;
+ BM_ITER_STRANDS(root, &iter, bm, BM_STRANDS_OF_MESH) {
+ int numverts = BM_strand_verts_count(root);
+ /* TODO if the root gets moved this would be non-zero */
+ Vector3f root_v = Vector3f(0.0f, 0.0f, 0.0f);
+
+ strand_solve(bm, root, orig, numverts, root_v);
+
+ orig += numverts;
+ }
}
+#endif
void BPH_strands_solve_constraints(Scene *scene, Object *ob, BMEditStrands *edit, float (*orig)[3])
{
@@ -451,12 +565,12 @@ void BPH_strands_solve_constraints(Scene *scene, Object *ob, BMEditStrands *edit
BKE_collision_cache_free(contacts);
}
- strands_apply_root_locations(edit);
-
#ifdef STRAND_CONSTRAINT_EDGERELAX
+ strands_apply_root_locations(edit);
strands_solve_edge_relaxation(edit);
#endif
#ifdef STRAND_CONSTRAINT_IK
+ strands_apply_root_locations(edit);
strands_solve_inverse_kinematics(ob, edit, orig);
#endif
#ifdef STRAND_CONSTRAINT_LAGRANGEMULT
More information about the Bf-blender-cvs
mailing list