[Big Feature] Add initial pybind interface

.gitignore

@@ -18,3 +18,6 @@ robotDescription.cpp
 log_files/
 
 *Example
+*.so
+*.egg-info
+__pycache__

CMakeLists.txt

@@ -50,6 +50,14 @@ else ()
     message(STATUS "Building with Magnum disabled")
 endif ()
 
+option(WITH_PYBIND OFF)
+if (WITH_PYBIND)
+    find_package(pybind11 REQUIRED)
+else ()
+    message(STATUS "Building with pybind disabled")
+endif ()
+
+
 if (${HAVE_SYMENGINE_LLVM})
     message(STATUS "Testing SymEngine LLVM & SBML support - found")
 else ()
@@ -199,6 +207,18 @@ else()
     message(WARNING "No robotDescription.cpp found. No executable will be created.")
 endif()
 
+if (WITH_PYBIND)
+    pybind11_add_module(py_dismech
+            src/app.cpp
+    )
+    target_link_libraries(py_dismech PRIVATE
+            common_sources
+    )
+
+    set_target_properties(py_dismech PROPERTIES
+    LIBRARY_OUTPUT_DIRECTORY "${CMAKE_SOURCE_DIR}/py_dismech"
+)
+endif()
 
 
 option(CREATE_EXAMPLES OFF)

README.md

