Add MSH writer.

src/wmtk/io/CMakeLists.txt

@@ -12,6 +12,8 @@ set(SRC_FILES
     MeshReader.cpp
     MshReader.hpp
     MshReader.cpp
+    MshWriter.hpp
+    MshWriter.cpp
     ParaviewWriter.hpp
     ParaviewWriter.cpp
 

src/wmtk/io/MshReader.cpp

@@ -297,6 +297,7 @@ void MshReader::extract_element_attribute(
         assert(data.entries.size() == tuples.size());
 
         const int64_t attr_dim = data.header.int_tags[1];
+        assert(data.header.int_tags[2] == tuples.size());
 
         auto tag_handle = m.register_attribute<double>(attr_name, primitive_type, attr_dim);
         auto acc = m.create_accessor<double>(tag_handle);
@@ -317,8 +318,9 @@ void MshReader::extract_element_attribute(
     }
 }
 
-auto MshReader::generate(const std::optional<std::vector<std::vector<std::string>>>&
-                             extra_attributes_opt) -> std::shared_ptr<Mesh>
+auto MshReader::generate(
+    const std::optional<std::vector<std::vector<std::string>>>& extra_attributes_opt)
+    -> std::shared_ptr<Mesh>
 {
     std::shared_ptr<Mesh> res;
     switch (get_mesh_dimension()) {

src/wmtk/io/MshWriter.cpp

@@ -0,0 +1,140 @@
+#include "MshWriter.hpp"
+
+#include <numeric>
+#include <wmtk/utils/Logger.hpp>
+
+namespace wmtk::io {
+
+void MshWriter::write(
+    const std::filesystem::path& m_name,
+    const Mesh& mesh,
+    const std::string& position_attribute_name,
+    const std::vector<std::string>& cell_attribute_names)
+{
+    auto pos_handle =
+        mesh.get_attribute_handle<double>(position_attribute_name, PrimitiveType::Vertex);
+    auto pos_acc = mesh.create_const_accessor<double>(pos_handle);
+
+    const auto vertices = mesh.get_all_id_simplex(PrimitiveType::Vertex);
+    const auto cells = mesh.get_all_id_simplex(mesh.top_simplex_type());
+
+    const int64_t dim = (int64_t)mesh.top_simplex_type();
+    const int64_t pos_dim = pos_handle.dimension();
+
+    mshio::MshSpec spec;
+
+    auto& format = spec.mesh_format;
+    format.version = "4.1"; // Only version "2.2" and "4.1" are supported.
+    format.file_type = 1; // 0: ASCII, 1: binary.
+    format.data_size = sizeof(size_t); // Size of data, defined as sizeof(size_t) = 8.
+
+    // nodes
+    {
+        auto& nodes = spec.nodes;
+        nodes.num_entity_blocks = 1; // Number of node blocks.
+        nodes.num_nodes = vertices.size(); // Total number of nodes.
+        nodes.min_node_tag = 1;
+        nodes.max_node_tag = vertices.size();
+        //nodes.entity_blocks = {...}; // A std::vector of node blocks.
+        nodes.entity_blocks.push_back({});
+
+        auto& block = nodes.entity_blocks[0];
+        block.entity_dim = pos_dim; // The dimension of the entity.
+        block.entity_tag = 1; // The entity these nodes belongs to.
+        block.parametric = 0; // 0: non-parametric, 1: parametric.
+        block.num_nodes_in_block = vertices.size(); // The number of nodes in block.
+        //block.tags = {...}; // A std::vector of unique, positive node tags.
+        block.tags.resize(vertices.size());
+        std::iota(block.tags.begin(), block.tags.end(), 1);
+        //block.data = {...}; // A std::vector of coordinates (x,y,z,<u>,<v>,<w>,...)
+
+        block.data.resize(vertices.size() * pos_dim);
+        for (const auto& v : vertices) {
+            const int64_t idx = mesh.id(v);
+            const auto p = pos_acc.const_vector_attribute(v);
+            for (int64_t i = 0; i < pos_dim; ++i) {
+                block.data[idx * pos_dim + i] = p(i);
+            }
+        }
+    }
+
+    // elements
+    {
+        auto& elements = spec.elements;
+        elements.num_entity_blocks = 1; // Number of element blocks.
+        elements.num_elements = cells.size(); // Total number of elmeents.
+        elements.min_element_tag = 1;
+        elements.max_element_tag = cells.size();
+        //elements.entity_blocks = {...}; // A std::vector of element blocks.
+        elements.entity_blocks.push_back({});
+
+        auto& block = elements.entity_blocks[0];
+        block.entity_dim = dim; // The dimension of the elements.
+        block.entity_tag = 1; // The entity these elements belongs to.
+        switch (dim) {
+        case 1: // edge
+            block.element_type = 1; // See element type table of mshio.
+            break;
+        case 2: // triangle
+            block.element_type = 2; // See element type table of mshio.
+            break;
+        case 3: // tet
+            block.element_type = 4; // See element type table of mshio.
+            break;
+        default: log_and_throw_error("Element type not supported for writing to MSH.");
+        }
+        block.num_elements_in_block = cells.size(); // The number of elements in this block.
+        //block.data = {...}; // See more detail below.
+
+        block.data.resize(cells.size() * (dim + 2));
+        for (const auto& c : cells) {
+            const int64_t idx = mesh.id(c);
+            const auto vs = mesh.orient_vertices(mesh.get_tuple_from_id_simplex(c));
+            assert(vs.size() == dim + 1);
+
+            // block data holds entries for every cell: element_tag node_tag ... node_tag
+            block.data[idx * (dim + 2) + 0] = idx + 1;
+            for (int64_t i = 0; i < vs.size(); ++i) {
+                const int64_t vid = mesh.id(vs[i], PrimitiveType::Vertex);
+                block.data[idx * (dim + 2) + 1 + i] = vid + 1;
+            }
+        }
+    }
+
+    // element attributes
+    for (const std::string& attribute_name : cell_attribute_names) {
+        const auto attr =
+            mesh.get_attribute_handle_typed<double>(attribute_name, mesh.top_simplex_type());
+        const auto acc = mesh.create_const_accessor(attr);
+
+        spec.element_data.push_back({});
+        auto& data = spec.element_data.back();
+
+        auto& header = data.header;
+        header.string_tags.push_back(attribute_name); // [field_name, <interpolation_scheme>, ...]
+        header.real_tags.push_back(0); // [<time value>, ...]
+        header.int_tags.resize(3); // [time step, num fields, num entities, <partition id>, ...]
+        header.int_tags[0] = 0;
+        header.int_tags[1] = acc.dimension();
+        header.int_tags[2] = cells.size();
+
+        data.entries.resize(cells.size());
+        for (const auto& c : cells) {
+            const int64_t idx = mesh.id(c);
+            auto& entry = data.entries[idx];
+
+            entry.tag = idx + 1;
+            // entry.num_nodes_per_element = dim + 1;
+            entry.data.resize(acc.dimension());
+
+            auto a = acc.const_vector_attribute(c);
+            for (int64_t i = 0; i < acc.dimension(); ++i) {
+                entry.data[i] = a(i);
+            }
+        }
+    }
+
+    mshio::save_msh(m_name.string(), spec);
+}
+
+} // namespace wmtk::io

src/wmtk/io/MshWriter.hpp

@@ -0,0 +1,23 @@
+#pragma once
+
+#include <mshio/mshio.h>
+#include <filesystem>
+#include <wmtk/Mesh.hpp>
+
+namespace wmtk::io {
+
+class MshWriter
+{
+public:
+    MshWriter() = delete;
+
+    static void write(
+        const std::filesystem::path& m_name,
+        const Mesh& mesh,
+        const std::string& position_attribute_name,
+        const std::vector<std::string>& cell_attribute_names = {});
+
+private:
+};
+
+} // namespace wmtk::io

tests/io/test_io.cpp

@@ -5,11 +5,12 @@
 #include <wmtk/io/HDF5Writer.hpp>
 #include <wmtk/io/MeshReader.hpp>
 #include <wmtk/io/MshReader.hpp>
+#include <wmtk/io/MshWriter.hpp>
 #include <wmtk/io/ParaviewWriter.hpp>
+#include <wmtk/multimesh/same_simplex_dimension_surjection.hpp>
 #include <wmtk/utils/Rational.hpp>
-#include <wmtk/utils/mesh_utils.hpp>
 #include <wmtk/utils/merkle_tree_diff.hpp>
-#include <wmtk/multimesh/same_simplex_dimension_surjection.hpp>
+#include <wmtk/utils/mesh_utils.hpp>
 
 #include <wmtk/operations/EdgeSplit.hpp>
 
@@ -97,8 +98,8 @@ TEST_CASE("paraview_2d", "[io]")
 TEST_CASE("paraview_2d_vtag", "[io]")
 {
     auto mesh = read_mesh(WMTK_DATA_DIR "/fan.msh");
-    mesh->register_attribute<int64_t>("tag", PrimitiveType::Triangle, 1);
-    mesh->register_attribute<int64_t>("tag1", PrimitiveType::Vertex, 1);
+    auto a1 = mesh->register_attribute<int64_t>("tag", PrimitiveType::Triangle, 1);
+    auto a2 = mesh->register_attribute<int64_t>("tag1", PrimitiveType::Vertex, 1);
 
     ParaviewWriter writer("paraview_2d_vtag", "vertices", *mesh, true, true, true, false);
     CHECK_NOTHROW(mesh->serialize(writer));
@@ -308,3 +309,56 @@ TEST_CASE("attribute_after_split", "[io][.]")
     ParaviewWriter writer("attribute_after_split", "vertices", m, true, true, true, false);
     CHECK_NOTHROW(m.serialize(writer));
 }
+
+TEST_CASE("msh_write", "[io]")
+{
+    std::shared_ptr<Mesh> m1_ptr;
+    REQUIRE_NOTHROW(m1_ptr = read_mesh(wmtk_data_dir / "sphere_delaunay.msh"));
+    Mesh& m1 = *m1_ptr;
+
+    auto a1 = m1.register_attribute<double>("a", m1.top_simplex_type(), 1);
+    {
+        auto pos_attr = m1.get_attribute_handle_typed<double>("vertices", PrimitiveType::Vertex);
+        auto pos_acc = m1.create_accessor(pos_attr);
+
+        auto a_acc = m1.create_accessor<double>(a1);
+
+        for (const Tuple& t : m1.get_all(m1.top_simplex_type())) {
+            auto p = pos_acc.const_vector_attribute(t);
+            a_acc.scalar_attribute(t) = p(0);
+        }
+    }
+
+    io::Cache cache("wmtk_test", ".");
+    const std::filesystem::path test_file = cache.get_cache_path() / "msh_write.msh";
+
+    REQUIRE_NOTHROW(io::MshWriter::write(test_file, m1, "vertices", {"a"}));
+
+    std::shared_ptr<Mesh> m2_ptr;
+    REQUIRE_NOTHROW(m2_ptr = read_mesh(test_file, false, {"a"}));
+    Mesh& m2 = *m2_ptr;
+
+    const auto v1 = m1.get_all(PrimitiveType::Vertex);
+    const auto v2 = m2.get_all(PrimitiveType::Vertex);
+    const bool v_same = v1 == v2;
+    CHECK(v_same);
+    const auto c1 = m1.get_all(m1.top_simplex_type());
+    const auto c2 = m2.get_all(m2.top_simplex_type());
+    const bool c_same = c1 == c2;
+    CHECK(c_same);
+
+    auto a2 = m2.get_attribute_handle<double>("a", m2.top_simplex_type());
+    REQUIRE(a1.dimension() == a2.dimension());
+    auto acc1 = m1.create_const_accessor<double>(a1);
+    auto acc2 = m2.create_const_accessor<double>(a2);
+    bool a_is_same = true;
+    for (const Tuple& t : c1) {
+        const auto x1 = acc1.const_vector_attribute(t);
+        const auto x2 = acc2.const_vector_attribute(t);
+        if (x1 != x2) {
+            a_is_same = false;
+            break;
+        }
+    }
+    CHECK(a_is_same);
+}