[Bf-blender-cvs] [09ceada01bc] soc-2020-soft-body: fast svd
over0219
noreply at git.blender.org
Wed Jul 22 04:23:49 CEST 2020
Commit: 09ceada01bc41d50ca167221136d13e151f808fe
Author: over0219
Date: Tue Jul 21 21:21:06 2020 -0500
Branches: soc-2020-soft-body
https://developer.blender.org/rB09ceada01bc41d50ca167221136d13e151f808fe
fast svd
===================================================================
M extern/softbody/CMakeLists.txt
M extern/softbody/src/admmpd_collision.cpp
M extern/softbody/src/admmpd_energy.cpp
M extern/softbody/src/admmpd_energy.h
M extern/softbody/src/admmpd_solver.cpp
M extern/softbody/src/admmpd_types.h
A extern/softbody/src/svd/ImplicitQRSVD.h
A extern/softbody/src/svd/Tools.h
M intern/softbody/admmpd_api.cpp
M source/blender/blenkernel/intern/softbody.c
===================================================================
diff --git a/extern/softbody/CMakeLists.txt b/extern/softbody/CMakeLists.txt
index 8120f1ec336..73e3c58b616 100644
--- a/extern/softbody/CMakeLists.txt
+++ b/extern/softbody/CMakeLists.txt
@@ -60,6 +60,8 @@ set(SRC
src/admmpd_mesh.h
src/admmpd_mesh.cpp
src/admmpd_types.h
+ src/svd/ImplicitQRSVD.h
+ src/svd/Tools.h
${SDFGEN_SRC}
)
diff --git a/extern/softbody/src/admmpd_collision.cpp b/extern/softbody/src/admmpd_collision.cpp
index e38df33df06..cf36c3ce009 100644
--- a/extern/softbody/src/admmpd_collision.cpp
+++ b/extern/softbody/src/admmpd_collision.cpp
@@ -107,10 +107,14 @@ int EmbeddedMeshCollision::detect(
return 0;
// Do we even need to process collisions?
- if (!this->settings.floor_collision &&
- (!this->settings.obs_collision || !obsdata.has_obs()) &&
+ if ((!this->settings.obs_collision || !obsdata.has_obs()) &&
!this->settings.self_collision)
- return 0;
+ {
+ if (x1->col(1).minCoeff() > this->settings.floor_z)
+ {
+ return 0;
+ }
+ }
if (!tet_tree.root())
throw std::runtime_error("EmbeddedMeshCollision::Detect: No tree");
@@ -127,7 +131,6 @@ int EmbeddedMeshCollision::detect(
typedef struct {
const Collision::Settings *settings;
const Collision *collision;
- const TetMeshData *tetmesh;
const EmbeddedMesh *embmesh;
const Collision::ObstacleData *obsdata;
const Eigen::MatrixXd *x0;
@@ -197,7 +200,6 @@ int EmbeddedMeshCollision::detect(
DetectThreadData thread_data = {
.settings = &settings,
.collision = this,
- .tetmesh = nullptr,
.embmesh = mesh.get(),
.obsdata = &obsdata,
.x0 = x0,
diff --git a/extern/softbody/src/admmpd_energy.cpp b/extern/softbody/src/admmpd_energy.cpp
index f8da41d11e5..4271b182ac3 100644
--- a/extern/softbody/src/admmpd_energy.cpp
+++ b/extern/softbody/src/admmpd_energy.cpp
@@ -2,6 +2,8 @@
// Distributed under the MIT License.
#include "admmpd_energy.h"
+#include "admmpd_types.h"
+#include "svd/ImplicitQRSVD.h"
#include <iostream>
#include <Eigen/Eigenvalues>
@@ -26,27 +28,49 @@ void EnergyTerm::signed_svd(
Eigen::Matrix<double,3,1> &S,
Eigen::Matrix<double,3,3> &V)
{
- JacobiSVD<Matrix3d> svd(A, ComputeFullU | ComputeFullV);
- S = svd.singularValues();
- U = svd.matrixU();
- V = svd.matrixV();
- Matrix3d J = Matrix3d::Identity();
- J(2,2) = -1.0;
- if (U.determinant() < 0.0)
- {
- U = U * J;
- S[2] = -S[2];
- }
- if (V.determinant() < 0.0)
- {
- Matrix3d Vt = V.transpose();
- Vt = J * Vt;
- V = Vt.transpose();
- S[2] = -S[2];
- }
+ JIXIE::singularValueDecomposition(A,U,S,V);
+ // Should replace this with
+// JacobiSVD<Matrix3d> svd(A, ComputeFullU | ComputeFullV);
+// S = svd.singularValues();
+// U = svd.matrixU();
+// V = svd.matrixV();
+// Matrix3d J = Matrix3d::Identity();
+// J(2,2) = -1.0;
+// if (U.determinant() < 0.0)
+// {
+// U = U * J;
+// S[2] = -S[2];
+// }
+// if (V.determinant() < 0.0)
+// {
+// Matrix3d Vt = V.transpose();
+// Vt = J * Vt;
+// V = Vt.transpose();
+// S[2] = -S[2];
+// }
}
void EnergyTerm::update(
+ int index,
+ int energyterm_type,
+ const Lame &lame,
+ double rest_volume,
+ double weight,
+ const Eigen::MatrixXd *x,
+ const Eigen::MatrixXd *Dx,
+ Eigen::MatrixXd *z,
+ Eigen::MatrixXd *u)
+{
+ switch (energyterm_type)
+ {
+ default:
+ case ENERGYTERM_TET: {
+ update_tet(index,lame,rest_volume,weight,x,Dx,z,u);
+ } break;
+ }
+} // end EnergyTerm::update
+
+void EnergyTerm::update_tet(
int index,
const Lame &lame,
double rest_volume,
@@ -147,10 +171,12 @@ void EnergyTerm::solve_prox(
const double eps = 1e-6;
int iter = 0;
- for(; iter<10; ++iter)
+ for(; iter<20; ++iter)
{
g.setZero();
energy_k = energy_density(lame,add_admm_pen,s0,s,&g); // e and g
+
+ // Converged because low gradient
if (g.norm() < eps)
break;
@@ -176,13 +202,17 @@ void EnergyTerm::solve_prox(
// Sometimes flattened tets will have a hard time
// uninverting, in which case they get linesearch
// blocked. There are ways around this, but for now
- // simply quitting is sufficient.
+ // simply quitting is fine. I've tried a few tricks
+ // to get past this, but in the end it doesn't really
+ // matter. The global step will smooth things over
+ // and get us in a better state next time ;)
if (ls_iter>=max_ls_iter)
{
s = s_prev;
break;
}
+ // Stopped making progress
if ((s-s_prev).norm() < eps)
break;
diff --git a/extern/softbody/src/admmpd_energy.h b/extern/softbody/src/admmpd_energy.h
index 55a15377a2e..953f1a7d81d 100644
--- a/extern/softbody/src/admmpd_energy.h
+++ b/extern/softbody/src/admmpd_energy.h
@@ -37,6 +37,17 @@ public:
// Updates the z and u variables for an element energy.
void update(
+ int index,
+ int energyterm_type,
+ const Lame &lame,
+ double rest_volume,
+ double weight,
+ const Eigen::MatrixXd *x,
+ const Eigen::MatrixXd *Dx,
+ Eigen::MatrixXd *z,
+ Eigen::MatrixXd *u);
+
+ void update_tet(
int index,
const Lame &lame,
double rest_volume,
diff --git a/extern/softbody/src/admmpd_solver.cpp b/extern/softbody/src/admmpd_solver.cpp
index a914aa9a480..599777a5c81 100644
--- a/extern/softbody/src/admmpd_solver.cpp
+++ b/extern/softbody/src/admmpd_solver.cpp
@@ -12,6 +12,7 @@
#include <stdio.h>
#include <iostream>
#include <unordered_map>
+#include <numeric>
#include "BLI_task.h" // threading
#include "BLI_assert.h"
@@ -109,7 +110,7 @@ void Solver::init_solve(
{
data->v.row(i) += dt*options->grav;
RowVector3d xbar_i = data->x.row(i) + dt*data->v.row(i);
- data->M_xbar.row(i) = data->m[i]*xbar_i;
+ data->M_xbar.row(i) = data->m[i]*xbar_i / (dt*dt);
data->x.row(i) = xbar_i; // initial guess
}
@@ -170,6 +171,7 @@ void Solver::solve_local_step(
lame.set_from_youngs_poisson(td->options->youngs,td->options->poisson);
EnergyTerm().update(
td->data->indices[i][0],
+ td->data->indices[i][2],
lame,
td->data->rest_volumes[i],
td->data->weights[i],
@@ -271,7 +273,7 @@ bool Solver::compute_matrices(
int ne = data->indices.size();
for (int i=0; i<ne; ++i)
{
- const Vector2i &idx = data->indices[i];
+ const Vector3i &idx = data->indices[i];
for (int j=0; j<idx[1]; ++j)
W2.coeffRef(idx[0]+j,idx[0]+j) = data->weights[i]*data->weights[i];
}
@@ -279,10 +281,10 @@ bool Solver::compute_matrices(
// Mass weighted Laplacian
data->D.resize(n_row_D,nx);
data->D.setFromTriplets(trips.begin(), trips.end());
- data->DtW2 = dt2 * data->D.transpose() * W2;
+ data->DtW2 = data->D.transpose() * W2;
data->A = data->DtW2 * data->D;
for (int i=0; i<nx; ++i)
- data->A.coeffRef(i,i) += data->m[i];
+ data->A.coeffRef(i,i) += data->m[i]/dt2;
data->ldltA.compute(data->A);
data->b.resize(nx,3);
data->b.setZero();
@@ -365,7 +367,7 @@ void Solver::append_energies(
data->energies_graph[ej].emplace(ek);
}
}
- data->indices.emplace_back(energy_index, energy_dim);
+ data->indices.emplace_back(energy_index, energy_dim, ENERGYTERM_TET);
energy_index += energy_dim;
}
diff --git a/extern/softbody/src/admmpd_types.h b/extern/softbody/src/admmpd_types.h
index f2bcfaacf90..bd12b1796d0 100644
--- a/extern/softbody/src/admmpd_types.h
+++ b/extern/softbody/src/admmpd_types.h
@@ -16,6 +16,10 @@
namespace admmpd {
template <typename T> using RowSparseMatrix = Eigen::SparseMatrix<T,Eigen::RowMajor>;
+enum EnergyTermType {
+ ENERGYTERM_TET = 0
+};
+
struct Options {
double timestep_s; // TODO: Figure out delta time from blender api
int max_admm_iters;
@@ -27,10 +31,10 @@ struct Options {
double min_res; // exit tolerance for global step
double youngs; // Young's modulus // TODO variable per-tet
double poisson; // Poisson ratio // TODO variable per-tet
- double density_kgm3; // density of deformables
+ double density_kgm3; // density of mesh
Eigen::Vector3d grav;
Options() :
- timestep_s(1.0/100.0),
+ timestep_s(1.0/24.0),
max_admm_iters(30),
max_cg_iters(10),
max_gs_iters(100),
@@ -40,22 +44,11 @@ struct Options {
min_res(1e-8),
youngs(10000000),
poisson(0.399),
- density_kgm3(1100),
+ density_kgm3(1522),
grav(0,0,-9.8)
{}
};
-// I think eventually I'll make the mesh
-// a virtual class with mapping functions.
-// That might clean up the API/interface a bit.
-// Will work out what we need for collisions and such first.
-
-struct TetMeshData {
- Eigen::MatrixXd x_rest; // verts at rest
- Eigen::MatrixXi faces; // surface elements, m x 3
- Eigen::MatrixXi tets; // internal elements, m x 4
-};
-
struct SolverData {
// Set from input
Eigen::MatrixXi tets; // elements t x 4, copy from mesh
@@ -95,7 +88,7 @@ struct SolverData {
} gsdata;
// Set in append_energies:
std::vector<std::set<int> > energies_graph; // per-vertex adjacency list (graph)
- std::vector<Eigen::Vector2i> indices; // per-energy index into D (row, num rows)
+ std::vector<Eigen::Vector3i> indices; // per-energy index into D (row, num rows, type)
std::vector<double> rest_volumes; // per-energy rest volume
std::vector<double> weights; // per-energy weights
};
diff --git a/extern/softbody/src/svd/ImplicitQRSVD.h b/extern/softbody/src/svd/ImplicitQRSVD.h
new file mode 100644
index 00000000000..98f5aa38d37
--- /dev/null
+++ b/extern/softbody/src/svd/ImplicitQRSVD.h
@@ -0,0 +1,878 @@
+/**
+Copyright (c) 2016 Theodore Gast, Chuyuan Fu, Chenfanfu Jiang, Joseph Teran
+
+Permission is hereby granted, free of charge, to any person obtaining a copy of
+this software and associated documentation files (the "Software"), to deal in
+the Software without restriction, including without limitation the rights to
+use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
+of the Software, and to permit persons to whom the Software is furnished to do
+so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+If the code is used in an article, the following paper shall be cited:
+ at techreport{qrsvd:2016,
+ title={Implicit-shifted Symmetric QR Singular Value Decomposition of 3x3 Matrices},
+
@@ Diff output truncated at 10240 characters. @@
More information about the Bf-blender-cvs
mailing list