Xpbd cloth and fascia simulation

build.sh

@@ -1,50 +1,50 @@
-[ "$1" ] && [ $1 = Debug ] && flag="Debug" || flag="Release"
-echo $flag
+# [ "$1" ] && [ $1 = Debug ] && flag="Debug" || flag="Release"
+# echo $flag
 
-cmake -B build -DCMAKE_BUILD_TYPE=$flag \
-    -DZENO_MULTIPROCESS=ON \
-    -DZENO_WITH_zenvdb:BOOL=ON \
-    -DZENO_SYSTEM_OPENVDB=OFF \
-    -DZENO_WITH_ZenoFX:BOOL=ON \
-    -DZENO_ENABLE_OPTIX:BOOL=ON \
-    -DZENO_WITH_FBX:BOOL=ON \
-    -DZENO_WITH_Alembic:BOOL=ON \
-    -DZENO_WITH_MeshSubdiv:BOOL=ON \
-    -DCMAKE_EXPORT_COMPILE_COMMANDS=ON \
-
-# cmake -B build -DCMAKE_BUILD_TYPE=Release \
+# cmake -B build -DCMAKE_BUILD_TYPE=$flag \
+#     -DZENO_MULTIPROCESS=ON \
+#     -DZENO_WITH_zenvdb:BOOL=ON \
+#     -DZENO_SYSTEM_OPENVDB=OFF \
 #     -DZENO_WITH_ZenoFX:BOOL=ON \
 #     -DZENO_ENABLE_OPTIX:BOOL=ON \
-#     -DZENO_SYSTEM_OPENVDB=OFF \
-#     -DZENO_WITH_zenvdb:BOOL=ON \
-#     -DZENOFX_ENABLE_OPENVDB:BOOL=ON \
-#     -DZENOFX_ENABLE_LBVH:BOOL=ON \
-#     -DZENO_WITH_FastFLIP:BOOL=ON \
-#     -DZENO_WITH_FEM:BOOL=ON \
-#     -DZENO_WITH_Rigid:BOOL=ON \
-#     -DZENO_WITH_cgmesh:BOOL=ON \
-#     -DZENO_WITH_oldzenbase:BOOL=ON \
-#     -DZENO_WITH_TreeSketch:BOOL=ON \
-#     -DZENO_WITH_Skinning:BOOL=ON \
-#     -DZENO_WITH_Euler:BOOL=ON \
-#     -DZENO_WITH_Functional:BOOL=ON \
-#     -DZENO_WITH_LSystem:BOOL=ON \
-#     -DZENO_WITH_mesher:BOOL=ON \
-#     -DZENO_WITH_Alembic:BOOL=ON \
 #     -DZENO_WITH_FBX:BOOL=ON \
-#     -DZENO_WITH_DemBones:BOOL=ON \
-#     -DZENO_WITH_SampleModel:BOOL=ON \
-#     -DZENO_WITH_CalcGeometryUV:BOOL=ON \
+#     -DZENO_WITH_Alembic:BOOL=ON \
 #     -DZENO_WITH_MeshSubdiv:BOOL=ON \
