Merge pull request #837 from wildmeshing/jiacheng/cdt-opt-component

Jiacheng/cdt opt application

applications/CMakeLists.txt

@@ -30,7 +30,7 @@ add_application(isotropic_remeshing OFF)
 add_application(tetwild_simplification ON)
 add_application(triwild ON)
 add_application(tetwild ON)
-add_application(cdt_sec ON)
+add_application(cdt_opt ON)
 add_application(shortest_edge_collapse ON)
 add_application(insertion ON)
 add_application(longest_edge_split ON)

applications/cdt_opt/CMakeLists.txt

@@ -0,0 +1,18 @@
+wmtk_add_application(cdt_opt_app
+    cdt_opt_main.cpp
+    cdt_opt_spec.hpp
+    )
+
+
+
+target_link_libraries(cdt_opt_app PRIVATE
+wmtk::input
+wmtk::multimesh
+wmtk::CDT
+wmtk::wildmeshing
+wmtk::output)
+
+wmtk_register_integration_test(EXEC_NAME cdt_opt_app
+    CONFIG_FILE ${CMAKE_CURRENT_SOURCE_DIR}/cdt_opt_test_config.json
+    GIT_REPOSITORY "https://github.com/wildmeshing/data.git"
+    GIT_TAG d1063ed50de45a1bcae9f8b6ae9b8b1d42885abe)

applications/cdt_sec/cdt_sec_main.cpp → applications/cdt_opt/cdt_opt_main.cpp

@@ -15,10 +15,10 @@
 #include <wmtk/components/input/input.hpp>
 #include <wmtk/components/multimesh/multimesh.hpp>
 #include <wmtk/components/output/output.hpp>
-#include <wmtk/components/shortest_edge_collapse/shortest_edge_collapse.hpp>
+#include <wmtk/components/wildmeshing/wildmeshing.hpp>
 
 
-#include "cdt_sec_spec.hpp"
+#include "cdt_opt_spec.hpp"
 
 using namespace wmtk;
 namespace fs = std::filesystem;
@@ -43,12 +43,12 @@ int main(int argc, char* argv[])
         j = nlohmann::json::parse(ifs);
 
         jse::JSE spec_engine;
-        bool r = spec_engine.verify_json(j, cdt_sec_spec);
+        bool r = spec_engine.verify_json(j, cdt_opt_spec);
         if (!r) {
             wmtk::logger().error("{}", spec_engine.log2str());
             return 1;
         } else {
-            j = spec_engine.inject_defaults(j, cdt_sec_spec);
+            j = spec_engine.inject_defaults(j, cdt_opt_spec);
         }
     }
 
@@ -57,7 +57,7 @@ int main(int argc, char* argv[])
     auto mesh = wmtk::components::input::input(input_file);
     wmtk::logger().info("mesh has {} vertices", mesh->get_all(PrimitiveType::Vertex).size());
 
-    auto mesh_after_cdt = wmtk::components::CDT(static_cast<const TriMesh&>(*mesh), true, false);
+    auto mesh_after_cdt = wmtk::components::CDT(static_cast<const TriMesh&>(*mesh), true, true);
 
     attribute::MeshAttributeHandle mesh_after_cdt_position_handle;
 
@@ -87,36 +87,51 @@ int main(int argc, char* argv[])
     }
 
     std::string output_file = j["output"];
-    wmtk::components::output::output(*parent_mesh, output_file + "_before_sec", "vertices");
-    wmtk::components::output::output(*child_mesh, output_file + "_surface_before_sec", "vertices");
+    wmtk::components::output::output(*parent_mesh, output_file + "_before_opt", "vertices");
+    wmtk::components::output::output(*child_mesh, output_file + "_surface_before_opt", "vertices");
 
     std::vector<attribute::MeshAttributeHandle> pass_through;
     auto boundary_handle =
         parent_mesh->get_attribute_handle<int64_t>("is_boundary", PrimitiveType::Triangle);
     pass_through.push_back(boundary_handle);
 
