save all primitives, including deleted (#463)

* save all primitives, including deleted
wildmeshing · Oct 27, 2023 · af3a3f6 · af3a3f6
1 parent 3542219
commit af3a3f6
Show file tree

Hide file tree

Showing 4 changed files with 116 additions and 22 deletions.
diff --git a/src/wmtk/Mesh.cpp b/src/wmtk/Mesh.cpp
@@ -38,16 +38,22 @@ PrimitiveType Mesh::top_simplex_type() const
 }
 
 std::vector<Tuple> Mesh::get_all(PrimitiveType type) const
+{
+    return get_all(type, false);
+}
+
+std::vector<Tuple> Mesh::get_all(PrimitiveType type, const bool include_deleted) const
 {
     ConstAccessor<char> flag_accessor = get_flag_accessor(type);
     const attribute::CachingAccessor<char>& flag_accessor_indices = flag_accessor.index_access();
     std::vector<Tuple> ret;
     long cap = capacity(type);
     ret.reserve(cap);
     for (size_t index = 0; index < cap; ++index) {
-        if ((flag_accessor_indices.const_scalar_attribute(index) & 1)) {
+        if (flag_accessor_indices.const_scalar_attribute(index) & 1)
             ret.emplace_back(tuple_from_id(type, index));
-        }
+        else if (include_deleted)
+            ret.emplace_back();
     }
     return ret;
 }

diff --git a/src/wmtk/Mesh.hpp b/src/wmtk/Mesh.hpp
@@ -414,6 +414,15 @@ class Mesh : public std::enable_shared_from_this<Mesh>
     // hashes for top level simplices (i.e cells) to identify whether tuples
     // are invalid or not
     MeshAttributeHandle<long> m_cell_hash_handle;
+
+
+    /**
+     * Generate a vector of Tuples from global vertex/edge/triangle/tetrahedron index
+     * @param type the type of tuple, can be vertex/edge/triangle/tetrahedron
+     * @param include_deleted if true returns also the deleted tuples (default false)
+     * @return vector of Tuples referring to each type
+     */
+    std::vector<Tuple> get_all(PrimitiveType type, const bool include_deleted) const;
 };
 
 

