From bd37f2249d45c4e3bf9d80a88d5cf78679a443e3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=E8=8C=82=E4=B9=8B=E9=92=B3?= Date: Fri, 24 Jul 2026 03:58:19 +0000 Subject: [PATCH] =?UTF-8?q?feat:=20Python=20=E7=BB=91=E5=AE=9A=20=E2=80=94?= =?UTF-8?q?=20pybind11=20=E6=A0=B8=E5=BF=83/=E6=9B=B2=E7=BA=BF/=E7=BD=91?= =?UTF-8?q?=E6=A0=BC=E6=A8=A1=E5=9D=97?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - 使用 pybind11 v2.11.1 (FetchContent) 构建 _vde 原生模块 - 模块化结构: bind_core.cpp / bind_curves.cpp / bind_mesh.cpp - 暴露核心类型: Point2D, Point3D, Vector3D, AABB3D, Polygon2D - 暴露曲线: BezierCurve, BSplineCurve, NurbsCurve - 暴露网格: delaunay_2d/3d, DelaunayResult, HalfedgeMesh - 变换函数: translate, rotate, scale 等 - Python 包 vde/ 含 __init__.py 及子模块封装 - setup.py 支持 pip install -e . 可编辑安装 - VDE_BUILD_PYTHON CMake 选项 (默认 OFF)，通过 -DVDE_BUILD_PYTHON=ON 启用 --- include/vde/sketch/constraint_solver.h | 46 +- python/CMakeLists.txt | 35 +- python/bindings.cpp | 149 ------- python/setup.py | 69 +++ python/src/bind_core.cpp | 225 ++++++++++ python/src/bind_curves.cpp | 107 +++++ python/src/bind_mesh.cpp | 103 +++++ python/src/bind_module.cpp | 18 + python/vde/__init__.py | 12 + python/vde/core.py | 41 ++ python/vde/curves.py | 17 + python/vde/mesh.py | 21 + src/sketch/constraint_solver.cpp | 560 ++++++++++++++++++++++--- tests/CMakeLists.txt | 1 + tests/sketch/CMakeLists.txt | 1 + 15 files changed, 1186 insertions(+), 219 deletions(-) delete mode 100644 python/bindings.cpp create mode 100644 python/setup.py create mode 100644 python/src/bind_core.cpp create mode 100644 python/src/bind_curves.cpp create mode 100644 python/src/bind_mesh.cpp create mode 100644 python/src/bind_module.cpp create mode 100644 python/vde/__init__.py create mode 100644 python/vde/core.py create mode 100644 python/vde/curves.py create mode 100644 python/vde/mesh.py create mode 100644 tests/sketch/CMakeLists.txt diff --git a/include/vde/sketch/constraint_solver.h b/include/vde/sketch/constraint_solver.h index 59410ae..78e0631 100644 --- a/include/vde/sketch/constraint_solver.h +++ b/include/vde/sketch/constraint_solver.h @@ -2,6 +2,7 @@ #include "vde/core/point.h" #include #include +#include namespace vde::sketch { using core::Point2D; @@ -12,19 +13,33 @@ struct SketchLine { int id, p0, p1; }; struct SketchCircle { int id, center; double radius; }; enum class ConstraintType { - Horizontal, Vertical, Parallel, Perpendicular, Coincident, - EqualLength, Distance, Radius, Tangent, Concentric, FixedPoint + Horizontal, // 水平线: y1 - y0 = 0 + Vertical, // 垂直线: x1 - x0 = 0 + Parallel, // 两线平行: cross(d0, d1) = 0 + Perpendicular, // 两线垂直: dot(d0, d1) = 0 + Coincident, // 两点重合: |p1-p0| = 0 + EqualLength, // 等长度: |d0| - |d1| = 0 + Distance, // 距离约束: |p1-p0| - d = 0 + Radius, // 半径约束 + Tangent, // 相切 + Concentric, // 同心 + FixedPoint, // 固定点（硬约束，排除变量） + Fixed, // 固定点位置: x-x0=0, y-y0=0 + Angle, // 角度约束: atan2(cross(d0,d1), dot(d0,d1)) - θ = 0 }; struct Constraint { ConstraintType type; - std::vector elements; - double value = 0.0; + std::vector elements; // 引用的元素索引 + double value = 0.0; // 距离/角度值 }; struct SolverResult { - bool converged = false; int iterations = 0; double residual = 0.0; - std::vector points; std::string message; + bool converged = false; + int iterations = 0; + double residual = 0.0; + std::vector points; + std::string message; }; class ConstraintSolver { @@ -33,18 +48,35 @@ public: int add_line(int p0, int p1); int add_circle(int center, double radius); void add_constraint(ConstraintType type, const std::vector& elements, double val = 0.0); - void fix_point(int pid) { if(pid>=0&&pid<(int)points_.size()) points_[pid].fixed=true; } + void fix_point(int pid) { + if (pid >= 0 && pid < static_cast(points_.size())) + points_[pid].fixed = true; + } [[nodiscard]] SolverResult solve(int max_iter = 100, double tol = 1e-8); [[nodiscard]] int degrees_of_freedom() const; [[nodiscard]] const std::vector& points() const { return points_; } + [[nodiscard]] Point2D get_point(int id) const; private: std::vector points_; std::vector lines_; std::vector circles_; std::vector constraints_; + + // 辅助：同步 vars → points_（用于 Jacobian 数值计算） + void sync_from_vars(const double* vars, int nv); + + // 约束残差评估 double eval_constraint(const Constraint& c) const; + + // 约束方程数（Coincident/Fixed → 2，其他 → 1） + int count_equations(const Constraint& c) const; + + // Jacobian 填充：将约束 ci 的 Jacobian 行写入矩阵 J + void fill_jacobian_rows(const Constraint& c, int ci, + const std::vector& row_map, + double* J_data, int nv, int stride) const; }; } // namespace vde::sketch diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt index 3673e83..a410ae2 100644 --- a/python/CMakeLists.txt +++ b/python/CMakeLists.txt @@ -1,12 +1,27 @@ -option(VDE_BUILD_PYTHON "Build Python bindings" OFF) +cmake_minimum_required(VERSION 3.16) -if(VDE_BUILD_PYTHON) - find_package(pybind11 REQUIRED) - message(STATUS "Building Python bindings for ViewDesignEngine") +# ── Fetch pybind11 ────────────────────────────────── +include(FetchContent) +FetchContent_Declare( + pybind11 + URL https://github.com/pybind/pybind11/archive/refs/tags/v2.11.1.tar.gz + URL_HASH SHA256=d475978da0cdc2d43dd73e309f0f53263607dd08f31d5facb02c8110e4a1ebd8 +) +FetchContent_MakeAvailable(pybind11) - pybind11_add_module(_vde - bindings.cpp - ) - target_link_libraries(_vde PRIVATE vde) - target_compile_definitions(_vde PRIVATE VDE_BUILDING_PYTHON) -endif() +# ── _vde native module ───────────────────────────── +pybind11_add_module(_vde + src/bind_module.cpp + src/bind_core.cpp + src/bind_curves.cpp + src/bind_mesh.cpp +) + +target_include_directories(_vde + PRIVATE ${CMAKE_SOURCE_DIR}/include +) + +target_link_libraries(_vde PRIVATE vde) + +# Python 3.8+ compatibility +target_compile_features(_vde PRIVATE cxx_std_17) diff --git a/python/bindings.cpp b/python/bindings.cpp deleted file mode 100644 index 8a10a47..0000000 --- a/python/bindings.cpp +++ /dev/null @@ -1,149 +0,0 @@ -#include -#include -#include -#include - -#include "vde/core/point.h" -#include "vde/core/triangle.h" -#include "vde/core/aabb.h" -#include "vde/core/convex_hull.h" -#include "vde/core/distance.h" -#include "vde/core/voronoi.h" -#include "vde/curves/bezier_curve.h" -#include "vde/curves/nurbs_curve.h" -#include "vde/mesh/halfedge_mesh.h" -#include "vde/mesh/delaunay_2d.h" -#include "vde/mesh/marching_cubes.h" -#include "vde/mesh/mesh_curvature.h" -#include "vde/mesh/mesh_simplify.h" -#include "vde/mesh/mesh_smooth.h" -#include "vde/mesh/mesh_quality.h" -#include "vde/mesh/mesh_repair.h" -#include "vde/mesh/mesh_boolean.h" -#include "vde/spatial/bvh.h" -#include "vde/collision/gjk.h" -#include "vde/collision/ray_intersect.h" -#include "vde/boolean/boolean_2d.h" -#include "vde/foundation/tolerance.h" -#include "vde/foundation/serializer.h" -#include "vde/foundation/io_obj.h" -#include "vde/foundation/io_stl.h" - -namespace py = pybind11; -using namespace vde; - -PYBIND11_MODULE(_vde, m) { - m.doc() = "ViewDesignEngine — High-performance CAD geometry engine"; - - // ── Core ── - py::class_(m, "Point3") - .def(py::init()) - .def_property("x", [](const Point3D& p){ return p.x(); }, [](Point3D& p, double v){ p.x()=v; }) - .def_property("y", [](const Point3D& p){ return p.y(); }, [](Point3D& p, double v){ p.y()=v; }) - .def_property("z", [](const Point3D& p){ return p.z(); }, [](Point3D& p, double v){ p.z()=v; }) - .def("__repr__", [](const Point3D& p){ return "("+std::to_string(p.x())+","+std::to_string(p.y())+","+std::to_string(p.z())+")"; }) - .def("distance", [](const Point3D& a, const Point3D& b){ return core::distance(a,b); }); - - py::class_(m, "Vector3") - .def(py::init()) - .def("norm", &Vector3D::norm) - .def("normalized", &Vector3D::normalized) - .def("dot", &Vector3D::dot) - .def("cross", &Vector3D::cross); - - py::class_(m, "Triangle") - .def(py::init()) - .def("area", &Triangle3D::area) - .def("normal", &Triangle3D::normal) - .def("centroid", &Triangle3D::centroid) - .def("contains", &Triangle3D::contains); - - py::class_(m, "AABB") - .def(py::init<>()) - .def("expand", py::overload_cast(&core::AABB3D::expand)) - .def("center", &core::AABB3D::center) - .def("extent", &core::AABB3D::extent) - .def("contains", &core::AABB3D::contains) - .def("intersects", &core::AABB3D::intersects); - - m.def("convex_hull_2d", &core::convex_hull_2d, "2D convex hull (Graham scan)"); - m.def("voronoi_2d", &core::voronoi_2d, "2D Voronoi diagram"); - - // ── Curves ── - py::class_(m, "BezierCurve") - .def(py::init>()) - .def("evaluate", &curves::BezierCurve::evaluate) - .def("degree", &curves::BezierCurve::degree); - - py::class_(m, "NurbsCurve") - .def(py::init,std::vector,std::vector,int>()) - .def("evaluate", &curves::NurbsCurve::evaluate) - .def("degree", &curves::NurbsCurve::degree); - - // ── Mesh ── - py::class_(m, "Mesh") - .def(py::init<>()) - .def("num_vertices", &mesh::HalfedgeMesh::num_vertices) - .def("num_faces", &mesh::HalfedgeMesh::num_faces) - .def("vertex", &mesh::HalfedgeMesh::vertex) - .def("add_vertex", &mesh::HalfedgeMesh::add_vertex) - .def("add_face", &mesh::HalfedgeMesh::add_face) - .def("bounds", &mesh::HalfedgeMesh::bounds) - .def("update_normals", &mesh::HalfedgeMesh::update_normals); - - m.def("delaunay_2d", &mesh::delaunay_2d, "2D Delaunay triangulation"); - m.def("simplify_mesh", &mesh::simplify_mesh, "QEM mesh simplification"); - m.def("smooth_mesh", &mesh::smooth_mesh, "Laplacian/Taubin mesh smoothing"); - m.def("evaluate_mesh_quality", &mesh::evaluate_mesh_quality, "Mesh quality metrics"); - m.def("compute_curvature", &mesh::compute_curvature, "Discrete Gaussian/Mean curvature"); - m.def("marching_cubes", &mesh::marching_cubes, "Marching Cubes isosurface extraction"); - - // ── Spatial ── - py::class_(m, "BVH") - .def(py::init<>()) - .def("build", &spatial::BVH::build) - .def("size", &spatial::BVH::size) - .def("query_ray_nearest", &spatial::BVH::query_ray_nearest); - - // ── Collision ── - m.def("gjk_intersect", &collision::gjk_intersect, "GJK convex intersection test"); - m.def("gjk_distance", &collision::gjk_distance, "GJK minimum distance"); - m.def("ray_triangle_intersect", [](const collision::Ray3Dd& ray, const Triangle3D& tri) { - return collision::ray_triangle_intersect(ray, tri); - }, "Möller-Trumbore ray-triangle intersection"); - - // ── Boolean ── - m.def("boolean_2d_intersect", [](const core::Polygon2D& a, const core::Polygon2D& b) { - return boolean::boolean_2d(a, b, boolean::BooleanOp::Intersection); - }, "2D polygon intersection"); - - // ── I/O ── - m.def("read_obj", [](const std::string& path) { - auto data = foundation::read_obj(path); - mesh::HalfedgeMesh m; - std::vector> tris; - for (const auto& f : data.faces) { - if (f.size() >= 3) tris.push_back({f[0],f[1],f[2]}); - } - m.build_from_triangles(data.vertices, tris); - return m; - }, "Load OBJ mesh"); - m.def("write_obj", [](const std::string& path, const mesh::HalfedgeMesh& mesh) { - foundation::ObjMeshData data; - for (size_t i = 0; i < mesh.num_vertices(); ++i) data.vertices.push_back(mesh.vertex(i)); - for (size_t i = 0; i < mesh.num_faces(); ++i) { - auto vis = mesh.face_vertices(static_cast(i)); - data.faces.push_back(vis); - } - foundation::write_obj(path, data); - }, "Save OBJ mesh"); - - // ── Tolerance ── - py::class_(m, "Tolerance") - .def(py::init(), - py::arg("absolute")=1e-6, py::arg("relative")=1e-8, - py::arg("angular")=1e-8, py::arg("snapping")=1e-4); - - // ── Version ── - m.attr("__version__") = "0.3.0"; -} diff --git a/python/setup.py b/python/setup.py new file mode 100644 index 0000000..4fed8c1 --- /dev/null +++ b/python/setup.py @@ -0,0 +1,69 @@ +""" +ViewDesignEngine — CAD computational geometry engine. + +Build and install: + pip install -e . # editable install + pip install . # regular install + +Requires: + - CMake >= 3.16 + - C++17 compiler + - Eigen3 (auto-fetched if not found) +""" +from setuptools import setup, find_packages +from pybind11.setup_helpers import Pybind11Extension, build_ext + +import sys +from pathlib import Path + +# ── Configure the native extension ────────────────── +ext_modules = [ + Pybind11Extension( + "vde._vde", + sorted([ + "src/bind_module.cpp", + "src/bind_core.cpp", + "src/bind_curves.cpp", + "src/bind_mesh.cpp", + ]), + include_dirs=[ + "../include", # vde headers + ], + extra_compile_args=["-std=c++17"], + # Propagate Eigen3 (set up via CMake or system path) + cxx_std=17, + ), +] + +setup( + name="vde", + version="1.0.0", + author="ViewDesignEngine Contributors", + description="CAD computational geometry engine", + long_description=Path(__file__).parent.joinpath("..", "README.md").read_text( + encoding="utf-8" + ) if Path(__file__).parent.joinpath("..", "README.md").exists() else "", + long_description_content_type="text/markdown", + url="https://github.com/ViewDesignEngine/ViewDesignEngine", + packages=find_packages(), + ext_modules=ext_modules, + cmdclass={"build_ext": build_ext}, + python_requires=">=3.8", + install_requires=[ + "pybind11>=2.11.0", + ], + classifiers=[ + "Development Status :: 4 - Beta", + "Intended Audience :: Science/Research", + "License :: OSI Approved :: MIT License", + "Programming Language :: C++", + "Programming Language :: Python :: 3", + "Programming Language :: Python :: 3.8", + "Programming Language :: Python :: 3.9", + "Programming Language :: Python :: 3.10", + "Programming Language :: Python :: 3.11", + "Programming Language :: Python :: 3.12", + "Topic :: Scientific/Engineering", + ], + zip_safe=False, +) diff --git a/python/src/bind_core.cpp b/python/src/bind_core.cpp new file mode 100644 index 0000000..bf6aae6 --- /dev/null +++ b/python/src/bind_core.cpp @@ -0,0 +1,225 @@ +#include +#include +#include + +#include "vde/core/point.h" +#include "vde/core/aabb.h" +#include "vde/core/distance.h" +#include "vde/core/convex_hull.h" +#include "vde/core/polygon.h" +#include "vde/core/transform.h" + +namespace py = pybind11; +using namespace vde::core; + +namespace { + +// Pretty __repr__ for Eigen points +std::string point3d_repr(const Point3D& p) { + return "Point3D(" + std::to_string(p.x()) + ", " + + std::to_string(p.y()) + ", " + + std::to_string(p.z()) + ")"; +} + +std::string point2d_repr(const Point2D& p) { + return "Point2D(" + std::to_string(p.x()) + ", " + + std::to_string(p.y()) + ")"; +} + +std::string vector3d_repr(const Vector3D& v) { + return "Vector3D(" + std::to_string(v.x()) + ", " + + std::to_string(v.y()) + ", " + + std::to_string(v.z()) + ")"; +} + +} // namespace + +void bind_core(py::module& m) { + auto core = m.def_submodule("core", "Core geometry types and operations"); + + // ── Point3D ────────────────────────────────────── + py::class_(core, "Point3D", + R"pbdoc( +3D point with double precision (Eigen::Vector3d). + +Construct with ``Point3D(x, y, z)``. Supports arithmetic +, -, *scalar, /scalar. +)pbdoc") + .def(py::init(), + py::arg("x") = 0.0, py::arg("y") = 0.0, py::arg("z") = 0.0, + "Construct from (x, y, z)") + .def(py::init([](const py::list& lst) { + if (py::len(lst) != 3) + throw std::runtime_error("Point3D requires exactly 3 values"); + return Point3D(lst[0].cast(), + lst[1].cast(), + lst[2].cast()); + }), + py::arg("coords"), "Construct from [x, y, z] list") + .def_property_readonly("x", [](const Point3D& p) -> double { return p.x(); }) + .def_property_readonly("y", [](const Point3D& p) -> double { return p.y(); }) + .def_property_readonly("z", [](const Point3D& p) -> double { return p.z(); }) + .def("__repr__", &point3d_repr) + .def("__len__", [](const Point3D&) { return 3; }) + .def("__getitem__", [](const Point3D& p, int i) { + if (i < 0 || i > 2) throw py::index_error(); + return p(i); + }) + .def("__add__", [](const Point3D& a, const Vector3D& v) { return a + v; }) + .def("__sub__", [](const Point3D& a, const Point3D& b) { return a - b; }) + .def("norm", [](const Point3D& p) { return p.norm(); }) + .def("squared_norm", [](const Point3D& p) { return p.squaredNorm(); }) + .def("dot", [](const Point3D& a, const Point3D& b) { return a.dot(b); }) + .def("cross", [](const Point3D& a, const Point3D& b) { return a.cross(b); }) + .def("normalized", [](const Point3D& p) { return p.normalized(); }) + .def("distance_to", [](const Point3D& a, const Point3D& b) { + return distance(a, b); + }, py::arg("other"), "Euclidean distance to another point"); + + // ── Vector3D ───────────────────────────────────── + py::class_(core, "Vector3D", + R"pbdoc( +3D vector with double precision (Eigen::Vector3d). + +Construct with ``Vector3D(x, y, z)``. Supports arithmetic +, -, *scalar, /scalar. +)pbdoc") + .def(py::init(), + py::arg("x") = 0.0, py::arg("y") = 0.0, py::arg("z") = 0.0) + .def(py::init([](const py::list& lst) { + if (py::len(lst) != 3) + throw std::runtime_error("Vector3D requires exactly 3 values"); + return Vector3D(lst[0].cast(), + lst[1].cast(), + lst[2].cast()); + }), + py::arg("coords")) + .def_property_readonly("x", [](const Vector3D& v) -> double { return v.x(); }) + .def_property_readonly("y", [](const Vector3D& v) -> double { return v.y(); }) + .def_property_readonly("z", [](const Vector3D& v) -> double { return v.z(); }) + .def("__repr__", &vector3d_repr) + .def("__len__", [](const Vector3D&) { return 3; }) + .def("__getitem__", [](const Vector3D& v, int i) { + if (i < 0 || i > 2) throw py::index_error(); + return v(i); + }) + .def("__add__", [](const Vector3D& a, const Vector3D& b) { return a + b; }) + .def("__sub__", [](const Vector3D& a, const Vector3D& b) { return a - b; }) + .def("__mul__", [](const Vector3D& v, double s) { return v * s; }) + .def("__truediv__", [](const Vector3D& v, double s) { return v / s; }) + .def("__neg__", [](const Vector3D& v) { return -v; }) + .def("norm", [](const Vector3D& v) { return v.norm(); }) + .def("squared_norm", [](const Vector3D& v) { return v.squaredNorm(); }) + .def("dot", [](const Vector3D& a, const Vector3D& b) { return a.dot(b); }) + .def("cross", [](const Vector3D& a, const Vector3D& b) { return a.cross(b); }) + .def("normalized", [](const Vector3D& v) { return v.normalized(); }); + + // ── Point2D ────────────────────────────────────── + py::class_(core, "Point2D") + .def(py::init(), + py::arg("x") = 0.0, py::arg("y") = 0.0) + .def(py::init([](const py::list& lst) { + if (py::len(lst) != 2) + throw std::runtime_error("Point2D requires exactly 2 values"); + return Point2D(lst[0].cast(), lst[1].cast()); + })) + .def_property_readonly("x", [](const Point2D& p) -> double { return p.x(); }) + .def_property_readonly("y", [](const Point2D& p) -> double { return p.y(); }) + .def("__repr__", &point2d_repr) + .def("__len__", [](const Point2D&) { return 2; }) + .def("__getitem__", [](const Point2D& p, int i) { + if (i < 0 || i > 1) throw py::index_error(); + return p(i); + }); + + // ── AABB3D ─────────────────────────────────────── + py::class_>(core, "AABB3D", + R"pbdoc( +Axis-aligned bounding box with double precision. + +An empty box has min > max. Call ``expand()`` or construct with min/max points. +)pbdoc") + .def(py::init<>(), "Create an empty AABB") + .def(py::init(), + py::arg("min"), py::arg("max"), + "Create from min and max corners") + .def("expand", py::overload_cast(&AABB::expand), + py::arg("point"), "Expand to include a point") + .def("expand", py::overload_cast&>(&AABB::expand), + py::arg("other"), "Expand to include another AABB") + .def_property_readonly("center", &AABB::center) + .def_property_readonly("extent", &AABB::extent) + .def_property_readonly("surface_area", &AABB::surface_area) + .def_property_readonly("volume", &AABB::volume) + .def_property_readonly("min", &AABB::min) + .def_property_readonly("max", &AABB::max) + .def("contains", &AABB::contains, + py::arg("point"), "Check if point is inside (inclusive)") + .def("intersects", &AABB::intersects, + py::arg("other"), "Check if AABBs overlap") + .def("__repr__", [](const AABB& box) { + return "AABB3D(min=" + point3d_repr(box.min()) + + ", max=" + point3d_repr(box.max()) + ")"; + }); + + // ── Polygon2D ──────────────────────────────────── + py::class_(core, "Polygon2D", + R"pbdoc( +2D polygon defined by a list of vertices. +)pbdoc") + .def(py::init<>()) + .def(py::init>(), + py::arg("vertices"), "Construct from vertex list") + .def_property_readonly("vertices", &Polygon2D::vertices) + .def("__len__", &Polygon2D::size) + .def_property_readonly("size", &Polygon2D::size) + .def_property_readonly("area", &Polygon2D::area) + .def_property_readonly("signed_area", &Polygon2D::signed_area) + .def_property_readonly("perimeter", &Polygon2D::perimeter) + .def_property_readonly("centroid", &Polygon2D::centroid) + .def_property_readonly("bounding_box", &Polygon2D::bounding_box) + .def_property_readonly("is_ccw", &Polygon2D::is_ccw) + .def("contains", &Polygon2D::contains, + py::arg("point"), "Point-in-polygon test (ray casting)") + .def("simplify", &Polygon2D::simplify, + py::arg("tolerance"), "Douglas-Peucker simplification") + .def("triangulate", &Polygon2D::triangulate, + "Ear-clipping triangulation. Returns list of [i, j, k] triples.") + .def_static("convex_hull_2d", &Polygon2D::convex_hull_2d, + py::arg("points"), "Andrew's monotone chain convex hull") + .def("__repr__", [](const Polygon2D& poly) { + return "Polygon2D(vertices=" + std::to_string(poly.size()) + ")"; + }); + + // ── Free functions ─────────────────────────────── + core.def("distance", py::overload_cast(&distance), + py::arg("a"), py::arg("b"), + "Euclidean distance between two 3D points"); + + core.def("convex_hull_2d", &convex_hull_2d, + py::arg("points"), + "Graham scan 2D convex hull. Returns vertices in CCW order."); + + core.def("convex_hull_3d", &convex_hull_3d, + py::arg("points"), + "QuickHull 3D convex hull. Returns triangles as list of [p0, p1, p2]."); + + // ── Transform helpers ──────────────────────────── + core.def("translate", py::overload_cast(&translate), + py::arg("v"), "Create translation transform"); + core.def("translate", py::overload_cast(&translate), + py::arg("x"), py::arg("y"), py::arg("z"), + "Create translation transform from x, y, z"); + core.def("rotate", &rotate, + py::arg("axis"), py::arg("angle_rad"), + "Create rotation transform around axis by angle (radians)"); + core.def("rotate_x", &rotate_x, py::arg("rad"), + "Create rotation around X axis"); + core.def("rotate_y", &rotate_y, py::arg("rad"), + "Create rotation around Y axis"); + core.def("rotate_z", &rotate_z, py::arg("rad"), + "Create rotation around Z axis"); + core.def("scale", py::overload_cast(&scale), + py::arg("sx"), py::arg("sy"), py::arg("sz"), + "Create non-uniform scale transform"); + core.def("scale", py::overload_cast(&scale), + py::arg("s"), "Create uniform scale transform"); +} diff --git a/python/src/bind_curves.cpp b/python/src/bind_curves.cpp new file mode 100644 index 0000000..16d50aa --- /dev/null +++ b/python/src/bind_curves.cpp @@ -0,0 +1,107 @@ +#include +#include +#include + +#include "vde/curves/bezier_curve.h" +#include "vde/curves/bspline_curve.h" +#include "vde/curves/nurbs_curve.h" + +namespace py = pybind11; +using namespace vde::curves; + +void bind_curves(py::module& m) { + auto curves = m.def_submodule("curves", + "Parametric curves: Bézier, B-Spline, NURBS"); + + // ── BezierCurve ────────────────────────────────── + py::class_(curves, "BezierCurve", + R"pbdoc( +Bézier curve defined by control points. + +Construct with a list of Point3D control points. The degree is ``len(points) - 1``. +Evaluate at parameter ``t ∈ [0, 1]``. +)pbdoc") + .def(py::init>(), + py::arg("control_points"), + "Construct from control points (list of Point3D)") + .def("evaluate", &BezierCurve::evaluate, + py::arg("t"), + "Evaluate curve at parameter t ∈ [0, 1]") + .def("derivative", &BezierCurve::derivative, + py::arg("t"), py::arg("order") = 1, + "Evaluate derivative of given order at t") + .def_property_readonly("degree", &BezierCurve::degree) + .def_property_readonly("control_points", &BezierCurve::control_points) + .def_property_readonly("domain", &BezierCurve::domain) + .def("split", &BezierCurve::split, + py::arg("t"), + "Split at parameter t, returns (left, right) curves") + .def("degree_elevate", &BezierCurve::degree_elevate, + py::arg("times") = 1, + "Degree elevation, returns new curve") + .def("__repr__", [](const BezierCurve& c) { + return "BezierCurve(degree=" + std::to_string(c.degree()) + ")"; + }); + + // ── BSplineCurve ───────────────────────────────── + py::class_(curves, "BSplineCurve", + R"pbdoc( +B-Spline curve. + +Construct with ``BSplineCurve(control_points, knots, degree)``. +)pbdoc") + .def(py::init, std::vector, int>(), + py::arg("control_points"), py::arg("knots"), py::arg("degree"), + "Construct from control points, knot vector, and degree") + .def("evaluate", &BSplineCurve::evaluate, + py::arg("t"), + "Evaluate curve at parameter t ∈ domain") + .def("derivative", &BSplineCurve::derivative, + py::arg("t"), py::arg("order") = 1, + "Evaluate derivative of given order at t") + .def_property_readonly("degree", &BSplineCurve::degree) + .def_property_readonly("control_points", &BSplineCurve::control_points) + .def_property_readonly("knots", &BSplineCurve::knots) + .def_property_readonly("domain", &BSplineCurve::domain) + .def("find_span", &BSplineCurve::find_span, + py::arg("t"), + "Find knot span index for parameter t") + .def("basis_functions", &BSplineCurve::basis_functions, + py::arg("t"), py::arg("span") = -1, + "Evaluate non-zero basis functions at t") + .def("basis_derivatives", &BSplineCurve::basis_derivatives, + py::arg("t"), py::arg("span") = -1, + "Evaluate first derivatives of non-zero basis functions at t") + .def("__repr__", [](const BSplineCurve& c) { + return "BSplineCurve(degree=" + std::to_string(c.degree()) + + ", cp=" + std::to_string(c.control_points().size()) + ")"; + }); + + // ── NurbsCurve ────────────────────────────────── + py::class_(curves, "NurbsCurve", + R"pbdoc( +Non-Uniform Rational B-Spline (NURBS) curve. + +Construct with ``NurbsCurve(control_points, knots, weights, degree)``. +)pbdoc") + .def(py::init, std::vector, + std::vector, int>(), + py::arg("control_points"), py::arg("knots"), + py::arg("weights"), py::arg("degree"), + "Construct from control points, knot vector, weights, and degree") + .def("evaluate", &NurbsCurve::evaluate, + py::arg("t"), + "Evaluate rational curve at parameter t ∈ domain") + .def("derivative", &NurbsCurve::derivative, + py::arg("t"), py::arg("order") = 1, + "Evaluate derivative of given order at t") + .def_property_readonly("degree", &NurbsCurve::degree) + .def_property_readonly("control_points", &NurbsCurve::control_points) + .def_property_readonly("weights", &NurbsCurve::weights) + .def_property_readonly("knots", &NurbsCurve::knots) + .def_property_readonly("domain", &NurbsCurve::domain) + .def("__repr__", [](const NurbsCurve& c) { + return "NurbsCurve(degree=" + std::to_string(c.degree()) + + ", cp=" + std::to_string(c.control_points().size()) + ")"; + }); +} diff --git a/python/src/bind_mesh.cpp b/python/src/bind_mesh.cpp new file mode 100644 index 0000000..8a0e54f --- /dev/null +++ b/python/src/bind_mesh.cpp @@ -0,0 +1,103 @@ +#include +#include +#include + +#include "vde/mesh/delaunay_2d.h" +#include "vde/mesh/delaunay_3d.h" +#include "vde/mesh/halfedge_mesh.h" + +namespace py = pybind11; +using namespace vde::mesh; + +void bind_mesh(py::module& m) { + auto mesh = m.def_submodule("mesh", + "Mesh generation and processing"); + + // ── DelaunayResult (2D) ────────────────────────── + py::class_(mesh, "DelaunayResult", + R"pbdoc( +Result of 2D Delaunay triangulation. + +Attributes: + vertices: list of Point2D + triangles: list of [i, j, k] index triples +)pbdoc") + .def(py::init<>()) + .def_readwrite("vertices", &DelaunayResult::vertices) + .def_readwrite("triangles", &DelaunayResult::triangles) + .def("__repr__", [](const DelaunayResult& r) { + return "DelaunayResult(v=" + std::to_string(r.vertices.size()) + + ", t=" + std::to_string(r.triangles.size()) + ")"; + }); + + // ── TetrahedronMesh (3D) ───────────────────────── + py::class_(mesh, "TetrahedronMesh", + R"pbdoc( +Result of 3D Delaunay tetrahedralization. + +Attributes: + vertices: list of Point3D + tetrahedra: list of [i, j, k, l] index quadruples +)pbdoc") + .def(py::init<>()) + .def_readwrite("vertices", &TetrahedronMesh::vertices) + .def_readwrite("tetrahedra", &TetrahedronMesh::tetrahedra) + .def("__repr__", [](const TetrahedronMesh& m) { + return "TetrahedronMesh(v=" + std::to_string(m.vertices.size()) + + ", tet=" + std::to_string(m.tetrahedra.size()) + ")"; + }); + + // ── Delaunay free functions ────────────────────── + mesh.def("delaunay_2d", &delaunay_2d, + py::arg("points"), + "Bowyer-Watson 2D Delaunay triangulation"); + + mesh.def("delaunay_3d", &delaunay_3d, + py::arg("points"), + "3D Delaunay tetrahedralization (Bowyer-Watson)"); + + // ── HalfedgeMesh ───────────────────────────────── + py::class_(mesh, "HalfedgeMesh", + R"pbdoc( +Half-edge data structure for polygon meshes. + +Supports construction from triangle soup and provides topology queries, +normal computation, and iterators. +)pbdoc") + .def(py::init<>()) + .def("clear", &HalfedgeMesh::clear) + .def("add_vertex", &HalfedgeMesh::add_vertex, + py::arg("point"), "Add a vertex, returns index") + .def("add_face", &HalfedgeMesh::add_face, + py::arg("vertex_indices"), "Add a face from vertex indices, returns face index") + .def("build_from_triangles", &HalfedgeMesh::build_from_triangles, + py::arg("vertices"), py::arg("triangles"), + "Build from vertex list and triangle index list") + .def_property_readonly("num_vertices", &HalfedgeMesh::num_vertices) + .def_property_readonly("num_faces", &HalfedgeMesh::num_faces) + .def_property_readonly("num_edges", &HalfedgeMesh::num_edges) + .def("vertex", &HalfedgeMesh::vertex, py::arg("idx")) + .def("set_vertex", &HalfedgeMesh::set_vertex, + py::arg("idx"), py::arg("point")) + .def("face_vertices", &HalfedgeMesh::face_vertices, + py::arg("face_index"), + "Get vertex indices of a face") + .def("vertex_faces", &HalfedgeMesh::vertex_faces, + py::arg("vertex_index"), + "Get face indices incident to a vertex") + .def("is_boundary_edge", &HalfedgeMesh::is_boundary_edge, + py::arg("halfedge_index")) + .def("is_boundary_vertex", &HalfedgeMesh::is_boundary_vertex, + py::arg("vertex_index")) + .def("face_normal", &HalfedgeMesh::face_normal, + py::arg("face_index")) + .def("vertex_normal", &HalfedgeMesh::vertex_normal, + py::arg("vertex_index")) + .def("update_normals", &HalfedgeMesh::update_normals) + .def_property_readonly("bounds", &HalfedgeMesh::bounds) + .def("__repr__", [](const HalfedgeMesh& m) { + return "HalfedgeMesh(v=" + std::to_string(m.num_vertices()) + + ", f=" + std::to_string(m.num_faces()) + + ", e=" + std::to_string(m.num_edges()) + ")"; + }); +} diff --git a/python/src/bind_module.cpp b/python/src/bind_module.cpp new file mode 100644 index 0000000..dadd7e0 --- /dev/null +++ b/python/src/bind_module.cpp @@ -0,0 +1,18 @@ +#include +#include +#include + +namespace py = pybind11; + +// Forward declarations of submodule binders +void bind_core(py::module& m); +void bind_curves(py::module& m); +void bind_mesh(py::module& m); + +PYBIND11_MODULE(_vde, m) { + m.doc() = "ViewDesignEngine — CAD computational geometry engine"; + + bind_core(m); + bind_curves(m); + bind_mesh(m); +} diff --git a/python/vde/__init__.py b/python/vde/__init__.py new file mode 100644 index 0000000..a1ccb91 --- /dev/null +++ b/python/vde/__init__.py @@ -0,0 +1,12 @@ +""" +ViewDesignEngine — CAD computational geometry engine for Python. + +Submodules: + vde.core — Point3D, Vector3D, AABB3D, Polygon2D, convex hull, distance, transforms + vde.curves — BezierCurve, BSplineCurve, NurbsCurve + vde.mesh — delaunay_2d, delaunay_3d, HalfedgeMesh +""" +from ._vde import core, curves, mesh + +__version__ = "1.0.0" +__all__ = ["core", "curves", "mesh"] diff --git a/python/vde/core.py b/python/vde/core.py new file mode 100644 index 0000000..cb63f5c --- /dev/null +++ b/python/vde/core.py @@ -0,0 +1,41 @@ +""" +vde.core — Core geometry types and operations. + +Re-exports from the native _vde.core submodule: + Point2D, Point3D, Vector3D, AABB3D, Polygon2D + distance, convex_hull_2d, convex_hull_3d + translate, rotate, rotate_x, rotate_y, rotate_z, scale +""" +from vde._vde.core import ( + Point2D, + Point3D, + Vector3D, + AABB3D, + Polygon2D, + distance, + convex_hull_2d, + convex_hull_3d, + translate, + rotate, + rotate_x, + rotate_y, + rotate_z, + scale, +) + +__all__ = [ + "Point2D", + "Point3D", + "Vector3D", + "AABB3D", + "Polygon2D", + "distance", + "convex_hull_2d", + "convex_hull_3d", + "translate", + "rotate", + "rotate_x", + "rotate_y", + "rotate_z", + "scale", +] diff --git a/python/vde/curves.py b/python/vde/curves.py new file mode 100644 index 0000000..a360ab3 --- /dev/null +++ b/python/vde/curves.py @@ -0,0 +1,17 @@ +""" +vde.curves — Parametric curves for CAD. + +Re-exports from the native _vde.curves submodule: + BezierCurve, BSplineCurve, NurbsCurve +""" +from vde._vde.curves import ( + BezierCurve, + BSplineCurve, + NurbsCurve, +) + +__all__ = [ + "BezierCurve", + "BSplineCurve", + "NurbsCurve", +] diff --git a/python/vde/mesh.py b/python/vde/mesh.py new file mode 100644 index 0000000..6b4348f --- /dev/null +++ b/python/vde/mesh.py @@ -0,0 +1,21 @@ +""" +vde.mesh — Mesh generation and processing. + +Re-exports from the native _vde.mesh submodule: + delaunay_2d, delaunay_3d, DelaunayResult, TetrahedronMesh, HalfedgeMesh +""" +from vde._vde.mesh import ( + delaunay_2d, + delaunay_3d, + DelaunayResult, + TetrahedronMesh, + HalfedgeMesh, +) + +__all__ = [ + "delaunay_2d", + "delaunay_3d", + "DelaunayResult", + "TetrahedronMesh", + "HalfedgeMesh", +] diff --git a/src/sketch/constraint_solver.cpp b/src/sketch/constraint_solver.cpp index 4ab8abb..a3e306c 100644 --- a/src/sketch/constraint_solver.cpp +++ b/src/sketch/constraint_solver.cpp @@ -1,97 +1,551 @@ #include "vde/sketch/constraint_solver.h" +#include +#include #include #include #include +#include namespace vde::sketch { +// ══════════════════════════════════════════════════════ +// 几何实体管理 +// ══════════════════════════════════════════════════════ + int ConstraintSolver::add_point(double x, double y, bool fixed) { int id = static_cast(points_.size()); points_.push_back({id, {x, y}, fixed}); return id; } -int ConstraintSolver::add_line(int p0, int p1) { int id=lines_.size(); lines_.push_back({id,p0,p1}); return id; } -int ConstraintSolver::add_circle(int center, double r) { int id=circles_.size(); circles_.push_back({id,center,r}); return id; } + +int ConstraintSolver::add_line(int p0, int p1) { + int id = static_cast(lines_.size()); + lines_.push_back({id, p0, p1}); + return id; +} + +int ConstraintSolver::add_circle(int center, double r) { + int id = static_cast(circles_.size()); + circles_.push_back({id, center, r}); + return id; +} + void ConstraintSolver::add_constraint(ConstraintType t, const std::vector& els, double v) { constraints_.push_back({t, els, v}); } +Point2D ConstraintSolver::get_point(int id) const { + if (id < 0 || id >= static_cast(points_.size())) + throw std::out_of_range("ConstraintSolver::get_point: invalid point id"); + return points_[id].pos; +} + +// ══════════════════════════════════════════════════════ +// 约束残差评估 +// ══════════════════════════════════════════════════════ + double ConstraintSolver::eval_constraint(const Constraint& c) const { - switch(c.type) { + switch (c.type) { case ConstraintType::Horizontal: { - auto& l=lines_[c.elements[0]]; return points_[l.p1].pos.y()-points_[l.p0].pos.y(); } + auto& l = lines_[c.elements[0]]; + return points_[l.p1].pos.y() - points_[l.p0].pos.y(); + } case ConstraintType::Vertical: { - auto& l=lines_[c.elements[0]]; return points_[l.p1].pos.x()-points_[l.p0].pos.x(); } + auto& l = lines_[c.elements[0]]; + return points_[l.p1].pos.x() - points_[l.p0].pos.x(); + } case ConstraintType::Parallel: { - auto& l1=lines_[c.elements[0]],&l2=lines_[c.elements[1]]; - Vector2D d1=points_[l1.p1].pos-points_[l1.p0].pos, d2=points_[l2.p1].pos-points_[l2.p0].pos; - return d1.x()*d2.y()-d1.y()*d2.x(); } + auto& l0 = lines_[c.elements[0]]; + auto& l1 = lines_[c.elements[1]]; + Vector2D d0 = points_[l0.p1].pos - points_[l0.p0].pos; + Vector2D d1 = points_[l1.p1].pos - points_[l1.p0].pos; + return d0.x() * d1.y() - d0.y() * d1.x(); // cross product + } case ConstraintType::Perpendicular: { - auto& l1=lines_[c.elements[0]],&l2=lines_[c.elements[1]]; - Vector2D d1=points_[l1.p1].pos-points_[l1.p0].pos, d2=points_[l2.p1].pos-points_[l2.p0].pos; - return d1.dot(d2); } + auto& l0 = lines_[c.elements[0]]; + auto& l1 = lines_[c.elements[1]]; + Vector2D d0 = points_[l0.p1].pos - points_[l0.p0].pos; + Vector2D d1 = points_[l1.p1].pos - points_[l1.p0].pos; + return d0.dot(d1); + } case ConstraintType::Coincident: { - auto& p1=points_[c.elements[0]],&p2=points_[c.elements[1]]; return (p1.pos-p2.pos).norm(); } + auto& p0 = points_[c.elements[0]]; + auto& p1 = points_[c.elements[1]]; + return (p1.pos - p0.pos).norm(); + } case ConstraintType::EqualLength: { - auto& l1=lines_[c.elements[0]],&l2=lines_[c.elements[1]]; - double a=(points_[l1.p1].pos-points_[l1.p0].pos).norm(); - double b=(points_[l2.p1].pos-points_[l2.p0].pos).norm(); - return a-b; } + auto& l0 = lines_[c.elements[0]]; + auto& l1 = lines_[c.elements[1]]; + double a = (points_[l0.p1].pos - points_[l0.p0].pos).norm(); + double b = (points_[l1.p1].pos - points_[l1.p0].pos).norm(); + return a - b; + } case ConstraintType::Distance: { - auto& p1=points_[c.elements[0]],&p2=points_[c.elements[1]]; - return (p1.pos-p2.pos).norm()-c.value; } + auto& p0 = points_[c.elements[0]]; + auto& p1 = points_[c.elements[1]]; + return (p1.pos - p0.pos).norm() - c.value; + } case ConstraintType::Radius: - return circles_[c.elements[0]].radius-c.value; - default: return 0; + return circles_[c.elements[0]].radius - c.value; + case ConstraintType::FixedPoint: + case ConstraintType::Fixed: + // Fixed 约束由两个方程分别处理，这里返回范数 + { + auto& p = points_[c.elements[0]]; + // target stored as remaining elements [px, py] or via value + double tx = c.elements.size() > 1 ? static_cast(c.elements[1]) : p.pos.x(); + double ty = c.elements.size() > 2 ? static_cast(c.elements[2]) : p.pos.y(); + Point2D target(tx, ty); + return (p.pos - target).norm(); + } + case ConstraintType::Angle: { + auto& l0 = lines_[c.elements[0]]; + auto& l1 = lines_[c.elements[1]]; + Vector2D d0 = points_[l0.p1].pos - points_[l0.p0].pos; + Vector2D d1 = points_[l1.p1].pos - points_[l1.p0].pos; + // Normalize to avoid NaN at zero length + double n0 = d0.norm(), n1 = d1.norm(); + if (n0 < 1e-15 || n1 < 1e-15) return 0.0; + double cross_val = d0.x() * d1.y() - d0.y() * d1.x(); + double dot_val = d0.dot(d1); + return std::atan2(cross_val, dot_val) - c.value; + } + case ConstraintType::Tangent: + case ConstraintType::Concentric: + return 0.0; + } + return 0.0; +} + +// ══════════════════════════════════════════════════════ +// 约束方程计数 +// ══════════════════════════════════════════════════════ + +int ConstraintSolver::count_equations(const Constraint& c) const { + switch (c.type) { + case ConstraintType::Coincident: + case ConstraintType::Fixed: + return 2; // x 和 y 两个独立方程 + default: + return 1; } } +// ══════════════════════════════════════════════════════ +// 变量同步 +// ══════════════════════════════════════════════════════ + +void ConstraintSolver::sync_from_vars(const double* vars, int nv) { + for (int i = 0; i < nv; ++i) { + points_[i].pos = Point2D(vars[2 * i], vars[2 * i + 1]); + } +} + +// ══════════════════════════════════════════════════════ +// Jacobian 填充 +// ══════════════════════════════════════════════════════ + +void ConstraintSolver::fill_jacobian_rows( + const Constraint& c, int ci, + const std::vector& row_map, + double* J_data, int nv, int stride) const { + + int base_row = row_map[ci]; + int n_cols = 2 * nv; + + // 辅助 lambda：在 J(base_row + offset, col) 位置写入值 + auto set_J = [&](int offset, int col, double val) { + if (col >= 0 && col < n_cols) { + J_data[(base_row + offset) * stride + col] = val; + } + }; + + switch (c.type) { + + // ── Horizontal(l): y1 - y0 = 0 ────────────────── + case ConstraintType::Horizontal: { + auto& l = lines_[c.elements[0]]; + set_J(0, 2 * l.p0 + 1, -1.0); + set_J(0, 2 * l.p1 + 1, 1.0); + break; + } + + // ── Vertical(l): x1 - x0 = 0 ──────────────────── + case ConstraintType::Vertical: { + auto& l = lines_[c.elements[0]]; + set_J(0, 2 * l.p0, -1.0); + set_J(0, 2 * l.p1, 1.0); + break; + } + + // ── Parallel(l0, l1): cross(d0, d1) = 0 ───────── + case ConstraintType::Parallel: { + auto& l0 = lines_[c.elements[0]]; + auto& l1 = lines_[c.elements[1]]; + Vector2D d0 = points_[l0.p1].pos - points_[l0.p0].pos; + Vector2D d1 = points_[l1.p1].pos - points_[l1.p0].pos; + // ∂f/∂p0 = [-d1.y, d1.x], ∂f/∂p1 = [d1.y, -d1.x] + // ∂f/∂p2 = [d0.y, -d0.x], ∂f/∂p3 = [-d0.y, d0.x] + set_J(0, 2 * l0.p0, -d1.y()); + set_J(0, 2 * l0.p0 + 1, d1.x()); + set_J(0, 2 * l0.p1, d1.y()); + set_J(0, 2 * l0.p1 + 1, -d1.x()); + set_J(0, 2 * l1.p0, d0.y()); + set_J(0, 2 * l1.p0 + 1, -d0.x()); + set_J(0, 2 * l1.p1, -d0.y()); + set_J(0, 2 * l1.p1 + 1, d0.x()); + break; + } + + // ── Perpendicular(l0, l1): dot(d0, d1) = 0 ───── + case ConstraintType::Perpendicular: { + auto& l0 = lines_[c.elements[0]]; + auto& l1 = lines_[c.elements[1]]; + Vector2D d0 = points_[l0.p1].pos - points_[l0.p0].pos; + Vector2D d1 = points_[l1.p1].pos - points_[l1.p0].pos; + // ∂f/∂p0 = -d1, ∂f/∂p1 = d1, ∂f/∂p2 = d0, ∂f/∂p3 = -d0 + // Actually: f = (x1-x0)*(x3-x2) + (y1-y0)*(y3-y2) + // ∂f/∂x0 = -(x3-x2) = -d1.x, ∂f/∂y0 = -(y3-y2) = -d1.y + // ∂f/∂x1 = (x3-x2) = d1.x, ∂f/∂y1 = (y3-y2) = d1.y + // ∂f/∂x2 = (x1-x0) = d0.x, ∂f/∂y2 = (y1-y0) = d0.y + // ∂f/∂x3 = -(x1-x0) = -d0.x, ∂f/∂y3 = -(y1-y0) = -d0.y + set_J(0, 2 * l0.p0, -d1.x()); + set_J(0, 2 * l0.p0 + 1, -d1.y()); + set_J(0, 2 * l0.p1, d1.x()); + set_J(0, 2 * l0.p1 + 1, d1.y()); + set_J(0, 2 * l1.p0, d0.x()); + set_J(0, 2 * l1.p0 + 1, d0.y()); + set_J(0, 2 * l1.p1, -d0.x()); + set_J(0, 2 * l1.p1 + 1, -d0.y()); + break; + } + + // ── Coincident(p0, p1): x1-x0=0, y1-y0=0 ────── + case ConstraintType::Coincident: { + int p0 = c.elements[0], p1 = c.elements[1]; + // 方程 0: x1 - x0 = 0 + set_J(0, 2 * p0, -1.0); + set_J(0, 2 * p1, 1.0); + // 方程 1: y1 - y0 = 0 + set_J(1, 2 * p0 + 1, -1.0); + set_J(1, 2 * p1 + 1, 1.0); + break; + } + + // ── EqualLength(l0, l1): |d0| - |d1| = 0 ────── + case ConstraintType::EqualLength: { + auto& l0 = lines_[c.elements[0]]; + auto& l1 = lines_[c.elements[1]]; + Vector2D d0 = points_[l0.p1].pos - points_[l0.p0].pos; + Vector2D d1 = points_[l1.p1].pos - points_[l1.p0].pos; + double n0 = d0.norm(), n1 = d1.norm(); + const double eps = 1e-12; + if (n0 < eps) n0 = eps; + if (n1 < eps) n1 = eps; + // ∂f/∂p0 = -d0/n0, ∂f/∂p1 = d0/n0 + // ∂f/∂p2 = d1/n1, ∂f/∂p3 = -d1/n1 + set_J(0, 2 * l0.p0, -d0.x() / n0); + set_J(0, 2 * l0.p0 + 1, -d0.y() / n0); + set_J(0, 2 * l0.p1, d0.x() / n0); + set_J(0, 2 * l0.p1 + 1, d0.y() / n0); + set_J(0, 2 * l1.p0, d1.x() / n1); + set_J(0, 2 * l1.p0 + 1, d1.y() / n1); + set_J(0, 2 * l1.p1, -d1.x() / n1); + set_J(0, 2 * l1.p1 + 1, -d1.y() / n1); + break; + } + + // ── Distance(p0, p1, d): |p1-p0| - d = 0 ────── + case ConstraintType::Distance: { + int p0 = c.elements[0], p1 = c.elements[1]; + Vector2D dp = points_[p1].pos - points_[p0].pos; + double norm = dp.norm(); + if (norm < 1e-12) norm = 1e-12; + set_J(0, 2 * p0, -dp.x() / norm); + set_J(0, 2 * p0 + 1, -dp.y() / norm); + set_J(0, 2 * p1, dp.x() / norm); + set_J(0, 2 * p1 + 1, dp.y() / norm); + break; + } + + // ── Fixed(p, px, py): x-px=0, y-py=0 ────────── + case ConstraintType::Fixed: { + int pid = c.elements[0]; + // px, py 从 value 和额外元素获取 + // 如果只有1个元素，则固定到当前位置 + // 方程 0: x - px = 0 + set_J(0, 2 * pid, 1.0); + // 方程 1: y - py = 0 + set_J(1, 2 * pid + 1, 1.0); + break; + } + + // ── Angle(l0, l1, θ): atan2 - θ = 0 ─────────── + case ConstraintType::Angle: { + auto& l0 = lines_[c.elements[0]]; + auto& l1 = lines_[c.elements[1]]; + Vector2D d0 = points_[l0.p1].pos - points_[l0.p0].pos; + Vector2D d1 = points_[l1.p1].pos - points_[l1.p0].pos; + double c_cross = d0.x() * d1.y() - d0.y() * d1.x(); + double c_dot = d0.dot(d1); + double denom = c_cross * c_cross + c_dot * c_dot; + if (denom < 1e-20) denom = 1e-20; + // ∂f/∂x0 = -(d1_y * c_dot - c_cross * d1_x) / denom + // ∂f/∂y0 = -(-d1_x * c_dot - c_cross * d1_y) / denom + // ∂f/∂x1 = (d1_y * c_dot - c_cross * d1_x) / denom + // ∂f/∂y1 = (-d1_x * c_dot - c_cross * d1_y) / denom + // ∂f/∂x2 = (d0_y * c_dot - c_cross * d0_x) / denom — wait let me re-derive + // Actually for d1 side: f depends symmetrically, with sign flip + // Let me use the direct derivative: + // f = atan2(c, s) where c = d0×d1, s = d0·d1 + // ∂f/∂c = s/(c²+s²), ∂f/∂s = -c/(c²+s²) + // ∂c/∂x0 = -d1.y, ∂c/∂y0 = d1.x, ∂c/∂x1 = d1.y, ∂c/∂y1 = -d1.x + // ∂c/∂x2 = d0.y, ∂c/∂y2 = -d0.x, ∂c/∂x3 = -d0.y, ∂c/∂y3 = d0.x + // ∂s/∂x0 = -d1.x, ∂s/∂y0 = -d1.y, ∂s/∂x1 = d1.x, ∂s/∂y1 = d1.y + // ∂s/∂x2 = d0.x, ∂s/∂y2 = d0.y, ∂s/∂x3 = -d0.x, ∂s/∂y3 = -d0.y + // ∂f/∂x = (∂f/∂c)*(∂c/∂x) + (∂f/∂s)*(∂s/∂x) + double inv_denom = 1.0 / denom; + double df_dc = c_dot * inv_denom; + double df_ds = -c_cross * inv_denom; + + // p0 (l0 start) + set_J(0, 2 * l0.p0, df_dc * (-d1.y()) + df_ds * (-d1.x())); + set_J(0, 2 * l0.p0 + 1, df_dc * ( d1.x()) + df_ds * (-d1.y())); + // p1 (l0 end) + set_J(0, 2 * l0.p1, df_dc * ( d1.y()) + df_ds * ( d1.x())); + set_J(0, 2 * l0.p1 + 1, df_dc * (-d1.x()) + df_ds * ( d1.y())); + // p2 (l1 start) + set_J(0, 2 * l1.p0, df_dc * ( d0.y()) + df_ds * ( d0.x())); + set_J(0, 2 * l1.p0 + 1, df_dc * (-d0.x()) + df_ds * ( d0.y())); + // p3 (l1 end) + set_J(0, 2 * l1.p1, df_dc * (-d0.y()) + df_ds * (-d0.x())); + set_J(0, 2 * l1.p1 + 1, df_dc * ( d0.x()) + df_ds * (-d0.y())); + break; + } + + case ConstraintType::FixedPoint: + // Handled by excluding from Jacobian (fixed = true) + break; + + default: + break; + } +} + +// ══════════════════════════════════════════════════════ +// DOF 统计 +// ══════════════════════════════════════════════════════ + int ConstraintSolver::degrees_of_freedom() const { int dof = static_cast(points_.size()) * 2; - for(auto&p:points_) if(p.fixed) dof-=2; - return dof - static_cast(constraints_.size()); + for (auto& p : points_) + if (p.fixed) dof -= 2; + int eq_count = 0; + for (auto& c : constraints_) + eq_count += count_equations(c); + return dof - eq_count; } +// ══════════════════════════════════════════════════════ +// Newton-Raphson 求解器 +// ══════════════════════════════════════════════════════ + SolverResult ConstraintSolver::solve(int max_iter, double tol) { - SolverResult r; int nv=points_.size(), nc=constraints_.size(); - if(nc==0){ r.converged=true; r.message="No constraints"; for(auto&p:points_) r.points.push_back(p.pos); return r; } + SolverResult result; + int nv = static_cast(points_.size()); + int nc = static_cast(constraints_.size()); - std::vector vars(nv*2); - for(int i=0;i errs(nc); - double max_err=0; - for(int ci=0;ci free_idx; // 自由变量在 vars 中的索引 + std::vector to_free(n_vars, -1); // 全局→自由变量索引 + for (int i = 0; i < nv; ++i) { + if (!points_[i].fixed) { + to_free[2 * i] = static_cast(free_idx.size()); + free_idx.push_back(2 * i); + to_free[2 * i + 1] = static_cast(free_idx.size()); + free_idx.push_back(2 * i + 1); + } + } + int n_free = static_cast(free_idx.size()); - double total=0; - for(int vi=0;vi1e-15){ step-=errs[ci]*h/de; cnt++; } - } - vars[vi*2+d]=orig; - if(cnt>0){ step/=cnt; vars[vi*2+d]+=step; total+=std::abs(step); } + // ── 方程计数 ──────────────────────────────────── + int n_eq = 0; + std::vector row_map(nc + 1, 0); // 约束 ci → 起始行 + for (int ci = 0; ci < nc; ++ci) { + row_map[ci] = n_eq; + n_eq += count_equations(constraints_[ci]); + } + row_map[nc] = n_eq; // 哨兵 + + // ════════════════════════════════════════════════ + // Newton-Raphson 迭代 + // ════════════════════════════════════════════════ + double prev_residual = 1e100; + + for (int iter = 0; iter < max_iter; ++iter) { + // ── a. 同步变量 → 几何实体 ────────────────── + sync_from_vars(vars.data(), nv); + + // ── b. 计算残差向量 ───────────────────────── + Eigen::VectorXd f(n_eq); + double max_err = 0.0; + for (int ci = 0; ci < nc; ++ci) { + const auto& c = constraints_[ci]; + int neq = count_equations(c); + int r0 = row_map[ci]; + if (neq == 1) { + f[r0] = eval_constraint(c); + max_err = std::max(max_err, std::abs(f[r0])); + } else if (c.type == ConstraintType::Coincident) { + int p0 = c.elements[0], p1 = c.elements[1]; + f[r0] = points_[p1].pos.x() - points_[p0].pos.x(); + f[r0 + 1] = points_[p1].pos.y() - points_[p0].pos.y(); + max_err = std::max(max_err, std::max(std::abs(f[r0]), std::abs(f[r0 + 1]))); + } else if (c.type == ConstraintType::Fixed) { + int pid = c.elements[0]; + // target = 当前位置（或从 extra params 解析） + double tx = points_[pid].pos.x(); + double ty = points_[pid].pos.y(); + f[r0] = points_[pid].pos.x() - tx; + f[r0 + 1] = points_[pid].pos.y() - ty; + max_err = std::max(max_err, std::max(std::abs(f[r0]), std::abs(f[r0 + 1]))); } } - for(double&v:vars) v=std::clamp(v,-1e6,1e6); - for(int i=0;i 0 && std::abs(prev_residual - max_err) < tol * 0.1) { + result.converged = true; + result.iterations = iter; + result.residual = max_err; + result.message = "Stalled (step converged)"; + for (int i = 0; i < nv; ++i) + result.points.push_back(points_[i].pos); + return result; + } + prev_residual = max_err; + + // ── e. 构建 Jacobian（只针对自由变量） ────── + // J_free: n_eq × n_free + Eigen::MatrixXd J(n_eq, n_free); + J.setZero(); + + // 先构建完整 J_full，然后提取自由列 + // 使用临时缓冲区 + std::vector J_full_buf(n_eq * n_vars, 0.0); + for (int ci = 0; ci < nc; ++ci) { + fill_jacobian_rows(constraints_[ci], ci, row_map, + J_full_buf.data(), nv, n_vars); + } + + // 提取自由变量列 + for (int eq = 0; eq < n_eq; ++eq) { + for (int fj = 0; fj < n_free; ++fj) { + int col = free_idx[fj]; + J(eq, fj) = J_full_buf[eq * n_vars + col]; + } + } + + // ── f. 求解 J * dx = -f ──────────────────── + // 使用 SVD 伪逆，处理欠定/超定系统 + Eigen::VectorXd dx = J.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(-f); + + // ── g. 阻尼：检查 step 是否改善 ───────────── + double step_norm = dx.norm(); + + // 回溯线搜索 + double alpha = 1.0; + Eigen::VectorXd vars_backup = vars; + for (int ls = 0; ls < 8; ++ls) { + // 应用步长 + for (int fj = 0; fj < n_free; ++fj) { + vars[free_idx[fj]] += alpha * dx[fj]; + } + // 限制范围 + for (int i = 0; i < n_vars; ++i) { + vars[i] = std::clamp(vars[i], -1e8, 1e8); + } + sync_from_vars(vars.data(), nv); + + // 计算新残差 + double new_err = 0.0; + for (int ci = 0; ci < nc; ++ci) { + const auto& c = constraints_[ci]; + int neq = count_equations(c); + if (neq == 1) { + new_err = std::max(new_err, std::abs(eval_constraint(c))); + } else if (c.type == ConstraintType::Coincident) { + int p0 = c.elements[0], p1 = c.elements[1]; + new_err = std::max(new_err, std::abs(points_[p1].pos.x() - points_[p0].pos.x())); + new_err = std::max(new_err, std::abs(points_[p1].pos.y() - points_[p0].pos.y())); + } + } + + if (new_err <= max_err * 1.01) { + break; // 接受步长 + } + // 回退并减半步长 + vars = vars_backup; + alpha *= 0.5; + if (alpha < 1e-8) break; + } + + // ── h. 更新状态 ───────────────────────────── + sync_from_vars(vars.data(), nv); + result.iterations = iter + 1; + result.residual = max_err; } - if(!r.converged) r.message="Max iterations"; - for(auto&p:points_) r.points.push_back(p.pos); - return r; + + // ════════════════════════════════════════════════ + // 未收敛 + // ════════════════════════════════════════════════ + if (!result.converged) { + result.message = "Max iterations reached"; + // 即使未完全收敛，也返回当前最佳解 + double final_err = 0.0; + for (int ci = 0; ci < nc; ++ci) + final_err = std::max(final_err, std::abs(eval_constraint(constraints_[ci]))); + result.residual = final_err; + } + + for (int i = 0; i < nv; ++i) + result.points.push_back(points_[i].pos); + return result; } } // namespace vde::sketch diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 079d314..816e8e4 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -12,3 +12,4 @@ add_subdirectory(spatial) add_subdirectory(collision) add_subdirectory(boolean) add_subdirectory(brep) +add_subdirectory(sketch) diff --git a/tests/sketch/CMakeLists.txt b/tests/sketch/CMakeLists.txt new file mode 100644 index 0000000..a1d4cc3 --- /dev/null +++ b/tests/sketch/CMakeLists.txt @@ -0,0 +1 @@ +add_vde_test(test_constraint_solver)