mesh_decimation

components/CMakeLists.txt

@@ -58,6 +58,7 @@ add_component(${WMTK_COMPONENT_PREFIX} "procedural")
 add_component(${WMTK_COMPONENT_PREFIX} "fusion")
 add_component(${WMTK_COMPONENT_PREFIX} "triangle_insertion")
 add_component(${WMTK_COMPONENT_PREFIX} "to_points")
+add_component(${WMTK_COMPONENT_PREFIX} "mesh_decimation")
 add_component(${WMTK_COMPONENT_PREFIX} "winding_number")
 
 string(LENGTH ${json_components} json_components_length)

components/tests/CMakeLists.txt

@@ -17,7 +17,6 @@ wmtk::components)
 target_include_directories(${WMTK_COMPONENT_TEST_TARGET} PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/..)
 target_compile_definitions(${WMTK_COMPONENT_TEST_TARGET} PRIVATE WMTK_TEST_DIR=\"${CMAKE_CURRENT_SOURCE_DIR}\")
 
-
 FetchContent_GetProperties(catch2)
 list(APPEND CMAKE_MODULE_PATH ${catch2_SOURCE_DIR}/extras)
 include(Catch)

components/tests/integration_test.cpp

@@ -88,9 +88,10 @@ int authenticate_json(const std::string& json_file, const bool compute_validatio
         auto tetrahedra = in_args["tests"]["tetrahedra"];
         if (meshes.size() != vertices.size() || meshes.size() != edges.size() ||
             meshes.size() != faces.size() || meshes.size() != tetrahedra.size()) {
-            spdlog::error("JSON size missmatch between meshes and vertices edges faces or "
-                          "tetrahedra or meshes key. Add a * "
-                          "to the beginning of filename to allow appends.");
+            spdlog::error(
+                "JSON size missmatch between meshes and vertices edges faces or "
+                "tetrahedra or meshes key. Set true for the parameter \"compute_validation\""
+                "to allow appends.");
             return 2;
         }
 
@@ -197,6 +198,7 @@ WMTK_INTEGRATION("isotropic_remeshing", false);
 WMTK_INTEGRATION("isotropic_remeshing_mm", false);
 WMTK_INTEGRATION("disk_fan_mm", false);
 // WMTK_INTEGRATION("grid",false);
+WMTK_INTEGRATION("mesh_decimation", false);
 // WMTK_INTEGRATION("wildmeshing_2d", false);
 // WMTK_INTEGRATION("wildmeshing_3d", false);
 WMTK_INTEGRATION("marching", false);

components/tests/test_component_mesh_decimation.cpp

@@ -0,0 +1,87 @@
+#include <catch2/catch_test_macros.hpp>
+#include <nlohmann/json.hpp>
+#include <tools/DEBUG_TetMesh.hpp>
+#include <tools/DEBUG_TriMesh.hpp>
+#include <tools/TetMesh_examples.hpp>
+#include <tools/TriMesh_examples.hpp>
+#include <wmtk/Mesh.hpp>
+#include <wmtk/TriMesh.hpp>
+#include <wmtk/components/mesh_decimation/internal/MeshDecimation.hpp>
+#include <wmtk/components/mesh_decimation/internal/MeshDecimationOptions.hpp>
+#include <wmtk/components/mesh_decimation/mesh_decimation.hpp>
+#include <wmtk/function/simplex/AMIPS.hpp>
+#include <wmtk/invariants/TodoInvariant.hpp>
+#include <wmtk/io/MeshReader.hpp>
+#include <wmtk/io/ParaviewWriter.hpp>
+#include <wmtk/operations/EdgeCollapse.hpp>
+#include <wmtk/operations/attribute_new/CollapseNewAttributeStrategy.hpp>
+#include <wmtk/simplex/link.hpp>
+#include <wmtk/utils/mesh_utils.hpp>
+
+using json = nlohmann::json;
+using namespace wmtk;
+
+const std::filesystem::path data_dir = WMTK_DATA_DIR;
+
+TEST_CASE("component_mesh_decimation_options", "[components][mesh_decimation]")
+{
+    using namespace components::internal;
+
+    json o = {
+        {"input", "input_mesh"},
+        {"output", "output_mesh"},
+        {"target_len", 1.0},
+        {"cell_constraint_tag_name", "tag"},
+        {"attributes", {"vertices", "tag"}},
+        {"pass_through", {"dummy"}}};
+
+    CHECK_NOTHROW(o.get<MeshDecimationOptions>());
+}
+
+TEST_CASE("component_mesh_decimation", "[components][2D][3D]")
+{
+    using namespace wmtk::components;
+    wmtk::io::Cache cache("wmtk_cache", ".");
+
+    SECTION("3D")
+    {
+        auto mesh_in = wmtk::read_mesh(data_dir / "unit_test/meshes/sphere_regularized.hdf5");
+        Mesh& mesh = *mesh_in;
+
+        std::vector<wmtk::attribute::MeshAttributeHandle> keep;
+        wmtk::attribute::MeshAttributeHandle constrait_cell_tag_handle =
+            mesh.get_attribute_handle<int64_t>("tag", mesh.top_simplex_type());
+        wmtk::attribute::MeshAttributeHandle pos_handle =
+            mesh.get_attribute_handle<double>("vertices", PrimitiveType::Vertex);
+        keep.emplace_back(constrait_cell_tag_handle);
+        keep.emplace_back(pos_handle);
+        mesh.clear_attributes(keep);
+
+        std::vector<wmtk::attribute::MeshAttributeHandle> pass_though;
+        internal::MeshDecimation md(mesh, constrait_cell_tag_handle, 5, pass_though);
+        md.process();
+
+        cache.write_mesh(mesh, "out3d");
+    }
+
+    SECTION("2D")
+    {
+        auto mesh_in = wmtk::read_mesh(data_dir / "2d/ellipse_layer/ellipse_01_substructure.hdf5");
+        Mesh& mesh = *mesh_in;
+
+        std::vector<wmtk::attribute::MeshAttributeHandle> keep;
+        wmtk::attribute::MeshAttributeHandle constrait_cell_tag_handle =
+            mesh.get_attribute_handle<int64_t>("tag", mesh.top_simplex_type());
+        wmtk::attribute::MeshAttributeHandle pos_handle =
+            mesh.get_attribute_handle<double>("vertices", PrimitiveType::Vertex);
+        keep.emplace_back(constrait_cell_tag_handle);
+        keep.emplace_back(pos_handle);
+        mesh.clear_attributes(keep);
+
+        std::vector<wmtk::attribute::MeshAttributeHandle> pass_though;
+        internal::MeshDecimation md(mesh, constrait_cell_tag_handle, 5, pass_though);
+        md.process();
+
+        cache.write_mesh(mesh, "out2d");
+    }
+}

components/wmtk_components/mesh_decimation/wmtk/components/mesh_decimation/CMakeLists.txt

@@ -0,0 +1,13 @@
+
+set(SRC_FILES
+    internal/MeshDecimationOptions.hpp
+    internal/MeshDecimationOptions.cpp
+    internal/MeshDecimation.hpp
+    internal/MeshDecimation.cpp
+	mesh_decimation.hpp
+    mesh_decimation.cpp
+    )
+
+
+#CURRENT_COMPONENT_LIB_NAME is set from the main cmake
+target_sources(${CURRENT_COMPONENT_LIB_NAME} PRIVATE ${SRC_FILES})

components/wmtk_components/mesh_decimation/wmtk/components/mesh_decimation/internal/MeshDecimation.cpp

@@ -0,0 +1,188 @@
+#include "MeshDecimation.hpp"
+#include "wmtk/Scheduler.hpp"
+#include "wmtk/function/simplex/AMIPS.hpp"
+#include "wmtk/invariants/MultiMeshLinkConditionInvariant.hpp"
+#include "wmtk/invariants/SimplexInversionInvariant.hpp"
+#include "wmtk/invariants/SmallerFunctionInvariant.hpp"
+#include "wmtk/invariants/TodoInvariant.hpp"
+#include "wmtk/operations/EdgeCollapse.hpp"
+#include "wmtk/operations/attribute_new/CollapseNewAttributeStrategy.hpp"
+#include "wmtk/operations/attribute_update/AttributeTransferStrategy.hpp"
+#include "wmtk/utils/Logger.hpp"
+
+namespace wmtk::components::internal {
+
+MeshDecimation::MeshDecimation(
+    Mesh& mesh,
+    attribute::MeshAttributeHandle constrainted_cell_tag_handle,
+    double target_len,
+    const std::vector<attribute::MeshAttributeHandle>& pass_through_attributes)
+    : m_mesh(mesh)
+    , m_constrainted_cell_tag_handle(constrainted_cell_tag_handle)
+    , m_target_len(target_len)
+    , m_pass_through_attributes(pass_through_attributes)
+{}
+
+void MeshDecimation::process()
+{
+    using namespace wmtk::attribute;
+    using namespace wmtk::invariants;
+    constexpr PrimitiveType PV = PrimitiveType::Vertex;
+    constexpr PrimitiveType PE = PrimitiveType::Edge;
+    constexpr PrimitiveType PF = PrimitiveType::Triangle;
+    constexpr PrimitiveType PT = PrimitiveType::Tetrahedron;
+
+    MeshAttributeHandle& cell_tag_handle = m_constrainted_cell_tag_handle;
+    MeshAttributeHandle position = m_mesh.get_attribute_handle<double>("vertices", PV);
+    MeshAttributeHandle todo_edge_handle =
+        m_mesh.register_attribute<int64_t>("mesh_decimation_todo_edge", PE, 1);
+    MeshAttributeHandle todo_vertex_handle =
+        m_mesh.register_attribute<int64_t>("mesh_decimation_todo_vertex", PV, 1);
+    MeshAttributeHandle edge_len_handle =
+        m_mesh.register_attribute<double>("mesh_decimation_edge_len", PE, 1);
+
+    Accessor<int64_t> acc_cell = m_mesh.create_accessor<int64_t>(cell_tag_handle);
+    Accessor<double> acc_pos = m_mesh.create_accessor<double>(position);
+    Accessor<int64_t> acc_edge = m_mesh.create_accessor<int64_t>(todo_edge_handle);
+    Accessor<int64_t> acc_vertex = m_mesh.create_accessor<int64_t>(todo_vertex_handle);
+    Accessor<double> acc_len = m_mesh.create_accessor<double>(edge_len_handle);
+
+    // build edge tag to protect to topology of the tagged input
+    switch (m_mesh.top_simplex_type()) {
+    case PF:
+        for (const Tuple& edge : m_mesh.get_all(PE)) {
+            if (m_mesh.is_boundary(PE, edge)) {
+                acc_vertex.scalar_attribute(edge) = 1;
+                acc_vertex.scalar_attribute(m_mesh.switch_tuple(edge, PV)) = 1;
+
+                acc_edge.scalar_attribute(edge) = 1;
+            } else if (
+                acc_cell.scalar_attribute(edge) !=
+                acc_cell.scalar_attribute(m_mesh.switch_tuple(edge, PF))) {
+                acc_vertex.scalar_attribute(edge) = 1;
+                acc_vertex.scalar_attribute(m_mesh.switch_tuple(edge, PV)) = 1;
+
+                acc_edge.scalar_attribute(edge) = 1;
+            }
+        }
+        break;
+    case PT: {
+        for (const Tuple& face : m_mesh.get_all(PF)) {
+            if (m_mesh.is_boundary(PF, face)) {
+                acc_vertex.scalar_attribute(face) = 1;
+                acc_vertex.scalar_attribute(m_mesh.switch_tuple(face, PV)) = 1;
+                acc_vertex.scalar_attribute(m_mesh.switch_tuples(face, {PE, PV})) = 1;
+
+                acc_edge.scalar_attribute(face) = 1;
+                acc_edge.scalar_attribute(m_mesh.switch_tuple(face, PE)) = 1;
+                acc_edge.scalar_attribute(m_mesh.switch_tuples(face, {PV, PE})) = 1;
+            } else if (
+                acc_cell.scalar_attribute(face) !=
+                acc_cell.scalar_attribute(m_mesh.switch_tuple(face, PT))) {
+                acc_vertex.scalar_attribute(face) = 1;
+                acc_vertex.scalar_attribute(m_mesh.switch_tuple(face, PV)) = 1;
+                acc_vertex.scalar_attribute(m_mesh.switch_tuples(face, {PE, PV})) = 1;
+
+                acc_edge.scalar_attribute(face) = 1;
+                acc_edge.scalar_attribute(m_mesh.switch_tuple(face, PE)) = 1;
+                acc_edge.scalar_attribute(m_mesh.switch_tuples(face, {PV, PE})) = 1;
+            }
+        }
+        break;
+    }
+    case PE:
+    case PV:
+    default:
+        log_and_throw_error("MeshDecimation.cpp: mesh_decimation component only supports tetmesh "
+                            "and trimesh for now!");
+    }
+
+    // stage 1: tag edges incident to the surface
+    for (const Tuple& edge : m_mesh.get_all(PE)) {
+        if (acc_vertex.scalar_attribute(edge) == 1 ||
+            acc_vertex.scalar_attribute(m_mesh.switch_tuple(edge, PV)) == 1) {
+            acc_edge.scalar_attribute(edge) = 1;
+        }
+    }
+    // stage 2: tag vertices incident to the tagged edge
+    for (const Tuple& edge : m_mesh.get_all(PE)) {
+        if (acc_edge.scalar_attribute(edge) == 1) {
+            acc_vertex.scalar_attribute(edge) = 1;
+            acc_vertex.scalar_attribute(m_mesh.switch_tuple(edge, PV)) = 1;
+        }
+    }
+    // stage 3: tag edges incident to the tagged vertices
+    for (const Tuple& edge : m_mesh.get_all(PE)) {
+        if (acc_vertex.scalar_attribute(edge) == 1 ||
+            acc_vertex.scalar_attribute(m_mesh.switch_tuple(edge, PV)) == 1) {
+            acc_edge.scalar_attribute(edge) = 1;
+        }
+    }
+
+    // Storing edge lengths and Edge length update
+    auto compute_edge_length = [](const Eigen::MatrixXd& P) -> Eigen::VectorXd {
+        assert(P.cols() == 2);
+        assert(P.rows() == 2 || P.rows() == 3);
+        return Eigen::VectorXd::Constant(1, (P.col(0) - P.col(1)).norm());
+    };
+    auto edge_length_update =
+        std::make_shared<wmtk::operations::SingleAttributeTransferStrategy<double, double>>(
+            edge_len_handle,
+            position,
+            compute_edge_length);
+    edge_length_update->run_on_all();
+
+    auto m_prio_short_edges_first = [&](const simplex::Simplex& s) {
+        assert(s.primitive_type() == PrimitiveType::Edge);
+        auto acc = m_mesh.create_accessor<double>(edge_len_handle);
+        return std::vector<double>({acc.scalar_attribute(s.tuple())});
+    };
+
+    auto op_collapse = std::make_shared<operations::EdgeCollapse>(m_mesh);
+
+    op_collapse->add_invariant(
+        std::make_shared<TodoInvariant>(m_mesh, todo_edge_handle.as<int64_t>(), 0));
+    op_collapse->add_invariant(
+        std::make_shared<TodoSmallerInvariant>(m_mesh, edge_len_handle.as<double>(), m_target_len));
+
+    auto m_amips = std::make_shared<function::AMIPS>(m_mesh, position);
+    auto m_link_conditions = std::make_shared<InvariantCollection>(m_mesh);
+    m_link_conditions->add(std::make_shared<MultiMeshLinkConditionInvariant>(m_mesh));
+    auto m_function_invariant =
+        std::make_shared<SmallerFunctionInvariant>(m_mesh.top_simplex_type(), m_amips, 30);
+    auto m_inversion_invariant =
+        std::make_shared<SimplexInversionInvariant>(m_mesh, position.as<double>());
+
+    op_collapse->add_invariant(m_link_conditions);
+    op_collapse->add_invariant(m_inversion_invariant);
+    op_collapse->add_invariant(m_function_invariant);
+
+    op_collapse->add_transfer_strategy(edge_length_update);
+
+    op_collapse->set_priority(m_prio_short_edges_first);
+
+    op_collapse->set_new_attribute_strategy(
+        position,
+        wmtk::operations::CollapseBasicStrategy::Mean);
+    op_collapse->set_new_attribute_strategy(todo_vertex_handle);
+
+    op_collapse->set_new_attribute_strategy(todo_edge_handle);
+    op_collapse->set_new_attribute_strategy(edge_len_handle);
+    op_collapse->set_new_attribute_strategy(cell_tag_handle);
+
+    // pass_through
+    for (const auto& attr : m_pass_through_attributes) {
+        op_collapse->set_new_attribute_strategy(attr);
+    }
+
+    while (true) {
+        Scheduler scheduler;
+        SchedulerStats pass_stats = scheduler.run_operation_on_all(*op_collapse);
+        m_mesh.consolidate();
+        if (pass_stats.number_of_successful_operations() == 0) {
+            break;
+        }
+    }
+}
+
+} // namespace wmtk::components::internal

components/wmtk_components/mesh_decimation/wmtk/components/mesh_decimation/internal/MeshDecimation.hpp

@@ -0,0 +1,27 @@
+#pragma once
+
+#include <wmtk/Mesh.hpp>
+
+namespace wmtk::components::internal {
+
+class MeshDecimation
+{
+public:
+    MeshDecimation(
+        Mesh& mesh,
+        attribute::MeshAttributeHandle constrainted_cell_tag_handle,
+        double target_len,
+        const std::vector<attribute::MeshAttributeHandle>& pass_through_attributes);
+
+    void process();
+
+private:
+    Mesh& m_mesh;
+
+    attribute::MeshAttributeHandle m_constrainted_cell_tag_handle;
+    double m_target_len;
+
+    std::vector<attribute::MeshAttributeHandle> m_pass_through_attributes;
+};
+
+} // namespace wmtk::components::internal

components/wmtk_components/mesh_decimation/wmtk/components/mesh_decimation/internal/MeshDecimationOptions.cpp

@@ -0,0 +1,28 @@
+#include "MeshDecimationOptions.hpp"
+
+#include <wmtk/utils/Logger.hpp>
+
+namespace wmtk::components::internal {
+
+void to_json(nlohmann::json& j, MeshDecimationOptions& o)
+{
+    j = {
+        {"input", o.input},
+        {"output", o.output},
+        {"target_len", o.target_len},
+        {"cell_constraint_tag_name", o.cell_constraint_tag_name},
+        {"attributes", o.attributes},
+        {"pass_through", o.pass_through}};
+}
+
+void from_json(const nlohmann::json& j, MeshDecimationOptions& o)
+{
+    o.input = j.at("input");
+    o.output = j.at("output");
+    o.target_len = j.at("target_len");
+    o.cell_constraint_tag_name = j.at("cell_constraint_tag_name");
+    j.at("attributes").get_to(o.attributes);
+    j.at("pass_through").get_to(o.pass_through);
+}
+
+} // namespace wmtk::components::internal