Implicit function clipping

.github/workflows/ci.yml

@@ -0,0 +1,68 @@
+name: CI
+
+on:
+  push:
+    branches: ["**"]
+  pull_request:
+    branches: ["**"]
+jobs:
+  tester:
+    name: Testing 
+    runs-on: ${{ matrix.os }}
+    strategy:
+      fail-fast: false
+      matrix:
+        os: [ubuntu-latest, macos-latest, macos-14,windows-latest]
+        python-version: ${{ fromJSON(vars.PYTHON_VERSIONS)}}
+    env:
+      MACOSX_DEPLOYMENT_TARGET: '15.0' 
+    steps:
+      - uses: actions/checkout@v4
+      - uses: actions/setup-python@v5
+        with:
+            python-version: ${{ matrix.python-version }}
+      - name: Install dependencies
+        run: |
+          python -m pip install --upgrade pip
+      - if: matrix.os == 'ubuntu-latest'
+        run: |
+          sudo apt-get install -y libgl1-mesa-dev libeigen3-dev libcgal-dev
+      - if: ${{ startsWith(matrix.os, 'macos') }}
+        run: |
+          brew install cgal
+      # Cache vcpkg on Windows
+      - name: Restore vcpkg cache
+        if: runner.os == 'Windows'
+        id: cache-vcpkg
+        uses: actions/cache@v4
+        with:
+          path: |
+            C:/vcpkg
+          key: vcpkg-${{ runner.os }}-${{ hashFiles('vcpkg.json') }}
+          restore-keys: |
+            vcpkg-${{ runner.os }}-
+      - name: Install vcpkg dependencies (Windows only)
+        if: runner.os == 'Windows'
+        shell: pwsh
+        run: |
+          $Env:VCPKG_ROOT = "C:/vcpkg"
+          if (!(Test-Path $Env:VCPKG_ROOT)) {
+              git clone https://github.com/microsoft/vcpkg $Env:VCPKG_ROOT
+              & "$Env:VCPKG_ROOT/bootstrap-vcpkg.bat"
+          }
+           # Only install if packages are missing
+          if (!(Test-Path "$Env:VCPKG_ROOT/installed/x64-windows/cgal")) {
+              & "$Env:VCPKG_ROOT/vcpkg.exe" install yasm-tool cgal
+          }
+          echo "VCPKG_ROOT=$Env:VCPKG_ROOT" >> $Env:GITHUB_ENV
+          # now inject your extra CMake flags:
+          $toolchain = "$Env:VCPKG_ROOT/scripts/buildsystems/vcpkg.cmake"
+          echo "CMAKE_ARGS=-DCMAKE_TOOLCHAIN_FILE=${toolchain} -DVCPKG_ROOT=${Env:VCPKG_ROOT}" >> $Env:GITHUB_ENV
+          echo $CMAKE_ARGS
+      - name: Pip install
+        run: |
+           pip install .[test]
+      - name: Run tests
+        run: |
+           pytest 
+

.gitignore