@@ -41,8 +41,10 @@ Based on the [Discrete Elastic Rods](https://www.cs.columbia.edu/cg/pdfs/143-rod
 ### TODO
 If you'd like DisMech to support a new feature, feel free create an issue and we'll add it to the list here.
 
+#### Ongoing
+- [ ] Expand Python bindings. PR [#13](https://github.com/StructuresComp/dismech-rods/pull/13)
+
 #### High priority
-- [ ] Migrate API to Python via pybind11
 - [ ] Add extension control via natural length manipulation.
 - [ ] Add shell functionality.
 - [ ] Improve robustness of friction.
@@ -177,10 +179,31 @@ Below is an example rendering.
     make -j$(nproc)
     ```
   - For those wishing to customize the rendering settings, users can take a look and adjust the source code in `magnumDERSimulationEnvironment.cpp` accordingly.
+- [Pybind11](https://github.com/pybind/pybind11)
+  - DisMech also offers Python bindings via pybind11. Users can either install by source or via `pip install pybind11`.
+  Afterward, the bindings can be built as follows:
+    ```bash
+    mkdir build && cd build
+    cmake -DWITH_PYBIND=on ..
+    make -j$(nproc)
+    cd ..
+    pip install -e .  # this will install py_dismech
+    ```
+  - Note that when building with older versions of `MKL` with Python, it may be possible to run into the following [static linkage bug](https://bugs.launchpad.net/ubuntu/+source/intel-mkl/+bug/1947626).
+  Users can either upgrade their `MKL` version or specify their `LD_PRELOAD` explicitly to include the following .so files.
+    ```bash
+    export LD_PRELOAD=$MKL_LIB/libmkl_def.so.2:\
+                      $MKL_LIB/libmkl_avx2.so.2:\
+                      $MKL_LIB/libmkl_core.so.2:\
+                      $MKL_LIB/libmkl_intel_lp64.so.2:\
+                      $MKL_LIB/libmkl_intel_thread.so.2:\
+                      /usr/lib/x86_64-linux-gnu/libiomp5.so
+     ```
+  - All available bindings can be seen in `app.cpp`. Users can expect heavy development for expanding the python bindings as time goes on.
 
 ***
 
-### Running Examples
+### Running Examples in C++
 DisMech is setup so that simulation environments can be instantiated using a single cpp file called `robotDescription.cpp`.
 
 Several example of working DisMech simulations can be seen in the `examples/` directory.
@@ -214,6 +237,12 @@ in `/examples`, their Jacobian matrices are small enough that any amount of para
 simulation. Therefore, it is recommended to set `OMP_NUM_THREADS=1` (`dismech.sh` does this as automatically) and see if parallelization is worth it
 for larger systems through profiling.
 
+### Running Examples in Python
+
+Python versions of certain examples can be found in `py_examples/`. To run, simply run the provided scripts:
+```bash
+OMP_NUM_THREADS=1 python py_examples/spider_case/spiderExample.py
+```
 ***
 
 ### Creating Custom Simulation Environments

py_dismech/__init__.py

@@ -0,0 +1 @@
+from .py_dismech import *

py_examples/helix_case/helixExample.py

@@ -0,0 +1,56 @@
+"""
+    Helix Under Gravity Example
+
+    This is a pybind version of helixExample.cpp
+
+"""
+import numpy as np
+import sys
+import py_dismech
+from pathlib import Path
+from functools import partial
+
+sim_manager = py_dismech.SimulationManager()
+
+soft_robots = sim_manager.soft_robots
+sim_params = sim_manager.sim_params
+render_params = sim_manager.render_params
+create_joint = sim_manager.soft_robots.createJoint
+add_to_joint = sim_manager.soft_robots.addToJoint
+add_force = sim_manager.forces.addForce
+
+############################
+
+sim_params.dt = 5e-3
+sim_params.sim_time = 10
+sim_params.dtol = 1e-3
+sim_params.integrator = py_dismech.IMPLICIT_MIDPOINT
+
+render_params.render_scale = 5.0
+render_params.render_per = 5
+
+# Read vertices describing helical shape from a file
+vertices = np.loadtxt(Path(__file__).parents[2] / 'examples/helix_case/helix_configuration.txt')
+
+# Create the helix limb with custom configuration
+radius = 5e-3
+young_mod = 1e7
+density = 1273.52
+poisson = 0.5
+
+add_limb = partial(sim_manager.soft_robots.addLimb, rho=density, rod_radius=radius,
+                   youngs_modulus=young_mod, poisson_ratio=poisson)
+
+# Add the helical structure as a sequential series of vertices
+add_limb(vertices)
+
+# Fix the top end of the helix (locking the first node)
+soft_robots.lockEdge(0, 0)
+
+# Add gravity
+gravity_force = py_dismech.GravityForce(soft_robots, np.array([0.0, 0.0, -9.8]))
+add_force(gravity_force)
+
+# Initialize and run the simulation
+sim_manager.initialize(sys.argv)
+sim_manager.run()

py_examples/spider_case/spiderExample.py

@@ -0,0 +1,88 @@
+"""
+    Spider on Incline Example
+
+    This is a pybind version of spiderExample.cpp
+"""
+
+import numpy as np
+import sys
+import py_dismech
+from functools import partial
+
+sim_manager = py_dismech.SimulationManager()
+
+soft_robots = sim_manager.soft_robots
+sim_params = sim_manager.sim_params
+render_params = sim_manager.render_params
+create_joint = sim_manager.soft_robots.createJoint
+add_to_joint = sim_manager.soft_robots.addToJoint
+add_force = sim_manager.forces.addForce
+
+SIM_FAST = True
+
+############################
+
+sim_params.sim_time = 2
+sim_params.ftol = 1e-3
+
+render_params.render_scale = 5.0
+render_params.show_mat_frames = True
+
+if SIM_FAST:
+    sim_params.dt = 2.5e-3
+    sim_params.max_iter.num_iters = 15
+    sim_params.max_iter.terminate_at_max = False
+    delta = 5e-3
+    nu = 1e-2
+else:
+    sim_params.dt = 1e-3
+    sim_params.adaptive_time_stepping = 7
+    delta = 5e-4
+    nu = 5e-3
+
+
+n = 25
+radius = 5e-3
+young_mod = 3e6
+density = 1180
+poisson = 0.5
+mu = 0.4
+add_limb = partial(sim_manager.soft_robots.addLimb, num_nodes=n, rho=density, rod_radius=radius,
+                   youngs_modulus=young_mod, poisson_ratio=poisson, mu=mu)
+
+add_limb(np.array([ 0.00,  0.00,  0.20]),    np.array([ 0.00,  0.00,  0.10]))
+add_limb(np.array([ 0.00,  0.00,  0.10]),    np.array([ 0.10,  0.00,  0.10]))
+add_limb(np.array([ 0.00,  0.00,  0.10]),    np.array([ 0.00,  0.10,  0.10]))
+add_limb(np.array([ 0.00,  0.00,  0.10]),    np.array([ 0.00, -0.10,  0.10]))
+add_limb(np.array([ 0.00,  0.00,  0.10]),    np.array([-0.10,  0.00,  0.10]))
+add_limb(np.array([ 0.10,  0.00,  0.10]),    np.array([ 0.10,  0.00,  0.00]))
+add_limb(np.array([ 0.00,  0.10,  0.10]),    np.array([ 0.00,  0.10,  0.00]))
+add_limb(np.array([ 0.00, -0.10,  0.10]),    np.array([ 0.00, -0.10,  0.00]))
+add_limb(np.array([-0.10,  0.00,  0.10]),    np.array([-0.10,  0.00,  0.00]))
+
+create_joint(0, -1)
+add_to_joint(0, 1, 0)
+add_to_joint(0, 2, 0)
+add_to_joint(0, 3, 0)
+add_to_joint(0, 4, 0)
+create_joint(1, -1)
+add_to_joint(1, 5, 0)
+create_joint(2, -1)
+add_to_joint(2, 6, 0)
+create_joint(3, -1)
+add_to_joint(3, 7, 0)
+create_joint(4, -1)
+add_to_joint(4, 8, 0)
+
+# Add gravity with a slight x-axis perturbation
+gravity_force = py_dismech.GravityForce(soft_robots, np.array([1.0, 0.0, -9.8]))
+add_force(gravity_force)
+
+# Add floor contact
+floor_z = -0.10
+floor_contact_force = py_dismech.FloorContactForce(soft_robots, delta, nu, floor_z)
+add_force(floor_contact_force)
+
+# Initialize and run the simulation
+sim_manager.initialize(sys.argv)
+sim_manager.run()

setup.py

@@ -0,0 +1,14 @@
+from setuptools import setup, Extension
+import glob
+import os
+
+# Use glob to find the actual .so file with cp310-... in its name
+so_file = glob.glob(os.path.join('py_dismech', 'py_dismech*.so'))[0]  # Adjust path as needed
+
+setup(
+    name='py_dismech',
+    version='0.1',
+    packages=['py_dismech'],
+    package_dir={'py_dismech': 'py_dismech'},
+    package_data={'py_dismech': [os.path.basename(so_file)]},  # Include the .so file
+)