From d48e17ef5e2902279ef6422f8a5c6c462bb3f911 Mon Sep 17 00:00:00 2001 From: ShuliangLu Date: Thu, 13 Jul 2023 17:52:12 +0800 Subject: [PATCH 1/4] ddg kernel --- projects/CuLagrange/CMakeLists.txt | 2 +- .../geometry/BiharmonicBoundedWeight.cu | 2 +- .../CuLagrange/geometry/SolveLaplacian.cu | 2 +- ...e_matrix.hpp => differential_geometry.hpp} | 33 +++++++++++++++++++ 4 files changed, 36 insertions(+), 3 deletions(-) rename projects/CuLagrange/geometry/kernel/{laplace_matrix.hpp => differential_geometry.hpp} (93%) diff --git a/projects/CuLagrange/CMakeLists.txt b/projects/CuLagrange/CMakeLists.txt index 4c79dbf03d..829fe30410 100644 --- a/projects/CuLagrange/CMakeLists.txt +++ b/projects/CuLagrange/CMakeLists.txt @@ -106,7 +106,7 @@ target_sources(zeno PRIVATE geometry/file_parser/vtk_types.hpp geometry/linear_system/mfcg.hpp geometry/linear_system/active_set.hpp - geometry/kernel/laplace_matrix.hpp + geometry/kernel/differential_geometry.hpp geometry/kernel/gradient_field.hpp geometry/kernel/bary_centric_weights.hpp geometry/kernel/compute_characteristic_length.hpp diff --git a/projects/CuLagrange/geometry/BiharmonicBoundedWeight.cu b/projects/CuLagrange/geometry/BiharmonicBoundedWeight.cu index 47d30a4a98..56e042b569 100644 --- a/projects/CuLagrange/geometry/BiharmonicBoundedWeight.cu +++ b/projects/CuLagrange/geometry/BiharmonicBoundedWeight.cu @@ -13,7 +13,7 @@ #include #include -#include "kernel/laplace_matrix.hpp" +#include "kernel/differential_geometry.hpp" #include "linear_system/active_set.hpp" namespace zeno { diff --git a/projects/CuLagrange/geometry/SolveLaplacian.cu b/projects/CuLagrange/geometry/SolveLaplacian.cu index 9125732a75..be790a2b54 100644 --- a/projects/CuLagrange/geometry/SolveLaplacian.cu +++ b/projects/CuLagrange/geometry/SolveLaplacian.cu @@ -13,7 +13,7 @@ #include #include -#include "kernel/laplace_matrix.hpp" +#include "kernel/differential_geometry.hpp" #include "linear_system/mfcg.hpp" namespace zeno { diff --git a/projects/CuLagrange/geometry/kernel/laplace_matrix.hpp b/projects/CuLagrange/geometry/kernel/differential_geometry.hpp similarity index 93% rename from projects/CuLagrange/geometry/kernel/laplace_matrix.hpp rename to projects/CuLagrange/geometry/kernel/differential_geometry.hpp index ebdb9f9cc2..17047467fa 100644 --- a/projects/CuLagrange/geometry/kernel/laplace_matrix.hpp +++ b/projects/CuLagrange/geometry/kernel/differential_geometry.hpp @@ -68,6 +68,39 @@ namespace zeno { theta(5) = zs::acos(cos_theta(5)); } + // template + // constexpr void compute_enclosed_volume_of_triangulate_mesh() { + + // } + + // template + // constexpr void compute_area_vector_of_triangulate_mesh() { + + // } + + // template + // constexpr void compute_total_area_of_surrface_mesh() { + + // } + + // template + // constexpr void compute_mean_curvature_normal_of_surface_mesh() { + + // } + + // template + // constexpr void compute_total_mean_curvature_of_surface_mesh() { + + // } + + // tempalte + // constexpr void compute_gauss_curvature_normal_of_surface_mesh() { + + // } + + + + template constexpr void compute_cotmatrix(Pol &pol,const ETileVec &eles, const VTileVec &verts, const zs::SmallString& xTag, From 0b2add70680c04cca81bd89d6fcb2aaa0dd23d84 Mon Sep 17 00:00:00 2001 From: ShuliangLu Date: Thu, 13 Jul 2023 17:52:30 +0800 Subject: [PATCH 2/4] graph coloring --- projects/CuLagrange/geometry/Topology.cu | 57 +++++ .../CuLagrange/geometry/kernel/topology.hpp | 208 +++++++++--------- 2 files changed, 160 insertions(+), 105 deletions(-) diff --git a/projects/CuLagrange/geometry/Topology.cu b/projects/CuLagrange/geometry/Topology.cu index 259cb39fc8..9e0bbe274f 100644 --- a/projects/CuLagrange/geometry/Topology.cu +++ b/projects/CuLagrange/geometry/Topology.cu @@ -687,4 +687,61 @@ ZENDEFNODE(BuildSurfFacetTetraNeighboring, {{{"zsparticles"}}, }, {"ZSGeometry"}}); +struct DoTopogicalColoring : zeno::INode { + virtual void apply() override { + using namespace zs; + using vec2i = zs::vec; + using vec3i = zs::vec; + using vec4i = zs::vec; + + auto cudaPol = zs::cuda_exec(); + constexpr auto space = zs::execspace_e::cuda; + + auto zsparticles = get_input("zsparticles"); + const auto& verts = zsparticles->getParticles(); + auto& elms = zsparticles->getQuadraturePoints(); + // auto& tris = (*zsparticles)[ZenoParticles::s_surfTriTag]; + // const auto& tets = zsparticles->getQuadraturePoints(); + auto cdim = elms.getPropertySize("inds"); + auto color_tag = get_param("colorTag"); + + if(!elms.hasProperty(color_tag)) + elms.append_channels(cudaPol,{{color_tag,1}}); + + zs::Vector topos{elms.get_allocator(),elms.size()}; + cudaPol(zs::range(elms.size()),[ + elms = proxy({},elms), + cdim, + topos = proxy(topos)] ZS_LAMBDA(int ti) mutable { + topos[ti] = vec4i::uniform(-1); + for(int i = 0;i != cdim;++i) { + topos[ti][i] = zs::reinterpret_bits(elms("inds",i,ti)); + } + }); + + zs::Vector colors{elms.get_allocator(),elms.size()}; + std::cout << "do topological coloring" << std::endl; + topological_coloring(cudaPol,topos,colors); + std::cout << "finish topological coloring" << std::endl; + + cudaPol(zs::range(elms.size()),[ + elms = proxy({},elms), + color_tag = zs::SmallString(color_tag), + colors = proxy(colors)] ZS_LAMBDA(int ei) mutable { + elms(color_tag,ei) = colors[ei]; + }); + + set_output("zsparticles",zsparticles); + } +}; + +ZENDEFNODE(DoTopogicalColoring, {{{"zsparticles"}}, + { + {"zsparticles"} + }, + { + {"string","colorTag","colorTag"} + }, + {"ZSGeometry"}}); + }; \ No newline at end of file diff --git a/projects/CuLagrange/geometry/kernel/topology.hpp b/projects/CuLagrange/geometry/kernel/topology.hpp index f37041a409..c26c48a28f 100644 --- a/projects/CuLagrange/geometry/kernel/topology.hpp +++ b/projects/CuLagrange/geometry/kernel/topology.hpp @@ -14,6 +14,7 @@ #include "zensim/graph/ConnectedComponents.hpp" #include "zensim/container/Bht.hpp" +#include "zensim/graph/Coloring.hpp" #include "compute_characteristic_length.hpp" @@ -1056,176 +1057,173 @@ namespace zeno { template= 2), (VecTI::etent <= 4)> = 0*/> void topological_incidence_matrix(Pol& pol, - int nm_points, + // size_t nm_points, const zs::Vector& topos, - zs::SparseMatrix& spmat) { + zs::SparseMatrix& spmat) { using namespace zs; using ICoord = zs::vec; constexpr auto CDIM = VecTI::extent; constexpr auto space = Pol::exec_tag::value; constexpr auto execTag = wrapv{}; + zs::Vector max_pi_vec{topos.get_allocator(),1}; + max_pi_vec.setVal(0); + pol(zs::range(topos),[max_pi_vec = proxy(max_pi_vec),execTag,CDIM] ZS_LAMBDA(const auto& topo) { + for(int i = 0;i != CDIM;++i) + if(topo[i] >= 0) + atomic_max(execTag,&max_pi_vec[0],(int)topo[i]); + }); + auto nm_points = max_pi_vec.getVal(0) + 1; - // auto cudaPol = cuda_exec(); - - zs::Vector exclusive_offsets{topos.get_allocator(),nm_points + 1}; + zs::Vector exclusive_offsets{topos.get_allocator(),(size_t)(nm_points)}; zs::Vector p2ts{topos.get_allocator(),0}; + zs::Vector max_tp_incidences{topos.get_allocator(),1}; + zs::Vector cnts{topos.get_allocator(),(size_t)nm_points}; - { - zs::Vector cnts{topos.get_allocator(),nm_points}; + { zs::Vector tab_buffer{topos.get_allocator(), topos.size() * CDIM}; bht tab{topos.get_allocator(), topos.size() * CDIM}; tab.reset(pol, true); - cnts.reset(0); + // cnts.reset(0); + pol(zs::range(cnts),[] ZS_LAMBDA(auto& cnt) {cnt = 0;}); pol(zs::range(topos.size()),[ topos = proxy(topos), tab = proxy(tab), tab_buffer = proxy(tab_buffer), cnts = proxy(cnts)] ZS_LAMBDA(int ti) mutable { - for(int i = 0;i != CDIM;++i) + for(int i = 0;i != CDIM;++i) { if(topos[ti][i] < 0) break; else{ - int local_offset = atomic_add(execTag,&cnts[topos[ti][i]], (int)1); - if(auto id = tab.insert(ICoord{topos[ti][i],local_offset}); id != bht::sentinel_v){ + auto local_offset = atomic_add(execTag,&cnts[topos[ti][i]], (int)1); + if(auto id = tab.insert(ICoord{topos[ti][i],(int)local_offset}); id != bht::sentinel_v){ tab_buffer[id] = ti; } } + } }); + // std::cout << "finish computing tab_buffer" << std::endl; + // pol(zs::range(cnts.size()),[cnts = proxy(cnts)] ZS_LAMBDA(int pi) mutable {printf("cnts[%d] = %d\n",pi,cnts[pi]);}); + pol(zs::range(exclusive_offsets),[] ZS_LAMBDA(auto& eoffset) {eoffset = 0;}); + exclusive_scan(pol,std::begin(cnts),std::end(cnts),std::begin(exclusive_offsets)); - auto nmPTEntries = exclusive_offsets.getVal(nm_points); + // pol(zs::range(exclusive_offsets.size()),[exclusive_offsets = proxy(exclusive_offsets)] ZS_LAMBDA(int pi) mutable {printf("eooffset[%d] = %d\n",pi,exclusive_offsets[pi]);}); + auto nmPTEntries = exclusive_offsets.getVal(nm_points - 1) + cnts.getVal(nm_points - 1); + // std::cout << "nmPTEntries " << nmPTEntries << std::endl; p2ts.resize(nmPTEntries); + max_tp_incidences.setVal(0); pol(zs::range(nm_points),[ topos = proxy(topos), tab = proxy(tab), cnts = proxy(cnts), + execTag, + max_tp_incidences = proxy(max_tp_incidences), p2ts = proxy(p2ts), tab_buffer = proxy(tab_buffer), exclusive_offsets = proxy(exclusive_offsets)] ZS_LAMBDA(int pi) mutable { auto pt_count = cnts[pi]; + atomic_max(execTag,&max_tp_incidences[0],pt_count); auto ex_offset = exclusive_offsets[pi]; for(int i = 0;i != pt_count;++i) if(auto id = tab.query(ICoord{pi,i}); id != bht::sentinel_v) { auto ti = tab_buffer[id]; p2ts[ex_offset + i] = ti; + // printf("p[%d] -> t[%d]\n",pi,ti); } }); - } - zs::Vector is{topos.get_allocator(),topos.size()}; - zs::Vector js{topos.get_allocator(),topos.size()}; - pol(enumerate(is, js), [] ZS_LAMBDA(int no, int &i, int &j) mutable { i = j = no; }); - auto reserveStorage = [&is, &js](std::size_t n) { - auto size = is.size(); - is.resize(size + n); - js.resize(size + n); - return size; - }; - // auto tets_entry_offset = reserveStorage(p2ts.size()); - - { - bool success = false; - zs::Vector cnts{topos.get_allocator(),topos.size()}; + // std::cout << "finish computing p2ts" << std::endl; + } - // the buffer size might need to be resized - bht tab{topos.get_allocator(),topos.size() * CDIM * 2}; - zs::Vector tab_buffer{topos.get_allocator(), topos.size() * CDIM * 2}; - cnts.reset(0); - pol(range(topos.size()),[ - topos = proxy(topos), - tab = proxy(tab), - tab_buffer = proxy(tab_buffer), - p2ts = proxy(p2ts), - cnts = proxy(cnts), - execTag, - CDIM, - exclusive_offsets = proxy(exclusive_offsets)] ZS_LAMBDA(int ti) mutable { - auto topo = topos[ti]; - for(int i = 0;i != CDIM;++i){ - auto vi = topo[i]; - if(vi < 0) - return; - auto ex_offset = exclusive_offsets[vi]; - auto nm_nts = exclusive_offsets[vi + 1] - exclusive_offsets[vi]; - for(int j = 0;j != nm_nts;++j) { - auto nti = p2ts[ex_offset + j]; - if(nti > ti) - continue; - if(auto id = tab.insert(ICoord{ti,atomic_add(execTag,&cnts[ti],(int)1)}); id != bht::sentinel_v){ - tab_buffer[id] = nti; - } - } + bht tij_tab{topos.get_allocator(), topos.size() * max_tp_incidences.getVal(0) * CDIM}; + tij_tab.reset(pol,true); + + pol(range(topos.size()),[ + topos = proxy(topos), + p2ts = proxy(p2ts), + tij_tab = proxy(tij_tab), + execTag, + CDIM, + cnts = proxy(cnts), + exclusive_offsets = proxy(exclusive_offsets)] ZS_LAMBDA(int ti) mutable { + auto topo = topos[ti]; + for(int i = 0;i != CDIM;++i){ + auto vi = topo[i]; + if(vi < 0) + return; + auto ex_offset = exclusive_offsets[vi]; + auto nm_nts = cnts[vi]; + for(int j = 0;j != nm_nts;++j) { + auto nti = p2ts[ex_offset + j]; + if(nti < ti) + continue; + tij_tab.insert(ICoord{ti,nti}); } - }); - exclusive_offsets.resize(topos.size() + 1); - exclusive_scan(pol,std::begin(cnts),std::end(cnts),std::begin(exclusive_offsets)); - int nm_topo_incidences = exclusive_offsets.getVal(topos.size()); - auto topo_conn_entry_offset = reserveStorage(nm_topo_incidences); + } + }); - pol(range(topos.size()),[ - topos = proxy(topos), - topo_conn_entry_offset = topo_conn_entry_offset, - exclusive_offsets = proxy(exclusive_offsets), - tab = proxy(tab), - is = proxy(is), - js = proxy(js), - tab_buffer = proxy(tab_buffer)] ZS_LAMBDA(int ti) mutable { - auto ex_offset = exclusive_offsets[ti]; - auto nm_ntopos = exclusive_offsets[ti + 1] - exclusive_offsets[ti]; - for(int i = 0;i != nm_ntopos;++i) { - if(auto id = tab.insert(ICoord{ti,i}); id != bht::sentinel_v){ - auto nti = tab_buffer[id]; - is[ex_offset + i] = ti; - js[ex_offset + i] = nti; - } - } - }); + // std::cout << "finish computing tij_tab" << std::endl; + zs::Vector is{topos.get_allocator(),tij_tab.size()}; + zs::Vector js{topos.get_allocator(),tij_tab.size()}; + pol(zip(zs::range(tij_tab.size()),zs::range(tij_tab._activeKeys)),[is = proxy(is),js = proxy(js)] ZS_LAMBDA(auto idx,const auto& pair) { + is[idx] = pair[0];js[idx] = pair[1]; + // printf("pair[%d] : %d %d\n",idx,pair[0],pair[1]); + }); - } + // pol(zs::range(is.size()),[is = proxy(is),js = proxy(js)] ZS_LAMBDA(int i) mutable {printf("ijs[%d] : %d %d\n",i,is[i],js[i]);}); + // std::cout << "topos.size() = " << topos.size() << std::endl; - spmat = zs::SparseMatrix{topos.get_allocator(),(int)topos.size(),(int)topos.size()}; - spmat.build(pol,(int)nm_points,(int)topos.size(),zs::range(is),zs::range(js)/*,zs::range(rs)*/,zs::false_c); - spmat.localOrdering(pol,zs::false_c); + // spmat = zs::SparseMatrix{topos.get_allocator(),(int)topos.size(),(int)topos.size()}; + spmat.build(pol,(int)topos.size(),(int)topos.size(),zs::range(is),zs::range(js)/*,zs::range(rs)*/,zs::true_c); + // spmat.localOrdering(pol,zs::false_c); spmat._vals.resize(spmat.nnz()); - spmat._vals.reset((int)1); + pol(spmat._vals, []ZS_LAMBDA(u32 &v) { v = 1; }); + // std::cout << "done connectivity graph build" << std::endl; + + // spmat._vals.reset((int)1); } template void topological_coloring(Pol& pol, - int nm_points, + // int nm_points, const zs::Vector& topo, - zs::Vector& coloring) { + zs::Vector& colors) { using namespace zs; constexpr auto space = Pol::exec_tag::value; - coloring.resize(topo.size()); - zs::SparseMatrix pt_incidence{}; - topo_nodal_incidence_matrix(pol,nm_points,topo,pt_incidence); + colors.resize(topo.size()); + zs::SparseMatrix topo_incidence_matrix{topo.get_allocator(),topo.size(),topo.size()}; + // std::cout << "compute incidence matrix " << std::endl; + topological_incidence_matrix(pol,topo,topo_incidence_matrix); + // std::cout << "finish compute incidence matrix " << std::endl; - union_find(pol,pt_incidence,range(coloring)); - zs::bcht, 16> vtab{coloring.get_allocator(),coloring.size()}; - pol(range(coloring.size()),[ - vtab = proxy(vtab), - coloring = proxy(coloring)] ZS_LAMBDA(int vi) mutable { - auto fa = coloring[vi]; - while(fa != coloring[fa]) - fa = coloring[fa]; - coloring[vi] = fa; - vtab.insert(fa); - }); + auto ompPol = omp_exec(); + constexpr auto omp_space = execspace_e::openmp; + zs::Vector weights(/*topo.get_allocator(),*/topo.size()); + { + bht tab{weights.get_allocator(),topo.size() * 2}; + tab.reset(ompPol, true); + ompPol(enumerate(weights), [tab1 = proxy(tab)] (int seed, u32 &w) mutable { + using tab_t = RM_CVREF_T(tab); + std::mt19937 rng; + rng.seed(seed); + u32 v = rng() % (u32)4294967291u; + // prevent weight duplications + while (tab1.insert(v) != tab_t::sentinel_v) + v = rng() % (u32)4294967291u; + w = v; + }); + } + weights = weights.clone(colors.memoryLocation()); - pol(range(coloring.size()),[ - coloring = proxy(coloring),vtab = proxy(vtab)] ZS_LAMBDA(int vi) mutable { - auto ancestor = coloring[vi]; - auto setNo = vtab.query(ancestor); - coloring[vi] = setNo; - }); + auto iterRef = maximum_independent_sets(pol, topo_incidence_matrix, weights, colors); + std::cout << "nm_colors : " << iterRef << std::endl; } template From 0be8c890ed22e71fbd645950b48dd459fad56f0e Mon Sep 17 00:00:00 2001 From: ShuliangLu Date: Thu, 13 Jul 2023 20:04:36 +0800 Subject: [PATCH 3/4] update previous fem pipeline --- projects/CuLagrange/CMakeLists.txt | 1 + .../CuLagrange/fem/FleshDynamicStepping.cu | 652 ++++++----------- .../collision_energy/evaluate_collision.hpp | 516 +++++++++++--- projects/CuLagrange/geometry/BasicGeoNodes.cu | 85 +++ projects/CuLagrange/geometry/CollisionVis.cu | 355 ++++----- projects/CuLagrange/geometry/SurfaceBinder.cu | 674 ++++++++++-------- projects/CuLagrange/geometry/VectorField.cu | 56 +- .../geometry/file_parser/read_vtk_mesh.hpp | 40 +- .../geometry/kernel/intersection.hpp | 147 ++-- 9 files changed, 1415 insertions(+), 1111 deletions(-) create mode 100755 projects/CuLagrange/geometry/BasicGeoNodes.cu diff --git a/projects/CuLagrange/CMakeLists.txt b/projects/CuLagrange/CMakeLists.txt index 829fe30410..f539068361 100644 --- a/projects/CuLagrange/CMakeLists.txt +++ b/projects/CuLagrange/CMakeLists.txt @@ -127,6 +127,7 @@ target_sources(zeno PRIVATE # geometry/VectorField.cu geometry/CollisionVis.cu # CHECK THIS # geometry/CalcSurfaceArea.cu + geometry/BasicGeoNodes.cu geometry/MarkSurfaceTag.cu geometry/SurfaceBinder.cu geometry/MarkInversion.cu diff --git a/projects/CuLagrange/fem/FleshDynamicStepping.cu b/projects/CuLagrange/fem/FleshDynamicStepping.cu index b3b0fe2c3f..df1f50a378 100644 --- a/projects/CuLagrange/fem/FleshDynamicStepping.cu +++ b/projects/CuLagrange/fem/FleshDynamicStepping.cu @@ -133,6 +133,7 @@ struct FleshDynamicStepping : INode { using tiles_t = typename ZenoParticles::particles_t; using vec2 = zs::vec; using vec3 = zs::vec; + using vec12 = zs::vec; using mat3 = zs::vec; using mat9 = zs::vec; using mat12 = zs::vec; @@ -291,11 +292,6 @@ struct FleshDynamicStepping : INode { ++nm_binders; } - // printf("binder_ids[%d] : %d : %d %d %d %d\n",ti,nm_binders, - // reinterpret_bits(tris(binderTag,0,ti)), - // reinterpret_bits(tris(binderTag,1,ti)), - // reinterpret_bits(tris(binderTag,2,ti)), - // reinterpret_bits(tris(binderTag,3,ti))); if(nm_binders == 0) return; @@ -351,26 +347,22 @@ struct FleshDynamicStepping : INode { cp[0] = kverts.pack(dim_c<3>,"x",idx); auto kstiffness = (T)1.0; if(kverts.hasProperty("binderStiffness")) - kstiffness = kverts("binderStiffness",idx); + kstiffness *= kverts("binderStiffness",idx); + auto average_vert_stiffness = (T)1.0; + if(verts.hasProperty("binderStiffness")){ + average_vert_stiffness = (T)0.0; + for(int j = 0;j != 3;++j) + average_vert_stiffness += verts("binderStiffness",tri[j]); + } + kstiffness *= average_vert_stiffness; auto alpha = binderStiffness * binder_weakness_param * kstiffness; auto beta = (T)1.0/(T)nm_binders; auto cgrad = -alpha * beta * VERTEX_FACE_SQRT_COLLISION::gradient(cp,mu,lam,ceps,from_inside); auto cH = alpha * beta * VERTEX_FACE_SQRT_COLLISION::hessian(cp,mu,lam,ceps,from_inside); - // printf("cgrad : %f cH %f params: %f %f %f %f\n",cgrad.norm(),cH.norm(), - // (float)kstiffness, - // (float)binderStiffness, - // (float)binder_weakness_param, - // (float)alpha); - - // if(isnan(cH.norm())) { - // printf("nan CH detected at Binder : %d from inside %d and ceps = \n",ti,from_inside,(float)ceps); - // printf("cp : \n%f %f %f\n%f %f %f\n%f %f %f\n%f %f %f\n", - // (float)cp[0][0],(float)cp[0][1],(float)cp[0][2], - // (float)cp[1][0],(float)cp[1][1],(float)cp[1][2], - // (float)cp[2][0],(float)cp[2][1],(float)cp[2][2], - // (float)cp[3][0],(float)cp[3][1],(float)cp[3][2]); - // } + if(isnan(cH.norm()) || isnan(cgrad.norm())) { + printf("nan cH and cgrad detected at [%d] [%d] : %f %f\n",ti,idx,(float)cH.norm(),(float)cgrad.norm()); + } for(int i = 3;i != 12;++i){ int d0 = i % 3; @@ -392,113 +384,10 @@ struct FleshDynamicStepping : INode { }); } - // template - // void computeCollisionGradientAndHessian(zs::CudaExecutionPolicy& cudaPol,const Model& model, - // dtiles_t& vtemp, - // dtiles_t& etemp, - // dtiles_t& sttemp, - // dtiles_t& setemp, - // // dtiles_t& ee_buffer, - // dtiles_t& fp_buffer, - // dtiles_t& kverts, - // dtiles_t& kc_buffer, - // dtiles_t& gh_buffer, - // T kd_theta = (T)0.0, - // bool explicit_collision = false, - // bool neglect_inverted = true) { - // using namespace zs; - // constexpr auto space = execspace_e::cuda; - - // int offset = eles.size(); - - // T lambda = model.lam; - // T mu = model.mu; - - // // auto stBvh = bvh_t{}; - // // auto bvs = retrieve_bounding_volumes(cudaPol,vtemp,tris,wrapv<3>{},(T)0.0,"xn"); - // // stBvh.build(cudaPol,bvs); - // // auto avgl = compute_average_edge_length(cudaPol,vtemp,"xn",tris); - // // auto bvh_thickness = 5 * avgl; - // // if(!calculate_facet_normal(cudaPol,vtemp,"xn",tris,sttemp,"nrm")){ - // // throw std::runtime_error("fail updating facet normal"); - // // } - // // if(!COLLISION_UTILS::calculate_cell_bisector_normal(cudaPol, - // // vtemp,"xn", - // // lines, - // // tris, - // // sttemp,"nrm", - // // setemp,"nrm")){ - // // throw std::runtime_error("fail calculate cell bisector normal"); - // // } - - - // COLLISION_UTILS::do_facet_point_collision_detection(cudaPol, - // vtemp,"xn", - // points, - // lines, - // tris, - // sttemp, - // setemp, - // fp_buffer, - // in_collisionEps,out_collisionEps); - - // COLLISION_UTILS::evaluate_fp_collision_grad_and_hessian(cudaPol, - // vtemp,"xn","vn",dt, - // fp_buffer, - // gh_buffer,offset, - // in_collisionEps,out_collisionEps, - // (T)collisionStiffness, - // (T)mu,(T)lambda,(T)kd_theta); - - - - // // COLLISION_UTILS::do_kinematic_point_collision_detection(cudaPol, - // // vtemp,"xn", - // // points, - // // lines, - // // tris, - // // setemp, - // // sttemp, - // // kverts, - // // kc_buffer, - // // (T)kine_in_collisionEps,(T)kine_out_collisionEps,false); - - // // offset = 0; - - // // COLLISION_UTILS::evaluate_kinematic_fp_collision_grad_and_hessian(cudaPol, - // // eles, - // // vtemp,"xn","vn",dt, - // // tris, - // // kverts, - // // kc_buffer, - // // gh_buffer,offset, - // // (T)kine_in_collisionEps,(T)kine_out_collisionEps, - // // (T)kineCollisionStiffness, - // // (T)mu,(T)lambda,(T)kd_theta); - - - // // adding collision damping on self collision - // // int offset = eles.size() + b_verts.size(); - // // cudaPol(zs::range(fp_buffer.size() + kc_buffer.size()), - // // [vtemp = proxy({},vtemp), - // // gh_buffer = proxy({},gh_buffer),offset,kd_theta] ZS_LAMBDA(int ci) mutable { - // // auto inds = gh_buffer.pack(dim_c<4>,"inds",ci).reinterpret_bits(int_c); - // // for(int i = 0;i != 4;++i) - // // if(inds[i] < 0) - // // return; - // // vec3 vs[4] = {}; - // // for(int i = 0;i = 4;++i) - // // vs[i] = vtemp.pack(dim_c<3>,"vn",inds[i]); - // // auto H = gh_buffer.pack(dim_c<12*12>,"H",ci); - // // gh_buffer.tuple(dim_c<12*12>,"H",ci) = H; - // // }); - - - // } void computePlaneConstraintGradientAndHessian2(zs::CudaExecutionPolicy& cudaPol, const dtiles_t& vtemp, - const dtiles_t& sttemp, + dtiles_t& sttemp, const dtiles_t& kverts, const dtiles_t& ktris, const std::string& planeConsBaryTag, @@ -520,7 +409,7 @@ struct FleshDynamicStepping : INode { plane_constraint_stiffness = plane_constraint_stiffness, use_sticky_condition = use_sticky_condition, nodal_gh_buffer = proxy({},nodal_gh_buffer)] ZS_LAMBDA(int vi) mutable { - return; + // return; auto idx = reinterpret_bits(verts(planeConsIDTag,vi)); if(idx < 0) return; @@ -530,7 +419,6 @@ struct FleshDynamicStepping : INode { if(is_inverted_vert) return; - #if 1 auto plane_root = kverts.pack(dim_c<3>,"x",ktri[0]); auto plane_nrm = ktris.pack(dim_c<3>,"nrm",idx); @@ -581,79 +469,135 @@ struct FleshDynamicStepping : INode { }); - cudaPol(zs::range(tris.size()),[ - vtemp = proxy({},vtemp), - sttemp = proxy({},sttemp), - verts = proxy({},verts), - tris = proxy({},tris), - kverts = proxy({},kverts), - ktris = proxy({},ktris), - cnorm = cnorm, - planeConsIDTag = zs::SmallString(planeConsIDTag), - kine_out_collisionEps = kine_out_collisionEps, - kine_in_collisionEps = kine_in_collisionEps, - plane_constraint_stiffness = plane_constraint_stiffness, - use_sticky_condition = use_sticky_condition, - tris_gh_buffer = proxy({},tris_gh_buffer)] ZS_LAMBDA(int ti) mutable { - auto kp_idx = reinterpret_bits(tris(planeConsIDTag,ti)); - if(kp_idx < 0) - return; - auto kp = kverts.pack(dim_c<3>,"x",kp_idx); - auto tri = tris.pack(dim_c<3>,"inds",ti).reinterpret_bits(int_c); - for(int i = 0;i != 3;++i){ - auto is_inverted_vert = vtemp("is_inverted",tri[i]) > (T)0.5; - if(is_inverted_vert) - return; - } + // cudaPol(zs::range(tris.size()),[ + // vtemp = proxy({},vtemp), + // sttemp = proxy({},sttemp), + // verts = proxy({},verts), + // tris = proxy({},tris), + // kverts = proxy({},kverts), + // ktris = proxy({},ktris), + // cnorm = cnorm, + // planeConsIDTag = zs::SmallString(planeConsIDTag), + // kine_out_collisionEps = kine_out_collisionEps, + // kine_in_collisionEps = kine_in_collisionEps, + // plane_constraint_stiffness = plane_constraint_stiffness, + // use_sticky_condition = use_sticky_condition, + // tris_gh_buffer = proxy({},tris_gh_buffer)] ZS_LAMBDA(int ti) mutable { + // auto kp_idx = reinterpret_bits(tris(planeConsIDTag,ti)); + // if(kp_idx < 0) + // return; + // auto kp = kverts.pack(dim_c<3>,"x",kp_idx); + // auto tri = tris.pack(dim_c<3>,"inds",ti).reinterpret_bits(int_c); + // for(int i = 0;i != 3;++i){ + // auto is_inverted_vert = vtemp("is_inverted",tri[i]) > (T)0.5; + // if(is_inverted_vert) + // return; + // } - // auto tnrm = sttemp.pack(dim_c<3>,"nrm",ti); + // // auto tnrm = sttemp.pack(dim_c<3>,"nrm",ti); - auto mu = verts("mu",tri[0]); - auto lam = verts("lam",tri[0]); + // auto mu = verts("mu",tri[0]); + // auto lam = verts("lam",tri[0]); - auto eps = kine_out_collisionEps; - vec3 vs[4] = {}; - vs[0] = kp; - for(int i = 0;i != 3;++i) - vs[i + 1] = vtemp.pack(dim_c<3>,"xn",tri[i]); + // auto eps = kine_out_collisionEps; + // vec3 vs[4] = {}; + // vs[0] = kp; + // for(int i = 0;i != 3;++i) + // vs[i + 1] = vtemp.pack(dim_c<3>,"xn",tri[i]); - vec3 e[3] = {}; - e[0] = vs[3] - vs[2]; - e[1] = vs[0] - vs[2]; - e[2] = vs[1] - vs[2]; - - auto n = e[2].cross(e[0]); - // if(n.norm() < 1e-4) - // return; - n = n/(n.norm() + 1e-6); - - T springLength = e[1].dot(n) - eps; - auto gvf = zs::vec::zeros(); - if(springLength < (T)0 || use_sticky_condition){ - auto gvf_v12 = COLLISION_UTILS::springLengthGradient(vs,e,n); - if(isnan(gvf_v12.norm())) - printf("nan gvf detected at %d %f %f\n",ti,gvf_v12.norm(),n.norm()); - for(int i = 0;i != 9;++i) - gvf[i] = gvf_v12[i + 3]; - } - cnorm = (T)1.0; - auto stiffness = plane_constraint_stiffness * cnorm; - // stiffness = (T)0; - auto g = -stiffness * (T)2.0 * mu * springLength * gvf; - auto H = stiffness * (T)2.0 * mu * zs::dyadic_prod(gvf, gvf); + // vec3 e[3] = {}; + // e[0] = vs[3] - vs[2]; + // e[1] = vs[0] - vs[2]; + // e[2] = vs[1] - vs[2]; + + // // auto seg = vs[0] - vs[2]; + // // auto n = e[2].cross(e[0]); + // auto n = LSL_GEO::facet_normal(vs[1],vs[2],vs[3]); + // // if(n.norm() < 1e-4) + // // return; + // // n = n/(n.norm() + 1e-6); + + // T springLength = e[1].dot(n) + eps; + // auto gvf = zs::vec::zeros(); + // if(springLength > (T)0 || use_sticky_condition){ + // auto gvf_v12 = COLLISION_UTILS::springLengthGradient(vs,e,n); + // if(isnan(gvf_v12.norm())) + // printf("nan gvf detected at %d %f %f\n",ti,gvf_v12.norm(),n.norm()); + // for(int i = 0;i != 9;++i) + // gvf[i] = gvf_v12[i + 3]; + // } + // cnorm = (T)1.0; + // auto stiffness = plane_constraint_stiffness * cnorm; + // // stiffness = (T)0; + // auto g = stiffness * (T)2.0 * mu * springLength * gvf; + // auto H = stiffness * (T)2.0 * mu * zs::dyadic_prod(gvf, gvf); - // if(springLength < (T)0) { - // auto springLengthH_M12 = COLLISION_UTILS::springLengthHessian(vs,e,n); - // auto springLengthH_M9 = mat9::zeros(); - // for(int r = 0;r != 9;++r) - // for(int c = 0;c != 9;++c) - // springLengthH_M9(r,c) = springLengthH_M12(r + 3,c+ 3); - // H += springLength * springLengthH_M9 * (T)2.0 * stiffness * mu; - // make_pd(H); - // } + // // if(springLength < (T)0) { + // // auto springLengthH_M12 = COLLISION_UTILS::springLengthHessian(vs,e,n); + // // auto springLengthH_M9 = mat9::zeros(); + // // for(int r = 0;r != 9;++r) + // // for(int c = 0;c != 9;++c) + // // springLengthH_M9(r,c) = springLengthH_M12(r + 3,c+ 3); + // // H += springLength * springLengthH_M9 * (T)2.0 * stiffness * mu; + // // make_pd(H); + // // } - tris_gh_buffer.tuple(dim_c<9>,"grad",ti) = g; - tris_gh_buffer.tuple(dim_c<9,9>,"H",ti) = H; + // tris_gh_buffer.tuple(dim_c<9>,"grad",ti) = g; + // tris_gh_buffer.tuple(dim_c<9,9>,"H",ti) = H; + // }); + + + cudaPol(zs::range(kverts.size()),[ + vtemp = proxy({},vtemp), + tris = proxy({},tris), + kverts = proxy({},kverts), + verts = proxy({},verts), + sttemp = proxy({},sttemp), + planeConsIDTag = zs::SmallString(planeConsIDTag), + kine_out_collisionEps = kine_out_collisionEps, + plane_constraint_stiffness = plane_constraint_stiffness, + use_sticky_condition = use_sticky_condition] ZS_LAMBDA(int kvi) mutable { + auto ti = reinterpret_bits(kverts(planeConsIDTag,kvi)); + if(ti < 0) + return; + auto kp = kverts.pack(dim_c<3>,"x",kvi); + auto tri = tris.pack(dim_c<3>,"inds",ti,int_c); + + auto plane_root = vtemp.pack(dim_c<3>,"xn",tri[0]); + vec3 tvs[3] = {}; + for(int i = 0;i != 3;++i) + tvs[i] = vtemp.pack(dim_c<3>,"xn",tri[i]); + auto plane_nrm = LSL_GEO::facet_normal(tvs[0],tvs[1],tvs[2]); + // auto plane_nrm = tris.pack(dim_c<3>,"nrm",ti); + + auto mu = verts("mu",tri[0]); + auto lam = verts("lam",tri[0]); + + auto eps = kine_out_collisionEps; + auto seg = kp - plane_root; + + auto fc = vec12::zeros(); + auto Hc = mat12::zeros(); + auto dist = seg.dot(plane_nrm) + eps; + if(dist > (T)0 || use_sticky_condition) { + // fc = -dist * mu * plane_constraint_stiffness * plane_nrm; + // Hc = mu * plane_constraint_stiffness * dyadic_prod(plane_nrm,plane_nrm); + // printf("detected kv2t pairs : %f %f\n",(T)fc.norm(),(T)Hc.norm()); + vec3 vs[4] = {}; + for(int i = 0;i != 3;++i) + vs[i + 1] = tvs[i]; + vs[0] = kp; + fc = -VERTEX_FACE_COLLISION::gradient(vs,mu,lam,-eps); + Hc = VERTEX_FACE_COLLISION::hessian(vs,mu,lam,-eps); + } + + + + for(int i = 0;i != 9;++i) + atomic_add(exec_cuda,&sttemp("grad",i,ti),fc[i + 3]); + for(int i = 0;i != 9;++i) + for(int j = 0;j != 9;++j) + atomic_add(exec_cuda,&sttemp("H",i * 9 + j,ti),Hc(i + 3,j + 3)); }); } @@ -707,56 +651,6 @@ struct FleshDynamicStepping : INode { nodal_gh_buffer.tuple(dim_c<3>,"grad",vi) = fc; nodal_gh_buffer.tuple(dim_c<3,3>,"H",vi) = Hc; }); - - // cudaPol(zs::range(tris.size()),[ - // verts = proxy({},verts), - // tris = proxy({},tris), - // vtemp = proxy({},vtemp), - // planeConsPosTag = zs::SmallString(planeConsPosTag), - // planeConsNrmTag = zs::SmallString(planeConsNrmTag), - // planeConsIDTag = zs::SmallString(planeConsIDTag), - // kine_out_collisionEps = kine_out_collisionEps, - // plane_constraint_stiffness = plane_constraint_stiffness, - // nodal_gh_buffer = proxy({},nodal_gh_buffer)] ZS_LAMBDA(int ti) mutable { - // auto idx = reinterpret_bits(tris(planeConsIDTag,ti)); - // if(idx < 0) - // return; - - // auto tri = tris.pack(dim_c<3>,"inds",ti).reinterpret_bits(int_c); - - // auto mu = verts("mu",tri[0]); - // auto lam = verts("lam",tri[0]); - - // auto eps = kine_out_collisionEps * 2.0; - // auto plane_nrm = tris.pack(dim_c<3>,planeConsNrmTag,ti); - // auto plane_root = tris.pack(dim_c<3>,planeConsPosTag,ti); - - // auto p = vec3::zeros(); - // for(int i = 0;i != 3;++i) - // p += vtemp.pack(dim_c<3>,"xn",tri[i])/(T)3.0; - // auto seg = p - plane_root; - - // auto fc = vec3::zeros(); - // auto Hc = mat3::zeros(); - // auto dist = seg.dot(plane_nrm) - eps; - // if(dist < (T)0){ - // fc = -dist * mu * plane_constraint_stiffness * plane_nrm; - // Hc = mu * plane_constraint_stiffness * dyadic_prod(plane_nrm,plane_nrm); - // } - - // // printf("apply plane constraint with force : %f %f\n",(float)dist,(float)fc.norm()); - // for(int i = 0;i != 3;++i) { - // auto vi = tri[i]; - // for(int d = 0;d != 3;++d) - // atomic_add(exec_cuda,&nodal_gh_buffer("grad",d,vi),fc[d]/(T)3.0); - // for(int r = 0;r != 3;++r) - // for(int c = 0;c != 3;++c) - // atomic_add(exec_cuda,&nodal_gh_buffer("H",r * 3 + c,vi),Hc(r,c)/(T)9.0); - // } - - // // nodal_gh_buffer.tuple(dim_c<3>,"grad",vi) = fc; - // // nodal_gh_buffer.tuple(dim_c<3,3>,"H",vi) = Hc; - // }); } template @@ -886,21 +780,6 @@ struct FleshDynamicStepping : INode { // add inertia hessian term auto H = dFdAct_dFdX.transpose() * Hq * dFdAct_dFdX * vole; - // if(isnan(H.norm())) { - // printf("nan CH detected at Elastic : %d %f %f %f %f\nFAct = \n%f %f %f\n%f %f %f\n%f %f %f\nF = \n%f %f %f\n%f %f %f\n%f %f %f\n",ei, - // (float)Hq.norm(), - // (float)dFdAct_dFdX.norm(), - // (float)P.norm(), - // (float)FAct.norm(), - // (float)FAct(0,0),(float)FAct(0,1),(float)FAct(0,2), - // (float)FAct(1,0),(float)FAct(1,1),(float)FAct(1,2), - // (float)FAct(2,0),(float)FAct(2,1),(float)FAct(2,2), - // (float)F(0,0),(float)F(0,1),(float)F(0,2), - // (float)F(1,0),(float)F(1,1),(float)F(1,2), - // (float)F(2,0),(float)F(2,1),(float)F(2,2) - // ); - // } - if(eles.hasProperty("Muscle_ID") && (int)eles("Muscle_ID",ei) >= 0) { auto fiber = eles.pack(dim_c<3>,"fiber",ei); if(zs::abs(fiber.norm() - 1.0) < 1e-3) { @@ -924,22 +803,6 @@ struct FleshDynamicStepping : INode { } } - - // adding rayleigh damping term - // vec3 v0[4] = {vtemp.pack(dim_c<3>,"vn", inds[0]), - // vtemp.pack(dim_c<3>,"vn", inds[1]), - // vtemp.pack(dim_c<3>,"vn", inds[2]), - // vtemp.pack(dim_c<3>,"vn", inds[3])}; - - // auto inertia = (T)1.0; - // if(eles.hasProperty("inertia")) - // inertia = eles("inertia",ei); - - // auto vel = COLLISION_UTILS::flatten(v0); - // auto m = eles("m",ei)/(T)4.0; - // auto C = kd_beta * H + kd_alpha * inertia * m * zs::vec::identity(); - // auto rdamping = C * vel; - gh_buffer.tuple(dim_c<12>,"grad",ei + offset) = gh_buffer.pack(dim_c<12>,"grad",ei + offset) + vf/* - rdamping*/; // gh_buffer.tuple(dim_c<12>,"grad",ei + offset) = gh_buffer.pack(dim_c<12>,"grad",ei + offset) - rdamping; // H += kd_beta*H/dt; @@ -952,14 +815,6 @@ struct FleshDynamicStepping : INode { auto nmEmbedVerts = b_verts.size(); - // TILEVEC_OPS::fill_range<4>(cudaPol,gh_buffer,"inds",zs::vec::uniform(-1).reinterpret_bits(float_c),eles.size() + offset,b_verts.size()); - // TILEVEC_OPS::fill_range<3>(cudaPol,gh_buffer,"grad",zs::vec::zeros(),eles.size() + offset,b_verts.size()); - // TILEVEC_OPS::fill_range<144>(cudaPol,gh_buffer,"H",zs::vec::zeros(),eles.size() + offset,b_verts.size()); - - // we should neglect the inverted element - // std::cout << "nmEmbedVerts : " << nmEmbedVerts << std::endl; - // std::cout << "bcwsize : " << b_bcws.size() << std::endl; - // return; cudaPol(zs::range(nmEmbedVerts), [ gh_buffer = proxy({},gh_buffer),model = model, bcws = proxy({},b_bcws),b_verts = proxy({},b_verts),vtemp = proxy({},vtemp),etemp = proxy({},etemp), @@ -973,23 +828,10 @@ struct FleshDynamicStepping : INode { if(b_verts.hasProperty("intersect")) if(b_verts("intersect",vi) > (T)0.5){ - printf("skip bverts[%d] due to intersection\n",vi); + // printf("skip bverts[%d] due to intersection\n",vi); return; } - // if(ei >= etemp.size()){ - // printf("ei too big for etemp\n"); - // return; - // } - // auto is_inverted = reinterpret_bits(etemp("is_inverted",ei)); - // if(is_inverted){ - // if(vi == 0) - // printf("inverted tet\n"); - // return; - // } - // auto FatID = eles("FatID",ei); - // if(FatID > 0) - // return; auto lambda = model.lam; auto mu = model.mu; @@ -1014,18 +856,6 @@ struct FleshDynamicStepping : INode { T stiffness = (2.0066 * mu + 1.0122 * lambda) * b_verts("strength",vi); - // zs::vec elm_grad{}; - // auto elm_H = zs::vec::zeros(); - - // if(vi == 0) { - // printf("stiff : %f dw : %f strength : %f cnorm : %f vol : %f bdw : %f\n", - // (float)stiffness, - // (float)bone_driven_weight, - // (float)bcws("strength",vi), - // (float)bcws("cnorm",vi), - // (float)eles("vol",ei), - // (float)eles("bdw",ei)); - // } auto alpha = stiffness * bone_driven_weight * bcws("strength",vi) * bcws("cnorm",vi) * eles("vol",ei) * eles("bdw",ei); @@ -1049,21 +879,9 @@ struct FleshDynamicStepping : INode { atomic_add(exec_cuda,&gh_buffer("H",(i*3 + d)*12 + j*3 + d,ei),beta); } } - - // for(int i = 0;i != 12;++i){ - // atomic_add(exec_cuda,&gh_buffer("grad",i,ei),elm_grad[i]); - // for(int j = 0;j != 12;++j) - // atomic_add(exec_cuda,&gh_buffer("H",i*12 + j,ei),elm_H(i,j)); - // } - // gh_buffer.tuple(dim_c<12>,"grad",vi + eles.size() + offset) = elm_grad; - // gh_buffer.tuple(dim_c<12*12>,"H",vi + eles.size() + offset) = elm_H; + }); - // cudaPol(zs::range(eles.size()), [gh_buffer = proxy({},gh_buffer)] ZS_LAMBDA (int ei) mutable { - // auto H = gh_buffer.template pack<12,12>("H",ei); - // make_pd(H); - // gh_buffer.template tuple<12*12>("H",ei) = H; - // }); } @@ -1499,6 +1317,11 @@ struct FleshDynamicStepping : INode { } } + auto planeConsPosTag = get_param("planeConsPosTag"); + auto planeConsNrmTag = get_param("planeConsNrmTag"); + auto planeConsIDTag = get_param("planeConsIDTag"); + auto planeConsBaryTag = get_param("planeConsBaryTag"); + // auto bone_driven_weight = (T)0.02; auto newton_res = get_input2("newton_res"); @@ -1513,14 +1336,9 @@ struct FleshDynamicStepping : INode { if(has_input("Acts")) { act_ = get_input("Acts")->getLiterial(); nm_acts = act_.size(); - - std::cout << "use Activation : "; - for(int i = 0;i != act_.size();++i) - std::cout << "[" << act_[i][0] << "\t" << act_[i][1] << "] "; - std::cout << std::endl; } - std::cout << "nmActs:" << nm_acts << std::endl; + // std::cout << "nmActs:" << nm_acts << std::endl; constexpr auto host_space = zs::execspace_e::openmp; auto ompExec = zs::omp_exec(); @@ -1590,6 +1408,7 @@ struct FleshDynamicStepping : INode { {"xp",3}, {"b_fail",1}, {"binderStiffness",1}, + {planeConsIDTag,1}, {"nrm",3}, {"area",1}},0,zs::memsrc_e::device,0); auto ktris = typename ZenoParticles::particles_t({ @@ -1597,28 +1416,28 @@ struct FleshDynamicStepping : INode { {"nrm",3}},0,zs::memsrc_e::device,0); - dtiles_t surf_tris_buffer{tris.get_allocator(),{ - {"inds",3}, - {"nrm",3}, - {"he_inds",1} - },tris.size()}; + // dtiles_t surf_tris_buffer{tris.get_allocator(),{ + // {"inds",3}, + // {"nrm",3}, + // {"he_inds",1} + // },tris.size()}; - dtiles_t surf_verts_buffer{points.get_allocator(),{ - {"inds",1}, - {"xn",3}, - {"is_corner",1}, - {"mustExclude",1} - },points.size()}; - TILEVEC_OPS::copy(cudaPol,points,"inds",surf_verts_buffer,"inds"); - TILEVEC_OPS::copy(cudaPol,tris,"inds",surf_tris_buffer,"inds"); - TILEVEC_OPS::copy(cudaPol,tris,"he_inds",surf_tris_buffer,"he_inds"); - reorder_topology(cudaPol,points,surf_tris_buffer); + // dtiles_t surf_verts_buffer{points.get_allocator(),{ + // {"inds",1}, + // {"xn",3}, + // {"is_loop_vertex",1}, + // {"mustExclude",1} + // },points.size()}; + // TILEVEC_OPS::copy(cudaPol,points,"inds",surf_verts_buffer,"inds"); + // TILEVEC_OPS::copy(cudaPol,tris,"inds",surf_tris_buffer,"inds"); + // TILEVEC_OPS::copy(cudaPol,tris,"he_inds",surf_tris_buffer,"he_inds"); + // reorder_topology(cudaPol,points,surf_tris_buffer); // zs::Vector nodal_colors{surf_verts_buffer.get_allocator(),surf_verts_buffer.size()}; dtiles_t gia_res{points.get_allocator(),{ {"ring_mask",1}, {"type_mask",1}, {"color_mask",1}, - {"is_corner",1} + {"is_loop_vertex",1} },points.size()}; dtiles_t tri_gia_res{tris.get_allocator(),{ {"ring_mask",1}, @@ -1637,6 +1456,11 @@ struct FleshDynamicStepping : INode { },tris.size() * 2}; + + bool use_plane_constraint = get_input2("use_plane_constraint"); + bool use_binder_constraint = get_input2("use_binder_constraint"); + bool use_kinematic_potential = get_input2("with_kinematic_potential"); + if(has_input("kinematic_boundary")){ auto kinematic_boundary = get_input("kinematic_boundary"); // if (kinematic_boundary.empty()) @@ -1668,6 +1492,11 @@ struct FleshDynamicStepping : INode { TILEVEC_OPS::copy(cudaPol,kb_verts,"binderStiffness",kverts,"binderStiffness"); else TILEVEC_OPS::fill(cudaPol,kverts,"binderStiffness",(T)1.0); + + if(kb_verts.hasProperty(planeConsIDTag)) + TILEVEC_OPS::copy(cudaPol,kb_verts,planeConsIDTag,kverts,planeConsIDTag); + else + TILEVEC_OPS::fill(cudaPol,kverts,planeConsIDTag,zs::reinterpret_bits((int)-1)); ktris.resize(kb_tris.size()); TILEVEC_OPS::copy<3>(cudaPol,kb_tris,"nrm",ktris,"nrm"); @@ -1958,10 +1787,7 @@ struct FleshDynamicStepping : INode { auto binderThicknessTag = get_param("binderThicknessTag"); auto binderInversionTag = get_param("binderInversionTag"); - auto planeConsPosTag = get_param("planeConsPosTag"); - auto planeConsNrmTag = get_param("planeConsNrmTag"); - auto planeConsIDTag = get_param("planeConsIDTag"); - auto planeConsBaryTag = get_param("planeConsBaryTag"); + auto planeConsStiffness = get_input2("planeConsStiffness"); @@ -2017,9 +1843,6 @@ struct FleshDynamicStepping : INode { auto max_cg_iters = get_param("max_cg_iters"); - bool use_plane_constraint = get_input2("use_plane_constraint"); - bool use_binder_constraint = get_input2("use_binder_constraint"); - bool use_kinematic_potential = get_input2("with_kinematic_potential"); bool use_line_search = get_param("use_line_search"); @@ -2076,6 +1899,21 @@ struct FleshDynamicStepping : INode { js[offset + stride * 3 + ei] = inds[2]; js[offset + stride * 4 + ei] = inds[3]; js[offset + stride * 5 + ei] = inds[3]; + + // js[offset + stride * 6 + ei] = inds[0]; + // js[offset + stride * 7 + ei] = inds[0]; + // js[offset + stride * 8 + ei] = inds[0]; + // js[offset + stride * 9 + ei] = inds[1]; + // js[offset + stride * 10 + ei] = inds[1]; + // js[offset + stride * 11 + ei] = inds[2]; + + // is[offset + stride * 6 + ei] = inds[1]; + // is[offset + stride * 7 + ei] = inds[2]; + // is[offset + stride * 8 + ei] = inds[3]; + // is[offset + stride * 9 + ei] = inds[2]; + // is[offset + stride * 10 + ei] = inds[3]; + // is[offset + stride * 11 + ei] = inds[3]; + }); spmat = spmat_t{verts.get_allocator(),(int)verts.size(),(int)verts.size()}; @@ -2174,7 +2012,7 @@ struct FleshDynamicStepping : INode { }else { std::cout << "apply no binder constraint" << std::endl; } - if(verts.hasProperty(planeConsPosTag) && verts.hasProperty(planeConsNrmTag) && verts.hasProperty(planeConsIDTag) && verts.hasProperty(planeConsBaryTag) && use_plane_constraint){ + if(verts.hasProperty(planeConsIDTag) && use_plane_constraint){ std::cout << "apply plane constraint" << std::endl; // A.computePlaneConstraintGradientAndHessian(cudaPol, @@ -2187,6 +2025,11 @@ struct FleshDynamicStepping : INode { planeConsIDTag, vtemp, sttemp,cnorm,use_sticky_condition); + + auto v2kt_force = TILEVEC_OPS::dot<3>(cudaPol,vtemp,"grad","grad"); + auto kv2t_force = TILEVEC_OPS::dot<9>(cudaPol,sttemp,"grad","grad"); + std::cout << "v2kt_force : " << v2kt_force << std::endl; + std::cout << "kv2t_force : " << kv2t_force << std::endl; } else{ std::cout << "apply no plane constraint : " << @@ -2218,7 +2061,7 @@ struct FleshDynamicStepping : INode { auto cforce = TILEVEC_OPS::dot<9>(cudaPol,sttemp,"grad","grad"); std::cout << "nm_csPT = " << nm_csPT << "\tkin_cforce : " << cforce << std::endl; #else - // std::cout << "apply kinematic collision" << std::endl; + std::cout << "apply kinematic collision" << std::endl; auto nm_csPT = COLLISION_UTILS::do_tetrahedra_surface_points_and_kinematic_boundary_collision_detection(cudaPol, kinematics[0], vtemp,"xn", @@ -2372,6 +2215,8 @@ struct FleshDynamicStepping : INode { out_collisionEps, in_collisionEps, csPT); + + // auto varea = TILEVEC_OPS::dot<1>(cudaPol,vtemp,"area","area"); // auto tarea = TILEVEC_OPS::dot<1>(cudaPol,tris,"area","area"); @@ -2395,8 +2240,8 @@ struct FleshDynamicStepping : INode { tris,"area", csPT,self_collision_fp_buffer,collisionStiffness,out_collisionEps); - // auto cforce = TILEVEC_OPS::dot<12>(cudaPol,fp_buffer,"grad","grad"); - // std::cout << "cforce " << cforce << std::endl; + auto cforce = TILEVEC_OPS::dot<12>(cudaPol,self_collision_fp_buffer,"grad","grad"); + std::cout << "cforce " << cforce << "\tnmcsPT : " << csPT.size() << std::endl; } timer.tock("eval hessian and gradient"); @@ -2620,103 +2465,4 @@ ZENDEFNODE(FleshDynamicStepping, {{"ZSParticles","kinematic_boundary", {"bool","use_line_search","0"} }, {"FEM"}}); - -// struct EvaluateElasticForce : zeno::INode { -// using T = float; -// using Ti = int; -// using dtiles_t = zs::TileVector; -// using tiles_t = typename ZenoParticles::particles_t; -// using vec2 = zs::vec; -// using vec3 = zs::vec; -// using mat3 = zs::vec; -// using mat9 = zs::vec; -// using mat12 = zs::vec; - -// using bvh_t = zs::LBvh<3,int,T>; -// using bv_t = zs::AABBBox<3, T>; - -// using pair3_t = zs::vec; -// using pair4_t = zs::vec; - -// virtual void apply() override { -// using namespace zs; -// auto zsparticles = get_input("ZSParticles"); -// auto models = zsparticles->getModel(); -// auto& verts = zsparticles->getParticles(); -// auto& eles = zsparticles->getQuadraturePoints(); - -// std::vector act_; -// std::size_t nm_acts = 0; - -// if(has_input("Acts")) { -// act_ = get_input("Acts")->getLiterial(); -// nm_acts = act_.size(); -// } - -// std::cout << "nmActs:" << nm_acts << std::endl; - -// constexpr auto host_space = zs::execspace_e::openmp; -// auto ompExec = zs::omp_exec(); -// auto act_buffer = dtiles_t{{{"act",2}},nm_acts,zs::memsrc_e::host}; -// ompExec(zs::range(act_buffer.size()), -// [act_buffer = proxy({},act_buffer),act_] (int i) mutable { -// act_buffer.tuple(dim_c<2>,"act",i) = vec2(act_[i][0],act_[i][1]); -// }); - -// act_buffer = act_buffer.clone({zs::memsrc_e::device, 0}); -// constexpr auto space = execspace_e::cuda; -// auto cudaPol = cuda_exec(); - -// auto forceTag = get_param("forceTag"); -// if(!verts.hasProperty(forceTag)){ -// verts.append_channels(cudaPol,{{forceTag,3}}); -// } -// TILEVEC_OPS::fill(cudaPol,verts,forceTag,(T)0.0); - - -// } -// }; - -// struct VisualizeBoneDrivenForce : zeno::INode { - -// }; - -// struct VisualizePlaneConstraintForce : zeno::INode { -// using T = float; -// using Ti = int; -// using dtiles_t = zs::TileVector; -// using tiles_t = typename ZenoParticles::particles_t; -// using vec2 = zs::vec; -// using vec3 = zs::vec; -// using mat3 = zs::vec; -// using mat9 = zs::vec; -// using mat12 = zs::vec; - -// using bvh_t = zs::LBvh<3,int,T>; -// using bv_t = zs::AABBBox<3, T>; - -// using pair3_t = zs::vec; -// using pair4_t = zs::vec; - -// virtual void apply() override { -// using namespace zs; -// auto zsparticles = get_input("ZSParticles"); -// auto& verts = zsparticles->getParticles(); -// auto& tris = (*zsparticles)[ZenoParticles::s_surfTriTag]; - -// auto kinematic_boundary = get_input("kinematic_boundary"); -// auto& kb_verts = kinematic_boundary->getParticles(); -// auto& kb_tris = kinematic_boundary->getQuadraturePoints(); - -// auto planeConsPosTag = get_param("planeConsPosTag"); -// auto planeConsNrmTag = get_param("planeConsNrmTag"); -// auto planeConsIDTag = get_param("planeConsIDTag"); -// auto planeConsBaryTag = get_param("planeConsBaryTag"); - -// auto planeConsStiffness = get_input2("planeConsStiffness"); - - -// } -// } - }; \ No newline at end of file diff --git a/projects/CuLagrange/fem/collision_energy/evaluate_collision.hpp b/projects/CuLagrange/fem/collision_energy/evaluate_collision.hpp index b0cb6e2564..858923649d 100644 --- a/projects/CuLagrange/fem/collision_energy/evaluate_collision.hpp +++ b/projects/CuLagrange/fem/collision_energy/evaluate_collision.hpp @@ -382,13 +382,14 @@ inline int do_tetrahedra_surface_points_and_kinematic_boundary_collision_detecti // std::cout << "finish write back gia res" << std::endl; auto ktris_cnorm = 2 * compute_average_edge_length(pol,kverts,"x",ktris); - auto bvh_thickness = ktris_cnorm > in_collisionEps ? ktris_cnorm : in_collisionEps; + // auto bvh_thickness = ktris_cnorm > in_collisionEps ? ktris_cnorm : in_collisionEps; + auto bvh_thickness = ktris_cnorm; auto ktriBvh = bvh_t{}; auto ktriBvs = retrieve_bounding_volumes(pol,kverts,ktris,wrapv<3>{},(T)0,"x"); ktriBvh.build(pol,ktriBvs); - // std::cout << "do csPT testing" << std::endl; + std::cout << "do csPT testing with k_active channel : " << kverts.hasProperty("k_active") << std::endl; bool colliding_from_inside = true; pol(zs::range(surf_verts_buffer.size()),[ @@ -410,18 +411,24 @@ inline int do_tetrahedra_surface_points_and_kinematic_boundary_collision_detecti pos_attr_tag = zs::SmallString(tet_pos_attr_tag), gia_res = proxy(gia_res), tris_gia_res = proxy(tris_gia_res)] ZS_LAMBDA(int vi) mutable { + if(verts_buffer("active",vi) < (T)0.5) + return; auto p = verts_buffer.pack(dim_c<3>,pos_attr_tag,vi); auto bv = bv_t{get_bounding_box(p - thickness,p + thickness)}; - T min_penertration_distance = (T)1e8; + T min_penertration_distance = limits::max(); int min_kti = -1; auto vnrm = verts_buffer.pack(dim_c<3>,"nrm",vi); auto process_vertex_kface_collision_pairs = [&](int kti) mutable { auto ktri = ktris.pack(dim_c<3>,"inds",kti,int_c); - if(kverts.hasProperty("k_active")) - for(int i = 0;i != 3;++i) - if(kverts("k_active",ktri[i]) < (T)0.5) - return; + // printf("testing pairs[%d %d] : %f %f %f\n",vi,kti, + // (float)kverts("k_active",ktri[0]), + // (float)kverts("k_active",ktri[1]), + // (float)kverts("k_active",ktri[2])); + // if(kverts.hasProperty("k_active")) + // for(int i = 0;i != 3;++i) + // if(kverts("k_active",ktri[i]) < (T)0.5) + // return; vec3 ktvs[3] = {}; for(int i = 0;i != 3;++i) ktvs[i] = kverts.pack(dim_c<3>,"x",ktri[i]); @@ -439,9 +446,30 @@ inline int do_tetrahedra_surface_points_and_kinematic_boundary_collision_detecti auto barySum = (T)1.0; T distance = LSL_GEO::pointTriangleDistance(ktvs[0],ktvs[1],ktvs[2],p,barySum); + // printf("testing pairs[%d %d] : %f %f\n",vi,kti,(float)distance,(float)collisionEps); + if(distance > collisionEps) return; + // printf("testing pairs[%d %d]\n",vi,kti); + int RING_MASK = 0; + for(int i = 0;i != ring_mask_width;++i) { + RING_MASK |= tris_gia_res[(kti + kt_offset) * ring_mask_width + i] & gia_res[vi * ring_mask_width + i]; + } + + // int RING_MASK = tris_gia_res[kti + kt_offset] & gia_res[vi]; + if(dist < 0 && RING_MASK == 0) { + // printf("negative distance but ring mask not matched\n"); + return; + } + if(dist > 0 && RING_MASK > 0) { + // printf("positive distance but ring mask matched\n"); + return; + } + + + + auto khi = zs::reinterpret_bits(ktris("he_inds",kti)); vec3 kbnrms[3] = {}; for(int i = 0;i != 3;++i) { @@ -470,22 +498,23 @@ inline int do_tetrahedra_surface_points_and_kinematic_boundary_collision_detecti } } - int RING_MASK = 0; - for(int i = 0;i != ring_mask_width;++i) { - RING_MASK |= tris_gia_res[(kti + kt_offset) * ring_mask_width + i] & gia_res[vi * ring_mask_width + i]; + if(dist < 0 && distance < min_penertration_distance){ + min_penertration_distance = distance; + min_kti = kti; } - // int RING_MASK = tris_gia_res[kti + kt_offset] & gia_res[vi]; - if(dist < 0 && RING_MASK == 0) - return; - if(dist > 0 && RING_MASK > 0) - return; - csPT.insert(zs::vec{vi,kti + ktoffset}); + if(dist > 0) + csPT.insert(zs::vec{vi,kti + ktoffset}); }; ktriBvh.iter_neighbors(bv,process_vertex_kface_collision_pairs); + + if(min_kti >= 0) + csPT.insert(zs::vec{vi,min_kti + ktoffset}); }); + + return csPT.size(); } @@ -785,8 +814,8 @@ inline void do_tetrahedra_surface_tris_and_points_self_collision_detection(Pol& PosTileVec& tet_verts,const zs::SmallString& pos_attr_tag, const TetTileVec& tets, const SurfPointTileVec& surf_points,const SurfTriTileVec& surf_tris, - const HalfEdgeTileVec& surf_halfedge, - const HalfFacetTileVec& tet_halffacet, + HalfEdgeTileVec& surf_halfedge, + HalfFacetTileVec& tet_halffacet, REAL outCollisionEps, REAL inCollisionEps, HashMap& csPT, @@ -798,28 +827,50 @@ inline void do_tetrahedra_surface_tris_and_points_self_collision_detection(Pol& PosTileVec surf_verts_buffer{surf_points.get_allocator(),{ {pos_attr_tag,3}, - {"ring_mask",1}, - {"color_mask",1}, - {"type_mask",1}, + // {"ring_mask",1}, + // {"color_mask",1}, + // {"type_mask",1}, {"embed_tet_id",1}, - {"mustExclude",1}, - {"is_corner",1}, + // {"mustExclude",1}, + // {"is_loop_vertex",1}, {"active",1} },surf_points.size()}; + // PosTileVec gia_res{surf_points.get_allocator(),{ + // {"color_mask",1} + // },(size_t)0}; SurfTriTileVec surf_tris_buffer{surf_tris.get_allocator(),{ {"inds",3}, {"he_inds",1}, + // {"ring_mask",1}, + // {"color_mask",1}, + // {"type_mask",1} + },surf_tris.size()}; + + PosTileVec gia_res{surf_points.get_allocator(),{ + // {pos_attr_tag,3}, + {"ring_mask",1}, + {"color_mask",1}, + {"type_mask",1}, + // {"mustExclude",1}, + {"is_loop_vertex",1}, + // {"active",1} + },(size_t)0}; + + SurfTriTileVec tris_gia_res{surf_tris.get_allocator(),{ + // {"inds",3}, + // {"he_inds",1}, {"ring_mask",1}, {"color_mask",1}, {"type_mask",1} - },surf_tris.size()}; + },(size_t)0}; + topological_sample(pol,surf_points,tet_verts,pos_attr_tag,surf_verts_buffer); - TILEVEC_OPS::fill(pol,surf_verts_buffer,"ring_mask",zs::reinterpret_bits((int)0)); - TILEVEC_OPS::fill(pol,surf_verts_buffer,"color_mask",zs::reinterpret_bits((int)0)); - TILEVEC_OPS::fill(pol,surf_verts_buffer,"type_mask",zs::reinterpret_bits((int)0)); + // TILEVEC_OPS::fill(pol,surf_verts_buffer,"ring_mask",zs::reinterpret_bits((int)0)); + // TILEVEC_OPS::fill(pol,surf_verts_buffer,"color_mask",zs::reinterpret_bits((int)0)); + // TILEVEC_OPS::fill(pol,surf_verts_buffer,"type_mask",zs::reinterpret_bits((int)0)); if(tet_verts.hasProperty("active")) { topological_sample(pol,surf_points,tet_verts,"active",surf_verts_buffer); }else { @@ -830,19 +881,19 @@ inline void do_tetrahedra_surface_tris_and_points_self_collision_detection(Pol& TILEVEC_OPS::copy(pol,surf_tris,"inds",surf_tris_buffer,"inds"); TILEVEC_OPS::copy(pol,surf_tris,"he_inds",surf_tris_buffer,"he_inds"); reorder_topology(pol,surf_points,surf_tris_buffer); - TILEVEC_OPS::fill(pol,surf_tris_buffer,"ring_mask",zs::reinterpret_bits((int)0)); - TILEVEC_OPS::fill(pol,surf_tris_buffer,"color_mask",zs::reinterpret_bits((int)0)); - TILEVEC_OPS::fill(pol,surf_tris_buffer,"type_mask",zs::reinterpret_bits((int)0)); + // TILEVEC_OPS::fill(pol,surf_tris_buffer,"ring_mask",zs::reinterpret_bits((int)0)); + // TILEVEC_OPS::fill(pol,surf_tris_buffer,"color_mask",zs::reinterpret_bits((int)0)); + // TILEVEC_OPS::fill(pol,surf_tris_buffer,"type_mask",zs::reinterpret_bits((int)0)); - auto cnorm = 3 * compute_average_edge_length(pol,surf_verts_buffer,pos_attr_tag,surf_tris_buffer); + auto cnorm = 5 * compute_average_edge_length(pol,surf_verts_buffer,pos_attr_tag,surf_tris_buffer); // evaluate_embed_tet_id auto tetBvh = bvh_t{}; auto tetBvs = retrieve_bounding_volumes(pol,tet_verts,tets,wrapv<4>{},(T)0,pos_attr_tag); tetBvh.build(pol,tetBvs); TILEVEC_OPS::fill(pol,surf_verts_buffer,"embed_tet_id",zs::reinterpret_bits((int)-1)); - zs::Vector nmExcludedPoints{tets.get_allocator(),1}; - nmExcludedPoints.setVal(0); + // zs::Vector nmExcludedPoints{tets.get_allocator(),1}; + // nmExcludedPoints.setVal(0); pol(zs::range(surf_verts_buffer.size()),[ pos_attr_tag = zs::SmallString(pos_attr_tag), surf_verts_buffer = proxy({},surf_verts_buffer), @@ -851,7 +902,7 @@ inline void do_tetrahedra_surface_tris_and_points_self_collision_detection(Pol& tets = proxy({},tets), thickness = cnorm, exec_tag, - nmExcludedPoints = proxy(nmExcludedPoints), + // nmExcludedPoints = proxy(nmExcludedPoints), tet_verts = proxy({},tet_verts)] ZS_LAMBDA(int pi) mutable { auto pv = surf_verts_buffer.pack(dim_c<3>,pos_attr_tag,pi); auto vi = zs::reinterpret_bits(surf_points("inds",pi)); @@ -867,44 +918,49 @@ inline void do_tetrahedra_surface_tris_and_points_self_collision_detection(Pol& tV[i] = tet_verts.pack(dim_c<3>,pos_attr_tag,tet[i]); if(LSL_GEO::is_inside_tet(tV[0],tV[1],tV[2],tV[3],pv)){ surf_verts_buffer("embed_tet_id",pi) = zs::reinterpret_bits((int)ei); - atomic_add(exec_tag,&nmExcludedPoints[0],(int)1); + // atomic_add(exec_tag,&nmExcludedPoints[0],(int)1); } }; tetBvh.iter_neighbors(bv,mark_interior_verts); }); - std::cout << "nm_excluded_points :" << nmExcludedPoints.getVal(0) << "\t" << surf_verts_buffer.size() << std::endl; + // std::cout << "nm_excluded_points :" << nmExcludedPoints.getVal(0) << "\t" << surf_verts_buffer.size() << std::endl; - pol(zs::range(surf_verts_buffer.size()),[ - surf_verts_buffer = proxy({},surf_verts_buffer)] ZS_LAMBDA(int pi) mutable { - auto embed_tet_id = zs::reinterpret_bits(surf_verts_buffer("embed_tet_id",pi)); - if(embed_tet_id >= 0) - surf_verts_buffer("mustExclude",pi) = (T)1.0; - else - surf_verts_buffer("mustExclude",pi) = (T)0.0; - }); - auto nm_rings = do_global_self_intersection_analysis(pol, + // pol(zs::range(surf_verts_buffer.size()),[ + // surf_verts_buffer = proxy({},surf_verts_buffer)] ZS_LAMBDA(int pi) mutable { + // auto embed_tet_id = zs::reinterpret_bits(surf_verts_buffer("embed_tet_id",pi)); + // if(embed_tet_id >= 0) + // surf_verts_buffer("mustExclude",pi) = (T)1.0; + // else + // surf_verts_buffer("mustExclude",pi) = (T)0.0; + // }); + auto ring_mask_width = do_global_self_intersection_analysis(pol, surf_verts_buffer,pos_attr_tag,surf_tris_buffer,surf_halfedge, - surf_verts_buffer,surf_tris_buffer); - std::cout << "nm_rings of GIA " << nm_rings << std::endl; + gia_res,tris_gia_res); + std::cout << "ring_mask_width of GIA : " << ring_mask_width << std::endl; if(write_back_gia_res) { - if(!tet_verts.hasProperty("ring_mask")) { - tet_verts.append_channels(pol,{{"ring_mask",1}}); - } + // if(!tet_verts.hasProperty("ring_mask")) { + // tet_verts.append_channels(pol,{{"ring_mask",1}}); + // } if(!tet_verts.hasProperty("flood")) { tet_verts.append_channels(pol,{{"flood",1}}); } - TILEVEC_OPS::fill(pol,tet_verts,"ring_mask",zs::reinterpret_bits((int)0)); + // TILEVEC_OPS::fill(pol,tet_verts,"ring_mask",zs::reinterpret_bits((int)0)); TILEVEC_OPS::fill(pol,tet_verts,"flood",(T)0); - pol(zs::range(surf_verts_buffer.size()),[ + pol(zs::range(surf_points.size()),[ tet_verts = proxy({},tet_verts), surf_points = proxy({},surf_points), - surf_verts_buffer = proxy({},surf_verts_buffer)] ZS_LAMBDA(int pi) mutable { + ring_mask_width = ring_mask_width, + gia_res = proxy({},gia_res)] ZS_LAMBDA(int pi) mutable { auto vi = zs::reinterpret_bits(surf_points("inds",pi)); - tet_verts("ring_mask",vi) = surf_verts_buffer("ring_mask",pi); - auto ring_mask = zs::reinterpret_bits(surf_verts_buffer("ring_mask",pi)); - if(ring_mask > 0) - tet_verts("flood",vi) = (T)1.0; + // tet_verts("ring_mask",vi) = surf_verts_buffer("ring_mask",pi); + for(int i = 0;i != ring_mask_width;++i) { + auto ring_mask = zs::reinterpret_bits(gia_res("ring_mask",pi * ring_mask_width + i)); + if(ring_mask > 0) { + tet_verts("flood",vi) = (T)1.0; + return; + } + } }); } @@ -922,6 +978,9 @@ inline void do_tetrahedra_surface_tris_and_points_self_collision_detection(Pol& surf_halfedge = proxy({},surf_halfedge), tets = proxy({},tets), tet_verts = proxy({},tet_verts), + gia_res = proxy({},gia_res), + tris_gia_res= proxy({},tris_gia_res), + ring_mask_width = ring_mask_width, thickness = cnorm, csPT = proxy(csPT), spBvh = proxy(spBvh)] ZS_LAMBDA(int ti) mutable { @@ -958,7 +1017,7 @@ inline void do_tetrahedra_surface_tris_and_points_self_collision_detection(Pol& hi = zs::reinterpret_bits(surf_halfedge("next_he",hi)); } - T min_penertration_depth = 1e8; + T min_penertration_depth = zs::limits::max(); int min_vi = -1; auto process_vertex_face_collision_pairs = [&, pos_attr_tag](int vi) { if(tri[0] == vi || tri[1] == vi || tri[2] == vi) @@ -995,74 +1054,309 @@ inline void do_tetrahedra_surface_tris_and_points_self_collision_detection(Pol& // if(dist < 0) { // auto neighbor_tet = zs::reinterpret_bits(halffacets("")) // } + bool is_valid_inverted_pair = false; + for(int bi = 0;bi != ring_mask_width;++bi) { + if(dist < 0) { + // do gia intersection test + int RING_MASK = zs::reinterpret_bits(gia_res("ring_mask",vi * ring_mask_width + bi)); + if(RING_MASK == 0) + continue; + // bool is_same_ring = false; + int TRING_MASK = ~0; + for(int i = 0;i != 3;++i) { + TRING_MASK &= zs::reinterpret_bits(gia_res("ring_mask",tri[i] * ring_mask_width + bi)); + } + RING_MASK = RING_MASK & TRING_MASK; + // the point and the tri should belong to the same ring + if(RING_MASK == 0) + continue; - if(dist < 0) { - // do gia intersection test - int RING_MASK = zs::reinterpret_bits(surf_verts_buffer("ring_mask",vi)); - if(RING_MASK == 0) - return; - // bool is_same_ring = false; - int TRING_MASK = 0; - for(int i = 0;i != 3;++i) { - TRING_MASK |= zs::reinterpret_bits(surf_verts_buffer("ring_mask",tri[i])); + // now the two pair belong to the same ring, check whether they belong black-white loop, and have different colors + auto COLOR_MASK = reinterpret_bits(gia_res("color_mask",vi * ring_mask_width + bi)); + auto TYPE_MASK = reinterpret_bits(gia_res("type_mask",vi * ring_mask_width + bi)); + + // only check the common type-1(white-black loop) rings + int TTYPE_MASK = 0; + for(int i = 0;i != 3;++i) + TTYPE_MASK |= reinterpret_bits(gia_res("type_mask",tri[i] * ring_mask_width + bi)); + + RING_MASK &= (TYPE_MASK & TTYPE_MASK); + // type-0 ring + if(RING_MASK == 0) { + is_valid_inverted_pair = true; + break; + } + // int nm_common_rings = 0; + // while() + // as long as there is one ring in which the pair have different colors, neglect the pair + int curr_ri_mask = 1; + bool is_color_same = false; + for(;RING_MASK > 0;RING_MASK = RING_MASK >> 1,curr_ri_mask = curr_ri_mask << 1) { + if(RING_MASK & 1) { + // int TCOLOR_MASK = ~0; + // for(int i = 0;i != 3;++i) + // TCOLOR_MASK &= reinterpret_bits(gia_res("color_mask",tri[i] * ring_mask_width + bi)) & curr_ri_mask; + // int auto VCOLOR_MASK = reinterpret_bits(gia_res("color_mask",vi * ring_mask_width + bi)) & curr_ri_mask; + for(int i = 0;i != 3;++i) { + auto TCOLOR_MASK = reinterpret_bits(gia_res("color_mask",tri[i] * ring_mask_width + bi)) & curr_ri_mask; + auto VCOLOR_MASK = reinterpret_bits(gia_res("color_mask",vi * ring_mask_width + bi)) & curr_ri_mask; + if(TCOLOR_MASK == VCOLOR_MASK) { + is_color_same = true; + break; + } + } + if(is_color_same) + break; + } + } + + if(!is_color_same) { + // type-1 ring with different color + is_valid_inverted_pair = true; + break; + } + + // break; + + // embed_tet_id >= 0 + // do shortest path detection + // auto ei = embed_tet_id; + + + + // is_valid_pair = true; } - RING_MASK = RING_MASK & TRING_MASK; - // the point and the tri should belong to the same ring - if(RING_MASK == 0) + // if(!is_valid_inverted_pair) + // return; + } + + + if(ring_mask_width == 0 && dist < 0) { + return; + } + + if(dist < 0 && is_valid_inverted_pair){ + min_vi = vi; + min_penertration_depth = distance; + } + if(dist > 0) + csPT.insert(zs::vec{vi,ti}); + }; + spBvh.iter_neighbors(bv,process_vertex_face_collision_pairs); + if(min_vi >= 0) + csPT.insert(zs::vec{min_vi,ti}); + }); + + auto stBvh = bvh_t{}; + auto stBvs = retrieve_bounding_volumes(pol,tet_verts,surf_tris,wrapv<3>{},(T)cnorm,pos_attr_tag); + stBvh.build(pol,stBvs); + + // csPT.reset(pol,true); + // for each verts, find the closest tri + pol(zs::range(surf_verts_buffer.size()),[ + outCollisionEps = outCollisionEps, + inCollisionEps = inCollisionEps, + surf_verts_buffer = proxy({},surf_verts_buffer,"surf_verts_buffer_problem"), + surf_tris_buffer = proxy({},surf_tris_buffer), + pos_attr_tag = zs::SmallString(pos_attr_tag), + surf_halfedge = proxy({},surf_halfedge), + tets = proxy({},tets), + tet_verts = proxy({},tet_verts), + gia_res = proxy({},gia_res), + tris_gia_res= proxy({},tris_gia_res), + ring_mask_width = ring_mask_width, + thickness = cnorm, + csPT = proxy(csPT), + stBvh = proxy(stBvh)] ZS_LAMBDA(int pi) mutable { + auto p = surf_verts_buffer.pack(dim_c<3>,pos_attr_tag,pi); + auto bv = bv_t{get_bounding_box(p - thickness,p + thickness)}; + + if(surf_verts_buffer("active",pi) < (T)0.5) + return; + T min_penertration_depth = zs::limits::max(); + int min_ti = -1; + auto process_vertex_face_collision_pairs = [&, pos_attr_tag](int ti) { + auto tri = surf_tris_buffer.pack(dim_c<3>,"inds",ti,int_c); + for(int i = 0;i != 3;++i) + if(surf_verts_buffer("active",tri[i]) < (T)0.5) return; - // now the two pair belong to the same ring, check whether they belong black-white loop, and have different colors - auto COLOR_MASK = reinterpret_bits(surf_verts_buffer("color_mask",vi)); - auto TYPE_MASK = reinterpret_bits(surf_verts_buffer("type_mask",vi)); + zs::vec tvs[3] = {}; + for(int i = 0;i != 3;++i) + tvs[i] = surf_verts_buffer.pack(dim_c<3>,pos_attr_tag,tri[i]); - // only check the common type-1(white-black loop) rings - int TTYPE_MASK = 0; - for(int i = 0;i != 3;++i) - TTYPE_MASK |= reinterpret_bits(surf_verts_buffer("type_mask",tri[i])); - - RING_MASK &= (TYPE_MASK & TTYPE_MASK); - // int nm_common_rings = 0; - // while() - // as long as there is one ring in which the pair have different colors, neglect the pair - int curr_ri_mask = 1; - for(;RING_MASK > 0;RING_MASK = RING_MASK >> 1,curr_ri_mask = curr_ri_mask << 1) { - if(RING_MASK & 1) { - for(int i = 0;i != 3;++i) { - auto TCOLOR_MASK = reinterpret_bits(surf_verts_buffer("color_mask",tri[i])) & curr_ri_mask; - auto VCOLOR_MASK = reinterpret_bits(surf_verts_buffer("color_mask",vi)) & curr_ri_mask; - if(TCOLOR_MASK == VCOLOR_MASK) - return; + if(tri[0] == pi || tri[1] == pi || tri[2] == pi) + return; + + // auto tri_center = zs::vec::zeros(); + // for(int i = 0;i != 3;++i) + // tri_center += tvs[i]/(T)3.0; + + auto embed_tet_id = zs::reinterpret_bits(surf_verts_buffer("embed_tet_id",pi)); + // auto p = surf_verts_buffer.pack(dim_c<3>,pos_attr_tag,vi); + auto tnrm = LSL_GEO::facet_normal(tvs[0],tvs[1],tvs[2]); + auto seg = p - tvs[0]; + auto dist = seg.dot(tnrm); + + auto collisionEps = dist > 0 ? outCollisionEps : inCollisionEps; + + auto barySum = (T)1.0; + T distance = LSL_GEO::pointTriangleDistance(tvs[0],tvs[1],tvs[2],p,barySum); + + if(distance > collisionEps) + return; + + + if(dist < 0 && distance > min_penertration_depth) + return; + + // if(dist < 0) { + // auto neighbor_tet = zs::reinterpret_bits(halffacets("")) + // } + if(dist < 0 && embed_tet_id < 0) + return; + // printf("testing pair %d %d\n",pi,ti); + + bool is_valid_inverted_pair = false; + for(int bi = 0;bi != ring_mask_width;++bi) { + if(dist < 0) { + // do gia intersection test + int V_RING_MASK = zs::reinterpret_bits(gia_res("ring_mask",pi * ring_mask_width + bi)); + if(V_RING_MASK == 0) + continue; + // bool is_same_ring = false; + // int TRING_MASK = ~0; + int TRING_MASK = ~0; + for(int i = 0;i != 3;++i) { + TRING_MASK &= zs::reinterpret_bits(gia_res("ring_mask",tri[i] * ring_mask_width + bi)); + } + auto RING_MASK = V_RING_MASK & TRING_MASK; + // the point and the tri should belong to the same ring + if(RING_MASK == 0) + continue; + // else if(V_RING_MASK > 0 && TRING_MASK > 0){ + // is_valid_inverted_pair = true; + // break; + // }else { + // continue; + // } + + // now the two pair belong to the same ring, check whether they belong black-white loop, and have different colors + auto COLOR_MASK = reinterpret_bits(gia_res("color_mask",pi * ring_mask_width + bi)); + auto TYPE_MASK = reinterpret_bits(gia_res("type_mask",pi * ring_mask_width + bi)); + + // only check the common type-1(white-black loop) rings + int TTYPE_MASK = 0; + for(int i = 0;i != 3;++i) + TTYPE_MASK |= reinterpret_bits(gia_res("type_mask",tri[i] * ring_mask_width + bi)); + + RING_MASK &= (TYPE_MASK & TTYPE_MASK); + // type-0 ring + if(RING_MASK == 0) { + is_valid_inverted_pair = true; + break; + } + // int nm_common_rings = 0; + // while() + // as long as there is one ring in which the pair have different colors, neglect the pair + int curr_ri_mask = 1; + bool is_color_same = false; + for(;RING_MASK > 0;RING_MASK = RING_MASK >> 1,curr_ri_mask = curr_ri_mask << 1) { + if(RING_MASK & 1) { + // int TCOLOR_MASK = ~0; + // for(int i = 0;i != 3;++i) + // TCOLOR_MASK &= reinterpret_bits(gia_res("color_mask",tri[i] * ring_mask_width + bi)) & curr_ri_mask; + // int auto VCOLOR_MASK = reinterpret_bits(gia_res("color_mask",vi * ring_mask_width + bi)) & curr_ri_mask; + for(int i = 0;i != 3;++i) { + auto TCOLOR_MASK = reinterpret_bits(gia_res("color_mask",tri[i] * ring_mask_width + bi)) & curr_ri_mask; + auto VCOLOR_MASK = reinterpret_bits(gia_res("color_mask",pi * ring_mask_width + bi)) & curr_ri_mask; + if(TCOLOR_MASK == VCOLOR_MASK) { + is_color_same = true; + break; + } + } + if(is_color_same) + break; } } + + if(!is_color_same) { + // type-1 ring with different color + is_valid_inverted_pair = true; + break; + } + // else { + // printf("skip with the same color %d %d\n",pi,ti); + // } + + // break; + + // embed_tet_id >= 0 + // do shortest path detection + // auto ei = embed_tet_id; + + + + // is_valid_pair = true; } + // if(!is_valid_inverted_pair) + // return; + } - // embed_tet_id >= 0 - // do shortest path detection - // auto ei = embed_tet_id; + if(barySum > (T)(1.0 + 0.1)) { + auto hi = zs::reinterpret_bits(surf_tris_buffer("he_inds",ti)); + vec3 bnrms[3] = {}; + for(int i = 0;i != 3;++i){ + auto edge_normal = tnrm; + auto opposite_hi = zs::reinterpret_bits(surf_halfedge("opposite_he",hi)); + if(opposite_hi >= 0){ + auto nti = zs::reinterpret_bits(surf_halfedge("to_face",opposite_hi)); + auto ntri = surf_tris_buffer.pack(dim_c<3>,"inds",nti,int_c); + auto ntnrm = LSL_GEO::facet_normal( + surf_verts_buffer.pack(dim_c<3>,pos_attr_tag,ntri[0]), + surf_verts_buffer.pack(dim_c<3>,pos_attr_tag,ntri[1]), + surf_verts_buffer.pack(dim_c<3>,pos_attr_tag,ntri[2])); + edge_normal = tnrm + ntnrm; + edge_normal = edge_normal/(edge_normal.norm() + (T)1e-6); + } + auto e01 = tvs[(i + 1) % 3] - tvs[(i + 0) % 3]; + bnrms[i] = edge_normal.cross(e01).normalized(); + hi = zs::reinterpret_bits(surf_halfedge("next_he",hi)); + } - // min_vi = vi; - // min_penertration_depth = distance; - - }else { // dist > 0 - int RING_MASK = zs::reinterpret_bits(surf_verts_buffer("ring_mask",vi)); - if(RING_MASK > 0) - return; - int TRING_MASK = 0; - for(int i = 0;i != 3;++i) { - TRING_MASK |= zs::reinterpret_bits(surf_verts_buffer("ring_mask",tri[i])); + for(int i = 0;i != 3;++i){ + seg = p - tvs[i]; + if(bnrms[i].dot(seg) < 0) { + // printf("skip due to bisector normal check\n"); + return; + } } - if(TRING_MASK > 0) - return; + } - + if(ring_mask_width == 0 && dist < 0) { + // printf("empty_ring_mask and negative dist\n"); + return; } - csPT.insert(zs::vec{vi,ti}); + + if(dist < 0 && is_valid_inverted_pair){ + min_ti = ti; + min_penertration_depth = distance; + } + else if(dist > 0){ + csPT.insert(zs::vec{pi,ti}); + // printf("find_new_positive pair %d %d\n",pi,ti); + } + // else { + // printf("invalid inverted pair\n"); + // } }; - spBvh.iter_neighbors(bv,process_vertex_face_collision_pairs); - // if(min_vi >= 0) - // csPT.insert(zs::vec{min_vi,ti}); + stBvh.iter_neighbors(bv,process_vertex_face_collision_pairs); + if(min_ti >= 0) { + // printf("find_new_negative pair %d %d\n",pi,min_ti); + csPT.insert(zs::vec{pi,min_ti}); + } }); } @@ -1381,7 +1675,7 @@ inline void do_kinematic_point_collision_detection(Pol& cudaPol, stBvh.build(cudaPol,bvs); auto avgl = compute_average_edge_length(cudaPol,verts,xtag,tris); - auto bvh_thickness = 5 * avgl; + auto bvh_thickness = 2 * avgl; if(update_normal) { if(!calculate_facet_normal(cudaPol,verts,xtag,tris,nrmTris,"nrm")){ diff --git a/projects/CuLagrange/geometry/BasicGeoNodes.cu b/projects/CuLagrange/geometry/BasicGeoNodes.cu new file mode 100755 index 0000000000..b5b4d12721 --- /dev/null +++ b/projects/CuLagrange/geometry/BasicGeoNodes.cu @@ -0,0 +1,85 @@ +#include "kernel/bary_centric_weights.hpp" +#include "zensim/io/MeshIO.hpp" +#include "zensim/math/bit/Bits.h" +#include "zensim/types/Property.h" +#include +#include +#include +#include +#include +#include + +#include "zensim/container/Bcht.hpp" +#include "kernel/tiled_vector_ops.hpp" +#include "kernel/geo_math.hpp" + +#include + +namespace zeno { + +struct ZSComputeSurfaceArea : zeno::INode { + using T = float; + virtual void apply() override { + using namespace zs; + constexpr auto cuda_space = execspace_e::cuda; + auto cudaPol = cuda_exec(); + constexpr auto exec_tag = wrapv{}; + + auto zsparticles = get_input("zsparticles"); + auto& verts = zsparticles->getParticles(); + bool is_tet_volume_mesh = zsparticles->category == ZenoParticles::category_e::tet; + auto &tris = is_tet_volume_mesh ? (*zsparticles)[ZenoParticles::s_surfTriTag] : zsparticles->getQuadraturePoints(); + + auto attrName = get_param("attrName"); + if(!verts.hasProperty(attrName)) { + verts.append_channels(cudaPol,{{attrName,1}}); + } + TILEVEC_OPS::fill(cudaPol,verts,attrName,(T)0.0); + + if(!tris.hasProperty(attrName)) { + tris.append_channels(cudaPol,{{attrName,1}}); + } + TILEVEC_OPS::fill(cudaPol,verts,attrName,(T)0.0); + + zs::Vector nmIncidentTris{verts.get_allocator(),verts.size()}; + cudaPol(zs::range(nmIncidentTris),[] ZS_LAMBDA(int& count) mutable {count = 0;}); + + cudaPol(zs::range(tris.size()),[ + exec_tag, + attrName = zs::SmallString(attrName), + tris = proxy({},tris), + nmIncidentTris = proxy(nmIncidentTris), + verts = proxy({},verts)] ZS_LAMBDA(int ti) mutable { + auto tri = tris.pack(dim_c<3>,"inds",ti,int_c); + zs::vec tV[3] = {}; + for(int i = 0;i != 3;++i) + tV[i] = verts.pack(dim_c<3>,"x",tri[i]); + auto A = LSL_GEO::area(tV[0],tV[1],tV[2]); + tris(attrName,ti) = A; + for(int i = 0;i != 3;++i) { + atomic_add(exec_tag,&verts(attrName,tri[i]),A); + atomic_add(exec_tag,&nmIncidentTris[0],(int)1); + } + }); + + cudaPol(zs::range(verts.size()),[ + verts = proxy({},verts), + attrName = zs::SmallString(attrName), + nmIncidentTris = proxy(nmIncidentTris)] ZS_LAMBDA(int vi) mutable { + if(nmIncidentTris[vi] > 0) + verts(attrName,vi) = verts(attrName,vi) / (T)nmIncidentTris[vi]; + }); + + set_output("zsparticles",zsparticles); + } +}; + + +ZENDEFNODE(ZSComputeSurfaceArea, {{{"zsparticles"}}, + {{"zsparticles"}}, + { + {"string","attrName","area"} + }, + {"ZSGeometry"}}); + +}; \ No newline at end of file diff --git a/projects/CuLagrange/geometry/CollisionVis.cu b/projects/CuLagrange/geometry/CollisionVis.cu index ce0937ec93..beb0d4f859 100644 --- a/projects/CuLagrange/geometry/CollisionVis.cu +++ b/projects/CuLagrange/geometry/CollisionVis.cu @@ -460,7 +460,7 @@ namespace zeno { {"ring_mask",1}, {"type_mask",1}, {"color_mask",1}, - {"is_corner",1} + {"is_loop_vertex",1} },verts_buffer.size()}; dtiles_t tris_gia_res{tri_buffer.get_allocator(),{ @@ -473,7 +473,7 @@ namespace zeno { auto nm_insts = do_global_self_intersection_analysis_on_surface_mesh_info( cudaPol,verts_buffer,"x",tri_buffer,halfedges,inst_buffer_info,gia_res,false); // zs::bht conn_of_first_ring{halfedges.get_allocator(),halfedges.size()}; - auto nm_rings = do_global_self_intersection_analysis(cudaPol, + auto ring_mask_width = do_global_self_intersection_analysis(cudaPol, verts_buffer,"x",tri_buffer,halfedges,gia_res,tris_gia_res); @@ -483,14 +483,18 @@ namespace zeno { verts.append_channels(cudaPol,{{markTag,1}}); } TILEVEC_OPS::fill(cudaPol,verts,markTag,(T)0.0); - cudaPol(zs::range(gia_res.size()),[ + cudaPol(zs::range(verts_buffer.size()),[ gia_res = proxy({},gia_res), verts = proxy({},verts), + ring_mask_width = ring_mask_width, verts_buffer = proxy({},verts_buffer), markTag = zs::SmallString(markTag) ] ZS_LAMBDA(int pi) mutable { auto vi = zs::reinterpret_bits(verts_buffer("inds",pi)); - auto ring_mask = zs::reinterpret_bits(gia_res("ring_mask",pi)); + int ring_mask = 0; + for(int i = 0;i != ring_mask_width;++i) { + ring_mask |= zs::reinterpret_bits(gia_res("ring_mask",pi * ring_mask_width + i)); + } verts(markTag,vi) = ring_mask == 0 ? (T)0.0 : (T)1.0; }); set_output("zsparticles",zsparticles); @@ -2221,18 +2225,18 @@ struct VisualizeSelfIntersections : zeno::INode { // zs::Vector nodal_colors{verts_buffer.get_allocator(),verts_buffer.size()}; // zs::Vector> instBuffer{tri_buffer.get_allocator(),tri_buffer.size() * 2}; - dtiles_t inst_buffer_info{tris.get_allocator(),{ - {"pair",2}, - {"type",1}, - {"its_edge_mark",6}, - {"int_points",6} - },tris.size() * 2}; + // dtiles_t inst_buffer_info{tris.get_allocator(),{ + // {"pair",2}, + // {"type",1}, + // {"its_edge_mark",6}, + // {"int_points",6} + // },tris.size() * 2}; dtiles_t gia_res{verts_buffer.get_allocator(),{ {"ring_mask",1}, {"type_mask",1}, {"color_mask",1}, - {"is_corner",1} + {"is_loop_vertex",1} },verts_buffer.size()}; dtiles_t tris_gia_res{tri_buffer.get_allocator(),{ {"ring_mask",1}, @@ -2241,10 +2245,10 @@ struct VisualizeSelfIntersections : zeno::INode { },tri_buffer.size()}; auto& halfedges = (*zsparticles)[ZenoParticles::s_surfHalfEdgeTag]; - auto nm_insts = do_global_self_intersection_analysis_on_surface_mesh_info( - cudaPol,verts_buffer,"x",tri_buffer,halfedges,inst_buffer_info,gia_res,false); + // auto nm_insts = do_global_self_intersection_analysis_on_surface_mesh_info( + // cudaPol,verts_buffer,"x",tri_buffer,halfedges,inst_buffer_info,gia_res,false); // zs::bht conn_of_first_ring{halfedges.get_allocator(),halfedges.size()}; - auto nm_rings = do_global_self_intersection_analysis(cudaPol, + auto ring_mask_width = do_global_self_intersection_analysis(cudaPol, verts_buffer,"x",tri_buffer,halfedges,gia_res,tris_gia_res); @@ -2380,177 +2384,186 @@ struct VisualizeSelfIntersections : zeno::INode { flood_region_vis->resize(verts.size()); auto& flood_region_verts = flood_region_vis->verts; auto& flood_region_mark = flood_region_vis->add_attr("flood"); - auto& is_corner_mark = flood_region_vis->add_attr("is_corner"); + auto& is_corner_mark = flood_region_vis->add_attr("is_loop_vertex"); ompPol(zs::range(verts_buffer.size()),[ &flood_region_verts, &flood_region_mark, &is_corner_mark, + ring_mask_width = ring_mask_width, flood_region = proxy({},flood_region), gia_res = proxy({},gia_res)] (int vi) mutable { auto p = flood_region.pack(dim_c<3>,"x",vi); flood_region_verts[vi] = p.to_array(); - auto ring_mask = zs::reinterpret_bits(gia_res("ring_mask",vi)); - flood_region_mark[vi] = ring_mask == 0 ? (float)0.0 : (float)1.0; - // auto is_corner = gia_res("is_corner",vi); - is_corner_mark[vi] = gia_res("is_corner",vi); - }); - set_output("flood_region",std::move(flood_region_vis)); + int ring_mask = 0; + bool is_corner = false; + for(int i = 0;i != ring_mask_width;++i) { + ring_mask |= zs::reinterpret_bits(gia_res("ring_mask",vi * ring_mask_width + i)); + if(gia_res("is_loop_vertex",vi) > (T)0.5) + is_corner = true; + } - dtiles_t self_intersect_buffer{tris.get_allocator(),{ - {"a0",3},{"A0",3}, - {"a1",3},{"A1",3}, - {"a2",3},{"A2",3}, - {"b0",3},{"B0",3}, - {"b1",3},{"B1",3}, - {"b2",3},{"B2",3}, - {"p0",3},{"p1",3} - },(size_t)nm_insts}; - cudaPol(zs::range(nm_insts),[ - // instBuffer = proxy(instBuffer), - inst_buffer_info = proxy({},inst_buffer_info), - verts_buffer = proxy({},verts_buffer), - self_intersect_buffer = proxy({},self_intersect_buffer), - tri_buffer = proxy({},tri_buffer)] ZS_LAMBDA(int sti) mutable { - auto tpair = inst_buffer_info.pack(dim_c<2>,"pair",sti,int_c); - auto ta = tpair[0]; - auto tb = tpair[1]; - - auto ints_p = inst_buffer_info.pack(dim_c<6>,"int_points",sti); - self_intersect_buffer.tuple(dim_c<3>,"p0",sti) = zs::vec{ints_p[0],ints_p[1],ints_p[2]}; - self_intersect_buffer.tuple(dim_c<3>,"p1",sti) = zs::vec{ints_p[3],ints_p[4],ints_p[5]}; - // auto ta = instBuffer[sti][0]; - // auto tb = instBuffer[sti][1]; - - auto triA = tri_buffer.pack(dim_c<3>,"inds",ta,int_c); - auto triB = tri_buffer.pack(dim_c<3>,"inds",tb,int_c); - self_intersect_buffer.tuple(dim_c<3>,"a0",sti) = verts_buffer.pack(dim_c<3>,"x",triA[0]); - self_intersect_buffer.tuple(dim_c<3>,"a1",sti) = verts_buffer.pack(dim_c<3>,"x",triA[1]); - self_intersect_buffer.tuple(dim_c<3>,"a2",sti) = verts_buffer.pack(dim_c<3>,"x",triA[2]); - - self_intersect_buffer.tuple(dim_c<3>,"b0",sti) = verts_buffer.pack(dim_c<3>,"x",triB[0]); - self_intersect_buffer.tuple(dim_c<3>,"b1",sti) = verts_buffer.pack(dim_c<3>,"x",triB[1]); - self_intersect_buffer.tuple(dim_c<3>,"b2",sti) = verts_buffer.pack(dim_c<3>,"x",triB[2]); - - self_intersect_buffer.tuple(dim_c<3>,"A0",sti) = verts_buffer.pack(dim_c<3>,"x",triA[0]); - self_intersect_buffer.tuple(dim_c<3>,"A1",sti) = verts_buffer.pack(dim_c<3>,"x",triA[1]); - self_intersect_buffer.tuple(dim_c<3>,"A2",sti) = verts_buffer.pack(dim_c<3>,"x",triA[2]); - - self_intersect_buffer.tuple(dim_c<3>,"B0",sti) = verts_buffer.pack(dim_c<3>,"x",triB[0]); - self_intersect_buffer.tuple(dim_c<3>,"B1",sti) = verts_buffer.pack(dim_c<3>,"x",triB[1]); - self_intersect_buffer.tuple(dim_c<3>,"B2",sti) = verts_buffer.pack(dim_c<3>,"x",triB[2]); + flood_region_mark[vi] = ring_mask == 0 ? (float)0.0 : (float)1.0; + // auto is_corner = gia_res("is_loop_vertex",vi); + is_corner_mark[vi] = is_corner ? (T)1.0 : (T)0.0; }); + set_output("flood_region",std::move(flood_region_vis)); - self_intersect_buffer = self_intersect_buffer.clone({zs::memsrc_e::host}); - - auto st_fact_vis = std::make_shared(); - auto& st_verts = st_fact_vis->verts; - auto& st_tris = st_fact_vis->tris; - st_verts.resize(self_intersect_buffer.size() * 6); - st_tris.resize(self_intersect_buffer.size() * 2); - - ompPol(zs::range(nm_insts),[ - &st_verts,&st_tris,self_intersect_buffer = proxy({},self_intersect_buffer)] (int sti) mutable { - st_verts[sti * 6 + 0] = self_intersect_buffer.pack(dim_c<3>,"a0",sti).to_array(); - st_verts[sti * 6 + 1] = self_intersect_buffer.pack(dim_c<3>,"a1",sti).to_array(); - st_verts[sti * 6 + 2] = self_intersect_buffer.pack(dim_c<3>,"a2",sti).to_array(); - st_verts[sti * 6 + 3] = self_intersect_buffer.pack(dim_c<3>,"b0",sti).to_array(); - st_verts[sti * 6 + 4] = self_intersect_buffer.pack(dim_c<3>,"b1",sti).to_array(); - st_verts[sti * 6 + 5] = self_intersect_buffer.pack(dim_c<3>,"b2",sti).to_array(); - - st_tris[sti * 2 + 0] = zeno::vec3i(sti * 6 + 0,sti * 6 + 1,sti * 6 + 2); - st_tris[sti * 2 + 1] = zeno::vec3i(sti * 6 + 3,sti * 6 + 4,sti * 6 + 5); - }); - - // std::cout << "nm_insts : " << nm_insts << std::endl; - set_output("st_facet_vis",std::move(st_fact_vis)); - - auto st_ring_vis = std::make_shared(); - auto& its_ring_verts = st_ring_vis->verts; - auto& its_ring_lines = st_ring_vis->lines; - its_ring_verts.resize(nm_insts * 2); - its_ring_lines.resize(nm_insts); - ompPol(zs::range(nm_insts),[ - &its_ring_verts,&its_ring_lines,self_intersect_buffer = proxy({},self_intersect_buffer)] (int sti) mutable { - auto p0 = self_intersect_buffer.pack(dim_c<3>,"p0",sti); - auto p1 = self_intersect_buffer.pack(dim_c<3>,"p1",sti); - its_ring_verts[sti * 2 + 0] = p0.to_array(); - its_ring_verts[sti * 2 + 1] = p1.to_array(); - its_ring_lines[sti] = zeno::vec2i{sti * 2 + 0,sti * 2 + 1}; - }); - set_output("st_ring_vis",std::move(st_ring_vis)); - - auto st_facet_rest_vis = std::make_shared(); - auto& st_rest_verts = st_facet_rest_vis->verts; - auto& st_rest_tris = st_facet_rest_vis->tris; - st_rest_verts.resize(self_intersect_buffer.size() * 6); - st_rest_tris.resize(self_intersect_buffer.size() * 2); - ompPol(zs::range(nm_insts),[ - &st_rest_verts,&st_rest_tris,self_intersect_buffer = proxy({},self_intersect_buffer)] (int sti) mutable { - st_rest_verts[sti * 6 + 0] = self_intersect_buffer.pack(dim_c<3>,"A0",sti).to_array(); - st_rest_verts[sti * 6 + 1] = self_intersect_buffer.pack(dim_c<3>,"A1",sti).to_array(); - st_rest_verts[sti * 6 + 2] = self_intersect_buffer.pack(dim_c<3>,"A2",sti).to_array(); - st_rest_verts[sti * 6 + 3] = self_intersect_buffer.pack(dim_c<3>,"B0",sti).to_array(); - st_rest_verts[sti * 6 + 4] = self_intersect_buffer.pack(dim_c<3>,"B1",sti).to_array(); - st_rest_verts[sti * 6 + 5] = self_intersect_buffer.pack(dim_c<3>,"B2",sti).to_array(); - - st_rest_tris[sti * 2 + 0] = zeno::vec3i(sti * 6 + 0,sti * 6 + 1,sti * 6 + 2); - st_rest_tris[sti * 2 + 1] = zeno::vec3i(sti * 6 + 3,sti * 6 + 4,sti * 6 + 5); - }); - set_output("st_facet_rest_vis",std::move(st_facet_rest_vis)); - - dtiles_t st_pair_buffer{tris.get_allocator(),{ - {"x0",3}, - {"x1",3} - },(std::size_t)nm_insts}; - cudaPol(zs::range(nm_insts),[ - inst_buffer_info = proxy({},inst_buffer_info), - // instBuffer = proxy(instBuffer), - st_pair_buffer = proxy({},st_pair_buffer), - verts = proxy({},verts_buffer), - tris = proxy({},tri_buffer)] ZS_LAMBDA(int sti) mutable { - auto tpair = inst_buffer_info.pack(dim_c<2>,"pair",sti,int_c); - auto ta = tpair[0]; - auto tb = tpair[1]; - // auto ta = instBuffer[sti][0]; - // auto tb = instBuffer[sti][1]; - - - auto triA = tris.pack(dim_c<3>,"inds",ta,int_c); - auto triB = tris.pack(dim_c<3>,"inds",tb,int_c); - - auto x0 = vec3::zeros(); - auto x1 = vec3::zeros(); + // dtiles_t self_intersect_buffer{tris.get_allocator(),{ + // {"a0",3},{"A0",3}, + // {"a1",3},{"A1",3}, + // {"a2",3},{"A2",3}, + // {"b0",3},{"B0",3}, + // {"b1",3},{"B1",3}, + // {"b2",3},{"B2",3}, + // // {"p0",3},{"p1",3} + // },(size_t)nm_insts}; + // cudaPol(zs::range(nm_insts),[ + // // instBuffer = proxy(instBuffer), + // inst_buffer_info = proxy({},inst_buffer_info), + // verts_buffer = proxy({},verts_buffer), + // self_intersect_buffer = proxy({},self_intersect_buffer), + // tri_buffer = proxy({},tri_buffer)] ZS_LAMBDA(int sti) mutable { + // auto tpair = inst_buffer_info.pack(dim_c<2>,"pair",sti,int_c); + // auto ta = tpair[0]; + // auto tb = tpair[1]; + + // // auto ints_p = inst_buffer_info.pack(dim_c<6>,"int_points",sti); + // // self_intersect_buffer.tuple(dim_c<3>,"p0",sti) = zs::vec{ints_p[0],ints_p[1],ints_p[2]}; + // // self_intersect_buffer.tuple(dim_c<3>,"p1",sti) = zs::vec{ints_p[3],ints_p[4],ints_p[5]}; + // // auto ta = instBuffer[sti][0]; + // // auto tb = instBuffer[sti][1]; + + // auto triA = tri_buffer.pack(dim_c<3>,"inds",ta,int_c); + // auto triB = tri_buffer.pack(dim_c<3>,"inds",tb,int_c); + // self_intersect_buffer.tuple(dim_c<3>,"a0",sti) = verts_buffer.pack(dim_c<3>,"x",triA[0]); + // self_intersect_buffer.tuple(dim_c<3>,"a1",sti) = verts_buffer.pack(dim_c<3>,"x",triA[1]); + // self_intersect_buffer.tuple(dim_c<3>,"a2",sti) = verts_buffer.pack(dim_c<3>,"x",triA[2]); + + // self_intersect_buffer.tuple(dim_c<3>,"b0",sti) = verts_buffer.pack(dim_c<3>,"x",triB[0]); + // self_intersect_buffer.tuple(dim_c<3>,"b1",sti) = verts_buffer.pack(dim_c<3>,"x",triB[1]); + // self_intersect_buffer.tuple(dim_c<3>,"b2",sti) = verts_buffer.pack(dim_c<3>,"x",triB[2]); + + // self_intersect_buffer.tuple(dim_c<3>,"A0",sti) = verts_buffer.pack(dim_c<3>,"x",triA[0]); + // self_intersect_buffer.tuple(dim_c<3>,"A1",sti) = verts_buffer.pack(dim_c<3>,"x",triA[1]); + // self_intersect_buffer.tuple(dim_c<3>,"A2",sti) = verts_buffer.pack(dim_c<3>,"x",triA[2]); + + // self_intersect_buffer.tuple(dim_c<3>,"B0",sti) = verts_buffer.pack(dim_c<3>,"x",triB[0]); + // self_intersect_buffer.tuple(dim_c<3>,"B1",sti) = verts_buffer.pack(dim_c<3>,"x",triB[1]); + // self_intersect_buffer.tuple(dim_c<3>,"B2",sti) = verts_buffer.pack(dim_c<3>,"x",triB[2]); + // }); - for(int i = 0;i != 3;++i) { - x0 += verts.pack(dim_c<3>,"x",triA[i]) / (T)3.0; - x1 += verts.pack(dim_c<3>,"x",triB[i]) / (T)3.0; - } + // self_intersect_buffer = self_intersect_buffer.clone({zs::memsrc_e::host}); + + // auto st_fact_vis = std::make_shared(); + // auto& st_verts = st_fact_vis->verts; + // auto& st_tris = st_fact_vis->tris; + // st_verts.resize(self_intersect_buffer.size() * 6); + // st_tris.resize(self_intersect_buffer.size() * 2); + + // ompPol(zs::range(nm_insts),[ + // &st_verts,&st_tris,self_intersect_buffer = proxy({},self_intersect_buffer)] (int sti) mutable { + // st_verts[sti * 6 + 0] = self_intersect_buffer.pack(dim_c<3>,"a0",sti).to_array(); + // st_verts[sti * 6 + 1] = self_intersect_buffer.pack(dim_c<3>,"a1",sti).to_array(); + // st_verts[sti * 6 + 2] = self_intersect_buffer.pack(dim_c<3>,"a2",sti).to_array(); + // st_verts[sti * 6 + 3] = self_intersect_buffer.pack(dim_c<3>,"b0",sti).to_array(); + // st_verts[sti * 6 + 4] = self_intersect_buffer.pack(dim_c<3>,"b1",sti).to_array(); + // st_verts[sti * 6 + 5] = self_intersect_buffer.pack(dim_c<3>,"b2",sti).to_array(); + + // st_tris[sti * 2 + 0] = zeno::vec3i(sti * 6 + 0,sti * 6 + 1,sti * 6 + 2); + // st_tris[sti * 2 + 1] = zeno::vec3i(sti * 6 + 3,sti * 6 + 4,sti * 6 + 5); + // }); + + // // std::cout << "nm_insts : " << nm_insts << std::endl; + // set_output("st_facet_vis",std::move(st_fact_vis)); + + // auto st_ring_vis = std::make_shared(); + // auto& its_ring_verts = st_ring_vis->verts; + // auto& its_ring_lines = st_ring_vis->lines; + // its_ring_verts.resize(nm_insts * 2); + // its_ring_lines.resize(nm_insts); + // ompPol(zs::range(nm_insts),[ + // &its_ring_verts,&its_ring_lines,self_intersect_buffer = proxy({},self_intersect_buffer)] (int sti) mutable { + // auto p0 = self_intersect_buffer.pack(dim_c<3>,"p0",sti); + // auto p1 = self_intersect_buffer.pack(dim_c<3>,"p1",sti); + // its_ring_verts[sti * 2 + 0] = p0.to_array(); + // its_ring_verts[sti * 2 + 1] = p1.to_array(); + // its_ring_lines[sti] = zeno::vec2i{sti * 2 + 0,sti * 2 + 1}; + // }); - st_pair_buffer.tuple(dim_c<3>,"x0",sti) = x0.to_array(); - st_pair_buffer.tuple(dim_c<3>,"x1",sti) = x1.to_array(); - }); + // set_output("st_ring_vis",std::move(st_ring_vis)); + + // auto st_facet_rest_vis = std::make_shared(); + // auto& st_rest_verts = st_facet_rest_vis->verts; + // auto& st_rest_tris = st_facet_rest_vis->tris; + // st_rest_verts.resize(self_intersect_buffer.size() * 6); + // st_rest_tris.resize(self_intersect_buffer.size() * 2); + // ompPol(zs::range(nm_insts),[ + // &st_rest_verts,&st_rest_tris,self_intersect_buffer = proxy({},self_intersect_buffer)] (int sti) mutable { + // st_rest_verts[sti * 6 + 0] = self_intersect_buffer.pack(dim_c<3>,"A0",sti).to_array(); + // st_rest_verts[sti * 6 + 1] = self_intersect_buffer.pack(dim_c<3>,"A1",sti).to_array(); + // st_rest_verts[sti * 6 + 2] = self_intersect_buffer.pack(dim_c<3>,"A2",sti).to_array(); + // st_rest_verts[sti * 6 + 3] = self_intersect_buffer.pack(dim_c<3>,"B0",sti).to_array(); + // st_rest_verts[sti * 6 + 4] = self_intersect_buffer.pack(dim_c<3>,"B1",sti).to_array(); + // st_rest_verts[sti * 6 + 5] = self_intersect_buffer.pack(dim_c<3>,"B2",sti).to_array(); + + // st_rest_tris[sti * 2 + 0] = zeno::vec3i(sti * 6 + 0,sti * 6 + 1,sti * 6 + 2); + // st_rest_tris[sti * 2 + 1] = zeno::vec3i(sti * 6 + 3,sti * 6 + 4,sti * 6 + 5); + // }); + // set_output("st_facet_rest_vis",std::move(st_facet_rest_vis)); + + // dtiles_t st_pair_buffer{tris.get_allocator(),{ + // {"x0",3}, + // {"x1",3} + // },(std::size_t)nm_insts}; + // cudaPol(zs::range(nm_insts),[ + // inst_buffer_info = proxy({},inst_buffer_info), + // // instBuffer = proxy(instBuffer), + // st_pair_buffer = proxy({},st_pair_buffer), + // verts = proxy({},verts_buffer), + // tris = proxy({},tri_buffer)] ZS_LAMBDA(int sti) mutable { + // auto tpair = inst_buffer_info.pack(dim_c<2>,"pair",sti,int_c); + // auto ta = tpair[0]; + // auto tb = tpair[1]; + // // auto ta = instBuffer[sti][0]; + // // auto tb = instBuffer[sti][1]; + + + // auto triA = tris.pack(dim_c<3>,"inds",ta,int_c); + // auto triB = tris.pack(dim_c<3>,"inds",tb,int_c); + + // auto x0 = vec3::zeros(); + // auto x1 = vec3::zeros(); + + // for(int i = 0;i != 3;++i) { + // x0 += verts.pack(dim_c<3>,"x",triA[i]) / (T)3.0; + // x1 += verts.pack(dim_c<3>,"x",triB[i]) / (T)3.0; + // } + + // st_pair_buffer.tuple(dim_c<3>,"x0",sti) = x0.to_array(); + // st_pair_buffer.tuple(dim_c<3>,"x1",sti) = x1.to_array(); + // }); - st_pair_buffer = st_pair_buffer.clone({zs::memsrc_e::host}); - auto st_pair_vis = std::make_shared(); - auto& st_pair_verts = st_pair_vis->verts; - auto& st_pair_lines = st_pair_vis->lines; - st_pair_verts.resize(st_pair_buffer.size() * 2); - st_pair_lines.resize(st_pair_buffer.size()); - - ompPol(zs::range(st_pair_buffer.size()),[ - st_pair_buffer = proxy({},st_pair_buffer), - &st_pair_verts,&st_pair_lines] (int spi) mutable { - auto x0 = st_pair_buffer.pack(dim_c<3>,"x0",spi); - auto x1 = st_pair_buffer.pack(dim_c<3>,"x1",spi); - st_pair_verts[spi * 2 + 0] = x0.to_array(); - st_pair_verts[spi * 2 + 1] = x1.to_array(); - st_pair_lines[spi] = zeno::vec2i{spi * 2 + 0,spi * 2 + 1}; - }); + // st_pair_buffer = st_pair_buffer.clone({zs::memsrc_e::host}); + // auto st_pair_vis = std::make_shared(); + // auto& st_pair_verts = st_pair_vis->verts; + // auto& st_pair_lines = st_pair_vis->lines; + // st_pair_verts.resize(st_pair_buffer.size() * 2); + // st_pair_lines.resize(st_pair_buffer.size()); + + // ompPol(zs::range(st_pair_buffer.size()),[ + // st_pair_buffer = proxy({},st_pair_buffer), + // &st_pair_verts,&st_pair_lines] (int spi) mutable { + // auto x0 = st_pair_buffer.pack(dim_c<3>,"x0",spi); + // auto x1 = st_pair_buffer.pack(dim_c<3>,"x1",spi); + // st_pair_verts[spi * 2 + 0] = x0.to_array(); + // st_pair_verts[spi * 2 + 1] = x1.to_array(); + // st_pair_lines[spi] = zeno::vec2i{spi * 2 + 0,spi * 2 + 1}; + // }); - set_output("st_pair_vis",std::move(st_pair_vis)); + // set_output("st_pair_vis",std::move(st_pair_vis)); // dtiles_t corner_verts_buffer{gia_res.get_allocator()} } @@ -2558,9 +2571,9 @@ struct VisualizeSelfIntersections : zeno::INode { ZENDEFNODE(VisualizeSelfIntersections, {{"zsparticles"}, { - "st_ring_vis", - "st_facet_rest_vis", - "st_facet_vis", + // "st_ring_vis", + // "st_facet_rest_vis", + // "st_facet_vis", "flood_region", // "be_vis" // "wire_fr_vis" @@ -2947,7 +2960,7 @@ struct VisualizeCollision2 : zeno::INode { const auto& tris = (*zsparticles)[ZenoParticles::s_surfTriTag]; const auto& points = (*zsparticles)[ZenoParticles::s_surfVertTag]; const auto& tets = zsparticles->getQuadraturePoints(); - const auto& halfedges = (*zsparticles)[ZenoParticles::s_surfHalfEdgeTag]; + auto& halfedges = (*zsparticles)[ZenoParticles::s_surfHalfEdgeTag]; const auto& halffacets = (*zsparticles)[ZenoParticles::s_tetHalfFacetTag]; zs::bht csPT{verts.get_allocator(),10000}; diff --git a/projects/CuLagrange/geometry/SurfaceBinder.cu b/projects/CuLagrange/geometry/SurfaceBinder.cu index 48c5485e99..66ec6f4f8b 100644 --- a/projects/CuLagrange/geometry/SurfaceBinder.cu +++ b/projects/CuLagrange/geometry/SurfaceBinder.cu @@ -18,6 +18,8 @@ #include "kernel/tiled_vector_ops.hpp" #include "kernel/geo_math.hpp" +#include "kernel/topology.hpp" +#include "kernel/intersection.hpp" #include "zensim/container/Bvh.hpp" #include "zensim/container/Bvs.hpp" @@ -27,6 +29,7 @@ #include "kernel/compute_characteristic_length.hpp" + namespace zeno { struct ZSSurfaceBind : zeno::INode { @@ -95,6 +98,7 @@ struct ZSSurfaceBind : zeno::INode { auto binder_tag = get_param("binder_tag"); auto thickness_tag = get_param("thickness_tag"); auto inversion_tag = get_param("inversion_tag"); + auto binder_bary_tag = get_param("binder_bary_tag"); auto align_direction = get_param("align_direction"); @@ -102,12 +106,14 @@ struct ZSSurfaceBind : zeno::INode { {binder_tag,max_nm_binders}, {thickness_tag,max_nm_binders}, {inversion_tag,max_nm_binders}, + {binder_bary_tag,max_nm_binders * 3}, {"nm_binders",1} }); TILEVEC_OPS::fill(cudaPol,tris,binder_tag,zs::reinterpret_bits((int)-1)); TILEVEC_OPS::fill(cudaPol,tris,thickness_tag,(T)0.0); TILEVEC_OPS::fill(cudaPol,tris,inversion_tag,(T)-1.0); TILEVEC_OPS::fill(cudaPol,tris,"nm_binders",reinterpret_bits((int)0)); + TILEVEC_OPS::fill(cudaPol,tris,binder_bary_tag,(T)0.0); auto kpBvh = bvh_t{}; auto bvs = retrieve_bounding_volumes(cudaPol,kverts,kverts,wrapv<1>{},(T)0.0,"x"); @@ -154,6 +160,7 @@ struct ZSSurfaceBind : zeno::INode { binder_tag = zs::SmallString(binder_tag), thickness_tag = zs::SmallString(thickness_tag), inversion_tag = zs::SmallString(inversion_tag), + binder_bary_tag = zs::SmallString(binder_bary_tag), max_nm_binders = max_nm_binders, align_angle_cosin = align_angle_cosin, kinInCollisionEps = kinInCollisionEps, @@ -192,15 +199,18 @@ struct ZSSurfaceBind : zeno::INode { auto t2 = verts.pack(dim_c<3>,"x",tri[2]); T barySum = (T)0.0; - T distance = LSL_GEO::pointTriangleDistance(t0,t1,t2,kp,barySum); + vec3 project_bary{}; + vec3 bary{}; + T distance = LSL_GEO::pointTriangleDistance(t0,t1,t2,kp,bary,project_bary); + barySum = fabs(bary[0]) + fabs(bary[0]) + fabs(bary[0]); auto nrm = tris.pack(dim_c<3>,"nrm",ti); auto knrm = kverts.pack(dim_c<3>,"nrm",kpi); // alignment - if(nrm.dot(knrm) < (T)align_angle_cosin && align_direction) + if(nrm.dot(knrm) < (T)align_angle_cosin && align_direction && align_angle_cosin > 0) return; - if(nrm.dot(knrm) > (T)-align_angle_cosin && !align_direction) + if(nrm.dot(knrm) > (T)align_angle_cosin && align_direction && align_angle_cosin < 0) return; auto dist = seg.dot(nrm); @@ -241,8 +251,10 @@ struct ZSSurfaceBind : zeno::INode { nm_tag++; tris(binder_tag,nm_binders,ti) = reinterpret_bits(kpi); - tris(thickness_tag,nm_binders,ti) = distance; + tris(thickness_tag,nm_binders,ti) = distance / (T)5.0; tris(inversion_tag,nm_binders,ti) = dist < 0 ? (T)1.0 : (T)-1.0; + for(int i = 0;i != 3;++i) + tris(binder_bary_tag,nm_binders * 3 + i,ti) = project_bary[i]; nm_binders++; kb_verts(markTag,kpi) = (T)1.0; }; @@ -273,6 +285,7 @@ ZENDEFNODE(ZSSurfaceBind, {{"zssurf","kboundary", {"string","thickness_tag","thicknessTag"}, {"string","inversion_tag","inversionTag"}, {"string","mark_tag","markTag"}, + {"string","binder_bary_tag","binderBaryTag"}, {"bool","align_direction","1"}, }, {"ZSGeometry"}}); @@ -766,9 +779,12 @@ struct ZSSurfaceClosestPoints : zeno::INode { virtual void apply() override { using namespace zs; constexpr auto space = execspace_e::cuda; + constexpr auto exec_tag = wrapv{}; auto cudaPol = cuda_exec(); auto zsparticles = get_input("zsparticles"); + + auto kboundary = get_input("kboundary"); auto& verts = zsparticles->getParticles(); @@ -776,250 +792,363 @@ struct ZSSurfaceClosestPoints : zeno::INode { (*zsparticles)[ZenoParticles::s_surfTriTag] : zsparticles->getQuadraturePoints(); + auto& halfedges = (*zsparticles)[ZenoParticles::s_surfHalfEdgeTag]; + const auto& points = (*zsparticles)[ZenoParticles::s_surfVertTag]; + // every vertex can only bind to one triangle auto& kverts = kboundary->getParticles(); auto& ktris = kboundary->getQuadraturePoints(); const auto& khalfedges = (*kboundary)[ZenoParticles::s_surfHalfEdgeTag]; // auto max_nm_binder = get_input2("nm_max_binders"); - auto project_pos_tag = get_param("project_pos_tag"); - auto project_nrm_tag = get_param("project_nrm_tag"); auto project_idx_tag = get_param("project_idx_tag"); - auto project_bary_tag = get_param("project_bary_tag"); auto align_direction = get_param("align_direction"); + auto xtag = get_param("xtag"); + auto kxtag = get_param("kxtag"); - // for each vertex of zsparticles, find a potential closest point on kboundary surface - // add a plane constraint - if(!verts.hasProperty(project_pos_tag) || !verts.hasProperty(project_nrm_tag) || !verts.hasProperty(project_idx_tag) || !verts.hasProperty(project_bary_tag)) { - verts.append_channels(cudaPol,{ - {project_pos_tag,3},// the idx of triangle of kboudary - {project_nrm_tag,3}, - {project_bary_tag,3}, - {project_idx_tag,1} - }); - } - + if(!verts.hasProperty(project_idx_tag)) + verts.append_channels(cudaPol,{{project_idx_tag,1}}); + if(!kverts.hasProperty(project_idx_tag)) + kverts.append_channels(cudaPol,{{project_idx_tag,1}}); TILEVEC_OPS::fill(cudaPol,verts,project_idx_tag,zs::reinterpret_bits((int)-1)); - - auto ktBvh = bvh_t{}; - auto bvs = retrieve_bounding_volumes(cudaPol,kverts,ktris,wrapv<3>{},(T)0.0,"x"); - ktBvh.build(cudaPol,bvs); + TILEVEC_OPS::fill(cudaPol,kverts,project_idx_tag,zs::reinterpret_bits((int)-1)); - auto kinInCollisionEps = get_input2("kinInColEps"); - auto kinOutCollisionEps = get_input2("kinOutColEps"); - auto thickness = kinInCollisionEps + kinOutCollisionEps; + dtiles_t surf_verts_buffer{points.get_allocator(),{ + {"x",3}, + {"nrm",3}, + {"X",3}, + {"k_active",1} + },points.size()}; - if(!tris.hasProperty("nrm")) - tris.append_channels(cudaPol,{{"nrm",3}}); - cudaPol(zs::range(tris.size()), - [tris = proxy({},tris), - verts = proxy({},verts)] ZS_LAMBDA(int ti) { - auto tri = tris.template pack<3>("inds",ti).reinterpret_bits(int_c); - auto v0 = verts.template pack<3>("x",tri[0]); - auto v1 = verts.template pack<3>("x",tri[1]); - auto v2 = verts.template pack<3>("x",tri[2]); + dtiles_t surf_tris_buffer{tris.get_allocator(),{ + {"nrm",3}, + {"inds",3}, + {"he_inds",1} + },tris.size()}; - auto e01 = v1 - v0; - auto e02 = v2 - v0; + dtiles_t kverts_buffer{kverts.get_allocator(),{ + {"x",3}, + {"nrm",3}, + {"X",3}, + {"k_active",1} + },kverts.size()}; - auto nrm = e01.cross(e02); - auto nrm_norm = nrm.norm(); - if(nrm_norm < 1e-8) - nrm = zs::vec::zeros(); - else - nrm = nrm / nrm_norm; + dtiles_t ktris_buffer{ktris.get_allocator(),{ + {"nrm",3}, + {"inds",3}, + {"he_inds",1}, + },ktris.size()}; - tris.tuple(dim_c<3>,"nrm",ti) = nrm; + topological_sample(cudaPol,points,verts,xtag,surf_verts_buffer,"x"); + if(verts.hasProperty("X")) + topological_sample(cudaPol,points,verts,"X",surf_verts_buffer,"X"); + else + topological_sample(cudaPol,points,verts,xtag,surf_verts_buffer,"X"); + if(verts.hasProperty("k_active")) + topological_sample(cudaPol,points,verts,"k_active",surf_verts_buffer,"k_active"); + else + TILEVEC_OPS::fill(cudaPol,surf_verts_buffer,"k_active",(T)1.0); + TILEVEC_OPS::copy(cudaPol,tris,"inds",surf_tris_buffer,"inds"); + TILEVEC_OPS::copy(cudaPol,tris,"he_inds",surf_tris_buffer,"he_inds"); + reorder_topology(cudaPol,points,surf_tris_buffer); + + TILEVEC_OPS::fill(cudaPol,surf_verts_buffer,"nrm",(T)0.0); + cudaPol(zs::range(surf_tris_buffer.size()),[exec_tag, + surf_verts_buffer = proxy({},surf_verts_buffer), + surf_tris_buffer = proxy({},surf_tris_buffer)] ZS_LAMBDA(int ti) mutable { + auto tri = surf_tris_buffer.pack(dim_c<3>,"inds",ti,int_c); + zs::vec tV[3] = {}; + for(int i = 0;i != 3;++i) + tV[i] = surf_verts_buffer.pack(dim_c<3>,"x",tri[i]); + auto tnrm = LSL_GEO::facet_normal(tV[0],tV[1],tV[2]); + surf_tris_buffer.tuple(dim_c<3>,"nrm",ti) = tnrm; + for(int i = 0;i != 3;++i) + for(int d = 0;d != 3;++d) + atomic_add(exec_tag,&surf_verts_buffer("nrm",d,tri[i]),tnrm[d]); }); - if(!verts.hasProperty("nrm")) - verts.append_channels(cudaPol,{{"nrm",3}}); - TILEVEC_OPS::fill(cudaPol,verts,"nrm",(T)0.0); - cudaPol(zs::range(tris.size()),[ - tris = proxy({},tris), - verts = proxy({},verts)] ZS_LAMBDA(int ti) mutable { - auto tri = tris.pack(dim_c<3>,"inds",ti).reinterpret_bits(int_c); - auto nrm = tris.pack(dim_c<3>,"nrm",ti); - for(int i = 0;i != 3;++i) - for(int d = 0;d != 3;++d) - atomic_add(exec_cuda,&verts("nrm",d,tri[i]),nrm[d]/*/(T)kverts("valence",ktri[i])*/); - }); - cudaPol(zs::range(verts.size()),[verts = proxy({},verts)] ZS_LAMBDA(int vi) mutable { - auto nrm = verts.pack(dim_c<3>,"nrm",vi); - nrm = nrm / (nrm.norm() + (T)1e-6); - verts.tuple(dim_c<3>,"nrm",vi) = nrm; - }); + TILEVEC_OPS::normalized_channel<3>(cudaPol,surf_verts_buffer,"nrm"); + + TILEVEC_OPS::copy(cudaPol,kverts,kxtag,kverts_buffer,"x"); + TILEVEC_OPS::copy(cudaPol,ktris,"inds",ktris_buffer,"inds"); + TILEVEC_OPS::copy(cudaPol,ktris,"he_inds",ktris_buffer,"he_inds"); + if(kverts.hasProperty("X")) + TILEVEC_OPS::copy(cudaPol,kverts,"X",kverts_buffer,"X"); + else + TILEVEC_OPS::copy(cudaPol,kverts,kxtag,kverts_buffer,"X"); + if(kverts.hasProperty("k_active")) + TILEVEC_OPS::copy(cudaPol,kverts,"k_active",kverts_buffer,"k_active"); + else + TILEVEC_OPS::fill(cudaPol,kverts_buffer,"k_active",(T)1.0); + + TILEVEC_OPS::fill(cudaPol,kverts_buffer,"nrm",(T)0.0); + cudaPol(zs::range(ktris_buffer.size()),[exec_tag, + kverts_buffer = proxy({},kverts_buffer), + ktris_buffer = proxy({},ktris_buffer)] ZS_LAMBDA(int kti) mutable { + auto ktri = ktris_buffer.pack(dim_c<3>,"inds",kti,int_c); + zs::vec ktV[3] = {}; + for(int i = 0;i != 3;++i) + ktV[i] = kverts_buffer.pack(dim_c<3>,"x",ktri[i]); + auto ktnrm = LSL_GEO::facet_normal(ktV[0],ktV[1],ktV[2]); + ktris_buffer.tuple(dim_c<3>,"nrm",kti) = ktnrm; + for(int i = 0;i != 3;++i) + for(int d = 0;d != 3;++d) + atomic_add(exec_tag,&kverts_buffer("nrm",d,ktri[i]),ktnrm[d]); + }); + TILEVEC_OPS::normalized_channel<3>(cudaPol,kverts_buffer,"nrm"); + + zs::Vector gia_res{surf_verts_buffer.get_allocator(),0}; + zs::Vector tris_gia_res{surf_tris_buffer.get_allocator(),0}; + + auto ring_mask_width = do_global_intersection_analysis_with_connected_manifolds(cudaPol, + surf_verts_buffer,"x",surf_tris_buffer,halfedges,false, + kverts_buffer,"x",ktris_buffer,khalfedges,true, + gia_res,tris_gia_res); + // for each vertex of zsparticles, find a potential closest point on kboundary surface + // add a plane constraint + + auto ktBvh = bvh_t{}; + auto kt_bvs = retrieve_bounding_volumes(cudaPol,kverts_buffer,ktris_buffer,wrapv<3>{},(T)0.0,"x"); + ktBvh.build(cudaPol,kt_bvs); + // auto cnorm = compute_average_edge_length(cudaPol,verts,"x",tris); + // cnorm *= 2; + auto tBvh = bvh_t{}; + auto tbvs = retrieve_bounding_volumes(cudaPol,surf_verts_buffer,surf_tris_buffer,wrapv<3>{},(T)0.0,"x"); + tBvh.build(cudaPol,tbvs); + + auto kinInCollisionEps = get_input2("kinInColEps"); + auto kinOutCollisionEps = get_input2("kinOutColEps"); + auto thickness = kinInCollisionEps + kinOutCollisionEps; auto align_angle_cosin = get_input2("align_angle_cosin"); + auto avg_norm = compute_average_edge_length(cudaPol,verts,xtag,tris); + auto avg_knorm = compute_average_edge_length(cudaPol,kverts,kxtag,ktris); + auto max_avg_norm = avg_norm > avg_knorm ? avg_norm : avg_knorm; + auto fail_distance = max_avg_norm * (T)5; - cudaPol(zs::range(verts.size()),[ + + cudaPol(zs::range(surf_verts_buffer.size()),[ verts = proxy({},verts), + surf_verts_buffer = proxy({},surf_verts_buffer), ktBvh = proxy(ktBvh), - kverts = proxy({},kverts), - ktris = proxy({},ktris), + fail_distance = fail_distance, + kverts_buffer = proxy({},kverts_buffer), + ktris_buffer = proxy({},ktris_buffer), + kt_offset = tris.size(), + kv_offset = points.size(), + points = proxy({},points), khalfedges = proxy({},khalfedges), align_angle_cosin = align_angle_cosin, - project_pos_tag = zs::SmallString(project_pos_tag), - project_nrm_tag = zs::SmallString(project_nrm_tag), project_idx_tag = zs::SmallString(project_idx_tag), - project_bary_tag = zs::SmallString(project_bary_tag), kinInCollisionEps = kinInCollisionEps, + gia_res = proxy(gia_res), + tris_gia_res = proxy(tris_gia_res), + ring_mask_width = ring_mask_width, align_direction = align_direction, kinOutCollisionEps = kinOutCollisionEps, - thickness = thickness] ZS_LAMBDA(int vi) mutable { - if(verts.hasProperty("is_surf")) - if(verts("is_surf",vi) < (T)0.5) - return; - if(verts.hasProperty("k_active"))// static unbind - if(verts("k_active",vi) < (T)0.5) - return; - // if(verts.hasProperty("k_fail"))// dynamic unbind - // if(verts("k_fail",vi) > (T)0.5) - // return; - auto p = verts.pack(dim_c<3>,"x",vi); + thickness = thickness] ZS_LAMBDA(int pi) mutable { + auto vi = zs::reinterpret_bits(points("inds",pi)); + + auto p = surf_verts_buffer.pack(dim_c<3>,"x",pi); + auto Xp = surf_verts_buffer.pack(dim_c<3>,"X",pi); + auto bv = bv_t{get_bounding_box(p - thickness,p + thickness)}; - auto min_dist = std::numeric_limits::infinity(); + auto min_dist = std::numeric_limits::max(); int min_tri_idx = -1; auto min_bary = vec3::zeros(); auto min_collision_eps = (T)0; - auto pnrm = verts.pack(dim_c<3>,"nrm",vi); + auto pnrm = surf_verts_buffer.pack(dim_c<3>,"nrm",pi); + + auto v_k_active = surf_verts_buffer("k_active",pi); + if(v_k_active < (T)0.5) + return; auto process_potential_closest_tris = [&](int kti) { - auto ktri = ktris.pack(dim_c<3>,"inds",kti).reinterpret_bits(int_c); - if(kverts.hasProperty("k_active")) - for(int i = 0;i != 3;++i) - if(kverts("k_active",ktri[i]) < (T)0.5) - return; - auto kv0 = kverts.pack(dim_c<3>,"x",ktri[0]); - auto kv1 = kverts.pack(dim_c<3>,"x",ktri[1]); - auto kv2 = kverts.pack(dim_c<3>,"x",ktri[2]); + auto ktri = ktris_buffer.pack(dim_c<3>,"inds",kti,int_c); + for(int i = 0;i != 3;++i) + if(kverts_buffer("k_active",ktri[i]) < (T)0.5) + return; - vec3 bary{}; - vec3 project_bary{}; - T distance = LSL_GEO::pointTriangleDistance(kv0,kv1,kv2,p,bary,project_bary); - if(distance > min_dist) + vec3 kvs[3] = {}; + for(int i = 0;i != 3;++i) + kvs[i] = kverts_buffer.pack(dim_c<3>,"x",ktri[i]); + + vec3 kXs[3] = {}; + for(int i = 0;i != 3;++i) + kXs[i] = kverts_buffer.pack(dim_c<3>,"X",ktri[i]); + + auto cX = vec3::zeros(); + for(int i = 0;i != 3;++i) + cX += kXs[i] / (T)3.0; + + auto test_fail_distance = (cX - Xp).norm(); + if(test_fail_distance > fail_distance) return; + vec3 bary{}; + T distance = LSL_GEO::pointTriangleDistance(kvs[0],kvs[1],kvs[2],p,bary); - auto seg = p - kv0; - auto knrm = ktris.pack(dim_c<3>,"nrm",kti); + auto seg = p - kvs[0]; + auto knrm = ktris_buffer.pack(dim_c<3>,"nrm",kti); auto dist = seg.dot(knrm); auto collisionEps = dist > 0 ? kinOutCollisionEps : kinInCollisionEps; if(distance > collisionEps) return; - + + // distance = dist > 0 ? distance : -distance; + if(distance > min_dist) + return; auto align = knrm.dot(pnrm); // TO RECOVER - if(align < align_angle_cosin && align_direction && dist < 0){ - // printf("failed of %d %d due to aligh = %f\n",vi,kti,(float)align); + if(align < align_angle_cosin && align_direction && dist < 0) return; - } - if(align > -align_angle_cosin && !align_direction && dist < 0){ - // printf("failed of %d %d due to aligh = %f\n",vi,kti,(float)align); + if(align > -align_angle_cosin && !align_direction && dist < 0) return; - } auto bary_sum = fabs(bary[0]) + fabs(bary[1]) + fabs(bary[2]); - if(bary_sum > 1.1) - return; - else { - auto khi = zs::reinterpret_bits(ktris("he_inds",kti)); + if(bary_sum > 1.1) { + auto khi = zs::reinterpret_bits(ktris_buffer("he_inds",kti)); for(int i = 0;i != 3;++i){ auto rkhi = zs::reinterpret_bits(khalfedges("opposite_he",khi)); - auto nti = zs::reinterpret_bits(khalfedges("to_face",rkhi)); - // auto nti = ntris[i]; auto edge_normal = vec3::zeros(); - if(nti < 0){ + if(rkhi < 0) { edge_normal = knrm; }else { - edge_normal = ktris.pack(dim_c<3>,"nrm",nti) + knrm; + auto nkti = zs::reinterpret_bits(khalfedges("to_face",rkhi)); + edge_normal = ktris_buffer.pack(dim_c<3>,"nrm",nkti) + knrm; edge_normal = edge_normal/(edge_normal.norm() + (T)1e-6); } - auto ke0 = kverts.pack(dim_c<3>,"x",ktri[(i + 0) % 3]); - auto ke1 = kverts.pack(dim_c<3>,"x",ktri[(i + 1) % 3]); + + auto ke0 = kverts_buffer.pack(dim_c<3>,"x",ktri[(i + 0) % 3]); + auto ke1 = kverts_buffer.pack(dim_c<3>,"x",ktri[(i + 1) % 3]); auto ke10 = ke1 - ke0; auto bisector_normal = edge_normal.cross(ke10).normalized(); - seg = p - kverts.pack(dim_c<3>,"x",ktri[(i + 0) % 3]); + seg = p - kverts_buffer.pack(dim_c<3>,"x",ktri[(i + 0) % 3]); if(bisector_normal.dot(seg) < 0) return; - khi = zs::reinterpret_bits(khalfedges("next_he",khi)); } } + int RING_MASK = 0; + for(int i = 0;i != ring_mask_width;++i) { + auto MASK = gia_res[pi * ring_mask_width + i] & tris_gia_res[(kt_offset + kti) * ring_mask_width + i]; + RING_MASK |= MASK; + } - + if(RING_MASK > 0 && dist > 0) + return; + if(RING_MASK == 0 && dist < 0) + return; min_dist = distance; min_tri_idx = kti; - min_bary = project_bary; - min_collision_eps = collisionEps; }; ktBvh.iter_neighbors(bv,process_potential_closest_tris); + verts(project_idx_tag,vi) = reinterpret_bits(min_tri_idx); + }); - if(min_tri_idx == -1) - return; - auto closest_ktri = ktris.pack(dim_c<3>,"inds",min_tri_idx).reinterpret_bits(int_c); - if(kverts.hasProperty("k_fail")) + cudaPol(zs::range(kverts_buffer.size()),[ + surf_verts_buffer = proxy({},surf_verts_buffer), + surf_tris_buffer = proxy({},surf_tris_buffer), + halfedges = proxy({},halfedges), + kv_offset = points.size(), + gia_res = proxy(gia_res), + tris_gia_res = proxy(tris_gia_res), + ring_mask_width = ring_mask_width, + tBvh = proxy(tBvh), + kverts = proxy({},kverts), + kverts_buffer = proxy({},kverts_buffer), + align_angle_cosin = align_angle_cosin, + project_idx_tag = zs::SmallString(project_idx_tag), + kinInCollisionEps = kinInCollisionEps, + kinOutCollisionEps = kinOutCollisionEps, + align_direction = align_direction, + thickness = thickness] ZS_LAMBDA(int kvi) mutable { + auto kp = kverts_buffer.pack(dim_c<3>,"x",kvi); + auto knrm = kverts_buffer.pack(dim_c<3>,"nrm",kvi); + auto min_dist = std::numeric_limits::max(); + int min_tri_idx = -1; + + auto bv = bv_t{get_bounding_box(kp - thickness,kp + thickness)}; + + auto process_potential_closest_tris = [&](int ti) mutable { + auto tri = surf_tris_buffer.pack(dim_c<3>,"inds",ti,int_c); + vec3 tvs[3] = {}; for(int i = 0;i != 3;++i) - if(kverts("k_fail",closest_ktri[i]) > (T)0.5) - return; + tvs[i] = surf_verts_buffer.pack(dim_c<3>,"x",tri[i]); + vec3 bnrms[3] = {}; + auto tnrm = surf_tris_buffer.pack(dim_c<3>,"nrm",ti); + auto hi = zs::reinterpret_bits(surf_tris_buffer("he_inds",ti)); + for(int i = 0;i != 3;++i) { + auto edge_normal = tnrm; + auto rhi = zs::reinterpret_bits(halfedges("opposite_he",hi)); + if(rhi >= 0) { + auto nti = zs::reinterpret_bits(halfedges("to_face",rhi)); + edge_normal = tnrm + surf_tris_buffer.pack(dim_c<3>,"nrm",nti); + edge_normal = edge_normal / (edge_normal.norm() + (T)1e-6); + } - // auto ori_bary = min_bary; - - min_bary[0] = min_bary[0] < 0 ? (T)0 : min_bary[0]; - min_bary[1] = min_bary[1] < 0 ? (T)0 : min_bary[1]; - min_bary[2] = min_bary[2] < 0 ? (T)0 : min_bary[2]; - min_bary = min_bary/min_bary.sum(); - if(min_bary[0] < 0 || min_bary[1] < 0 || min_bary[2] < 0) - printf("invalid min_bary[%f %f %f]\n", - (float)min_bary[0], - (float)min_bary[1], - (float)min_bary[2]); - if((zs::abs(min_bary[0]) + zs::abs(min_bary[1]) + zs::abs(min_bary[2])) > 1.1) - printf("invalid min_bary[%f %f %f]\n", - (float)min_bary[0], - (float)min_bary[1], - (float)min_bary[2]); - - if((zs::abs(min_bary[0] + min_bary[1] + min_bary[2] - 1.0)) > 0.1) - printf("invalid min_bary[%f %f %f]\n", - (float)min_bary[0], - (float)min_bary[1], - (float)min_bary[2]); + auto e01 = tvs[(i + 1) % 3] - tvs[(i + 0) % 3]; + bnrms[i] = edge_normal.cross(e01).normalized(); + hi = zs::reinterpret_bits(halfedges("next_he",hi)); + } - auto project_kv = vec3::zeros(); - for(int i = 0;i != 3;++i) - project_kv += kverts.pack(dim_c<3>,"x",closest_ktri[i]) * min_bary[i]; - auto project_knrm = vec3::zeros(); - for(int i = 0;i != 3;++i) - project_knrm += kverts.pack(dim_c<3>,"nrm",closest_ktri[i]) * min_bary[i]; - project_knrm /= (project_knrm.norm() + 1e-6); + vec3 bary{}; + T distance = LSL_GEO::pointTriangleDistance(tvs[0],tvs[1],tvs[2],kp,bary); - // printf("find closest pairs : %d %d\n",vi,min_tri_idx); + auto seg = kp - tvs[0]; + auto dist = -seg.dot(tnrm); - // printf("vert[%d] = %f %f %f closest to ktri[%d] = %f %f %f bary = %f %f %f\n",vi, - // (float)p[0],(float)p[1],(float)p[2], - // min_tri_idx, - // (float)project_kv[0],(float)project_kv[1],(float)project_kv[2], - // (float)ori_bary[0],(float)ori_bary[1],(float)ori_bary[2] - // ); + auto collisionEps = dist > 0 ? kinOutCollisionEps : kinInCollisionEps; + if(distance > collisionEps) + return; + + // distance = dist > 0 ? distance : -distance; + if(distance > min_dist) + return; - verts.tuple(dim_c<3>,project_pos_tag,vi) = project_kv; + auto align = knrm.dot(tnrm); + // TO RECOVER + if(align < align_angle_cosin && align_direction && dist < 0){ + // printf("failed of %d %d due to aligh = %f\n",vi,kti,(float)align); + return; + } + if(align > -align_angle_cosin && !align_direction && dist < 0){ + // printf("failed of %d %d due to aligh = %f\n",vi,kti,(float)align); + return; + } - auto distance = (verts.pack(dim_c<3>,"x",vi) - project_kv).norm(); - if(distance > kinInCollisionEps) - printf("find invalid distance %f > %f : %f\n",(float)distance,(float)kinInCollisionEps,(float)min_collision_eps); + auto bary_sum = fabs(bary[0]) + fabs(bary[1]) + fabs(bary[2]); + if(bary_sum > 1.1) { + for(int i = 0;i != 3;++i){ + seg = kp - surf_verts_buffer.pack(dim_c<3>,"x",tri[(i + 0) % 3]); + if(bnrms[i].dot(seg) < 0) + return; + } + } - verts.tuple(dim_c<3>,project_nrm_tag,vi) = project_knrm; - verts.tuple(dim_c<3>,project_bary_tag,vi) = min_bary; - verts(project_idx_tag,vi) = reinterpret_bits(min_tri_idx); + int RING_MASK = 0; + for(int i = 0;i != ring_mask_width;++i) { + auto MASK = gia_res[(kvi + kv_offset) * ring_mask_width + i] & tris_gia_res[ti * ring_mask_width + i]; + RING_MASK |= MASK; + } + if(RING_MASK > 0 && dist > 0) + return; + if(RING_MASK == 0 && dist < 0) + return; + + min_dist = distance; + min_tri_idx = ti; + }; + tBvh.iter_neighbors(bv,process_potential_closest_tris); + kverts(project_idx_tag,kvi) = reinterpret_bits(min_tri_idx); }); set_output("zsparticles",zsparticles); @@ -1038,10 +1167,9 @@ ZENDEFNODE(ZSSurfaceClosestPoints, { {"zsparticles","kboundary"}, { {"bool","align_direction","1"}, - {"string","project_pos_tag","project_pos_tag"}, - {"string","project_nrm_tag","project_nrm_tag"}, {"string","project_idx_tag","project_idx_tag"}, - {"string","project_bary_tag","project_bary_tag"} + {"string","xtag","x"}, + {"string","kxtag","x"} }, {"ZSGeometry"}}); @@ -1415,15 +1543,16 @@ struct ZSSurfaceClosestTris : zeno::INode { auto hi = zs::reinterpret_bits(tris("he_inds",ti)); vec3 bnrms[3] = {}; for(int i = 0;i != 3;++i){ + auto edge_normal = tnrm; auto rhi = zs::reinterpret_bits(halfedges("opposite_he",hi)); - auto nti = zs::reinterpret_bits(halfedges("to_face",rhi)); - auto edge_normal = vec3::zeros(); - if(nti < 0) + if(rhi < 0) edge_normal = tnrm; else{ + auto nti = zs::reinterpret_bits(halfedges("to_face",rhi)); edge_normal = tnrm + tris.pack(dim_c<3>,"nrm",nti); edge_normal = edge_normal/(edge_normal.norm() + (T)1e-6); } + auto e01 = tvs[(i + 1) % 3] - tvs[(i + 0) % 3]; bnrms[i] = edge_normal.cross(e01).normalized(); hi = zs::reinterpret_bits(halfedges("next_he",hi)); @@ -1533,144 +1662,117 @@ struct VisualizeClosestTris : zeno::INode { // auto project_nrm_tag = get_param("project_nrm_tag"); // the id of the vertex on kboundary auto project_idx_tag = get_param("project_idx_tag"); - auto kinOutCollisionEps = get_input2("kin_out_collision_eps"); + // auto kinOutCollisionEps = get_input2("kin_out_collision_eps"); - dtiles_t verts_buffer{tris.get_allocator(),{ + dtiles_t verts_buffer{verts.get_allocator(),{ {"x",3}, - {"xp",3}, - {"nrm",3}, - {"inds",3}, - {"grad",9} - },tris.size()}; + {"xp",3} + // {"nrm",3}, + // {"inds",3}, + // {"grad",9} + },verts.size()}; - dtiles_t force_buffer{verts.get_allocator(),{ - {"force",3}, - {"x",3} - },verts.size()}; + dtiles_t kverts_buffer{kverts.get_allocator(),{ + {"x",3}, + {"xp",3} + },kverts.size()}; - TILEVEC_OPS::copy(cudaPol,tris,"inds",verts_buffer,"inds"); - TILEVEC_OPS::copy(cudaPol,verts,"x",force_buffer,"x"); - TILEVEC_OPS::fill(cudaPol,verts_buffer,"grad",(T)0.0); - TILEVEC_OPS::fill(cudaPol,force_buffer,"force",(T)0.0); - cudaPol(zs::range(tris.size()),[ + // TILEVEC_OPS::copy(cudaPol,tris,"inds",verts_buffer,"inds"); + // TILEVEC_OPS::copy(cudaPol,verts,"x",force_buffer,"x"); + // TILEVEC_OPS::fill(cudaPol,verts_buffer,"grad",(T)0.0); + // TILEVEC_OPS::fill(cudaPol,force_buffer,"force",(T)0.0); + + TILEVEC_OPS::copy(cudaPol,verts,"x",verts_buffer,"x"); + TILEVEC_OPS::copy(cudaPol,verts,"x",verts_buffer,"xp"); + TILEVEC_OPS::copy(cudaPol,kverts,"x",kverts_buffer,"x"); + TILEVEC_OPS::copy(cudaPol,kverts,"x",kverts_buffer,"xp"); + + std::cout << "evaluate verts_buffer" << std::endl; + + cudaPol(zs::range(verts.size()),[ verts_buffer = proxy({},verts_buffer), verts = proxy({},verts), - tris = proxy({},tris), kverts = proxy({},kverts), ktris = proxy({},ktris), - kinOutCollisionEps = kinOutCollisionEps, - project_idx_tag = zs::SmallString(project_idx_tag)] ZS_LAMBDA(int ti) mutable { - auto tri = tris.pack(dim_c<3>,"inds",ti).reinterpret_bits(int_c); - auto tp = vec3::zeros(); + // kinOutCollisionEps = kinOutCollisionEps, + project_idx_tag = zs::SmallString(project_idx_tag)] ZS_LAMBDA(int vi) mutable { + auto kti = zs::reinterpret_bits(verts(project_idx_tag,vi)); + if(kti < 0) + return; + auto ktri = ktris.pack(dim_c<3>,"inds",kti).reinterpret_bits(int_c); + auto ktc = vec3::zeros(); for(int i = 0;i != 3;++i) - tp += verts.pack(dim_c<3>,"x",tri[i]) / (T)3.0; + ktc += kverts.pack(dim_c<3>,"x",ktri[i]) / (T)3.0; - verts_buffer.tuple(dim_c<3>,"x",ti) = tp; - verts_buffer.tuple(dim_c<3>,"nrm",ti) = vec3::zeros(); - auto kp_idx = reinterpret_bits(tris(project_idx_tag,ti)); - if(kp_idx < 0) - verts_buffer.tuple(dim_c<3>,"xp",ti) = tp; - else{ - auto kp = kverts.pack(dim_c<3>,"x",kp_idx); - auto tnrm = tris.pack(dim_c<3>,"nrm",ti); - verts_buffer.tuple(dim_c<3>,"xp",ti) = kp; - verts_buffer.tuple(dim_c<3>,"nrm",ti) = tnrm; + verts_buffer.tuple(dim_c<3>,"xp",vi) = ktc; + }); - vec3 vs[4] = {}; - auto nrm = tris.pack(dim_c<3>,"nrm",ti); - vs[0] = kp; - for(int i = 0;i != 3;++i) - vs[i + 1] = verts.pack(dim_c<3>,"x",tri[i]); - - vec3 e[3] = {}; - e[0] = vs[3] - vs[2]; - e[1] = vs[0] - vs[2]; - e[2] = vs[1] - vs[2]; - - auto n = e[2].cross(e[0]); - n = n/(n.norm() + 1e-6); - - T springLength = e[1].dot(n) - kinOutCollisionEps; - auto gvf = zs::vec::zeros(); - if(springLength < (T)0){ - auto gvf_v12 = COLLISION_UTILS::springLengthGradient(vs,e,n); - if(isnan(gvf_v12.norm())) - printf("nan gvf detected at %d %f %f\n",ti,gvf_v12.norm(),n.norm()); - for(int i = 0;i != 9;++i) - gvf[i] = gvf_v12[i + 3]; - } - auto stiffness = (T)1.0; - auto g = -stiffness * springLength * gvf; - // auto H = stiffness * zs::dyadic_prod(gvf, gvf); - verts_buffer.tuple(dim_c<9>,"grad",ti) = g; - } + std::cout << "evaluate kverts_buffer" << std::endl; + cudaPol(zs::range(kverts.size()),[ + kverts_buffer = proxy({},kverts_buffer), + kverts = proxy({},kverts), + verts = proxy({},verts), + tris = proxy({},tris), + // kinOutCollisionEps = kinOutCollisionEps, + project_idx_tag = zs::SmallString(project_idx_tag)] ZS_LAMBDA(int kvi) mutable { + auto ti = zs::reinterpret_bits(kverts(project_idx_tag,kvi)); + if(ti < 0) + return; + auto tri = tris.pack(dim_c<3>,"inds",ti).reinterpret_bits(int_c); + auto tc = vec3::zeros(); + for(int i = 0;i != 3;++i) + tc += verts.pack(dim_c<3>,"x",tri[i]) / (T)3.0; + + kverts_buffer.tuple(dim_c<3>,"xp",kvi) = tc; }); - TILEVEC_OPS::assemble(cudaPol,verts_buffer,"grad","inds",force_buffer,"force"); - constexpr auto omp_space = execspace_e::openmp; auto ompPol = omp_exec(); verts_buffer = verts_buffer.clone({zs::memsrc_e::host}); - force_buffer = force_buffer.clone({zs::memsrc_e::host}); - auto closest_points_vis = std::make_shared(); - auto& pverts = closest_points_vis->verts; - auto& plines = closest_points_vis->lines; - pverts.resize(verts_buffer.size() * 2); - plines.resize(verts_buffer.size()); + kverts_buffer = kverts_buffer.clone({zs::memsrc_e::host}); - ompPol(zs::range(verts_buffer.size()),[ - verts_buffer = proxy({},verts_buffer), - &pverts,&plines] (int vi) mutable { - pverts[vi * 2 + 0] = verts_buffer.pack(dim_c<3>,"x",vi).to_array(); - pverts[vi * 2 + 1] = verts_buffer.pack(dim_c<3>,"xp",vi).to_array(); - plines[vi] = zeno::vec2i{vi * 2 + 0,vi * 2 + 1}; - }); + std::cout << "evaluate closest_ktris_vis" << std::endl; + + auto closest_ktris_vis = std::make_shared(); + auto& c_verts = closest_ktris_vis->verts; + auto& c_lines = closest_ktris_vis->lines; + c_verts.resize(verts_buffer.size() * 2); + c_lines.resize(verts_buffer.size()); - auto nrm_scale = get_input2("nrm_scale"); - auto normal_vis = std::make_shared(); - auto& nverts = normal_vis->verts; - auto& nlines = normal_vis->lines; - nverts.resize(verts_buffer.size() * 2); - nlines.resize(verts_buffer.size()); ompPol(zs::range(verts_buffer.size()),[ verts_buffer = proxy({},verts_buffer), - &nverts,&nlines,nrm_scale = nrm_scale] (int vi) mutable { - nverts[vi * 2 + 0] = verts_buffer.pack(dim_c<3>,"x",vi).to_array(); - auto ep = verts_buffer.pack(dim_c<3>,"nrm",vi) * nrm_scale + verts_buffer.pack(dim_c<3>,"x",vi); - nverts[vi * 2 + 1] = ep.to_array(); - nlines[vi] = zeno::vec2i{vi * 2 + 0,vi * 2 + 1}; + &c_verts,&c_lines] (int vi) mutable { + c_verts[vi * 2 + 0] = verts_buffer.pack(dim_c<3>,"x",vi).to_array(); + c_verts[vi * 2 + 1] = verts_buffer.pack(dim_c<3>,"xp",vi).to_array(); + c_lines[vi] = zeno::vec2i{vi * 2 + 0,vi * 2 + 1}; }); - auto force_scale = get_input2("force_scale"); - auto force_vis = std::make_shared(); - auto& fverts = force_vis->verts; - auto& flines = force_vis->lines; - fverts.resize(2 * verts.size()); - flines.resize(verts.size()); - ompPol(zs::range(verts.size()),[ - force_buffer = proxy({},force_buffer), - &fverts,&flines,force_scale = force_scale] (int vi) mutable { - fverts[vi * 2 + 0] = force_buffer.pack(dim_c<3>,"x",vi).to_array(); - auto ep = force_buffer.pack(dim_c<3>,"force",vi) * force_scale + force_buffer.pack(dim_c<3>,"x",vi); - fverts[vi * 2 + 1] = ep.to_array(); - flines[vi] = zeno::vec2i{vi * 2 + 0,vi * 2 + 1}; + std::cout << "evaluate closest_tris_vis" << std::endl; + + auto closest_tris_vis = std::make_shared(); + auto& kc_verts = closest_tris_vis->verts; + auto& kc_lines = closest_tris_vis->lines; + kc_verts.resize(kverts_buffer.size() * 2); + kc_lines.resize(kverts_buffer.size()); + + ompPol(zs::range(kverts_buffer.size()),[ + kverts_buffer = proxy({},kverts_buffer), + &kc_verts,&kc_lines] (int kvi) mutable { + kc_verts[kvi * 2 + 0] = kverts_buffer.pack(dim_c<3>,"x",kvi).to_array(); + kc_verts[kvi * 2 + 1] = kverts_buffer.pack(dim_c<3>,"xp",kvi).to_array(); + kc_lines[kvi] = zeno::vec2i{kvi * 2 + 0,kvi * 2 + 1}; }); - set_output("closest_vis",std::move(closest_points_vis)); - set_output("normal_vis",std::move(normal_vis)); - set_output("force_vis",std::move(force_vis)); + set_output("closest_v2kt_vis",std::move(closest_ktris_vis)); + set_output("closest_kv2t_vis",std::move(closest_tris_vis)); } }; -ZENDEFNODE(VisualizeClosestTris, {{"zsparticles","kboundary", - {"float","nrm_scale","1.0"}, - {"float","force_scale","1.0"}, - {"float","kin_out_collision_eps","0.001"} - }, - {"closest_vis","normal_vis","force_vis"}, +ZENDEFNODE(VisualizeClosestTris, {{"zsparticles","kboundary"}, + {"closest_v2kt_vis","closest_kv2t_vis"}, { {"string","project_idx_tag","project_idx_tag"}, }, @@ -2022,14 +2124,18 @@ struct ZSSurfaceClosestPointsGrp : zeno::INode { auto hi = zs::reinterpret_bits(ktris("he_inds",kti)); for(int i = 0;i != 3;++i){ auto rhi = zs::reinterpret_bits(khalfedges("opposite_he",hi)); - auto nti = zs::reinterpret_bits(khalfedges("to_face",rhi)); auto edge_normal = vec3::zeros(); - if(nti < 0){ + + if(rhi < 0){ edge_normal = knrm; }else { + auto nti = zs::reinterpret_bits(khalfedges("to_face",rhi)); edge_normal = ktris.pack(dim_c<3>,"nrm",nti) + knrm; edge_normal = edge_normal/(edge_normal.norm() + (T)1e-6); } + + + auto ke0 = kverts.pack(dim_c<3>,"x",ktri[(i + 0) % 3]); auto ke1 = kverts.pack(dim_c<3>,"x",ktri[(i + 1) % 3]); auto ke10 = ke1 - ke0; @@ -2189,13 +2295,15 @@ struct ZSSurfaceClosestPointsGrp : zeno::INode { auto khi = zs::reinterpret_bits(ktris("he_inds",kti)); for(int i = 0;i != 3;++i){ + auto edge_normal = vec3::zeros(); auto rkhi = zs::reinterpret_bits(khalfedges("opposite_he",khi)); - auto nti = zs::reinterpret_bits(khalfedges("to_face",rkhi)); + // auto nti = ntris[i]; - auto edge_normal = vec3::zeros(); - if(nti < 0){ + + if(rkhi < 0){ edge_normal = knrm; }else { + auto nti = zs::reinterpret_bits(khalfedges("to_face",rkhi)); edge_normal = ktris.pack(dim_c<3>,"nrm",nti) + knrm; edge_normal = edge_normal/(edge_normal.norm() + (T)1e-6); } diff --git a/projects/CuLagrange/geometry/VectorField.cu b/projects/CuLagrange/geometry/VectorField.cu index e633c3a2bd..228d50947f 100644 --- a/projects/CuLagrange/geometry/VectorField.cu +++ b/projects/CuLagrange/geometry/VectorField.cu @@ -167,14 +167,14 @@ struct ZSRetrieveVectorField : zeno::INode { bool on_elm = (type == "quad" || type == "tri"); - if((type == "quad" || type == "tri") && (!eles.hasProperty(gtag) || !eles.hasProperty(color_tag))){ + if((type == "quad" || type == "tri") && (!eles.hasProperty(gtag))){ if(!eles.hasProperty(gtag)) fmt::print("the elements does not contain element-wise gradient field : {}\n",gtag); - if(!eles.hasProperty(color_tag)) - fmt::print("the elements does not contain element-wise color_tag field : {}\n",color_tag); + // if(!eles.hasProperty(color_tag)) + // fmt::print("the elements does not contain element-wise color_tag field : {}\n",color_tag); throw std::runtime_error("the volume does not contain element-wise gradient field"); } - if(type == "vert" && !verts.hasProperty(gtag) && !verts.hasProperty(color_tag)){ + if(type == "vert" && !verts.hasProperty(gtag)){ fmt::print("the volume does not contain nodal-wize gradient field : {}\n",gtag); throw std::runtime_error("the volume does not contain nodal-wize gradient field"); } @@ -210,13 +210,18 @@ struct ZSRetrieveVectorField : zeno::INode { } vec_buffer.tuple<3>("x",i) = bx; vec_buffer.tuple<3>("vec",i) = scale * eles.pack<3>(gtag,i)/* / eles.pack<3>(gtag,i).norm()*/; - zsvec_buffer[i] = eles(color_tag,i); - // vec_buffer(color_tag,i) = eles(color_tag,i); + if(eles.hasProperty(color_tag)) + zsvec_buffer[i] = eles(color_tag,i); + else + zsvec_buffer[i] = (T)1.0; }else{ vec_buffer.tuple<3>("x",i) = verts.pack<3>(xtag,i); vec_buffer.tuple<3>("vec",i) = scale * verts.pack<3>(gtag,i)/* / verts.pack<3>(gtag,i).norm()*/; // vec_buffer(color_tag,i) = verts(color_tag,i); - zsvec_buffer[i] = verts(color_tag,i); + if(verts.hasProperty(color_tag)) + zsvec_buffer[i] = verts(color_tag,i); + else + zsvec_buffer[i] = (T)1.0; } }); @@ -315,26 +320,30 @@ struct ZSSampleQuadratureAttr2Vert : zeno::INode { fmt::print("the verts' {} attr[{}] and quads' {} attr[{}] not matched\n",attr,verts.getPropertySize(attr),attr,attr_dim); } cudaPol(range(verts.size()), - [verts = proxy({},verts),attr_dim,attr = SmallString(attr)] + [skip_bou,bou_tag = zs::SmallString(bou_tag),verts = proxy({},verts),attr_dim,attr = SmallString(attr)] __device__(int vi) mutable { + if(skip_bou && verts(bou_tag,vi) > 1e-6) + return; for(int i = 0;i != attr_dim;++i) verts(attr,i,vi) = 0.; }); - static dtiles_t vtemp(verts.get_allocator(),{{"wsum",1}},verts.size()); - vtemp.resize(verts.size()); - cudaPol(range(vtemp.size()), - [vtemp = proxy({},vtemp)] ZS_LAMBDA(int vi) mutable { - vtemp("wsum",vi) = 0; - }); + // static dtiles_t vtemp(verts.get_allocator(),{{"wsum",1}},verts.size()); + // vtemp.resize(verts.size()); + // cudaPol(range(vtemp.size()), + // [vtemp = proxy({},vtemp)] ZS_LAMBDA(int vi) mutable { + // vtemp("wsum",vi) = 0; + // }); + zs::Vector wsum_vec{verts.get_allocator(),verts.size()}; + cudaPol(zs::range(wsum_vec),[] ZS_LAMBDA(auto& wsum){wsum = (float)0;}); // std::cout << "check here 2" << std::endl; cudaPol(range(quads.size()), [verts = proxy({},verts),quads = proxy({},quads),attr_dim,attr = SmallString(attr),simplex_size,weight = SmallString(weight), - execTag = wrapv{},skip_bou,bou_tag = zs::SmallString(bou_tag),vtemp = proxy({},vtemp)] + execTag = wrapv{},skip_bou = skip_bou,bou_tag = zs::SmallString(bou_tag),wsum_vec = proxy(wsum_vec)] __device__(int ei) mutable { - float w = quads(weight,ei); + auto w = quads(weight,ei); // if(ei == 0) // printf("w : %f\n",(float)w); // w = 1.0;// cancel out the specified weight info @@ -348,23 +357,28 @@ struct ZSSampleQuadratureAttr2Vert : zeno::INode { // verts(attr,j,idx) += w * quads(attr,j,ei) / (float)simplex_size; atomic_add(execTag,&verts(attr,j,idx),alpha * quads(attr,j,ei)); } - atomic_add(execTag,&vtemp("wsum",idx),alpha); + atomic_add(execTag,&wsum_vec[idx],alpha); } }); // std::cout << "check here 3 aaaa" << std::endl; // std::cout << "attr_dim = " << attr_dim << std::endl; - cudaPol(range(verts.size()), - [ + cudaPol(range(verts.size()),[skip_bou,bou_tag = zs::SmallString(bou_tag),simplex_size, verts = proxy({},verts),attr = SmallString(attr), - attr_dim,vtemp = proxy({},vtemp)] ZS_LAMBDA(int vi) mutable { + attr_dim,wsum_vec = proxy(wsum_vec)] ZS_LAMBDA(int vi) mutable { // if(vi == 0) // printf("wsum : %f\n",(float)vtemp("wsum",vi)); + if(skip_bou && verts(bou_tag,vi) > 1e-6) + return; for(int j = 0;j != attr_dim;++j) { // verts(attr,j,idx) += w * quads(attr,j,ei) / (float)simplex_size; - verts(attr,j,vi) = verts(attr,j,vi) / vtemp("wsum",vi); + verts(attr,j,vi) = verts(attr,j,vi) / (wsum_vec[vi] + (float)1e-6); + if(zs::isnan(verts(attr,j,vi) )) + printf("nan verts attr %s detected at %d %d\n",attr.asChars(),(int)j,(int)vi); } + + }); // std::cout << "check here 4" << std::endl; diff --git a/projects/CuLagrange/geometry/file_parser/read_vtk_mesh.hpp b/projects/CuLagrange/geometry/file_parser/read_vtk_mesh.hpp index 0e6f6f2d57..00634c7f80 100644 --- a/projects/CuLagrange/geometry/file_parser/read_vtk_mesh.hpp +++ b/projects/CuLagrange/geometry/file_parser/read_vtk_mesh.hpp @@ -126,9 +126,12 @@ namespace zeno { if(j != 2 || nm_points_read != (numberofpoints - 1)) bufferp = find_next_numeric(bufferp,buffer,fp,&line_count); + } - // printf("\n"); + // printf("points[%d] at line[%d] : (%f %f %f)\n",nm_points_read,line_count,(float)verts[nm_points_read][0],(float)verts[nm_points_read][1],(float)verts[nm_points_read][2]); nm_points_read++; + // printf("\n"); + } return true; } @@ -143,7 +146,7 @@ namespace zeno { int nm_cells_read = 0; - // printf("numberofcells : %d\n",numberofcells); + printf("numberofcells : %d\n",numberofcells); while(nm_cells_read < numberofcells){ bufferp = readline(buffer,fp,&line_count); @@ -592,19 +595,21 @@ namespace zeno { continue; } if(!strcmp(id,"POINTS")){ - printf("reading points\n"); + printf("reading points %d\n",line_count); int numberofpoints = 0; sscanf(line,"%s %d %s",id,&numberofpoints,dummy_str); printf("number of points %d\n",numberofpoints); parsing_verts_coord(fp,prim->verts,numberofpoints,line_count); + printf("finish reading points %d\n",line_count); continue; } if(!strcmp(id,"CELLS")){ - printf("reading cells\n"); + printf("reading cells %d\n",line_count); int numberofcells = 0; int numberofdofs = 0; sscanf(line,"%s %d %d",id,&numberofcells,&numberofdofs); simplex_size = numberofdofs/numberofcells - 1; + printf("simplex_size %d\n",simplex_size); if(simplex_size == 4) parsing_cells_topology<4>(fp,prim->quads,numberofcells,line_count); else if(simplex_size == 3) @@ -616,25 +621,28 @@ namespace zeno { continue; } if(!strcmp(id,"CELL_TYPES")){ - printf("reading cell types\n"); + printf("reading cell types %d\n",line_count); int numberofcells = 0; sscanf(line,"%s %d",id,&numberofcells); printf("number of cell types : %d\n",numberofcells); bufferp = readline(line,fp,&line_count); if(numberofcells > 0){ - int type = strtol(bufferp,&bufferp,0); - if(type != VTK_TETRA && simplex_size == 4){ - printf("non-tetra cell detected on line %d parsing cell types with simplex size = 4\n",line_count); - fclose(fp); - return false; - }else if(type != VTK_TRIANGLE && simplex_size == 3) { - printf("non-triangle cell detected on line %d parsing cell types with simplex size = 3\n",line_count); - fclose(fp); - return false; + for(int i = 0;i != numberofcells;++i) { + int type = strtol(bufferp,&bufferp,0); + if(type != VTK_TETRA && simplex_size == 4){ + printf("non-tetra cell detected on line %d parsing cell types with simplex size = 4\n",line_count); + fclose(fp); + return false; + }else if(type != VTK_TRIANGLE && simplex_size == 3) { + printf("non-triangle cell detected on line %d parsing cell types with simplex size = 3\n",line_count); + fclose(fp); + return false; + } + if(i+1 != numberofcells) + bufferp = find_next_numeric(bufferp,line,fp,&line_count); } - bufferp = find_next_numeric(bufferp,line,fp,&line_count); } - printf("finish cell type check\n"); + printf("finish cell type check at %d\n",line_count); continue; } diff --git a/projects/CuLagrange/geometry/kernel/intersection.hpp b/projects/CuLagrange/geometry/kernel/intersection.hpp index f52d229f6b..0cb88d5ebb 100644 --- a/projects/CuLagrange/geometry/kernel/intersection.hpp +++ b/projects/CuLagrange/geometry/kernel/intersection.hpp @@ -574,9 +574,12 @@ size_t retrieve_self_intersection_tri_halfedge_list_info(Pol& pol, auto cnorm = compute_average_edge_length(pol,verts,xtag,tris); cnorm *= 3; + auto max_intersections = intersect_buffers.size(); + pol(zs::range(halfedges.size()),[ exec_tag, nmIts = proxy(nmIts), + max_intersections = max_intersections, halfedges = proxy({},halfedges),/*'to_vertex' 'to_face' 'opposite_he' 'next_he'*/ verts = proxy({},verts), nm_verts = verts.size(), @@ -641,10 +644,12 @@ size_t retrieve_self_intersection_tri_halfedge_list_info(Pol& pol, // LSL_GEO::tri_ray_intersect_d(eV[0],eV[1],tV[0],tV[1],tV[2],r); if(LSL_GEO::tri_ray_intersect_d(eV[0],eV[1],tV[0],tV[1],tV[2],r)) { auto offset = atomic_add(exec_tag,&nmIts[0],(int)1); + if(offset >= max_intersections) + return; auto intp = r * dir + eV[0]; intersect_buffers.tuple(dim_c<2>,"pair",offset) = zs::vec{hei,ti}.reinterpret_bits(float_c); intersect_buffers.tuple(dim_c<3>,"int_points",offset) = intp; - intersect_buffers("r",offset) = r; + intersect_buffers("r",offset) = (T)r; // make sure the opposite he - tri pairs are also inserted // auto opposite_hei = zs::reinterpret_bits(halfedges("opposite_he",hei)); @@ -662,6 +667,10 @@ size_t retrieve_self_intersection_tri_halfedge_list_info(Pol& pol, // std::cout << "initialize corner_idx : " << nmIts.getVal(0) << std::endl; + if(nmIts.getVal(0) >= max_intersections) { + throw std::runtime_error("max_size_of_intersections buffer reach"); + } + pol(zs::range(nmIts.getVal(0)),[ intersect_buffers = proxy({},intersect_buffers), tris = proxy({},tris), @@ -702,11 +711,11 @@ int do_global_self_intersection_analysis(Pol& pol, const PosTileVec& verts, const zs::SmallString& xtag, const TriTileVec& tris, - const HalfEdgeTileVec& halfedges, + HalfEdgeTileVec& halfedges, GIA_TILEVEC& gia_res, GIA_TILEVEC& tris_gia_res, // zs::bht& conn_of_first_ring, - size_t max_nm_intersections = 10000) { + size_t max_nm_intersections = 50000) { using namespace zs; using T = typename PosTileVec::value_type; using index_type = std::make_signed_t; @@ -733,9 +742,9 @@ int do_global_self_intersection_analysis(Pol& pol, {"r",1}, },max_nm_intersections}; - // if(!halfedges.hasProperty("broken")) - // halfedges.append_channels(pol,{{"broken",1}}); - // TILEVEC_OPS::fill(pol,halfedges,"broken",(T)0.0); + if(!halfedges.hasProperty("broken")) + halfedges.append_channels(pol,{{"broken",1}}); + TILEVEC_OPS::fill(pol,halfedges,"broken",(T)0.0); auto nm_insts = retrieve_self_intersection_tri_halfedge_list_info(pol,verts,xtag,tris,halfedges,ints_buffer); table_vec2i_type cftab{ints_buffer.get_allocator(),(size_t)nm_insts}; @@ -862,7 +871,26 @@ int do_global_self_intersection_analysis(Pol& pol, auto nm_rings = mark_disconnected_island(pol,conn_topo,ringTag); - std::cout << "finish Mark disconnected island with nm_rings : " << nm_rings << std::endl; + // std::cout << "finish Mark disconnected island with nm_rings : " << nm_rings << std::endl; + + auto ring_mask_width = (nm_rings + 31) / 32; + + gia_res.resize(verts.size() * ring_mask_width); + pol(zs::range(gia_res.size()),[gia_res = proxy({},gia_res)] ZS_LAMBDA(int mi) mutable { + // nodal_colors[ni] = 0; + gia_res("ring_mask",mi) = zs::reinterpret_bits((int)0); + gia_res("color_mask",mi) = zs::reinterpret_bits((int)0); + gia_res("type_mask",mi) = zs::reinterpret_bits((int)0); + gia_res("is_loop_vertex",mi) = zs::reinterpret_bits((int)0); + }); + tris_gia_res.resize(tris.size() * ring_mask_width); + pol(zs::range(tris_gia_res.size()),[tris_gia_res = proxy({},tris_gia_res)] ZS_LAMBDA(int mi) mutable { + // nodal_colors[ni] = 0; + tris_gia_res("ring_mask",mi) = zs::reinterpret_bits((int)0); + tris_gia_res("color_mask",mi) = zs::reinterpret_bits((int)0); + tris_gia_res("type_mask",mi) = zs::reinterpret_bits((int)0); + // tris_gia_res("is_loop_vertex",ti) = zs::reinterpret_bits((int)0); + }); // return nm_rings; @@ -881,22 +909,6 @@ int do_global_self_intersection_analysis(Pol& pol, zs::Vector island_buffer{verts.get_allocator(),verts.size()}; - gia_res.resize(verts.size()); - pol(zs::range(verts.size()),[gia_res = proxy({},gia_res)] ZS_LAMBDA(int ni) mutable { - // nodal_colors[ni] = 0; - gia_res("ring_mask",ni) = zs::reinterpret_bits((int)0); - gia_res("color_mask",ni) = zs::reinterpret_bits((int)0); - gia_res("type_mask",ni) = zs::reinterpret_bits((int)0); - gia_res("is_corner",ni) = zs::reinterpret_bits((int)0); - }); - tris_gia_res.resize(tris.size()); - pol(zs::range(tris.size()),[tris_gia_res = proxy({},tris_gia_res)] ZS_LAMBDA(int ti) mutable { - // nodal_colors[ni] = 0; - tris_gia_res("ring_mask",ti) = zs::reinterpret_bits((int)0); - tris_gia_res("color_mask",ti) = zs::reinterpret_bits((int)0); - tris_gia_res("type_mask",ti) = zs::reinterpret_bits((int)0); - // tris_gia_res("is_corner",ti) = zs::reinterpret_bits((int)0); - }); zs::Vector> edge_topos{tris.get_allocator(),tris.size() * 3}; pol(range(tris.size()),[ @@ -917,8 +929,13 @@ int do_global_self_intersection_analysis(Pol& pol, for(int ri = 0;ri != nm_rings;++ri) { auto rsize = (size_t)ringSize.getVal(ri); + // if(output_intermediate_information) - printf("ring[%d] Size : %d\n",ri,rsize); + // printf("ring[%d] Size : %d\n",ri,rsize); + + + int cur_ri_mask = 1 << (ri % 32); + int ri_offset = ri / 32; // edge_topo_type dc_edge_topos{tris.get_allocator(),rsize * 6}; table_int_type disable_points{tris.get_allocator(),rsize * 8}; @@ -931,6 +948,9 @@ int do_global_self_intersection_analysis(Pol& pol, ringTag = proxy(ringTag), // output_intermediate_information, ri, + cur_ri_mask = cur_ri_mask, + ri_offset = ri_offset, + ring_mask_width = ring_mask_width, halfedges = proxy({},halfedges), // topo_tag = zs::SmallString(topo_tag), // dc_edge_topos = proxy(dc_edge_topos), @@ -949,11 +969,11 @@ int do_global_self_intersection_analysis(Pol& pol, auto ti = pair[1]; int cur_ri_mask = 1 << ri; - auto tring_mask = zs::reinterpret_bits(tris_gia_res("ring_mask",ti)); + auto tring_mask = zs::reinterpret_bits(tris_gia_res("ring_mask",ti * ring_mask_width + ri_offset)); tring_mask |= cur_ri_mask; - tris_gia_res("ring_mask",ti) = zs::reinterpret_bits(tring_mask); + tris_gia_res("ring_mask",ti * ring_mask_width + ri_offset) = zs::reinterpret_bits(tring_mask); - // halfedges("broken",hi) = (T)1.0; + halfedges("broken",hi) = (T)1.0; // auto ti = pair[1]; // auto type = zs::reinterpret_bits(ints_buffer("type",isi)); @@ -990,7 +1010,7 @@ int do_global_self_intersection_analysis(Pol& pol, // } auto corner_idx = zs::reinterpret_bits(ints_buffer("corner_idx",isi)); if(corner_idx >= 0){ - gia_res("is_corner",corner_idx) = (T)1.0; + gia_res("is_loop_vertex",corner_idx * ring_mask_width + ri_offset) = (T)1.0; disable_points.insert(corner_idx); } @@ -1029,7 +1049,7 @@ int do_global_self_intersection_analysis(Pol& pol, // }); // std::cout << "size of conn_of_first_ring : " << conn_of_first_ring.size() << std::endl; // } - std::cout << "ring[" << ri << "] : " << nm_islands << "\tnm_broken_edges : " << disable_lines.size() << "\tnm_broken_corners : " << disable_points.size() << std::endl; + // std::cout << "ring[" << ri << "] : " << nm_islands << "\tnm_broken_edges : " << disable_lines.size() << "\tnm_broken_corners : " << disable_points.size() << std::endl; zs::Vector nm_cmps_every_island_count{verts.get_allocator(),(size_t)nm_islands}; @@ -1038,15 +1058,17 @@ int do_global_self_intersection_analysis(Pol& pol, nm_cmps_every_island_count = proxy(nm_cmps_every_island_count)] ZS_LAMBDA(int i) mutable { nm_cmps_every_island_count[i] = 0; }); + + // it is a really bad idea to use mustExclude here, as this tag might no be locally significant pol(zs::range(verts.size()),[ exec_tag, verts = proxy({},verts), island_buffer = proxy(island_buffer), nm_cmps_every_island_count = proxy(nm_cmps_every_island_count)] ZS_LAMBDA(int vi) mutable { auto island_idx = island_buffer[vi]; - if(verts.hasProperty("mustExclude")) - if(verts("mustExclude",vi) > (T)0.5) - return; + // if(verts.hasProperty("mustExclude")) + // if(verts("mustExclude",vi) > (T)0.5) + // return; atomic_add(exec_tag,&nm_cmps_every_island_count[island_idx],(int)1); }); @@ -1070,18 +1092,22 @@ int do_global_self_intersection_analysis(Pol& pol, break; } } - auto cur_ri_mask = (int)1 << ri; + + + // auto cur_ri_mask = (int)1 << ri; for(int i = 0;i != nm_islands;++i) std::cout << nm_cmps_every_island_count.getVal(i) << "\t"; - std::cout << "max_island = " << max_island_idx << std::endl; + // std::cout << "max_island = " << max_island_idx << std::endl; // std::cout << std::endl; pol(zs::range(verts.size()),[ gia_res = proxy({},gia_res), nm_islands, cur_ri_mask, + ri_offset = ri_offset, + ring_mask_width = ring_mask_width, black_island_idx, exec_tag, // ints_types = proxy(ints_types), @@ -1091,9 +1117,9 @@ int do_global_self_intersection_analysis(Pol& pol, if(island_idx == max_island_idx) return; // might exceed the integer range - auto ring_mask = zs::reinterpret_bits(gia_res("ring_mask",vi)); - auto color_mask = zs::reinterpret_bits(gia_res("color_mask",vi)); - auto type_mask = zs::reinterpret_bits(gia_res("type_mask",vi)); + auto ring_mask = zs::reinterpret_bits(gia_res("ring_mask",vi * ring_mask_width + ri_offset)); + auto color_mask = zs::reinterpret_bits(gia_res("color_mask",vi * ring_mask_width + ri_offset)); + auto type_mask = zs::reinterpret_bits(gia_res("type_mask",vi * ring_mask_width + ri_offset)); // ring_mask += ((int) << ri) // if(island_idx != max_island_idx/* || ints_types[island_idx] == 1*/){ @@ -1103,20 +1129,23 @@ int do_global_self_intersection_analysis(Pol& pol, type_mask |= cur_ri_mask; if(nm_islands == 3 && island_idx == black_island_idx) color_mask |= cur_ri_mask; - gia_res("ring_mask",vi) = zs::reinterpret_bits(ring_mask); - gia_res("color_mask",vi) = zs::reinterpret_bits(color_mask); - gia_res("type_mask",vi) = zs::reinterpret_bits(type_mask); + gia_res("ring_mask",vi * ring_mask_width + ri_offset) = zs::reinterpret_bits(ring_mask); + gia_res("color_mask",vi * ring_mask_width + ri_offset) = zs::reinterpret_bits(color_mask); + gia_res("type_mask",vi * ring_mask_width + ri_offset) = zs::reinterpret_bits(type_mask); }); } pol(zs::range(gia_res.size()),[ - gia_res = proxy({},gia_res)] ZS_LAMBDA(int vi) mutable { - auto is_corner = gia_res("is_corner",vi); - if(is_corner > (T)0.5) - gia_res("ring_mask",vi) = zs::reinterpret_bits((int)0); + // ring_mask_width = ring_mask_width, + gia_res = proxy({},gia_res)] ZS_LAMBDA(int mi) mutable { + // for(int i = 0;i != ring_mask_width;++i) { + auto is_corner = gia_res("is_loop_vertex",mi); + if(is_corner > (T)0.5) + gia_res("ring_mask",mi) = zs::reinterpret_bits((int)0); + // } }); - return nm_rings; + return ring_mask_width; } template @@ -2231,10 +2260,16 @@ int retrieve_intersection_tri_halfedge_info_of_two_meshes(Pol& pol, edge_A[0] = tri_A[(local_vert_id_A + 0) % 3]; edge_A[1] = tri_A[(local_vert_id_A + 1) % 3]; - // auto ohei_A = zs::reinterpret_bits(halfedges_A("opposite_he",hei_A)); + auto ohei_A = zs::reinterpret_bits(halfedges_A("opposite_he",hei_A)); + + if(edge_A[0] > edge_A[1] && ohei_A >= 0) + return; - // if(edge_A[0] > edge_A[1] && ohei_A >= 0) - // return; + // if(edge_A[0] > edge_A[1]) { + // auto tmp = edge_A[0]; + // edge_A[0] = edge_A[1]; + // edge_A[1] = tmp; + // } vec3 eV_A[2] = {}; for(int i = 0;i != 2;++i) @@ -2261,12 +2296,12 @@ int retrieve_intersection_tri_halfedge_info_of_two_meshes(Pol& pol, intersect_buffers("r",offset) = (T)r; // make sure the opposite he - tri pairs are also inserted // auto opposite_hei_A = zs::reinterpret_bits(halfedges_A("opposite_he",hei_A)); - // if(opposite_hei_A >= 0) { - // offset = atomic_add(exec_tag,&nmIts[0],(int)1); - // intersect_buffers.tuple(dim_c<2>,"pair",offset) = zs::vec{opposite_hei_A,ti_B}.reinterpret_bits(float_c); - // intersect_buffers.tuple(dim_c<3>,"int_points",offset) = intp; - // intersect_buffers("r",offset) = (T)(1 - r); - // } + if(ohei_A >= 0) { + offset = atomic_add(exec_tag,&nmIts[0],(int)1); + intersect_buffers.tuple(dim_c<2>,"pair",offset) = zs::vec{ohei_A,ti_B}.reinterpret_bits(float_c); + intersect_buffers.tuple(dim_c<3>,"int_points",offset) = intp; + intersect_buffers("r",offset) = (T)(1 - r); + } } } }; @@ -3175,7 +3210,7 @@ int do_global_intersection_analysis_with_connected_manifolds(Pol& pol, if(a2b_isi + _A_offset > na2b_isi + _A_offset) incidentItsTab.insert(vec2i{a2b_isi + _A_offset,na2b_isi + _A_offset}); }else { - printf("impossible reaching here, the hi and ohi should both have been inserted\n"); + printf("do_global_intersection_analysis_with_connected_manifolds_new::impossible reaching here, the hi and ohi should both have been inserted\n"); atomic_add(exec_tag,&nmInvalid[0],(int)1); } } @@ -3203,7 +3238,7 @@ int do_global_intersection_analysis_with_connected_manifolds(Pol& pol, } hb = zs::reinterpret_bits(_B_halfedges("next_he",hb)); } - printf("impossible reaching here, the intersection ring seems to be broken\n"); + printf("do_global_intersection_analysis_with_connected_manifolds_new::impossible reaching here, the intersection ring seems to be broken\n"); atomic_add(exec_tag,&nmInvalid[0],(int)1); }); }; From 82a2d638f78675bb1424c0b5cbb18acbc1d053c2 Mon Sep 17 00:00:00 2001 From: ShuliangLu Date: Wed, 26 Jul 2023 17:05:10 +0800 Subject: [PATCH 4/4] fix gia crash encountering broken rings --- .../geometry/kernel/intersection.hpp | 38 ++- .../CuLagrange/geometry/kernel/topology.hpp | 263 ++++++++++++++++-- 2 files changed, 272 insertions(+), 29 deletions(-) diff --git a/projects/CuLagrange/geometry/kernel/intersection.hpp b/projects/CuLagrange/geometry/kernel/intersection.hpp index 0cb88d5ebb..ba17ace0aa 100644 --- a/projects/CuLagrange/geometry/kernel/intersection.hpp +++ b/projects/CuLagrange/geometry/kernel/intersection.hpp @@ -740,6 +740,7 @@ int do_global_self_intersection_analysis(Pol& pol, {"pair",2}, {"int_points",3}, {"r",1}, + {"is_broken",1} },max_nm_intersections}; if(!halfedges.hasProperty("broken")) @@ -747,6 +748,7 @@ int do_global_self_intersection_analysis(Pol& pol, TILEVEC_OPS::fill(pol,halfedges,"broken",(T)0.0); auto nm_insts = retrieve_self_intersection_tri_halfedge_list_info(pol,verts,xtag,tris,halfedges,ints_buffer); + TILEVEC_OPS::fill(pol,ints_buffer,"is_broken",(T)0); table_vec2i_type cftab{ints_buffer.get_allocator(),(size_t)nm_insts}; cftab.reset(pol,true); zs::Vector cfbuffer{ints_buffer.get_allocator(),(size_t)nm_insts}; @@ -821,7 +823,9 @@ int do_global_self_intersection_analysis(Pol& pol, printf("do_global_self_intersection_analysis::impossible reaching here, the hi and ohi should both have been inserted %f %f %f\n",(float)hr,(float)ohr,(float)ints_buffer("r",isi)); + ints_buffer("is_broken",isi) = (T)1.0; atomic_add(exec_tag,&nmInvalid[0],(int)1); + return; } } auto corner_idx = zs::reinterpret_bits(ints_buffer("corner_idx",isi)); @@ -854,12 +858,15 @@ int do_global_self_intersection_analysis(Pol& pol, } printf("do_global_self_intersection_analysis::impossible reaching here with broken insertion ring %f\n",(float)ints_buffer("r",isi)); + ints_buffer("is_broken",isi) = (T)1.0; atomic_add(exec_tag,&nmInvalid[0],(int)1); }); auto nmInvalidCount = nmInvalid.getVal(0); if(nmInvalidCount > 0) - throw std::runtime_error("SELF GIA invalid state detected"); + printf("SELF GIA invalid state detected\n"); + // there might be some broken rings + auto nmEntries = incidentItsTab.size(); zs::Vector> conn_topo{tris.get_allocator(),nmEntries}; @@ -871,6 +878,20 @@ int do_global_self_intersection_analysis(Pol& pol, auto nm_rings = mark_disconnected_island(pol,conn_topo,ringTag); + zs::Vector is_broken_rings{ringTag.get_allocator(),(size_t)nm_rings}; + pol(zs::range(is_broken_rings),[] ZS_LAMBDA(auto& is_br) mutable {is_br = 0;}); + pol(zs::range(nm_insts),[ringTag = proxy(ringTag),is_broken_rings = proxy(is_broken_rings),ints_buffer = proxy({},ints_buffer)] ZS_LAMBDA(int isi) mutable { + if(ints_buffer("is_broken",isi) > (T)0.5) { + auto ring_id = ringTag[isi]; + is_broken_rings[ring_id] = 1; + } + }); + + std::cout << "broken_ring_tag : "; + for(int i = 0;i != nm_rings;++i) + std::cout << is_broken_rings.getVal(i) << "\t"; + std::cout << std::endl; + // std::cout << "finish Mark disconnected island with nm_rings : " << nm_rings << std::endl; auto ring_mask_width = (nm_rings + 31) / 32; @@ -907,6 +928,11 @@ int do_global_self_intersection_analysis(Pol& pol, atomic_add(exec_tag,&ringSize[ringTag[isi]],(int)1); }); + // pol(zs::range(nm_rings),[ringSize = proxy(ringSize),is_broken_rings = proxy(is_broken_rings)] ZS_LAMBDA(int ri) mutable { + // if(is_broken_rings[ri]) + // ringSize[ri] = 0; + // }); + zs::Vector island_buffer{verts.get_allocator(),verts.size()}; @@ -929,7 +955,11 @@ int do_global_self_intersection_analysis(Pol& pol, for(int ri = 0;ri != nm_rings;++ri) { auto rsize = (size_t)ringSize.getVal(ri); - + // if(rsize == 0) + // continue; + auto is_broken_ring = is_broken_rings.getVal(ri); + if(is_broken_ring) + continue; // if(output_intermediate_information) // printf("ring[%d] Size : %d\n",ri,rsize); @@ -2079,8 +2109,8 @@ int do_global_self_intersection_analysis_on_surface_mesh_info(Pol& pol, tri_pairs[1] = tris.pack(dim_c<3>,topo_tag,tb,int_c); for(int t = 0;t != 2;++t) { - zs::vec out_edges[3] = {}; - elm_to_edges(tri_pairs[t],out_edges); + // zs::vec out_edges[3] = {}; + auto out_edges = elm_to_edges(tri_pairs[t]); for(int i = 0;i != 3;++i) { auto a = out_edges[i][0]; auto b = out_edges[i][1]; diff --git a/projects/CuLagrange/geometry/kernel/topology.hpp b/projects/CuLagrange/geometry/kernel/topology.hpp index c26c48a28f..6f993ca80c 100644 --- a/projects/CuLagrange/geometry/kernel/topology.hpp +++ b/projects/CuLagrange/geometry/kernel/topology.hpp @@ -383,25 +383,82 @@ namespace zeno { return true; } - constexpr void elm_to_edges(const zs::vec& single_edge,zs::vec out_edges[1]) { - out_edges[0] = single_edge; - } + // constexpr auto elm_to_edges(const zs::vec& single_edge) { + // zs::vec out_edges[1] = {}; + // out_edges[0] = single_edge; + // return out_edges; + // } - constexpr void elm_to_edges(const zs::vec& tri,zs::vec out_edges[3]) { - out_edges[0] = zs::vec{tri[0],tri[1]}; - out_edges[1] = zs::vec{tri[1],tri[2]}; - out_edges[2] = zs::vec{tri[2],tri[0]}; - } + // constexpr auto elm_to_edges(const zs::vec& tri) { + // zs::vec out_edges[3] {}; + // out_edges[0] = zs::vec{tri[0],tri[1]}; + // out_edges[1] = zs::vec{tri[1],tri[2]}; + // out_edges[2] = zs::vec{tri[2],tri[0]}; + // return out_edges; + // } + + template 1)> = 0> + constexpr auto elm_to_edges(const zs::VecInterface& elm) { + using Ti = typename VecTi::value_type; + constexpr auto CODIM = VecTi::extent; + constexpr auto NM_EDGES = (CODIM - 1) * (CODIM) / 2; + + zs::vec out_edges[NM_EDGES] = {}; + int nm_out_edges = 0; + for(int i = 0;i != CODIM;++i) + for(int j = i + 1;j != CODIM;++j) + out_edges[nm_out_edges++] = zs::vec{elm[i],elm[j]}; - constexpr void elm_to_edges(const zs::vec& tet,zs::vec out_edges[6]) { - out_edges[0] = zs::vec{tet[0],tet[1]}; - out_edges[1] = zs::vec{tet[0],tet[2]}; - out_edges[2] = zs::vec{tet[0],tet[3]}; - out_edges[3] = zs::vec{tet[1],tet[2]}; - out_edges[4] = zs::vec{tet[1],tet[3]}; - out_edges[5] = zs::vec{tet[2],tet[3]}; + return out_edges; } + // template + // constexpr auto elm_to_edges(const VecI& tet,zs::wrapv) { + // constexpr int NM_EDGES = codim * (codim - 1) / 2; + // zs::vec out_edges[NM_EDGES] = {}; + // if(codim == 4) { + + // out_edges[0] = zs::vec{tet[0],tet[1]}; + // out_edges[1] = zs::vec{tet[0],tet[2]}; + // out_edges[2] = zs::vec{tet[0],tet[3]}; + // out_edges[3] = zs::vec{tet[1],tet[2]}; + // out_edges[4] = zs::vec{tet[1],tet[3]}; + // out_edges[5] = zs::vec{tet[2],tet[3]}; + // return out_edges; + // } + // if(codim == 3) { + // zs::vec out_edges[3] {}; + // out_edges[0] = zs::vec{tri[0],tri[1]}; + // out_edges[1] = zs::vec{tri[1],tri[2]}; + // out_edges[2] = zs::vec{tri[2],tri[0]}; + // return out_edges; + // } + // } + + // constexpr auto elm_to_edges(const zs::vec& single_edge) { + // zs::vec out_edges[1] = {}; + // out_edges[0] = single_edge; + // return out_edges; + // } + + // constexpr auto elm_to_edges(const zs::vec& tri) { + // zs::vec out_edges[3] {}; + // out_edges[0] = zs::vec{tri[0],tri[1]}; + // out_edges[1] = zs::vec{tri[1],tri[2]}; + // out_edges[2] = zs::vec{tri[2],tri[0]}; + // return out_edges; + // } + + // constexpr auto elm_to_edges(const zs::vec& tet) { + // zs::vec out_edges[6] = {}; + // out_edges[0] = zs::vec{tet[0],tet[1]}; + // out_edges[1] = zs::vec{tet[0],tet[2]}; + // out_edges[2] = zs::vec{tet[0],tet[3]}; + // out_edges[3] = zs::vec{tet[1],tet[2]}; + // out_edges[4] = zs::vec{tet[1],tet[3]}; + // out_edges[5] = zs::vec{tet[2],tet[3]}; + // return out_edges; + // } // template // bool topo_to_incident_matrix(Pol& pol, // const TopoTileVec& topo, @@ -1055,20 +1112,21 @@ namespace zeno { } - template= 2), (VecTI::etent <= 4)> = 0*/> + template= 2), (VecTI::etent <= 4)> = 0*/> void topological_incidence_matrix(Pol& pol, // size_t nm_points, - const zs::Vector& topos, + const TopoRangT& topos, zs::SparseMatrix& spmat) { using namespace zs; using ICoord = zs::vec; - constexpr auto CDIM = VecTI::extent; + // constexpr auto CDIM = VecTI::extent; + constexpr auto CDIM = RM_CVREF_T(topos[0])::extent; constexpr auto space = Pol::exec_tag::value; constexpr auto execTag = wrapv{}; zs::Vector max_pi_vec{topos.get_allocator(),1}; max_pi_vec.setVal(0); - pol(zs::range(topos),[max_pi_vec = proxy(max_pi_vec),execTag,CDIM] ZS_LAMBDA(const auto& topo) { + pol(zs::range(topos),[max_pi_vec = proxy(max_pi_vec),execTag,CDIM] ZS_LAMBDA(const auto& topo) mutable { for(int i = 0;i != CDIM;++i) if(topo[i] >= 0) atomic_max(execTag,&max_pi_vec[0],(int)topo[i]); @@ -1135,8 +1193,6 @@ namespace zeno { // printf("p[%d] -> t[%d]\n",pi,ti); } }); - - // std::cout << "finish computing p2ts" << std::endl; } @@ -1178,6 +1234,8 @@ namespace zeno { // pol(zs::range(is.size()),[is = proxy(is),js = proxy(js)] ZS_LAMBDA(int i) mutable {printf("ijs[%d] : %d %d\n",i,is[i],js[i]);}); // std::cout << "topos.size() = " << topos.size() << std::endl; + // for(int i = 0;i != topos.size();++i) + // std::cout << topos.getVal(i)[0] << "\t" << topos.getVal(i)[1] << std::endl; // spmat = zs::SparseMatrix{topos.get_allocator(),(int)topos.size(),(int)topos.size()}; spmat.build(pol,(int)topos.size(),(int)topos.size(),zs::range(is),zs::range(js)/*,zs::range(rs)*/,zs::true_c); @@ -1189,16 +1247,18 @@ namespace zeno { // spmat._vals.reset((int)1); } - template + template void topological_coloring(Pol& pol, // int nm_points, - const zs::Vector& topo, - zs::Vector& colors) { + const TopoRangeT& topo, + ColorRangeT& colors) { using namespace zs; constexpr auto space = Pol::exec_tag::value; + using Ti = RM_CVREF_T(colors[0]); + colors.resize(topo.size()); - zs::SparseMatrix topo_incidence_matrix{topo.get_allocator(),topo.size(),topo.size()}; + zs::SparseMatrix topo_incidence_matrix{topo.get_allocator(),(int)topo.size(),(int)topo.size()}; // std::cout << "compute incidence matrix " << std::endl; topological_incidence_matrix(pol,topo,topo_incidence_matrix); // std::cout << "finish compute incidence matrix " << std::endl; @@ -1220,10 +1280,163 @@ namespace zeno { w = v; }); } + + // pol(zs::range()) weights = weights.clone(colors.memoryLocation()); + // for(int i = 0;i != weights.size();++i) + // printf("w[%d] : %u\n",i,weights.getVal(i)); auto iterRef = maximum_independent_sets(pol, topo_incidence_matrix, weights, colors); std::cout << "nm_colors : " << iterRef << std::endl; + pol(zs::range(colors),[] ZS_LAMBDA(auto& clr) mutable {clr = clr - (Ti)1;}); + + } + + template + void sort_topology_by_coloring_tag(Pol& pol, + const COLOR_RANGE& colors, + REORDERED_MAP_RANGE& reordered_map, + EXCLUSIVE_OFFSET_RANGE& offset_out) { + using namespace zs; + constexpr auto space = Pol::exec_tag::value; + constexpr auto exec_tag = wrapv{}; + + // zs::Vector reordered_map{colors.get_allocator(),colors.size()}; + reordered_map.resize(colors.size()); + zs::Vector max_color{colors.get_allocator(),1}; + max_color.setVal(0); + + pol(zs::range(colors.size()),[ + colors = proxy(colors),\ + exec_tag = exec_tag, + max_color = proxy(max_color)] ZS_LAMBDA(int ci) mutable { + auto color = (int)colors[ci]; + atomic_max(exec_tag,&max_color[0],color); + }); + + int nm_total_colors = max_color.getVal(0) + 1; + // zs::bht color_buffer{} + zs::Vector nm_colors{colors.get_allocator(),nm_total_colors}; + pol(zs::range(nm_colors),[] ZS_LAMBDA(auto& nclr) mutable {nclr = 0;}); + pol(zs::range(colors),[nm_colors = proxy(nm_colors),exec_tag] ZS_LAMBDA(const auto& clrf) mutable { + auto clr = (int)clrf; + atomic_add(exec_tag,&nm_colors[clr],1); + }); + + zs::Vector exclusive_offsets{colors.get_allocator(),nm_total_colors}; + pol(zs::range(exclusive_offsets),[] ZS_LAMBDA(auto& eoffset) {eoffset = 0;}); + exclusive_scan(pol,std::begin(nm_colors),std::end(nm_colors),std::begin(exclusive_offsets)); + pol(zs::range(nm_colors),[] ZS_LAMBDA(auto& nclr) {nclr = 0;}); + + offset_out.resize(nm_total_colors); + + pol(zip(zs::range(exclusive_offsets.size()),exclusive_offsets),[offset_out = proxy(offset_out)] ZS_LAMBDA(auto i,auto offset) mutable {offset_out[i] = offset;}); + pol(zs::range(colors.size()),[ + nm_colors = proxy(nm_colors), + colors = proxy(colors), + exec_tag, + exclusive_offsets = proxy(exclusive_offsets), + reordered_map = proxy(reordered_map)] ZS_LAMBDA(auto ci) mutable { + auto clr = (int)colors[ci]; + auto offset = atomic_add(exec_tag,&nm_colors[clr],1); + auto eoffset = exclusive_offsets[clr]; + + reordered_map[eoffset + offset] = ci; + }); + + // zs::Vector topos_copy{topos.get_allocator(),topos.size()}; + // pol(zip(zs::range(topos.size()),topos),[topos_copy = proxy(topos_copy)] ZS_LAMBDA(auto ti,const auto& topo) mutable {topos_copy[ti] = topo;}); + + // pol(zip(zs::range(topos.size()),topos),[ + // topos_copy = proxy(topos_copy), + // reordered_map = proxy(reordered_map)] ZS_LAMBDA(auto ti,auto& topo) mutable {topo = topos_copy[reordered_map[ti]];}); + } + + template> + zs::Vector tilevec_topo_to_zsvec_topo(Pol& pol,const TopoTileVec& source,zs::wrapv) { + zs::Vector out_topo{source.get_allocator(),source.size()}; + auto sr = zs::range(source, "inds", zs::dim_c, zs::int_c); + pol(zip(sr, out_topo), []ZS_LAMBDA(auto id, VecTi& dst) mutable { + if constexpr (std::is_integral_v) + dst[0] = id; + else + dst = id; + }); + return out_topo; + } + + template + void retrieve_tri_bending_topology(Pol& pol, + const TriTileVec& tris, + const HalfEdgeTileVec& halfedges, + zs::Vector>& tb_topos) { + using namespace zs; + constexpr auto space = RM_CVREF_T(pol)::exec_tag::value; + constexpr auto exec_tag = wrapv{}; + + // zs::Vector nm_interior_edges{halfedges.get_allocator(),1}; + // nm_interior_edges.setVal(0); + + zs::bht interior_edges{halfedges.get_allocator(),halfedges.size()}; + interior_edges.reset(pol,true); + + pol(zs::range(halfedges.size()),[ + halfedges = proxy({},halfedges), + exec_tag, + interior_edges = proxy(interior_edges)] ZS_LAMBDA(int hi) mutable { + auto ohi = zs::reinterpret_bits(halfedges("opposite_he",hi)); + // the boundary halfedge will return -1 for opposite_he here, so it is automatically neglected + if(ohi < hi) + return; + interior_edges.insert(hi); + }); + + tb_topos.resize(interior_edges.size()); + pol(zs::zip(zs::range(interior_edges.size()),interior_edges._activeKeys),[ + tb_topos = proxy(tb_topos), + halfedges = proxy({},halfedges), + tris = proxy({},tris)] ZS_LAMBDA(auto id,auto hi_vec) mutable { + auto hi = hi_vec[0]; + auto ti = zs::reinterpret_bits(halfedges("to_face",hi)); + auto vid = zs::reinterpret_bits(halfedges("local_vertex_id",hi)); + auto ohi = zs::reinterpret_bits(halfedges("opposite_he",hi)); + auto oti = zs::reinterpret_bits(halfedges("to_face",ohi)); + auto ovid = zs::reinterpret_bits(halfedges("local_vertex_id",ohi)); + + auto tri = tris.pack(dim_c<3>,"inds",ti,int_c); + auto otri = tris.pack(dim_c<3>,"inds",oti,int_c); + + tb_topos[id] = zs::vec(tri[(vid + 0) % 3],tri[(vid + 1) % 3],tri[(vid + 2) % 3],otri[(ovid + 2) % 3]); + }); + } + + template + void retrieve_edges_topology(Pol& pol, + const zs::Vector& src_topos, + zs::Vector>& edges_topos) { + using namespace zs; + constexpr auto space = RM_CVREF_T(pol)::exec_tag::value; + // constexpr auto CDIM = VecTi::extent; + // constexpr auto NM_EDGES = CDIM * (CDIM - 1) / 2; + + zs::bht edges_tab{src_topos.get_allocator(),src_topos.size() * 6}; + edges_tab.reset(pol,true); + pol(zs::range(src_topos.size()),[ + src_topos = proxy(src_topos), + edges_tab = proxy(edges_tab)] ZS_LAMBDA(int ei) mutable { + auto elm_edges = elm_to_edges(src_topos[ei]); + for(int i = 0;i != NM_EDGES;++i) { + auto edge = elm_edges[i]; + if(edge[0] < edge[1]) + edges_tab.insert(zs::vec{edge[0],edge[1]}); + else + edges_tab.insert(zs::vec{edge[1],edge[0]}); + } + }); + + edges_topos.resize(edges_tab.size()); + pol(zip(zs::range(edges_tab.size()),zs::range(edges_tab._activeKeys)),[ + edges_topos = proxy(edges_topos)] ZS_LAMBDA(auto ei,const auto& edge){edges_topos[ei] = edge;}); } template