Python interface

.devcontainer/Dockerfile

@@ -14,7 +14,10 @@ RUN apt-get update && apt-get install -y \
     libtbb-dev \
     libboost-all-dev \
     libilmbase-dev \
-    libopenexr-dev 
+    libopenexr-dev \
+    python3-dev \
+    python3-pytest \
+    python3-numpy
 
 # Set working directory
 WORKDIR /workspace

.github/workflows/cmake-single-platform.yml

@@ -23,7 +23,7 @@ jobs:
     - uses: actions/checkout@v4
 
     - name: Install deps
-      run: sudo apt-get install libeigen3-dev libgtest-dev libopenvdb-dev libtbb-dev libboost-all-dev libilmbase-dev libopenexr-dev
+      run: sudo apt-get install libeigen3-dev libgtest-dev libopenvdb-dev libtbb-dev libboost-all-dev libilmbase-dev libopenexr-dev python3-dev python3-pytest python3-numpy
 
     - name: Configure CMake
       # Configure CMake in a 'build' subdirectory. `CMAKE_BUILD_TYPE` is only required if you are using a single-configuration generator such as make.

.gitignore

@@ -3,3 +3,6 @@ build/
 postprocess/backup/
 *.vscode/
 deps/
+.pytest_cache/
+__pycache__/
+*.pyc

CMakeLists.txt

@@ -27,17 +27,41 @@ if(CMAKE_BUILD_TYPE STREQUAL "Release")
     endif()
 endif()
 
-# Check LTO/IPO support — result used by src/CMakeLists.txt to enable it on the mpm target only
+# Check LTO/IPO support — result used below to enable it on both the mpm executable
+# and the matter Python module
 include(CheckIPOSupported)
 check_ipo_supported(RESULT ipo_supported OUTPUT ipo_output)
 if(NOT ipo_supported)
     message(STATUS "IPO/LTO not supported: ${ipo_output}")
 endif()
 # =================================================
 
+# Compile settings shared between the mpm executable (src/) and the matter Python
+# module (python/), so simulation code is compiled identically either way — single
+# source of truth, referenced by both src/CMakeLists.txt and python/CMakeLists.txt.
+set(MATTER_EIGEN_DEFS
+    EIGEN_NO_DEBUG                    # Disable debug checks in Eigen for speed
+    EIGEN_UNROLLING_LIMIT=256         # Unroll small fixed-size matrix loops (raised from 128)
+    EIGEN_FAST_MATH=1                 # Use fast math approximations where possible
+    EIGEN_VECTORIZE_FMA               # Enable fused multiply-add instructions
+)
+set(MATTER_RELEASE_FLAGS
+    # SPECIFIC RELEASE OPTIMIZATIONS (safer fast-math subset — preserves NaN/Inf detection)
+    $<$<CONFIG:Release>:-fno-math-errno>
+    $<$<CONFIG:Release>:-fno-trapping-math>
+    $<$<CONFIG:Release>:-fno-signed-zeros>
+    $<$<CONFIG:Release>:-fassociative-math>
+    $<$<CONFIG:Release>:-freciprocal-math>
+    $<$<CONFIG:Release>:-funroll-loops>
+    $<$<CONFIG:Release>:-fomit-frame-pointer>
+)
+
 option (USE_VDB "Use OpenVDB part of the code" ON)
 message(STATUS "Use of OpenVDB: ${USE_VDB}")
 
+option (BUILD_PYTHON_BINDINGS "Build pybind11 Python bindings for Simulation" ON)
+message(STATUS "Build Python bindings: ${BUILD_PYTHON_BINDINGS}")
+
 enable_testing()
 
 # Custom target to clean CMake cache (like make clean for cache)
