Commit: 90ac7b5741c4be4a7923aa95bfb6e8c59896caf1
Author: mattoverby
Date: Tue Aug 11 14:06:26 2020 -0500
Branches: soc-2020-soft-body
https://developer.blender.org/rB90ac7b5741c4be4a7923aa95bfb6e8c59896caf1
cleaning up api
===================================================================
M extern/softbody/src/admmpd_energy.cpp
M extern/softbody/src/admmpd_energy.h
M extern/softbody/src/admmpd_linsolve.cpp
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/BKE_softbody.h
M source/blender/blenkernel/intern/object.c
M source/blender/blenkernel/intern/softbody.c
===================================================================
diff --git a/extern/softbody/src/admmpd_energy.cpp
b/extern/softbody/src/admmpd_energy.cpp
index eb4cb9c9534..ff0e31c3a5c 100644
--- a/extern/softbody/src/admmpd_energy.cpp
+++ b/extern/softbody/src/admmpd_energy.cpp
@@ -10,7 +10,7 @@
namespace admmpd {
using namespace Eigen;
-Lame::Lame() : m_material(ELASTIC_ARAP)
+Lame::Lame()
{
m_limit = Vector2d(-999,999);
set_from_youngs_poisson(10000000,0.399);
@@ -52,9 +52,9 @@ void EnergyTerm::signed_svd(
}
void EnergyTerm::update(
+ const Options *options,
int index,
int energyterm_type,
- const Lame &lame,
double rest_volume,
double weight,
const Eigen::MatrixXd *x,
@@ -66,17 +66,17 @@ void EnergyTerm::update(
{
default: break;
case ENERGYTERM_TET: {
- update_tet(index,lame,rest_volume,weight,x,Dx,z,u);
+ update_tet(options,index,rest_volume,weight,x,Dx,z,u);
} break;
case ENERGYTERM_TRIANGLE: {
- update_tri(index,lame,rest_volume,weight,x,Dx,z,u);
+ update_tri(options,index,rest_volume,weight,x,Dx,z,u);
} break;
}
} // end EnergyTerm::update
void EnergyTerm::update_tet(
+ const Options *options,
int index,
- const Lame &lame,
double rest_volume,
double weight,
const Eigen::MatrixXd *x,
@@ -84,6 +84,9 @@ void EnergyTerm::update_tet(
Eigen::MatrixXd *z,
Eigen::MatrixXd *u)
{
+ Lame lame;
+ lame.set_from_youngs_poisson(options->youngs,options->poisson);
+
(void)(x);
Matrix3d Dix = Dx->block<3,3>(index,0);
Matrix3d ui = u->block<3,3>(index,0);
@@ -94,7 +97,7 @@ void EnergyTerm::update_tet(
signed_svd(zi, U, S, V);
Vector3d s0 = S;
- switch (lame.m_material)
+ switch (options->elastic_material)
{
default:
case ELASTIC_ARAP:
@@ -109,7 +112,7 @@ void EnergyTerm::update_tet(
case ELASTIC_NH:
{
S = Vector3d::Ones();
- solve_prox(index,lame,s0,S);
+ solve_prox(options,index,lame,s0,S);
zi = U * S.asDiagonal() * V.transpose();
} break;
}
@@ -121,8 +124,8 @@ void EnergyTerm::update_tet(
} // end EnergyTerm::update
void EnergyTerm::update_tri(
+ const Options *options,
int index,
- const Lame &lame,
double rest_area,
double weight,
const Eigen::MatrixXd *x,
@@ -130,6 +133,10 @@ void EnergyTerm::update_tri(
Eigen::MatrixXd *z,
Eigen::MatrixXd *u)
{
+ Lame lame;
+ lame.set_from_youngs_poisson(options->youngs,options->poisson);
+
+ (void)(options);
typedef Matrix<double,2,3> Matrix23d;
// For we'll assume ARAP energy. If we need nonlinear energies
// in the future that's totally doable.
@@ -165,14 +172,17 @@ void EnergyTerm::update_tri(
}
int EnergyTerm::init_tet(
+ const Options *options,
int index,
- const Lame &lame,
const Eigen::Matrix<int,1,4> &prim,
const Eigen::MatrixXd *x,
double &volume,
double &weight,
std::vector< Eigen::Triplet<double> > &triplets )
{
+ Lame lame;
+ lame.set_from_youngs_poisson(options->youngs,options->poisson);
+
Matrix<double,3,3> edges;
edges.col(0) = x->row(prim[1]) - x->row(prim[0]);
edges.col(1) = x->row(prim[2]) - x->row(prim[0]);
@@ -203,14 +213,17 @@ int EnergyTerm::init_tet(
}
int EnergyTerm::init_triangle(
+ const Options *options,
int index,
- const Lame &lame,
const Eigen::RowVector3i &prim,
const Eigen::MatrixXd *x,
double &area,
double &weight,
std::vector< Eigen::Triplet<double> > &triplets)
{
+ Lame lame;
+ lame.set_from_youngs_poisson(options->youngs,options->poisson);
+
Matrix<double,3,2> edges;
edges.col(0) = x->row(prim[1]) - x->row(prim[0]);
edges.col(1) = x->row(prim[2]) - x->row(prim[0]);
@@ -244,10 +257,11 @@ int EnergyTerm::init_triangle(
}
void EnergyTerm::solve_prox(
- int index,
- const Lame &lame,
- const Eigen::Vector3d &s0,
- Eigen::Vector3d &s)
+ const Options *options,
+ int index,
+ const Lame &lame,
+ const Eigen::Vector3d &s0,
+ Eigen::Vector3d &s)
{
(void)(index);
Vector3d g; // gradient
@@ -262,18 +276,18 @@ void EnergyTerm::solve_prox(
for(; iter<20; ++iter)
{
g.setZero();
- energy_k = energy_density(lame,add_admm_pen,s0,s,&g); // e and g
+ energy_k = energy_density(options,lame,add_admm_pen,s0,s,&g);
// e and g
// Converged because low gradient
if (g.norm() < eps)
break;
- hessian_spd(lame,add_admm_pen,s,H);
+ hessian_spd(options,lame,add_admm_pen,s,H);
p = H.ldlt().solve(-g); // Newton step direction
s_prev = s;
s = s_prev + p;
- energy_k1 = energy_density(lame,add_admm_pen,s0,s);
+ energy_k1 = energy_density(options,lame,add_admm_pen,s0,s);
// Backtracking line search
double alpha = 1;
@@ -283,7 +297,7 @@ void EnergyTerm::solve_prox(
{
alpha *= 0.5;
s = s_prev + alpha*p;
- energy_k1 = energy_density(lame,add_admm_pen,s0,s);
+ energy_k1 =
energy_density(options,lame,add_admm_pen,s0,s);
ls_iter++;
}
@@ -312,6 +326,7 @@ void EnergyTerm::solve_prox(
} // end solve prox
double EnergyTerm::energy_density(
+ const Options *options,
const Lame &lame,
bool add_admm_penalty,
const Eigen::Vector3d &s0,
@@ -323,7 +338,7 @@ double EnergyTerm::energy_density(
// Compute energy and gradient
//
https://github.com/mattoverby/admm-elastic/blob/master/src/TetEnergyTerm.cpp
// I suppose I should add ARAP for completeness even though it's not
used
- switch (lame.m_material)
+ switch (options->elastic_material)
{
default: {
if (g != nullptr) g->setZero();
@@ -362,15 +377,16 @@ double EnergyTerm::energy_density(
} // end energy
void EnergyTerm::hessian_spd(
- const Lame &lame,
- bool add_admm_penalty,
- const Eigen::Vector3d &s,
- Eigen::Matrix3d &H)
+ const Options *options,
+ const Lame &lame,
+ bool add_admm_penalty,
+ const Eigen::Vector3d &s,
+ Eigen::Matrix3d &H)
{
static const Matrix3d I = Matrix3d::Identity();
// Compute specific Hessian
- switch (lame.m_material)
+ switch (options->elastic_material)
{
default:
{
@@ -421,5 +437,4 @@ void EnergyTerm::hessian_spd(
} // end hessian
-
} // end namespace mcl
diff --git a/extern/softbody/src/admmpd_energy.h
b/extern/softbody/src/admmpd_energy.h
index c26596c5ebf..1b88efa9f54 100644
--- a/extern/softbody/src/admmpd_energy.h
+++ b/extern/softbody/src/admmpd_energy.h
@@ -13,7 +13,6 @@ namespace admmpd {
class Lame {
public:
- int m_material;
double m_mu;
double m_lambda;
double m_bulk_mod;
@@ -33,9 +32,9 @@ public:
// Updates the z and u variables for an element energy.
void update(
+ const Options *options,
int index,
int energyterm_type,
- const Lame &lame,
double rest_volume,
double weight,
const Eigen::MatrixXd *x,
@@ -45,8 +44,8 @@ public:
// Updates the z and u variables for a tet
void update_tet(
+ const Options *options,
int index,
- const Lame &lame,
double rest_volume,
double weight,
const Eigen::MatrixXd *x,
@@ -56,8 +55,8 @@ public:
// Updates the z and u variables for a tri
void update_tri(
+ const Options *options,
int index,
- const Lame &lame,
double rest_area,
double weight,
const Eigen::MatrixXd *x,
@@ -67,8 +66,8 @@ public:
// Initializes tet energy, returns num rows for D
int init_tet(
+ const Options *options,
int index,
- const Lame &lame,
const Eigen::RowVector4i &prim,
const Eigen::MatrixXd *x,
double &volume,
@@ -77,8 +76,8 @@ public:
// Initializes facet energy, returns num rows for D
int init_triangle(
+ const Options *options,
int index,
- const Lame &lame,
const Eigen::RowVector3i &prim,
const Eigen::MatrixXd *x,
double &area,
@@ -87,6 +86,7 @@ public:
// Solves proximal energy function
void solve_prox(
+ const Options *options,
int index,
const Lame &lame,
const Eigen::Vector3d &s0,
@@ -95,6 +95,7 @@ public:
// Returns gradient and energy (+ADMM penalty)
// of a material evaluated at s (singular values).
double energy_density(
+ const Options *options,
const Lame &lame,
bool add_admm_penalty,
const Eigen::Vector3d &s0,
@@ -104,6 +105,7 @@ public:
// Returns the Hessian of a material at S
// projected to nearest SPD
void hessian_spd(
+ const Options *options,
const Lame &lame,
bool add_admm_penalty,
const Eigen::Vector3d &s,
diff --git a/extern/softbody/src/admmpd_linsolve.cpp
b/extern/softbody/src/admmpd_linsolve.cpp
index c54ad16f89f..bda806b2e2a 100644
--- a/extern/softbody/src/admmpd_linsolve.cpp
+++ b/extern/softbody/src/admmpd_linsolve.cpp
@@ -257,7 +257,6 @@ void ConjugateGradients::solve(
C.resize(1,nx*3);
d = VectorXd::Zero(1);
}
-
// Compute RHS
data->ls.rhs.noalias() =
@@ -328,17 +327,20 @@ void ConjugateGradients::apply_preconditioner(
Eigen::MatrixXd &x,
const Eigen::MatrixXd &b)
{
- x = data->ls.ldlt_A_PtP->solve(b);
-// x = m_ldlt.cholesky()->solve(b);
-/*
+ // Only bother with parallel solve if we have large dof
+ if (x.rows()<10000)
+ {
+ x = data->ls.ldlt_A_PtP->solve(b);
+ return;
+ }
+
BLI_assert(b.cols()==3);
if (x.rows() != b.rows())
x.resize(b.rows(),3);
- Cholesky *chol = m_ldlt.cholesky();
+ Cholesky *chol = data->ls.ldlt_A_PtP.get();
const auto & linsolve = [&x,&b,&chol](int col)
{
- if (col < 0 || col > 2) { return; }
x.col(col) = chol->solve(b.col(col));
};
@@ -351,8 +353,8 @@ void ConjugateGradients::apply_preconditioner(
if (pool[i].joinable())
pool[i].join();
}
-*/
-}
+
+} // end apply preconditioner
#if 0
void GaussSeidel::solve(
@@ -360,7 +362,6 @@ void GaussSeidel::solve(
SolverData *data,
Collision *collision)
{
- init_solve(options,data,collision);
MatrixXd dx(data->x.rows(),3);
dx.setZero();
@@ -377,7 +378,7 @@ void GaussSeidel::solve(
GaussSeidelThreadData thread_data = {
.iter = 0,
.color = 0,
- .colors = &data->gsdata.A3_plus_CtC_colors,
+ .colors = &data->ls.A_colors,
.options = options,
.data = data,
.collision = collision,
@@ -398,7 +399,7 @@ void GaussSeidel::solve(
Vector3d inv_aii(0,0,0);
for (int j=0; j<3; ++j)
{
- InnerIter rit(td->data->gsdata.A3_CtC_PtP, idx*3+j);
+ InnerIter rit(td->data->ls.A3_CtC_PtP, idx*3+j);
for (; rit; ++rit)
{
double v = rit.value();
@@ -420,7 +421,7 @@ void GaussSeidel::solve(
} // end loop segment
// Update x
- Vector3d bi = td->data->gsdata.b3_Ctd_Ptx.segment<3>(idx*3);
+ Vector3d bi = ...;
//Vector3d bi = td->data->b.
@@ Diff output truncated at 10240 characters. @@
_______________________________________________
Bf-blender-cvs mailing list
[email protected]
https://lists.blender.org/mailman/listinfo/bf-blender-cvs
