[Bf-blender-cvs] [67bc4d23b27] soc-2020-soft-body: updated embedded mesh
over0219
noreply at git.blender.org
Wed Jul 22 04:23:48 CEST 2020
Commit: 67bc4d23b2730aaf7fbe7278322e42947ff4c719
Author: over0219
Date: Tue Jul 21 16:56:34 2020 -0500
Branches: soc-2020-soft-body
https://developer.blender.org/rB67bc4d23b2730aaf7fbe7278322e42947ff4c719
updated embedded mesh
===================================================================
M extern/softbody/src/admmpd_bvh.h
M extern/softbody/src/admmpd_geom.cpp
M extern/softbody/src/admmpd_geom.h
M extern/softbody/src/admmpd_mesh.cpp
M extern/softbody/src/admmpd_mesh.h
M intern/softbody/admmpd_api.cpp
===================================================================
diff --git a/extern/softbody/src/admmpd_bvh.h b/extern/softbody/src/admmpd_bvh.h
index 88ff5393dc4..5905b676adb 100644
--- a/extern/softbody/src/admmpd_bvh.h
+++ b/extern/softbody/src/admmpd_bvh.h
@@ -100,6 +100,7 @@ public:
// (only create children if it contains prims) and does not
// create a cell if it is outside the mesh.
// ** Assumes a closed mesh and only defined for 3D
+ // eps is added to the min and max of the boxes.
void init(const MatrixXT *V, const Eigen::MatrixXi *F, int stopdepth);
// Returns bounding box of the root node
diff --git a/extern/softbody/src/admmpd_geom.cpp b/extern/softbody/src/admmpd_geom.cpp
index c03c5718383..f6a545276a6 100644
--- a/extern/softbody/src/admmpd_geom.cpp
+++ b/extern/softbody/src/admmpd_geom.cpp
@@ -44,7 +44,6 @@ bool geom::point_in_tet(
const Vector3d &c,
const Vector3d &d)
{
- using namespace Eigen;
auto check_face = [](
const Vector3d &point,
const Vector3d &p0,
@@ -64,6 +63,43 @@ bool geom::point_in_tet(
check_face(p, d, a, b, c);
}
+void geom::create_tets_from_box(
+ const Eigen::Vector3d &bmin,
+ const Eigen::Vector3d &bmax,
+ std::vector<Eigen::Vector3d> &verts,
+ std::vector<Eigen::RowVector4i> &tets)
+{
+ std::vector<Vector3d> v = {
+ // Top plane, clockwise looking down
+ bmax,
+ Vector3d(bmin[0], bmax[1], bmax[2]),
+ Vector3d(bmin[0], bmax[1], bmin[2]),
+ Vector3d(bmax[0], bmax[1], bmin[2]),
+ // Bottom plane, clockwise looking down
+ Vector3d(bmax[0], bmin[1], bmax[2]),
+ Vector3d(bmin[0], bmin[1], bmax[2]),
+ bmin,
+ Vector3d(bmax[0], bmin[1], bmin[2])
+ };
+ // Add vertices and get indices of the box
+ std::vector<int> b;
+ for(int i=0; i<8; ++i)
+ {
+ b.emplace_back(verts.size());
+ verts.emplace_back(v[i]);
+ }
+ // From the box, create five new tets
+ std::vector<RowVector4i> new_tets = {
+ RowVector4i( b[0], b[5], b[7], b[4] ),
+ RowVector4i( b[5], b[7], b[2], b[0] ),
+ RowVector4i( b[5], b[0], b[2], b[1] ),
+ RowVector4i( b[7], b[2], b[0], b[3] ),
+ RowVector4i( b[5], b[2], b[7], b[6] )
+ };
+ for(int i=0; i<5; ++i)
+ tets.emplace_back(new_tets[i]);
+}
+
// From Real-Time Collision Detection by Christer Ericson
template<typename T>
Eigen::Matrix<T,3,1> geom::point_triangle_barys(
@@ -342,6 +378,47 @@ bool geom::ray_triangle(
return true;
} // end ray - triangle test
+void geom::merge_close_vertices(
+ std::vector<Vector3d> &verts,
+ std::vector<RowVector4i> &tets,
+ double eps)
+{
+ int nv = verts.size();
+ std::vector<Vector3d> new_v(nv); // new verts
+ std::vector<int> idx(nv,0); // index mapping
+ std::vector<int> visited(nv,0);
+ int count = 0;
+ for (int i=0; i<nv; ++i)
+ {
+ if(!visited[i])
+ {
+ visited[i] = 1;
+ new_v[count] = verts[i];
+ idx[i] = count;
+ Vector3d vi = verts[i];
+ for (int j = i+1; j<nv; ++j)
+ {
+ if((verts[j]-vi).norm() < eps)
+ {
+ visited[j] = 1;
+ idx[j] = count;
+ }
+ }
+ count++;
+ }
+ }
+ new_v.resize(count);
+ verts = new_v;
+ int nt = tets.size();
+ for (int i=0; i<nt; ++i)
+ {
+ for (int j=0; j<4; ++j)
+ {
+ tets[i][j] = idx[tets[i][j]];
+ }
+ }
+}
+
//
// Compile template types
//
diff --git a/extern/softbody/src/admmpd_geom.h b/extern/softbody/src/admmpd_geom.h
index 200cdf4c68a..46da52e683c 100644
--- a/extern/softbody/src/admmpd_geom.h
+++ b/extern/softbody/src/admmpd_geom.h
@@ -5,6 +5,7 @@
#define ADMMPD_GEOM_H_
#include <Eigen/Geometry>
+#include <vector>
// Common geometry kernels
namespace admmpd {
@@ -25,6 +26,12 @@ static bool point_in_tet(
const Eigen::Vector3d &c,
const Eigen::Vector3d &d);
+static void create_tets_from_box(
+ const Eigen::Vector3d &bmin,
+ const Eigen::Vector3d &bmax,
+ std::vector<Eigen::Vector3d> &verts,
+ std::vector<Eigen::RowVector4i> &tets);
+
template<typename T>
static Eigen::Matrix<T,3,1> point_triangle_barys(
const Eigen::Matrix<T,3,1> &p,
@@ -72,6 +79,11 @@ static bool ray_triangle(
T &t_max,
Eigen::Matrix<T,3,1> *bary=nullptr);
+static void merge_close_vertices(
+ std::vector<Eigen::Vector3d> &verts,
+ std::vector<Eigen::RowVector4i> &tets,
+ double eps = 1e-12);
+
}; // class geom
} // namespace admmpd
diff --git a/extern/softbody/src/admmpd_mesh.cpp b/extern/softbody/src/admmpd_mesh.cpp
index c4225fa6dc1..1d4e6ee3d0a 100644
--- a/extern/softbody/src/admmpd_mesh.cpp
+++ b/extern/softbody/src/admmpd_mesh.cpp
@@ -7,10 +7,356 @@
#include <set>
#include <iostream>
#include "BLI_assert.h"
+#include "BLI_task.h"
namespace admmpd {
using namespace Eigen;
+static void gather_octree_tets(
+ Octree<double,3>::Node *node,
+ const MatrixXd *emb_V,
+ const MatrixXi *emb_F,
+ Eigen::VectorXi *emb_v_to_tet,
+ Eigen::MatrixXd *emb_barys,
+ AABBTree<double,3> *face_tree,
+ std::vector<Vector3d> &verts,
+ std::vector<RowVector4i> &tets
+ )
+{
+ if (node == nullptr)
+ return;
+
+ bool is_leaf = node->is_leaf();
+ bool has_prims = (int)node->prims.size()>0;
+ if (is_leaf)
+ {
+ Vector3d bmin = node->center-Vector3d::Ones()*node->halfwidth;
+ Vector3d bmax = node->center+Vector3d::Ones()*node->halfwidth;
+
+ // If we have primitives in the cell,
+ // create tets and compute embedding
+ if (has_prims)
+ {
+// int prev_tets = tets.size();
+ geom::create_tets_from_box(bmin,bmax,verts,tets);
+// int num_box_tets = (int)tets.size()-prev_tets;
+/*
+ // Loop over faces embedded in this box.
+ // Loop over tets created for this box.
+ // Find which tet the face vertex is embedded in.
+ int np = node->prims.size();
+ for (int i=0; i<np; ++i)
+ {
+ for (int j=0; j<3; ++j)
+ {
+ int v_idx = emb_F->operator()(i,j);
+ Vector3d v = emb_V->row(v_idx);
+ for (int k=0; k<num_box_tets; ++k)
+ {
+ const RowVector4i &t = tets[prev_tets+k];
+ if (geom::point_in_tet(v,verts[t[0]],verts[t[1]],verts[t[2]],verts[t[3]]))
+ {
+ emb_v_to_tet->operator[](v_idx) = prev_tets+k;
+ emb_barys->row(v_idx) = geom::point_tet_barys(v,
+ verts[t[0]],verts[t[1]],verts[t[2]],verts[t[3]]);
+ break;
+ }
+
+ }
+ }
+ }
+*/
+ }
+ else
+ {
+ // Otherwise, launch a ray
+ // to determine if we are inside or outside
+ // the mesh. If we're outside, don't create tets.
+ PointInTriangleMeshTraverse<double> pt_in_mesh(node->center,emb_V,emb_F);
+ face_tree->traverse(pt_in_mesh);
+ if (pt_in_mesh.output.is_inside())
+ geom::create_tets_from_box(bmin,bmax,verts,tets);
+ }
+ return;
+ }
+ for (int i=0; i<8; ++i)
+ {
+ gather_octree_tets(node->children[i],emb_V,emb_F,emb_v_to_tet,emb_barys,face_tree,verts,tets);
+ }
+
+} // end gather octree tets
+
+bool EmbeddedMesh::create(
+ const float *verts, // size nv*3
+ int nv,
+ const unsigned int *faces, // size nf*3
+ int nf,
+ const unsigned int *tets, // ignored
+ int nt)
+{
+ if (nv<=0 || verts == nullptr)
+ return false;
+ if (nf<=0 || faces == nullptr)
+ return false;
+
+ emb_V0.resize(nv,3);
+ for (int i=0; i<nv; ++i)
+ {
+ emb_V0(i,0) = verts[i*3+0];
+ emb_V0(i,1) = verts[i*3+1];
+ emb_V0(i,2) = verts[i*3+2];
+ }
+
+ emb_F.resize(nf,3);
+ std::vector<AlignedBox<double,3> > emb_leaves(nf);
+ for (int i=0; i<nf; ++i)
+ {
+ emb_F(i,0) = faces[i*3+0];
+ emb_F(i,1) = faces[i*3+1];
+ emb_F(i,2) = faces[i*3+2];
+ emb_leaves.emplace_back();
+ AlignedBox<double,3> &box = emb_leaves.back();
+ box.extend(emb_V0.row(emb_F(i,0)).transpose());
+ box.extend(emb_V0.row(emb_F(i,1)).transpose());
+ box.extend(emb_V0.row(emb_F(i,2)).transpose());
+ box.extend(box.min()-Vector3d::Ones()*1e-8);
+ box.extend(box.max()+Vector3d::Ones()*1e-8);
+ }
+
+ Octree<double,3> octree;
+ octree.init(&emb_V0,&emb_F,options.max_subdiv_levels);
+ emb_rest_facet_tree.init(emb_leaves);
+
+ emb_v_to_tet.resize(nv);
+ emb_v_to_tet.array() = -1;
+ emb_barys.resize(nv,4);
+ emb_barys.setZero();
+
+ std::vector<Vector3d> lat_verts;
+ std::vector<RowVector4i> lat_tets;
+ Octree<double,3>::Node *root = octree.root().get();
+ gather_octree_tets(
+ root,
+ &emb_V0,
+ &emb_F,
+ &emb_v_to_tet,
+ &emb_barys,
+ &emb_rest_facet_tree,
+ lat_verts,
+ lat_tets);
+ geom::merge_close_vertices(lat_verts,lat_tets);
+
+ int nlv = lat_verts.size();
+ lat_V0.resize(nlv,3);
+ for (int i=0; i<nlv; ++i)
+ {
+ for(int j=0; j<3; ++j){
+ lat_V0(i,j) = lat_verts[i][j];
+ }
+ }
+ int nlt = lat_tets.size();
+ lat_T.resize(nlt,4);
+ for(int i=0; i<nlt; ++i){
+ for(int j=0; j<4; ++j){
+ lat_T(i,j) = lat_tets[i][j];
+ }
+ }
+
+ compute_embedding();
+
+ auto return_error = [](const std::string &msg)
+ {
+ printf("EmbeddedMesh::generate create: %s\n", msg.c_str());
+ return false;
+ };
+
+ // Verify embedding is correct
+ for (int i=0; i<nv; ++i)
+ {
+ if (emb_v_to_tet[i]<0)
+ return return_error("Failed embedding");
+ if (std::abs(emb_barys.row(i).sum()-1.0)>1e-6)
+ {
+ std::stringstream ss; ss << emb_barys.row(i);
+ return return_error("Bad embedding barys: "+ss.str());
+ }
+ }
+
+ if (!emb_rest_facet_tree.root())
+ return return_error("Failed to create tree");
+ if (lat_V0.rows()==0)
+ return return_error("Failed to create verts");
+ if (lat_T.rows()==0)
+ return return_error("Failed to create tets");
+ if (emb_F.rows()==0)
+ return return_error("Did not set faces");
+ if (emb_V0.rows()==0)
+ return return_error("Did not set verts");
+
+ return true;
+}
+
+bool EmbeddedMesh::compute_embedding()
+{
+ struct FindTetThreadData {
+ AABBTree<double,3> *tree;
+ EmbeddedMesh *emb_mesh; // thread sets vtx_to_tet and barys
+ MatrixXd *lat_V0;
+ MatrixXi *lat_T;
+ MatrixXd *emb_barys;
+ VectorXi *emb_v_to_tet;
+ };
+
+ auto parallel_point_in_tet = [](
+ void *__restrict userdata,
+ const int i,
+ const TaskParallelTLS *__restrict tls)->void
+ {
+ (void)(tls);
+ FindTetThreadData *td = (FindTetThreadData*)userdata;
+ Vector3d pt = td->emb_mesh->rest_facet_verts().row(i);
+ PointInTetMeshTraverse<double> traverser(
+ pt,
+ td->lat_V0,
+ td->lat_T);
+ bool success = td->tree->traverse(traverser);
+ int tet_idx = traverser.output.prim;
+ if (success && tet_idx >= 0)
+ {
+ RowVector4i tet = td->emb_mesh->prims().row(tet_idx);
+ Vector3d t[4] = {
+ td->emb_mesh->rest_prim_verts().row(tet[0]),
+ td->emb_mesh->rest_prim_verts().row(tet[1])
@@ Diff output truncated at 10240 characters. @@
More information about the Bf-blender-cvs
mailing list