From a19d398c7878ada59e6f3e020fa946b5227599a1 Mon Sep 17 00:00:00 2001 From: Thomas Krijnen Date: Tue, 7 Apr 2026 15:47:58 +0200 Subject: [PATCH] Vibe code an implementation that uses manifold --- cmake/CMakeLists.txt | 17 + src/ifcconvert/IfcConvert.cpp | 5 +- src/ifcgeom/hybrid_kernel.h | 23 +- .../manifold/ManifoldConversionResult.cpp | 534 +++++++++ .../manifold/ManifoldConversionResult.h | 80 ++ .../kernels/manifold/ManifoldKernel.cpp | 1028 +++++++++++++++++ src/ifcgeom/kernels/manifold/ManifoldKernel.h | 38 + win/build-deps.cmd | 35 +- 8 files changed, 1757 insertions(+), 3 deletions(-) create mode 100644 src/ifcgeom/kernels/manifold/ManifoldConversionResult.cpp create mode 100644 src/ifcgeom/kernels/manifold/ManifoldConversionResult.h create mode 100644 src/ifcgeom/kernels/manifold/ManifoldKernel.cpp create mode 100644 src/ifcgeom/kernels/manifold/ManifoldKernel.h diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt index 65ce6e1d0e..86ae1d6143 100644 --- a/cmake/CMakeLists.txt +++ b/cmake/CMakeLists.txt @@ -74,6 +74,7 @@ option(BUILD_PACKAGE "" OFF) option(WITH_OPENCASCADE "Enable geometry interpretation using Open CASCADE" ON) option(WITH_CGAL "Enable geometry interpretation using CGAL" ON) +option(WITH_MANIFOLD "Enable geometry interpretation using Manifold" OFF) option(COLLADA_SUPPORT "Build IfcConvert with COLLADA support (requires OpenCOLLADA)." ON) option(GLTF_SUPPORT "Build IfcConvert with glTF support (requires json.hpp)." OFF) option(HDF5_SUPPORT "Enable HDF5 support (requires HDF5, zlib)" ON) @@ -215,6 +216,22 @@ if(BUILD_IFCGEOM AND WITH_OPENCASCADE) list(APPEND GEOMETRY_KERNELS opencascade) endif() +message(STATUS "BUILD_IFCGEOM WITH_MANIFOLD: ${BUILD_IFCGEOM} ${WITH_MANIFOLD}") + +if(BUILD_IFCGEOM AND WITH_MANIFOLD) + find_package(manifold CONFIG REQUIRED) + if(TARGET manifold::manifold) + set(MANIFOLD_LIBRARIES manifold::manifold) + elseif(TARGET manifold) + set(MANIFOLD_LIBRARIES manifold) + else() + message(FATAL_ERROR "Unable to determine manifold target") + endif() + add_definitions(-DIFOPSH_WITH_MANIFOLD) + set(SWIG_DEFINES ${SWIG_DEFINES} -DIFOPSH_WITH_MANIFOLD) + list(APPEND GEOMETRY_KERNELS manifold) +endif() + if(GLTF_SUPPORT) UNIFY_ENVVARS_AND_CACHE(JSON_INCLUDE_DIR) if(NOT JSON_INCLUDE_DIR) diff --git a/src/ifcconvert/IfcConvert.cpp b/src/ifcconvert/IfcConvert.cpp index 8d3713500b..5907c988b3 100644 --- a/src/ifcconvert/IfcConvert.cpp +++ b/src/ifcconvert/IfcConvert.cpp @@ -277,6 +277,9 @@ int main(int argc, char** argv) { std::string offset_str, rotation_str; std::string default_kernel; +#ifdef IFOPSH_WITH_MANIFOLD + default_kernel = "manifold"; +#endif #ifdef IFOPSH_WITH_CGAL default_kernel = "cgal"; #endif @@ -292,7 +295,7 @@ int main(int argc, char** argv) { po::options_description geom_options("Geometry options"); geom_options.add_options() ("kernel", po::value(&geometry_kernel)->default_value(default_kernel), - "Geometry kernel to use (opencascade, cgal, cgal-simple, hybrid-cgal-simple-opencascade).") + "Geometry kernel to use (opencascade, cgal, cgal-simple, manifold, hybrid-cgal-simple-opencascade).") ("threads,j", po::value(&num_threads)->default_value(1), "Number of parallel processing threads for geometry interpretation.") ("center-model", diff --git a/src/ifcgeom/hybrid_kernel.h b/src/ifcgeom/hybrid_kernel.h index a2bf5ae829..e7e40cd3fd 100644 --- a/src/ifcgeom/hybrid_kernel.h +++ b/src/ifcgeom/hybrid_kernel.h @@ -36,6 +36,10 @@ #undef CgalKernel #endif +#ifdef IFOPSH_WITH_MANIFOLD +#include "../ifcgeom/kernels/manifold/ManifoldKernel.h" +#endif + namespace { inline bool is_valid_for_kernel(const ifcopenshell::geometry::kernels::AbstractKernel* k, const IfcGeom::ConversionResult& shp) { #ifdef IFOPSH_WITH_OPENCASCADE @@ -50,6 +54,11 @@ namespace { if (k->geometry_library() == "cgal") { return dynamic_cast(shp.Shape().get()) != nullptr; } +#endif +#ifdef IFOPSH_WITH_MANIFOLD + if (k->geometry_library() == "manifold") { + return dynamic_cast(shp.Shape().get()) != nullptr; + } #endif return false; } @@ -186,6 +195,12 @@ namespace ifcopenshell { } #endif +#ifdef IFOPSH_WITH_MANIFOLD + if (geometry_library_lower == "manifold") { + return std::make_unique(conv_settings); + } +#endif + if (geometry_library_lower.rfind("hybrid-", 0) == 0) { geometry_library_lower = geometry_library_lower.substr(strlen("hybrid")); std::vector> kernels; @@ -213,6 +228,12 @@ namespace ifcopenshell { kernels.emplace_back(new CgalKernel(conv_settings)); geometry_library_lower = geometry_library_lower.substr(strlen("cgal")); } +#endif +#ifdef IFOPSH_WITH_MANIFOLD + if (geometry_library_lower.find("manifold", 0) == 0) { + kernels.emplace_back(new ManifoldKernel(conv_settings)); + geometry_library_lower = geometry_library_lower.substr(strlen("manifold")); + } #endif if (kernels.size() != n + 1) { throw ifcopenshell::exception("Invalid hybrid kernel " + geometry_library); @@ -236,4 +257,4 @@ namespace ifcopenshell { } } -#endif \ No newline at end of file +#endif diff --git a/src/ifcgeom/kernels/manifold/ManifoldConversionResult.cpp b/src/ifcgeom/kernels/manifold/ManifoldConversionResult.cpp new file mode 100644 index 0000000000..7ac1234335 --- /dev/null +++ b/src/ifcgeom/kernels/manifold/ManifoldConversionResult.cpp @@ -0,0 +1,534 @@ +#include "ManifoldConversionResult.h" + +#include "../../../ifcgeom/IfcGeomRepresentation.h" + +#include + +#include +#include +#include +#include +#include + +using IfcGeom::ConversionResultShape; +using IfcGeom::NumberNativeDouble; +using IfcGeom::OpaqueCoordinate; +using IfcGeom::OpaqueNumber; + +namespace { + using Mesh = manifold::MeshGL64; + + Mesh transform_mesh(const Mesh& mesh, const ifcopenshell::geometry::taxonomy::matrix4& place) { + Mesh result = mesh; + const auto& m = place.ccomponents(); + const bool flip = m.block<3, 3>(0, 0).determinant() < 0.; + for (size_t i = 0; i < mesh.NumVert(); ++i) { + Eigen::Vector4d v( + mesh.vertProperties[i * mesh.numProp + 0], + mesh.vertProperties[i * mesh.numProp + 1], + mesh.vertProperties[i * mesh.numProp + 2], + 1.); + auto v2 = m * v; + result.vertProperties[i * result.numProp + 0] = v2(0); + result.vertProperties[i * result.numProp + 1] = v2(1); + result.vertProperties[i * result.numProp + 2] = v2(2); + } + if (flip) { + for (size_t i = 0; i < mesh.NumTri(); ++i) { + std::swap(result.triVerts[i * 3 + 1], result.triVerts[i * 3 + 2]); + } + } + result.runTransform.clear(); + return result; + } + + std::optional make_manifold(const Mesh& mesh) { + manifold::Manifold solid(mesh); + if (solid.Status() == manifold::Manifold::Error::NoError) { + return solid; + } + return std::nullopt; + } + + manifold::Box mesh_bbox(const Mesh& mesh) { + if (!mesh.NumVert()) { + return {}; + } + manifold::Box box( + manifold::vec3( + mesh.vertProperties[0], + mesh.vertProperties[1], + mesh.vertProperties[2]), + manifold::vec3( + mesh.vertProperties[0], + mesh.vertProperties[1], + mesh.vertProperties[2])); + for (size_t i = 1; i < mesh.NumVert(); ++i) { + box.Union(manifold::vec3( + mesh.vertProperties[i * mesh.numProp + 0], + mesh.vertProperties[i * mesh.numProp + 1], + mesh.vertProperties[i * mesh.numProp + 2])); + } + return box; + } + + double bbox_volume(const manifold::Box& box) { + if (!box.IsFinite()) { + return 0.; + } + const auto size = box.Size(); + return size[0] * size[1] * size[2]; + } + + double triangle_area(const Mesh& mesh, size_t tri) { + auto idx = [&](int corner) { return mesh.triVerts[tri * 3 + corner]; }; + auto point = [&](uint32_t i) { + return Eigen::Vector3d( + mesh.vertProperties[i * mesh.numProp + 0], + mesh.vertProperties[i * mesh.numProp + 1], + mesh.vertProperties[i * mesh.numProp + 2]); + }; + const auto a = point(idx(0)); + const auto b = point(idx(1)); + const auto c = point(idx(2)); + return 0.5 * ((b - a).cross(c - a)).norm(); + } + + double mesh_area(const Mesh& mesh) { + double area = 0.; + for (size_t i = 0; i < mesh.NumTri(); ++i) { + area += triangle_area(mesh, i); + } + return area; + } + + double mesh_volume(const Mesh& mesh) { + double volume = 0.; + for (size_t i = 0; i < mesh.NumTri(); ++i) { + auto idx = [&](int corner) { return mesh.triVerts[i * 3 + corner]; }; + auto point = [&](uint32_t v) { + return Eigen::Vector3d( + mesh.vertProperties[v * mesh.numProp + 0], + mesh.vertProperties[v * mesh.numProp + 1], + mesh.vertProperties[v * mesh.numProp + 2]); + }; + const auto a = point(idx(0)); + const auto b = point(idx(1)); + const auto c = point(idx(2)); + volume += a.dot(b.cross(c)) / 6.; + } + return std::abs(volume); + } + + struct EdgeHash { + size_t operator()(const std::pair& edge) const { + return std::hash()((uint64_t(edge.first) << 32) ^ uint64_t(edge.second)); + } + }; + + std::unordered_map, int, EdgeHash> count_edges(const Mesh& mesh) { + std::unordered_map, int, EdgeHash> edges; + for (size_t i = 0; i < mesh.NumTri(); ++i) { + uint32_t tri[3] = { + mesh.triVerts[i * 3 + 0], + mesh.triVerts[i * 3 + 1], + mesh.triVerts[i * 3 + 2] + }; + for (int j = 0; j < 3; ++j) { + auto a = tri[j]; + auto b = tri[(j + 1) % 3]; + if (a > b) { + std::swap(a, b); + } + edges[{a, b}]++; + } + } + return edges; + } + + double mesh_length(const Mesh& mesh) { + double length = 0.; + auto edges = count_edges(mesh); + for (const auto& pair : edges) { + const auto a = pair.first.first; + const auto b = pair.first.second; + Eigen::Vector3d p( + mesh.vertProperties[a * mesh.numProp + 0], + mesh.vertProperties[a * mesh.numProp + 1], + mesh.vertProperties[a * mesh.numProp + 2]); + Eigen::Vector3d q( + mesh.vertProperties[b * mesh.numProp + 0], + mesh.vertProperties[b * mesh.numProp + 1], + mesh.vertProperties[b * mesh.numProp + 2]); + length += (q - p).norm(); + } + return length; + } + + int mesh_edges(const Mesh& mesh) { + return (int)count_edges(mesh).size(); + } + + ifcopenshell::geometry::ManifoldPart make_part(const manifold::Manifold& solid) { + return { solid.GetMeshGL64(), solid }; + } + + ifcopenshell::geometry::ManifoldPart make_box_part(const manifold::Box& box) { + const auto size = box.Size(); + auto solid = manifold::Manifold::Cube(manifold::vec3(size[0], size[1], size[2]), false).Translate(box.min); + return make_part(solid); + } +} + +ifcopenshell::geometry::ManifoldShape::ManifoldShape(const ManifoldPart& part) + : parts_{ part } {} + +ifcopenshell::geometry::ManifoldShape::ManifoldShape(ManifoldPart&& part) + : parts_{ std::move(part) } {} + +ifcopenshell::geometry::ManifoldShape::ManifoldShape(const std::vector& parts) + : parts_(parts) {} + +ifcopenshell::geometry::ManifoldShape::ManifoldShape(std::vector&& parts) + : parts_(std::move(parts)) {} + +std::optional ifcopenshell::geometry::ManifoldShape::as_manifold() const { + if (parts_.empty()) { + return std::nullopt; + } + std::vector solids; + solids.reserve(parts_.size()); + for (const auto& part : parts_) { + if (!part.solid) { + return std::nullopt; + } + solids.push_back(*part.solid); + } + if (solids.size() == 1) { + return solids.front(); + } + return manifold::Manifold::BatchBoolean(solids, manifold::OpType::Add); +} + +void ifcopenshell::geometry::ManifoldShape::Triangulate(ifcopenshell::geometry::Settings, const ifcopenshell::geometry::taxonomy::matrix4& place, IfcGeom::Representation::Triangulation* t, int item_id, int surface_style_id) const { + for (const auto& part : parts_) { + auto mesh = transform_mesh(part.mesh, place); + std::vector indices(mesh.NumVert()); + for (size_t i = 0; i < mesh.NumVert(); ++i) { + indices[i] = t->addVertex( + item_id, + surface_style_id, + mesh.vertProperties[i * mesh.numProp + 0], + mesh.vertProperties[i * mesh.numProp + 1], + mesh.vertProperties[i * mesh.numProp + 2]); + } + auto edges = count_edges(mesh); + for (size_t i = 0; i < mesh.NumTri(); ++i) { + t->addFace( + item_id, + surface_style_id, + indices[mesh.triVerts[i * 3 + 0]], + indices[mesh.triVerts[i * 3 + 1]], + indices[mesh.triVerts[i * 3 + 2]]); + } + for (const auto& edge : edges) { + if (edge.second == 1) { + t->registerEdge(item_id, indices[edge.first.first], indices[edge.first.second]); + } + } + } +} + +void ifcopenshell::geometry::ManifoldShape::Serialize(const ifcopenshell::geometry::taxonomy::matrix4& place, std::string& result) const { + std::stringstream stream; + stream << std::setprecision(17); + size_t offset = 0; + for (const auto& part : parts_) { + auto mesh = transform_mesh(part.mesh, place); + for (size_t i = 0; i < mesh.NumVert(); ++i) { + stream << "v " + << mesh.vertProperties[i * mesh.numProp + 0] << " " + << mesh.vertProperties[i * mesh.numProp + 1] << " " + << mesh.vertProperties[i * mesh.numProp + 2] << "\n"; + } + for (size_t i = 0; i < mesh.NumTri(); ++i) { + stream << "f " + << mesh.triVerts[i * 3 + 0] + 1 + offset << " " + << mesh.triVerts[i * 3 + 1] + 1 + offset << " " + << mesh.triVerts[i * 3 + 2] + 1 + offset << "\n"; + } + offset += mesh.NumVert(); + } + result = stream.str(); +} + +int ifcopenshell::geometry::ManifoldShape::surface_genus() const { + int genus = 0; + for (const auto& part : parts_) { + if (!part.solid) { + return 0; + } + genus += part.solid->Genus(); + } + return genus; +} + +bool ifcopenshell::geometry::ManifoldShape::is_manifold() const { + return std::all_of(parts_.begin(), parts_.end(), [](const auto& part) { return part.solid.has_value(); }); +} + +int ifcopenshell::geometry::ManifoldShape::num_vertices() const { + size_t total = 0; + for (const auto& part : parts_) { + total += part.mesh.NumVert(); + } + return (int)total; +} + +int ifcopenshell::geometry::ManifoldShape::num_edges() const { + int total = 0; + for (const auto& part : parts_) { + total += part.solid ? (int)part.solid->NumEdge() : mesh_edges(part.mesh); + } + return total; +} + +int ifcopenshell::geometry::ManifoldShape::num_faces() const { + size_t total = 0; + for (const auto& part : parts_) { + total += part.solid ? part.solid->NumTri() : part.mesh.NumTri(); + } + return (int)total; +} + +double ifcopenshell::geometry::ManifoldShape::bounding_box(void*& box_ptr) const { + bool initialized = false; + auto* box = static_cast(box_ptr); + if (!box) { + box = new manifold::Box(); + box_ptr = box; + } + for (const auto& part : parts_) { + auto bbox = part.solid ? part.solid->BoundingBox() : mesh_bbox(part.mesh); + if (!bbox.IsFinite()) { + continue; + } + if (!initialized) { + *box = bbox; + initialized = true; + } else { + box->Union(bbox.min); + box->Union(bbox.max); + } + } + return initialized ? bbox_volume(*box) : 0.; +} + +std::pair, OpaqueCoordinate<3>> ifcopenshell::geometry::ManifoldShape::bounding_box() const { + void* box_ptr = nullptr; + bounding_box(box_ptr); + auto* box = static_cast(box_ptr); + if (!box || !box->IsFinite()) { + delete box; + throw std::runtime_error("Invalid shape"); + } + auto result = std::make_pair( + OpaqueCoordinate<3>( + new NumberNativeDouble(box->min[0]), + new NumberNativeDouble(box->min[1]), + new NumberNativeDouble(box->min[2])), + OpaqueCoordinate<3>( + new NumberNativeDouble(box->max[0]), + new NumberNativeDouble(box->max[1]), + new NumberNativeDouble(box->max[2]))); + delete box; + return result; +} + +void ifcopenshell::geometry::ManifoldShape::set_box(void* box_ptr) { + auto* box = static_cast(box_ptr); + if (!box || !box->IsFinite()) { + throw std::runtime_error("Invalid shape"); + } + parts_ = { make_box_part(*box) }; +} + +OpaqueNumber* ifcopenshell::geometry::ManifoldShape::length() { + double total = 0.; + for (const auto& part : parts_) { + total += mesh_length(part.mesh); + } + return new NumberNativeDouble(total); +} + +OpaqueNumber* ifcopenshell::geometry::ManifoldShape::area() { + double total = 0.; + for (const auto& part : parts_) { + total += part.solid ? part.solid->SurfaceArea() : mesh_area(part.mesh); + } + return new NumberNativeDouble(total); +} + +OpaqueNumber* ifcopenshell::geometry::ManifoldShape::volume() { + double total = 0.; + for (const auto& part : parts_) { + total += part.solid ? part.solid->Volume() : mesh_volume(part.mesh); + } + return new NumberNativeDouble(total); +} + +OpaqueCoordinate<3> ifcopenshell::geometry::ManifoldShape::position() { + throw std::runtime_error("Invalid shape"); +} + +OpaqueCoordinate<3> ifcopenshell::geometry::ManifoldShape::axis() { + throw std::runtime_error("Invalid shape"); +} + +OpaqueCoordinate<4> ifcopenshell::geometry::ManifoldShape::plane_equation() { + throw std::runtime_error("Invalid shape"); +} + +std::vector ifcopenshell::geometry::ManifoldShape::convex_decomposition() { + throw std::runtime_error("Not implemented"); +} + +ConversionResultShape* ifcopenshell::geometry::ManifoldShape::halfspaces() { + throw std::runtime_error("Not implemented"); +} + +ConversionResultShape* ifcopenshell::geometry::ManifoldShape::box() { + void* box_ptr = nullptr; + bounding_box(box_ptr); + auto* box = static_cast(box_ptr); + if (!box || !box->IsFinite()) { + delete box; + throw std::runtime_error("Invalid shape"); + } + auto* result = new ManifoldShape(make_box_part(*box)); + delete box; + return result; +} + +ConversionResultShape* ifcopenshell::geometry::ManifoldShape::solid() { + if (!is_manifold()) { + throw std::runtime_error("Invalid shape"); + } + return new ManifoldShape(parts_); +} + +ConversionResultShape* ifcopenshell::geometry::ManifoldShape::wrap_in_compound() { + return new ManifoldShape(parts_); +} + +std::vector ifcopenshell::geometry::ManifoldShape::vertices() { + throw std::runtime_error("Not implemented"); +} + +std::vector ifcopenshell::geometry::ManifoldShape::edges() { + throw std::runtime_error("Not implemented"); +} + +std::vector ifcopenshell::geometry::ManifoldShape::facets() { + throw std::runtime_error("Not implemented"); +} + +ConversionResultShape* ifcopenshell::geometry::ManifoldShape::add(ConversionResultShape* other) { + auto* rhs = dynamic_cast(other); + if (!rhs) { + throw std::runtime_error("Invalid shape"); + } + auto a = as_manifold(); + auto b = rhs->as_manifold(); + if (!a || !b) { + throw std::runtime_error("Invalid shape"); + } + return new ManifoldShape(make_part(*a + *b)); +} + +ConversionResultShape* ifcopenshell::geometry::ManifoldShape::subtract(ConversionResultShape* other) { + auto* rhs = dynamic_cast(other); + if (!rhs) { + throw std::runtime_error("Invalid shape"); + } + auto a = as_manifold(); + auto b = rhs->as_manifold(); + if (!a || !b) { + throw std::runtime_error("Invalid shape"); + } + return new ManifoldShape(make_part(*a - *b)); +} + +ConversionResultShape* ifcopenshell::geometry::ManifoldShape::intersect(ConversionResultShape* other) { + auto* rhs = dynamic_cast(other); + if (!rhs) { + throw std::runtime_error("Invalid shape"); + } + auto a = as_manifold(); + auto b = rhs->as_manifold(); + if (!a || !b) { + throw std::runtime_error("Invalid shape"); + } + return new ManifoldShape(make_part(*a ^ *b)); +} + +ConversionResultShape* ifcopenshell::geometry::ManifoldShape::concat(ConversionResultShape* other) { + auto* rhs = dynamic_cast(other); + if (!rhs) { + throw std::runtime_error("Invalid shape"); + } + auto parts = parts_; + parts.insert(parts.end(), rhs->parts_.begin(), rhs->parts_.end()); + return new ManifoldShape(std::move(parts)); +} + +void ifcopenshell::geometry::ManifoldShape::map(OpaqueCoordinate<4>&, OpaqueCoordinate<4>&) { + throw std::runtime_error("Not implemented"); +} + +void ifcopenshell::geometry::ManifoldShape::map(const std::vector>&, const std::vector>&) { + throw std::runtime_error("Not implemented"); +} + +ConversionResultShape* ifcopenshell::geometry::ManifoldShape::moved(ifcopenshell::geometry::taxonomy::matrix4::ptr place) const { + std::vector moved_parts; + moved_parts.reserve(parts_.size()); + for (const auto& part : parts_) { + auto mesh = transform_mesh(part.mesh, *place); + auto solid = part.solid ? make_manifold(mesh) : std::nullopt; + if (part.solid && !solid) { + throw std::runtime_error("Failed to transform shape"); + } + moved_parts.push_back({ std::move(mesh), std::move(solid) }); + } + return new ManifoldShape(std::move(moved_parts)); +} + +bool ifcopenshell::geometry::ManifoldShape::surface_area_along_direction(double, const ifcopenshell::geometry::taxonomy::matrix4::ptr& place, double& along_x, double& along_y, double& along_z) const { + along_x = along_y = along_z = 0.; + for (const auto& part : parts_) { + auto mesh = transform_mesh(part.mesh, *place); + for (size_t i = 0; i < mesh.NumTri(); ++i) { + auto point = [&](uint32_t v) { + return Eigen::Vector3d( + mesh.vertProperties[v * mesh.numProp + 0], + mesh.vertProperties[v * mesh.numProp + 1], + mesh.vertProperties[v * mesh.numProp + 2]); + }; + const auto a = point(mesh.triVerts[i * 3 + 0]); + const auto b = point(mesh.triVerts[i * 3 + 1]); + const auto c = point(mesh.triVerts[i * 3 + 2]); + auto n = (b - a).cross(c - a); + const auto norm = n.norm(); + if (norm < 1.e-12) { + continue; + } + const auto area = 0.5 * norm; + n /= norm; + along_x += area * std::abs(n(0)); + along_y += area * std::abs(n(1)); + along_z += area * std::abs(n(2)); + } + } + return true; +} diff --git a/src/ifcgeom/kernels/manifold/ManifoldConversionResult.h b/src/ifcgeom/kernels/manifold/ManifoldConversionResult.h new file mode 100644 index 0000000000..423da07ab6 --- /dev/null +++ b/src/ifcgeom/kernels/manifold/ManifoldConversionResult.h @@ -0,0 +1,80 @@ +#ifndef IFCGEOMMANIFOLDREPRESENTATION_H +#define IFCGEOMMANIFOLDREPRESENTATION_H + +#include + +#include "../../../ifcgeom/ConversionResult.h" +#include "../../../ifcgeom/kernels/ifc_geomlibrary_api.h" + +#include + +namespace ifcopenshell { +namespace geometry { + +struct IFC_GEOMLIBRARY_API ManifoldPart { + manifold::MeshGL64 mesh; + std::optional solid; +}; + +class IFC_GEOMLIBRARY_API ManifoldShape : public IfcGeom::ConversionResultShape { +public: + ManifoldShape() = default; + explicit ManifoldShape(const ManifoldPart& part); + explicit ManifoldShape(ManifoldPart&& part); + explicit ManifoldShape(const std::vector& parts); + explicit ManifoldShape(std::vector&& parts); + + const std::vector& parts() const { return parts_; } + std::optional as_manifold() const; + + virtual void Triangulate(ifcopenshell::geometry::Settings settings, const ifcopenshell::geometry::taxonomy::matrix4& place, IfcGeom::Representation::Triangulation* t, int item_id, int surface_style_id) const; + virtual void Serialize(const ifcopenshell::geometry::taxonomy::matrix4& place, std::string&) const; + + virtual int surface_genus() const; + virtual bool is_manifold() const; + + virtual int num_vertices() const; + virtual int num_edges() const; + virtual int num_faces() const; + + virtual double bounding_box(void*&) const; + virtual std::pair, IfcGeom::OpaqueCoordinate<3>> bounding_box() const; + virtual void set_box(void* b); + + virtual IfcGeom::OpaqueNumber* length(); + virtual IfcGeom::OpaqueNumber* area(); + virtual IfcGeom::OpaqueNumber* volume(); + + virtual IfcGeom::OpaqueCoordinate<3> position(); + virtual IfcGeom::OpaqueCoordinate<3> axis(); + virtual IfcGeom::OpaqueCoordinate<4> plane_equation(); + + virtual std::vector convex_decomposition(); + virtual IfcGeom::ConversionResultShape* halfspaces(); + virtual IfcGeom::ConversionResultShape* box(); + virtual IfcGeom::ConversionResultShape* solid(); + virtual IfcGeom::ConversionResultShape* wrap_in_compound(); + + virtual std::vector vertices(); + virtual std::vector edges(); + virtual std::vector facets(); + + virtual IfcGeom::ConversionResultShape* add(IfcGeom::ConversionResultShape*); + virtual IfcGeom::ConversionResultShape* subtract(IfcGeom::ConversionResultShape*); + virtual IfcGeom::ConversionResultShape* intersect(IfcGeom::ConversionResultShape*); + virtual IfcGeom::ConversionResultShape* concat(IfcGeom::ConversionResultShape*); + + virtual void map(IfcGeom::OpaqueCoordinate<4>& from, IfcGeom::OpaqueCoordinate<4>& to); + virtual void map(const std::vector>& from, const std::vector>& to); + virtual IfcGeom::ConversionResultShape* moved(ifcopenshell::geometry::taxonomy::matrix4::ptr) const; + + virtual bool surface_area_along_direction(double tol, const ifcopenshell::geometry::taxonomy::matrix4::ptr&, double& along_x, double& along_y, double& along_z) const; + +private: + std::vector parts_; +}; + +} +} + +#endif diff --git a/src/ifcgeom/kernels/manifold/ManifoldKernel.cpp b/src/ifcgeom/kernels/manifold/ManifoldKernel.cpp new file mode 100644 index 0000000000..fa865ed75e --- /dev/null +++ b/src/ifcgeom/kernels/manifold/ManifoldKernel.cpp @@ -0,0 +1,1028 @@ +#include "ManifoldKernel.h" + +#include "../../../ifcparse/logger.h" + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace ifcopenshell::geometry; +using namespace ifcopenshell::geometry::kernels; + +namespace { + using Mesh = manifold::MeshGL64; + using Part = ifcopenshell::geometry::ManifoldPart; + + std::string manifold_error_string(manifold::Manifold::Error error) { + switch (error) { + case manifold::Manifold::Error::NoError: + return "no error"; + case manifold::Manifold::Error::NonFiniteVertex: + return "non-finite vertex"; + case manifold::Manifold::Error::NotManifold: + return "not manifold"; + case manifold::Manifold::Error::VertexOutOfBounds: + return "vertex out of bounds"; + case manifold::Manifold::Error::PropertiesWrongLength: + return "properties wrong length"; + case manifold::Manifold::Error::MissingPositionProperties: + return "missing position properties"; + case manifold::Manifold::Error::MergeVectorsDifferentLengths: + return "merge vectors different lengths"; + case manifold::Manifold::Error::MergeIndexOutOfBounds: + return "merge index out of bounds"; + case manifold::Manifold::Error::TransformWrongLength: + return "transform wrong length"; + case manifold::Manifold::Error::RunIndexWrongLength: + return "run index wrong length"; + case manifold::Manifold::Error::FaceIDWrongLength: + return "face id wrong length"; + case manifold::Manifold::Error::InvalidConstruction: + return "invalid construction"; + case manifold::Manifold::Error::ResultTooLarge: + return "result too large"; + } + return "unknown error"; + } + + struct VertexKey { + long long x; + long long y; + long long z; + + bool operator==(const VertexKey& other) const { + return x == other.x && y == other.y && z == other.z; + } + }; + + struct VertexKeyHash { + size_t operator()(const VertexKey& key) const { + auto h = std::hash()(key.x); + h ^= std::hash()(key.y) + 0x9e3779b97f4a7c15ull + (h << 6) + (h >> 2); + h ^= std::hash()(key.z) + 0x9e3779b97f4a7c15ull + (h << 6) + (h >> 2); + return h; + } + }; + + struct MeshBuilder { + double precision; + std::vector vertices; + std::unordered_map vertex_map; + std::vector tri_verts; + std::vector face_ids; + + explicit MeshBuilder(double p) : precision(p > 0. ? p : 1.e-9) {} + + VertexKey key(const Eigen::Vector3d& p) const { + return { + (long long)std::llround(p(0) / precision), + (long long)std::llround(p(1) / precision), + (long long)std::llround(p(2) / precision) + }; + } + + uint64_t add_vertex(const Eigen::Vector3d& p) { + auto entry = vertex_map.find(key(p)); + if (entry != vertex_map.end()) { + return entry->second; + } + auto idx = (uint64_t)vertices.size(); + vertices.push_back(p); + vertex_map.insert({ key(p), idx }); + return idx; + } + + void add_triangle(uint64_t a, uint64_t b, uint64_t c, uint64_t face_id) { + if (a == b || b == c || c == a) { + return; + } + tri_verts.push_back(a); + tri_verts.push_back(b); + tri_verts.push_back(c); + face_ids.push_back(face_id); + } + + Mesh build() const { + Mesh mesh; + mesh.numProp = 3; + mesh.vertProperties.reserve(vertices.size() * 3); + for (const auto& vertex : vertices) { + mesh.vertProperties.push_back(vertex(0)); + mesh.vertProperties.push_back(vertex(1)); + mesh.vertProperties.push_back(vertex(2)); + } + mesh.triVerts = tri_verts; + mesh.faceID = face_ids; + mesh.tolerance = precision; + return mesh; + } + }; + + struct LoopPoint { + Eigen::Vector3d xyz; + manifold::vec2 uv; + }; + + using LoopPolygon = std::vector; + + struct EdgeKey { + uint64_t a; + uint64_t b; + + bool operator==(const EdgeKey& other) const { + return a == other.a && b == other.b; + } + }; + + struct EdgeKeyHash { + size_t operator()(const EdgeKey& key) const { + auto h = std::hash()(key.a); + h ^= std::hash()(key.b) + 0x9e3779b97f4a7c15ull + (h << 6) + (h >> 2); + return h; + } + }; + + struct FaceKey { + uint64_t a; + uint64_t b; + uint64_t c; + + bool operator==(const FaceKey& other) const { + return a == other.a && b == other.b && c == other.c; + } + }; + + struct FaceKeyHash { + size_t operator()(const FaceKey& key) const { + auto h = std::hash()(key.a); + h ^= std::hash()(key.b) + 0x9e3779b97f4a7c15ull + (h << 6) + (h >> 2); + h ^= std::hash()(key.c) + 0x9e3779b97f4a7c15ull + (h << 6) + (h >> 2); + return h; + } + }; + + struct EdgeUseCount { + size_t forward = 0; + size_t reverse = 0; + }; + + struct MeshDiagnostics { + size_t vertices = 0; + size_t triangles = 0; + size_t unique_edges = 0; + size_t invalid_indices = 0; + size_t nonfinite_vertices = 0; + size_t degenerate_triangles = 0; + size_t zero_area_triangles = 0; + size_t duplicate_faces = 0; + size_t boundary_edges = 0; + size_t nonmanifold_edges = 0; + size_t orientation_conflicts = 0; + long long euler_characteristic = 0; + bool has_bounds = false; + Eigen::Vector3d bounds_min = Eigen::Vector3d::Zero(); + Eigen::Vector3d bounds_max = Eigen::Vector3d::Zero(); + double min_edge = std::numeric_limits::infinity(); + double max_edge = 0.; + double min_area = std::numeric_limits::infinity(); + double max_area = 0.; + }; + + struct ShellDiagnostics { + size_t faces = 0; + size_t loops = 0; + size_t edges = 0; + size_t faces_with_inner_loops = 0; + size_t max_loops_per_face = 0; + size_t non_planar_faces = 0; + size_t non_polygonal_edges = 0; + size_t implicit_vertices = 0; + }; + + std::optional explicit_point(const taxonomy::edge::ptr& edge) { + if (edge->start.index() != 1) { + return std::nullopt; + } + return std::get(edge->start)->ccomponents(); + } + + bool face_supported(const taxonomy::face::ptr& face) { + if (face->basis && face->basis->kind() != taxonomy::PLANE) { + return false; + } + for (const auto& loop : face->children) { + if (!loop->is_polyhedron()) { + return false; + } + for (const auto& edge : loop->children) { + if (edge->start.index() != 1) { + return false; + } + } + } + return true; + } + + bool extrusion_face_supported(const taxonomy::face::ptr& face) { + return face && face->children.size() == 1 && face_supported(face); + } + + bool face_basis(const taxonomy::face::ptr& face, Eigen::Vector3d& origin, Eigen::Vector3d& x, Eigen::Vector3d& y) { + double best_score = -1.; + std::vector best_points; + for (const auto& loop : face->children) { + std::vector points; + points.reserve(loop->children.size()); + for (const auto& edge : loop->children) { + auto point = explicit_point(edge); + if (!point) { + return false; + } + if (!points.empty()) { + const auto d = points.back() - *point; + if (d.squaredNorm() <= 1.e-24) { + continue; + } + } + points.push_back(*point); + } + if (points.size() > 1) { + const auto d = points.front() - points.back(); + if (d.squaredNorm() <= 1.e-24) { + points.pop_back(); + } + } + if (points.size() < 3) { + continue; + } + Eigen::Vector3d normal = Eigen::Vector3d::Zero(); + for (size_t i = 0; i < points.size(); ++i) { + const auto& a = points[i]; + const auto& b = points[(i + 1) % points.size()]; + normal(0) += (a(1) - b(1)) * (a(2) + b(2)); + normal(1) += (a(2) - b(2)) * (a(0) + b(0)); + normal(2) += (a(0) - b(0)) * (a(1) + b(1)); + } + const auto score = normal.squaredNorm(); + if (score <= 1.e-24 || score <= best_score) { + continue; + } + best_score = score; + best_points = std::move(points); + } + if (best_points.size() < 3) { + return false; + } + Eigen::Vector3d normal = Eigen::Vector3d::Zero(); + for (size_t i = 0; i < best_points.size(); ++i) { + const auto& a = best_points[i]; + const auto& b = best_points[(i + 1) % best_points.size()]; + normal(0) += (a(1) - b(1)) * (a(2) + b(2)); + normal(1) += (a(2) - b(2)) * (a(0) + b(0)); + normal(2) += (a(0) - b(0)) * (a(1) + b(1)); + } + if (normal.norm() < 1.e-12) { + return false; + } + origin = best_points.front(); + x = best_points[1] - best_points.front(); + x -= normal.normalized() * x.dot(normal.normalized()); + if (x.norm() < 1.e-12) { + return false; + } + x.normalize(); + y = normal.normalized().cross(x).normalized(); + return true; + } + + double signed_area(const manifold::SimplePolygonIdx& polygon) { + double area = 0.; + for (size_t i = 0; i < polygon.size(); ++i) { + const auto& a = polygon[i].pos; + const auto& b = polygon[(i + 1) % polygon.size()].pos; + area += a[0] * b[1] - b[0] * a[1]; + } + return 0.5 * area; + } + + Eigen::Vector3d mesh_vertex(const Mesh& mesh, size_t index) { + return Eigen::Vector3d( + mesh.vertProperties[index * mesh.numProp + 0], + mesh.vertProperties[index * mesh.numProp + 1], + mesh.vertProperties[index * mesh.numProp + 2]); + } + + std::string format_number(double value) { + std::ostringstream ss; + ss << std::setprecision(6) << value; + return ss.str(); + } + + std::string format_vector(const Eigen::Vector3d& value) { + return "(" + format_number(value(0)) + ", " + format_number(value(1)) + ", " + format_number(value(2)) + ")"; + } + + MeshDiagnostics diagnose_mesh(const Mesh& mesh, double precision) { + MeshDiagnostics diagnostics; + diagnostics.vertices = mesh.NumVert(); + diagnostics.triangles = mesh.NumTri(); + std::unordered_map edge_use_count; + std::unordered_map face_use_count; + for (size_t i = 0; i < diagnostics.vertices; ++i) { + auto p = mesh_vertex(mesh, i); + if (!std::isfinite(p(0)) || !std::isfinite(p(1)) || !std::isfinite(p(2))) { + diagnostics.nonfinite_vertices++; + continue; + } + if (!diagnostics.has_bounds) { + diagnostics.has_bounds = true; + diagnostics.bounds_min = p; + diagnostics.bounds_max = p; + } else { + diagnostics.bounds_min = diagnostics.bounds_min.cwiseMin(p); + diagnostics.bounds_max = diagnostics.bounds_max.cwiseMax(p); + } + } + for (size_t i = 0; i < diagnostics.triangles; ++i) { + auto a = mesh.triVerts[i * 3 + 0]; + auto b = mesh.triVerts[i * 3 + 1]; + auto c = mesh.triVerts[i * 3 + 2]; + if (a >= diagnostics.vertices || b >= diagnostics.vertices || c >= diagnostics.vertices) { + diagnostics.invalid_indices++; + continue; + } + auto pa = mesh_vertex(mesh, (size_t)a); + auto pb = mesh_vertex(mesh, (size_t)b); + auto pc = mesh_vertex(mesh, (size_t)c); + const auto ab = (pb - pa).norm(); + const auto bc = (pc - pb).norm(); + const auto ca = (pa - pc).norm(); + if (std::isfinite(ab) && ab > 0.) { + diagnostics.min_edge = std::min(diagnostics.min_edge, ab); + diagnostics.max_edge = std::max(diagnostics.max_edge, ab); + } + if (std::isfinite(bc) && bc > 0.) { + diagnostics.min_edge = std::min(diagnostics.min_edge, bc); + diagnostics.max_edge = std::max(diagnostics.max_edge, bc); + } + if (std::isfinite(ca) && ca > 0.) { + diagnostics.min_edge = std::min(diagnostics.min_edge, ca); + diagnostics.max_edge = std::max(diagnostics.max_edge, ca); + } + if (a == b || b == c || c == a) { + diagnostics.degenerate_triangles++; + continue; + } + const auto area = 0.5 * ((pb - pa).cross(pc - pa)).norm(); + if (std::isfinite(area)) { + diagnostics.min_area = std::min(diagnostics.min_area, area); + diagnostics.max_area = std::max(diagnostics.max_area, area); + if (area <= precision * precision) { + diagnostics.zero_area_triangles++; + } + } + std::array face = { a, b, c }; + std::sort(face.begin(), face.end()); + const FaceKey face_key{ face[0], face[1], face[2] }; + auto face_it = face_use_count.find(face_key); + if (face_it == face_use_count.end()) { + face_use_count.insert({ face_key, 1 }); + } else { + face_it->second++; + diagnostics.duplicate_faces++; + } + std::array edges = { + EdgeKey{ std::min(a, b), std::max(a, b) }, + EdgeKey{ std::min(b, c), std::max(b, c) }, + EdgeKey{ std::min(c, a), std::max(c, a) } + }; + std::array forward = { + a < b, + b < c, + c < a + }; + for (const auto& edge : edges) { + if (edge.a == edge.b) { + continue; + } + } + for (size_t j = 0; j < edges.size(); ++j) { + const auto& edge = edges[j]; + auto& use = edge_use_count[edge]; + if (forward[j]) { + use.forward++; + } else { + use.reverse++; + } + } + } + for (const auto& entry : edge_use_count) { + const auto total = entry.second.forward + entry.second.reverse; + if (total == 1) { + diagnostics.boundary_edges++; + } else if (total > 2) { + diagnostics.nonmanifold_edges++; + } else if (entry.second.forward != 1 || entry.second.reverse != 1) { + diagnostics.orientation_conflicts++; + } + } + diagnostics.unique_edges = edge_use_count.size(); + diagnostics.euler_characteristic = (long long)diagnostics.vertices - (long long)diagnostics.unique_edges + (long long)diagnostics.triangles; + return diagnostics; + } + + ShellDiagnostics diagnose_shell(const taxonomy::shell::ptr& shell) { + ShellDiagnostics diagnostics; + diagnostics.faces = shell->children.size(); + for (const auto& face : shell->children) { + diagnostics.max_loops_per_face = std::max(diagnostics.max_loops_per_face, face->children.size()); + if (face->children.size() > 1) { + diagnostics.faces_with_inner_loops++; + } + diagnostics.loops += face->children.size(); + if (face->basis && face->basis->kind() != taxonomy::PLANE) { + diagnostics.non_planar_faces++; + } + for (const auto& loop : face->children) { + diagnostics.edges += loop->children.size(); + for (const auto& edge : loop->children) { + if (edge->basis && edge->basis->kind() != taxonomy::LINE) { + diagnostics.non_polygonal_edges++; + } + if (edge->start.index() != 1) { + diagnostics.implicit_vertices++; + } + } + } + } + return diagnostics; + } + + std::string mesh_diagnostics_string(const MeshDiagnostics& diagnostics) { + std::ostringstream ss; + ss << "verts=" << diagnostics.vertices + << " tris=" << diagnostics.triangles + << " unique_edges=" << diagnostics.unique_edges + << " euler=" << diagnostics.euler_characteristic + << " invalid_idx=" << diagnostics.invalid_indices + << " nonfinite_verts=" << diagnostics.nonfinite_vertices + << " degenerate_tris=" << diagnostics.degenerate_triangles + << " zero_area_tris=" << diagnostics.zero_area_triangles + << " duplicate_faces=" << diagnostics.duplicate_faces + << " boundary_edges=" << diagnostics.boundary_edges + << " nonmanifold_edges=" << diagnostics.nonmanifold_edges + << " orientation_conflicts=" << diagnostics.orientation_conflicts; + if (diagnostics.has_bounds) { + ss << " bbox_min=" << format_vector(diagnostics.bounds_min) + << " bbox_max=" << format_vector(diagnostics.bounds_max); + } + if (std::isfinite(diagnostics.min_edge)) { + ss << " min_edge=" << format_number(diagnostics.min_edge); + } + if (diagnostics.max_edge > 0.) { + ss << " max_edge=" << format_number(diagnostics.max_edge); + } + if (std::isfinite(diagnostics.min_area)) { + ss << " min_area=" << format_number(diagnostics.min_area); + } + if (diagnostics.max_area > 0.) { + ss << " max_area=" << format_number(diagnostics.max_area); + } + return ss.str(); + } + + std::string shell_diagnostics_string(const ShellDiagnostics& diagnostics) { + std::ostringstream ss; + ss << "faces=" << diagnostics.faces + << " loops=" << diagnostics.loops + << " edges=" << diagnostics.edges + << " faces_with_inner_loops=" << diagnostics.faces_with_inner_loops + << " max_loops_per_face=" << diagnostics.max_loops_per_face + << " non_planar_faces=" << diagnostics.non_planar_faces + << " non_polygonal_edges=" << diagnostics.non_polygonal_edges + << " implicit_vertices=" << diagnostics.implicit_vertices; + return ss.str(); + } + + std::string matrix_diagnostics_string(const taxonomy::matrix4::ptr& place) { + const auto& m = place->ccomponents(); + const auto linear = m.block<3, 3>(0, 0); + const auto c0 = linear.col(0); + const auto c1 = linear.col(1); + const auto c2 = linear.col(2); + std::ostringstream ss; + ss << "det=" << format_number(linear.determinant()) + << " scale=(" << format_number(c0.norm()) << ", " << format_number(c1.norm()) << ", " << format_number(c2.norm()) << ")" + << " dot=(" << format_number(c0.dot(c1)) << ", " << format_number(c0.dot(c2)) << ", " << format_number(c1.dot(c2)) << ")" + << " translation=" << format_vector(m.col(3).head<3>()); + return ss.str(); + } + + std::string solid_shell_failure_diagnosis(const Part& part, const MeshDiagnostics& before, const MeshDiagnostics& after, manifold::Manifold::Error before_status, manifold::Manifold::Error after_status) { + const bool before_problematic = + !part.solid || + before_status != manifold::Manifold::Error::NoError || + before.invalid_indices != 0 || + before.nonfinite_vertices != 0 || + before.degenerate_triangles != 0 || + before.zero_area_triangles != 0 || + before.boundary_edges != 0 || + before.nonmanifold_edges != 0; + const bool after_problematic = + after_status != manifold::Manifold::Error::NoError || + after.invalid_indices != 0 || + after.nonfinite_vertices != 0 || + after.degenerate_triangles != 0 || + after.zero_area_triangles != 0 || + after.boundary_edges != 0 || + after.nonmanifold_edges != 0; + if (before_problematic) { + return "shell is already problematic before transform"; + } + if (after_problematic) { + return "shell is valid before transform, failure is likely introduced by transform or precision collapse"; + } + return "shell looks clean before and after mesh inspection, issue may be in manifold validation details"; + } + + void log_solid_shell_transform_failure(const taxonomy::shell::ptr& shell, const Part& before_part, const Mesh& after_mesh, const taxonomy::matrix4::ptr& place, double precision, manifold::Manifold::Error before_status, manifold::Manifold::Error after_status) { + const auto shell_info = diagnose_shell(shell); + const auto before = diagnose_mesh(before_part.mesh, precision); + const auto after = diagnose_mesh(after_mesh, precision); + logger::warning( + "Manifold kernel: solid shell manifold validation failed; before_transform=" + + std::string(before_part.solid ? "solid" : "mesh-only") + + " (" + manifold_error_string(before_status) + "), after_transform=(" + manifold_error_string(after_status) + ")", + shell->instance); + logger::warning("Manifold kernel: solid shell diagnosis: " + solid_shell_failure_diagnosis(before_part, before, after, before_status, after_status), shell->instance); + logger::warning("Manifold kernel: solid shell input: " + shell_diagnostics_string(shell_info), shell->instance); + logger::warning("Manifold kernel: solid shell mesh before transform: " + mesh_diagnostics_string(before), shell->instance); + logger::warning("Manifold kernel: solid shell transform: " + matrix_diagnostics_string(place), shell->instance); + logger::warning("Manifold kernel: solid shell mesh after transform: " + mesh_diagnostics_string(after), shell->instance); + } + + double signed_area(const manifold::SimplePolygon& polygon) { + double area = 0.; + for (size_t i = 0; i < polygon.size(); ++i) { + const auto& a = polygon[i]; + const auto& b = polygon[(i + 1) % polygon.size()]; + area += a[0] * b[1] - b[0] * a[1]; + } + return 0.5 * area; + } + + double signed_area(const LoopPolygon& polygon) { + double area = 0.; + for (size_t i = 0; i < polygon.size(); ++i) { + const auto& a = polygon[i].uv; + const auto& b = polygon[(i + 1) % polygon.size()].uv; + area += a[0] * b[1] - b[0] * a[1]; + } + return 0.5 * area; + } + + bool append_simple_loop(const taxonomy::loop::ptr& loop, const Eigen::Vector3d& origin, const Eigen::Vector3d& x, const Eigen::Vector3d& y, double precision, LoopPolygon& polygon) { + polygon.clear(); + polygon.reserve(loop->children.size()); + for (const auto& edge : loop->children) { + if (edge->basis && edge->basis->kind() != taxonomy::LINE) { + return false; + } + auto point = explicit_point(edge); + if (!point) { + return false; + } + if (!polygon.empty()) { + const auto d = polygon.back().xyz - *point; + if (d.squaredNorm() <= precision * precision) { + continue; + } + } + auto v = *point - origin; + polygon.push_back({ *point, manifold::vec2(v.dot(x), v.dot(y)) }); + } + if (polygon.size() > 1) { + const auto d = polygon.front().xyz - polygon.back().xyz; + if (d.squaredNorm() <= precision * precision) { + polygon.pop_back(); + } + } + if (polygon.size() < 3) { + return false; + } + return true; + } + + void reverse_loop(LoopPolygon& polygon) { + std::reverse(polygon.begin(), polygon.end()); + } + + void append_loop(const LoopPolygon& loop_polygon, MeshBuilder& builder, manifold::PolygonsIdx& polygons) { + manifold::SimplePolygonIdx polygon; + polygon.reserve(loop_polygon.size()); + for (const auto& point : loop_polygon) { + manifold::PolyVert poly_vert; + poly_vert.pos = point.uv; + poly_vert.idx = (int)builder.add_vertex(point.xyz); + polygon.push_back(poly_vert); + } + polygons.push_back(std::move(polygon)); + } + + bool append_face(const taxonomy::face::ptr& face, MeshBuilder& builder, uint64_t face_id) { + if (!face_supported(face)) { + return false; + } + Eigen::Vector3d origin; + Eigen::Vector3d x; + Eigen::Vector3d y; + if (!face_basis(face, origin, x, y)) { + return false; + } + std::vector loops; + loops.reserve(face->children.size()); + size_t outer_index = 0; + double outer_area = 0.; + manifold::PolygonsIdx polygons; + polygons.reserve(face->children.size()); + for (const auto& loop : face->children) { + LoopPolygon loop_polygon; + if (!append_simple_loop(loop, origin, x, y, builder.precision, loop_polygon)) { + return false; + } + const auto area = signed_area(loop_polygon); + if (std::abs(area) > std::abs(outer_area)) { + outer_area = area; + outer_index = loops.size(); + } + loops.push_back(std::move(loop_polygon)); + } + if (loops.empty() || std::abs(outer_area) < 1.e-12) { + return false; + } + for (size_t i = 0; i < loops.size(); ++i) { + if (i != outer_index && signed_area(loops[i]) * outer_area > 0.) { + reverse_loop(loops[i]); + } + append_loop(loops[i], builder, polygons); + } + auto triangles = manifold::TriangulateIdx(polygons, builder.precision, true); + for (const auto& tri : triangles) { + builder.add_triangle((uint32_t)tri[0], (uint32_t)tri[1], (uint32_t)tri[2], face_id); + } + return !triangles.empty(); + } + + bool shell_to_mesh(const taxonomy::shell::ptr& shell, double precision, Mesh& mesh) { + MeshBuilder builder(precision); + uint64_t face_id = 0; + bool any = false; + for (const auto& face : shell->children) { + if (!append_face(face, builder, face_id++)) { + return false; + } + any = true; + } + if (!any) { + return false; + } + mesh = builder.build(); + return mesh.NumTri() > 0; + } + + std::optional part_from_mesh(const Mesh& mesh, bool require_manifold, manifold::Manifold::Error* status_ptr = nullptr) { + auto solid = std::optional{}; + manifold::Manifold candidate(mesh); + auto status = candidate.Status(); + if (status_ptr) { + *status_ptr = status; + } + if (status == manifold::Manifold::Error::NoError) { + solid = candidate; + } + if (!solid && require_manifold) { + return std::nullopt; + } + if (solid) { + return Part{ solid->GetMeshGL64(), solid }; + } + return Part{ mesh, std::nullopt }; + } + + std::optional part_from_shell(const taxonomy::shell::ptr& shell, double precision, manifold::Manifold::Error* status_ptr = nullptr) { + Mesh mesh; + if (!shell_to_mesh(shell, precision, mesh)) { + return std::nullopt; + } + return part_from_mesh(mesh, false, status_ptr); + } + + Mesh transform_mesh(const Mesh& mesh, const taxonomy::matrix4::ptr& place); + + std::optional part_from_extrusion(const taxonomy::extrusion::ptr& extrusion, double precision) { + if (extrusion->depth < precision) { + return std::nullopt; + } + auto face = taxonomy::dcast(extrusion->basis); + if (!extrusion_face_supported(face)) { + return std::nullopt; + } + Eigen::Vector3d origin; + Eigen::Vector3d x; + Eigen::Vector3d y; + if (!face_basis(face, origin, x, y)) { + return std::nullopt; + } + auto dir = extrusion->direction->ccomponents(); + if (dir.norm() < 1.e-12) { + return std::nullopt; + } + dir.normalize(); + auto normal = x.cross(y); + if (normal.norm() < 1.e-12) { + return std::nullopt; + } + normal.normalize(); + auto direction_sign = normal.dot(dir); + if (std::abs(direction_sign) < 1.e-9) { + return std::nullopt; + } + LoopPolygon polygon; + if (!append_simple_loop(face->children.front(), origin, x, y, precision, polygon)) { + return std::nullopt; + } + manifold::SimplePolygonIdx polygon_idx; + polygon_idx.reserve(polygon.size()); + for (size_t i = 0; i < polygon.size(); ++i) { + polygon_idx.push_back({ polygon[i].uv, (int)i }); + } + manifold::PolygonsIdx polygons = { polygon_idx }; + auto triangles = manifold::TriangulateIdx(polygons, precision, true); + if (triangles.empty()) { + return std::nullopt; + } + auto offset = dir * extrusion->depth; + MeshBuilder builder(precision); + std::vector bottom; + std::vector top; + bottom.reserve(polygon.size()); + top.reserve(polygon.size()); + for (const auto& point : polygon) { + bottom.push_back(builder.add_vertex(point.xyz)); + top.push_back(builder.add_vertex(point.xyz + offset)); + } + for (const auto& tri : triangles) { + auto a = (size_t)tri[0]; + auto b = (size_t)tri[1]; + auto c = (size_t)tri[2]; + if (direction_sign > 0.) { + builder.add_triangle(bottom[c], bottom[b], bottom[a], 0); + builder.add_triangle(top[a], top[b], top[c], 1); + } else { + builder.add_triangle(bottom[a], bottom[b], bottom[c], 0); + builder.add_triangle(top[c], top[b], top[a], 1); + } + } + for (size_t i = 0; i < polygon.size(); ++i) { + auto j = (i + 1) % polygon.size(); + if (direction_sign > 0.) { + builder.add_triangle(bottom[i], bottom[j], top[j], 2 + (uint64_t)i); + builder.add_triangle(bottom[i], top[j], top[i], 2 + (uint64_t)i); + } else { + builder.add_triangle(bottom[i], top[j], bottom[j], 2 + (uint64_t)i); + builder.add_triangle(bottom[i], top[i], top[j], 2 + (uint64_t)i); + } + } + return part_from_mesh(builder.build(), true); + } + + Mesh transform_mesh(const Mesh& mesh, const taxonomy::matrix4::ptr& place) { + const auto& m = place->ccomponents(); + Mesh result = mesh; + const bool flip = m.block<3, 3>(0, 0).determinant() < 0.; + for (size_t i = 0; i < mesh.NumVert(); ++i) { + Eigen::Vector4d v( + mesh.vertProperties[i * mesh.numProp + 0], + mesh.vertProperties[i * mesh.numProp + 1], + mesh.vertProperties[i * mesh.numProp + 2], + 1.); + auto v2 = m * v; + result.vertProperties[i * result.numProp + 0] = v2(0); + result.vertProperties[i * result.numProp + 1] = v2(1); + result.vertProperties[i * result.numProp + 2] = v2(2); + } + if (flip) { + for (size_t i = 0; i < mesh.NumTri(); ++i) { + std::swap(result.triVerts[i * 3 + 1], result.triVerts[i * 3 + 2]); + } + } + result.runTransform.clear(); + return result; + } + + taxonomy::style::ptr fallback_style(const taxonomy::geom_item::ptr& item, const IfcGeom::ConversionResults& results) { + if (item->surface_style) { + return item->surface_style; + } + for (const auto& result : results) { + if (result.hasStyle()) { + return result.StylePtr(); + } + } + return nullptr; + } + + std::optional result_to_manifold(const IfcGeom::ConversionResult& result) { + auto moved = std::unique_ptr(result.apply_transform()); + auto* shape = dynamic_cast(moved.get()); + if (!shape) { + return std::nullopt; + } + return shape->as_manifold(); + } + + std::optional results_to_operand(const IfcGeom::ConversionResults& results) { + std::vector operands; + for (const auto& result : results) { + auto operand = result_to_manifold(result); + if (operand) { + operands.push_back(*operand); + } + } + if (operands.empty()) { + return std::nullopt; + } + if (operands.size() == 1) { + return operands.front(); + } + return manifold::Manifold::BatchBoolean(operands, manifold::OpType::Add); + } + + std::optional boolean_result_from_operands(const std::vector& operands, taxonomy::boolean_result::operation_t operation) { + if (operands.empty()) { + return std::nullopt; + } + if (operands.size() == 1) { + return operands.front(); + } + switch (operation) { + case taxonomy::boolean_result::UNION: + return manifold::Manifold::BatchBoolean(operands, manifold::OpType::Add); + case taxonomy::boolean_result::INTERSECTION: + return manifold::Manifold::BatchBoolean(operands, manifold::OpType::Intersect); + case taxonomy::boolean_result::SUBTRACTION: + return manifold::Manifold::BatchBoolean(operands, manifold::OpType::Subtract); + } + return std::nullopt; + } +} + +bool ManifoldKernel::convert_impl(const taxonomy::extrusion::ptr extrusion, IfcGeom::ConversionResults& results) { + auto part = part_from_extrusion(extrusion, settings_.get().get()); + if (!part) { + logger::warning("Manifold kernel: failed to convert extrusion, only simple extrusions with a single polygonal outer bound are supported", extrusion->instance); + return false; + } + results.emplace_back(IfcGeom::ConversionResult( + extrusion->instance.id(), + extrusion->matrix, + new ifcopenshell::geometry::ManifoldShape(std::move(*part)), + extrusion->surface_style)); + return true; +} + +bool ManifoldKernel::convert_impl(const taxonomy::shell::ptr shell, IfcGeom::ConversionResults& results) { + manifold::Manifold::Error status = manifold::Manifold::Error::NoError; + auto part = part_from_shell(shell, settings_.get().get(), &status); + if (!part) { + logger::warning("Manifold kernel: failed to convert shell, requires planar polygonal faces with explicit vertices", shell->instance); + return false; + } + if (!part->solid) { + logger::notice("Manifold kernel: shell converted as mesh only (" + manifold_error_string(status) + ")", shell->instance); + } + results.emplace_back(IfcGeom::ConversionResult( + shell->instance.id(), + shell->matrix, + new ifcopenshell::geometry::ManifoldShape(std::move(*part)), + shell->surface_style)); + return true; +} + +bool ManifoldKernel::convert_impl(const taxonomy::solid::ptr solid, IfcGeom::ConversionResults& results) { + std::vector shells; + for (const auto& shell : solid->children) { + const auto precision = settings_.get().get(); + manifold::Manifold::Error before_status = manifold::Manifold::Error::NoError; + auto part = part_from_shell(shell, precision, &before_status); + if (!part) { + logger::warning("Manifold kernel: failed to convert solid shell, requires planar polygonal faces with explicit vertices", shell->instance); + return false; + } + auto place = shell->matrix ? shell->matrix : taxonomy::make(); + auto transformed_mesh = transform_mesh(part->mesh, place); + manifold::Manifold::Error after_status = manifold::Manifold::Error::NoError; + auto transformed = part_from_mesh(transformed_mesh, true, &after_status); + if (!transformed || !transformed->solid) { + log_solid_shell_transform_failure(shell, *part, transformed_mesh, place, precision, before_status, after_status); + return false; + } + shells.push_back(*transformed->solid); + } + if (shells.empty()) { + return false; + } + auto result = shells.front(); + for (size_t i = 1; i < shells.size(); ++i) { + result -= shells[i]; + } + results.emplace_back(IfcGeom::ConversionResult( + solid->instance.id(), + solid->matrix, + new ifcopenshell::geometry::ManifoldShape(Part{ result.GetMeshGL64(), result }), + solid->surface_style)); + return true; +} + +bool ManifoldKernel::convert_impl(const taxonomy::boolean_result::ptr br, IfcGeom::ConversionResults& results) { + std::vector operands; + taxonomy::style::ptr style; + for (const auto& child : br->children) { + IfcGeom::ConversionResults converted; + if (!AbstractKernel::convert(child, converted)) { + logger::warning("Manifold kernel: failed to convert boolean operand", child->instance); + return false; + } + auto operand = results_to_operand(converted); + if (!operand) { + logger::warning("Manifold kernel: boolean operand is not a valid manifold solid", child->instance); + return false; + } + operands.push_back(*operand); + if (!style) { + style = fallback_style(child, converted); + } + } + auto result = boolean_result_from_operands(operands, br->operation); + if (!result || result->IsEmpty()) { + logger::warning("Manifold kernel: boolean operation produced no result", br->instance); + return false; + } + results.emplace_back(IfcGeom::ConversionResult( + br->instance.id(), + br->matrix, + new ifcopenshell::geometry::ManifoldShape(Part{ result->GetMeshGL64(), *result }), + br->surface_style ? br->surface_style : style)); + return true; +} + +bool ManifoldKernel::convert_openings(const express::Base&, const std::vector>& openings, const IfcGeom::ConversionResults& entity_shapes, const taxonomy::matrix4& entity_trsf, IfcGeom::ConversionResults& cut_shapes) { + std::vector opening_operands; + for (const auto& opening : openings) { + IfcGeom::ConversionResults converted; + if (!AbstractKernel::convert(opening.first, converted)) { + logger::warning("Manifold kernel: failed to convert opening operand", opening.first->instance); + return false; + } + const auto relative = taxonomy::make(entity_trsf.ccomponents().inverse() * opening.second.ccomponents()); + for (const auto& result : converted) { + auto moved = std::unique_ptr(result.Shape()->moved(taxonomy::make(relative->ccomponents() * result.Placement()->ccomponents()))); + auto* shape = dynamic_cast(moved.get()); + if (!shape) { + logger::warning("Manifold kernel: opening result is not a manifold shape"); + return false; + } + auto operand = shape->as_manifold(); + if (!operand) { + logger::warning("Manifold kernel: opening result is not a valid manifold solid", opening.first->instance); + return false; + } + opening_operands.push_back(*operand); + } + } + if (opening_operands.empty()) { + return false; + } + auto opening_union = manifold::Manifold::BatchBoolean(opening_operands, manifold::OpType::Add); + for (const auto& entity_shape : entity_shapes) { + auto operand = result_to_manifold(entity_shape); + if (!operand) { + logger::warning("Manifold kernel: host shape is not a valid manifold solid"); + return false; + } + auto result = *operand - opening_union; + cut_shapes.emplace_back(IfcGeom::ConversionResult( + entity_shape.ItemId(), + new ifcopenshell::geometry::ManifoldShape(Part{ result.GetMeshGL64(), result }), + entity_shape.StylePtr())); + } + return !cut_shapes.empty(); +} diff --git a/src/ifcgeom/kernels/manifold/ManifoldKernel.h b/src/ifcgeom/kernels/manifold/ManifoldKernel.h new file mode 100644 index 0000000000..70c4dbfeae --- /dev/null +++ b/src/ifcgeom/kernels/manifold/ManifoldKernel.h @@ -0,0 +1,38 @@ +#ifndef MANIFOLD_KERNEL_H +#define MANIFOLD_KERNEL_H + +#include + +#include "../../../ifcgeom/AbstractKernel.h" +#include "../../../ifcgeom/kernels/ifc_geomlibrary_api.h" +#include "../../../ifcgeom/kernels/manifold/ManifoldConversionResult.h" + +namespace ifcopenshell { +namespace geometry { +namespace kernels { + +class IFC_GEOMLIBRARY_API ManifoldKernel : public AbstractKernel { +public: + ManifoldKernel(const Settings& settings) + : AbstractKernel("manifold", settings) {} + + virtual AbstractKernel* clone() const { + return new ManifoldKernel(settings()); + } + + virtual bool supports_boolean_operations() const { return true; } + + virtual bool convert_impl(const taxonomy::extrusion::ptr, IfcGeom::ConversionResults&); + virtual bool convert_impl(const taxonomy::shell::ptr, IfcGeom::ConversionResults&); + virtual bool convert_impl(const taxonomy::solid::ptr, IfcGeom::ConversionResults&); + virtual bool convert_impl(const taxonomy::boolean_result::ptr, IfcGeom::ConversionResults&); + + virtual bool convert_openings(const express::Base& entity, const std::vector>& openings, + const IfcGeom::ConversionResults& entity_shapes, const ifcopenshell::geometry::taxonomy::matrix4& entity_trsf, IfcGeom::ConversionResults& cut_shapes); +}; + +} +} +} + +#endif diff --git a/win/build-deps.cmd b/win/build-deps.cmd index f490f4a20f..2ad44b8444 100644 --- a/win/build-deps.cmd +++ b/win/build-deps.cmd @@ -200,7 +200,7 @@ IF "%IFCOS_INSTALL_PYTHON%"=="TRUE" ( echo PYTHONHOME=%PYTHONHOME%>>"%~dp0\%BUILD_DEPS_CACHE_PATH%" ) -goto :SWIG +goto :manifold :nuget set DEPENDENCY_NAME=nuget @@ -732,6 +732,39 @@ IF NOT %ERRORLEVEL%==0 GOTO :Error call :InstallCMakeProject "%DEPENDENCY_DIR%\%BUILD_DIR%" %BUILD_CFG% IF NOT %ERRORLEVEL%==0 GOTO :Error +:manifold +set DEPENDENCY_NAME=manifold +set MANIFOLD_VERSION=3.2.1 +set DEPENDENCY_DIR=%DEPS_DIR%\manifold-%MANIFOLD_VERSION% +set DEPENDENCY_INSTALL_DIR=%INSTALL_DIR%\manifold-%MANIFOLD_VERSION% +echo MANIFOLD_ROOT=%DEPENDENCY_INSTALL_DIR%>>"%~dp0\%BUILD_DEPS_CACHE_PATH%" + +IF EXIST "%DEPENDENCY_INSTALL_DIR%" ( + echo Found existing "%DEPENDENCY_INSTALL_DIR%", skipping + goto :Eigen +) + +call :GitCloneAndCheckoutRevision https://github.com/elalish/manifold.git "%DEPENDENCY_DIR%" v%MANIFOLD_VERSION% +IF NOT %ERRORLEVEL%==0 GOTO :Error +cd "%DEPENDENCY_DIR%" +git reset --hard + +call :RunCMake -DCMAKE_INSTALL_PREFIX="%DEPENDENCY_INSTALL_DIR%" ^ + -DBUILD_SHARED_LIBS=OFF ^ + -DMANIFOLD_PAR=OFF ^ + -DMANIFOLD_CROSS_SECTION=OFF ^ + -DMANIFOLD_PYBIND=OFF ^ + -DMANIFOLD_JSBIND=OFF ^ + -DMANIFOLD_CBIND=OFF ^ + -DMANIFOLD_TEST=OFF ^ + -DMANIFOLD_EXPORT=OFF ^ + -DMANIFOLD_DOWNLOADS=OFF +IF NOT %ERRORLEVEL%==0 GOTO :Error +call :BuildSolution "%DEPENDENCY_DIR%\%BUILD_DIR%\manifold.sln" %BUILD_CFG% +IF NOT %ERRORLEVEL%==0 GOTO :Error +call :InstallCMakeProject "%DEPENDENCY_DIR%\%BUILD_DIR%" %BUILD_CFG% +IF NOT %ERRORLEVEL%==0 GOTO :Error + :: :tbb :: set DEPENDENCY_NAME=tbb :: set DEPENDENCY_DIR=%DEPS_DIR%\tbb