Merge remote-tracking branch 'origin/main' into hackathon/edge_mesh

components/wmtk_components/CMakeLists.txt

@@ -1,4 +1,4 @@
-
+add_subdirectory(delaunay)
 add_subdirectory(input)
 add_subdirectory(isotropic_remeshing)
 add_subdirectory(mesh_info)

components/wmtk_components/delaunay/CMakeLists.txt

@@ -0,0 +1,13 @@
+
+set(SRC_FILES
+    internal/DelaunayOptions.hpp
+    internal/delaunay_2d.hpp
+    internal/delaunay_2d.cpp
+    internal/delaunay_3d.hpp
+    internal/delaunay_3d.cpp
+    internal/delaunay_geogram.hpp
+    internal/delaunay_geogram.cpp
+	delaunay.hpp
+    delaunay.cpp
+    )
+target_sources(wildmeshing_components PRIVATE ${SRC_FILES})

components/wmtk_components/delaunay/delaunay.cpp

@@ -0,0 +1,114 @@
+#include "delaunay.hpp"
+
+#include <wmtk/PointMesh.hpp>
+#include <wmtk/TetMesh.hpp>
+#include <wmtk/TriMesh.hpp>
+#include <wmtk/io/HDF5Writer.hpp>
+#include <wmtk/io/MeshReader.hpp>
+#include <wmtk/utils/mesh_utils.hpp>
+
+#include "internal/DelaunayOptions.hpp"
+#include "internal/delaunay_2d.hpp"
+#include "internal/delaunay_3d.hpp"
+
+namespace wmtk {
+namespace components {
+
+template <int D>
+RowVectors<double, D> points_to_rowvectors(PointMesh& point_cloud)
+{
+    auto pts_attr = point_cloud.get_attribute_handle<double>("position", PrimitiveType::Vertex);
+    auto pts_acc = point_cloud.create_accessor(pts_attr);
+
+    const auto vertices = point_cloud.get_all(PrimitiveType::Vertex);
+
+    RowVectors<double, D> vec(vertices.size(), D);
+    size_t i = 0;
+    for (const Tuple& t : vertices) {
+        const auto p = pts_acc.vector_attribute(t);
+        vec.row(i) = p.transpose();
+        ++i;
+    }
+
+    return vec;
+}
+
+template <int D, typename MeshT>
+void delaunay_exec(
+    const internal::DelaunayOptions& options,
+    std::map<std::string, std::filesystem::path>& files)
+{
+    // 2d --> TriMesh
+    // 3d --> TetMesh
+    static_assert(
+        (D == 2 && std::is_same<MeshT, TriMesh>()) || (D == 3 && std::is_same<MeshT, TetMesh>()));
+
+    // input
+    PointMesh point_cloud;
+    {
+        const std::filesystem::path& file = files[options.input];
+        MeshReader reader(file);
+        reader.read(point_cloud);
+    }
+
+    // make sure dimensions fit
+    {
+        auto pts_attr = point_cloud.get_attribute_handle<double>("position", PrimitiveType::Vertex);
+        auto pts_acc = point_cloud.create_accessor(pts_attr);
+        assert(pts_acc.dimension() == options.cell_dimension);
+    }
+
+    if constexpr (D == 2) {
+        throw "not tested for 2d";
+    }
+
+    MeshT mesh;
+    Eigen::MatrixXd vertices;
+    Eigen::MatrixXi faces;
+    const auto pts_vec = points_to_rowvectors<D>(point_cloud);
+    if constexpr (D == 2) {
+        std::tie(vertices, faces) = internal::delaunay_2d(pts_vec);
+    } else if constexpr (D == 3) {
+        std::tie(vertices, faces) = internal::delaunay_3d(pts_vec);
+    } else {
+        throw "unsupported cell dimension in delaunay component";
+    }
+
+    mesh.initialize(faces.cast<long>());
+    mesh_utils::set_matrix_attribute(vertices, "position", PrimitiveType::Vertex, mesh);
+
+    // output
+    {
+        const std::filesystem::path cache_dir = "cache";
+        const std::filesystem::path cached_mesh_file = cache_dir / (options.output + ".hdf5");
+
+        HDF5Writer writer(cached_mesh_file);
+        mesh.serialize(writer);
+
+        files[options.output] = cached_mesh_file;
+    }
+}
+
+void delaunay(const nlohmann::json& j, std::map<std::string, std::filesystem::path>& files)
+{
+    using namespace internal;
+
+    DelaunayOptions options = j.get<DelaunayOptions>();
+
+    // delaunay
+    switch (options.cell_dimension) {
+    case 2: {
+        delaunay_exec<2, TriMesh>(options, files);
+        break;
+    }
+    case 3: {
+        delaunay_exec<3, TetMesh>(options, files);
+        break;
+    }
+    default: {
+        throw "unsupported cell dimension in delaunay component";
+    }
+    }
+}
+} // namespace components
+} // namespace wmtk

components/wmtk_components/delaunay/delaunay.hpp

@@ -0,0 +1,9 @@
+#pragma once
+#include <map>
+#include <nlohmann/json.hpp>
+
+namespace wmtk::components {
+
+void delaunay(const nlohmann::json& j, std::map<std::string, std::filesystem::path>& files);
+
+} // namespace wmtk::components

components/wmtk_components/delaunay/internal/DelaunayOptions.hpp

@@ -0,0 +1,22 @@
+#pragma once
+
+#include <spdlog/spdlog.h>
+#include <nlohmann/json.hpp>
+
+namespace wmtk {
+namespace components {
+namespace internal {
+
+struct DelaunayOptions
+{
+    std::string type;
+    std::string input;
+    std::string output;
+    long cell_dimension = -1;
+};
+
+NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(DelaunayOptions, type, input, output, cell_dimension);
+
+} // namespace internal
+} // namespace components
+} // namespace wmtk

