fixed a few of bugs with solver;

now it can work on a variety of cases.
Known bug: due to current design of primitive, some primitives with sharp rounding may not be detected, so that their integrand may not be correct, as well as their output mesh.

frontend/src/implicit_surface_network_processor.cpp

@@ -2,6 +2,7 @@
 
 #include <extract_patch.hpp>
 #include <implicit_arrangement.hpp>
+#include <optional>
 #include <patch_connectivity.hpp>
 #include <pair_faces.hpp>
 #include <topology_ray_shooting.hpp>
@@ -129,7 +130,7 @@ solve_result_t ImplicitSurfaceNetworkProcessor::run(const virtual_node_t& tree_n
 
     // compute arrangement in each tet
     // HINT: we skip robust test for this part for now
-    stl_vector_mp<std::shared_ptr<arrangement_t>> cut_results{};
+    stl_vector_mp<std::optional<arrangement_t>> cut_results{};
     uint32_t                                      num_1_func    = 0;
     uint32_t                                      num_2_func    = 0;
     uint32_t                                      num_more_func = 0;
@@ -143,7 +144,7 @@ solve_result_t ImplicitSurfaceNetworkProcessor::run(const virtual_node_t& tree_n
             const auto start_index              = start_index_of_tet[i];
             const auto active_funcs_in_curr_tet = start_index_of_tet[i + 1] - start_index;
             if (active_funcs_in_curr_tet == 0) {
-                cut_results.emplace_back(nullptr);
+                cut_results.emplace_back(std::nullopt);
                 continue;
             }
             g_timers_manager.push_timer(active_funcs_in_curr_tet == 1
@@ -160,7 +161,7 @@ solve_result_t ImplicitSurfaceNetworkProcessor::run(const virtual_node_t& tree_n
                                             -vertex_scalar_values[tet[2]][fid],
                                             -vertex_scalar_values[tet[3]][fid]});
             }
-            cut_results.emplace_back(std::make_shared<arrangement_t>(std::move(compute_arrangement(planes))));
+            cut_results.emplace_back(std::move(compute_arrangement(planes)));
             switch (active_funcs_in_curr_tet) {
                 case 1:
                     g_timers_manager.pop_timer("implicit arrangements calculation (1 func)");

frontend/src/patch_propagator.cpp

@@ -215,6 +215,7 @@ dynamic_bitset_mp<> PatchPropagator::filter_cells_by_boolean(const virtual_node_
                                                              const stl_vector_mp<uint32_t>&            leaf_index_of_primitive,
                                                              const stl_vector_mp<dynamic_bitset_mp<>>& function_cell_labels)
 {
+    // system("pause");
     const auto                       num_funcs             = get_primitive_count();
     const auto                       num_cells             = function_cell_labels[0].size();
     const auto                       max_parent_node_index = tree_root.inner_index;

implicit_arrangements/src/ia_cut_1_face.cpp

@@ -33,8 +33,9 @@ std::array<uint32_t, 3> ia_cut_1_face(ia_complex_t&                ia_complex,
         return static_cast<uint32_t>(vertices.size() - 1);
     };
 
+    if (o0 == 0 && o1 == 0) {
         // do nothing if edge is coplanar with plane (o0 == 0 && o1 == 0)
-    if (o0 >= 0 && o1 >= 0) {
+    } else if (o0 >= 0 && o1 >= 0) {
         positive_subedge_id = eid;
     } else if (o0 <= 0 && o1 <= 0) {
         negative_subedge_id = eid;

network_process/interface/extract_patch.hpp

@@ -6,7 +6,7 @@
 ISNP_API void extract_iso_mesh(uint32_t                                           num_1_func,
                                uint32_t                                           num_2_func,
                                uint32_t                                           num_more_func,
-                               const stl_vector_mp<std::shared_ptr<arrangement_t>>& cut_results,
+                               const stl_vector_mp<std::optional<arrangement_t>>& cut_results,
                                const stl_vector_mp<uint32_t>&                     func_in_tet,
                                const stl_vector_mp<uint32_t>&                     start_index_of_tet,
                                const tetrahedron_mesh_t&                          background_mesh,

network_process/interface/pair_faces.hpp

@@ -12,7 +12,7 @@ ISNP_API void compute_patch_order(const iso_edge_t
                                   const stl_vector_mp<tetrahedron_vertex_indices_t>         &tets,
                                   const stl_vector_mp<iso_vertex_t>                         &iso_verts,
                                   const stl_vector_mp<polygon_face_t>                       &iso_faces,
-                                  const stl_vector_mp<std::shared_ptr<arrangement_t>>       &cut_results,
+                                  const stl_vector_mp<std::optional<arrangement_t>>         &cut_results,
                                   const stl_vector_mp<uint32_t>                             &func_in_tet,
                                   const stl_vector_mp<uint32_t>                             &start_index_of_tet,
                                   const flat_hash_map_mp<uint32_t, stl_vector_mp<uint32_t>> &incident_tets,
@@ -34,7 +34,7 @@ ISNP_API void pair_patches_in_tets(const iso_edge_t
                                    const stl_vector_mp<uint32_t>                     &containing_tetIds,
                                    const stl_vector_mp<tetrahedron_vertex_indices_t> &tets,
                                    const stl_vector_mp<polygon_face_t>               &iso_faces,
-                                   const stl_vector_mp<std::shared_ptr<arrangement_t>> &cut_results,
+                                   const stl_vector_mp<std::optional<arrangement_t>> &cut_results,
                                    const stl_vector_mp<uint32_t>                     &func_in_tet,
                                    const stl_vector_mp<uint32_t>                     &start_index_of_tet,
                                    const stl_vector_mp<uint32_t>                     &patch_of_face_mapping,

network_process/interface/topology_ray_shooting.hpp

@@ -75,7 +75,7 @@ struct equal_to<face_with_orient_t> {
 
 // topological ray shooting for implicit arrangement
 ISNP_API void topo_ray_shooting(const tetrahedron_mesh_t                          &tet_mesh,
-                                const stl_vector_mp<std::shared_ptr<arrangement_t>> &cut_results,
+                                const stl_vector_mp<std::optional<arrangement_t>> &cut_results,
                                 const stl_vector_mp<iso_vertex_t>                 &iso_verts,
                                 const stl_vector_mp<polygon_face_t>               &iso_faces,
                                 const stl_vector_mp<stl_vector_mp<uint32_t>>      &patches,

network_process/src/extract_patch.cpp

@@ -7,7 +7,7 @@
 ISNP_API void extract_iso_mesh(uint32_t                                           num_1_func,
                                uint32_t                                           num_2_func,
                                uint32_t                                           num_more_func,
-                               const stl_vector_mp<std::shared_ptr<arrangement_t>>& cut_results,
+                               const stl_vector_mp<std::optional<arrangement_t>>& cut_results,
                                const stl_vector_mp<uint32_t>&                     func_in_tet,
                                const stl_vector_mp<uint32_t>&                     start_index_of_tet,
                                const tetrahedron_mesh_t&                          background_mesh,
@@ -56,8 +56,8 @@ ISNP_API void extract_iso_mesh(uint32_t
 
     //
     for (uint32_t i = 0; i < n_tets; i++) {
-        if (cut_results[i]) {
+        if (cut_results[i].has_value()) {
-            const auto& arrangement = *cut_results[i].get();
+            const auto& arrangement = cut_results[i].value();
             const auto& vertices    = arrangement.vertices;
             const auto& faces       = arrangement.faces;
             auto        start_index = start_index_of_tet[i];

network_process/src/pair_faces.cpp

@@ -9,7 +9,7 @@ ISNP_API void compute_patch_order(const iso_edge_t
                                   const stl_vector_mp<tetrahedron_vertex_indices_t>         &tets,
                                   const stl_vector_mp<iso_vertex_t>                         &iso_verts,
                                   const stl_vector_mp<polygon_face_t>                       &iso_faces,
-                                  const stl_vector_mp<std::shared_ptr<arrangement_t>>       &cut_results,
+                                  const stl_vector_mp<std::optional<arrangement_t>>         &cut_results,
                                   const stl_vector_mp<uint32_t>                             &func_in_tet,
                                   const stl_vector_mp<uint32_t>                             &start_index_of_tet,
                                   const flat_hash_map_mp<uint32_t, stl_vector_mp<uint32_t>> &incident_tets,
@@ -31,7 +31,7 @@ ISNP_API void compute_patch_order(const iso_edge_t
     if (containing_tets.size() == 1) {
         //        std::cout << ">>>>>>>> iso-edge in tet" << std::endl;
         auto tet_id = *containing_tets.begin();
-        pair_patches_in_one_tet(*cut_results[tet_id].get(), iso_faces, iso_edge, patch_of_face_mapping, half_patch_adj_list);
+        pair_patches_in_one_tet(cut_results[tet_id].value(), iso_faces, iso_edge, patch_of_face_mapping, half_patch_adj_list);
     } else {
         const auto  v1          = iso_edge.v1;
         const auto  v2          = iso_edge.v2;
@@ -206,7 +206,7 @@ ISNP_API void pair_patches_in_tets(const iso_edge_t
                                    const stl_vector_mp<uint32_t>                     &containing_tetIds,
                                    const stl_vector_mp<tetrahedron_vertex_indices_t> &tets,
                                    const stl_vector_mp<polygon_face_t>               &iso_faces,
-                                   const stl_vector_mp<std::shared_ptr<arrangement_t>> &cut_results,
+                                   const stl_vector_mp<std::optional<arrangement_t>> &cut_results,
                                    const stl_vector_mp<uint32_t>                     &func_in_tet,
                                    const stl_vector_mp<uint32_t>                     &start_index_of_tet,
                                    const stl_vector_mp<uint32_t>                     &patch_of_face_mapping,
@@ -321,7 +321,7 @@ ISNP_API void pair_patches_in_tets(const iso_edge_t
             }
         } else {
             // non-empty tet i
-            const auto &arrangement = *cut_results[i].get();
+            const auto &arrangement = cut_results[i].value();
             const auto &vertices    = arrangement.vertices;
             const auto &faces       = arrangement.faces;
             auto        start_index = start_index_of_tet[i];
@@ -442,7 +442,7 @@ ISNP_API void pair_patches_in_tets(const iso_edge_t
     // the orientation of an iso-face is defined by the smallest-index implicit function passing the iso-face
     auto get_half_iso_face = [&](face_header_t tet_face, int8_t orient, uint32_t &iso_face_id, int8_t &iso_orient) {
         iso_face_id              = iso_face_Id_of_face[tet_face];
-        const auto &cell_complex = *cut_results[tet_face.volume_index].get();
+        const auto &cell_complex = cut_results[tet_face.volume_index].value();
         const auto &faces        = cell_complex.faces;
         auto        supp_pId     = faces[tet_face.local_face_index].supporting_plane;
         if (supp_pId > 3) { // plane 0,1,2,3 are tet boundary planes
@@ -486,7 +486,7 @@ ISNP_API void pair_patches_in_tets(const iso_edge_t
             }
         } else {
             // non-empty tet
-            const auto &cell_complex = *cut_results[tet_id].get();
+            const auto &cell_complex = cut_results[tet_id].value();
             uint32_t    cell_id =
                 (orient == 1 ? cell_complex.faces[tet_face_id].positive_cell : cell_complex.faces[tet_face_id].negative_cell);
             if (cell_id != invalid_index) {

network_process/src/topology_ray_shooting.cpp

@@ -2,7 +2,7 @@
 #include <patch_connectivity.hpp>
 
 ISNP_API void topo_ray_shooting(const tetrahedron_mesh_t                          &tet_mesh,
-                                const stl_vector_mp<std::shared_ptr<arrangement_t>> &cut_results,
+                                const stl_vector_mp<std::optional<arrangement_t>> &cut_results,
                                 const stl_vector_mp<iso_vertex_t>                 &iso_verts,
                                 const stl_vector_mp<polygon_face_t>               &iso_faces,
                                 const stl_vector_mp<stl_vector_mp<uint32_t>>      &patches,
@@ -50,7 +50,7 @@ ISNP_API void topo_ray_shooting(const tetrahedron_mesh_t
             const auto  extreme_v2     = extremal_edge_of_component[2 * i + 1];
             const auto  iso_vId        = iso_vert_on_v_v_next[extreme_v1];
             const auto  tetId          = iso_verts[iso_vId].header.volume_index;
-            const auto &tet_cut_result = *cut_results[tetId].get();
+            const auto &tet_cut_result = cut_results[tetId].value();
             // get local index of v1 and v2 in the tet
             uint32_t    local_v1, local_v2;
             for (uint32_t j = 0; j < 4; ++j) {
@@ -107,7 +107,7 @@ ISNP_API void topo_ray_shooting(const tetrahedron_mesh_t
                     // reached iso-vert at end of the ray
                     const auto  iso_vId_end        = iso_vert_on_v_v_next[v_curr];
                     const auto  end_tetId          = iso_verts[iso_vId_end].header.volume_index;
-                    const auto &end_tet_cut_result = *cut_results[end_tetId].get();
+                    const auto &end_tet_cut_result = cut_results[end_tetId].value();
                     auto        v_next             = next_vert[v_curr];
                     // find local vertex indices in the end tetrahedron
                     for (uint32_t j = 0; j < 4; ++j) {
@@ -307,8 +307,8 @@ ISNP_API void compute_edge_intersection_order(const arrangement_t     &tet_cut_r
         const auto  num_vert = face.vertices.size();
         for (uint32_t i = 0; i < num_vert; ++i) {
             const auto i_next  = (i + 1) % num_vert;
-            const auto vi      = tet_cut_result.vertices[fId][i];
+            const auto vi      = face.vertices[i];
-            const auto vi_next = tet_cut_result.vertices[fId][i_next];
+            const auto vi_next = face.vertices[i_next];
             // add fId to edge (vi, vi_next)
             auto iter_inserted = faces_of_edge.try_emplace(std::make_pair(vi, vi_next), std::make_pair(fId, invalid_index));
             if (!iter_inserted.second) { // inserted before
@@ -359,8 +359,8 @@ ISNP_API void compute_edge_intersection_order(const arrangement_t     &tet_cut_r
         // visit all edges of face, find edge_prev, edge_next and edge_on_vu
         for (uint32_t i = 0; i < num_vert; ++i) {
             const auto i_next  = (i + 1) % num_vert;
-            const auto vi      = tet_cut_result.vertices[f_curr][i];
+            const auto vi      = face.vertices[i];
-            const auto vi_next = tet_cut_result.vertices[f_curr][i_next];
+            const auto vi_next = face.vertices[i_next];
             if (is_on_edge_vu[vi] && !is_on_edge_vu[vi_next]) {
                 auto &two_faces  = faces_of_edge[std::make_pair(vi, vi_next)];
                 auto  other_face = (two_faces.first == f_curr) ? two_faces.second : two_faces.first;