diff --git a/src/wmtk/io/ParaviewWriter.cpp b/src/wmtk/io/ParaviewWriter.cpp
@@ -84,30 +84,35 @@ ParaviewWriter::ParaviewWriter(
     for (int i = 0; i < 4; ++i) {
         auto pt = PrimitiveType(i);
         if (m_enabled[i]) {
-            const auto tuples = mesh.get_all(pt);
+            // include deleted tuples so that attributes are aligned
+            const auto tuples = mesh.get_all(pt, true);
             cells[i].resize(tuples.size(), i + 1);
 
             for (size_t j = 0; j < tuples.size(); ++j) {
                 const auto& t = tuples[j];
-                long vid = mesh.id(t, PrimitiveType::Vertex);
-                cells[i](j, 0) = vid;
-                if (i > 0) {
-                    auto t1 = mesh.switch_tuple(t, PrimitiveType::Vertex);
-
-                    cells[i](j, 1) = mesh.id(t1, PrimitiveType::Vertex);
-                }
-                if (i > 1) {
-                    auto t1 = mesh.switch_tuple(t, PrimitiveType::Edge);
-                    auto t2 = mesh.switch_tuple(t1, PrimitiveType::Vertex);
-
-                    cells[i](j, 2) = mesh.id(t2, PrimitiveType::Vertex);
-                }
-                if (i > 2) {
-                    auto t1 = mesh.switch_tuple(t, PrimitiveType::Face);
-                    auto t2 = mesh.switch_tuple(t1, PrimitiveType::Edge);
-                    auto t3 = mesh.switch_tuple(t2, PrimitiveType::Vertex);
-
-                    cells[i](j, 3) = mesh.id(t3, PrimitiveType::Vertex);
+                if (t.is_null()) {
+                    for (int d = 0; d < i; ++d) cells[i](j, d) = 0;
+                } else {
+                    long vid = mesh.id(t, PrimitiveType::Vertex);
+                    cells[i](j, 0) = vid;
+                    if (i > 0) {
+                        auto t1 = mesh.switch_tuple(t, PrimitiveType::Vertex);
+
+                        cells[i](j, 1) = mesh.id(t1, PrimitiveType::Vertex);
+                    }
+                    if (i > 1) {
+                        auto t1 = mesh.switch_tuple(t, PrimitiveType::Edge);
+                        auto t2 = mesh.switch_tuple(t1, PrimitiveType::Vertex);
+
+                        cells[i](j, 2) = mesh.id(t2, PrimitiveType::Vertex);
+                    }
+                    if (i > 2) {
+                        auto t1 = mesh.switch_tuple(t, PrimitiveType::Face);
+                        auto t2 = mesh.switch_tuple(t1, PrimitiveType::Edge);
+                        auto t3 = mesh.switch_tuple(t2, PrimitiveType::Vertex);
+
+                        cells[i](j, 3) = mesh.id(t3, PrimitiveType::Vertex);
+                    }
                 }
             }
         }

diff --git a/tests/test_io.cpp b/tests/test_io.cpp
@@ -6,11 +6,22 @@
 #include <wmtk/io/ParaviewWriter.hpp>
 #include <wmtk/utils/mesh_utils.hpp>
 
+#include <wmtk/operations/OperationFactory.hpp>
+#include <wmtk/operations/tri_mesh/EdgeSplit.hpp>
+
+#include "tools/DEBUG_TriMesh.hpp"
+#include "tools/TriMesh_examples.hpp"
+
 #include <catch2/catch_test_macros.hpp>
 
 #include <igl/read_triangle_mesh.h>
 
 using namespace wmtk;
+using namespace wmtk::tests;
+
+
+constexpr PrimitiveType PV = PrimitiveType::Vertex;
+constexpr PrimitiveType PE = PrimitiveType::Edge;
 
 TEST_CASE("hdf5_2d", "[io]")
 {
@@ -86,6 +97,69 @@ TEST_CASE("paraview_3d", "[io]")
     mesh.serialize(writer);
 }
 
+TEST_CASE("attribute_after_split", "[io]")
+{
+    DEBUG_TriMesh m = single_triangle_with_position();
+    wmtk::MeshAttributeHandle<long> attribute_handle =
+        m.register_attribute<long>(std::string("test_attribute"), PE, 1);
+    wmtk::MeshAttributeHandle<double> pos_handle =
+        m.get_attribute_handle<double>(std::string("position"), PV);
+
+    {
+        Accessor<long> acc_attribute = m.create_accessor<long>(attribute_handle);
+        Accessor<double> acc_pos = m.create_accessor<double>(pos_handle);
+
+        const Tuple edge = m.edge_tuple_between_v1_v2(0, 1, 0);
+
+        Eigen::Vector3d p_mid;
+        {
+            const Eigen::Vector3d p0 = acc_pos.vector_attribute(edge);
+            const Eigen::Vector3d p1 = acc_pos.vector_attribute(m.switch_vertex(edge));
+            p_mid = 0.5 * (p0 + p1);
+        }
+
+        {
+            // set edge(0,1)'s tag as 1
+            acc_attribute.scalar_attribute(edge) = 1;
+
+            // all edges hold 0 besides "edge"
+            for (const Tuple& t : m.get_all(PE)) {
+                if (m.simplices_are_equal(Simplex::edge(edge), Simplex::edge(t))) {
+                    CHECK(acc_attribute.scalar_attribute(t) == 1);
+                } else {
+                    CHECK(acc_attribute.scalar_attribute(t) == 0);
+                }
+            }
+
+            wmtk::operations::OperationSettings<operations::tri_mesh::EdgeSplit> op_settings;
+            op_settings.split_boundary_edges = true;
+            op_settings.initialize_invariants(m);
+
+            operations::tri_mesh::EdgeSplit op(m, edge, op_settings);
+            REQUIRE(op());
+
+            // set new vertex position
+            acc_pos.vector_attribute(op.return_tuple()) = p_mid;
+        }
+
+        // since the default value is 0, there should be no other value in this triangle
+        for (const Tuple& t : m.get_all(PE)) {
+            CHECK(acc_attribute.scalar_attribute(t) == 0);
+        }
+    } // end of scope for the accessors
+
+    {
+        Accessor<long> acc_attribute = m.create_accessor<long>(attribute_handle);
+        for (const Tuple& t : m.get_all(PE)) {
+            CHECK(acc_attribute.scalar_attribute(t) == 0);
+        }
+    }
+
+    // attribute_after_split_edges.hdf contains a 1 in the "test_attribute"
+    ParaviewWriter writer("attribute_after_split", "position", m, true, true, true, false);
+    m.serialize(writer);
+}
+
 // TEST_CASE("io", "[io][mshio]")
 // {
 //     using namespace wmtk;