-    auto child_mesh_position_handle =
-        child_mesh->get_attribute_handle<double>("vertices", PrimitiveType::Vertex);
-
-    attribute::MeshAttributeHandle parent_mesh_position_handle =
-        parent_mesh->get_attribute_handle<double>("vertices", PrimitiveType::Vertex);
-
-    {
-        using namespace components::shortest_edge_collapse;
-
-        ShortestEdgeCollapseOptions options;
-        options.position_handle = child_mesh_position_handle;
-        options.other_position_handles.emplace_back(parent_mesh_position_handle);
-        options.length_rel = j["length_rel"];
-        options.envelope_size = j["envelope_size"];
-        options.check_inversions = true;
-        options.pass_through_attributes = pass_through;
-
-        shortest_edge_collapse(*parent_mesh, options);
-    }
-
-    wmtk::components::output::output(*parent_mesh, output_file, "vertices");
-    wmtk::components::output::output(*child_mesh, output_file + "_surface", "vertices");
+    std::vector<wmtk::components::EnvelopeOptions> enves;
+
+    wmtk::components::EnvelopeOptions e_surface;
+    e_surface.envelope_name = "surface";
+    e_surface.envelope_constrained_mesh = child_mesh;
+    e_surface.envelope_geometry_mesh = child_mesh;
+    e_surface.constrained_position_name = "vertices";
+    e_surface.geometry_position_name = "vertices";
+    e_surface.thickness = j["envelope_size"];
+
+    enves.push_back(e_surface);
+
+    wmtk::components::WildMeshingOptions wmo;
+    wmo.input_mesh = parent_mesh;
+    wmo.input_mesh_position = "vertices";
+    wmo.target_edge_length = j["length_rel"];
+    wmo.target_max_amips = j["target_max_amips"];
+    wmo.max_passes = j["max_passes"];
+    wmo.intermediate_output = false;
+    wmo.replace_double_coordinate = false;
+    wmo.scheduler_update_frequency = 0;
+    wmo.intermediate_output_path = "";
+    wmo.intermediate_output_name = j["output"];
+    wmo.envelopes = enves;
+    wmo.pass_through = pass_through;
+    wmo.skip_split = false;
+    wmo.skip_collapse = false;
+    wmo.skip_swap = false;
+    wmo.skip_smooth = true;
+
+    auto meshes_after_wildmeshing = wildmeshing(wmo);
+
+    wmtk::components::output::output(*meshes_after_wildmeshing[0].first, output_file, "vertices");
+    wmtk::components::output::output(
+        *meshes_after_wildmeshing[1].first,
+        output_file + "_surface",
+        "vertices");
 
     const std::string report = j["report"];
     if (!report.empty()) {

applications/cdt_sec/cdt_sec_spec.hpp → applications/cdt_opt/cdt_opt_spec.hpp

@@ -2,7 +2,7 @@
 #include <nlohmann/json.hpp>
 namespace {
 
-nlohmann::json cdt_sec_spec = R"(
+nlohmann::json cdt_opt_spec = R"(
 [
     {
         "pointer": "/",
@@ -12,7 +12,9 @@ nlohmann::json cdt_sec_spec = R"(
           "root",
           "report",
           "envelope_size",
-          "length_rel"
+          "length_rel",
+          "target_max_amips",
+          "max_passes"
         ]
 },
 {
@@ -29,6 +31,16 @@ nlohmann::json cdt_sec_spec = R"(
     "type": "float",
     "default": 0.1
 },
+{
+    "pointer": "/target_max_amips",
+    "type": "float",
+    "default": 50
+},
+{
+    "pointer": "/max_passes",
+    "type": "int",
+    "default": 10
+},
 {
     "pointer": "/output",
     "type": "string"

applications/cdt_sec/cdt_sec_test_config.json → applications/cdt_opt/cdt_opt_test_config.json

@@ -2,7 +2,7 @@
     "test_directory": "unit_test",
     "tests": [],
     "release_only_tests": [
-        "cdt_sec_out.json"
+        "cdt_opt_out.json"
     ],
     "input_tag": "input",
     "oracle_tag": "report",

components/CDT/wmtk/components/CDT/internal/CDT.cpp

@@ -2,9 +2,14 @@
 
 //
 #include "CDT.hpp"
+#include <wmtk/Scheduler.hpp>
+#include <wmtk/operations/Rounding.hpp>
 #include <wmtk/utils/Logger.hpp>
 #include <wmtk/utils/mesh_utils.hpp>
 #include <wmtk/utils/orient.hpp>
+
+#include <wmtk/invariants/SimplexInversionInvariant.hpp>
+
 #include "cdt_lib.hpp"
 #include "get_vf.hpp"
 
@@ -79,6 +84,27 @@ std::shared_ptr<wmtk::TetMesh> CDT_internal(
         }
 
         mesh_utils::set_matrix_attribute(V, "vertices", PrimitiveType::Vertex, *tm);
+
+        auto rounding_pt_attribute =
+            tm->get_attribute_handle_typed<Rational>("vertices", PrimitiveType::Vertex);
+
+
+        auto rounding = std::make_shared<wmtk::operations::Rounding>(*tm, rounding_pt_attribute);
+        rounding->add_invariant(
+            std::make_shared<SimplexInversionInvariant<Rational>>(*tm, rounding_pt_attribute));
+
+        Scheduler scheduler;
+        auto stats = scheduler.run_operation_on_all(*rounding);
+
+        logger().info(
+            "Executed rounding, {} ops (S/F) {}/{}. Time: collecting: {}, sorting: {}, "
+            "executing: {}",
+            stats.number_of_performed_operations(),
+            stats.number_of_successful_operations(),
+            stats.number_of_failed_operations(),
+            stats.collecting_time,
+            stats.sorting_time,
+            stats.executing_time);
     } else {
         MatrixX<double> V_double;
         V_double.resize(V_final_str.size(), 3);

components/wildmeshing/wmtk/components/wildmeshing/internal/wildmeshing3d.cpp

@@ -226,7 +226,7 @@ std::vector<std::pair<std::shared_ptr<Mesh>, std::string>> wildmeshing3d(
         PrimitiveType::Edge,
         1,
         false,
-        options.target_edge_length); // defaults to target edge length
+        target_edge_length); // defaults to target edge length
 
     // Target edge length update
     const double min_edge_length = [&]() -> double {

src/wmtk/invariants/CMakeLists.txt

@@ -93,6 +93,8 @@ set(SRC_FILES
     EnergyFilterInvariant.hpp
     EnergyFilterInvariant.cpp
 
+    Swap2dEdgeLengthInvariant.hpp
+    Swap2dEdgeLengthInvariant.cpp
     FrozenVertexInvariant.hpp
     FrozenVertexInvariant.cpp
     Swap2dUnroundedVertexInvariant.hpp

src/wmtk/invariants/Swap2dEdgeLengthInvariant.cpp

@@ -0,0 +1,63 @@
+#include "Swap2dEdgeLengthInvariant.hpp"
+#include <wmtk/Mesh.hpp>
+#include <wmtk/function/utils/amips.hpp>
+#include <wmtk/utils/orient.hpp>
+
+namespace wmtk {
+Swap2dEdgeLengthInvariant::Swap2dEdgeLengthInvariant(
+    const Mesh& m,
+    const attribute::TypedAttributeHandle<double>& coordinate)
+    : Invariant(m, true, false, false)
+    , m_coordinate_handle(coordinate)
+{}
+
+bool Swap2dEdgeLengthInvariant::before(const simplex::Simplex& t) const
+{
+    assert(mesh().top_cell_dimension() == 2);
+    constexpr static PrimitiveType PV = PrimitiveType::Vertex;
+    constexpr static PrimitiveType PE = PrimitiveType::Edge;
+    constexpr static PrimitiveType PF = PrimitiveType::Triangle;
+
+    auto accessor = mesh().create_const_accessor(m_coordinate_handle);
+
+    // get the coords of the vertices
+    // input face end points
+    const Tuple v0 = t.tuple();
+    const Tuple v1 = mesh().switch_tuple(v0, PV);
+    // other 2 vertices
+    const Tuple v2 = mesh().switch_tuples(v0, {PE, PV});
+    const Tuple v3 = mesh().switch_tuples(v0, {PF, PE, PV});
+
+    // use length
+    // const double length_old =
+    //     (accessor.const_vector_attribute(v0) - accessor.const_vector_attribute(v1)).norm();
+    // const double length_new =
+    //     (accessor.const_vector_attribute(v2) - accessor.const_vector_attribute(v3)).norm();
+
+    // return length_new > length_old;
+
+    // use height
+
+    Eigen::Vector3d p0 = accessor.const_vector_attribute(v0);
+    Eigen::Vector3d p1 = accessor.const_vector_attribute(v1);
+    Eigen::Vector3d p2 = accessor.const_vector_attribute(v2);
+    Eigen::Vector3d p3 = accessor.const_vector_attribute(v3);
+
+    const double h0 = ((p0 - p2).cross(p0 - p3)).norm() / (p2 - p3).norm();
+    const double h1 = ((p1 - p2).cross(p1 - p3)).norm() / (p2 - p3).norm();
+    const double h2 = ((p2 - p0).cross(p2 - p1)).norm() / (p0 - p1).norm();
+    const double h3 = ((p3 - p0).cross(p3 - p1)).norm() / (p0 - p1).norm();
+
+    const double min_old = std::min(h2, h3);
+    const double min_new = std::min(h0, h1);
+
+    if (((p0 - p1).norm() * h2 / 2) < 1e-10 || ((p0 - p1).norm() * h3 / 2) < 1e-10) {
+        if (((p2 - p3).norm() * h0 / 2) > 1e-10 && ((p2 - p3).norm() * h1 / 2) > 1e-10) {
+            return true;
+        }
+    }
+
+    return min_old < min_new;
+}
+
+} // namespace wmtk

src/wmtk/invariants/Swap2dEdgeLengthInvariant.hpp

@@ -0,0 +1,20 @@
+#pragma once
+
+#include <wmtk/attribute/TypedAttributeHandle.hpp>
+#include "Invariant.hpp"
+
+namespace wmtk {
+class Swap2dEdgeLengthInvariant : public Invariant
+{
+public:
+    Swap2dEdgeLengthInvariant(
+        const Mesh& m,
+        const attribute::TypedAttributeHandle<double>& coordinate);
+    using Invariant::Invariant;
+
+    bool before(const simplex::Simplex& t) const override;
+
+private:
+    const attribute::TypedAttributeHandle<double> m_coordinate_handle;
+};
+} // namespace wmtk