-#     -DZENO_WITH_Audio:BOOL=ON \
-#     -DZENO_WITH_PBD:BOOL=ON \
-#     -DZENO_WITH_GUI:BOOL=ON \
-#     -DZENO_WITH_ImgCV:BOOL=ON \
-#     -DCMAKE_EXPORT_COMPILE_COMMANDS=ON
+#     -DCMAKE_EXPORT_COMPILE_COMMANDS=ON \
+
+cmake -B build -DCMAKE_BUILD_TYPE=Release \
+    -DZENO_WITH_ZenoFX:BOOL=ON \
+    -DZENO_ENABLE_OPTIX:BOOL=ON \
+    -DZENO_SYSTEM_OPENVDB=OFF \
+    -DZENO_WITH_zenvdb:BOOL=ON \
+    -DZENOFX_ENABLE_OPENVDB:BOOL=ON \
+    -DZENOFX_ENABLE_LBVH:BOOL=ON \
+    -DZENO_WITH_FastFLIP:BOOL=ON \
+    -DZENO_WITH_FEM:BOOL=ON \
+    -DZENO_WITH_Rigid:BOOL=ON \
+    -DZENO_WITH_cgmesh:BOOL=ON \
+    -DZENO_WITH_oldzenbase:BOOL=ON \
+    -DZENO_WITH_TreeSketch:BOOL=ON \
+    -DZENO_WITH_Skinning:BOOL=ON \
+    -DZENO_WITH_Euler:BOOL=ON \
+    -DZENO_WITH_Functional:BOOL=ON \
+    -DZENO_WITH_LSystem:BOOL=ON \
+    -DZENO_WITH_mesher:BOOL=ON \
+    -DZENO_WITH_Alembic:BOOL=ON \
+    -DZENO_WITH_FBX:BOOL=ON \
+    -DZENO_WITH_DemBones:BOOL=ON \
+    -DZENO_WITH_SampleModel:BOOL=ON \
+    -DZENO_WITH_CalcGeometryUV:BOOL=ON \
+    -DZENO_WITH_MeshSubdiv:BOOL=ON \
+    -DZENO_WITH_Audio:BOOL=ON \
+    -DZENO_WITH_PBD:BOOL=ON \
+    -DZENO_WITH_GUI:BOOL=ON \
+    -DZENO_WITH_ImgCV:BOOL=ON \
+    -DCMAKE_EXPORT_COMPILE_COMMANDS=ON
 
 cmake --build build --parallel $(nproc) \
 
 # ln -s ./build/compile_commands.json ./
-cp ./build/compile_commands.json ./
+# cp ./build/compile_commands.json ./
 
-rm -rf /var/tmp/OptixCache_$USER/
+# rm -rf /var/tmp/OptixCache_$USER/

projects/CUDA/Structures.hpp

