[Bf-blender-cvs] [14a76718e5e] soc-2020-soft-body: added goal positions
over0219
noreply at git.blender.org
Thu Jul 9 21:42:24 CEST 2020
Commit: 14a76718e5e0d16c87e45f27f20a31aebe1e19e6
Author: over0219
Date: Thu Jul 9 14:42:20 2020 -0500
Branches: soc-2020-soft-body
https://developer.blender.org/rB14a76718e5e0d16c87e45f27f20a31aebe1e19e6
added goal positions
===================================================================
M extern/softbody/CMakeLists.txt
M extern/softbody/src/admmpd_bvh.h
M extern/softbody/src/admmpd_linsolve.cpp
A extern/softbody/src/admmpd_pin.cpp
A extern/softbody/src/admmpd_pin.h
M extern/softbody/src/admmpd_solver.cpp
M extern/softbody/src/admmpd_solver.h
M extern/softbody/src/admmpd_types.h
M intern/softbody/admmpd_api.cpp
M intern/softbody/admmpd_api.h
M source/blender/blenkernel/intern/softbody.c
===================================================================
diff --git a/extern/softbody/CMakeLists.txt b/extern/softbody/CMakeLists.txt
index cadd9ef51ee..7f631400ce0 100644
--- a/extern/softbody/CMakeLists.txt
+++ b/extern/softbody/CMakeLists.txt
@@ -43,6 +43,8 @@ set(SRC
src/admmpd_linsolve.cpp
src/admmpd_math.h
src/admmpd_math.cpp
+ src/admmpd_pin.h
+ src/admmpd_pin.cpp
src/admmpd_solver.h
src/admmpd_solver.cpp
src/admmpd_tetmesh.h
diff --git a/extern/softbody/src/admmpd_bvh.h b/extern/softbody/src/admmpd_bvh.h
index dca84e6803f..19f47545435 100644
--- a/extern/softbody/src/admmpd_bvh.h
+++ b/extern/softbody/src/admmpd_bvh.h
@@ -49,6 +49,8 @@ protected:
AABB aabb;
Node *left, *right;
std::vector<int> prims;
+ VecType normal;
+ T angle;
bool is_leaf() const { return prims.size()>0; }
Node() : left(nullptr), right(nullptr) {}
~Node()
diff --git a/extern/softbody/src/admmpd_linsolve.cpp b/extern/softbody/src/admmpd_linsolve.cpp
index 1ad3fe9a78d..7d5a29eb964 100644
--- a/extern/softbody/src/admmpd_linsolve.cpp
+++ b/extern/softbody/src/admmpd_linsolve.cpp
@@ -54,10 +54,12 @@ void ConjugateGradients::solve(
BLI_assert(data->C.cols() == nx*3);
BLI_assert(data->d.rows() > 0);
BLI_assert(data->C.rows() == data->d.rows());
+ BLI_assert(data->PtP.cols() == nx*3);
+ BLI_assert(data->PtP.rows() == nx*3);
// If we don't have any constraints, we don't need to perform CG
data->b.noalias() = data->M_xbar + data->DtW2*(data->z-data->u);
- if (data->C.nonZeros()==0)
+ if (data->C.nonZeros()==0 && data->PtP.nonZeros()==0)
{
data->x = data->ldltA.solve(data->b);
return;
@@ -71,28 +73,30 @@ void ConjugateGradients::solve(
gsdata->z.resize(nx*3);
gsdata->p.resize(nx*3);
gsdata->A3p.resize(nx*3);
- gsdata->b3_plus_Ctd.resize(nx*3);
+ gsdata->b3_Ctd_Ptx.resize(nx*3);
make_n3(data->A, gsdata->A3);
}
- gsdata->Ctd = data->spring_k * data->C.transpose()*data->d;
gsdata->CtC = data->spring_k * data->C.transpose()*data->C;
- gsdata->A3_plus_CtC = gsdata->A3 + gsdata->CtC;
+ gsdata->Ctd = data->spring_k * data->C.transpose()*data->d;
+ gsdata->A3_CtC_PtP = gsdata->A3 + gsdata->CtC + data->PtP;
VectorXd x3(nx*3);
for (int i=0; i<nx; ++i)
{
- Vector3d bi = data->b.row(i);
- gsdata->b3_plus_Ctd.segment<3>(i*3) = bi+gsdata->Ctd.segment<3>(i*3);
+ gsdata->b3_Ctd_Ptx.segment<3>(i*3) =
+ data->b.row(i).transpose() +
+ gsdata->Ctd.segment<3>(i*3) +
+ data->Ptq.segment<3>(i*3);
x3.segment<3>(i*3) = data->x.row(i);
}
- gsdata->r.noalias() = gsdata->b3_plus_Ctd - gsdata->A3_plus_CtC*x3;
+ gsdata->r.noalias() = gsdata->b3_Ctd_Ptx - gsdata->A3_CtC_PtP*x3;
solve_Ax_b(data,&gsdata->z,&gsdata->r);
gsdata->p = gsdata->z;
for (int iter=0; iter<options->max_cg_iters; ++iter)
{
- gsdata->A3p.noalias() = gsdata->A3_plus_CtC*gsdata->p;
+ gsdata->A3p.noalias() = gsdata->A3_CtC_PtP*gsdata->p;
double p_dot_Ap = gsdata->p.dot(gsdata->A3p);
if (p_dot_Ap==0.0)
break;
@@ -196,7 +200,7 @@ void GaussSeidel::solve(
Vector3d inv_aii(0,0,0);
for (int j=0; j<3; ++j)
{
- InnerIter rit(td->data->gsdata.A3_plus_CtC, idx*3+j);
+ InnerIter rit(td->data->gsdata.A3_CtC_PtP, idx*3+j);
for (; rit; ++rit)
{
double v = rit.value();
@@ -218,7 +222,7 @@ void GaussSeidel::solve(
} // end loop segment
// Update x
- Vector3d bi = td->data->gsdata.b3_plus_Ctd.segment<3>(idx*3);
+ Vector3d bi = td->data->gsdata.b3_Ctd_Ptx.segment<3>(idx*3);
//Vector3d bi = td->data->b.row(idx);
Vector3d xi = td->data->x.row(idx);
@@ -277,6 +281,15 @@ void GaussSeidel::init_solve(
SolverData *data,
Collision *collision)
{
+
+ // TODO:
+ //
+ // When it comes to improving run time of Gauss-Seidel after
+ // the instability issues have been addressed, reducing the
+ // matrix-vector mults that occur here should be a priority.
+ // Many of them are unnecessary and can be done
+ // within the Gauss-Seidel sweeps!
+
BLI_assert(options != nullptr);
BLI_assert(data != nullptr);
BLI_assert(collision != nullptr);
@@ -308,13 +321,13 @@ void GaussSeidel::init_solve(
{
data->gsdata.CtC = data->spring_k * data->C.transpose()*data->C;
data->gsdata.Ctd.noalias() = data->spring_k * data->C.transpose()*data->d;
- data->gsdata.A3_plus_CtC = data->gsdata.A3 + data->gsdata.CtC;
- data->gsdata.b3_plus_Ctd.resize(nx*3);
+ data->gsdata.A3_CtC_PtP = data->gsdata.A3 + data->gsdata.CtC;
+ data->gsdata.b3_Ctd_Ptx.resize(nx*3);
for (int i=0; i<nx; ++i)
{
- data->gsdata.b3_plus_Ctd[i*3+0] = data->b(i,0)+data->gsdata.Ctd[i*3+0];
- data->gsdata.b3_plus_Ctd[i*3+1] = data->b(i,1)+data->gsdata.Ctd[i*3+1];
- data->gsdata.b3_plus_Ctd[i*3+2] = data->b(i,2)+data->gsdata.Ctd[i*3+2];
+ data->gsdata.b3_Ctd_Ptx[i*3+0] = data->b(i,0)+data->gsdata.Ctd[i*3+0];
+ data->gsdata.b3_Ctd_Ptx[i*3+1] = data->b(i,1)+data->gsdata.Ctd[i*3+1];
+ data->gsdata.b3_Ctd_Ptx[i*3+2] = data->b(i,2)+data->gsdata.Ctd[i*3+2];
}
std::vector<std::set<int> > c_graph;
collision->graph(c_graph);
@@ -328,13 +341,13 @@ void GaussSeidel::init_solve(
if (data->gsdata.Ctd.rows() != nx*3)
data->gsdata.Ctd.resize(nx*3);
data->gsdata.Ctd.setZero();
- data->gsdata.A3_plus_CtC = data->gsdata.A3;
- data->gsdata.b3_plus_Ctd.resize(nx*3);
+ data->gsdata.A3_CtC_PtP = data->gsdata.A3;
+ data->gsdata.b3_Ctd_Ptx.resize(nx*3);
for (int i=0; i<nx; ++i)
{
- data->gsdata.b3_plus_Ctd[i*3+0] = data->b(i,0);
- data->gsdata.b3_plus_Ctd[i*3+1] = data->b(i,1);
- data->gsdata.b3_plus_Ctd[i*3+2] = data->b(i,2);
+ data->gsdata.b3_Ctd_Ptx[i*3+0] = data->b(i,0);
+ data->gsdata.b3_Ctd_Ptx[i*3+1] = data->b(i,1);
+ data->gsdata.b3_Ctd_Ptx[i*3+2] = data->b(i,2);
}
data->gsdata.A3_plus_CtC_colors = data->gsdata.A_colors;
}
diff --git a/extern/softbody/src/admmpd_pin.cpp b/extern/softbody/src/admmpd_pin.cpp
new file mode 100644
index 00000000000..ff983a08e87
--- /dev/null
+++ b/extern/softbody/src/admmpd_pin.cpp
@@ -0,0 +1,74 @@
+// Copyright Matt Overby 2020.
+// Distributed under the MIT License.
+
+#include "admmpd_pin.h"
+#include "BLI_assert.h"
+
+namespace admmpd {
+using namespace Eigen;
+
+void EmbeddedMeshPin::clear()
+{
+ pin_k.clear();
+ pin_pos.clear();
+}
+
+void EmbeddedMeshPin::set_pin(
+ int idx,
+ const Eigen::Vector3d &qi,
+ const Eigen::Vector3d &ki_)
+{
+ if (!mesh)
+ return;
+
+ if (idx<0 || idx>=mesh->emb_rest_x.rows())
+ return;
+
+ // Clamp
+ Vector3d ki = ki_;
+ for (int i=0; i<3; ++i)
+ ki[i] = std::max(0.0, ki[i]);
+
+ pin_k[idx] = ki;
+ pin_pos[idx] = qi;
+}
+
+void EmbeddedMeshPin::linearize(
+ const Eigen::MatrixXd *x, // not used yet
+ std::vector<Eigen::Triplet<double> > *trips,
+ std::vector<double> *q)
+{
+
+ (void)(x);
+ int np = pin_k.size();
+ trips->reserve((int)trips->size() + np*3*4);
+ q->reserve((int)q->size() + np*3);
+
+ std::unordered_map<int,Eigen::Vector3d>::const_iterator it_k = pin_k.begin();
+ for (; it_k != pin_k.end(); ++it_k)
+ {
+ int emb_idx = it_k->first;
+ const Vector3d &qi = pin_pos[emb_idx];
+ const Vector3d &ki = it_k->second;
+
+ int tet_idx = mesh->emb_vtx_to_tet[emb_idx];
+ RowVector4d bary = mesh->emb_barys.row(emb_idx);
+ RowVector4i tet = mesh->tets.row(tet_idx);
+
+ for (int i=0; i<3; ++i)
+ {
+ int p_idx = q->size();
+ q->emplace_back(qi[i]*ki[i]);
+ for (int j=0; j<4; ++j)
+ trips->emplace_back(p_idx, tet[j]*3+i, bary[j]*ki[i]);
+ }
+ }
+
+} // end linearize
+
+//bool EmbeddedMeshPin::has_pin(int idx) const
+//{
+// return pin_k.count(idx);
+//}
+
+} // namespace admmpd
diff --git a/extern/softbody/src/admmpd_pin.h b/extern/softbody/src/admmpd_pin.h
new file mode 100644
index 00000000000..8080d9d2c85
--- /dev/null
+++ b/extern/softbody/src/admmpd_pin.h
@@ -0,0 +1,88 @@
+// Copyright Matt Overby 2020.
+// Distributed under the MIT License.
+
+#ifndef ADMMPD_PIN_H_
+#define ADMMPD_PIN_H_
+
+#include "admmpd_bvh.h"
+#include "admmpd_types.h"
+#include <unordered_map>
+#include <vector>
+
+namespace admmpd {
+
+class Pin {
+public:
+
+ virtual ~Pin() {}
+
+ // Clears all pins
+ virtual void clear() = 0;
+
+ // Set the pin location (q) and per-axis stiffness (k)
+ // Stiffness should be 0 to 1. It can go larger, but
+ // the resulting matrix will be poorly conditioned.
+ virtual void set_pin(
+ int idx,
+ const Eigen::Vector3d &q,
+ const Eigen::Vector3d &k) = 0;
+
+ // Returns true if the vert is pinned
+// virtual bool has_pin(int idx) const = 0;
+
+ // Creates linearization for constraint:
+ // Px=q with stiffnesses baked in
+ virtual void linearize(
+ const Eigen::MatrixXd *x, // not used yet
+ std::vector<Eigen::Triplet<double> > *trips,
+ std::vector<double> *q) = 0;
+
+ // Returns per-axis pin stiffness
+// virtual Eigen::Vector3d get_pin_k(int idx) const = 0;
+
+ // Returns pin location, or zero vector if not set
+ // virtual Eigen::Vector3d get_pin_pos(int idx) const = 0;
+};
+
+class EmbeddedMeshPin : public Pin {
+public:
+ EmbeddedMeshPin(const EmbeddedMeshData *mesh_) : mesh(mesh_){}
+
+ // Clears all pins
+ void clear();
+
+ // Set the pin location of the embedded vertex
+ void set_pin(
+ int idx,
+ const Eigen::Vector3d &p,
+ const Eigen::Vector3d &k);
+
+ // Returns true if the deforming vertex is pinned
+// bool has_pin(int idx) const;
+
+ void linearize(
+ const Eigen::MatrixXd *x, // not used yet
+ std::vector<Eigen::Triplet<double> > *trips,
+ std::vector<double> *q);
+
+ // Returns per-axis pin stiffness of the deforming vertex (not embedded)
+ // or zero if not pinned
+ // Baryweights are included.
+// Eigen::Vector3d get_pin_k(int idx) const;
+
+ // Returns pin location of the deforming vertex (not embedded)
+ // or zero vector if not set
+// Eigen::Vector3d get_pin_pos(int idx) const;
+
+protected:
+ // A ptr to the embedded mesh data
+ const EmbeddedMeshData *mesh;
+
+ // Pins for embedded vertices:
+ std::unordered_map<int,Eigen::Vector3d> pin_k;
+ std::unordered_map<int,Eigen::Vector3d> pin_pos;
+};
+
+} // namespace admmpd
+
+#endif // ADMMPD_PIN_H_
diff --git a/extern/softbody/src/admmpd_solver.cpp b/extern/softbody/src/admmpd_solver.cpp
index eeb3163ccd9..e1408b4d53e 1
@@ Diff output truncated at 10240 characters. @@
More information about the Bf-blender-cvs
mailing list