[Bf-blender-cvs] [2c2243c012d] soc-2020-soft-body: cleanup
mattoverby
noreply at git.blender.org
Wed Jul 29 02:14:10 CEST 2020
Commit: 2c2243c012d1844e459631593839ccd70d3f51be
Author: mattoverby
Date: Tue Jul 28 19:14:06 2020 -0500
Branches: soc-2020-soft-body
https://developer.blender.org/rB2c2243c012d1844e459631593839ccd70d3f51be
cleanup
===================================================================
M extern/softbody/src/admmpd_linsolve.cpp
M extern/softbody/src/admmpd_linsolve.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/makesrna/intern/rna_object_force.c
===================================================================
diff --git a/extern/softbody/src/admmpd_linsolve.cpp b/extern/softbody/src/admmpd_linsolve.cpp
index 95751957a39..e926c36db15 100644
--- a/extern/softbody/src/admmpd_linsolve.cpp
+++ b/extern/softbody/src/admmpd_linsolve.cpp
@@ -23,6 +23,7 @@ void ConjugateGradients::init_solve(
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);
@@ -45,7 +46,6 @@ void ConjugateGradients::solve(
BLI_assert(options != NULL);
int nx = data->x.rows();
BLI_assert(nx > 0);
- BLI_assert(data->b.rows() == nx);
BLI_assert(data->C.cols() == nx*3);
BLI_assert(data->d.rows() > 0);
BLI_assert(data->C.rows() == data->d.rows());
@@ -68,11 +68,11 @@ void ConjugateGradients::solve(
// Compute RHS
VectorXd x3(nx*3);
- data->b.noalias() = data->M_xbar + data->DtW2*(data->z-data->u);
+ b.noalias() = data->M_xbar + data->DtW2*(data->z-data->u);
for (int i=0; i<nx; ++i)
{
b3_Ptq_Ctd.segment<3>(i*3) =
- data->b.row(i).transpose() +
+ b.row(i).transpose() +
Ptq.segment<3>(i*3) +
Ctd.segment<3>(i*3);
x3.segment<3>(i*3) = data->x.row(i);
diff --git a/extern/softbody/src/admmpd_linsolve.h b/extern/softbody/src/admmpd_linsolve.h
index eaed804afa2..1df2d4b35a3 100644
--- a/extern/softbody/src/admmpd_linsolve.h
+++ b/extern/softbody/src/admmpd_linsolve.h
@@ -39,6 +39,7 @@ public:
protected:
RowSparseMatrix<double> A3_PtP_CtC;
RowSparseMatrix<double> CtC;
+ Eigen::MatrixXd b; // M xbar + DtW2(z-u)
Eigen::VectorXd Ctd; // ck * Ct d
Eigen::VectorXd Ptq; // pk * Pt q
Eigen::VectorXd b3_Ptq_Ctd; // M xbar + DtW2(z-u) + ck Ct d + pk Pt q
diff --git a/extern/softbody/src/admmpd_solver.cpp b/extern/softbody/src/admmpd_solver.cpp
index 6302384f388..2f8214821d2 100644
--- a/extern/softbody/src/admmpd_solver.cpp
+++ b/extern/softbody/src/admmpd_solver.cpp
@@ -39,11 +39,12 @@ bool Solver::init(
data->x = mesh->rest_prim_verts();
data->v.resize(data->x.rows(), 3);
data->v.setZero();
- data->tets = mesh->prims();
mesh->compute_masses(data->x, options->density_kgm3, data->m);
- init_matrices(options,data);
+ init_matrices(mesh,options,data);
- printf("Solver::init:\n\tNum tets: %d\n\tNum verts: %d\n",(int)data->tets.rows(),(int)data->x.rows());
+ int nt = mesh->prims().rows();
+ int nx = data->x.rows();
+ printf("Solver::init:\n\tNum tets: %d\n\tNum verts: %d\n",nt,nx);
return true;
} // end init
@@ -60,8 +61,8 @@ int Solver::solve(
BLI_assert(data->x.rows() > 0);
BLI_assert(options->max_admm_iters > 0);
- update_pin_matrix(mesh,options,data);
- update_global_matrix(options,data);
+ update_pin_matrix(mesh,options,data); // P,q
+ update_global_matrix(options,data); // A+PtP
ConjugateGradients cg;
cg.init_solve(options,data);
@@ -83,8 +84,16 @@ int Solver::solve(
update_collisions(options,data,collision);
// Solve Ax=b s.t. Px=q and Cx=d
+ data->x_prev = data->x;
cg.solve(options,data,collision);
+ // Check convergence
+ if (options->min_res>0)
+ {
+ if (residual_norm(options,data) < options->min_res)
+ break;
+ }
+
} // end solver iters
// Update velocity (if not static solve)
@@ -94,6 +103,16 @@ int Solver::solve(
return iters;
} // end solve
+double Solver::residual_norm(
+ const Options *options,
+ SolverData *data)
+{
+ (void)(options);
+ double ra = ((data->D*data->x) - data->z).norm();
+ double rx = (data->D*(data->x-data->x_prev)).norm();
+ return ra + rx;
+}
+
void Solver::init_solve(
const Mesh *mesh,
const Options *options,
@@ -220,6 +239,7 @@ void Solver::update_collisions(
} // end update constraints
void Solver::init_matrices(
+ const Mesh *mesh,
const Options *options,
SolverData *data)
{
@@ -243,15 +263,15 @@ void Solver::init_matrices(
// Add per-element energies to data
std::vector<Triplet<double> > trips;
- append_energies(options,data,trips);
+ append_energies(mesh,options,data,trips);
if (trips.size()==0)
{
throw_err("compute_matrices","No reduction coeffs");
}
int n_row_D = trips.back().row()+1;
- RowSparseMatrix<double> W2;
- compute_weight_matrix_squared(options,data,n_row_D,&W2);
+ update_weight_matrix(options,data,n_row_D);
+ RowSparseMatrix<double> W2 = data->W*data->W;
// Constraint data
data->C.resize(1,nx*3);
@@ -263,8 +283,6 @@ void Solver::init_matrices(
data->D.resize(n_row_D,nx);
data->D.setFromTriplets(trips.begin(), trips.end());
data->DtW2 = data->D.transpose() * W2;
- data->b.resize(nx,3);
- data->b.setZero();
update_global_matrix(options,data);
// Perform factorization
@@ -278,30 +296,29 @@ void Solver::init_matrices(
} // end compute matrices
-void Solver::compute_weight_matrix_squared(
+void Solver::update_weight_matrix(
const Options *options,
SolverData *data,
- int rows,
- RowSparseMatrix<double> *W2) const
+ int rows)
{
(void)(options);
- W2->resize(rows,rows);
+ data->W.resize(rows,rows);
VectorXi W_nnz = VectorXi::Ones(rows);
- W2->reserve(W_nnz);
+ data->W.reserve(W_nnz);
int ne = data->indices.size();
if (ne != (int)data->weights.size())
- throw_err("compute_weight_matrix","bad num indices/weights");
+ throw_err("update_weight_matrix","bad num indices/weights");
for (int i=0; i<ne; ++i)
{
const Vector3i &idx = data->indices[i];
if (idx[0]+idx[1] > rows)
- throw_err("compute_weight_matrix","bad matrix dim");
+ throw_err("update_weight_matrix","bad matrix dim");
for (int j=0; j<idx[1]; ++j)
- W2->coeffRef(idx[0]+j,idx[0]+j) = data->weights[i]*data->weights[i];
+ data->W.coeffRef(idx[0]+j,idx[0]+j) = data->weights[i];
}
- W2->finalize();
+ data->W.finalize();
}
void Solver::update_pin_matrix(
@@ -354,30 +371,32 @@ void Solver::update_global_matrix(
if (dt2 < 0) // static solve
dt2 = 1.0;
- SparseMatrix<double> A = data->DtW2 * data->D;
+ data->A = data->DtW2 * data->D;
data->A_diag_max = 0;
for (int i=0; i<nx; ++i)
{
- A.coeffRef(i,i) += data->m[i]/dt2;
- double Aii = A.coeff(i,i);
+ data->A.coeffRef(i,i) += data->m[i]/dt2;
+ double Aii = data->A.coeff(i,i);
if (Aii>data->A_diag_max)
data->A_diag_max = Aii;
}
SparseMatrix<double> A3;
- geom::make_n3<double>(A,A3);
+ geom::make_n3<double>(data->A,A3);
double pk = options->mult_pk * data->A_diag_max;
data->A3_plus_PtP = A3 + pk * data->P.transpose()*data->P;
}
void Solver::append_energies(
+ const Mesh *mesh,
const Options *options,
SolverData *data,
std::vector<Triplet<double> > &D_triplets)
{
BLI_assert(data != NULL);
BLI_assert(options != NULL);
- int nt = data->tets.rows();
+ Ref<const MatrixXi> tets = mesh->prims();
+ int nt = tets.rows();
BLI_assert(nt > 0);
int nx = data->x.rows();
@@ -398,7 +417,7 @@ void Solver::append_energies(
int energy_index = 0;
for (int i=0; i<nt; ++i)
{
- RowVector4i ele = data->tets.row(i);
+ RowVector4i ele = tets.row(i);
// Initialize the energy
data->rest_volumes.emplace_back();
diff --git a/extern/softbody/src/admmpd_solver.h b/extern/softbody/src/admmpd_solver.h
index 3514d768d50..d3b1a99c329 100644
--- a/extern/softbody/src/admmpd_solver.h
+++ b/extern/softbody/src/admmpd_solver.h
@@ -32,31 +32,45 @@ public:
protected:
- void update_collisions(
+ // Returns the combined residual norm
+ double residual_norm(
const Options *options,
- SolverData *data,
- Collision *collision);
+ SolverData *data);
+ // Computes start-of-solve quantites
void init_solve(
const Mesh *mesh,
const Options *options,
SolverData *data,
Collision *collision);
+ // Performs collision detection
+ // and updates C, d
+ void update_collisions(
+ const Options *options,
+ SolverData *data,
+ Collision *collision);
+
+ // Update z and u in parallel
void solve_local_step(
const Options *options,
SolverData *data);
+ // Called once at start of simulation.
+ // Computes constant quantities
void init_matrices(
+ const Mesh *mesh,
const Options *options,
SolverData *data);
- void compute_weight_matrix_squared(
+ // Computes the W matrix
+ // from the current weights
+ void update_weight_matrix(
const Options *options,
SolverData *data,
- int rows,
- RowSparseMatrix<double> *W2) const;
+ int rows);
+ // Linearizes pin constraints P, q
void update_pin_matrix(
const Mesh *mesh,
const Options *options,
@@ -67,7 +81,9 @@ protected:
const Options *options,
SolverData *data);
+ // Generates energies from the mesh
void append_energies(
+ const Mesh *mesh,
const Options *options,
SolverData *data,
std::vector<Eigen::Triplet<double> > &D_triplets);
diff --git a/extern/softbody/src/admmpd_types.h b/extern/softbody/src/admmpd_types.h
index 55a15f3c44f..4f9ea607b53 100644
--- a/extern/softbody/src/admmpd_types.h
+++ b/extern/softbody/src/admmpd_types.h
@@ -58,20 +58,19 @@ struct Options {
};
struct SolverData {
- // Set from input
- Eigen::MatrixXi tets; // elements t x 4, copy from mesh
Eigen::MatrixXd x; // vertices, n x 3
Eigen::MatrixXd v; // velocity, n x 3
-
Eigen::MatrixXd x_start; // x at t=0 (and goal if k>0), n x 3
+ Eigen::MatrixXd x_prev; // x at k-1
Eigen::VectorXd m; // masses, n x 1
Eigen::MatrixXd z; // ADMM z variable
Eigen::MatrixXd u; // ADMM u aug lag with W inv
Eigen::MatrixXd M_xbar; // M*(x + dt v)
Eigen::MatrixXd Dx; // D * x
- Eigen::MatrixXd b; // M xbar + DtW2(z-u)
RowSparseMatrix<double> D; // reduction matrix
RowSparseMatrix<double> DtW2; // D'W'W
+ RowSparseMatrix<double> A; // M + DtW'WD
+ RowSparseMatrix<double> W; // weight matrix
double A_diag_max; // Max coeff of diag of A
RowSparseMatrix<double> A3_plus_PtP; // A + pk PtP replicated
diff --git a/intern/softbody/admmpd_api.cpp b/intern/softbody/admmpd_api.cpp
index c599f59391d..c4be926aa88 100644
--- a/intern/softbody/admmpd_api.cpp
+++ b/intern/softbody/admmpd_api.cpp
@@ -178,6 +178,7 @@ static inline int admmpd_init_with_tetgen(ADMMPDInterfaceData *iface, Object *ob
if (di
@@ Diff output truncated at 10240 characters. @@
More information about the Bf-blender-cvs
mailing list