@@ -163,6 +163,7 @@ struct ZenoParticles : IObjectClone<ZenoParticles> {
     static constexpr auto s_surfVertTag = "surfVerts";
     static constexpr auto s_surfHalfEdgeTag = "surfHalfEdges";
     static constexpr auto s_tetHalfFacetTag = "tetHalfFacets";
+    static constexpr auto s_surfBoundaryEdgeTag = "surfBoundaryEdges";
 
     ZenoParticles() = default;
     ~ZenoParticles() = default;

projects/CUDA/Utils.hpp

@@ -65,6 +65,165 @@ retrieve_bounding_volumes(Pol &pol, const TileVecT &vtemp,
   return ret;
 }
 
+
+template <int codim, typename Pol, typename TileVecT>
+zs::Vector<typename ZenoParticles::lbvh_t::Box>
+retrieve_bounding_volumes(Pol &pol, const TileVecT &vtemp,
+                          float thickness = 0.f,
+                          const zs::SmallString &xTag = "xn") {
+  using namespace zs;
+  using bv_t = typename ZenoParticles::lbvh_t::Box;
+  // static_assert(codim >= 1 && codim <= 4, "invalid co-dimension!\n");
+  constexpr auto space = Pol::exec_tag::value;
+#if ZS_ENABLE_CUDA && defined(__CUDACC__)
+  // ZS_LAMBDA -> __device__
+  static_assert(space == execspace_e::cuda,
+                "specialized policy and compiler not match");
+#else
+  static_assert(space != execspace_e::cuda,
+                "specialized policy and compiler not match");
+#endif
+  size_t nm_elms = vtemp.size() / codim;
+
+  Vector<bv_t> ret{vtemp.get_allocator(), nm_elms};
+  pol(zs::range(nm_elms), [bvs = proxy<space>(ret),
+                               vtemp = proxy<space>({}, vtemp), xTag,
+                               thickness] ZS_LAMBDA(int ei) mutable {
+    auto x0 = vtemp.template pack<3>(xTag, ei * codim);
+    bv_t bv{x0, x0};
+    for (int d = 1; d != codim; ++d)
+      merge(bv, vtemp.template pack<3>(xTag, ei * codim + d));
+    bv._min -= thickness / 2;
+    bv._max += thickness / 2;
+    bvs[ei] = bv;
+  });
+  return ret;
+}
+
+template <typename Pol, typename TileVecT>
+zs::Vector<typename ZenoParticles::lbvh_t::Box>
+retrieve_bounding_volumes(Pol &pol, const TileVecT &vtemp,
+                          float radius = 0.f,
+                          const zs::SmallString &xTag = "xn") {
+  using namespace zs;
+  using bv_t = typename ZenoParticles::lbvh_t::Box;
+  constexpr auto space = Pol::exec_tag::value;
+#if ZS_ENABLE_CUDA && defined(__CUDACC__)
+  // ZS_LAMBDA -> __device__
+  static_assert(space == execspace_e::cuda,
+                "specialized policy and compiler not match");
+#else
+  static_assert(space != execspace_e::cuda,
+                "specialized policy and compiler not match");
+#endif
+  Vector<bv_t> ret{vtemp.get_allocator(), vtemp.size()};
+  pol(zs::range(vtemp.size()), [bvs = proxy<space>(ret),
+                               vtemp = proxy<space>({}, vtemp),
+                               xTag = xTag,
+                               radius = radius] ZS_LAMBDA(int vi) mutable {
+    auto x0 = vtemp.template pack<3>(xTag, vi);
+    bv_t bv{x0, x0};
+    bv._min -= radius;
+    bv._max += radius;
+    bvs[vi] = bv;
+  });
+  return ret;
+}
+
+template <typename Pol, typename TileVecT>
+void retrieve_bounding_volumes(Pol &pol, const TileVecT &vtemp,
+                          zs::Vector<ZenoParticles::lbvh_t::Box>& ret,
+                          float radius = 0.f,
+                          const zs::SmallString &xTag = "xn") {
+  using namespace zs;
+  using bv_t = typename ZenoParticles::lbvh_t::Box;
+  constexpr auto space = Pol::exec_tag::value;
+#if ZS_ENABLE_CUDA && defined(__CUDACC__)
+  // ZS_LAMBDA -> __device__
+  static_assert(space == execspace_e::cuda,
+                "specialized policy and compiler not match");
+#else
+  static_assert(space != execspace_e::cuda,
+                "specialized policy and compiler not match");
+#endif
+  // Vector<bv_t> ret{vtemp.get_allocator(), vtemp.size()};
+  ret.resize(vtemp.size());
+  pol(zs::range(vtemp.size()), [bvs = proxy<space>(ret),
+                               vtemp = proxy<space>({}, vtemp),
+                               xTag = xTag,
+                               radius = radius] ZS_LAMBDA(int vi) mutable {
+    auto x0 = vtemp.template pack<3>(xTag, vi);
+    bv_t bv{x0, x0};
+    bv._min -= radius;
+    bv._max += radius;
+    bvs[vi] = bv;
+  });
+  // return ret;
+}
+
+template <typename Pol, typename TileVecT>
+zs::Vector<typename ZenoParticles::lbvh_t::Box>
+retrieve_bounding_volumes(Pol &pol, const TileVecT &vtemp,
+                          const zs::SmallString &rTag = "r",
+                          const zs::SmallString &xTag = "xn") {
+  using namespace zs;
+  using bv_t = typename ZenoParticles::lbvh_t::Box;
+  constexpr auto space = Pol::exec_tag::value;
+#if ZS_ENABLE_CUDA && defined(__CUDACC__)
+  // ZS_LAMBDA -> __device__
+  static_assert(space == execspace_e::cuda,
+                "specialized policy and compiler not match");
+#else
+  static_assert(space != execspace_e::cuda,
+                "specialized policy and compiler not match");
+#endif
+  Vector<bv_t> ret{vtemp.get_allocator(), vtemp.size()};
+  pol(zs::range(vtemp.size()), [bvs = proxy<space>(ret),
+                               vtemp = proxy<space>({}, vtemp),
+                               rTag,
+                               xTag] ZS_LAMBDA(int vi) mutable {
+    auto x0 = vtemp.template pack<3>(xTag, vi);
+    bv_t bv{x0, x0};
+    auto radius = vtemp(rTag,vi);
+    bv._min -= radius;
+    bv._max += radius;
+    bvs[vi] = bv;
+  });
+  return ret;
+}
+
+template <typename Pol, typename TileVecT>
+void retrieve_bounding_volumes(Pol &pol, const TileVecT &vtemp,
+                          zs::Vector<typename ZenoParticles::lbvh_t::Box>& ret,
+                          const zs::SmallString &rTag = "r",
+                          const zs::SmallString &xTag = "xn") {
+  using namespace zs;
+  using bv_t = typename ZenoParticles::lbvh_t::Box;
+  constexpr auto space = Pol::exec_tag::value;
+#if ZS_ENABLE_CUDA && defined(__CUDACC__)
+  // ZS_LAMBDA -> __device__
+  static_assert(space == execspace_e::cuda,
+                "specialized policy and compiler not match");
+#else
+  static_assert(space != execspace_e::cuda,
+                "specialized policy and compiler not match");
+#endif
+  // Vector<bv_t> ret{vtemp.get_allocator(), vtemp.size()};
+  ret.resize(vtemp.size());
+  pol(zs::range(vtemp.size()), [bvs = proxy<space>(ret),
+                               vtemp = proxy<space>({}, vtemp),
+                               rTag,
+                               xTag] ZS_LAMBDA(int vi) mutable {
+    auto x0 = vtemp.template pack<3>(xTag, vi);
+    bv_t bv{x0, x0};
+    auto radius = vtemp(rTag,vi);
+    bv._min -= radius;
+    bv._max += radius;
+    bvs[vi] = bv;
+  });
+  // return ret;
+}
+
 // for ccd
 template <typename Pol, typename TileVecT0, typename TileVecT1, int codim = 3>
 zs::Vector<typename ZenoParticles::lbvh_t::Box> retrieve_bounding_volumes(

projects/CuLagrange/CMakeLists.txt

@@ -35,10 +35,9 @@ target_sources(zeno PRIVATE
 
 # pbd
 target_sources(zeno PRIVATE
-  pbd/Constraints.cu
-  # pbd/SelfCollision.cu
+  pbd/ConstraintsBuilder.cu
+  pbd/ConstraintsSolver.cu
   pbd/KinematicCollision.cu
-  # pbd/Pipeline.cu
 )
 
 # fem-cloth
@@ -114,7 +113,7 @@ target_sources(zeno PRIVATE
   geometry/kernel/compute_characteristic_length.hpp
   geometry/kernel/topology.hpp
   geometry/kernel/calculate_facet_normal.hpp
-  geometry/kernel/calculate_bisector_normal.hpp
+  # geometry/kernel/calculate_bisector_normal.hpp
   geometry/kernel/tiled_vector_ops.hpp
   geometry/kernel/intersection.hpp
   geometry/SolveLaplacian.cu

projects/CuLagrange/fem/Ccds.hpp

@@ -92,6 +92,10 @@ ptccd(VecT p, VecT t0, VecT t1, VecT t2, VecT dp, VecT dt0, VecT dt1, VecT dt2,
   T gap = eta * (dist2_cur - thickness * thickness) / (dist_cur + thickness);
   T toc_prev = toc;
   toc = tStart;
+
+  if(gap < 0)
+    return true;
+
   int iter = 0;
   while (++iter < 1000000) {
     // while (true) {
@@ -137,6 +141,9 @@ eeccd(VecT ea0, VecT ea1, VecT eb0, VecT eb1, VecT dea0, VecT dea1, VecT deb0,
 
   T dist2_cur = dist2_ee_unclassified(ea0, ea1, eb0, eb1);
   T dFunc = dist2_cur - thickness * thickness;
+// suspecious code, check twice
+  if(dFunc <= 0)
+    return true;
   if (dFunc <= 0) {
     // since we ensured other place that all dist smaller than dHat are
     // positive, this must be some far away nearly parallel edges
@@ -151,10 +158,16 @@ eeccd(VecT ea0, VecT ea1, VecT eb0, VecT eb1, VecT dea0, VecT dea1, VecT deb0,
     }
     dFunc = dist2_cur - thickness * thickness;
   }
+
+
   T dist_cur = zs::sqrt(dist2_cur);
   T gap = eta * dFunc / (dist_cur + thickness);
   T toc_prev = toc;
   toc = tStart;
+
+  if(dFunc <= 0)
+    return true;
+
   int iter = 0;
   while (++iter < 1000000) {
     // while (true) {

projects/CuLagrange/fem/ClothCollision.cu

@@ -107,7 +107,7 @@ void ClothSystem::findBoundaryCellCollisionConstraints(zs::CudaExecutionPolicy &
                         biNrm *= -1;
                     return biNrm;
                 };
-                if (!COLLISION_UTILS::pointProjectsInsideTriangle(t0, t1, t2, p))
+                if (!COLLISION_UTILS::vertex_projects_inside_triangle(t0, t1, t2, p))
                     for (int i = 0; i != 3; ++i) {
                         auto bisector_normal = get_bisector_orient(i);
                         if (bisector_normal.dot(seg) < 0)