@@ -49,6 +73,10 @@ add_custom_target(clean-cache
 add_subdirectory(src)
 add_subdirectory(tests)
 
+if (BUILD_PYTHON_BINDINGS)
+    add_subdirectory(python)
+endif()
+
 #### Download header files from external Github repos ####
 
 # Create deps folder

README.md

@@ -71,57 +71,72 @@ With _Matter_, you can simulate granular flow on various simple and complex terr
 
 * Initial particle positions can be sampled using the **Poisson disk sampling** scheme by R. Bridson, ACM SIGGRAPH 2007, here based on the implementation by [Tommy Hinks](https://github.com/thinks/poisson-disk-sampling)
 
+* Python interface using [pybind11](https://github.com/pybind/pybind11)   
+
 
 ## Get started
 
 ### Installing dependencies
 
 The required non-header-only dependencies are **[CMake](https://cmake.org/)**, **[OpenMP](https://www.openmp.org/)**. The standard C++ linear algebra template library **[Eigen](https://eigen.tuxfamily.org/)** is required. **[OpenVDB](https://www.openvdb.org/)** is optional, however, needed if vdb-files are used as input (either as terrain/objects or as particle sampling domain). 
 
-##### Linux
+<details>
+<summary><b>Linux</b></summary>
+
 Install all the above dependencies with this simple command:        
 `sudo apt-get install -y cmake libeigen3-dev libopenvdb-dev libtbb-dev libboost-all-dev libilmbase-dev libopenexr-dev`    
 where only `cmake` and `libeigen3-dev` are strictly required. 
+</details>
+
+<details>
+<summary><b>MacOS</b></summary>
 
-##### MacOS
 You can download/install the required dependencies with [Homebrew](https://brew.sh/) through   
 `brew install cmake eigen libomp openvdb`    
 although OpenVDB may require other dependecies, see https://formulae.brew.sh/formula/openvdb. 
+</details>
 
 
-However, if not using Linux, we recommend instead using [VSCode Dev Containers](https://code.visualstudio.com/docs/devcontainers/tutorial). 
+We recommend in any case using [VSCode Dev Containers](https://code.visualstudio.com/docs/devcontainers/tutorial). 
 Download [VSCode](https://code.visualstudio.com/) and [Docker](https://www.docker.com/products/docker-desktop/). In VSCode, install the "Dev Containers" extension. When you open the Matter repository in VSCode, you will be asked if you want to run it in a container. When you accept, a container will be set up and you are ready to simulate!
 
 
 ### Build and run the code
 
-1. Set up your simulation parameters and initial state in `mpm.cpp`.  In `tools.hpp`, the user can specify the dimension of the simulation (default: 2D) and order of interpolation (default: quadratic). 
+1. Specify the dimension of the simulation (default: 2D) and order of interpolation (default: quadratic) in `src/tools.hpp`
 
-2. Create a build directory:     
+2. From the top-level `matter/`, create a build directory:     
    `mkdir build`
 
 3. Enter the build directory:    
    `cd build`
 
 4. Specify CMake options. Use default options (release mode and use OpenVDB):     
    `cmake ..`     
-   or optionally specify     
+   or optionally specify particular build-types, e.g., debug mode and without OpenVDB    
    `cmake -DCMAKE_BUILD_TYPE=Debug -DUSE_VDB=OFF ..`       
 
 5. Compile (NB: the number of threads for the _simulation_ is specified in `mpm.cpp`)      
    `make -j <number of cores for compilation>` 
 
-6. Run the sim you set up in `mpm.cpp`:      
-   `./src/mpm`
+6. Run the simulation:       
+   * If using C++ interface: Run the sim you set up in `mpm.cpp`:      
+   `./src/mpm`     
+   * If using Python interface: Run the sim in any python file:     
+   `PYTHONPATH=<path_to_matter>/build/python python3 my_script.py`
 
 7. [Optional] Tests:       
    Run all tests with `make test` or `ctest`. Run single test with `ctest -R <name of test>`
 
 You can "clean" using `make clean` and `make clean-cache`, e.g., if you want to change CMake options.
 
+
 ### Example of setup file
 
-Here is a minimal setup file `mpm.cpp` for a simple granular collapse.    
+Here is a minimal setup for a simple granular collapse, shown both as a C++ setup file (`mpm.cpp`) and as the equivalent Python script (see [Python bindings](#python-bindings) above for how to build/run it).
+
+<details>
+<summary><b>C++ (<code>mpm.cpp</code>)</b></summary>
 
 ```cpp
 #include "tools.hpp"
@@ -167,7 +182,45 @@ int main(){
 }
 ```
 
-In the `examples` folder, other examples will be archived (to use one of these examples, simply copy it into the `src` folder and rename it `mpm.cpp`).
+</details>
+
+<details>
+<summary><b>Python</b></summary>
+
+```python
+import matter
+
+sim = matter.Simulation()
+sim.initialize(True, "output/", "collapse")
+
+sim.end_frame = 20  # last frame to simulate
+sim.fps = 10         # frames per second
+sim.n_threads = 8    # number of threads in parallel
+
+sim.Lx = 1
+sim.Ly = 1
+# sim.Lz = 1  # only 3D, otherwise remove
+
+matter.sample_particles(sim, 0.01)  # sampling radius
+
+sim.add_plate(matter.ObjectPlate(0, matter.PlateType.bottom, matter.BC.NoSlip, 0.5))
+
+sim.rho = 1000              # density (kg/m3)
+sim.gravity = [0, -9.81]    # gravity
+
+sim.E = 1e6    # Young's modulus (Pa)
+sim.nu = 0.3   # Poisson's ratio (-)
+
+sim.plastic_model = matter.PlasticModel.DP  # Drucker-Prager yield
+sim.M = 0.5          # internal friction
+sim.q_cohesion = 0   # cohesion
+
+sim.simulate()
+```
+
+</details>
+
+In the `examples` and `python/examples/` folder, other examples will be archived (if you are using C++, copy the example `.cpp` files into the `src` folder and rename it `mpm.cpp`).
 
 ### Objects and terrains
 
@@ -177,12 +230,15 @@ Rigid objects and terrains (boundaries) are either
 
 The objects already available can be found in the `objects` directory. The forces acting on these objects can be measured by setting their parameter `force_calc = true`, see the example `collapse.cpp`. Be aware that if the force measurement is turned on for an object which is contantly in collision with a large part of the material, this may slightly slow down the simulation.
 
-##### Analytical objects
+<details>
+<summary><b>Analytical objects</b></summary>
 Analytical objects can be specified as a derived class from the general `ObjectGeneral` class. An example of this is `ObjectBump` which provides the terrain of a smooth bump used in the flow experiments in [Viroulet et al. (2017)](https://doi.org/10.1017/jfm.2017.41). For the very common case of an axis-aligned plate, an `ObjectPlate` class has been made separate from `ObjectGeneral` class for efficiency and convenience. In `ObjectPlate`, you can also assign a speed to the plate, as well as controls on the time-evolution of the speed. Any plate must either a `PlateType::top`, `PlateType::bottom`, `PlateType::front`, `PlateType::back`, `PlateType::left` or `PlateType::right`. 
-
-##### OpenVDB objects
+</details>
+<details>
+<summary><b>OpenVDB objects</b></summary>
 A terrain/object from a `.vdb` is stored in an instance of the `ObjectVdb` class which is derived from `ObjectGeneral`.
 Examples of `.vdb`-files are found in the folder `levelsets`.
+</details>
 
 Multiple objects in a simulation are possible. Note that all `ObjectGeneral` instances must be added to the vector `objects` and `ObjectPlate` instances are added to the vector `plates`. 
 
@@ -200,7 +256,11 @@ PLY files can also be easily read by Python. This is shown in the file `load_ply
 If forces are measured on objects, this will be saved as `force_<name-of-object>.csv` where each line contains (time, Fx, Fy(, Fz)) for each frame.
 
 
-### Key parameters and options
+### Key parameters and models
+
+<details>
+<summary><b>Click to expand: Main parameters and options</b></summary>
+
 This is a non-exhaustive list of parameters and options (of the `Simulation` class) to be specified in the input file `mpm.cpp`. See `simulation.hpp` for the complete list, and take advantage of the current `mpm.cpp` example file. Other example files are found in the `examples` directory.
 
 | Parameter  | Default value  | Description  |
@@ -227,8 +287,10 @@ This is a non-exhaustive list of parameters and options (of the `Simulation` cla
 | `nu`                   | 0.3            | The 3D Poisson's ratio (-)
 | `q_prefac` | sqrt(1/2) | Prefactor used in definition of q, default is q = sqrt(1/2 s:s). 
 
+</details>
 
-Here is a list of the various plastic models and their parameters:
+<details>
+<summary><b>Click to expand: Constitutive plastic models</b></summary>
 
 | Model                                | Name                  | Parameters                | Default value   |
 |  ----                                | ----     |    ----      |          ---    |
@@ -285,6 +347,7 @@ In the MCC-based models, one must also choose a corresponding hardening law:
 * exponential implicit law `HardeningLaw::ExpoImpl`     
 * hyperbolic sine implicit law `HardeningLaw::SinhImpl`     
 
+</details>
 
 ## Limitations
 
@@ -306,10 +369,15 @@ Please attribute this software by citing this article:
 
 * Blatny, L. & Gaume, J. (2025). Matter (v1): an open-source MPM solver for granular matter, _Geoscientific Model Development_, 18, 9149–9166. [DOI: 10.5194/gmd-18-9149-2025](https://doi.org/10.5194/gmd-18-9149-2025).
 
-If you use the $\mu(I)$-rheology models, please also cite this article: 
+If you use the $\mu(I)$-rheology models, please cite this article: 
 
 * Blatny, L., Gray, J.M.N.T., & Gaume, J. (2024). A critical state μ(I)-rheology model for cohesive granular flows. _Journal of Fluid Mechanics_, 997, A67. [DOI: 10.1017/jfm.2024.643](https://doi.org/10.1017/jfm.2024.643)    
 
+If you use the viscous models, please cite this article: 
+
+* Blatny, L. & Pellet, A. (2026). Modeling elasto-viscoplastic free-surface flows with different yield surfaces
+. _arXiv preprint_ arXiv:2607.17380. [DOI: 10.48550/arXiv.2607.17380](https://doi.org/10.48550/arXiv.2607.17380)
+
 If you use the sparse background grid feature (by setting `use_sparse=true`), please cite this article:
 
 * Zhao, Y., Blatny, L., Feng, X., Juel, M. M., Jiang, C., & Gaume, J. (2026). Unified sparse framework for large-scale material point method simulations. _arXiv preprint_ arXiv:2605.28525. [DOI: 10.48550/arXiv.2605.28525](https://doi.org/10.48550/arXiv.2605.28525)

examples/README.md

@@ -1,4 +1,13 @@
 # Examples
 
+Here you will find some example setup files. The setup can be done using either C++ or using Python.
+
+### Using C++
+
 An example file can be copied to the `src` directory and renamed `mpm.cpp`.
 (Alternatively, change the `CMakeLists.txt` files)
+
+### Python
+
+Simply run any Python file like this:
+`PYTHONPATH=<path_to_matter>/build/python python3 my_example.py`

examples/block_incline.py

@@ -0,0 +1,82 @@
+# Copyright (C) 2026 Lars Blatny. Released under GPL-3.0 license.
+
+import numpy as np
+import matter
+
+DIM = matter.Simulation().dim
+def vec(x=0.0, y=0.0, z=0.0):
+    return [x, y, z][:DIM]
+
+sim = matter.Simulation()
+
+sim.initialize(True, "output/", "block_incline")
+
+sim.save_grid = True
+sim.end_frame = 10       # last frame to simulate
+sim.fps = 20              # frames per second
+sim.n_threads = 8         # number of threads in parallel
+sim.cfl = 0.5              # CFL constant, typically around 0.5
+sim.flip_ratio = -0.95    # (A)PIC-(A)FLIP ratio in [-1,1].
+sim.reduce_verbose = True
+
+# INITIALIZE ELASTICITY 
+sim.elastic_model = matter.ElasticModel.Hencky
+sim.plastic_model = matter.PlasticModel.NoPlasticity
+sim.E = 1e6      # Young's modulus (Pa)
+sim.nu = 0.3     # Poisson's ratio (-)
+sim.rho = 1000   # Density (kg/m3)
+
+# INCLINE ANGLE AND FRICTION 
+theta = np.radians(24)
+sim.gravity = vec(9.81 * np.sin(theta), -9.81 * np.cos(theta))
+
+# PARTICLES
+block_size = vec(0.1, 0.05, 0.05 if DIM == 3 else 0.1)
+sim.Lx = block_size[0]
+sim.Ly = block_size[1]
+if DIM == 3:
+    sim.Lz = block_size[2]
+
+k_rad = 0.005
+matter.sample_particles(sim, k_rad, ppc=4 if DIM == 2 else 8)
+
+positions = sim.particles.get_positions()
+positions[:, 1] += 0.5 * sim.dx
+sim.particles.set_positions(positions)
+
+# FRICTION DATA
+mu = np.tan(np.radians(15)) 
+basal_friction = matter.BasalFrictionField()
+if DIM == 3:
+    friction_grid = np.full((2, 2), mu)  # of shape (nz, nx)
+    basal_friction.set(friction_grid, x0=0.0, z0=0.0, d=max(sim.Lx, sim.Lz))
+else:
+    friction_line = np.full(2, mu)
+    basal_friction.set(friction_line, x0=0, d=sim.Lx)
+sim.basal_friction_field = basal_friction
+sim.use_basal_friction_field = True
+
+# sim.add_plate(matter.ObjectPlate(0, matter.PlateType.bottom, matter.BC.SlipFree, mu))
+sim.add_object(matter.ObjectGround(matter.BC.SlipFree, mu, "ground", False, 0))
+
+sim.simulate()
+
+# VERIFICATION:
+# under Coulomb friction: a = g*(sin(theta) - mu*cos(theta))
+final_positions = sim.particles.get_positions()
+final_velocities = sim.particles.get_velocities()
+
+x_sim = final_positions[:, 0].mean()
+vx_sim = final_velocities[:, 0].mean()
+
+T_total = sim.end_frame / sim.fps
+a_analytical = 9.81 * (np.sin(theta) - mu * np.cos(theta))
+vx_analytical = a_analytical * T_total
+dx_analytical = 0.5 * a_analytical * T_total**2
+
+x0_mean = block_size[0] / 2   # initial mean x position (block starts at x=0)
+print("Simulated slide distance:  ", x_sim - x0_mean, "m")
+print("Analytical slide distance: ", dx_analytical, "m")
+print("Simulated velocity:  ", vx_sim, "m/s")
+print("Analytical velocity: ", vx_analytical, "m/s")
+print("Relative velocity difference: ", 100 * abs(vx_sim - vx_analytical) / vx_analytical, "%")