[Bf-blender-cvs] [8af51750b24] soc-2020-soft-body: working on mcgs
over0219
noreply at git.blender.org
Wed Jul 1 00:07:33 CEST 2020
Commit: 8af51750b24ea037c5ad64641d79b05a60448d11
Author: over0219
Date: Tue Jun 30 17:07:29 2020 -0500
Branches: soc-2020-soft-body
https://developer.blender.org/rB8af51750b24ea037c5ad64641d79b05a60448d11
working on mcgs
===================================================================
M extern/softbody/src/admmpd_collision.cpp
M extern/softbody/src/admmpd_collision.h
M extern/softbody/src/admmpd_embeddedmesh.cpp
M extern/softbody/src/admmpd_embeddedmesh.h
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
===================================================================
diff --git a/extern/softbody/src/admmpd_collision.cpp b/extern/softbody/src/admmpd_collision.cpp
index bbda72f2428..b7e835a9185 100644
--- a/extern/softbody/src/admmpd_collision.cpp
+++ b/extern/softbody/src/admmpd_collision.cpp
@@ -20,8 +20,7 @@ VFCollisionPair::VFCollisionPair() :
p_is_obs(0), // 0 or 1
q_idx(-1), // face
q_is_obs(0), // 0 or 1
- pt_on_q(0,0,0),
- q_n(0,0,0)
+ pt_on_q(0,0,0)
{}
void Collision::set_obstacles(
@@ -89,7 +88,6 @@ void Collision::detect_discrete_vf(
pair.q_idx = -1;
pair.q_is_obs = 1;
pair.pt_on_q = Vector3d(pt[0],pt[1],floor_z);
- pair.q_n = Vector3d(0,0,1);
}
// Any faces to detect against?
@@ -133,18 +131,17 @@ void Collision::detect_discrete_vf(
pair.pt_on_q = find_nearest_tri.output.pt_on_tri;
// Compute face normal
- BLI_assert(pair.q_idx >= 0);
- BLI_assert(pair.q_idx < mesh_tris->rows());
- RowVector3i tri_inds = mesh_tris->row(pair.q_idx);
- BLI_assert(tri_inds[0]>=0 && tri_inds[0]<mesh_x->rows());
- BLI_assert(tri_inds[1]>=0 && tri_inds[1]<mesh_x->rows());
- BLI_assert(tri_inds[2]>=0 && tri_inds[2]<mesh_x->rows());
- Vector3d tri[3] = {
- mesh_x->row(tri_inds[0]),
- mesh_x->row(tri_inds[1]),
- mesh_x->row(tri_inds[2])
- };
- pair.q_n = ((tri[1]-tri[0]).cross(tri[2]-tri[0])).normalized();
+// BLI_assert(pair.q_idx >= 0);
+// BLI_assert(pair.q_idx < mesh_tris->rows());
+// RowVector3i tri_inds = mesh_tris->row(pair.q_idx);
+// BLI_assert(tri_inds[0]>=0 && tri_inds[0]<mesh_x->rows());
+// BLI_assert(tri_inds[1]>=0 && tri_inds[1]<mesh_x->rows());
+// BLI_assert(tri_inds[2]>=0 && tri_inds[2]<mesh_x->rows());
+// Vector3d tri[3] = {
+// mesh_x->row(tri_inds[0]),
+// mesh_x->row(tri_inds[1]),
+// mesh_x->row(tri_inds[2])
+// };
// std::stringstream ss;
// ss << "\nhit:" <<
@@ -247,12 +244,10 @@ int EmbeddedMeshCollision::detect(
return vf_pairs.size();
} // end detect
-void EmbeddedMeshCollision::jacobian(
+void EmbeddedMeshCollision::linearize(
const Eigen::MatrixXd *x,
- std::vector<Eigen::Triplet<double> > *trips_x,
- std::vector<Eigen::Triplet<double> > *trips_y,
- std::vector<Eigen::Triplet<double> > *trips_z,
- std::vector<double> *l)
+ std::vector<Eigen::Triplet<double> > *trips,
+ std::vector<double> *d)
{
BLI_assert(x != NULL);
BLI_assert(x->cols() == 3);
@@ -261,16 +256,45 @@ void EmbeddedMeshCollision::jacobian(
if (np==0)
return;
- l->reserve((int)l->size() + np);
- trips_x->reserve((int)trips_x->size() + np*4);
- trips_y->reserve((int)trips_y->size() + np*4);
- trips_z->reserve((int)trips_z->size() + np*4);
+ d->reserve((int)d->size() + np);
+ trips->reserve((int)trips->size() + np*3*4);
for (int i=0; i<np; ++i)
{
VFCollisionPair &pair = vf_pairs[i];
int emb_p_idx = pair.p_idx;
+
+ //p_idx(-1), // point
+ //p_is_obs(0), // 0 or 1
+ //q_idx(-1), // face
+ //q_is_obs(0), // 0 or 1
+// pt_on_q(0,0,0)
+
+ // Compute normal of triangle at x
+ Vector3d q_n(0,0,0);
+ RowVector3i q_inds(-1,-1,-1);
+ std::vector<Vector3d> q_tris;
+ if (pair.q_is_obs)
+ {
+ // Special case, floor
+ if (pair.q_idx == -1)
+ {
+ q_n = Vector3d(0,0,1);
+ }
+ else
+ {
+ q_inds = obsdata.F.row(pair.q_idx);
+ q_tris = {
+ obsdata.V1.row(q_inds[0]),
+ obsdata.V1.row(q_inds[1]),
+ obsdata.V1.row(q_inds[2])
+ };
+ q_n = ((q_tris[1]-q_tris[0]).cross(q_tris[2]-q_tris[0]));
+ q_n.normalize();
+ }
+ }
+
// TODO: obstacle point intersecting mesh
if (pair.p_is_obs)
{
@@ -283,24 +307,14 @@ void EmbeddedMeshCollision::jacobian(
RowVector4d bary = mesh->barys.row(emb_p_idx);
int tet_idx = mesh->vtx_to_tet[emb_p_idx];
RowVector4i tet = mesh->tets.row(tet_idx);
-
- // Okay this is pretty ugly. I'm going to have to think about
- // whether or not this is reasonable.
- for (int j=0; j<3; ++j)
+ int c_idx = d->size();
+ d->emplace_back(q_n.dot(pair.pt_on_q));
+ for (int j=0; j<4; ++j)
{
- int c_idx = l->size();
- for (int k=0; k<4; ++k)
- {
- if (j==0)
- trips_x->emplace_back(c_idx, tet[k], bary[k]*pair.q_n[0]);
- if (j==1)
- trips_y->emplace_back(c_idx, tet[k], bary[k]*pair.q_n[1]);
- if (j==2)
- trips_z->emplace_back(c_idx, tet[k], bary[k]*pair.q_n[2]);
- }
- l->emplace_back(pair.pt_on_q[j]*pair.q_n[j]);
+ trips->emplace_back(c_idx, tet[j]*3+0, bary[j]*q_n[0]);
+ trips->emplace_back(c_idx, tet[j]*3+1, bary[j]*q_n[1]);
+ trips->emplace_back(c_idx, tet[j]*3+2, bary[j]*q_n[2]);
}
-
continue;
}
diff --git a/extern/softbody/src/admmpd_collision.h b/extern/softbody/src/admmpd_collision.h
index f9c6363b89f..2ab7d343c05 100644
--- a/extern/softbody/src/admmpd_collision.h
+++ b/extern/softbody/src/admmpd_collision.h
@@ -52,15 +52,11 @@ public:
// Special case for floor since it's common.
virtual void set_floor(double z) = 0;
- // Linearize the constraints and return Jacobian tensor.
- // Constraints are linearized about x for constraint
- // K x = l
- virtual void jacobian(
+ // Linearize the constraints and return Jacobian.
+ virtual void linearize(
const Eigen::MatrixXd *x,
- std::vector<Eigen::Triplet<double> > *trips_x,
- std::vector<Eigen::Triplet<double> > *trips_y,
- std::vector<Eigen::Triplet<double> > *trips_z,
- std::vector<double> *l) = 0;
+ std::vector<Eigen::Triplet<double> > *trips,
+ std::vector<double> *d) = 0;
// Given a point and a surface mesh,
// perform discrete collision and create
@@ -83,7 +79,8 @@ class EmbeddedMeshCollision : public Collision {
public:
EmbeddedMeshCollision(const EmbeddedMeshData *mesh_) :
mesh(mesh_),
- floor_z(-std::numeric_limits<double>::max())
+// floor_z(-std::numeric_limits<double>::max())
+ floor_z(0)
{}
// A floor is so common that it makes sense to hard
@@ -100,12 +97,10 @@ public:
// Linearizes the collision pairs about x
// for the constraint Kx=l
- void jacobian(
+ void linearize(
const Eigen::MatrixXd *x,
- std::vector<Eigen::Triplet<double> > *trips_x,
- std::vector<Eigen::Triplet<double> > *trips_y,
- std::vector<Eigen::Triplet<double> > *trips_z,
- std::vector<double> *l);
+ std::vector<Eigen::Triplet<double> > *trips,
+ std::vector<double> *d);
protected:
// A ptr to the embedded mesh data
diff --git a/extern/softbody/src/admmpd_embeddedmesh.cpp b/extern/softbody/src/admmpd_embeddedmesh.cpp
index 68cd38b27d5..22ac59fef0c 100644
--- a/extern/softbody/src/admmpd_embeddedmesh.cpp
+++ b/extern/softbody/src/admmpd_embeddedmesh.cpp
@@ -67,7 +67,8 @@ bool EmbeddedMesh::generate(
const Eigen::MatrixXd &V, // embedded verts
const Eigen::MatrixXi &F, // embedded faces
EmbeddedMeshData *emb_mesh, // where embedding is stored
- Eigen::MatrixXd *x_tets) // lattice vertices, n x 3
+ Eigen::MatrixXd *x_tets, // lattice vertices, n x 3
+ bool trim_lattice)
{
// How big the grid cells are as a fraction
// of the total mesh.
@@ -191,9 +192,12 @@ bool EmbeddedMesh::generate(
.tets = &tets,
.keep_tet = &keep_tet
};
- TaskParallelSettings settings;
- BLI_parallel_range_settings_defaults(&settings);
- BLI_task_parallel_range(0, nt0, &thread_data, parallel_keep_tet, &settings);
+ if (trim_lattice)
+ {
+ TaskParallelSettings settings;
+ BLI_parallel_range_settings_defaults(&settings);
+ BLI_task_parallel_range(0, nt0, &thread_data, parallel_keep_tet, &settings);
+ }
// Loop over tets and remove as needed.
// Mark referenced vertices to compute a mapping.
diff --git a/extern/softbody/src/admmpd_embeddedmesh.h b/extern/softbody/src/admmpd_embeddedmesh.h
index 85c64c7b0a5..eb7583b93b4 100644
--- a/extern/softbody/src/admmpd_embeddedmesh.h
+++ b/extern/softbody/src/admmpd_embeddedmesh.h
@@ -15,7 +15,8 @@ public:
const Eigen::MatrixXd &V, // embedded verts
const Eigen::MatrixXi &F, // embedded faces
EmbeddedMeshData *emb_mesh, // where embedding is stored
- Eigen::MatrixXd *x_tets); // lattice vertices, n x 3
+ Eigen::MatrixXd *x_tets, // lattice vertices, n x 3
+ bool trim_lattice = true); // remove elements outside embedded volume
// Returns the vtx mapped from x/v and tets
Eigen::Vector3d get_mapped_vertex(
diff --git a/extern/softbody/src/admmpd_linsolve.cpp b/extern/softbody/src/admmpd_linsolve.cpp
index 5b3a62dc7bb..0d40c49ef4d 100644
--- a/extern/softbody/src/admmpd_linsolve.cpp
+++ b/extern/softbody/src/admmpd_linsolve.cpp
@@ -4,12 +4,154 @@
#include "admmpd_linsolve.h"
#include <numeric>
#include <iostream>
+#include <set>
#include <unordered_set>
+#include <random>
#include "BLI_assert.h"
+#include "BLI_task.h"
namespace admmpd {
using namespace Eigen;
+// Makes full n3 x n3 matrix
+static inline void make_n3(
+ const RowSparseMatrix<double> &A,
+ RowSparseMatrix<double> &A3)
+{
+ int na = A.rows();
+ SparseMatrix<double> A_cm = A; // A to CM
+ SparseMatrix<double> A3_cm(na*3, na*3);
+ A3_cm.reserve(A_cm.nonZeros()*3);
+ int cols = A_cm.rows();
+ for(int c=0; c<cols; ++c)
+ {
+ for(int dim=0; dim<3; ++dim)
+ {
+ int col = c*3+dim;
+ A3_cm.startVec(col);
+ for(SparseMatrix<double>::InnerIterator itA(A_cm,c); itA; ++itA)
+ A3_cm.insertBack((itA.row()*3+dim), col) = itA.value();
+ }
+ }
+ A3_cm.finalize();
+ A3_cm.makeCompressed();
+ A3 = A3_cm; // to RowMajor
+} // end make_n3
+
+void ConjugateGradients::solve(
+ const Options *options,
+ SolverData *data)
+{
+ BLI_assert(data != NULL);
+ BLI_assert(options != NULL);
+ int nx = data->x.rows();
+ BLI_assert(nx > 0);
+ BLI_assert(data->A.rows() == nx);
+ BLI_assert(data->A.cols() == nx);
+ 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());
+
+ // 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)
+ {
+ data->x = data->ldltA.solve(data->b);
+ return;
+ }
+
+ // Resize data associated with Conjugate-Gradient
+ admmpd::SolverData::GlobalStepData *gsdata = &data->gsdata;
+ if (gsdata->r.rows() != nx*3)
+ {
+ gsdata->
@@ Diff output truncated at 10240 characters. @@
More information about the Bf-blender-cvs
mailing list