Commit: 97b08dcdbb39b0cba2573d4b709b0aeb53c66e00
Author: mattoverby
Date: Wed Aug 5 14:27:40 2020 -0500
Branches: soc-2020-soft-body
https://developer.blender.org/rB97b08dcdbb39b0cba2573d4b709b0aeb53c66e00
update cg solver
===================================================================
M extern/softbody/src/admmpd_collision.cpp
M extern/softbody/src/admmpd_collision.h
M extern/softbody/src/admmpd_energy.cpp
M extern/softbody/src/admmpd_linsolve.cpp
M extern/softbody/src/admmpd_linsolve.h
M extern/softbody/src/admmpd_mesh.cpp
M extern/softbody/src/admmpd_mesh.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 release/scripts/startup/bl_ui/properties_physics_softbody.py
M source/blender/blenkernel/intern/softbody.c
M source/blender/makesdna/DNA_object_force_types.h
M source/blender/makesrna/intern/rna_object_force.c
===================================================================
diff --git a/extern/softbody/src/admmpd_collision.cpp
b/extern/softbody/src/admmpd_collision.cpp
index fee3456e022..1d37ffc688e 100644
--- a/extern/softbody/src/admmpd_collision.cpp
+++ b/extern/softbody/src/admmpd_collision.cpp
@@ -459,7 +459,7 @@ void EmbeddedMeshCollision::graph(
void EmbeddedMeshCollision::linearize(
const Eigen::MatrixXd *x,
std::vector<Eigen::Triplet<double> > *trips,
- std::vector<double> *d)
+ std::vector<double> *d) const
{
BLI_assert(x != NULL);
BLI_assert(x->cols() == 3);
@@ -475,7 +475,7 @@ void EmbeddedMeshCollision::linearize(
for (int i=0; i<np; ++i)
{
const Vector2i &pair_idx = vf_pairs[i];
- VFCollisionPair &pair =
per_vertex_pairs[pair_idx[0]][pair_idx[1]];
+ const VFCollisionPair &pair =
per_vertex_pairs[pair_idx[0]][pair_idx[1]];
int emb_p_idx = pair.p_idx;
// Vector3d p_pt =
meshdata.mesh->get_mapped_facet_vertex(x,emb_p_idx);
diff --git a/extern/softbody/src/admmpd_collision.h
b/extern/softbody/src/admmpd_collision.h
index d4a3b57f20c..4228a53f8f4 100644
--- a/extern/softbody/src/admmpd_collision.h
+++ b/extern/softbody/src/admmpd_collision.h
@@ -25,8 +25,6 @@ struct VFCollisionPair {
class Collision {
public:
- typedef Discregrid::CubicLagrangeDiscreteGrid SDFType;
-
struct ObstacleData {
bool has_obs() const { return F.rows()>0; }
Eigen::MatrixXd V;
@@ -65,11 +63,12 @@ public:
const unsigned int *faces,
int nf);
- // Linearize the constraints about x and return Jacobian.
+ // Linearizes active collision pairs about x
+ // for the constraint Cx=d
virtual void linearize(
const Eigen::MatrixXd *x,
std::vector<Eigen::Triplet<double> > *trips,
- std::vector<double> *d) = 0;
+ std::vector<double> *d) const = 0;
// Given a point and a mesh, perform
// discrete collision detection.
@@ -79,6 +78,7 @@ public:
const Eigen::Vector3d &pt_t1,
const ObstacleData *obs) const;
+ // Perform self collision detection
virtual std::pair<bool,VFCollisionPair>
detect_against_self(
int pt_idx,
@@ -101,13 +101,13 @@ public:
void graph(
std::vector<std::set<int> > &g);
-
+
// Linearizes the collision pairs about x
- // for the constraint Kx=l
+ // for the constraint Cx=d
void linearize(
const Eigen::MatrixXd *x,
std::vector<Eigen::Triplet<double> > *trips,
- std::vector<double> *d);
+ std::vector<double> *d) const;
// Updates the tetmesh BVH for self collisions.
void update_bvh(
diff --git a/extern/softbody/src/admmpd_energy.cpp
b/extern/softbody/src/admmpd_energy.cpp
index 394769d6291..eb4cb9c9534 100644
--- a/extern/softbody/src/admmpd_energy.cpp
+++ b/extern/softbody/src/admmpd_energy.cpp
@@ -147,6 +147,18 @@ void EnergyTerm::update_tri(
double w2 = weight*weight;
zi_T = (kv*p + w2*zi_T) / (w2 + kv);
+ // Apply strain limiting
+ bool check_strain = lame.m_limit[0] > 0.0 || lame.m_limit[1] < 99.0;
+ if (check_strain)
+ {
+ double l_col0 = zi_T.col(0).norm();
+ double l_col1 = zi_T.col(1).norm();
+ if (l_col0 < lame.m_limit[0]) { zi_T.col(0) *= (
lame.m_limit[0]/l_col0 ); }
+ if (l_col1 < lame.m_limit[0]) { zi_T.col(1) *= (
lame.m_limit[0]/l_col1 ); }
+ if (l_col0 > lame.m_limit[1]) { zi_T.col(0) *= (
lame.m_limit[1]/l_col0 ); }
+ if (l_col1 > lame.m_limit[1]) { zi_T.col(1) *= (
lame.m_limit[1]/l_col1 ); }
+ }
+
ui += (Dix - zi_T.transpose());
u->block<2,3>(index,0) = ui;
z->block<2,3>(index,0) = zi_T.transpose();
diff --git a/extern/softbody/src/admmpd_linsolve.cpp
b/extern/softbody/src/admmpd_linsolve.cpp
index e926c36db15..5330e07cfaf 100644
--- a/extern/softbody/src/admmpd_linsolve.cpp
+++ b/extern/softbody/src/admmpd_linsolve.cpp
@@ -8,109 +8,327 @@
#include <set>
#include <unordered_set>
#include <random>
+#include <thread>
#include "BLI_assert.h"
#include "BLI_task.h"
namespace admmpd {
using namespace Eigen;
-void ConjugateGradients::init_solve(
+// Throws an exception for a given function with message
+static inline void throw_err(const std::string &f, const std::string &m)
+{
+ throw std::runtime_error("LinearSolver::"+f+": "+m);
+}
+
+LDLT::LDLT() :
+ last_pk(-1)
+ {}
+
+void LDLT::init_solve(
+ const Mesh *mesh,
+ const Options *options,
+ const Collision *collision,
+ SolverData *data)
+{
+ (void)(collision);
+ int nx = std::max((int)data->x.rows(),1);
+
+ // Get the P matrix
+ std::set<int> pin_inds;
+ std::vector<Triplet<double> > trips;
+ std::vector<double> q_coeffs;
+ bool replicate = false;
+ bool new_P = mesh->linearize_pins(trips, q_coeffs, pin_inds, replicate);
+
+ // If we've changed the stiffness but not the pins,
+ // the P matrix is still changing
+ double pk = options->mult_pk * data->A_diag_max;
+ if (std::abs(pk-last_pk) > 1e-8 && trips.size()>0)
+ new_P = true;
+
+ // Compute P
+ last_pk = pk;
+ int np = q_coeffs.size();
+ SparseMatrix<double> P;
+ MatrixXd q;
+ if (m_Ptq.rows() != nx) { m_Ptq.resize(nx,3); }
+ if (np==0)
+ { // no springs
+ P.resize(1,nx);
+ P.setZero();
+ q = MatrixXd::Zero(1,3);
+ m_Ptq.setZero();
+ }
+ else
+ {
+ P.resize(np,nx);
+ P.setFromTriplets(trips.begin(), trips.end());
+ q.resize(np,3);
+ for (int i=0; i<np; ++i)
+ {
+ q(i,0) = q_coeffs[i*3+0];
+ q(i,1) = q_coeffs[i*3+1];
+ q(i,2) = q_coeffs[i*3+2];
+ }
+ m_Ptq = pk * P.transpose() * q;
+ }
+
+ // Compute A + P'P and factorize:
+ // 1) A not computed
+ // 2) P has changed
+ // 3) factorization not set
+ if (m_A_PtP.nonZeros()==0 || new_P || !m_ldlt_A_PtP)
+ {
+ if (!m_ldlt_A_PtP)
+ m_ldlt_A_PtP = std::make_unique<Cholesky>();
+
+ m_A_PtP = SparseMatrix<double>(data->A) + pk * P.transpose()*P;
+ m_ldlt_A_PtP->compute(m_A_PtP);
+ if(m_ldlt_A_PtP->info() != Eigen::Success) {
+ throw_err("init_solve","facorization failed");
+ }
+ geom::make_n3<double>(m_A_PtP,m_A_PtP_3); // replicate
+ }
+} // end init solve
+
+void LDLT::solve(
+ const Mesh *mesh,
const Options *options,
+ const Collision *collision,
SolverData *data)
{
- int nx3 = 3*data->x.rows();
- if (nx3==0)
+ if (!m_ldlt_A_PtP)
+ throw_err("solve","you must call init_solve before solve");
+
+ int nx = std::max((int)data->x.rows(),1);
+
+ // Linearize collision constraints
+ RowSparseMatrix<double> C;
+ VectorXd d;
+ if (collision != nullptr)
+ {
+ std::vector<double> d_coeffs;
+ std::vector<Eigen::Triplet<double> > trips;
+ collision->linearize(&data->x, &trips, &d_coeffs);
+
+ int nc = d_coeffs.size();
+ if (nc>0)
+ {
+ d = Map<VectorXd>(d_coeffs.data(), d_coeffs.size());
+ C.resize(nc,nx*3);
+ C.setFromTriplets(trips.begin(),trips.end());
+ }
+ }
+
+ // Compute RHS
+ m_rhs.noalias() = data->M_xbar + data->DtW2*(data->z-data->u) + m_Ptq;
+
+ // If there are no collision constraints,
+ // we can use our initial factorization.
+ if (C.nonZeros()==0)
+ {
+ data->x.noalias() = m_ldlt_A_PtP->solve(m_rhs);
return;
+ }
- A3_PtP_CtC.resize(nx3,nx3);
- b.resize(data->x.rows(),data->x.cols());
- CtC.resize(nx3,nx3);
- Ctd.resize(nx3);
- Ptq.resize(nx3);
- b3_Ptq_Ctd.resize(nx3);
- r.resize(nx3);
- z.resize(nx3);
- p.resize(nx3);
- Ap.resize(nx3);
- double pk = options->mult_pk * data->A_diag_max;
- Ptq = pk * data->P.transpose()*data->q;
+ // Otherwise we have to solve the full system:
+ // (A + PtP + CtC) x = b + Ptq + Ctd
+ double ck = options->mult_ck * data->A_diag_max;
+ SparseMatrix<double> A3 = m_A_PtP_3 + ck * C.transpose()*C;
+ VectorXd Ctd = ck * C.transpose() * d;
+ VectorXd rhs3(nx*3);
+ for (int i=0; i<nx; ++i)
+ {
+ rhs3.segment<3>(i*3) =
+ m_rhs.row(i).transpose() +
+ Ctd.segment<3>(i*3);
+ }
+
+ Cholesky ldlt_full(A3);
+ VectorXd x3 = ldlt_full.solve(rhs3);
+ for (int i=0; i<nx; ++i)
+ data->x.row(i) = x3.segment<3>(i*3);
+
+} // end solve
+
+ConjugateGradients::ConjugateGradients() :
+ factor_A_PtP(true),
+ last_pk(-1)
+ {}
+
+void ConjugateGradients::init_solve(
+ const Mesh *mesh,
+ const Options *options,
+ const Collision *collision,
+ SolverData *data)
+{
+ int nx = data->x.rows();
+ if (nx==0) { return; }
+
+ if (!factor_A_PtP)
+ throw_err("init_solve","todo: not factor A+P'P");
+
+ if (!m_ldlt)
+ m_ldlt = std::make_unique<LDLT>();
+
+ // We'll just use our LDLT implementation
+ // to decide when to refactor the matrix
+ m_ldlt->init_solve(mesh,options,collision,data);
+
+ rhs.resize(nx,3);
+ Ctd.resize(nx,3);
+ r.resize(nx,3);
+ z.resize(nx,3);
+ p.resize(nx,3);
+ p3.resize(nx*3);
+ Ap.resize(nx,3);
}
void ConjugateGradients::solve(
- const Options *options,
- SolverData *data,
- Collision *collision)
+ const Mesh *mesh,
+ const Options *options,
+ const Collision *collision,
+ SolverData *data)
{
- (void)(collision); // unused
+ auto map_vector_to_matrix = [](const VectorXd &x3, MatrixXd &x)
+ {
+ int nx = x3.rows()/3;
+ if (x.rows() != nx) { x.resize(nx,3); }
+ for (int i=0; i<nx; ++i)
+ x.row(i) = x3.segment<3>(i*3);
+ };
+
+ auto map_matrix_to_vector = [](const MatrixXd &x, VectorXd &x3)
+ {
+ int nx = x.rows();
+ if (x3.rows() != nx*3) { x3.resize(nx*3); }
+ for (int i=0; i<nx; ++i)
+ x3.segment<3>(i*3) = x.row(i);
+ };
+
+ auto mat_inner = [](const MatrixXd &A, const MatrixXd &B)
+ {
+ double sum = 0.0;
+ int rows = std::min(A.rows(), B.rows());
+ int dim = std::min(A.cols(), B.cols());
+ for (int i=0; i<rows; ++i)
+ {
+ for (int j=0; j<dim; ++j)
+ {
+ sum += A(i,j)*B(i,j);
+ }
+ }
+ return sum;
+ };
+
BLI_assert(data != NULL);
BLI_assert(options != NULL);
int nx = data->x.rows();
BLI_assert(nx > 0);
- BLI_assert(data->C.cols() == nx*3);
- BLI_assert(data->d.rows() > 0);
- BLI_assert(data->C.rows() == data->d.rows());
- // Update A matrix with collision constraints
- double ck = options->mult_ck * data->A_diag_max;
- RowSparseMatrix<double> *A;
- if (data->C.nonZeros()>0)
- {
- Ctd = ck * data->C.transpose()*data->d;
- CtC = ck * data->C.transpose()*data->C;
- A3_PtP_CtC = data->A3_plus_PtP + CtC;
- A = &A3_PtP_CtC;
- }
- else
+ if (!m_ldlt)
+ throw_err("solve","you must call init_solve before solve");
+
+ // Linearize collision constraints
+ RowSparseMatrix<double> C;
+ VectorXd d;
+ if (collision != nullptr)
{
- Ctd.setZero();
- A = &data->A3_plus_PtP;
+ std::vector<double> d_coeffs;
@@ Diff output truncated at 10240 characters. @@
_______________________________________________
Bf-blender-cvs mailing list
[email protected]
https://lists.blender.org/mailman/listinfo/bf-blender-cvs
