diff --git a/loop_cgal/__init__.py b/loop_cgal/__init__.py index fcbbe20..d33e977 100644 --- a/loop_cgal/__init__.py +++ b/loop_cgal/__init__.py @@ -1,5 +1,6 @@ from __future__ import annotations +import copy from typing import Tuple import numpy as np @@ -20,12 +21,20 @@ class TriMesh(_TriMesh): Inherits from the base TriMesh class and provides additional functionality. """ - def __init__(self, surface: pv.PolyData): + def __init__(self, surface): + if isinstance(surface, _TriMesh): + # Copy-construct directly from another TriMesh at the C++ level. + # This is used by read_from_file and anywhere a TriMesh-to-TriMesh + # copy is needed without a pyvista round-trip. + super().__init__(surface) + return + # Validate input surface validate_pyvista_polydata(surface, "input surface") # Triangulate to ensure we have triangular faces - surface = surface.triangulate() + if not surface.is_all_triangles: + surface = surface.triangulate() # Extract vertices and triangles verts = np.array(surface.points, dtype=np.float64).copy() @@ -34,6 +43,7 @@ def __init__(self, surface: pv.PolyData): raise ValueError("Invalid surface geometry") super().__init__(verts, faces) + @classmethod def from_vertices_and_triangles( cls, vertices: np.ndarray, triangles: np.ndarray @@ -56,7 +66,10 @@ def from_vertices_and_triangles( # Create a temporary PyVista PolyData object for validation if not validate_vertices_and_faces(vertices, triangles): raise ValueError("Invalid vertices or triangles") - surface = pv.PolyData(vertices, np.hstack((np.full((triangles.shape[0], 1), 3), triangles)).flatten()) + surface = pv.PolyData( + vertices, + np.hstack((np.full((triangles.shape[0], 1), 3), triangles)).flatten(), + ) return cls(surface) def get_vertices_and_triangles( @@ -107,13 +120,62 @@ def vtk( """ return self.to_pyvista(area_threshold, duplicate_vertex_threshold) - def copy(self) -> TriMesh: + @property + def area(self) -> float: + """Surface area computed directly from the CGAL mesh (no pyvista conversion).""" + return self._cgal_area() + + @property + def points(self) -> np.ndarray: + """Vertex coordinates as (n, 3) numpy array (no pyvista conversion).""" + return self._cgal_points() + + @property + def n_cells(self) -> int: + """Number of faces in the CGAL mesh (no pyvista conversion).""" + return self._cgal_n_faces() + + @property + def n_points(self) -> int: + """Number of vertices in the CGAL mesh (no pyvista conversion).""" + return self._cgal_n_vertices() + + @classmethod + def read_from_file(cls, path: str) -> "TriMesh": + """Read a mesh from a binary file written by :meth:`write_to_file`. + + Bypasses the pyvista/VTK stack entirely — much faster for temp files + passed between CGAL operations. """ - Create a copy of the TriMesh. + return cls(_TriMesh._read_from_file(path)) + + def clone(self) -> TriMesh: + """ + Return a deep copy of this TriMesh as a Python ``TriMesh`` instance. + + Uses the C++ copy-constructor binding to clone the CGAL mesh directly, + with no numpy array roundtrip. + """ + return TriMesh(self) + + def copy(self, deep: bool = True) -> TriMesh: + """ + Return a deep copy of this TriMesh. + + Uses the C++-level ``clone()`` to copy the CGAL mesh directly without + a pyvista round-trip, preserving all vertices, faces, and fixed edges. Returns ------- TriMesh - A copy of the TriMesh object. + A deep copy of this TriMesh. """ - return TriMesh(self.to_pyvista()) + return self.clone() + + def __copy__(self) -> TriMesh: + return self.clone() + + def __deepcopy__(self, memo: dict) -> TriMesh: + result = self.clone() + memo[id(self)] = result + return result diff --git a/loop_cgal/bindings.cpp b/loop_cgal/bindings.cpp index cd9fa21..4eee326 100644 --- a/loop_cgal/bindings.cpp +++ b/loop_cgal/bindings.cpp @@ -22,6 +22,15 @@ PYBIND11_MODULE(_loop_cgal, m) py::class_(m, "TriMesh") .def(py::init &, const pybind11::array_t &>(), py::arg("vertices"), py::arg("triangles")) + .def(py::init([](const TriMesh& other) { return other.clone(); }), + py::arg("other"), + "Copy-construct a TriMesh as a deep clone of another (no array roundtrip).") + .def("clip_with_plane", &TriMesh::clipWithPlane, + py::arg("a"), py::arg("b"), py::arg("c"), py::arg("d"), + py::arg("use_exact_kernel") = true, + "Clip the mesh with the halfspace ax+by+cz+d <= 0. " + "Uses PMP::clip(mesh, Plane_3) directly — no corefinement, no skirt construction. " + "Returns the number of faces removed (0 = no-op).") .def("cut_with_surface", &TriMesh::cutWithSurface, py::arg("surface"), py::arg("preserve_intersection") = false, py::arg("preserve_intersection_clipper") = false, @@ -42,6 +51,21 @@ PYBIND11_MODULE(_loop_cgal, m) "Vertex index pairs defining edges to be fixed in mesh when remeshing.") .def("cut_with_implicit_function", &TriMesh::cut_with_implicit_function, py::arg("property"), py::arg("value"),py::arg("cutmode") = ImplicitCutMode::KEEP_POSITIVE_SIDE, - "Cut the mesh with an implicit function defined by vertex properties."); + "Cut the mesh with an implicit function defined by vertex properties.") + .def("_cgal_area", &TriMesh::area, "Surface area computed directly from CGAL mesh.") + .def("_cgal_n_faces", &TriMesh::n_faces, "Number of faces in the CGAL mesh.") + .def("_cgal_n_vertices", &TriMesh::n_vertices, "Number of vertices in the CGAL mesh.") + .def("_cgal_points", &TriMesh::get_points, "Vertex coordinates as (n, 3) numpy array.") + .def("overlaps", &TriMesh::overlaps, py::arg("other"), py::arg("bbox_tol") = 1e-6, + "Return True if this mesh intersects another TriMesh (AABB fast-reject + triangle-level check).") + .def("clone", &TriMesh::clone, + "Return a deep copy of this TriMesh at the C++ level (no pyvista round-trip). " + "Preserves the full fixed-edge set including any user-added edges.") + .def("write_to_file", &TriMesh::write_to_file, py::arg("path"), + "Write the mesh to a compact binary file (LCMESH format). " + "Much faster than converting to pyvista and saving as VTK.") + .def_static("_read_from_file", &TriMesh::read_from_file, py::arg("path"), + "Read a mesh from a binary file written by write_to_file. " + "Returns a _TriMesh — use TriMesh.read_from_file() from Python instead."); } // End of PYBIND11_MODULE \ No newline at end of file diff --git a/loop_cgal/utils.py b/loop_cgal/utils.py index 011395c..06c57d2 100644 --- a/loop_cgal/utils.py +++ b/loop_cgal/utils.py @@ -1,6 +1,8 @@ import pyvista as pv import numpy as np import scipy.sparse as sp + + def validate_pyvista_polydata( surface: pv.PolyData, surface_name: str = "surface" ) -> None: @@ -80,19 +82,10 @@ def validate_vertices_and_faces(verts, faces): # build a ntris x nverts matrix # populate with true for vertex in each triangle # sum rows and if not equal to 3 then it is degenerate - face_idx = np.arange(faces.shape[0]) - face_idx = np.tile(face_idx, (3, 1)).T.flatten() - faces_flat = faces.flatten() - m = sp.coo_matrix( - (np.ones(faces_flat.shape[0]), (faces_flat, face_idx)), - shape=(verts.shape[0], faces.shape[0]), - dtype=bool, - ) - # coo duplicates entries so just make sure its boolean - m = m > 0 - if not np.all(m.sum(axis=0) == 3): - degen_idx = np.where(m.sum(axis=0) != 3)[1] + a, b, c = faces[:, 0], faces[:, 1], faces[:, 2] + if np.any((a == b) | (b == c) | (a == c)): + degen_idx = np.where((a == b) | (b == c) | (a == c)) raise ValueError( f"Surface contains degenerate triangles: {degen_idx} (each triangle must have exactly 3 vertices)" ) - return True \ No newline at end of file + return True diff --git a/src/mesh.cpp b/src/mesh.cpp index 74f83a4..7940391 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -2,7 +2,12 @@ #include "meshutils.h" #include "globals.h" #include "sizet.h" +#include +#include +#include #include +#include +#include #include #include #include @@ -116,6 +121,42 @@ void TriMesh::init() _edge_is_constrained_map = CGAL::make_boolean_property_map(_fixedEdges); } +// Move constructor: _edge_is_constrained_map stores a raw pointer into +// _fixedEdges, so after relocating _fixedEdges we must re-bind the map. +TriMesh::TriMesh(TriMesh&& other) noexcept + : _mesh(std::move(other._mesh)) + , _fixedEdges(std::move(other._fixedEdges)) + , _edge_is_constrained_map(CGAL::make_boolean_property_map(_fixedEdges)) +{} + +TriMesh& TriMesh::operator=(TriMesh&& other) noexcept +{ + if (this != &other) { + _mesh = std::move(other._mesh); + _fixedEdges = std::move(other._fixedEdges); + _edge_is_constrained_map = CGAL::make_boolean_property_map(_fixedEdges); + } + return *this; +} + +// Copy constructor: same issue — copy creates a new _fixedEdges at a new +// address, so the property map must be re-bound to *this* object's set. +TriMesh::TriMesh(const TriMesh& other) + : _mesh(other._mesh) + , _fixedEdges(other._fixedEdges) + , _edge_is_constrained_map(CGAL::make_boolean_property_map(_fixedEdges)) +{} + +TriMesh& TriMesh::operator=(const TriMesh& other) +{ + if (this != &other) { + _mesh = other._mesh; + _fixedEdges = other._fixedEdges; + _edge_is_constrained_map = CGAL::make_boolean_property_map(_fixedEdges); + } + return *this; +} + void TriMesh::add_fixed_edges(const pybind11::array_t &pairs) { if (!CGAL::is_valid_polygon_mesh(_mesh, LoopCGAL::verbose)) @@ -348,46 +389,47 @@ int TriMesh::cutWithSurface(TriMesh &clipper, CGAL::square(target_bb.ymax() - target_bb.ymin()) + CGAL::square(target_bb.zmax() - target_bb.zmin())); - // Pass 1: accumulate per-ver tex outward directions from each adjacent border face. - // For a border halfedge h (source→target), the outward direction is fn × d, - // where fn is the adjacent face normal and d is the normalised edge direction. - // This lies in the face's tangent plane and points away from its interior. + // Copy the clipper and stitch any collocated border vertices first. + // stitch_borders can remove/merge vertices, so all subsequent passes must + // operate on extended_mesh (not clipper._mesh) to keep vertex descriptors valid. + TriangleMesh extended_mesh = clipper._mesh; + PMP::stitch_borders(extended_mesh); + + // Pass 1: accumulate per-vertex outward directions from each border halfedge + // in extended_mesh (post-stitch). For a border halfedge h (source→target), + // the outward direction is fn × d, where fn is the adjacent face normal and + // d is the normalised edge direction. This lies in the face's tangent plane + // and points away from the interior. std::map outward_sum; - for (auto he : clipper._mesh.halfedges()) + for (auto he : extended_mesh.halfedges()) { - if (!clipper._mesh.is_border(he)) continue; + if (!extended_mesh.is_border(he)) continue; - const Point &ps = clipper._mesh.point(clipper._mesh.source(he)); - const Point &pt = clipper._mesh.point(clipper._mesh.target(he)); + const Point &ps = extended_mesh.point(extended_mesh.source(he)); + const Point &pt = extended_mesh.point(extended_mesh.target(he)); Vector d = pt - ps; const double d_len = std::sqrt(d.squared_length()); if (d_len < 1e-10) continue; d = d / d_len; - const auto adj_face = clipper._mesh.face(clipper._mesh.opposite(he)); - const Vector fn = PMP::compute_face_normal(adj_face, clipper._mesh); + const auto adj_face = extended_mesh.face(extended_mesh.opposite(he)); + const Vector fn = PMP::compute_face_normal(adj_face, extended_mesh); const Vector out = CGAL::cross_product(fn, d); const double out_len = std::sqrt(out.squared_length()); if (out_len < 1e-10) continue; - auto vs = clipper._mesh.source(he); - auto vt = clipper._mesh.target(he); + auto vs = extended_mesh.source(he); + auto vt = extended_mesh.target(he); outward_sum.try_emplace(vs, 0.0, 0.0, 0.0); outward_sum.try_emplace(vt, 0.0, 0.0, 0.0); outward_sum[vs] = outward_sum[vs] + out / out_len; outward_sum[vt] = outward_sum[vt] + out / out_len; } - // Pass 2: copy the clipper, merge any collocated border vertices (they - // produce degenerate faces whose normals are near-zero, which would cause - // skirt edges to be silently skipped and leave bridge gaps), then add one - // new vertex per boundary vertex pushed outward by target_diag, and stitch - // a skirt quad (two triangles) per boundary edge. The winding order - // (source, target, target_new) produces normals consistent with the - // adjacent interior face. - TriangleMesh extended_mesh = clipper._mesh; - PMP::stitch_borders(extended_mesh); - + // Pass 2: add one new vertex per boundary vertex pushed outward by target_diag, + // and stitch a skirt quad (two triangles) per boundary edge of extended_mesh. + // The winding order (source, target, target_new) produces normals consistent + // with the adjacent interior face. std::map skirt_vertex; for (auto &[v, dir_sum] : outward_sum) { @@ -401,11 +443,11 @@ int TriMesh::cutWithSurface(TriMesh &clipper, p.z() + target_diag * dir.z())); } - for (auto he : clipper._mesh.halfedges()) + for (auto he : extended_mesh.halfedges()) { - if (!clipper._mesh.is_border(he)) continue; - auto vs = clipper._mesh.source(he); - auto vt = clipper._mesh.target(he); + if (!extended_mesh.is_border(he)) continue; + auto vs = extended_mesh.source(he); + auto vt = extended_mesh.target(he); if (!skirt_vertex.count(vs) || !skirt_vertex.count(vt)) continue; auto vs_new = skirt_vertex[vs]; auto vt_new = skirt_vertex[vt]; @@ -488,6 +530,51 @@ int TriMesh::cutWithSurface(TriMesh &clipper, return faces_before - faces_after; } +int TriMesh::clipWithPlane(double a, double b, double c, double d, bool use_exact_kernel) +{ + if (_mesh.number_of_vertices() == 0 || _mesh.number_of_faces() == 0) + { + std::cerr << "Error: Source mesh is empty!" << std::endl; + return 0; + } + + const int faces_before = static_cast(_mesh.number_of_faces()); + + try + { + if (use_exact_kernel) + { + Exact_K::Plane_3 plane(a, b, c, d); + Exact_Mesh exact = convert_to_exact(*this); + PMP::clip(exact, plane, CGAL::parameters::clip_volume(false)); + set_mesh(convert_to_double_mesh(exact)); + } + else + { + Plane plane(a, b, c, d); + PMP::clip(_mesh, plane, CGAL::parameters::clip_volume(false)); + } + } + catch (const std::exception &e) + { + std::cerr << "Plane clip failed: " << e.what() << std::endl; + } + + const int faces_after = static_cast(_mesh.number_of_faces()); + if (LoopCGAL::verbose && faces_after >= faces_before) + { + std::cerr << "Warning: clipWithPlane removed no faces (before=" + << faces_before << ", after=" << faces_after << ")." << std::endl; + } + // Compact the mesh: PMP::clip marks removed elements as tombstones without + // freeing them. Iterating _mesh.vertices() when has_garbage()==true requires + // scanning every tombstoned slot, so a mesh that has been clipped 6+ times + // can accumulate enough tombstones to make get_points() take seconds. + if (_mesh.has_garbage()) + _mesh.collect_garbage(); + return faces_before - faces_after; +} + NumpyMesh TriMesh::save(double area_threshold, double duplicate_vertex_threshold) { @@ -805,3 +892,145 @@ void TriMesh::cut_with_implicit_function(const std::vector &property, do // Replace internal mesh _mesh = std::move(newmesh); } + +double TriMesh::area() const +{ + return PMP::area(_mesh); +} + +std::size_t TriMesh::n_faces() const +{ + return _mesh.number_of_faces(); +} + +std::size_t TriMesh::n_vertices() const +{ + return _mesh.number_of_vertices(); +} + +pybind11::array_t TriMesh::get_points() const +{ + std::size_t n = _mesh.number_of_vertices(); + pybind11::array_t result({n, std::size_t(3)}); + auto r = result.mutable_unchecked<2>(); + std::size_t i = 0; + for (auto v : _mesh.vertices()) { + const auto& p = _mesh.point(v); + r(i, 0) = p.x(); + r(i, 1) = p.y(); + r(i, 2) = p.z(); + ++i; + } + return result; +} + +TriMesh TriMesh::clone() const +{ + TriMesh result(_mesh); // copies _mesh via CGAL Surface_mesh copy-ctor, then init() sets border edges + // Overwrite with the full fixed-edge set (border edges + any user-added edges). + // CGAL Surface_mesh indices are integer handles and remain valid across copies of the same topology. + result._fixedEdges = _fixedEdges; + result._edge_is_constrained_map = CGAL::make_boolean_property_map(result._fixedEdges); + return result; +} + +// --------------------------------------------------------------------------- +// Native binary file I/O — avoids the pyvista round-trip for temp files. +// +// Format (little-endian, no padding): +// [6 bytes] magic "LCMESH" +// [4 bytes] n_vertices (uint32) +// [4 bytes] n_faces (uint32) +// [nv×24] vertices (3 × float64 each: x, y, z) +// [nf×12] faces (3 × uint32 each: i0, i1, i2) +// --------------------------------------------------------------------------- +void TriMesh::write_to_file(const std::string& path) const +{ + std::ofstream out(path, std::ios::binary); + if (!out) + throw std::runtime_error("TriMesh::write_to_file — cannot open: " + path); + + // Build a compact vertex index map (mesh may have tombstoned entries + // left by CGAL operations that have not been garbage-collected yet). + std::vector valid_verts; + std::unordered_map vmap; + valid_verts.reserve(_mesh.number_of_vertices()); + vmap.reserve(_mesh.number_of_vertices()); + for (auto v : _mesh.vertices()) { + vmap[static_cast(v.idx())] = + static_cast(valid_verts.size()); + valid_verts.push_back(v); + } + + uint32_t nv = static_cast(valid_verts.size()); + uint32_t nf = static_cast(_mesh.number_of_faces()); + + out.write("LCMESH", 6); + out.write(reinterpret_cast(&nv), 4); + out.write(reinterpret_cast(&nf), 4); + + for (auto v : valid_verts) { + const auto& p = _mesh.point(v); + double x = p.x(), y = p.y(), z = p.z(); + out.write(reinterpret_cast(&x), 8); + out.write(reinterpret_cast(&y), 8); + out.write(reinterpret_cast(&z), 8); + } + + for (auto f : _mesh.faces()) { + auto h = _mesh.halfedge(f); + for (int i = 0; i < 3; ++i) { + uint32_t idx = vmap.at(static_cast(_mesh.target(h).idx())); + out.write(reinterpret_cast(&idx), 4); + h = _mesh.next(h); + } + } +} + +TriMesh TriMesh::read_from_file(const std::string& path) +{ + std::ifstream in(path, std::ios::binary); + if (!in) + throw std::runtime_error("TriMesh::read_from_file — cannot open: " + path); + + char magic[6]; + in.read(magic, 6); + if (std::string(magic, 6) != "LCMESH") + throw std::runtime_error("TriMesh::read_from_file — bad magic in: " + path); + + uint32_t nv, nf; + in.read(reinterpret_cast(&nv), 4); + in.read(reinterpret_cast(&nf), 4); + + TriangleMesh mesh; + std::vector verts(nv); + for (uint32_t i = 0; i < nv; ++i) { + double x, y, z; + in.read(reinterpret_cast(&x), 8); + in.read(reinterpret_cast(&y), 8); + in.read(reinterpret_cast(&z), 8); + verts[i] = mesh.add_vertex(Point(x, y, z)); + } + + for (uint32_t i = 0; i < nf; ++i) { + uint32_t i0, i1, i2; + in.read(reinterpret_cast(&i0), 4); + in.read(reinterpret_cast(&i1), 4); + in.read(reinterpret_cast(&i2), 4); + mesh.add_face(verts[i0], verts[i1], verts[i2]); + } + + return TriMesh(std::move(mesh)); // private ctor calls init() → border edges +} + +bool TriMesh::overlaps(const TriMesh& other, double bbox_tol) const +{ + CGAL::Bbox_3 b1 = PMP::bbox(_mesh); + CGAL::Bbox_3 b2 = PMP::bbox(other._mesh); + if (b1.xmax() < b2.xmin() - bbox_tol || b2.xmax() < b1.xmin() - bbox_tol || + b1.ymax() < b2.ymin() - bbox_tol || b2.ymax() < b1.ymin() - bbox_tol || + b1.zmax() < b2.zmin() - bbox_tol || b2.zmax() < b1.zmin() - bbox_tol) { + return false; + } + return PMP::do_intersect(_mesh, other._mesh); +} diff --git a/src/mesh.h b/src/mesh.h index a231bc7..e153e33 100644 --- a/src/mesh.h +++ b/src/mesh.h @@ -28,6 +28,13 @@ class TriMesh TriMesh(const pybind11::array_t &vertices, const pybind11::array_t &triangles); + // Move/copy constructors — must re-bind _edge_is_constrained_map after + // relocation, because it stores a raw pointer into _fixedEdges. + TriMesh(TriMesh&& other) noexcept; + TriMesh& operator=(TriMesh&& other) noexcept; + TriMesh(const TriMesh& other); + TriMesh& operator=(const TriMesh& other); + // Method to cut the mesh with another surface object. // Returns the number of faces removed (0 = no-op / bad cut). int cutWithSurface(TriMesh &surface, @@ -35,6 +42,12 @@ class TriMesh bool preserve_intersection_clipper = false, bool use_exact_kernel = true); + // Clip the mesh with a halfspace defined by the plane ax+by+cz+d=0. + // The negative side (ax+by+cz+d < 0) is kept. + // Returns the number of faces removed (0 = no-op / bad cut). + int clipWithPlane(double a, double b, double c, double d, + bool use_exact_kernel = true); + // Method to remesh the triangle mesh void remesh(bool split_long_edges, double target_edge_length, int number_of_iterations, bool protect_constraints, @@ -45,6 +58,14 @@ class TriMesh void reverseFaceOrientation(); NumpyMesh save(double area_threshold, double duplicate_vertex_threshold); void add_fixed_edges(const pybind11::array_t &pairs); + double area() const; + std::size_t n_faces() const; + std::size_t n_vertices() const; + pybind11::array_t get_points() const; + bool overlaps(const TriMesh& other, double bbox_tol = 1e-6) const; + TriMesh clone() const; + void write_to_file(const std::string& path) const; + static TriMesh read_from_file(const std::string& path); const TriangleMesh& get_mesh() const { return _mesh; } void set_mesh(const TriangleMesh& mesh) { _mesh = mesh; } private: diff --git a/tests/test_file_io.py b/tests/test_file_io.py new file mode 100644 index 0000000..945ee00 --- /dev/null +++ b/tests/test_file_io.py @@ -0,0 +1,255 @@ +"""Tests for TriMesh native binary file I/O (write_to_file / read_from_file). + +These tests isolate the LCMESH round-trip and the clone/copy paths that +depend on the pybind11 copy-constructor, both of which were implicated in +a segfault caused by calling _TriMesh.__init__ on an uninitialized holder. +""" +from __future__ import annotations + +import copy + +import numpy as np +import pyvista as pv +import pytest + +import loop_cgal +from loop_cgal import TriMesh + + +# --------------------------------------------------------------------------- +# Fixtures +# --------------------------------------------------------------------------- + + +@pytest.fixture +def simple_plane(): + """1×1 unit square as pv.PolyData.""" + return pv.Plane(center=(0, 0, 0), direction=(0, 0, 1), i_size=1.0, j_size=1.0) + + +@pytest.fixture +def unit_mesh(simple_plane): + return TriMesh(simple_plane) + + +@pytest.fixture +def large_plane(): + """Denser mesh for round-trip fidelity checks.""" + return pv.Plane( + center=(0, 0, -500), + direction=(0, 1, 0), + i_size=2000.0, + j_size=1000.0, + i_resolution=10, + j_resolution=5, + ).triangulate() + + +# --------------------------------------------------------------------------- +# write_to_file / read_from_file round-trip +# --------------------------------------------------------------------------- + + +def test_write_creates_file(tmp_path, unit_mesh): + path = str(tmp_path / "mesh.lcm") + unit_mesh.write_to_file(path) + import os + assert os.path.exists(path) + assert os.path.getsize(path) > 0 + + +def test_roundtrip_n_points(tmp_path, unit_mesh): + path = str(tmp_path / "mesh.lcm") + unit_mesh.write_to_file(path) + loaded = TriMesh.read_from_file(path) + assert loaded.n_points == unit_mesh.n_points + + +def test_roundtrip_n_cells(tmp_path, unit_mesh): + path = str(tmp_path / "mesh.lcm") + unit_mesh.write_to_file(path) + loaded = TriMesh.read_from_file(path) + assert loaded.n_cells == unit_mesh.n_cells + + +def test_roundtrip_area(tmp_path, unit_mesh): + path = str(tmp_path / "mesh.lcm") + unit_mesh.write_to_file(path) + loaded = TriMesh.read_from_file(path) + assert abs(loaded.area - unit_mesh.area) < 1e-9 + + +def test_roundtrip_vertices(tmp_path, unit_mesh): + path = str(tmp_path / "mesh.lcm") + unit_mesh.write_to_file(path) + loaded = TriMesh.read_from_file(path) + orig_pts = np.sort(unit_mesh.points, axis=0) + load_pts = np.sort(loaded.points, axis=0) + np.testing.assert_allclose(orig_pts, load_pts, atol=1e-10) + + +def test_roundtrip_returns_trimesh_subclass(tmp_path, unit_mesh): + """read_from_file must return a Python TriMesh, not a bare _TriMesh.""" + path = str(tmp_path / "mesh.lcm") + unit_mesh.write_to_file(path) + loaded = TriMesh.read_from_file(path) + assert type(loaded) is TriMesh + + +def test_roundtrip_large_mesh(tmp_path, large_plane): + mesh = TriMesh(large_plane) + path = str(tmp_path / "large.lcm") + mesh.write_to_file(path) + loaded = TriMesh.read_from_file(path) + assert loaded.n_points == mesh.n_points + assert loaded.n_cells == mesh.n_cells + assert abs(loaded.area - mesh.area) < 1e-4 + + +def test_read_bad_magic(tmp_path): + bad = tmp_path / "bad.lcm" + bad.write_bytes(b"BADMAG\x00\x00\x00\x00\x00\x00\x00\x00") + with pytest.raises(Exception, match="[Bb]ad magic|cannot open|LCMESH"): + TriMesh.read_from_file(str(bad)) + + +def test_read_missing_file(tmp_path): + with pytest.raises(Exception): + TriMesh.read_from_file(str(tmp_path / "nonexistent.lcm")) + + +# --------------------------------------------------------------------------- +# Loaded mesh is fully usable (no dangling state) +# --------------------------------------------------------------------------- + + +def test_loaded_mesh_to_pyvista(tmp_path, unit_mesh): + path = str(tmp_path / "mesh.lcm") + unit_mesh.write_to_file(path) + loaded = TriMesh.read_from_file(path) + pv_mesh = loaded.to_pyvista() + assert isinstance(pv_mesh, pv.PolyData) + assert pv_mesh.n_points > 0 + + +def test_loaded_mesh_remesh(tmp_path, large_plane): + """read_from_file result can be remeshed without crashing.""" + mesh = TriMesh(large_plane) + path = str(tmp_path / "mesh.lcm") + mesh.write_to_file(path) + loaded = TriMesh.read_from_file(path) + loaded.remesh(target_edge_length=200.0, number_of_iterations=1) + assert loaded.n_cells > 0 + + +def test_loaded_mesh_cut_with_surface(tmp_path, large_plane): + """read_from_file result can be cut without crashing.""" + mesh = TriMesh(large_plane) + path = str(tmp_path / "mesh.lcm") + mesh.write_to_file(path) + loaded = TriMesh.read_from_file(path) + + clipper_pv = pv.Plane( + center=(0, 0, -500), + direction=(1, 0, 0), + i_size=3000.0, + j_size=1500.0, + i_resolution=4, + j_resolution=4, + ).triangulate() + clipper = TriMesh(clipper_pv) + loaded.cut_with_surface(clipper) + assert loaded.n_cells >= 0 # just must not crash + + +# --------------------------------------------------------------------------- +# clone / copy after file round-trip +# --------------------------------------------------------------------------- + + +def test_clone_of_loaded_mesh(tmp_path, unit_mesh): + path = str(tmp_path / "mesh.lcm") + unit_mesh.write_to_file(path) + loaded = TriMesh.read_from_file(path) + cloned = loaded.clone() + assert type(cloned) is TriMesh + assert cloned.n_cells == loaded.n_cells + assert cloned.n_points == loaded.n_points + + +def test_clone_of_loaded_is_independent(tmp_path, large_plane): + mesh = TriMesh(large_plane) + path = str(tmp_path / "mesh.lcm") + mesh.write_to_file(path) + loaded = TriMesh.read_from_file(path) + cloned = loaded.clone() + n_before = loaded.n_cells + + clipper_pv = pv.Plane( + center=(0, 0, -500), + direction=(1, 0, 0), + i_size=3000.0, + j_size=1500.0, + ).triangulate() + clipper = TriMesh(clipper_pv) + cloned.cut_with_surface(clipper) + + assert loaded.n_cells == n_before, "clone mutation leaked into loaded mesh" + + +def test_deepcopy_of_loaded_mesh(tmp_path, unit_mesh): + path = str(tmp_path / "mesh.lcm") + unit_mesh.write_to_file(path) + loaded = TriMesh.read_from_file(path) + deep = copy.deepcopy(loaded) + assert type(deep) is TriMesh + assert deep.n_cells == loaded.n_cells + + +def test_write_then_read_then_write(tmp_path, unit_mesh): + """Double round-trip must be stable.""" + p1 = str(tmp_path / "a.lcm") + p2 = str(tmp_path / "b.lcm") + unit_mesh.write_to_file(p1) + loaded = TriMesh.read_from_file(p1) + loaded.write_to_file(p2) + reloaded = TriMesh.read_from_file(p2) + assert reloaded.n_points == unit_mesh.n_points + assert reloaded.n_cells == unit_mesh.n_cells + assert abs(reloaded.area - unit_mesh.area) < 1e-9 + + +# --------------------------------------------------------------------------- +# write_to_file after mesh mutation (tombstoned vertices/faces) +# --------------------------------------------------------------------------- + + +def test_write_after_remesh(tmp_path, large_plane): + """write_to_file must handle tombstoned CGAL entries after remesh.""" + mesh = TriMesh(large_plane) + mesh.remesh(target_edge_length=150.0, number_of_iterations=1) + path = str(tmp_path / "remeshed.lcm") + mesh.write_to_file(path) + loaded = TriMesh.read_from_file(path) + assert loaded.n_cells == mesh.n_cells + assert loaded.n_points == mesh.n_points + assert abs(loaded.area - mesh.area) < 1e-3 + + +def test_write_after_cut(tmp_path, large_plane): + """write_to_file must handle tombstoned entries after cut_with_surface.""" + mesh = TriMesh(large_plane) + clipper_pv = pv.Plane( + center=(0, 0, -500), + direction=(1, 0, 0), + i_size=3000.0, + j_size=1500.0, + ).triangulate() + clipper = TriMesh(clipper_pv) + mesh.cut_with_surface(clipper) + + path = str(tmp_path / "cut.lcm") + mesh.write_to_file(path) + loaded = TriMesh.read_from_file(path) + assert loaded.n_cells == mesh.n_cells + assert loaded.n_points == mesh.n_points diff --git a/tests/test_mesh_operations.py b/tests/test_mesh_operations.py index daf6fa6..922bf6c 100644 --- a/tests/test_mesh_operations.py +++ b/tests/test_mesh_operations.py @@ -1,4 +1,5 @@ from __future__ import annotations +import copy import numpy as np import pyvista as pv import pytest @@ -6,6 +7,10 @@ from loop_cgal._loop_cgal import ImplicitCutMode +# --------------------------------------------------------------------------- +# Fixtures +# --------------------------------------------------------------------------- + @pytest.fixture def square_surface(): # Unit square made of two triangles @@ -18,6 +23,180 @@ def clipper_surface(): return pv.Plane(center=(0,0,0),direction=(1,0,0),i_size=2.0,j_size=2.0) +@pytest.fixture +def unit_trimesh(): + """TriMesh built from a known 1×1 unit square.""" + return loop_cgal.TriMesh(pv.Plane(center=(0, 0, 0), direction=(0, 0, 1), i_size=1.0, j_size=1.0)) + + +# --------------------------------------------------------------------------- +# Properties: area, points, n_cells, n_points +# --------------------------------------------------------------------------- + +def test_area_unit_square(unit_trimesh): + """area property should return the surface area directly from CGAL.""" + assert abs(unit_trimesh.area - 1.0) < 1e-6 + + +def test_points_shape(unit_trimesh): + """points property returns (n, 3) float array.""" + pts = unit_trimesh.points + assert isinstance(pts, np.ndarray) + assert pts.ndim == 2 + assert pts.shape[1] == 3 + assert pts.shape[0] == unit_trimesh.n_points + + +def test_n_cells_positive(unit_trimesh): + """n_cells should equal the number of triangles (≥1).""" + assert unit_trimesh.n_cells > 0 + saved_tris = np.array(unit_trimesh.save().triangles) + assert unit_trimesh.n_cells == saved_tris.shape[0] + + +def test_n_points_positive(unit_trimesh): + """n_points should equal the number of vertices (≥3).""" + assert unit_trimesh.n_points >= 3 + saved_verts = np.array(unit_trimesh.save().vertices) + assert unit_trimesh.n_points == saved_verts.shape[0] + + +def test_n_cells_and_n_points_consistent(unit_trimesh): + """n_cells and n_points should both report non-zero consistent values.""" + assert unit_trimesh.n_cells > 0 + assert unit_trimesh.n_points > 0 + + +# --------------------------------------------------------------------------- +# Constructors: from_vertices_and_triangles, get_vertices_and_triangles +# --------------------------------------------------------------------------- + +def test_from_vertices_and_triangles_roundtrip(): + """from_vertices_and_triangles should produce an equivalent mesh.""" + surface = pv.Plane(center=(0, 0, 0), direction=(0, 0, 1), i_size=2.0, j_size=2.0) + original = loop_cgal.TriMesh(surface) + verts, tris = original.get_vertices_and_triangles() + + rebuilt = loop_cgal.TriMesh.from_vertices_and_triangles(verts, tris) + + assert rebuilt.n_cells == original.n_cells + assert rebuilt.n_points == original.n_points + assert abs(rebuilt.area - original.area) < 1e-6 + + +def test_from_vertices_and_triangles_rejects_invalid(): + """from_vertices_and_triangles should raise ValueError for bad input.""" + with pytest.raises((ValueError, Exception)): + loop_cgal.TriMesh.from_vertices_and_triangles( + np.zeros((2, 3)), # only 2 vertices — can't form a triangle + np.array([[0, 1, 2]]), + ) + + +def test_get_vertices_and_triangles_shapes(unit_trimesh): + """get_vertices_and_triangles returns correctly shaped arrays.""" + verts, tris = unit_trimesh.get_vertices_and_triangles() + assert verts.ndim == 2 and verts.shape[1] == 3 + assert tris.ndim == 2 and tris.shape[1] == 3 + + +# --------------------------------------------------------------------------- +# Conversions: to_pyvista, vtk +# --------------------------------------------------------------------------- + +def test_to_pyvista_returns_polydata(unit_trimesh): + pv_mesh = unit_trimesh.to_pyvista() + assert isinstance(pv_mesh, pv.PolyData) + assert pv_mesh.n_points > 0 + assert pv_mesh.n_cells > 0 + + +def test_vtk_alias_matches_to_pyvista(unit_trimesh): + """vtk() is an alias for to_pyvista() and must return identical geometry.""" + via_to_pyvista = unit_trimesh.to_pyvista() + via_vtk = unit_trimesh.vtk() + np.testing.assert_array_equal(via_to_pyvista.points, via_vtk.points) + np.testing.assert_array_equal(via_to_pyvista.faces, via_vtk.faces) + + +# --------------------------------------------------------------------------- +# Deep copy: clone, copy, __copy__, __deepcopy__ +# --------------------------------------------------------------------------- + +def test_clone_returns_trimesh(unit_trimesh): + cloned = unit_trimesh.clone() + assert isinstance(cloned, loop_cgal.TriMesh) + + +def test_clone_has_same_geometry(unit_trimesh): + cloned = unit_trimesh.clone() + assert cloned.n_cells == unit_trimesh.n_cells + assert cloned.n_points == unit_trimesh.n_points + assert abs(cloned.area - unit_trimesh.area) < 1e-10 + + +def test_clone_is_independent(square_surface, clipper_surface): + """Mutating the clone must not affect the original.""" + original = loop_cgal.TriMesh(square_surface) + cloned = original.clone() + n_cells_before = original.n_cells + + clipper = loop_cgal.TriMesh(clipper_surface) + cloned.cut_with_surface(clipper) + + assert original.n_cells == n_cells_before, "clone mutation leaked into original" + + +def test_copy_method_returns_trimesh(unit_trimesh): + c = unit_trimesh.copy() + assert isinstance(c, loop_cgal.TriMesh) + assert c.n_cells == unit_trimesh.n_cells + + +def test_copy_module_shallow(unit_trimesh): + """copy.copy() should use __copy__ and return an independent TriMesh.""" + c = copy.copy(unit_trimesh) + assert isinstance(c, loop_cgal.TriMesh) + assert c.n_cells == unit_trimesh.n_cells + + +def test_copy_module_deep(unit_trimesh): + """copy.deepcopy() should use __deepcopy__ and return an independent TriMesh.""" + c = copy.deepcopy(unit_trimesh) + assert isinstance(c, loop_cgal.TriMesh) + assert c.n_cells == unit_trimesh.n_cells + + +def test_deepcopy_is_independent(square_surface, clipper_surface): + """Mutating the deepcopy must not affect the original.""" + original = loop_cgal.TriMesh(square_surface) + deep = copy.deepcopy(original) + n_cells_before = original.n_cells + + clipper = loop_cgal.TriMesh(clipper_surface) + deep.cut_with_surface(clipper) + + assert original.n_cells == n_cells_before, "deepcopy mutation leaked into original" + + +# --------------------------------------------------------------------------- +# overlaps +# --------------------------------------------------------------------------- + +def test_overlaps_intersecting(): + """Two overlapping planes should report overlaps=True.""" + a = loop_cgal.TriMesh(pv.Plane(center=(0, 0, 0), direction=(0, 0, 1), i_size=2.0, j_size=2.0)) + b = loop_cgal.TriMesh(pv.Plane(center=(0, 0, 0), direction=(1, 0, 0), i_size=2.0, j_size=2.0)) + assert a.overlaps(b) + + +def test_overlaps_separated(): + """Two well-separated meshes must not report overlap.""" + a = loop_cgal.TriMesh(pv.Plane(center=(0, 0, 0), direction=(0, 0, 1), i_size=1.0, j_size=1.0)) + b = loop_cgal.TriMesh(pv.Plane(center=(100, 0, 0), direction=(0, 0, 1), i_size=1.0, j_size=1.0)) + assert not a.overlaps(b) + + # @pytest.mark.parametrize("remesh_kwargs", [ # {"split_long_edges": True, "target_edge_length": 0.2, "number_of_iterations": 1, "protect_constraints": True, "relax_constraints": False}, # {"split_long_edges": False, "target_edge_length": 0.02, "number_of_iterations": 2, "protect_constraints": True, "relax_constraints": True},