@@ -43,4 +43,5 @@ Thumbs.db
 loop_cgal/__pycache__/__init__.cpython-311.pyc
 build/*
 *.vtk
-*.ply
+*.ply
+build_win64/*

CMakeLists.txt

@@ -25,26 +25,30 @@ add_library(_loop_cgal MODULE
 target_link_libraries(_loop_cgal PRIVATE pybind11::module CGAL::CGAL)
 target_include_directories(_loop_cgal PRIVATE ${CMAKE_SOURCE_DIR}/src)
 
-set_target_properties(_loop_cgal PROPERTIES PREFIX "" SUFFIX ".pyd")
+if(WIN32)
+    set_target_properties(_loop_cgal PROPERTIES PREFIX "" SUFFIX ".pyd")
+else()
+    set_target_properties(_loop_cgal PROPERTIES PREFIX "" SUFFIX ".so")
+endif()
 
 # Windows: copy required DLLs from VCPKG_ROOT
+
+
 if(WIN32)
-    if(NOT DEFINED ENV{VCPKG_ROOT})
-        message(FATAL_ERROR "VCPKG_ROOT environment variable is not set!")
-    endif()
-    set(VCPKG_ROOT $ENV{VCPKG_ROOT})
-    message(STATUS "Using VCPKG_ROOT: ${VCPKG_ROOT}")
-
-    add_custom_command(TARGET _loop_cgal POST_BUILD
-        COMMAND ${CMAKE_COMMAND} -E copy_if_different
-            "${VCPKG_ROOT}/packages/gmp_x64-windows/bin/gmp-10.dll"
-            "$<TARGET_FILE_DIR:_loop_cgal>"
-        COMMAND ${CMAKE_COMMAND} -E copy_if_different
-            "${VCPKG_ROOT}/packages/mpfr_x64-windows/bin/mpfr-6.dll"
-            "$<TARGET_FILE_DIR:_loop_cgal>"
-    )
+    install(CODE "
+        file(INSTALL
+            DESTINATION \"\${CMAKE_INSTALL_PREFIX}/loop_cgal\"
+            TYPE FILE
+            FILES
+                \"${VCPKG_ROOT}/packages/gmp_x64-windows/bin/gmp-10.dll\"
+                \"${VCPKG_ROOT}/packages/mpfr_x64-windows/bin/mpfr-6.dll\"
+        )
+        message(\"[CMake Install] Copied gmp-10.dll and mpfr-6.dll to install directory.\")
+    ")
 endif()
 
+
+
 # Install for Python packages
 install(TARGETS _loop_cgal
     LIBRARY DESTINATION ${CMAKE_INSTALL_PREFIX}/loop_cgal

examples/cut_example.py

@@ -0,0 +1,34 @@
+import numpy as np
+from loop_cgal import TriMesh, set_verbose
+from LoopStructural.datatypes import BoundingBox
+set_verbose(True)
+
+
+bb = BoundingBox(np.zeros(3), np.ones(3))
+grid = bb.structured_grid().vtk()
+grid["scalars"] = grid.points[:, 0]
+surface = grid.contour([0.5])
+property = surface.points[:, 1].copy()
+surface['prop'] = property
+surface.save("input_mesh.vtk")
+
+mesh = TriMesh(surface)
+
+# Define a scalar property that is the Y coordinate of each vertex
+
+# Cut value at y = 0.6; this will split triangles crossing y=0.6
+cut_value = property.mean()
+# mesh.remesh(0.10, protect_constraints=False)
+
+mesh.vtk().save('before_cut.vtk')
+# Invoke the new method (property is python list or numpy array)
+mesh.cut_with_implicit_function(property.tolist(), cut_value)
+mesh.reverse_face_orientation()
+# Convert back to pyvista for visualization
+out = mesh.to_pyvista()
+out.save('cut_mesh.vtk')
+# # Visualize the result
+# pl = pv.Plotter()
+# pl.add_mesh(out, show_edges=True, color='lightgray')
+# pl.add_mesh(pv.PolyData(np.array([[cut_value, -1e3, 0.0],[cut_value,1e3,0.0]])), color='red')
+# pl.show()

loop_cgal/__init__.py

@@ -7,104 +7,86 @@
 import pyvista as pv
 
 from ._loop_cgal import TriMesh as _TriMesh
-from ._loop_cgal import verbose # noqa: F401
+from ._loop_cgal import verbose  # noqa: F401
 from ._loop_cgal import set_verbose as set_verbose
 
-
-
-
-
-def validate_pyvista_polydata(surface: pv.PolyData, surface_name: str = "surface") -> None:
-    """Validate a PyVista PolyData object.
-
-    Parameters
-    ----------
-    surface : pv.PolyData
-        The surface to validate
-    surface_name : str
-        Name of the surface for error messages
-
-    Raises
-    ------
-    ValueError
-        If the surface is invalid
-    """
-    if not isinstance(surface, pv.PolyData):
-        raise ValueError(f"{surface_name} must be a pyvista.PolyData object")
-
-    if surface.n_points == 0:
-        raise ValueError(f"{surface_name} has no points")
-
-    if surface.n_cells == 0:
-        raise ValueError(f"{surface_name} has no cells")
-
-    points = np.asarray(surface.points)
-    if not np.isfinite(points).all():
-        raise ValueError(f"{surface_name} points contain NaN or infinite values")
-
+from .utils import validate_pyvista_polydata, validate_vertices_and_faces
 
 
 class TriMesh(_TriMesh):
     """
     A class for handling triangular meshes using CGAL.
-    
+
     Inherits from the base TriMesh class and provides additional functionality.
     """
+
     def __init__(self, surface: pv.PolyData):
         # Validate input surface
         validate_pyvista_polydata(surface, "input surface")
-        
+
         # Triangulate to ensure we have triangular faces
         surface = surface.triangulate()
-        
+
         # Extract vertices and triangles
         verts = np.array(surface.points, dtype=np.float64).copy()
         faces = surface.faces.reshape(-1, 4)[:, 1:].copy().astype(np.int32)
-
-        # Additional validation on extracted data
-        if verts.size == 0:
-            raise ValueError("Surface has no vertices after triangulation")
-
-        if faces.size == 0:
-            raise ValueError("Surface has no triangular faces after triangulation")
-
-        if not np.isfinite(verts).all():
-            raise ValueError("Surface vertices contain NaN or infinite values")
-
-        # Check triangle indices
-        max_vertex_index = verts.shape[0] - 1
-        if faces.min() < 0:
-            raise ValueError("Surface has negative triangle indices")
-
-        if faces.max() > max_vertex_index:
-            raise ValueError(f"Surface triangle indices exceed vertex count (max index: {faces.max()}, vertex count: {verts.shape[0]})")
-        # Check for degenerate triangles
-        # 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]
-            raise ValueError(f"Surface contains degenerate triangles: {degen_idx} (each triangle must have exactly 3 vertices)")
-
+        if not validate_vertices_and_faces(verts, faces):
+            raise ValueError("Invalid surface geometry")
 
         super().__init__(verts, faces)
-
-    def to_pyvista(self, area_threshold: float = 1e-6,  # this is the area threshold for the faces, if the area is smaller than this it will be removed
-            duplicate_vertex_threshold: float = 1e-4,  # this is the threshold for duplicate vertices
-            ) -> pv.PolyData:
+    @classmethod
+    def from_vertices_and_triangles(
+        cls, vertices: np.ndarray, triangles: np.ndarray
+    ) -> TriMesh:
+        """
+        Create a TriMesh from vertices and triangle indices.
+
+        Parameters
+        ----------
+        vertices : np.ndarray
+            An array of shape (n_vertices, 3) containing the vertex coordinates.
+        triangles : np.ndarray
+            An array of shape (n_triangles, 3) containing the triangle vertex indices.
+
+        Returns
+        -------
+        TriMesh
+            The created TriMesh object.
+        """
+        # 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())
+        return cls(surface)
+
+    def get_vertices_and_triangles(
+        self,
+        area_threshold: float = 1e-6,  # this is the area threshold for the faces, if the area is smaller than this it will be removed
+        duplicate_vertex_threshold: float = 1e-4,
+    ) -> Tuple[np.ndarray, np.ndarray]:
+        """
+        Get the vertices and triangle indices of the TriMesh.
+
+        Returns
+        -------
+        Tuple[np.ndarray, np.ndarray]
+            A tuple containing:
+            - An array of shape (n_vertices, 3) with the vertex coordinates.
+            - An array of shape (n_triangles, 3) with the triangle vertex indices.
+        """
+        np_mesh = self.save(area_threshold, duplicate_vertex_threshold)
+        vertices = np.array(np_mesh.vertices).copy()
+        triangles = np.array(np_mesh.triangles).copy()
+        return vertices, triangles
+
+    def to_pyvista(
+        self,
+        area_threshold: float = 1e-6,  # this is the area threshold for the faces, if the area is smaller than this it will be removed
+        duplicate_vertex_threshold: float = 1e-4,  # this is the threshold for duplicate vertices
+    ) -> pv.PolyData:
         """
         Convert the TriMesh to a pyvista PolyData object.
-        
+
         Returns
         -------
         pyvista.PolyData
@@ -115,3 +97,23 @@ def to_pyvista(self, area_threshold: float = 1e-6,  # this is the area threshold
         triangles = np.array(np_mesh.triangles).copy()
         return pv.PolyData.from_regular_faces(vertices, triangles)
 
+    def vtk(
+        self,
+        area_threshold: float = 1e-6,
+        duplicate_vertex_threshold: float = 1e-4,
+    ) -> pv.PolyData:
+        """
+        Alias for to_pyvista method.
+        """
+        return self.to_pyvista(area_threshold, duplicate_vertex_threshold)
+
+    def copy(self) -> TriMesh:
+        """
+        Create a copy of the TriMesh.
+
+        Returns
+        -------
+        TriMesh
+            A copy of the TriMesh object.
+        """
+        return TriMesh(self.to_pyvista())

loop_cgal/bindings.cpp

@@ -14,12 +14,18 @@ PYBIND11_MODULE(_loop_cgal, m)
          .def(py::init<>())
          .def_readwrite("vertices", &NumpyMesh::vertices)
          .def_readwrite("triangles", &NumpyMesh::triangles);
+     py::enum_<ImplicitCutMode>(m, "ImplicitCutMode")
+         .value("PRESERVE_INTERSECTION", ImplicitCutMode::PRESERVE_INTERSECTION)
+         .value("KEEP_POSITIVE_SIDE", ImplicitCutMode::KEEP_POSITIVE_SIDE)
+         .value("KEEP_NEGATIVE_SIDE", ImplicitCutMode::KEEP_NEGATIVE_SIDE)
+         .export_values();
      py::class_<TriMesh>(m, "TriMesh")
          .def(py::init<const pybind11::array_t<double> &, const pybind11::array_t<int> &>(),
               py::arg("vertices"), py::arg("triangles"))
          .def("cut_with_surface", &TriMesh::cutWithSurface, py::arg("surface"),
               py::arg("preserve_intersection") = false,
-              py::arg("preserve_intersection_clipper") = false)
+              py::arg("preserve_intersection_clipper") = false,
+              py::arg("use_exact_kernel") = true)
          .def("remesh", &TriMesh::remesh, py::arg("split_long_edges") = true,
               py::arg("target_edge_length") = 10.0,
               py::arg("number_of_iterations") = 3,
@@ -31,6 +37,9 @@ PYBIND11_MODULE(_loop_cgal, m)
               "Reverse the face orientation of the mesh.")
          .def("add_fixed_edges", &TriMesh::add_fixed_edges,
               py::arg("pairs"),
-              "Vertex index pairs defining edges to be fixed in mesh when remeshing.");
+              "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.");
 
 } // End of PYBIND11_MODULE