components/wmtk_components/delaunay/internal/delaunay_2d.cpp

@@ -0,0 +1,13 @@
+#include "delaunay_2d.hpp"
+#include <wmtk/Types.hpp>
+#include "delaunay_geogram.hpp"
+
+
+namespace wmtk::components::internal {
+
+std::tuple<Eigen::MatrixXd, Eigen::MatrixXi> delaunay_2d(Eigen::Ref<const RowVectors2d> points)
+{
+    return delaunay_geogram(points);
+}
+
+} // namespace wmtk::components::internal

components/wmtk_components/delaunay/internal/delaunay_2d.hpp

@@ -0,0 +1,7 @@
+#pragma once
+
+#include <wmtk/Types.hpp>
+
+namespace wmtk::components::internal {
+std::tuple<Eigen::MatrixXd, Eigen::MatrixXi> delaunay_2d(Eigen::Ref<const RowVectors2d> points);
+} // namespace wmtk::components::internal

components/wmtk_components/delaunay/internal/delaunay_3d.cpp

@@ -0,0 +1,12 @@
+#include "delaunay_3d.hpp"
+#include "delaunay_geogram.hpp"
+
+
+namespace wmtk::components::internal {
+
+std::tuple<Eigen::MatrixXd, Eigen::MatrixXi> delaunay_3d(Eigen::Ref<const RowVectors3d> points)
+{
+    return delaunay_geogram(points);
+}
+
+} // namespace wmtk::components::internal

components/wmtk_components/delaunay/internal/delaunay_3d.hpp

@@ -0,0 +1,7 @@
+#pragma once
+
+#include <wmtk/Types.hpp>
+
+namespace wmtk::components::internal {
+std::tuple<Eigen::MatrixXd, Eigen::MatrixXi> delaunay_3d(Eigen::Ref<const RowVectors3d> points);
+} // namespace wmtk::components::internal

components/wmtk_components/delaunay/internal/delaunay_geogram.cpp

@@ -0,0 +1,89 @@
+#include "delaunay_geogram.hpp"
+
+#include <cassert>
+#include <mutex>
+#include <wmtk/utils/Logger.hpp>
+
+#include <geogram/delaunay/delaunay.h>
+
+namespace wmtk::components::internal {
+
+std::tuple<Eigen::MatrixXd, Eigen::MatrixXi> delaunay_geogram(
+    Eigen::Ref<const Eigen::MatrixXd> points)
+{
+    Eigen::MatrixXd vertices;
+    Eigen::MatrixXi tetrahedra;
+
+    if (points.rows() == 0) {
+        wmtk::logger().warn("Cannot run delaunay on an empty point set");
+        return {vertices, tetrahedra};
+    }
+
+    long dim = points.cols();
+    assert(dim == 2 || dim == 3);
+
+    static std::once_flag once_flag;
+    std::call_once(once_flag, []() { GEO::initialize(); });
+
+    GEO::Delaunay_var engine;
+    switch (dim) {
+    case 2: {
+        engine = GEO::Delaunay::create(2, "BDEL2d");
+        break;
+    }
+    case 3: {
+        engine = GEO::Delaunay::create(3, "BDEL");
+        break;
+    }
+    default: throw "unexpected vector type in delaunay_geogram";
+    }
+
+    assert(engine);
+
+    // Some settings
+    engine->set_reorder(false);
+    engine->set_stores_cicl(false); // Incident tetrahedral list.
+    engine->set_stores_neighbors(false); // Vertex neighbors.
+    engine->set_refine(false);
+    engine->set_keeps_infinite(false);
+
+    // Run!
+    geo_assert(points.size() > 0);
+    Eigen::MatrixXd points_transposed = points.transpose();
+    engine->set_vertices(static_cast<GEO::index_t>(points.rows()), points_transposed.data());
+
+    // Extract output.
+    vertices.resize(engine->nb_vertices(), dim);
+    tetrahedra.resize(engine->nb_cells(), dim + 1);
+
+    for (GEO::index_t i = 0; i < vertices.rows(); ++i) {
+        //   depending on branch prediction to make these calls cost nothing
+        switch (dim) {
+        case 2: {
+            vertices.row(i) = Eigen::Map<const Eigen::Vector2d>(engine->vertex_ptr(i)).transpose();
+            break;
+        }
+        case 3: {
+            vertices.row(i) = Eigen::Map<const Eigen::Vector3d>(engine->vertex_ptr(i)).transpose();
+            break;
+        }
+        default: {
+            // for generality allowing for arbitrary dimensions
+            vertices.row(i) =
+                Eigen::Map<const Eigen::VectorXd>(engine->vertex_ptr(i), dim).transpose();
+            break;
+        }
+        }
+    }
+
+    for (GEO::index_t i = 0; i < tetrahedra.rows(); ++i) {
+        for (GEO::index_t j = 0; j < tetrahedra.cols(); ++j) {
+            tetrahedra(i, j) = engine->cell_vertex(i, j);
+        }
+    }
+
+    return {vertices, tetrahedra};
+}
+
+
+} // namespace wmtk::components::internal

components/wmtk_components/delaunay/internal/delaunay_geogram.hpp

@@ -0,0 +1,9 @@
+#pragma once
+
+#include <wmtk/Types.hpp>
+
+namespace wmtk::components::internal {
+
+std::tuple<Eigen::MatrixXd, Eigen::MatrixXi> delaunay_geogram(
+    Eigen::Ref<const Eigen::MatrixXd> points);
+} // namespace wmtk::components::internal