diff --git a/.github/dependabot.yml b/.github/dependabot.yml
index b47ad1889..8717de9c9 100644
--- a/.github/dependabot.yml
+++ b/.github/dependabot.yml
@@ -10,5 +10,5 @@ updates:
schedule:
interval: "daily"
open-pull-requests-limit: 99
- target-branch: "update-v1"
+ target-branch: "master"
versioning-strategy: lockfile-only
diff --git a/.github/workflows/elasticapp.yml b/.github/workflows/elasticapp.yml
new file mode 100644
index 000000000..8fbc24d55
--- /dev/null
+++ b/.github/workflows/elasticapp.yml
@@ -0,0 +1,61 @@
+name: elasticapp (Elastica++ based backend) tests
+
+# trigger run only on changes to the backend folder.
+on:
+ push:
+ paths:
+ - backend/**
+ pull_request:
+ paths:
+ - backend/**
+
+jobs:
+ build:
+ runs-on: ${{ matrix.os }}
+ strategy:
+ matrix:
+ python-version: ["3.11"] #, "3.12"]
+ os: [ubuntu-latest] # , macos-latest]
+ include:
+ - os: ubuntu-latest
+ path: ~/.cache/pip
+ defaults:
+ run:
+ shell: bash
+
+ steps:
+ - uses: actions/checkout@v4
+ - name: Install uv
+ uses: astral-sh/setup-uv@v5
+ with:
+ uv-version: latest
+ - name: Compile OpenMP
+ env:
+ OMP_NUM_THREADS: 2
+ run: |
+ sudo apt-get update; sudo apt-get install -y libomp5 libomp-dev
+ - name: Set up Python ${{ matrix.python-version }}
+ uses: actions/setup-python@v5
+ with:
+ python-version: ${{ matrix.python-version }}
+ - name: Set up cache
+ uses: actions/cache@v5
+ with:
+ path: ${{ matrix.path }}
+ key: ${{ runner.os }}-pip-${{ matrix.python-version }}-${{ hashFiles('pyproject.toml') }}-${{ hashFiles('uv.lock') }}
+ restore-keys: |
+ ${{ runner.os }}-pip-${{ matrix.python-version }}-${{ hashFiles('pyproject.toml') }}-${{ hashFiles('uv.lock') }}
+ ${{ runner.os }}-pip-${{ matrix.python-version }}-${{ hashFiles('pyproject.toml') }}
+
+ - name: Install PyElastica and dependencies
+ run: |
+ make install-dev-deps PYTHON_VERSION=${{ matrix.python-version }}
+ uv cache prune --ci
+
+ - name: Run elasticapp tests
+ run: |
+ source .venv/bin/activate
+ cd backend
+ make clean-build
+ make test
+ # pytest backend/tests/py
diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml
index 989f99642..8acf2e14a 100644
--- a/.github/workflows/main.yml
+++ b/.github/workflows/main.yml
@@ -3,14 +3,6 @@ name: CI
# Controls when the action will run.
on: [push, pull_request]
-# Older settings:
-# Triggers the workflow on push request events for the master branch,
-# and pull request events for all branches.
-#on:
-# push:
-# branches: [ master ]
-# pull_request:
-# branches: [ '**' ]
# A workflow run is made up of one or more jobs that can run sequentially or in parallel
jobs:
@@ -52,7 +44,7 @@ jobs:
run: python -c "import sys; print(sys.version)"
# Set up cache
- name: Set up cache
- uses: actions/cache@v4
+ uses: actions/cache@v5
with:
path: ${{ matrix.path }}
key: ${{ runner.os }}-pip-${{ matrix.python-version }}-${{ hashFiles('pyproject.toml') }}-${{ hashFiles('uv.lock') }}
@@ -62,7 +54,7 @@ jobs:
# Install dependencies
- name: Install dependencies
run: |
- make install-dev-deps
+ make install-dev-deps PYTHON_VERSION=${{ matrix.python-version }}
uv cache prune --ci
source .venv/bin/activate
# Runs a single command using the runners shell
@@ -102,7 +94,7 @@ jobs:
with:
python-version: ${{ matrix.python-version }}
- name: Set up cache
- uses: actions/cache@v3
+ uses: actions/cache@v4
with:
path: ${{ matrix.path }}
key: ${{ runner.os }}-pip-${{ matrix.python-version }}-${{ hashFiles('pyproject.toml') }}-${{ hashFiles('uv.lock') }}
@@ -111,7 +103,7 @@ jobs:
${{ runner.os }}-pip-${{ matrix.python-version }}-${{ hashFiles('pyproject.toml') }}
- name: Install dependencies
run: |
- make install-dev-deps
+ make install-dev-deps PYTHON_VERSION=${{ matrix.python-version }}
uv cache prune --ci
source .venv/bin/activate
- name: Test PyElastica using pytest
@@ -119,6 +111,6 @@ jobs:
run: |
make test_coverage_xml
- name: Upload coverage to Codecov
- uses: codecov/codecov-action@v3
+ uses: codecov/codecov-action@v5
with:
token: ${{ secrets.CODECOV_TOKEN }}
diff --git a/.github/workflows/publish-to-pypi.yml b/.github/workflows/publish-to-pypi.yml
deleted file mode 100644
index 86f70c90c..000000000
--- a/.github/workflows/publish-to-pypi.yml
+++ /dev/null
@@ -1,27 +0,0 @@
-name: publish
-
-on:
- release:
- types: [created]
-
-jobs:
- deploy:
- runs-on: ubuntu-latest
- steps:
- - uses: actions/checkout@v4
- - name: Install uv
- uses: astral-sh/setup-uv@v5
- with:
- uv-version: latest
- - name: Set up Python
- uses: actions/setup-python@v5
- with:
- python-version: "3.x"
- - name: Build
- run: |
- make build
- - name: Publish distribution 📦 to PyPI
- if: startsWith(github.ref, 'refs/tags')
- uses: pypa/gh-action-pypi-publish@release/v1
- with:
- password: ${{ secrets.PYPI_API_TOKEN }}
diff --git a/.gitignore b/.gitignore
index da20c6d3e..cba0be5fe 100644
--- a/.gitignore
+++ b/.gitignore
@@ -7,6 +7,10 @@ __pycache__/
*.so
*.swp
+.vscode
+docs/_gallery
+docs/gen_modules
+docs/sg_execution_times.rst
# Distribution / packaging
.Python
@@ -215,6 +219,7 @@ sample_prog.py
# txt files
*.txt
+!CMakelists.txt
# movie or video file formats
*.mp4
@@ -237,3 +242,6 @@ outcmaes/*
# csv files
*.csv
+
+# ./backend dependencies
+deps
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 43ebcfe21..deeafcc14 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -1,17 +1,20 @@
default_language_version:
python: python3
-default_stages: [commit, push]
+default_stages: [pre-commit, pre-push]
repos:
- repo: https://github.com/pre-commit/pre-commit-hooks
- rev: v2.5.0
+ rev: v6.0.0
hooks:
- id: trailing-whitespace
- id: check-json
- id: check-yaml
+ - id: check-toml
- id: end-of-file-fixer
exclude: LICENSE
+ - id: check-added-large-files
+ - id: mixed-line-ending
- repo: local
hooks:
diff --git a/.readthedocs.yaml b/.readthedocs.yaml
index e0a0cc28c..cb878e965 100644
--- a/.readthedocs.yaml
+++ b/.readthedocs.yaml
@@ -4,6 +4,8 @@ version: 2
# Set the version of Python and other tools you might need
build:
os: ubuntu-24.04
+ apt_packages:
+ - ffmpeg
tools: { python: "3.10" }
jobs:
create_environment:
@@ -14,6 +16,7 @@ build:
install:
- "true"
+
# Build documentation in the docs/ directory with Sphinx
sphinx:
builder: html
diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md
index e4108c526..9358a501b 100644
--- a/CONTRIBUTING.md
+++ b/CONTRIBUTING.md
@@ -121,7 +121,7 @@ Please don't hesitate improving [documentation](https://github.com/GazzolaLab/Py
We also have many related projects in separate repositories that utilize the PyElastica as a core library.
Since the package is often used for research purpose, many on-going projects are typically not published.
-If you are interested in hearing more, please contact one of our the maintainer.
+If you are interested in hearing more, please contact the maintainer.
### Pull requests
diff --git a/Makefile b/Makefile
index 523ab386a..7be633efc 100644
--- a/Makefile
+++ b/Makefile
@@ -2,15 +2,24 @@
PYTHON := python3
PYTHONPATH := `pwd`
AUTOFLAKE_ARGS := -r
+PYTHON_VERSION :=
#* Installation
.PHONY: install
install:
+ifdef PYTHON_VERSION
+ uv sync --python $(PYTHON_VERSION)
+else
uv sync
+endif
.PHONY: install-dev-deps
install-dev-deps:
+ifdef PYTHON_VERSION
+ uv sync --all-groups --all-extras --python $(PYTHON_VERSION)
+else
uv sync --all-groups --all-extras
+endif
.PHONY: build
@@ -25,7 +34,7 @@ pre-commit-install:
.PHONY: black
black:
uv run black --version
- uv run black --config pyproject.toml elastica tests examples
+ uv run black --config pyproject.toml elastica tests examples backend
.PHONY: black-check
black-check:
@@ -35,7 +44,7 @@ black-check:
.PHONY: flake8
flake8:
uv run flake8 --version
- uv run flake8 elastica tests
+ uv run flake8 elastica
.PHONY: autoflake-check
autoflake-check:
@@ -45,7 +54,7 @@ autoflake-check:
.PHONY: autoflake-format
autoflake-format:
uv run autoflake --version
- uv run autoflake --in-place $(AUTOFLAKE_ARGS) elastica tests examples
+ uv run autoflake --in-place $(AUTOFLAKE_ARGS) elastica tests examples backend
.PHONY: format-codestyle
format-codestyle: black autoflake-format
@@ -62,15 +71,15 @@ mypy:
.PHONY: test
test:
- uv run pytest -c pyproject.toml
+ uv run pytest -c pyproject.toml tests
.PHONY: test_coverage
test_coverage:
- NUMBA_DISABLE_JIT=1 uv run pytest --cov=elastica -c pyproject.toml
+ NUMBA_DISABLE_JIT=1 uv run pytest --cov=elastica -c pyproject.toml tests
.PHONY: test_coverage_xml
test_coverage_xml:
- NUMBA_DISABLE_JIT=1 uv run pytest --cov=elastica --cov-report=xml -c pyproject.toml
+ NUMBA_DISABLE_JIT=1 uv run pytest --cov=elastica --cov-report=xml -c pyproject.toml tests
.PHONY: check-codestyle
check-codestyle: black-check flake8 autoflake-check
@@ -103,8 +112,12 @@ pytestcache-remove:
build-remove:
rm -rf build/ dist/
+.PHONY: doc-remove
+doc-remove:
+ rm -rf docs/_build docs/gen_modules/ docs/sg_execution_times.rst docs/_gallery/
+
.PHONY: cleanup
-cleanup: pycache-remove dsstore-remove ipynbcheckpoints-remove pytestcache-remove mypycache-remove build-remove
+cleanup: pycache-remove dsstore-remove ipynbcheckpoints-remove pytestcache-remove mypycache-remove build-remove doc-remove
all: format-codestyle cleanup test
diff --git a/README.md b/README.md
index 8112480b6..837e6d6e1 100644
--- a/README.md
+++ b/README.md
@@ -1,12 +1,12 @@
PyElastica
-[![CI][badge-CI]][link-CI] [![Documentation Status][badge-docs-status]][link-docs-status] [![codecov][badge-codecov]][link-codecov] [![Downloads][badge-pepy-download-count]][link-pepy-download-count] [![DOI][badge-doi]][link-doi] [![Binder][badge-binder]][link-binder] [![Gitter][badge-gitter]][link-gitter]
+[![CI][badge-CI]][link-CI] [![Documentation Status][badge-docs-status]][link-docs-status] [![codecov][badge-codecov]][link-codecov] [![Downloads][badge-pepy-download-count]][link-pepy-download-count] [![DOI][badge-doi]][link-doi] [![Gitter][badge-gitter]][link-gitter]
PyElastica is the python implementation of **Elastica**: an *open-source* project for simulating assemblies of slender, one-dimensional structures using Cosserat Rod theory.
-[![gallery][link-readme-gallary]][link-project-website]
+[![gallery][link-readme-gallery]][link-project-website]
Visit [www.cosseratrods.org][link-project-website] for more information and learn about Elastica and Cosserat rod theory.
@@ -84,13 +84,6 @@ We ask that any publications which use Elastica cite as following:
- [Controlling a CyberOctopus soft arm with muscle-like actuation](https://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=9683318) (UIUC, 2020) (IEEE CDC 2021)
- [Energy shaping control of a CyberOctopus soft arm](https://ieeexplore.ieee.org/document/9304408) (UIUC, 2020) (IEEE CDC 2020)
-## Tutorials
-[![Binder][badge-binder-tutorial]][link-binder]
-
-We have created several Jupyter notebooks and Python scripts to help users get started with PyElastica. The Jupyter notebooks are available on Binder, allowing you to try out some of the tutorials without having to install PyElastica.
-
-We have also included an example script for visualizing PyElastica simulations using POVray. This script is located in the examples folder ([`examples/Visualization`](examples/Visualization)).
-
## Contribution
If you would like to participate, please read our [contribution guideline](CONTRIBUTING.md). Private development branches are moved to [elastica-python](https://github.com/GazzolaLab/elastica-python) repository; access is limited to the core developers, collaborators, and maintainers.
@@ -113,7 +106,7 @@ _Names arranged alphabetically_
[//]: # (Collection of URLs.)
-[link-readme-gallary]: https://github.com/skim0119/PyElastica/blob/assets_logo/assets/alpha_gallery.gif
+[link-readme-gallery]: https://github.com/skim0119/PyElastica/blob/assets_logo/assets/alpha_gallery.gif
[link-project-website]: https://cosseratrods.org
[link-lab-website]: http://mattia-lab.com/
@@ -122,7 +115,6 @@ _Names arranged alphabetically_
[badge-pypi]: https://badge.fury.io/py/pyelastica.svg
[badge-CI]: https://github.com/GazzolaLab/PyElastica/workflows/CI/badge.svg
[badge-docs-status]: https://readthedocs.org/projects/pyelastica/badge/?version=latest
-[badge-binder]: https://mybinder.org/badge_logo.svg
[badge-pepy-download-count]: https://pepy.tech/badge/pyelastica
[badge-codecov]: https://codecov.io/gh/GazzolaLab/PyElastica/branch/master/graph/badge.svg
[badge-gitter]: https://badges.gitter.im/PyElastica/community.svg
@@ -133,7 +125,5 @@ _Names arranged alphabetically_
[link-pepy-download-count]: https://pepy.tech/project/pyelastica
[link-codecov]: https://codecov.io/gh/GazzolaLab/PyElastica
-[badge-binder-tutorial]: https://img.shields.io/badge/Launch-PyElastica%20Tutorials-579ACA.svg?logo=data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAFkAAABZCAMAAABi1XidAAAB8lBMVEX///9XmsrmZYH1olJXmsr1olJXmsrmZYH1olJXmsr1olJXmsrmZYH1olL1olJXmsr1olJXmsrmZYH1olL1olJXmsrmZYH1olJXmsr1olL1olJXmsrmZYH1olL1olJXmsrmZYH1olL1olL0nFf1olJXmsrmZYH1olJXmsq8dZb1olJXmsrmZYH1olJXmspXmspXmsr1olL1olJXmsrmZYH1olJXmsr1olL1olJXmsrmZYH1olL1olLeaIVXmsrmZYH1olL1olL1olJXmsrmZYH1olLna31Xmsr1olJXmsr1olJXmsrmZYH1olLqoVr1olJXmsr1olJXmsrmZYH1olL1olKkfaPobXvviGabgadXmsqThKuofKHmZ4Dobnr1olJXmsr1olJXmspXmsr1olJXmsrfZ4TuhWn1olL1olJXmsqBi7X1olJXmspZmslbmMhbmsdemsVfl8ZgmsNim8Jpk8F0m7R4m7F5nLB6jbh7jbiDirOEibOGnKaMhq+PnaCVg6qWg6qegKaff6WhnpKofKGtnomxeZy3noG6dZi+n3vCcpPDcpPGn3bLb4/Mb47UbIrVa4rYoGjdaIbeaIXhoWHmZYHobXvpcHjqdHXreHLroVrsfG/uhGnuh2bwj2Hxk17yl1vzmljzm1j0nlX1olL3AJXWAAAAbXRSTlMAEBAQHx8gICAuLjAwMDw9PUBAQEpQUFBXV1hgYGBkcHBwcXl8gICAgoiIkJCQlJicnJ2goKCmqK+wsLC4usDAwMjP0NDQ1NbW3Nzg4ODi5+3v8PDw8/T09PX29vb39/f5+fr7+/z8/Pz9/v7+zczCxgAABC5JREFUeAHN1ul3k0UUBvCb1CTVpmpaitAGSLSpSuKCLWpbTKNJFGlcSMAFF63iUmRccNG6gLbuxkXU66JAUef/9LSpmXnyLr3T5AO/rzl5zj137p136BISy44fKJXuGN/d19PUfYeO67Znqtf2KH33Id1psXoFdW30sPZ1sMvs2D060AHqws4FHeJojLZqnw53cmfvg+XR8mC0OEjuxrXEkX5ydeVJLVIlV0e10PXk5k7dYeHu7Cj1j+49uKg7uLU61tGLw1lq27ugQYlclHC4bgv7VQ+TAyj5Zc/UjsPvs1sd5cWryWObtvWT2EPa4rtnWW3JkpjggEpbOsPr7F7EyNewtpBIslA7p43HCsnwooXTEc3UmPmCNn5lrqTJxy6nRmcavGZVt/3Da2pD5NHvsOHJCrdc1G2r3DITpU7yic7w/7Rxnjc0kt5GC4djiv2Sz3Fb2iEZg41/ddsFDoyuYrIkmFehz0HR2thPgQqMyQYb2OtB0WxsZ3BeG3+wpRb1vzl2UYBog8FfGhttFKjtAclnZYrRo9ryG9uG/FZQU4AEg8ZE9LjGMzTmqKXPLnlWVnIlQQTvxJf8ip7VgjZjyVPrjw1te5otM7RmP7xm+sK2Gv9I8Gi++BRbEkR9EBw8zRUcKxwp73xkaLiqQb+kGduJTNHG72zcW9LoJgqQxpP3/Tj//c3yB0tqzaml05/+orHLksVO+95kX7/7qgJvnjlrfr2Ggsyx0eoy9uPzN5SPd86aXggOsEKW2Prz7du3VID3/tzs/sSRs2w7ovVHKtjrX2pd7ZMlTxAYfBAL9jiDwfLkq55Tm7ifhMlTGPyCAs7RFRhn47JnlcB9RM5T97ASuZXIcVNuUDIndpDbdsfrqsOppeXl5Y+XVKdjFCTh+zGaVuj0d9zy05PPK3QzBamxdwtTCrzyg/2Rvf2EstUjordGwa/kx9mSJLr8mLLtCW8HHGJc2R5hS219IiF6PnTusOqcMl57gm0Z8kanKMAQg0qSyuZfn7zItsbGyO9QlnxY0eCuD1XL2ys/MsrQhltE7Ug0uFOzufJFE2PxBo/YAx8XPPdDwWN0MrDRYIZF0mSMKCNHgaIVFoBbNoLJ7tEQDKxGF0kcLQimojCZopv0OkNOyWCCg9XMVAi7ARJzQdM2QUh0gmBozjc3Skg6dSBRqDGYSUOu66Zg+I2fNZs/M3/f/Grl/XnyF1Gw3VKCez0PN5IUfFLqvgUN4C0qNqYs5YhPL+aVZYDE4IpUk57oSFnJm4FyCqqOE0jhY2SMyLFoo56zyo6becOS5UVDdj7Vih0zp+tcMhwRpBeLyqtIjlJKAIZSbI8SGSF3k0pA3mR5tHuwPFoa7N7reoq2bqCsAk1HqCu5uvI1n6JuRXI+S1Mco54YmYTwcn6Aeic+kssXi8XpXC4V3t7/ADuTNKaQJdScAAAAAElFTkSuQmCC
-[link-binder]: https://mybinder.org/v2/gh/GazzolaLab/PyElastica/master?filepath=examples%2FBinder%2F0_PyElastica_Tutorials_Overview.ipynb
[link-gitter]: https://gitter.im/PyElastica/community?utm_source=badge&utm_medium=badge&utm_campaign=pr-badge&utm_content=badge
[link-doi]: https://zenodo.org/badge/latestdoi/254172891
diff --git a/backend/.clang-format b/backend/.clang-format
new file mode 100644
index 000000000..48d4ec71d
--- /dev/null
+++ b/backend/.clang-format
@@ -0,0 +1,17 @@
+---
+# We'll use defaults from the Google style.
+# See http://clang.llvm.org/docs/ClangFormat.html for help.
+Language: Cpp
+BasedOnStyle: Google
+AllowShortIfStatementsOnASingleLine: false
+AllowShortLoopsOnASingleLine: false
+PointerAlignment: Left
+DerivePointerAlignment: false
+FixNamespaceComments: true
+IncludeCategories:
+ - Regex: "^<.*"
+ Priority: 1
+ - Regex: ".*"
+ Priority: 2
+NamespaceIndentation: All
+SortIncludes: false
diff --git a/backend/.gitignore b/backend/.gitignore
new file mode 100644
index 000000000..4689dcf2c
--- /dev/null
+++ b/backend/.gitignore
@@ -0,0 +1,12 @@
+CMakeLists.txt.user
+CMakeCache.txt
+CMakeFiles
+CMakeScripts
+Testing
+Makefile
+cmake_install.cmake
+install_manifest.txt
+compile_commands.json
+CTestTestfile.cmake
+_deps
+
diff --git a/backend/CMakeLists.txt b/backend/CMakeLists.txt
new file mode 100644
index 000000000..be937721d
--- /dev/null
+++ b/backend/CMakeLists.txt
@@ -0,0 +1,290 @@
+cmake_minimum_required(VERSION 3.20)
+project(elasticapp VERSION 0.0.4 LANGUAGES CXX)
+
+# Set C++ standard
+set(CMAKE_CXX_STANDARD 20)
+set(CMAKE_CXX_STANDARD_REQUIRED ON)
+set(CMAKE_CXX_EXTENSIONS OFF)
+
+# Compiler flags
+add_compile_options(-Wno-unused -fdiagnostics-color=always -Wall -Wextra -pedantic -flto)
+
+# Find Python and pybind11 (temp)
+# pybind11 will find Python and NumPy automatically
+# Set Python interpreter to ../.venv
+set(Python3_ROOT_DIR "${CMAKE_CURRENT_SOURCE_DIR}/../.venv")
+set(PYTHON_EXECUTABLE "${Python3_ROOT_DIR}/bin/python")
+
+# Find pybind11 and ensure it uses the Python interpreter from our venv
+set(PYBIND11_FINDPYTHON "ON")
+set(PYBIND11_PYTHON_EXECUTABLE "${PYTHON_EXECUTABLE}" CACHE FILEPATH "" FORCE)
+find_package(pybind11 REQUIRED)
+
+# ------------------------------------------------------------ #
+# External dependency management #
+# Include FetchContent for Catch2 and Eigen3
+include(FetchContent)
+
+# Save current BUILD_TESTING state and disable it for dependencies
+set(BUILD_TESTING_SAVED ${BUILD_TESTING})
+set(BUILD_TESTING OFF)
+
+# Fetch Eigen3 for numerical computations
+FetchContent_Declare(
+ Eigen3
+ GIT_REPOSITORY https://gitlab.com/libeigen/eigen.git
+ GIT_TAG 3.4.0
+)
+# Disable Eigen3 tests and documentation
+set(EIGEN_BUILD_DOC OFF)
+# Set Eigen's default storage order to row-major
+# This ensures Eigen's default behavior matches our explicit MatrixType template parameter
+set(EIGEN_DEFAULT_TO_ROW_MAJOR ON CACHE BOOL "Use row-major as default matrix storage order" FORCE)
+FetchContent_MakeAvailable(Eigen3)
+
+# Restore BUILD_TESTING state for our tests
+set(BUILD_TESTING ${BUILD_TESTING_SAVED})
+
+# ------------------------------------------------------------ #
+
+# Define directories #
+set(TEST_DIR "${CMAKE_CURRENT_SOURCE_DIR}/tests")
+set(PYTHON_PACKAGE_DIR "${CMAKE_CURRENT_SOURCE_DIR}/src/py")
+set(SRC_DIR "${CMAKE_CURRENT_SOURCE_DIR}/src")
+
+option(ELASTICAPP_VECTORIZATION_REPORTS "Enable compiler vectorization reports" OFF)
+option(ELASTICAPP_BUILD_TESTS "Build C++ unit tests" OFF)
+option(ELASTICAPP_GIL_RELEASE "Release Python GIL around compute-heavy bindings" ON )
+
+# Global CMake option to enable/disable for unittests
+if (DEFINED ENV{ELASTICAPP_BUILD_TESTS})
+ string(TOUPPER "$ENV{ELASTICAPP_BUILD_TESTS}" ELASTICAPP_BUILD_TESTS_ENV_VALUE)
+ if ("${ELASTICAPP_BUILD_TESTS_ENV_VALUE}" MATCHES "^(ON|TRUE|1)$")
+ set(ELASTICAPP_BUILD_TESTS ON)
+ message(STATUS "ELASTICAPP_BUILD_TESTS set to TRUE from environment variable.")
+ elseif ("${ELASTICAPP_BUILD_TESTS_ENV_VALUE}" MATCHES "^(OFF|FALSE|0)$")
+ set(ELASTICAPP_BUILD_TESTS OFF)
+ message(STATUS "ELASTICAPP_BUILD_TESTS set to OFF from environment variable.")
+ else()
+ message(WARNING "Environment variable ELASTICAPP_BUILD_TESTS has an unrecognized value. It will not be used to set the CMake option.")
+ endif()
+endif()
+
+# Option to set number of threads per component in update_dynamics parallel sections
+# This controls how many threads are used for each spatial component (x, y, z) in the parallel sections
+# Default is 3 (one thread per component), but can be adjusted for different hardware configurations
+set(ELASTICAPP_COMPONENT_THREADS "2" CACHE STRING "Number of threads per component in update_dynamics parallel sections (default: 3)")
+mark_as_advanced(ELASTICAPP_COMPONENT_THREADS)
+
+# Configure threading support (OpenMP) - check early but apply per-target
+# On macOS with Homebrew, help CMake find libomp
+# This is needed for both Makefile builds and scikit-build-core builds
+if(APPLE)
+ # Try to find Homebrew-installed libomp
+ execute_process(
+ COMMAND brew --prefix libomp
+ OUTPUT_VARIABLE HOMEBREW_LIBOMP_PREFIX
+ OUTPUT_STRIP_TRAILING_WHITESPACE
+ ERROR_QUIET
+ )
+ if(HOMEBREW_LIBOMP_PREFIX)
+ message(STATUS "Found Homebrew libomp at: ${HOMEBREW_LIBOMP_PREFIX}")
+ # Set OpenMP variables for clang on macOS
+ if(CMAKE_CXX_COMPILER_ID MATCHES "Clang")
+ set(OpenMP_CXX_FLAGS "-Xpreprocessor -fopenmp -I${HOMEBREW_LIBOMP_PREFIX}/include")
+ set(OpenMP_CXX_LIB_NAMES "omp")
+ set(OpenMP_omp_LIBRARY "${HOMEBREW_LIBOMP_PREFIX}/lib/libomp.dylib")
+ set(OpenMP_CXX_LIBRARIES "${HOMEBREW_LIBOMP_PREFIX}/lib/libomp.dylib")
+ message(STATUS "Configured OpenMP for clang with Homebrew libomp")
+ endif()
+ endif()
+endif()
+
+find_package(OpenMP COMPONENTS CXX REQUIRED)
+
+# Helper function to apply optimization flags to a target
+# This function must be defined before targets are created
+function(apply_optimization_flags target_name)
+ if(CMAKE_CXX_COMPILER_ID MATCHES "GNU|Clang")
+ target_compile_options(${target_name} PRIVATE -O3 -march=native)
+ target_compile_options(${target_name} PRIVATE -ffunction-sections -fdata-sections)
+ elseif(CMAKE_CXX_COMPILER_ID MATCHES "MSVC")
+ target_compile_options(${target_name} PRIVATE /O2)
+ endif()
+ # Add NDEBUG flag to disable assertions in release builds
+ target_compile_definitions(${target_name} PRIVATE NDEBUG)
+
+ # Enable Eigen vectorization
+ target_compile_definitions(${target_name} PRIVATE EIGEN_VECTORIZE)
+
+ # Add architecture-specific optimization flags
+ if(CMAKE_CXX_COMPILER_ID MATCHES "GNU|Clang")
+ # AVX2 and FMA are only available on x86_64, not ARM
+ if(CMAKE_SYSTEM_PROCESSOR MATCHES "x86_64|AMD64")
+ target_compile_options(${target_name} PRIVATE -mavx2 -mfma)
+ endif()
+
+ # Add vectorization report flags if requested
+ if(ELASTICAPP_VECTORIZATION_REPORTS)
+ if(CMAKE_CXX_COMPILER_ID MATCHES "Clang")
+ # Clang vectorization reports
+ target_compile_options(${target_name} PRIVATE
+ -Rpass=loop-vectorize
+ -Rpass-missed=loop-vectorize
+ -Rpass-analysis=loop-vectorize
+ )
+ message(STATUS "Clang vectorization reports enabled for ${target_name}")
+ elseif(CMAKE_CXX_COMPILER_ID MATCHES "GNU")
+ # GCC vectorization reports
+ target_compile_options(${target_name} PRIVATE
+ -fopt-info-vec
+ -fopt-info-vec-missed
+ -fopt-info-vec-all
+ )
+ message(STATUS "GCC vectorization reports enabled for ${target_name}")
+ endif()
+ endif()
+ elseif(CMAKE_CXX_COMPILER_ID MATCHES "MSVC")
+ # MSVC flags for AVX2
+ target_compile_options(${target_name} PRIVATE /arch:AVX2)
+ endif()
+
+ if(OpenMP_CXX_FOUND)
+ target_link_libraries(${target_name} PRIVATE OpenMP::OpenMP_CXX)
+ message(STATUS "Threading support enabled for ${target_name}")
+ endif()
+endfunction()
+
+# Extension module: _memory_block (Python module name)
+pybind11_add_module(_memory_block
+ "${SRC_DIR}/_api.cpp"
+)
+target_compile_definitions(_memory_block PRIVATE VERSION_INFO=${PROJECT_VERSION})
+target_include_directories(_memory_block PRIVATE ${SRC_DIR})
+target_link_libraries(_memory_block PRIVATE Eigen3::Eigen)
+
+if(APPLE)
+ target_link_options(_memory_block PRIVATE -Wl,-dead_strip)
+elseif(CMAKE_CXX_COMPILER_ID MATCHES "GNU")
+ target_link_options(_memory_block PRIVATE -Wl,--gc-sections)
+elseif(CMAKE_CXX_COMPILER_ID MATCHES "Clang")
+ target_link_options(_memory_block PRIVATE -Wl,--gc-sections)
+endif()
+
+
+# Apply component threads option to _memory_block (only if explicitly set and non-empty)
+# If not set, the code will use OpenMP's default thread count
+if(ELASTICAPP_COMPONENT_THREADS AND NOT ELASTICAPP_COMPONENT_THREADS STREQUAL "" AND NOT ELASTICAPP_COMPONENT_THREADS STREQUAL "0")
+ target_compile_definitions(_memory_block PRIVATE ELASTICAPP_COMPONENT_THREADS=${ELASTICAPP_COMPONENT_THREADS})
+endif()
+
+if(ELASTICAPP_GIL_RELEASE)
+ target_compile_definitions(_memory_block PRIVATE ELASTICAPP_GIL_RELEASE)
+endif()
+
+# Apply optimization flags to _memory_block
+apply_optimization_flags(_memory_block)
+
+# Install _memory_block module
+install(TARGETS _memory_block
+ LIBRARY DESTINATION "elasticapp"
+ COMPONENT python)
+
+# Extension module: version (Python module name)
+pybind11_add_module(version
+ "${SRC_DIR}/_version.cpp"
+)
+target_compile_definitions(version PRIVATE VERSION_INFO=${PROJECT_VERSION})
+target_include_directories(version PRIVATE ${SRC_DIR})
+
+# Install version module
+install(TARGETS version
+ LIBRARY DESTINATION "elasticapp"
+ COMPONENT python)
+
+# Extension module: _collision (Python module name)
+pybind11_add_module(_collision
+ "${SRC_DIR}/environment/_api.cpp"
+ "${SRC_DIR}/environment/collision/course_detection/hash_grid.cpp"
+ "${SRC_DIR}/environment/collision/fine_detection/max_contacts.cpp"
+ "${SRC_DIR}/environment/collision/batching/union_find.cpp"
+ "${SRC_DIR}/environment/collision/physics/linear_spring_dashpot.cpp"
+ "${SRC_DIR}/environment/collision/physics/no_interaction.cpp"
+ "${SRC_DIR}/math/node_element_mapping.cpp"
+)
+target_compile_definitions(_collision PRIVATE VERSION_INFO=${PROJECT_VERSION})
+target_include_directories(_collision PRIVATE ${SRC_DIR})
+target_link_libraries(_collision PRIVATE Eigen3::Eigen)
+
+if(APPLE)
+ target_link_options(_collision PRIVATE -Wl,-dead_strip)
+elseif(CMAKE_CXX_COMPILER_ID MATCHES "GNU")
+ target_link_options(_collision PRIVATE -Wl,--gc-sections)
+elseif(CMAKE_CXX_COMPILER_ID MATCHES "Clang")
+ target_link_options(_collision PRIVATE -Wl,--gc-sections)
+endif()
+
+if(ELASTICAPP_GIL_RELEASE)
+ target_compile_definitions(_collision PRIVATE ELASTICAPP_GIL_RELEASE)
+endif()
+
+# Apply optimization flags to _collision
+apply_optimization_flags(_collision)
+
+# Install _collision module
+install(TARGETS _collision
+ LIBRARY DESTINATION "elasticapp"
+ COMPONENT python)
+
+# Install Python package files
+install(FILES
+ "${PYTHON_PACKAGE_DIR}/__init__.py"
+ "${PYTHON_PACKAGE_DIR}/typing_alias.py"
+ "${PYTHON_PACKAGE_DIR}/memory_block_rod.py"
+ "${PYTHON_PACKAGE_DIR}/collision_physics.py"
+ "${PYTHON_PACKAGE_DIR}/module_collision.py"
+ DESTINATION "elasticapp"
+ COMPONENT python)
+
+# C++ unit tests with Catch2
+if(ELASTICAPP_BUILD_TESTS)
+ enable_testing()
+ include(CTest)
+
+ # Fetch Catch2 for C++ unit testing
+ FetchContent_Declare(
+ Catch2
+ GIT_REPOSITORY https://github.com/catchorg/Catch2.git
+ GIT_TAG v3.5.4
+ )
+ set(CATCH_BUILD_TESTING OFF CACHE BOOL "" FORCE)
+ set(CATCH_BUILD_EXAMPLES OFF CACHE BOOL "" FORCE)
+ FetchContent_MakeAvailable(Catch2)
+
+ # Test executable
+ add_executable(cpp_tests
+ "${TEST_DIR}/cpp/test_version.cpp"
+ "${TEST_DIR}/cpp/test_block.cpp"
+ "${TEST_DIR}/cpp/test_traits.cpp"
+ "${TEST_DIR}/cpp/test_system.cpp"
+ "${TEST_DIR}/cpp/test_system_variable_validation.cpp"
+ "${TEST_DIR}/cpp/test_block_view.cpp"
+ "${TEST_DIR}/cpp/test_block_operations.cpp"
+ "${TEST_DIR}/cpp/collision/test_collision_system.cpp"
+ )
+
+ target_compile_definitions(cpp_tests PRIVATE VERSION_INFO=${PROJECT_VERSION})
+ target_include_directories(cpp_tests PRIVATE ${SRC_DIR} ${TEST_DIR}/cpp/mock)
+ target_link_libraries(cpp_tests PRIVATE Catch2::Catch2WithMain Eigen3::Eigen)
+
+ # Apply optimization flags to cpp_tests
+ apply_optimization_flags(cpp_tests)
+
+ # Use Catch2's test discovery to automatically register tests with CTest
+ # Add Catch2 extras directory to module path (FetchContent uses uppercase)
+ list(APPEND CMAKE_MODULE_PATH ${Catch2_SOURCE_DIR}/extras)
+ include(Catch)
+ catch_discover_tests(cpp_tests)
+else()
+ message(STATUS "C++ tests disabled (ELASTICAPP_BUILD_TESTS=OFF)")
+endif()
diff --git a/backend/Makefile b/backend/Makefile
new file mode 100644
index 000000000..9e94afa7d
--- /dev/null
+++ b/backend/Makefile
@@ -0,0 +1,62 @@
+#* Variables
+PYTHON := python3
+PYTHONPATH := `pwd`/..
+PYTHON_VERSION :=
+
+#* Installation
+.PHONY: install
+install:
+ uv pip install "." -v
+
+#* C++ Build (direct CMake, no pip)
+.PHONY: build-cpp
+build-cpp:
+ @echo "Building C++ code with CMake..." && \
+ VENV_PYTHON=$$(which python) && \
+ PYBIND11_DIR=$$(python -c "import pybind11; import os; print(os.path.join(os.path.dirname(pybind11.__file__), 'share', 'cmake', 'pybind11'))") && \
+ mkdir -p build && \
+ cd build && \
+ LIBOMP_PREFIX=$$(brew --prefix libomp 2>/dev/null || echo "") && \
+ LLVM_PREFIX=$$(brew --prefix llvm 2>/dev/null || echo "") && \
+ cmake .. \
+ -DCMAKE_BUILD_TYPE=Release \
+ -DCMAKE_C_COMPILER=$$LLVM_PREFIX/bin/clang \
+ -DCMAKE_CXX_COMPILER=$$LLVM_PREFIX/bin/clang++ \
+ -DPython3_ROOT_DIR=$$(cd .. && pwd)/../.venv \
+ -DPython3_EXECUTABLE="$$VENV_PYTHON" \
+ -DPYTHON_EXECUTABLE="$$VENV_PYTHON" \
+ -DPython_EXECUTABLE="$$VENV_PYTHON" \
+ -Dpybind11_DIR="$$PYBIND11_DIR" \
+ -DOpenMP_CXX_FLAGS="-Xpreprocessor -fopenmp -I$$LIBOMP_PREFIX/include" \
+ -DOpenMP_CXX_LIB_NAMES="omp" \
+ -DOpenMP_omp_LIBRARY="$$LIBOMP_PREFIX/lib/libomp.dylib" \
+ -DOpenMP_CXX_LIBRARIES="$$LIBOMP_PREFIX/lib/libomp.dylib" \
+ -DELASTICAPP_VECTORIZATION_REPORTS=ON \
+ -DELASTICAPP_SAVE_ASSEMBLY=ON \
+ -DELASTICAPP_BUILD_TESTS=OFF; \
+ cmake --build . --parallel 2>&1 | tee ../vectorization_report.txt
+
+#* C++ Tests (build and run, no pip)
+.PHONY: test-cpp-direct
+test-cpp-direct: build-cpp
+ @echo "Running C++ tests..."
+ @cd build && ./cpp_tests
+
+#* Testing
+.PHONY: test
+test:
+ @echo "Running C++ and Python tests..."
+ ELASTICAPP_BUILD_TESTS=ON make install && cd build && ./cpp_tests
+ pytest tests/py -v
+
+#* Cleaning
+.PHONY: clean
+clean: clean-build
+ rm -rf dist/ *.egg-info
+ find . -type d -name __pycache__ -exec rm -r {} + 2>/dev/null || true
+ find . -type f -name "*.pyc" -delete
+ find . -type f -name "*.pyo" -delete
+
+.PHONY: clean-build
+clean-build:
+ rm -rf build/
diff --git a/backend/README.md b/backend/README.md
new file mode 100644
index 000000000..a8828bc6e
--- /dev/null
+++ b/backend/README.md
@@ -0,0 +1,61 @@
+# C++ for PyElastica
+
+This directory includes basic C++ backend for experimental purposes.
+The full functionality is not yet moved to this directory.
+The purpose of this update is to provide a way to implement functionality in C++, such as multithreading and more controlled vectorization.
+For large-number of rods and complex environment, computation with contact and self-interactions requires higher control over parallelization.
+
+> Note: manual tuning of vectorization and threading is required for optimal performance, which may vary depending on the hardware.
+
+> Note: expected performance improvement depends on the vector size. It is expected to be faster for large rod systems, mostly due to the overhead created by the python-binding and interpreter. In most of the case with less than 1e5 elements, the performance of `numba` could be better.
+
+## Usage
+
+```python
+import elastica as ea
+import elasticapp as epp
+
+class SystemSimulator(
+ ea.BaseSystemCollection,
+ ...
+):
+ pass
+
+simulator = SystemSimulator()
+# Replace C++ block module for CosseratRod
+simulator.enable_block_supports(ea.CosseratRod, epp.MemoryBlockCosseratRod)
+```
+
+## Installation
+
+From the `backend` directory, run:
+
+```bash
+cd backend
+make install # or run pip install "."
+```
+> Make sure you install the package from _PyElastica source tree_.
+
+## Testing
+
+Test includes both `cpp` and `python` testing.
+
+```bash
+make test
+```
+
+## File structure
+
+- All cpp files are in `src` directory. `src/py` directory includes python bindings for C++ classes.
+- All test files are in `tests` directory. `cpp` tests are in `tests/cpp` directory, and `python` tests are in `tests/py` directory.
+- All benchmark files are in `benchmarking` directory.
+
+
+## Contributed By
+
+- Tejaswin Parthasarathy (Teja)
+- [Seung Hyun Kim](https://github.com/skim0119)
+- Ankith Pai
+- [Yashraj Bhosale](https://github.com/bhosale2)
+- Arman Tekinalp
+- Songyuan Cui
diff --git a/backend/benchmarking/PDE_benchmark.ipynb b/backend/benchmarking/PDE_benchmark.ipynb
new file mode 100644
index 000000000..f76de55b6
--- /dev/null
+++ b/backend/benchmarking/PDE_benchmark.ipynb
@@ -0,0 +1,231 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "id": "ea3e0b99",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from typing import Any\n",
+ "import timeit\n",
+ "import numpy as np\n",
+ "import matplotlib.pyplot as plt\n",
+ "import time\n",
+ "from tqdm import tqdm\n",
+ "import os\n",
+ "# os.environ[\"NUMBA_DISABLE_JIT\"] = \"1\"\n",
+ "\n",
+ "import elastica as ea\n",
+ "from elasticapp import MemoryBlockCosseratRod\n",
+ "\n",
+ "# Try to import set_num_threads if threading is enabled\n",
+ "try:\n",
+ " from elasticapp._memory_block import set_num_threads, get_max_threads\n",
+ " THREADING_AVAILABLE = True\n",
+ "except ImportError:\n",
+ " THREADING_AVAILABLE = False\n",
+ " raise RuntimeError(\"Threading not available (ELASTICAPP_USE_THREADING not enabled)\")\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "id": "bc221bc9",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "100%|██████████| 50000/50000 [00:24<00:00, 2016.30it/s]\n"
+ ]
+ }
+ ],
+ "source": [
+ "# Create 50000 rods with 6 elements each, iterate 100 operations\n",
+ "N = 100\n",
+ "n_rods = 50000\n",
+ "n_elems_per_rod = 6\n",
+ "\n",
+ "# Create rods\n",
+ "rods = [\n",
+ " ea.CosseratRod.straight_rod(\n",
+ " n_elements=n_elems_per_rod,\n",
+ " start=np.zeros(3),\n",
+ " direction=np.array([0.0, 0.0, 1.0]),\n",
+ " normal=np.array([1.0, 0.0, 0.0]),\n",
+ " base_length=1.0,\n",
+ " base_radius=0.01,\n",
+ " density=3000,\n",
+ " youngs_modulus=1e6,\n",
+ " )\n",
+ " for _ in tqdm(range(n_rods))\n",
+ "]\n",
+ "\n",
+ "def jitter_rod(rod):\n",
+ " rng = np.random.default_rng(42)\n",
+ " rod.external_forces[:] = rng.uniform(-1e-6, 1e-6, rod.external_forces.shape)\n",
+ " rod.external_torques[:] = rng.uniform(-1e-6, 1e-6, rod.external_torques.shape)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "id": "35e3d633",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Time per call: 42.549 ms (compute_internal_forces_and_torques - inverse rotation, equations)\n",
+ "Time per call: 5.615 ms (update_accelerations - block addition)\n",
+ "Time per call: 15.882 ms (update_kinematics - rodrigues)\n",
+ "Time per call: 0.974 ms (update_dynamics)\n"
+ ]
+ }
+ ],
+ "source": [
+ "block_rod = ea.MemoryBlockCosseratRod(rods, range(n_rods))\n",
+ "jitter_rod(block_rod)\n",
+ "\n",
+ "py_result = {}\n",
+ "\n",
+ "for _ in range(10): # Pre-run\n",
+ " block_rod.compute_internal_forces_and_torques(0.0)\n",
+ "T = timeit.timeit(lambda: block_rod.compute_internal_forces_and_torques(0.0), number=N)\n",
+ "py_result['compute_internal_forces_and_torques'] = T/N*1000\n",
+ "print(f\"Time per call: {T/N*1000:.3f} ms (compute_internal_forces_and_torques - inverse rotation, equations)\")\n",
+ "for _ in range(10): # Pre-run\n",
+ " block_rod.update_accelerations(0.0)\n",
+ "T = timeit.timeit(lambda: block_rod.update_accelerations(0.0), number=N)\n",
+ "py_result['update_accelerations'] = T/N*1000\n",
+ "print(f\"Time per call: {T/N*1000:.3f} ms (update_accelerations - block addition)\")\n",
+ "for _ in range(10): # Pre-run\n",
+ " block_rod.update_kinematics(0.0, 1e-4)\n",
+ "T = timeit.timeit(lambda: block_rod.update_kinematics(0.0, 1e-4), number=N)\n",
+ "py_result['update_kinematics'] = T/N*1000\n",
+ "print(f\"Time per call: {T/N*1000:.3f} ms (update_kinematics - rodrigues)\")\n",
+ "for _ in range(10): # Pre-run\n",
+ " block_rod.update_dynamics(0.0, 1e-4)\n",
+ "T = timeit.timeit(lambda: block_rod.update_dynamics(0.0, 1e-4), number=N)\n",
+ "py_result['update_dynamics'] = T/N*1000\n",
+ "print(f\"Time per call: {T/N*1000:.3f} ms (update_dynamics)\")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "id": "31e6d7a1",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "100%|██████████| 50000/50000 [00:24<00:00, 2057.62it/s]\n"
+ ]
+ },
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Time per call: 21.076 ms (compute_internal_forces_and_torques - inverse rotation, equations)\n",
+ "Time per call: 3.041 ms (update_accelerations - block addition)\n",
+ "Time per call: 1.601 ms (update_kinematics - rodrigues)\n",
+ "Time per call: 0.756 ms (update_dynamics)\n",
+ "------------------------------------------------------------\n"
+ ]
+ }
+ ],
+ "source": [
+ "assert THREADING_AVAILABLE\n",
+ "\n",
+ "\n",
+ "# Benchmark parameters\n",
+ "# thread_counts = [1, 2, 4, 8, 16]\n",
+ "\n",
+ "# Storage for results\n",
+ "cpp_results_for_threads = []\n",
+ "\n",
+ "# Run benchmark for each thread count\n",
+ "# for num_threads in thread_counts:\n",
+ " # set_num_threads(num_threads)\n",
+ " # print(f\"Threading available: {num_threads}\")\n",
+ "\n",
+ "rods = [\n",
+ " ea.CosseratRod.straight_rod(\n",
+ " n_elements=n_elems_per_rod,\n",
+ " start=np.zeros(3),\n",
+ " direction=np.array([0.0, 0.0, 1.0]),\n",
+ " normal=np.array([1.0, 0.0, 0.0]),\n",
+ " base_length=1.0,\n",
+ " base_radius=0.01,\n",
+ " density=3000,\n",
+ " youngs_modulus=1e6,\n",
+ " )\n",
+ " for _ in tqdm(range(n_rods))\n",
+ "]\n",
+ "\n",
+ "block_rod = MemoryBlockCosseratRod(rods, range(n_rods))\n",
+ "jitter_rod(block_rod)\n",
+ "\n",
+ "cpp_result = {}\n",
+ "cpp_results_for_threads.append(cpp_result)\n",
+ "\n",
+ "for _ in range(10): # Pre-run\n",
+ " block_rod.compute_internal_forces_and_torques(0.0)\n",
+ "T = timeit.timeit(lambda: block_rod.compute_internal_forces_and_torques(0.0), number=N)\n",
+ "cpp_result['compute_internal_forces_and_torques'] = T/N*1000\n",
+ "print(f\"Time per call: {T/N*1000:.3f} ms (compute_internal_forces_and_torques - inverse rotation, equations)\")\n",
+ "for _ in range(10): # Pre-run\n",
+ " block_rod.update_accelerations(0.0)\n",
+ "T = timeit.timeit(lambda: block_rod.update_accelerations(0.0), number=N)\n",
+ "cpp_result['update_accelerations'] = T/N*1000\n",
+ "print(f\"Time per call: {T/N*1000:.3f} ms (update_accelerations - block addition)\")\n",
+ "for _ in range(10): # Pre-run\n",
+ " block_rod.update_kinematics(0.0, 1e-4)\n",
+ "T = timeit.timeit(lambda: block_rod.update_kinematics(0.0, 1e-4), number=N)\n",
+ "cpp_result['update_kinematics'] = T/N*1000\n",
+ "print(f\"Time per call: {T/N*1000:.3f} ms (update_kinematics - rodrigues)\")\n",
+ "for _ in range(10): # Pre-run\n",
+ " block_rod.update_dynamics(0.0, 1e-4)\n",
+ "T = timeit.timeit(lambda: block_rod.update_dynamics(0.0, 1e-4), number=N)\n",
+ "cpp_result['update_dynamics'] = T/N*1000\n",
+ "print(f\"Time per call: {T/N*1000:.3f} ms (update_dynamics)\")\n",
+ "\n",
+ "print(\"-\"*60)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "bb8950dc",
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": ".venv",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.11.14"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/backend/benchmarking/memory_block_integrity_check.py b/backend/benchmarking/memory_block_integrity_check.py
new file mode 100644
index 000000000..98e4e6964
--- /dev/null
+++ b/backend/benchmarking/memory_block_integrity_check.py
@@ -0,0 +1,111 @@
+"""
+Example script to benchmark C++ operations with timing done entirely in C++.
+
+This script demonstrates how to use the benchmark_cpp function to measure
+performance without Python loop or pybind11 call overhead.
+"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+import elastica as epy
+import elasticapp as epp
+
+# %%
+# Create test rods
+n_rods = 50
+n_elems_per_rod = 200
+
+print(f"Creating {n_rods} rods with {n_elems_per_rod} elements each...")
+keys = list(epp.memory_block_rod.PY2CPP_VARNAMES.keys())
+rods_cpp = [
+ epy.CosseratRod.straight_rod(
+ n_elements=n_elems_per_rod,
+ start=np.zeros(3),
+ direction=np.array([0.0, 0.0, 1.0]),
+ normal=np.array([1.0, 0.0, 0.0]),
+ base_length=1.0,
+ base_radius=0.01,
+ density=3000,
+ youngs_modulus=1e6,
+ )
+ for _ in range(n_rods)
+]
+rods_py = [
+ epy.CosseratRod.straight_rod(
+ n_elements=n_elems_per_rod,
+ start=np.zeros(3),
+ direction=np.array([0.0, 0.0, 1.0]),
+ normal=np.array([1.0, 0.0, 0.0]),
+ base_length=1.0,
+ base_radius=0.01,
+ density=3000,
+ youngs_modulus=1e6,
+ )
+ for _ in range(n_rods)
+]
+rng = np.random.default_rng(43)
+for i in range(n_rods):
+ for key in ["rest_kappa", "rest_sigma"]:
+ shape = getattr(rods_py[i], key).shape
+ values = rng.random(size=shape)
+ getattr(rods_py[i], key)[...] = values.copy()
+ getattr(rods_cpp[i], key)[...] = values.copy()
+
+# Create block rod system
+block_cpp = epp.MemoryBlockCosseratRod(rods_cpp, range(n_rods))
+block_py = epy.MemoryBlockCosseratRod(rods_py, range(n_rods))
+
+# %%
+# Cross-check ghost indices
+# -------------------------
+np.testing.assert_array_equal(block_cpp.ghost_nodes_idx, block_py.ghost_nodes_idx)
+np.testing.assert_array_equal(block_cpp.ghost_elems_idx[:-1], block_py.ghost_elems_idx)
+np.testing.assert_array_equal(
+ block_cpp.ghost_voronoi_idx[:-2], block_py.ghost_voronoi_idx
+)
+
+
+# %%
+# Cross-check block memory
+# ------------------------
+def cross_check_block_memory():
+ for key in keys:
+ cpp_value = getattr(block_cpp, key)
+ py_value = getattr(block_py, key)
+ assert cpp_value.shape == py_value.shape
+ assert np.allclose(cpp_value, py_value), f"{key} is not equal"
+ print(f"{key} | values and shapes checked")
+
+
+cross_check_block_memory()
+
+# %%
+# Cross-check computing internal forces and torques
+# -------------------------------------------------
+block_cpp.compute_internal_forces_and_torques(0.0)
+block_py.compute_internal_forces_and_torques(0.0)
+
+cross_check_block_memory()
+
+# %%
+# Cross-check updating accelerations
+# ----------------------------------
+block_cpp.update_accelerations(0.0)
+block_py.update_accelerations(0.0)
+
+cross_check_block_memory()
+# %%
+# Cross-check updating kinematics
+# -----------------------------
+block_cpp.update_kinematics(0.0, 1.4)
+block_py.update_kinematics(0.0, 1.4)
+
+cross_check_block_memory()
+# %%
+# Cross-check updating dynamics
+# -----------------------------
+block_cpp.update_dynamics(0.0, 1.6)
+block_py.update_dynamics(0.0, 1.6)
+
+cross_check_block_memory()
+# %%
diff --git a/backend/examples/run_timoshenko.py b/backend/examples/run_timoshenko.py
new file mode 100644
index 000000000..c4603bcc9
--- /dev/null
+++ b/backend/examples/run_timoshenko.py
@@ -0,0 +1,176 @@
+"""
+Test case for many-rod simulation with CPP memory block in the back.
+"""
+
+import numpy as np
+
+pass
+from tqdm import tqdm
+
+import elastica as ea
+import elasticapp as epp
+
+# %%
+# Now that we have imported all the necessary classes, we want to create our beam system. We do this by combining all the modules we need to represent the physics that we to include in the simulation. In this case, that is the ``BaseSystemCollection``, ``Constraint``, ``Forcings`` and ``Damping`` because the simulation will consider a rod that is fixed in place on one end, and subject to an applied force on the other end.
+
+
+class TimoshenkoBeamSimulator(
+ ea.BaseSystemCollection, ea.Constraints, ea.Forcing, ea.CallBacks, ea.Damping
+):
+ pass
+
+
+timoshenko_sim = TimoshenkoBeamSimulator()
+timoshenko_sim.enable_block_supports(ea.CosseratRod, epp.MemoryBlockCosseratRod)
+
+# %%
+# Creating Rods
+# -------------
+# With our simulator set up, we can now define the numerical, material, and geometric properties.
+#
+# First we define the number of elements in the rod. Next, the material properties are defined for every rod. These are the Young's modulus, the Poisson ratio, the density and the viscous damping coefficient. Finally, the geometry of the rod also needs to be defined by specifying the location of the rod and its orientation, length and radius.
+#
+# All of the values defined here are done in SI units, though this is not strictly necessary. You can rescale properties however you want, as long as you use consistent units throughout the simulation. See `here `_ for an example of consistent units.
+#
+# In order to make the difference between a shearable and unshearable rod more clear, we are using a Poisson ratio of 99. This is an unphysical value, as Poisson ratios can not exceed 0.5, however, it is used here for demonstration purposes.
+
+# setting up test params
+simulation_time = 5
+
+n_elem = 100
+start = np.zeros((3,))
+direction = np.array([0.0, 0.0, 1.0])
+normal = np.array([0.0, 1.0, 0.0])
+base_length = 3.0
+base_radius = 0.25
+base_area = np.pi * base_radius**2
+density = 5000
+nu = 0.1 / 7 / density / base_area
+E = 1e6
+# For shear modulus of 1e4, nu is 99!
+poisson_ratio = 99
+shear_modulus = E / (poisson_ratio + 1.0)
+
+# %%
+# With all of the rod's parameters set, we can now create a rod with the specificed properties and add the rod to the simulator system. **Important:** Make sure that any rods you create get added to the simulator system (``timoshenko_sim``), otherwise they will not be included in your simulation.
+
+n_rods = 10_000_000
+for _ in range(n_rods):
+ shearable_rod = ea.CosseratRod.straight_rod(
+ n_elem,
+ start,
+ direction,
+ normal,
+ base_length,
+ base_radius,
+ density,
+ youngs_modulus=E,
+ shear_modulus=shear_modulus,
+ )
+ timoshenko_sim.append(shearable_rod)
+
+# %%
+# Adding Damping
+# --------------
+# With the rod added to the simulator, we can add damping to the rod. We do this using the ``.dampen()`` option and the ``AnalyticalLinearDamper``. We are modifying ``timoshenko_sim`` simulator to ``dampen`` the ``shearable_rod`` object using ``AnalyticalLinearDamper`` type of dissipation (damping) model.
+#
+# We also need to define ``damping_constant`` and simulation ``time_step`` and pass in ``.using()`` method.
+
+dl = base_length / n_elem
+dt = 0.07 * dl
+timoshenko_sim.dampen(shearable_rod).using(
+ ea.AnalyticalLinearDamper,
+ damping_constant=nu,
+ time_step=dt,
+)
+
+# %%
+# Adding Boundary Conditions
+# --------------------------
+# With the rod added to the system, we need to apply boundary conditions. The first condition we will apply is fixing the location of one end of the rod. We do this using the ``.constrain()`` option and the ``OneEndFixedRod`` boundary condition. We are modifying the ``timoshenko_sim`` simulator to ``constrain`` the ``shearable_rod`` object using the ``OneEndFixedRod`` type of constraint.
+#
+# We also need to define which node of the rod is being constrained. We do this by passing the index of the nodes that we want to constrain to ``constrained_position_idx``. Here we are fixing the first node in the rod. In order to keep the rod from rotating around the fixed node, we also need to constrain an element between two nodes. This fixes the orientation of the rod. We do this by passing the index of the element that we want to fix to ``constrained_director_idx``. Like with the position, we are fixing the first element of the rod. Together, this constrains the position and orientation of the rod at the origin.
+
+# One end of the rod is now fixed in place
+timoshenko_sim.constrain(shearable_rod).using(
+ ea.OneEndFixedBC, constrained_position_idx=(0,), constrained_director_idx=(0,)
+)
+
+# %%
+# The next boundary condition that we want to apply is the endpoint force. Similarly to how we constrained one of the points, we want the ``timoshenko_sim`` simulator to ``add_forcing_to`` the ``shearable_rod`` object using the ``EndpointForces`` type of forcing. This ``EndpointForces`` applies forces to both ends of the rod. We want to apply a negative force in the :math:`d_1` direction, but only at the end of the rod. We do this by specifying the force vector to be applied at each end as ``origin_force`` and ``end_force``. We also want to ramp up the force over time, so we make the force take some ``ramp_up_time`` to reach its steady-state value. This helps avoid numerical errors due to discontinuities in the applied force.
+
+# Forces added to the rod
+end_force = np.array([-15.0, 0.0, 0.0])
+timoshenko_sim.add_forcing_to(shearable_rod).using(
+ ea.EndpointForces, 0.0 * end_force, end_force, ramp_up_time=simulation_time / 2.0
+)
+
+# %%
+# Add Unshearable Rod
+# -------------------
+#
+# Along with the shearable rod, we also want to add an unshearable rod to be able to compare the difference between the two. We do this the same way we did for the first rod, however, because this rod is unsherable, we need to change the Poisson ratio to make the rod unsherable. For a truely unsheraable rod, you would need a Poisson ratio of -1.0, however, this causes the system to be numerically unstable, so instead we make the system nearly unshearable by using a Poisson ratio of -0.85.
+
+# Start into the plane
+unshearable_start = np.array([0.0, -1.0, 0.0])
+shear_modulus = E / (-0.7 + 1.0)
+unshearable_rod = ea.CosseratRod.straight_rod(
+ n_elem,
+ unshearable_start,
+ direction,
+ normal,
+ base_length,
+ base_radius,
+ density,
+ youngs_modulus=E,
+ # Unshearable rod needs G -> inf, which is achievable with -ve poisson ratio
+ shear_modulus=shear_modulus,
+)
+
+timoshenko_sim.append(unshearable_rod)
+
+# add damping
+timoshenko_sim.dampen(unshearable_rod).using(
+ ea.AnalyticalLinearDamper,
+ damping_constant=nu,
+ time_step=dt,
+)
+# add boundary conditions
+timoshenko_sim.constrain(unshearable_rod).using(
+ ea.OneEndFixedBC, constrained_position_idx=(0,), constrained_director_idx=(0,)
+)
+timoshenko_sim.add_forcing_to(unshearable_rod).using(
+ ea.EndpointForces, 0.0 * end_force, end_force, ramp_up_time=simulation_time / 2.0
+)
+
+# %%
+# System Finalization
+# -------------------
+# We have now added all the necessary rods and boundary conditions to our system. The last thing we need to do is finalize the system. This goes through the system, rearranges things, and precomputes useful quantities to prepare the system for simulation.
+#
+# As a note, if you make any changes to the rod after calling finalize, you will need to re-setup the system. This requires rerunning all cells above this point.
+
+timoshenko_sim.finalize()
+
+# %%
+# Define Simulation Time
+# ----------------------
+# The last thing we need to do decide how long we want the simulation to run for and what timestepping method to use. Currently, the PositionVerlet algorithim is suggested default method.
+#
+# In this example, we are trying to match a steady-state solution by temporally evolving our system to reach equilibrium. As such, there is a tradeoff between letting the simulation run long enough to reach the equilibrium and waiting around for the simulation to be done. Here we are running the simulation for 10 seconds, this produces reasonable agreement with the analytical solution without taking to long to finish. If you run the simulation for longer, you will get better agreement with the analytical solution.
+
+timestepper = ea.PositionVerlet()
+# timestepper = PEFRL()
+
+total_steps = int(simulation_time / dt)
+print("Total steps", total_steps)
+
+# %%
+# Run Simulation
+# --------------
+#
+# We are now ready to perform the simulation. To run the simulation, we ``integrate`` the ``timoshenko_sim`` system using the ``timestepper`` method until ``final_time`` by taking ``total_steps``. As currently setup, the beam simulation takes about 1 minute to run.
+
+time = 0.0
+for i in tqdm(range(total_steps)):
+ time = timestepper.step(timoshenko_sim, time, dt)
diff --git a/backend/pyproject.toml b/backend/pyproject.toml
new file mode 100644
index 000000000..5ebedb6e0
--- /dev/null
+++ b/backend/pyproject.toml
@@ -0,0 +1,78 @@
+[project]
+name = "elasticapp"
+version = "0.0.4"
+description = "A CPP accelerated backend for PyElastica kernels"
+requires-python = ">=3.10"
+license = {text = "MIT License"}
+dependencies = [
+ "pybind11",
+ "pytest",
+ "pytest-rng",
+ "pytest-mock",
+]
+
+# this is a list that can be expanded
+authors = [{name = "Seung Hyun Kim", email = "skim0119@gmail.com"}]
+
+# classifiers can be added here, later
+classifiers = [
+]
+
+[project.urls]
+"Homepage" = "https://github.com/GazzolaLab/PyElastica/"
+"Bug Reports" = "https://github.com/GazzolaLab/PyElastica/issues"
+"Source" = "https://github.com/GazzolaLab/PyElastica/"
+"Documentation" = "https://docs.cosseratrods.org/"
+
+[build-system]
+requires = ["scikit-build-core", "pybind11"]
+build-backend = "scikit_build_core.build"
+
+[tool.scikit-build]
+build-dir = "build"
+
+[tool.scikit-build.cmake]
+version = ">=3.20"
+build-type = "Release"
+args = []
+
+# [tool.mypy]
+# files = "setup.py"
+# python_version = "3.7"
+# strict = true
+# show_error_codes = true
+# enable_error_code = ["ignore-without-code", "redundant-expr", "truthy-bool"]
+# warn_unreachable = true
+#
+# [[tool.mypy.overrides]]
+# module = ["ninja"]
+# ignore_missing_imports = true
+
+[tool.pytest.ini_options]
+minversion = "6.0"
+addopts = ["-ra", "--showlocals", "--strict-markers", "--strict-config"]
+xfail_strict = true
+filterwarnings = [
+ "error",
+ "ignore:(ast.Str|Attribute s|ast.NameConstant|ast.Num) is deprecated:DeprecationWarning:_pytest",
+]
+testpaths = ["tests/py"]
+
+[tool.cibuildwheel]
+test-command = "pytest {project}/tests"
+test-extras = ["test"]
+test-skip = ["*universal2:arm64"]
+# Setuptools bug causes collision between pypy and cpython artifacts
+before-build = "rm -rf {project}/build"
+
+[tool.ruff]
+target-version = "py37"
+
+[tool.ruff.lint]
+extend-select = [
+ "B", # flake8-bugbear
+ "I", # isort
+ "PGH", # pygrep-hooks
+ "RUF", # Ruff-specific
+ "UP", # pyupgrade
+]
diff --git a/backend/src/_api.cpp b/backend/src/_api.cpp
new file mode 100644
index 000000000..2e3ba0b58
--- /dev/null
+++ b/backend/src/_api.cpp
@@ -0,0 +1,603 @@
+#include
+#include
+#include
+#include "block.h"
+#include "block_view.h"
+#include "traits.h"
+#include "api.h"
+#include "operations.h"
+#include "cosserat_rod_system.h"
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+#include
+
+namespace py = pybind11;
+
+namespace elasticapp {
+
+#ifdef ELASTICAPP_GIL_RELEASE
+#define ELASTICAPP_GIL_RELEASE_SCOPE() py::gil_scoped_release gil_release;
+#else
+#define ELASTICAPP_GIL_RELEASE_SCOPE()
+#endif
+
+// Helper to find a variable type by name at runtime
+// Iterates through SystemType::Variables tuple and finds matching name
+template
+auto get_variable_by_name_impl(BlockView& view, const std::string& var_name) {
+ using CurrentVar = std::tuple_element_t;
+
+ // Check if current variable's name matches
+ if (var_name == std::string(CurrentVar::name)) {
+ return view.template get();
+ }
+
+ // Recurse to next variable if not last
+ if constexpr (Index + 1 < std::tuple_size_v) {
+ return get_variable_by_name_impl(view, var_name);
+ }
+
+ // If we've exhausted all variables, throw error
+ throw std::runtime_error("Unknown variable name: " + var_name);
+}
+
+// Helper function to map string variable names to types and call get()
+// This allows Python to use block.at(index).get("position") syntax
+template
+auto get_variable_by_name(BlockView& view, const std::string& var_name) {
+ using VariablesTuple = typename SystemType::Variables;
+
+ if constexpr (std::tuple_size_v > 0) {
+ return get_variable_by_name_impl(view, var_name);
+ } else {
+ throw std::runtime_error("System has no variables");
+ }
+}
+
+// Helper to find a variable type by name at runtime for Block
+// Iterates through SystemType::Variables tuple and finds matching name
+template
+auto get_block_variable_by_name_impl(BlockType& block, const std::string& var_name) {
+ using CurrentVar = std::tuple_element_t;
+
+ // Check if current variable's name matches
+ if (var_name == std::string(CurrentVar::name)) {
+ return block.template get();
+ }
+
+ // Recurse to next variable if not last
+ if constexpr (Index + 1 < std::tuple_size_v) {
+ return get_block_variable_by_name_impl(block, var_name);
+ }
+
+ // If we've exhausted all variables, throw error
+ throw std::runtime_error("Unknown variable name: " + var_name);
+}
+
+// Helper function to map string variable names to types and call get() for Block
+// This allows Python to use block.get("position") syntax
+template
+auto get_block_variable_by_name(BlockType& block, const std::string& var_name) {
+ using VariablesTuple = typename BlockType::Variables;
+
+ if constexpr (std::tuple_size_v > 0) {
+ return get_block_variable_by_name_impl(block, var_name);
+ } else {
+ throw std::runtime_error("System has no variables");
+ }
+}
+
+// Helper to reset ghost for a variable by name
+template
+void reset_ghost_for_variable_by_name_impl(BlockType& block, const std::string& var_name) {
+ using CurrentVar = std::tuple_element_t;
+
+ // Check if current variable's name matches
+ if (var_name == std::string(CurrentVar::name)) {
+ block.template reset_ghost_for_variable();
+ return;
+ }
+
+ // Recurse to next variable if not last
+ if constexpr (Index + 1 < std::tuple_size_v) {
+ reset_ghost_for_variable_by_name_impl(block, var_name);
+ } else {
+ throw std::runtime_error("Unknown variable name: " + var_name);
+ }
+}
+
+// Helper function to reset ghost for a variable by name
+template
+void reset_ghost_for_variable_by_name(BlockType& block, const std::string& var_name) {
+ using VariablesTuple = typename BlockType::Variables;
+
+ if constexpr (std::tuple_size_v > 0) {
+ reset_ghost_for_variable_by_name_impl(block, var_name);
+ } else {
+ throw std::runtime_error("System has no variables");
+ }
+}
+
+// Helper to convert Eigen Block view to numpy array
+// For Scalar variables (rows == 1), returns a 1D array instead of 2D
+// For Matrix variables (rows == 9), returns a 3D array (3, 3, N) directly
+template
+py::object block_to_numpy(BlockExpr&& block_expr, py::object parent) {
+ // Evaluate the expression to get dimensions
+ auto rows = static_cast(block_expr.rows());
+ auto cols = static_cast(block_expr.cols());
+
+ // Get actual strides from the Eigen Block expression
+ // For Eigen Blocks, innerStride() is the stride between elements in the same row/column
+ // and outerStride() is the stride between rows/columns depending on storage order
+ // For column-major: innerStride() = 1 (between rows), outerStride() = underlying_rows (between columns)
+ // For row-major: innerStride() = 1 (between columns), outerStride() = underlying_cols (between rows)
+ auto inner_stride = static_cast(block_expr.innerStride() * sizeof(double));
+ auto outer_stride = static_cast(block_expr.outerStride() * sizeof(double));
+
+ // For Scalar variables (rows == 1), return as 1D array
+ if (rows == 1) {
+ // For 1D array, stride is the column stride
+ py::ssize_t stride;
+ if constexpr (IsColMajor) {
+ // Column-major: stride between columns is outer_stride
+ stride = outer_stride;
+ } else {
+ // Row-major: stride between columns is inner_stride
+ stride = inner_stride;
+ }
+
+ return py::array_t(
+ {cols},
+ {stride},
+ const_cast(block_expr.data()),
+ parent // Keep parent object alive
+ );
+ }
+
+ // For Matrix variables (rows == 9), return as 3D array (3, 3, N)
+ // Matrix variables represent 3x3 matrices stored as flattened 9-element vectors
+ // Storage order: [m00, m10, m20, m01, m11, m21, m02, m12, m22] (column-major)
+ if (rows == 9) {
+ // For (3, 3, N) view from (9, N) column-major:
+ // Mapping: arr_3d[a, b, c] -> arr_2d[a*3 + b, c]
+ // Strides:
+ // - Page stride (a dimension): 3 * row_stride (to skip 3 rows)
+ // - Row stride (b dimension): row_stride (to skip 1 row)
+ // - Col stride (c dimension): col_stride (to skip 1 column)
+ py::ssize_t page_stride, row_stride_3d, col_stride_3d;
+ if constexpr (IsColMajor) {
+ // Column-major: row_stride = inner_stride, col_stride = outer_stride
+ page_stride = 3 * inner_stride; // Stride for first 3x3 dimension (skip 3 rows)
+ row_stride_3d = inner_stride; // Stride between rows in 3x3 (skip 1 row)
+ col_stride_3d = outer_stride; // Stride between columns (N dimension)
+ } else {
+ // Row-major: row_stride = outer_stride, col_stride = inner_stride
+ page_stride = 3 * outer_stride;
+ row_stride_3d = outer_stride;
+ col_stride_3d = inner_stride;
+ }
+
+ // Use py::buffer_info for 3D arrays with custom strides
+ std::vector shape = {3, 3, cols};
+ std::vector strides = {page_stride, row_stride_3d, col_stride_3d};
+ py::buffer_info buf_info(
+ const_cast(block_expr.data()),
+ sizeof(double),
+ py::format_descriptor::format(),
+ 3,
+ shape,
+ strides
+ );
+ return py::array_t(buf_info, parent);
+ }
+
+ // For Vector variables (rows == 3), return as 2D array
+ // For numpy, strides are in bytes and represent the step size for each dimension
+ // For column-major (Eigen default): row_stride = inner_stride, col_stride = outer_stride
+ // For row-major: row_stride = outer_stride, col_stride = inner_stride
+ py::ssize_t row_stride, col_stride;
+ if constexpr (IsColMajor) {
+ // Column-major: stride between rows is inner_stride, between columns is outer_stride
+ row_stride = inner_stride;
+ col_stride = outer_stride;
+ } else {
+ // Row-major: stride between rows is outer_stride, between columns is inner_stride
+ row_stride = outer_stride;
+ col_stride = inner_stride;
+ }
+
+ // Create numpy array view (non-owning) with correct strides
+ return py::array_t(
+ {rows, cols},
+ {row_stride, col_stride},
+ const_cast(block_expr.data()),
+ parent // Keep parent object alive
+ );
+}
+
+PYBIND11_MODULE(_memory_block, m) {
+ m.doc() = R"pbdoc(
+ Elasticapp BlockCosseratRod module
+ ----------------------------------------
+
+ Provides BlockCosseratRod class for 2D array management with Eigen backend.
+ )pbdoc";
+
+ // Forward declare BlockRodSystemView so it can be used as a return type
+ using BlockRodSystemViewType = BlockRodSystem::View;
+
+ // Thread management functions (only available when threading is enabled)
+ // Disable Eigen's internal threading to prevent oversubscription with OpenMP
+ // We use OpenMP for explicit parallelization, so Eigen should use single-threaded operations
+ Eigen::setNbThreads(1);
+
+ m.def("set_num_threads", [](int num_threads) {
+ if (num_threads > 0) {
+ omp_set_num_threads(num_threads);
+ // Ensure Eigen uses single thread to avoid oversubscription
+ Eigen::setNbThreads(1);
+ }
+ },
+ R"pbdoc(
+ Set the number of OpenMP threads to use for parallel operations.
+
+ Args:
+ num_threads: Number of threads to use (must be > 0).
+ If 0 or negative, OpenMP will use its default
+ (typically all available CPU cores).
+
+ Note:
+ This affects all subsequent parallel regions in the code.
+ The environment variable OMP_NUM_THREADS can also be used
+ to control thread count.
+ )pbdoc",
+ py::arg("num_threads"));
+
+ m.def("get_num_threads", []() {
+ return omp_get_num_threads();
+ },
+ R"pbdoc(
+ Get the current number of threads in the current parallel region.
+
+ Returns:
+ int: Number of threads (returns 1 if called outside a parallel region).
+ )pbdoc");
+
+ m.def("get_max_threads", []() {
+ return omp_get_max_threads();
+ },
+ R"pbdoc(
+ Get the maximum number of threads that can be used.
+
+ Returns:
+ int: Maximum number of threads available.
+ )pbdoc");
+
+ m.def("get_thread_num", []() {
+ return omp_get_thread_num();
+ },
+ R"pbdoc(
+ Get the current thread number (0 to num_threads-1).
+
+ Returns:
+ int: Current thread number (returns 0 if called outside a parallel region).
+ )pbdoc");
+
+ // BlockRodSystem class
+ py::class_(m, "BlockRodSystem")
+ .def(py::init([](py::object n_elems_per_rod_obj) {
+ std::vector n_elems_per_rod;
+
+ // Handle list, tuple, or numpy array
+ if (py::isinstance(n_elems_per_rod_obj)) {
+ py::list lst = n_elems_per_rod_obj.cast();
+ for (auto item : lst) {
+ n_elems_per_rod.push_back(item.cast());
+ }
+ } else if (py::isinstance(n_elems_per_rod_obj)) {
+ py::tuple tup = n_elems_per_rod_obj.cast();
+ for (auto item : tup) {
+ n_elems_per_rod.push_back(item.cast());
+ }
+ } else if (py::isinstance(n_elems_per_rod_obj)) {
+ py::array arr = n_elems_per_rod_obj.cast();
+ auto buf = arr.request();
+ if (buf.ndim != 1) {
+ throw std::runtime_error("numpy array must be 1-dimensional");
+ }
+ // Properly convert numpy array elements to std::size_t
+ // Handle different integer types safely
+ n_elems_per_rod.reserve(buf.size);
+ if (buf.itemsize == sizeof(std::int32_t)) {
+ auto* ptr = static_cast(buf.ptr);
+ for (py::ssize_t i = 0; i < buf.size; ++i) {
+ n_elems_per_rod.push_back(static_cast(ptr[i]));
+ }
+ } else if (buf.itemsize == sizeof(std::int64_t)) {
+ auto* ptr = static_cast(buf.ptr);
+ for (py::ssize_t i = 0; i < buf.size; ++i) {
+ n_elems_per_rod.push_back(static_cast(ptr[i]));
+ }
+ } else if (buf.itemsize == sizeof(std::size_t)) {
+ auto* ptr = static_cast(buf.ptr);
+ n_elems_per_rod.assign(ptr, ptr + buf.size);
+ } else {
+ // Fallback: iterate and cast each element
+ for (py::ssize_t i = 0; i < buf.size; ++i) {
+ py::object item = arr.attr("__getitem__")(i);
+ n_elems_per_rod.push_back(item.cast());
+ }
+ }
+ } else {
+ throw std::runtime_error("n_elems_per_rod must be a list, tuple, or numpy array");
+ }
+
+ return new BlockRodSystem(n_elems_per_rod);
+ }),
+ R"pbdoc(
+ Create a BlockRodSystem from list of element counts per rod.
+
+ Args:
+ n_elems_per_rod: List, tuple, or numpy array of integers representing
+ number of elements in each rod
+ )pbdoc",
+ py::arg("n_elems_per_rod"))
+ .def_property_readonly("n_systems", [](const BlockRodSystem& block) {
+ return block.n_systems();
+ },
+ R"pbdoc(
+ Number of systems (rods) in the block.
+ )pbdoc")
+ .def_property_readonly("shape", [](const BlockRodSystem& block) {
+ auto shape = block.shape();
+ return py::make_tuple(shape.first, shape.second);
+ },
+ R"pbdoc(
+ Get the shape of the block as (depth, width).
+
+ Returns:
+ tuple: (depth, width) tuple representing the block dimensions
+ )pbdoc")
+ .def("as_ref", [](const BlockRodSystem& block) {
+ return block_to_numpy(block.data(), py::cast(block));
+ },
+ R"pbdoc(
+ Get a numpy array view of the entire block data.
+
+ Returns:
+ numpy.ndarray: A writable numpy array view into the block's data.
+ The array does not own the data.
+ )pbdoc",
+ py::keep_alive<0, 1>())
+ .def("system_start_index", [](const BlockRodSystem& block, std::size_t index) {
+ return block.system_start_index(index);
+ },
+ R"pbdoc(
+ Get the starting column index for a specific rod.
+
+ Args:
+ index: Index of the rod
+
+ Returns:
+ int: Starting column index for the rod in the block
+ )pbdoc",
+ py::arg("index"))
+ .def("at", [](BlockRodSystem& block, std::size_t index) -> BlockRodSystemViewType {
+ return block.at(index);
+ },
+ R"pbdoc(
+ Get a view for a specific rod.
+
+ Args:
+ index: Index of the rod
+
+ Returns:
+ BlockRodSystemView: View object for accessing variables of this rod
+ )pbdoc",
+ py::arg("index"), py::return_value_policy::reference_internal)
+ .def("get", [](BlockRodSystem& block, const std::string& var_name) {
+ auto block_expr = get_block_variable_by_name(block, var_name);
+ return block_to_numpy(block_expr, py::cast(block));
+ },
+ R"pbdoc(
+ Get a variable by name as a numpy array view across all rods.
+
+ Args:
+ var_name: Name of the variable (e.g., "position", "velocity", "director")
+
+ Returns:
+ numpy.ndarray: A writable numpy array view into the variable's data
+ across all rods. The array does not own the data.
+ )pbdoc",
+ py::arg("var_name"), py::keep_alive<0, 1>())
+ .def("compute_internal_forces_and_torques", [](BlockRodSystem& block, double time) {
+ ELASTICAPP_GIL_RELEASE_SCOPE();
+ block.compute_internal_forces_and_torques(time);
+ },
+ R"pbdoc(
+ Compute internal forces and torques for all rods in the block.
+
+ This operation computes the internal forces and torques based on the
+ current state of the rods (positions, velocities, etc.).
+
+ Args:
+ time: Current simulation time.
+ )pbdoc",
+ py::arg("time"))
+ .def("compute_strains", [](BlockRodSystem& block, double time) {
+ ELASTICAPP_GIL_RELEASE_SCOPE();
+ block.compute_strains(time);
+ },
+ R"pbdoc(
+ Compute strains for all rods in the block.
+
+ This operation computes the shear/stretch strains and bending/twist strains
+ based on the current state of the rods.
+
+ Args:
+ time: Current simulation time (included for API compatibility, not used in implementation).
+ )pbdoc",
+ py::arg("time"))
+ .def("update_accelerations", [](BlockRodSystem& block, double time) {
+ ELASTICAPP_GIL_RELEASE_SCOPE();
+ block.update_accelerations(time);
+ },
+ R"pbdoc(
+ Update accelerations based on forces and torques.
+
+ This operation updates the acceleration variables based on the
+ computed forces and torques.
+
+ Args:
+ time: Current simulation time (included for API compatibility, not used in implementation).
+ )pbdoc",
+ py::arg("time"))
+ .def("zeroed_out_external_forces_and_torques", [](BlockRodSystem& block, double time) {
+ ELASTICAPP_GIL_RELEASE_SCOPE();
+ block.zeroed_out_external_forces_and_torques(time);
+ },
+ R"pbdoc(
+ Zero out external forces and torques for all rods.
+
+ This operation sets all external forces and torques to zero,
+ typically called at the beginning of each time step.
+
+ Args:
+ time: Current simulation time (included for API compatibility, not used in implementation).
+ )pbdoc",
+ py::arg("time"))
+ .def("update_kinematics", [](BlockRodSystem& block, double time, double prefac) {
+ ELASTICAPP_GIL_RELEASE_SCOPE();
+ block.update_kinematics(prefac);
+ },
+ R"pbdoc(
+ Update kinematics (position, director) for all rods.
+
+ This operation updates the kinematic variables based on velocity and omega.
+ Updates: position += prefac * velocity, director = R(prefac * omega) @ director
+
+ Args:
+ time: Current time (for compatibility with Python interface, not used in C++)
+ prefac: Integration prefactor (e.g., time step dt)
+ )pbdoc",
+ py::arg("time"), py::arg("prefac"))
+ .def("update_dynamics", [](BlockRodSystem& block, double time, double prefac) {
+ ELASTICAPP_GIL_RELEASE_SCOPE();
+ block.update_dynamics(prefac);
+ },
+ R"pbdoc(
+ Update dynamics (velocity, omega) for all rods.
+
+ This operation updates the dynamic variables based on acceleration and alpha.
+ Updates: velocity += prefac * acceleration, omega += prefac * alpha
+
+ Args:
+ time: Current time (for compatibility with Python interface, not used in C++)
+ prefac: Integration prefactor (e.g., time step dt)
+ )pbdoc",
+ py::arg("time"), py::arg("prefac"))
+ .def_property_readonly("ghost_nodes_idx", [](const BlockRodSystem& block) {
+ auto indices = block.ghost_nodes_idx();
+ // Convert to numpy array (pybind11 will handle the conversion automatically)
+ return py::cast(indices);
+ },
+ R"pbdoc(
+ Get indices of ghost nodes between rods.
+
+ Returns:
+ numpy.ndarray: An array of ghost node indices (length: n_rods - 1).
+ The array does not own the data.
+ )pbdoc")
+ .def_property_readonly("ghost_elems_idx", [](const BlockRodSystem& block) {
+ auto indices = block.ghost_elems_idx();
+ // Convert to numpy array (pybind11 will handle the conversion automatically)
+ return py::cast(indices);
+ },
+ R"pbdoc(
+ Get indices of ghost elements between rods.
+
+ Returns:
+ numpy.ndarray: An array of ghost element indices (length: 2 * (n_rods - 1)).
+ The array does not own the data.
+ )pbdoc")
+ .def_property_readonly("ghost_voronoi_idx", [](const BlockRodSystem& block) {
+ auto indices = block.ghost_voronoi_idx();
+ // Convert to numpy array (pybind11 will handle the conversion automatically)
+ return py::cast(indices);
+ },
+ R"pbdoc(
+ Get indices of ghost voronoi nodes between rods.
+
+ Returns:
+ numpy.ndarray: An array of ghost voronoi indices (length: 3 * (n_rods - 1)).
+ The array does not own the data.
+ )pbdoc")
+ .def("reset_ghost_for_variable", [](BlockRodSystem& block, const std::string& var_name) {
+ // Helper to reset ghost for a variable by name
+ reset_ghost_for_variable_by_name(block, var_name);
+ },
+ R"pbdoc(
+ Reset ghost values for a specific variable by name.
+
+ Args:
+ var_name: Name of the variable (e.g., "position", "velocity", "director")
+ )pbdoc",
+ py::arg("var_name"))
+ .def("reset_ghost", [](BlockRodSystem& block) {
+ block.reset_ghost();
+ },
+ R"pbdoc(
+ Reset ghost values for all variables.
+
+ This operation sets all ghost node/element/voronoi values to their
+ default ghost_value as defined in each variable type.
+ )pbdoc");
+
+ // BlockRodSystemView class
+ py::class_(m, "BlockRodSystemView")
+ .def_property_readonly("shape", [](const BlockRodSystemViewType& view) {
+ auto shape = view.shape();
+ return py::make_tuple(shape.first, shape.second);
+ },
+ R"pbdoc(
+ Get the shape of the view as (depth, width).
+
+ Returns:
+ tuple: (depth, width) tuple representing the view dimensions
+ )pbdoc")
+ .def("as_ref", [](const BlockRodSystemViewType& view) {
+ return block_to_numpy(view.data(), py::cast(view));
+ },
+ R"pbdoc(
+ Get a numpy array view of the entire view data.
+
+ Returns:
+ numpy.ndarray: A writable numpy array view into the view's data.
+ The array does not own the data.
+ )pbdoc",
+ py::keep_alive<0, 1>())
+ .def("get", [](BlockRodSystemViewType& view, const std::string& var_name) {
+ auto block_expr = get_variable_by_name(view, var_name);
+ return block_to_numpy(block_expr, py::cast(view));
+ },
+ R"pbdoc(
+ Get a variable by name as a numpy array view.
+
+ Args:
+ var_name: Name of the variable (e.g., "position", "velocity", "director")
+
+ Returns:
+ numpy.ndarray: A writable numpy array view into the variable's data.
+ The array does not own the data.
+ )pbdoc",
+ py::arg("var_name"), py::keep_alive<0, 1>());
+
+}
+} // namespace elasticapp
diff --git a/backend/src/_version.cpp b/backend/src/_version.cpp
new file mode 100644
index 000000000..e4e99d909
--- /dev/null
+++ b/backend/src/_version.cpp
@@ -0,0 +1,27 @@
+#include
+#include "version.h"
+
+namespace py = pybind11;
+
+PYBIND11_MODULE(version, m) {
+ m.doc() = R"pbdoc(
+ Elasticapp version module
+ -------------------------
+
+ Provides version information for the elasticapp package.
+ )pbdoc";
+
+ // Version function
+ m.def("version", &elasticapp::version, R"pbdoc(
+ Return the current version of elasticapp.
+
+ Returns:
+ str: The version string
+ )pbdoc");
+
+#ifdef VERSION_INFO
+ m.attr("__version__") = elasticapp::version();
+#else
+ m.attr("__version__") = "dev";
+#endif
+}
diff --git a/backend/src/api.h b/backend/src/api.h
new file mode 100644
index 000000000..e17670b42
--- /dev/null
+++ b/backend/src/api.h
@@ -0,0 +1,21 @@
+#pragma once
+
+#include "block.h"
+#include "cosserat_rod_system.h"
+#include "operations.h"
+
+namespace elasticapp {
+
+// BlockRodSystem is a CRTP mix of CosseratRodSystem, Block, and CosseratRodOperations
+// It combines System functionality with Block data storage and rod-specific operations
+// This allows Block to have access to all CosseratRodSystem variables and operations
+// Block inherits from:
+// - CosseratRodSystem (system variables and methods)
+// - CosseratRodOperations> (rod-specific operations)
+// So it has:
+// - All CosseratRodSystem methods
+// - All Block methods
+// - All CosseratRodOperations methods (compute_internal_forces_and_torques, etc.)
+// - Access to CosseratRodSystem::Variables for template metaprogramming
+using BlockRodSystem = Block;
+} // namespace elasticapp
diff --git a/backend/src/block.h b/backend/src/block.h
new file mode 100644
index 000000000..0f9cb4514
--- /dev/null
+++ b/backend/src/block.h
@@ -0,0 +1,280 @@
+#pragma once
+
+#include
+#include
+#include
+#include
+#include "traits.h"
+#include "system.h"
+#include "block_get_impl.h"
+#include "block_view.h"
+#include "operations.h"
+
+namespace elasticapp {
+
+// Block class using CRTP with SystemModel and Operations
+// Block inherits from SystemModel, giving it access to all system variables and methods
+// Block also inherits from OperationsType using CRTP, allowing for extensible operations
+// SystemModel must be a System type
+// OperationsType must be a template class that takes the Block type as a template parameter
+template class OperationsType = DefaultOperations>
+class Block : public SystemType, public OperationsType> {
+public:
+
+ using System = SystemType;
+ using Variables = typename SystemType::Variables;
+ using View = BlockView;
+
+ constexpr static std::size_t GHOST_NODE_WIDTH = 1; // Size of ghost for node variables.
+
+ // Constructor: takes list of element counts per rod
+ explicit Block(const std::vector& n_elems_per_rod) {
+ // Validate input: reject empty list or less than 6 total elements
+ if (n_elems_per_rod.empty()) {
+ throw std::invalid_argument("n_elems_per_rod cannot be empty");
+ }
+
+ std::size_t total_elements = 0;
+ for (std::size_t n_elems : n_elems_per_rod) {
+ total_elements += n_elems;
+ }
+ if (total_elements < 6) {
+ throw std::invalid_argument("Total number of elements must be at least 6, got " +
+ std::to_string(total_elements));
+ }
+
+ compute_width_and_indices(n_elems_per_rod);
+ depth_ = SystemType::get_depth();
+ data_ = MatrixType(static_cast(depth_), static_cast(width_));
+ data_.setZero();
+ initialize_ghost_indices();
+ reset_ghost(); // Initialize all ghost values
+ }
+
+ std::pair shape() const {
+ return std::make_pair(
+ static_cast(data_.rows()),
+ static_cast(data_.cols())
+ );
+ }
+
+ // Accessors
+ std::size_t width() const { return width_; }
+ std::size_t depth() const { return depth_; }
+ std::size_t n_systems() const { return rod_start_indices_.size(); }
+
+ std::size_t system_start_index(std::size_t rod_index) const {
+ if (rod_index >= rod_start_indices_.size()) {
+ throw std::out_of_range("System index out of range");
+ }
+ return rod_start_indices_[rod_index];
+ }
+
+ MatrixType& data() { return data_; }
+ const MatrixType& data() const { return data_; }
+
+ // Get number of elements for a specific rod
+ std::size_t rod_n_elems(std::size_t rod_index) const {
+ if (rod_index >= rod_n_elems_.size()) {
+ throw std::out_of_range("System index out of range");
+ }
+ return rod_n_elems_[rod_index];
+ }
+
+ std::size_t rod_n_nodes(std::size_t rod_index) const {
+ return rod_n_elems(rod_index) + 1;
+ }
+
+ std::size_t rod_n_voronoi(std::size_t rod_index) const {
+ return rod_n_elems(rod_index) - 1;
+ }
+
+ // Get the n_elems_per_rod vector (for BlockView construction)
+ const std::vector& get_n_elems_per_rod() const { return rod_n_elems_; }
+
+ // Get ghost indices (return const references for efficiency)
+ inline const std::vector& ghost_nodes_idx() const {
+ return ghost_nodes_idx_;
+ }
+ inline const std::vector& ghost_elems_idx() const {
+ return ghost_elems_idx_;
+ }
+ inline const std::vector& ghost_voronoi_idx() const {
+ return ghost_voronoi_idx_;
+ }
+
+ // Reset ghost values for a specific variable
+ // Uses VariableTag::ghost_value and appropriate ghost indices based on placement
+ template
+ void reset_ghost_for_variable() {
+ static_assert(tuple_contains_v>,
+ "VariableTag is not a valid member of tuple SystemType::Variables");
+
+ // Compute row offset for this variable
+ constexpr std::size_t row_offset = compute_variable_offset();
+ constexpr std::size_t var_dimension = get_dimension_v;
+
+ // Get appropriate ghost indices based on placement
+ std::vector ghost_indices;
+ if constexpr (std::is_base_of_v) {
+ ghost_indices = ghost_nodes_idx();
+ } else if constexpr (std::is_base_of_v) {
+ ghost_indices = ghost_elems_idx();
+ } else if constexpr (std::is_base_of_v) {
+ ghost_indices = ghost_voronoi_idx();
+ }
+
+ // Set ghost values at each ghost index
+ // Note: ghost indices are in the full width coordinate system
+ // For OnElement and OnVoronoi variables, we need to ensure ghost indices
+ // are within the adjusted width (since get() returns adjusted width views)
+ const auto& ghost_val = VariableTag::ghost_value;
+
+ // Optimized: Add simd hint for inner loop (var_dimension is compile-time constant)
+ for (std::size_t ghost_col : ghost_indices) {
+ // Only reset ghost values that are within the adjusted width
+ // (ghost indices beyond adjusted width are not accessible via get())
+ IndexType data_col = static_cast(ghost_col);
+ #pragma omp simd
+ for (std::size_t row = 0; row < var_dimension; ++row) {
+ IndexType data_row = static_cast(row_offset + row);
+ data_(data_row, data_col) = ghost_val(static_cast(row), 0);
+ }
+ }
+ }
+
+ // Reset ghost values for all variables
+ // Iterates over all variables and calls reset_ghost_for_variable for each
+ void reset_ghost() {
+ reset_ghost_impl, 0>();
+ }
+
+ // Get a view for a specific variable across all rods
+ // Returns a view into the variable's data (rows) and adjusted columns based on placement
+ // - OnNode: full width
+ // - OnElement: width - 1
+ // - OnVoronoi: width - 2
+ // No data is copied - returns a reference to the same matrix
+ template
+ auto get() {
+ // Assert VariableTag is a valid member of tuple SystemType::Variables
+ static_assert(tuple_contains_v>,
+ "VariableTag is not a valid member of tuple SystemType::Variables");
+ // Compute row offset for this variable
+ constexpr std::size_t row_offset = compute_variable_offset();
+ constexpr std::size_t var_dimension = get_dimension_v;
+
+ // Adjust width based on placement type
+ std::size_t adjusted_width = width_;
+ if constexpr (std::is_base_of_v) {
+ adjusted_width = width_ > 0 ? width_ - 1 : 0;
+ } else if constexpr (std::is_base_of_v) {
+ adjusted_width = width_ > 1 ? width_ - 2 : 0;
+ }
+ // OnNode: no adjustment needed (full width)
+
+ // Return a view into the specific rows and adjusted columns
+ return get_block_slice(data_, row_offset, var_dimension, 0, adjusted_width);
+ }
+
+ template
+ auto get() const {
+ // Assert VariableTag is a valid member of tuple SystemType::Variables
+ static_assert(tuple_contains_v>,
+ "VariableTag is not a valid member of tuple SystemType::Variables");
+ // Compute row offset for this variable
+ constexpr std::size_t row_offset = compute_variable_offset();
+ constexpr std::size_t var_dimension = get_dimension_v;
+
+ // Adjust width based on placement type
+ std::size_t adjusted_width = width_;
+ if constexpr (std::is_base_of_v) {
+ adjusted_width = width_ > 0 ? width_ - 1 : 0;
+ } else if constexpr (std::is_base_of_v) {
+ adjusted_width = width_ > 1 ? width_ - 2 : 0;
+ }
+ // OnNode: no adjustment needed (full width)
+
+ // Return a view into the specific rows and adjusted columns
+ return get_block_slice(data_, row_offset, var_dimension, 0, adjusted_width);
+ }
+
+ // Create a view for a specific rod
+ View at(std::size_t rod_index);
+
+private:
+ MatrixType data_;
+ std::size_t width_;
+ std::size_t depth_;
+ std::vector rod_start_indices_;
+ std::vector rod_n_elems_; // Store n_elems for each rod
+
+ std::vector ghost_nodes_idx_; // Indices of ghost nodes between rods.
+ std::vector ghost_elems_idx_; // Indices of ghost elements between rods.
+ std::vector ghost_voronoi_idx_; // Indices of ghost voronoi nodes between rods.
+
+ void compute_width_and_indices(const std::vector& n_elems_per_rod) {
+ width_ = 0;
+ rod_start_indices_.clear();
+ rod_n_elems_.clear();
+ rod_start_indices_.reserve(n_elems_per_rod.size());
+ rod_n_elems_.reserve(n_elems_per_rod.size());
+
+ for (std::size_t n_elems : n_elems_per_rod) {
+ rod_start_indices_.push_back(width_);
+ rod_n_elems_.push_back(n_elems);
+ // Each rod has n_elems + 1 nodes
+ width_ += n_elems + 1 + GHOST_NODE_WIDTH;
+ }
+ width_ -= GHOST_NODE_WIDTH; // Remove the last ghost node width.
+ }
+
+ void initialize_ghost_indices() {
+ ghost_nodes_idx_.reserve(rod_n_elems_.size() * GHOST_NODE_WIDTH);
+ ghost_elems_idx_.reserve(rod_n_elems_.size() * (1 + GHOST_NODE_WIDTH));
+ ghost_voronoi_idx_.reserve(rod_n_elems_.size() * (2 + GHOST_NODE_WIDTH));
+
+ std::size_t cumulative_nodes = 0;
+ for (std::size_t i = 0; i < rod_n_elems_.size(); ++i) {
+ cumulative_nodes += rod_n_elems_[i] + 1; // n_elems + 1 = n_nodes
+ ghost_elems_idx_.push_back(cumulative_nodes - 1);
+ ghost_voronoi_idx_.push_back(cumulative_nodes - 2);
+ ghost_voronoi_idx_.push_back(cumulative_nodes - 1);
+ for (std::size_t j = 0; j < GHOST_NODE_WIDTH; ++j) {
+ ghost_nodes_idx_.push_back(cumulative_nodes + j);
+ ghost_elems_idx_.push_back(cumulative_nodes + j);
+ ghost_voronoi_idx_.push_back(cumulative_nodes + j);
+ }
+ cumulative_nodes += GHOST_NODE_WIDTH;
+ }
+
+ ghost_nodes_idx_.pop_back();
+ ghost_elems_idx_.pop_back();
+ ghost_voronoi_idx_.pop_back();
+ }
+
+ // Helper to iterate over all variables and reset ghost values
+ template
+ void reset_ghost_impl() {
+ using CurrentVar = std::tuple_element_t;
+ reset_ghost_for_variable();
+
+ // Recurse to next variable if not last
+ if constexpr (Index + 1 < std::tuple_size_v) {
+ reset_ghost_impl();
+ }
+ }
+
+};
+
+
+// Implement at() after BlockView is fully defined
+template class OperationsType>
+typename Block::View Block::at(std::size_t rod_index) {
+ std::size_t rod_start_col = system_start_index(rod_index);
+ std::size_t rod_n_elems = rod_n_elems_[rod_index];
+ return BlockView(data_, rod_index, rod_start_col, rod_n_elems);
+}
+
+
+} // namespace elasticapp
diff --git a/backend/src/block_get_impl.h b/backend/src/block_get_impl.h
new file mode 100644
index 000000000..3b70341d4
--- /dev/null
+++ b/backend/src/block_get_impl.h
@@ -0,0 +1,70 @@
+#pragma once
+
+#include
+#include
+#include
+#include "system.h"
+#include "variable_offsets.h"
+#include "traits.h"
+
+namespace elasticapp {
+
+// Uses C++17 fold expressions for a concise implementation
+// Helper to check if a type is in a tuple
+template
+struct tuple_contains_impl;
+
+template
+struct tuple_contains_impl> {
+ static constexpr bool value = (std::is_same_v || ...);
+};
+
+template
+constexpr bool tuple_contains_v = tuple_contains_impl::value;
+
+
+// Helper function to get number of columns for a variable based on placement
+template
+inline constexpr std::size_t get_variable_num_cols(std::size_t rod_n_nodes,
+ std::size_t rod_n_elems,
+ std::size_t rod_n_voronoi) {
+ if constexpr (std::is_base_of_v) {
+ return rod_n_nodes;
+ } else if constexpr (std::is_base_of_v) {
+ return rod_n_elems;
+ } else if constexpr (std::is_base_of_v) {
+ return rod_n_voronoi;
+ } else {
+ static_assert(std::is_base_of_v ||
+ std::is_base_of_v ||
+ std::is_base_of_v,
+ "VariableTag must have a Placement tag");
+ return 0; // Should never reach here
+ }
+}
+
+// Generic implementation helper for Block::get() and BlockRodSystemView::get()
+// Extracted to reduce code duplication between const and non-const versions
+// This can be used by both Block and BlockRodSystemView classes
+template
+inline auto get_impl(MatrixRef&& matrix,
+ std::size_t rod_n_nodes,
+ std::size_t rod_n_elems,
+ std::size_t rod_n_voronoi,
+ std::size_t rod_start_col) {
+ // Assert VariableTag is a valid member of tuple SystemType::Variables
+ static_assert(tuple_contains_v>,
+ "VariableTag is not a valid member of tuple SystemType::Variables");
+ // Compute row offset for this variable
+ constexpr std::size_t row_offset = compute_variable_offset();
+ constexpr std::size_t var_dimension = get_dimension_v;
+
+ // Determine number of columns based on placement
+ std::size_t num_cols = get_variable_num_cols(
+ rod_n_nodes, rod_n_elems, rod_n_voronoi);
+
+ // Return a view into the specific rows and columns
+ return get_block_slice(matrix, row_offset, var_dimension, rod_start_col, num_cols);
+}
+
+} // namespace elasticapp
diff --git a/backend/src/block_view.h b/backend/src/block_view.h
new file mode 100644
index 000000000..2e98627fd
--- /dev/null
+++ b/backend/src/block_view.h
@@ -0,0 +1,74 @@
+#pragma once
+
+#include
+#include
+#include
+#include "system.h"
+#include "block_get_impl.h"
+#include "traits.h"
+
+namespace elasticapp {
+
+// BlockView provides access to variables for a specific rod
+// This is templated on SystemType and uses template metaprogramming to expose variables
+template
+class BlockView {
+public:
+ // using Variables = typename SystemType::Variables; // tuple of variables
+ constexpr static std::size_t depth = SystemType::get_depth();
+
+ BlockView(MatrixType& data, std::size_t rod_index,
+ std::size_t rod_start_col, std::size_t rod_n_elems);
+
+ std::pair shape() const {
+ return std::make_pair(depth, rod_n_nodes_);
+ }
+
+ // Return a view into the submatrix (columns for this rod)
+ // Uses traits helper function to encapsulate Eigen-specific operations
+ auto data() {
+ return get_column_slice(data_, rod_start_col_, rod_n_nodes_);
+ }
+
+ auto data() const {
+ return get_column_slice(data_, rod_start_col_, rod_n_nodes_);
+ }
+
+ // Get a view for a specific variable
+ // Returns a view into the variable's data (rows) and this rod's columns
+ // No data is copied - returns a reference to the same matrix
+ template
+ auto get() {
+ return get_impl(
+ data_, rod_n_nodes_, rod_n_elems_, rod_n_voronoi_, rod_start_col_);
+ }
+
+ template
+ auto get() const {
+ return get_impl(
+ data_, rod_n_nodes_, rod_n_elems_, rod_n_voronoi_, rod_start_col_);
+ }
+
+protected:
+ MatrixType& data_;
+ std::size_t rod_index_;
+ std::size_t rod_n_elems_;
+ std::size_t rod_n_nodes_;
+ std::size_t rod_n_voronoi_;
+ std::size_t rod_start_col_;
+};
+
+// Note: block.h is included after BlockView definition to break circular dependency
+// The BlockView constructor implementation is in block.h after Block is fully defined
+template
+BlockView::BlockView(MatrixType& data, std::size_t rod_index,
+ std::size_t rod_start_col, std::size_t rod_n_elems)
+ : data_(data), rod_index_(rod_index),
+ rod_n_elems_(rod_n_elems),
+ rod_n_nodes_(rod_n_elems + 1),
+ rod_n_voronoi_((rod_n_elems > 0) ? rod_n_elems - 1 : 0),
+ rod_start_col_(rod_start_col)
+{
+}
+
+} // namespace elasticapp
diff --git a/backend/src/cosserat_equations.h b/backend/src/cosserat_equations.h
new file mode 100644
index 000000000..75a62f149
--- /dev/null
+++ b/backend/src/cosserat_equations.h
@@ -0,0 +1,607 @@
+#pragma once
+
+#include
+#include
+#include "cosserat_rod_system.h"
+#include "traits.h"
+#include "math/eigen_detail/eigen_linear_algebra.hpp"
+#include "math/eigen_detail/eigen_calculus.hpp"
+#include
+
+namespace elasticapp {
+
+// Forward declaration
+template class OperationsType>
+class Block;
+
+// Compute geometry from state (lengths, tangents, radius)
+// Updates: lengths, tangents, radius
+template
+inline void compute_geometry_from_state(BlockType& block) {
+ // Get variable views
+ auto&& position = block.template get();
+ auto&& volume = block.template get();
+ auto&& lengths = block.template get();
+ auto&& tangents = block.template get();
+ auto&& radius = block.template get();
+
+ // Compute position differences using difference_kernel
+ // position_diff = position[:, 1:] - position[:, :-1]
+ // This is equivalent to: difference_kernel(position)
+ auto position_diff = difference_kernel(position);
+
+ // Compute lengths = norm(position_diff) + 1e-14
+ // FIXME: 1e-14 is added to fix ghost lengths, which is 0, and causes division by zero error!
+ auto lengths_vec = batch_norm(position_diff);
+ const IndexType n_elems = lengths.cols();
+ #pragma omp parallel for simd schedule(static)
+ for (IndexType k = 0; k < n_elems; ++k) {
+ lengths(0, k) = lengths_vec(k) + 1e-14;
+ }
+
+ // Compute tangents = position_diff / lengths
+ #pragma omp parallel for simd schedule(static)
+ for (IndexType k = 0; k < n_elems; ++k) {
+ const double len = lengths(0, k);
+ tangents(0, k) = position_diff(0, k) / len;
+ tangents(1, k) = position_diff(1, k) / len;
+ tangents(2, k) = position_diff(2, k) / len;
+ }
+
+ // Compute radius from volume conservation: radius = sqrt(volume / (lengths * pi))
+ const double pi = 3.14159265358979323846;
+ #pragma omp parallel for simd schedule(static)
+ for (IndexType k = 0; k < n_elems; ++k) {
+ radius(0, k) = std::sqrt(volume(0, k) / (lengths(0, k) * pi));
+ }
+}
+
+// Compute all dilatations (element and voronoi)
+// Updates: lengths, tangents, radius, dilatation, voronoi_dilatation
+template
+inline void compute_all_dilatations(BlockType& block) {
+ // Get variable views
+ auto&& lengths = block.template get();
+ auto&& rest_lengths = block.template get();
+ auto&& dilatation = block.template get();
+ auto&& rest_voronoi_lengths = block.template get();
+ auto&& voronoi_dilatation = block.template get();
+
+ // Compute dilatation = lengths / rest_lengths
+ const IndexType n_elems = lengths.cols();
+ #pragma omp parallel for simd schedule(static)
+ for (IndexType k = 0; k < n_elems; ++k) {
+ dilatation(0, k) = lengths(0, k) / rest_lengths(0, k);
+ }
+
+ // Compute voronoi_lengths from average of lengths
+ // voronoi_lengths = 0.5 * (lengths[k+1] + lengths[k]) for k in [0, n_voronoi-1]
+ auto voronoi_lengths = average_kernel(lengths);
+
+ // Compute voronoi_dilatation = voronoi_lengths / rest_voronoi_lengths
+ const IndexType n_voronoi = voronoi_dilatation.cols();
+ #pragma omp parallel for simd schedule(static)
+ for (IndexType k = 0; k < n_voronoi; ++k) {
+ voronoi_dilatation(0, k) = voronoi_lengths(0, k) / rest_voronoi_lengths(0, k);
+ }
+}
+
+// Compute shear/stretch strains (sigma)
+// Updates: lengths, tangents, radius, dilatation, voronoi_dilatation, sigma
+template
+inline void compute_shear_stretch_strains(BlockType& block) {
+ // Get variable views
+ auto&& dilatation = block.template get();
+ auto&& director = block.template get();
+ auto&& tangents = block.template get();
+ auto&& sigma = block.template get();
+
+ // Compute sigma = dilatation * batch_matvec(director_collection, tangents) - z_vector
+ // director is stored as (9, n_elems) - flattened 3x3 matrices
+ // Storage order: [d00, d10, d20, d01, d11, d21, d02, d12, d22] (column-major)
+ // tangents is (3, n_elems)
+ // sigma is (3, n_elems)
+ const IndexType n_elems = sigma.cols();
+
+ #pragma omp parallel for simd schedule(static)
+ for (IndexType k = 0; k < n_elems; ++k) {
+ // Extract 3x3 director matrix for element k
+ // director[:, k] is (row-wise flattened 3x3 matrix)
+ // [d00, d01, d02]
+ // [d10, d11, d12]
+ // [d20, d21, d22]
+ double d00 = director(0, k);
+ double d01 = director(1, k);
+ double d02 = director(2, k);
+ double d10 = director(3, k);
+ double d11 = director(4, k);
+ double d12 = director(5, k);
+ double d20 = director(6, k);
+ double d21 = director(7, k);
+ double d22 = director(8, k);
+
+ // Compute director @ tangents for element k
+ // result = director_matrix * tangents[:, k]
+ double result0 = d00 * tangents(0, k) + d01 * tangents(1, k) + d02 * tangents(2, k);
+ double result1 = d10 * tangents(0, k) + d11 * tangents(1, k) + d12 * tangents(2, k);
+ double result2 = d20 * tangents(0, k) + d21 * tangents(1, k) + d22 * tangents(2, k);
+
+ // Compute sigma = dilatation * result - z_vector
+ const double dil = dilatation(0, k);
+ sigma(0, k) = dil * result0;
+ sigma(1, k) = dil * result1;
+ sigma(2, k) = dil * result2 - 1.0;
+ }
+}
+
+// Compute internal shear/stretch stresses from model
+// Updates: lengths, tangents, radius, dilatation, voronoi_dilatation, sigma, internal_stress
+template
+inline void compute_internal_shear_stretch_stresses_from_model(BlockType& block) {
+ // Get variable views
+ auto&& shear_matrix = block.template get();
+ auto&& sigma = block.template get();
+ auto&& rest_sigma = block.template get();
+ auto&& internal_stress = block.template get();
+
+ // Compute sigma_diff = sigma - rest_sigma
+ // sigma is (3, n_elems), rest_sigma is (3, n_elems)
+ const IndexType n_elems = sigma.cols();
+
+ // Compute internal_stress = batch_matvec(shear_matrix, sigma - rest_sigma)
+ // shear_matrix is stored as (9, n_elems) - flattened 3x3 matrices
+ // Storage order: [m00, m10, m20, m01, m11, m21, m02, m12, m22] (column-major)
+ #pragma omp parallel for simd schedule(static)
+ for (IndexType k = 0; k < n_elems; ++k) {
+ // Compute sigma_diff = sigma - rest_sigma for element k
+ double sigma_diff0 = sigma(0, k) - rest_sigma(0, k);
+ double sigma_diff1 = sigma(1, k) - rest_sigma(1, k);
+ double sigma_diff2 = sigma(2, k) - rest_sigma(2, k);
+
+ // Extract 3x3 shear_matrix for element k
+ // shear_matrix[:, k] is [m00, m10, m20, m01, m11, m21, m02, m12, m22]
+ double m00 = shear_matrix(0, k);
+ double m10 = shear_matrix(1, k);
+ double m20 = shear_matrix(2, k);
+ double m01 = shear_matrix(3, k);
+ double m11 = shear_matrix(4, k);
+ double m21 = shear_matrix(5, k);
+ double m02 = shear_matrix(6, k);
+ double m12 = shear_matrix(7, k);
+ double m22 = shear_matrix(8, k);
+
+ // Compute shear_matrix @ sigma_diff for element k
+ // result = shear_matrix * sigma_diff
+ internal_stress(0, k) = m00 * sigma_diff0 + m01 * sigma_diff1 + m02 * sigma_diff2;
+ internal_stress(1, k) = m10 * sigma_diff0 + m11 * sigma_diff1 + m12 * sigma_diff2;
+ internal_stress(2, k) = m20 * sigma_diff0 + m21 * sigma_diff1 + m22 * sigma_diff2;
+ }
+}
+
+// Compute internal forces for Cosserat rod
+// This is a template function that will be implemented later
+template
+inline void compute_internal_forces(BlockType& block) {
+ // Get variable views
+ auto&& director = block.template get();
+ auto&& internal_stress = block.template get();
+ auto&& dilatation = block.template get();
+ auto&& internal_forces = block.template get();
+ auto&& cosserat_internal_stress = block.template get();
+
+ // Compute cosserat_internal_stress = director^T @ internal_stress / dilatation
+ // director is stored as (9, n_elems) - flattened 3x3 matrices
+ // Storage order: [d00, d10, d20, d01, d11, d21, d02, d12, d22] (column-major)
+ // internal_stress is (3, n_elems) where n_elems includes ghost elements
+ // cosserat_internal_stress is (3, n_elems)
+ // Note: internal_stress has adjusted width (width - 1 for OnElement), but we need full width
+ // Get the full width from director which also has adjusted width
+ const IndexType n_elems = director.cols();
+
+ // Temporary matrix for cosserat_internal_stress
+ // MatrixType cosserat_internal_stress(3, n_elems);
+
+ // Compute cosserat_internal_stress = director^T @ internal_stress
+ // Python: cosserat_internal_stress[i, k] = sum_j(director_collection[j, i, k] * internal_stress[j, k])
+ // In C++: director_collection[j, i, k] = director(j*3 + i, k)
+ for (IndexType i = 0; i < 3; ++i) {
+ #pragma omp parallel for simd schedule(static)
+ for (IndexType k = 0; k < n_elems; ++k) {
+ double sum = 0.0;
+ for (IndexType j = 0; j < 3; ++j) {
+ // director_collection[j, i, k] = director(j*3 + i, k)
+ IndexType director_idx = j * 3 + i;
+ sum += director(director_idx, k) * internal_stress(j, k);
+ }
+ cosserat_internal_stress(i, k) = sum / dilatation(0, k);
+ }
+ }
+
+ // Reset ghost values for cosserat_internal_stress (OnElement)
+ block.template reset_ghost_for_variable();
+
+ // Compute internal_forces = two_point_difference_kernel(cosserat_internal_stress)
+ // internal_forces is OnNode (3, n_nodes) where n_nodes = n_elems + 1
+ two_point_difference_kernel(internal_forces, cosserat_internal_stress);
+}
+
+// Compute bending/twist strains (kappa)
+// Updates: kappa
+template
+inline void compute_bending_twist_strains(BlockType& block) {
+ // Get variable views
+ auto&& director = block.template get();
+ auto&& rest_voronoi_lengths = block.template get();
+ auto&& kappa = block.template get();
+
+ // director is stored as (9, n_elems) - flattened 3x3 matrices
+ // Storage order: [d00, d10, d20, d01, d11, d21, d02, d12, d22] (column-major)
+ // rest_voronoi_lengths is (1, n_voronoi) where n_voronoi = n_elems - 1
+ // kappa is (3, n_voronoi)
+ const IndexType n_voronoi = kappa.cols();
+
+ // Compute temp = inv_rotate(director_collection)
+ // inv_rotate computes relative rotation between consecutive director frames
+ // Python: temp = _inv_rotate(director_collection)
+ // temp has shape (3, n_voronoi) where n_voronoi = n_elems - 1
+ double temp_x, temp_y, temp_z;
+
+ // Compute inv_rotate manually: relative rotation between consecutive directors
+ // Python implementation computes cross products between consecutive director rows
+ #pragma omp parallel for schedule(static)
+ for (IndexType k = 0; k < n_voronoi; ++k) {
+ // Extract director matrices for k and k+1
+ // director[:, k] is (row-wise flattened 3x3 matrix)
+ // [d00, d01, d02]
+ // [d10, d11, d12]
+ // [d20, d21, d22]
+
+ // For director[k]:
+ double d00_k = director(0, k);
+ double d01_k = director(1, k);
+ double d02_k = director(2, k);
+ double d10_k = director(3, k);
+ double d11_k = director(4, k);
+ double d12_k = director(5, k);
+ double d20_k = director(6, k);
+ double d21_k = director(7, k);
+ double d22_k = director(8, k);
+
+ // For director[k+1]:
+ double d00_k1 = director(0, k + 1);
+ double d01_k1 = director(1, k + 1);
+ double d02_k1 = director(2, k + 1);
+ double d10_k1 = director(3, k + 1);
+ double d11_k1 = director(4, k + 1);
+ double d12_k1 = director(5, k + 1);
+ double d20_k1 = director(6, k + 1);
+ double d21_k1 = director(7, k + 1);
+ double d22_k1 = director(8, k + 1);
+
+ // Compute inv_rotate: cross product between consecutive director frames
+ temp_x = (d20_k1 * d10_k + d21_k1 * d11_k + d22_k1 * d12_k) -
+ (d10_k1 * d20_k + d11_k1 * d21_k + d12_k1 * d22_k);
+ temp_y = (d00_k1 * d20_k + d01_k1 * d21_k + d02_k1 * d22_k) -
+ (d20_k1 * d00_k + d21_k1 * d01_k + d22_k1 * d02_k);
+ temp_z = (d10_k1 * d00_k + d11_k1 * d01_k + d12_k1 * d02_k) -
+ (d00_k1 * d10_k + d01_k1 * d11_k + d02_k1 * d12_k);
+
+ double trace = (d00_k1 * d00_k + d01_k1 * d01_k + d02_k1 * d02_k) +
+ (d10_k1 * d10_k + d11_k1 * d11_k + d12_k1 * d12_k) +
+ (d20_k1 * d20_k + d21_k1 * d21_k + d22_k1 * d22_k);
+
+ // Clip the trace to between -1 and 3.
+ // Any deviation beyond this is numerical error
+ trace = std::clamp(trace, -1.0, 3.0);
+ double cos = 0.5 * trace - 0.5;
+ double theta = std::acos(cos) + 1e-14;
+ double magnitude = -0.5 * theta / std::sin(theta) / rest_voronoi_lengths(0, k);
+ kappa(0, k) = temp_x * magnitude;
+ kappa(1, k) = temp_y * magnitude;
+ kappa(2, k) = temp_z * magnitude;
+ }
+}
+
+// Compute internal bending/twist stresses from model
+// Updates: kappa, internal_couple
+template
+inline void compute_internal_bending_twist_stresses_from_model(BlockType& block) {
+ // Get variable views
+ auto&& kappa = block.template get();
+ auto&& rest_kappa = block.template get();
+ auto&& bend_matrix = block.template get();
+ auto&& internal_couple = block.template get();
+
+ // kappa is (3, n_voronoi), rest_kappa is (3, n_voronoi)
+ // bend_matrix is stored as (9, n_voronoi) - flattened 3x3 matrices
+ // internal_couple is (3, n_voronoi)
+ const IndexType n_voronoi = kappa.cols();
+
+ // Compute diff_kappa = kappa - rest_kappa
+ // MatrixType diff_kappa(3, n_voronoi);
+ auto&& diff_kappa = block.template get();
+
+ for (IndexType i = 0; i < 3; ++i) {
+ #pragma omp parallel for simd schedule(static)
+ for (IndexType k = 0; k < n_voronoi; ++k) {
+ diff_kappa(i, k) = kappa(i, k) - rest_kappa(i, k);
+ }
+ }
+
+ // Compute internal_couple = batch_matvec(bend_matrix, diff_kappa)
+ // bend_matrix is stored as (9, n_voronoi) - flattened 3x3 matrices
+ // Storage order: [m00, m10, m20, m01, m11, m21, m02, m12, m22] (column-major)
+ #pragma omp parallel for simd schedule(static)
+ for (IndexType k = 0; k < n_voronoi; ++k) {
+ // Extract 3x3 bend_matrix for voronoi k
+ // bend_matrix[:, k] is [m00, m10, m20, m01, m11, m21, m02, m12, m22]
+ double m00 = bend_matrix(0, k);
+ double m10 = bend_matrix(1, k);
+ double m20 = bend_matrix(2, k);
+ double m01 = bend_matrix(3, k);
+ double m11 = bend_matrix(4, k);
+ double m21 = bend_matrix(5, k);
+ double m02 = bend_matrix(6, k);
+ double m12 = bend_matrix(7, k);
+ double m22 = bend_matrix(8, k);
+
+ // Compute bend_matrix @ diff_kappa for voronoi k
+ // result = bend_matrix * diff_kappa[:, k]
+ internal_couple(0, k) = m00 * diff_kappa(0, k) + m01 * diff_kappa(1, k) + m02 * diff_kappa(2, k);
+ internal_couple(1, k) = m10 * diff_kappa(0, k) + m11 * diff_kappa(1, k) + m12 * diff_kappa(2, k);
+ internal_couple(2, k) = m20 * diff_kappa(0, k) + m21 * diff_kappa(1, k) + m22 * diff_kappa(2, k);
+ }
+}
+
+// Compute dilatation rate
+// Updates: dilatation_rate
+template
+inline void compute_dilatation_rate(BlockType& block) {
+ // Get variable views
+ auto&& position = block.template get();
+ auto&& velocity = block.template get();
+ auto&& lengths = block.template get();
+ auto&& rest_lengths = block.template get();
+ auto&& dilatation_rate = block.template get();
+
+ // position is (3, n_nodes), velocity is (3, n_nodes)
+ // lengths is (1, n_elems), rest_lengths is (1, n_elems)
+ // dilatation_rate is (1, n_elems)
+ const IndexType n_nodes = position.cols();
+ const IndexType n_elems = lengths.cols();
+
+ // Compute r_dot_v = batch_dot(position, velocity)
+ // This is the dot product of position and velocity at each node
+ // MatrixType r_dot_v(1, n_nodes);
+ auto&& r_dot_v = block.template get();
+ #pragma omp parallel for simd schedule(static)
+ for (IndexType k = 0; k < n_nodes; ++k) {
+ r_dot_v(0, k) = position(0, k) * velocity(0, k) +
+ position(1, k) * velocity(1, k) +
+ position(2, k) * velocity(2, k);
+ }
+
+ // Compute r_plus_one_dot_v = batch_dot(position[..., 1:], velocity[..., :-1])
+ // Dot product of position[1:] and velocity[:-1] (both have n_elems elements)
+ // MatrixType r_plus_one_dot_v(1, n_elems);
+ auto&& r_plus_one_dot_v = block.template get();
+ #pragma omp parallel for simd schedule(static)
+ for (IndexType k = 0; k < n_elems; ++k) {
+ r_plus_one_dot_v(0, k) = position(0, k + 1) * velocity(0, k) +
+ position(1, k + 1) * velocity(1, k) +
+ position(2, k + 1) * velocity(2, k);
+ }
+
+ // Compute r_dot_v_plus_one = batch_dot(position[..., :-1], velocity[..., 1:])
+ // Dot product of position[:-1] and velocity[1:] (both have n_elems elements)
+ // MatrixType r_dot_v_plus_one(1, n_elems);
+ auto&& r_dot_v_plus_one = block.template get();
+ #pragma omp parallel for simd schedule(static)
+ for (IndexType k = 0; k < n_elems; ++k) {
+ r_dot_v_plus_one(0, k) = position(0, k) * velocity(0, k + 1) +
+ position(1, k) * velocity(1, k + 1) +
+ position(2, k) * velocity(2, k + 1);
+ }
+
+ // Compute dilatation_rate for each element
+ // dilatation_rate[k] = (r_dot_v[k] + r_dot_v[k + 1] - r_dot_v_plus_one[k] - r_plus_one_dot_v[k])
+ // / lengths[k] / rest_lengths[k]
+ #pragma omp parallel for simd schedule(static)
+ for (IndexType k = 0; k < n_elems; ++k) {
+ dilatation_rate(0, k) = (r_dot_v(0, k) + r_dot_v(0, k + 1) -
+ r_dot_v_plus_one(0, k) - r_plus_one_dot_v(0, k)) /
+ lengths(0, k) / rest_lengths(0, k);
+ }
+}
+
+// Compute internal torques for Cosserat rod
+// This is a template function that will be implemented later
+template
+inline void compute_internal_torques(BlockType& block) {
+ // Get variable views
+ auto&& voronoi_dilatation = block.template get();
+ auto&& internal_couple = block.template get();
+ auto&& kappa = block.template get();
+ auto&& rest_voronoi_lengths = block.template get();
+ auto&& director = block.template get();
+ auto&& tangents = block.template get();
+ auto&& internal_stress = block.template get();
+ auto&& rest_lengths = block.template get();
+ auto&& mass_second_moment_of_inertia = block.template get();
+ auto&& omega = block.template get();
+ auto&& dilatation = block.template get();
+ auto&& dilatation_rate = block.template get();
+ auto&& internal_torques = block.template get();
+
+ // voronoi_dilatation is (1, n_voronoi), internal_couple is (3, n_voronoi)
+ // kappa is (3, n_voronoi), rest_voronoi_lengths is (1, n_voronoi)
+ // director is (9, n_elems), tangents is (3, n_elems)
+ // internal_stress is (3, n_elems), rest_lengths is (1, n_elems)
+ // mass_second_moment_of_inertia is (9, n_elems), omega is (3, n_elems)
+ // dilatation is (1, n_elems), dilatation_rate is (1, n_elems)
+ // internal_torques is (3, n_elems)
+ const IndexType n_voronoi = voronoi_dilatation.cols();
+ const IndexType n_elems = internal_torques.cols();
+ const IndexType n_nodes = n_elems + 1;
+
+ // Scratch buffers
+ auto&& voronoi_dilatation_inv_cube_cached = block.template get();
+ auto&& bend_twist_couple_2D = block.template get();
+ auto&& J_omega_upon_e = block.template get();
+ auto&& scratch_vec_b_voronoi = block.template get();
+ auto&& bend_twist_couple_3D = block.template get();
+ auto&& shear_stretch_couple = block.template get();
+ auto&& lagrangian_transport = block.template get();
+ auto&& unsteady_dilatation = block.template get();
+
+ // Compute voronoi_dilatation_inv_cube_cached = 1.0 / voronoi_dilatation^3
+ // MatrixType voronoi_dilatation_inv_cube_cached(1, n_voronoi);
+ #pragma omp parallel for simd schedule(static)
+ for (IndexType k = 0; k < n_voronoi; ++k) {
+ double voronoi_dil = voronoi_dilatation(0, k);
+ voronoi_dilatation_inv_cube_cached(0, k) = 1.0 / (voronoi_dil * voronoi_dil * voronoi_dil);
+ }
+
+ // Compute bend_twist_couple_2D = difference_kernel(internal_couple * voronoi_dilatation_inv_cube_cached, ghost_voronoi_idx)
+ // First compute the product
+ // MatrixType internal_couple_scaled(3, n_voronoi);
+ auto internal_couple_scaled = scratch_vec_b_voronoi;
+ for (IndexType i = 0; i < 3; ++i) {
+ #pragma omp parallel for simd schedule(static)
+ for (IndexType k = 0; k < n_voronoi; ++k) {
+ internal_couple_scaled(i, k) = internal_couple(i, k) * voronoi_dilatation_inv_cube_cached(0, k);
+ }
+ }
+
+ // Reset ghost values for internal_couple_scaled (OnVoronoi)
+ // reset_ghost :: internal_couple_scaled
+ block.template reset_ghost_for_variable();
+
+ // Apply difference_kernel (two_point_difference_kernel)
+ // MatrixType bend_twist_couple_2D(3, n_nodes);
+ two_point_difference_kernel(bend_twist_couple_2D, internal_couple_scaled);
+
+ // Compute bend_twist_couple_3D = quadrature_kernel((kappa x internal_couple) * rest_voronoi_lengths * voronoi_dilatation_inv_cube_cached, ghost_voronoi_idx)
+ // First compute kappa x internal_couple
+ // MatrixType kappa_cross_internal_couple(3, n_voronoi);
+ auto kappa_cross_internal_couple = scratch_vec_b_voronoi;
+ batch_cross(kappa_cross_internal_couple, kappa, internal_couple);
+
+ // Multiply by rest_voronoi_lengths and voronoi_dilatation_inv_cube_cached
+ // MatrixType bend_twist_couple_3D_input(3, n_voronoi);
+ auto bend_twist_couple_3D_input = scratch_vec_b_voronoi;
+ for (IndexType i = 0; i < 3; ++i) {
+ #pragma omp parallel for simd schedule(static)
+ for (IndexType k = 0; k < n_voronoi; ++k) {
+ bend_twist_couple_3D_input(i, k) = kappa_cross_internal_couple(i, k) *
+ rest_voronoi_lengths(0, k) *
+ voronoi_dilatation_inv_cube_cached(0, k);
+ }
+ }
+
+ // Reset ghost values
+ // reset_ghost:: bend_twist_couple_3D_input;
+ block.template reset_ghost_for_variable();
+
+ // Apply quadrature_kernel (trapezoidal)
+ // MatrixType bend_twist_couple_3D(3, n_nodes);
+ quadrature_kernel(bend_twist_couple_3D, bend_twist_couple_3D_input);
+
+ // Compute shear_stretch_couple = (Q^T * tangents) x internal_stress * rest_lengths
+ // First compute Q^T * tangents (same as director^T @ tangents)
+ // { // DEPRECATED BLOCK: REPLACED WITH BELOW OPS
+ // MatrixType director_tangents(3, n_elems);
+ // for (IndexType i = 0; i < 3; ++i) {
+ // #pragma omp parallel for simd schedule(static)
+ // for (IndexType k = 0; k < n_elems; ++k) {
+ // double sum = 0.0;
+ // for (IndexType j = 0; j < 3; ++j) {
+ // IndexType director_idx = j * 3 + i;
+ // sum += director(director_idx, k) * tangents(j, k);
+ // }
+ // director_tangents(i, k) = sum;
+ // }
+ // }
+
+ // // Compute cross product: (Q^T * tangents) x internal_stress
+ // // MatrixType shear_stretch_couple(3, n_elems);
+ // batch_cross(shear_stretch_couple, director_tangents, internal_stress);
+
+ // // Multiply by rest_lengths
+ // for (IndexType i = 0; i < 3; ++i) {
+ // #pragma omp parallel for simd schedule(static)
+ // for (IndexType k = 0; k < n_elems; ++k) {
+ // shear_stretch_couple(i, k) *= rest_lengths(0, k);
+ // }
+ // }
+ // }
+
+ #pragma omp parallel for simd schedule(static)
+ for (IndexType k = 0; k < n_elems; ++k) {
+ // Compute Q^T * tangents
+
+ double a0 = director(0, k) * tangents(0, k) + director(3, k) * tangents(1, k) + director(6, k) * tangents(2, k);
+ double a1 = director(1, k) * tangents(0, k) + director(4, k) * tangents(1, k) + director(7, k) * tangents(2, k);
+ double a2 = director(2, k) * tangents(0, k) + director(5, k) * tangents(1, k) + director(8, k) * tangents(2, k);
+
+ // internal_stress alias
+ double i0 = internal_stress(0, k);
+ double i1 = internal_stress(1, k);
+ double i2 = internal_stress(2, k);
+
+ // cross-product (Q*t x n) * l_hat
+ shear_stretch_couple(0, k) += (a1 * i2 - a2 * i1) * rest_lengths(0, k);
+ shear_stretch_couple(1, k) += (a2 * i0 - a0 * i2) * rest_lengths(0, k);
+ shear_stretch_couple(2, k) += (a0 * i1 - a1 * i0) * rest_lengths(0, k);
+ }
+
+ // Compute J_omega_upon_e = batch_matvec(mass_second_moment_of_inertia, omega) / dilatation
+ // MatrixType J_omega_upon_e(3, n_elems);
+ #pragma omp parallel for simd schedule(static)
+ for (IndexType k = 0; k < n_elems; ++k) {
+ // Extract 3x3 mass_second_moment_of_inertia for element k
+ double m00 = mass_second_moment_of_inertia(0, k);
+ double m10 = mass_second_moment_of_inertia(1, k);
+ double m20 = mass_second_moment_of_inertia(2, k);
+ double m01 = mass_second_moment_of_inertia(3, k);
+ double m11 = mass_second_moment_of_inertia(4, k);
+ double m21 = mass_second_moment_of_inertia(5, k);
+ double m02 = mass_second_moment_of_inertia(6, k);
+ double m12 = mass_second_moment_of_inertia(7, k);
+ double m22 = mass_second_moment_of_inertia(8, k);
+
+ // Compute mass_second_moment_of_inertia @ omega
+ J_omega_upon_e(0, k) = (m00 * omega(0, k) + m01 * omega(1, k) + m02 * omega(2, k)) / dilatation(0, k);
+ J_omega_upon_e(1, k) = (m10 * omega(0, k) + m11 * omega(1, k) + m12 * omega(2, k)) / dilatation(0, k);
+ J_omega_upon_e(2, k) = (m20 * omega(0, k) + m21 * omega(1, k) + m22 * omega(2, k)) / dilatation(0, k);
+ }
+
+ // Compute lagrangian_transport = (J * omega / dilatation) x omega
+ // MatrixType lagrangian_transport(3, n_elems);
+ batch_cross(lagrangian_transport, J_omega_upon_e, omega);
+
+ // Compute unsteady_dilatation = J_omega_upon_e * dilatation_rate / dilatation
+ // MatrixType unsteady_dilatation(3, n_elems);
+ for (IndexType i = 0; i < 3; ++i) {
+ #pragma omp parallel for simd schedule(static)
+ for (IndexType k = 0; k < n_elems; ++k) {
+ unsteady_dilatation(i, k) = J_omega_upon_e(i, k) * dilatation_rate(0, k) / dilatation(0, k);
+ }
+ }
+
+ // Compute internal_torques = sum of all components
+ // Note: bend_twist_couple_2D and bend_twist_couple_3D are (3, n_nodes), but we need (3, n_elems)
+ // So we need to take the element values (columns 0 to n_elems-1)
+ for (IndexType i = 0; i < 3; ++i) {
+ #pragma omp parallel for simd schedule(static)
+ for (IndexType k = 0; k < n_elems; ++k) {
+ internal_torques(i, k) += // bend_twist_couple_2D(i, k) + // Note: internal_torques share bend_twist_couple_2D
+ bend_twist_couple_3D(i, k) +
+ shear_stretch_couple(i, k) +
+ lagrangian_transport(i, k) +
+ unsteady_dilatation(i, k);
+ }
+ }
+}
+
+} // namespace elasticapp
diff --git a/backend/src/cosserat_rod_system.h b/backend/src/cosserat_rod_system.h
new file mode 100644
index 000000000..11d6875ed
--- /dev/null
+++ b/backend/src/cosserat_rod_system.h
@@ -0,0 +1,181 @@
+#pragma once
+
+#include "system.h"
+
+namespace elasticapp {
+
+// CosseratRod-specific variable tags
+// These variable types are now made internal to this translation unit
+namespace system::cosserat_rod {
+ // Node variables
+ struct Position : Placement::OnNode, DataType::Vector {
+ static constexpr std::string_view name = "position";
+ };
+ struct Velocity : Placement::OnNode, DataType::Vector {
+ static constexpr std::string_view name = "velocity";
+ };
+ struct Acceleration : Placement::OnNode, DataType::Vector {
+ static constexpr std::string_view name = "acceleration";
+ };
+ struct Mass : Placement::OnNode, DataType::Scalar {
+ static constexpr std::string_view name = "mass";
+ inline static MatrixType ghost_value = MatrixType::Constant(1, 1, 1.0);
+ };
+ struct InternalForces : Placement::OnNode, DataType::Vector {
+ static constexpr std::string_view name = "internal_forces";
+ };
+ struct ExternalForces : Placement::OnNode, DataType::Vector {
+ static constexpr std::string_view name = "external_forces";
+ };
+
+ // Element variables
+ struct Director : Placement::OnElement, DataType::Matrix {
+ static constexpr std::string_view name = "director";
+ };
+ struct Omega : Placement::OnElement, DataType::Vector {
+ static constexpr std::string_view name = "omega";
+ };
+ struct Alpha : Placement::OnElement, DataType::Vector {
+ static constexpr std::string_view name = "alpha";
+ };
+ struct RestLengths : Placement::OnElement, DataType::Scalar {
+ static constexpr std::string_view name = "rest_lengths";
+ // Note: Not constexpr because Eigen dynamic matrices don't have constexpr constructors
+ // For dynamic matrices, Constant requires dimensions: Constant(rows, cols, value)
+ inline static MatrixType ghost_value = MatrixType::Constant(1, 1, 1.0);
+ };
+ struct Density : Placement::OnElement, DataType::Scalar {
+ static constexpr std::string_view name = "density";
+ };
+ struct Volume : Placement::OnElement, DataType::Scalar {
+ static constexpr std::string_view name = "volume";
+ };
+ struct MassSecondMomentOfInertia : Placement::OnElement, DataType::Matrix {
+ static constexpr std::string_view name = "mass_second_moment_of_inertia";
+ };
+ struct InvMassSecondMomentOfInertia : Placement::OnElement, DataType::Matrix {
+ static constexpr std::string_view name = "inv_mass_second_moment_of_inertia";
+ };
+ struct InternalTorques : Placement::OnElement, DataType::Vector {
+ static constexpr std::string_view name = "internal_torques";
+ };
+ struct ExternalTorques : Placement::OnElement, DataType::Vector {
+ static constexpr std::string_view name = "external_torques";
+ };
+ struct Lengths : Placement::OnElement, DataType::Scalar {
+ static constexpr std::string_view name = "lengths";
+ };
+ struct Tangents : Placement::OnElement, DataType::Vector {
+ static constexpr std::string_view name = "tangents";
+ };
+ struct Radius : Placement::OnElement, DataType::Scalar {
+ static constexpr std::string_view name = "radius";
+ };
+ struct Dilatation : Placement::OnElement, DataType::Scalar {
+ static constexpr std::string_view name = "dilatation";
+ };
+ struct DilatationRate : Placement::OnElement, DataType::Scalar {
+ static constexpr std::string_view name = "dilatation_rate";
+ };
+ struct Sigma : Placement::OnElement, DataType::Vector {
+ static constexpr std::string_view name = "sigma";
+ };
+ struct RestSigma : Placement::OnElement, DataType::Vector {
+ static constexpr std::string_view name = "rest_sigma";
+ };
+ struct InternalStress : Placement::OnElement, DataType::Vector {
+ static constexpr std::string_view name = "internal_stress";
+ };
+ struct ShearMatrix : Placement::OnElement, DataType::Matrix {
+ static constexpr std::string_view name = "shear_matrix";
+ };
+
+ // Voronoi variables
+ struct RestVoronoiLengths : Placement::OnVoronoi, DataType::Scalar {
+ static constexpr std::string_view name = "rest_voronoi_lengths";
+ // Note: Not constexpr because Eigen dynamic matrices don't have constexpr constructors
+ // For dynamic matrices, Constant requires dimensions: Constant(rows, cols, value)
+ inline static MatrixType ghost_value = MatrixType::Constant(1, 1, 1.0);
+ };
+ struct VoronoiDilatation : Placement::OnVoronoi, DataType::Scalar {
+ static constexpr std::string_view name = "voronoi_dilatation";
+ };
+ struct Kappa : Placement::OnVoronoi, DataType::Vector {
+ static constexpr std::string_view name = "kappa";
+ };
+ struct RestKappa : Placement::OnVoronoi, DataType::Vector {
+ static constexpr std::string_view name = "rest_kappa";
+ };
+ struct InternalCouple : Placement::OnVoronoi, DataType::Vector {
+ static constexpr std::string_view name = "internal_couple";
+ };
+ struct BendMatrix : Placement::OnVoronoi, DataType::Matrix {
+ static constexpr std::string_view name = "bend_matrix";
+ };
+
+ // Dummy variables for computation
+ struct ScratchVectorA : Placement::OnElement, DataType::Vector { static constexpr std::string_view name = "scratch_vector_a"; };
+ struct ScratchVectorB : Placement::OnVoronoi, DataType::Vector { static constexpr std::string_view name = "scratch_vector_b"; };
+ struct ScratchVectorC : Placement::OnNode, DataType::Vector { static constexpr std::string_view name = "scratch_vector_c"; };
+ struct ScratchVectorD : Placement::OnNode, DataType::Vector { static constexpr std::string_view name = "scratch_vector_d"; };
+ struct ScratchVectorE : Placement::OnNode, DataType::Vector { static constexpr std::string_view name = "scratch_vector_e"; };
+ struct ScratchVectorF : Placement::OnNode, DataType::Vector { static constexpr std::string_view name = "scratch_vector_f"; };
+
+ struct ScratchScalarA : Placement::OnNode, DataType::Scalar { static constexpr std::string_view name = "scratch_scalar_a"; };
+ struct ScratchScalarB : Placement::OnNode, DataType::Scalar { static constexpr std::string_view name = "scratch_scalar_b"; };
+ struct ScratchScalarC : Placement::OnNode, DataType::Scalar { static constexpr std::string_view name = "scratch_scalar_c"; };
+}
+
+// CosseratRodSystem is a System with all variables from CosseratRod
+// Variables are organized by placement (Node, Element, Voronoi)
+using CosseratRodSystem = System<
+ // Node variables
+ system::cosserat_rod::Position,
+ system::cosserat_rod::Velocity,
+ system::cosserat_rod::Acceleration,
+ system::cosserat_rod::Mass,
+ system::cosserat_rod::InternalForces,
+ system::cosserat_rod::ExternalForces,
+
+ // Element variables
+ system::cosserat_rod::Omega,
+ system::cosserat_rod::Alpha,
+ system::cosserat_rod::Director,
+ system::cosserat_rod::RestLengths,
+ system::cosserat_rod::Density,
+ system::cosserat_rod::Volume,
+ system::cosserat_rod::MassSecondMomentOfInertia,
+ system::cosserat_rod::InvMassSecondMomentOfInertia,
+ system::cosserat_rod::InternalTorques,
+ system::cosserat_rod::ExternalTorques,
+ system::cosserat_rod::Lengths,
+ system::cosserat_rod::Tangents,
+ system::cosserat_rod::Radius,
+ system::cosserat_rod::Dilatation,
+ system::cosserat_rod::DilatationRate,
+ system::cosserat_rod::Sigma,
+ system::cosserat_rod::RestSigma,
+ system::cosserat_rod::InternalStress,
+ system::cosserat_rod::ShearMatrix,
+
+ // Voronoi variables
+ system::cosserat_rod::RestVoronoiLengths,
+ system::cosserat_rod::VoronoiDilatation,
+ system::cosserat_rod::Kappa,
+ system::cosserat_rod::RestKappa,
+ system::cosserat_rod::InternalCouple,
+ system::cosserat_rod::BendMatrix,
+
+ // Dummy variables
+ system::cosserat_rod::ScratchVectorA,
+ system::cosserat_rod::ScratchVectorB,
+ system::cosserat_rod::ScratchVectorC,
+ system::cosserat_rod::ScratchVectorD,
+ system::cosserat_rod::ScratchVectorE,
+ system::cosserat_rod::ScratchVectorF,
+ system::cosserat_rod::ScratchScalarA,
+ system::cosserat_rod::ScratchScalarB,
+ system::cosserat_rod::ScratchScalarC
+>;
+
+} // namespace elasticapp
diff --git a/backend/src/environment/_api.cpp b/backend/src/environment/_api.cpp
new file mode 100644
index 000000000..419df9813
--- /dev/null
+++ b/backend/src/environment/_api.cpp
@@ -0,0 +1,142 @@
+#include
+#include
+#include
+#include "../api.h"
+#include "api.h"
+#include "collision/collision_system.h"
+#include "collision/physics/no_interaction.h"
+#include "collision/physics/linear_spring_dashpot.h"
+
+namespace py = pybind11;
+
+PYBIND11_MODULE(_collision, m) {
+ m.doc() = R"pbdoc(
+ Elasticapp Collision module
+ ---------------------------
+
+ Provides collision detection and resolution for Cosserat rods using
+ the Discrete Element Method (DEM).
+ )pbdoc";
+
+ using namespace elasticapp;
+ using namespace elasticapp::environment::collision;
+ using namespace elasticapp::environment::collision::physics;
+
+ // NoInteraction model class (for testing)
+ py::class_(m, "NoInteraction")
+ .def(py::init<>(),
+ R"pbdoc(
+ Initialize a NoInteraction collision physics model.
+
+ This model returns zero force for all contacts, useful for testing
+ collision detection without applying forces.
+
+ Args:
+ None (no parameters required)
+ )pbdoc");
+
+ // LinearSpringDashpot model class
+ py::class_(m, "LinearSpringDashpot")
+ .def(py::init(),
+ R"pbdoc(
+ Initialize a LinearSpringDashpot collision physics model.
+
+ Args:
+ k_normal: Normal spring constant for repulsion force
+ eta_normal: Normal damping coefficient (also used for tangential damping)
+ friction: Static friction coefficient
+ )pbdoc",
+ py::arg("k_normal") = 1.0,
+ py::arg("eta_normal") = 0.1,
+ py::arg("friction") = 0.5)
+ .def(py::init(),
+ R"pbdoc(
+ Initialize a LinearSpringDashpot collision physics model with explicit tangential damping.
+
+ Args:
+ k_normal: Normal spring constant for repulsion force
+ eta_normal: Normal damping coefficient
+ eta_tangential: Tangential damping coefficient
+ friction: Static friction coefficient
+ )pbdoc",
+ py::arg("k_normal"),
+ py::arg("eta_normal"),
+ py::arg("eta_tangential"),
+ py::arg("friction"))
+ .def(py::init(),
+ R"pbdoc(
+ Initialize a LinearSpringDashpot collision physics model with explicit tangential spring and damping.
+
+ Args:
+ k_normal: Normal spring constant for repulsion force
+ eta_normal: Normal damping coefficient
+ k_tangential: Tangential spring constant
+ eta_tangential: Tangential damping coefficient
+ friction: Static friction coefficient
+ )pbdoc",
+ py::arg("k_normal"),
+ py::arg("eta_normal"),
+ py::arg("k_tangential"),
+ py::arg("eta_tangential"),
+ py::arg("friction"));
+
+ // CollisionSystem class (using DefaultCollisionSystem)
+ // Support both LinearSpringDashpot and NoInteraction models
+ py::class_(m, "CollisionSystem")
+ .def(py::init(),
+ R"pbdoc(
+ Initialize a CollisionSystem with a LinearSpringDashpot physics model.
+
+ Uses default policies: HashGrid (coarse), MaxContacts (fine), UnionFind (batching).
+
+ Args:
+ model: The LinearSpringDashpot collision physics model
+ detect_every: Perform coarse detection every N steps (default: 1, meaning every step)
+ )pbdoc",
+ py::arg("model"),
+ py::arg("detect_every") = 1)
+ .def(py::init(),
+ R"pbdoc(
+ Initialize a CollisionSystem with a NoInteraction physics model.
+
+ Uses default policies: HashGrid (coarse), MaxContacts (fine), UnionFind (batching).
+ Useful for testing collision detection without applying forces.
+
+ Args:
+ model: The NoInteraction collision physics model
+ detect_every: Perform coarse detection every N steps (default: 1, meaning every step)
+ )pbdoc",
+ py::arg("model"),
+ py::arg("detect_every") = 1)
+ .def("detect_every", &DefaultCollisionSystem::detect_every,
+ R"pbdoc(
+ Get the detect_every parameter.
+
+ Returns:
+ int: Number of steps between coarse detection calls
+ )pbdoc")
+ .def("set_detect_every", &DefaultCollisionSystem::set_detect_every,
+ R"pbdoc(
+ Set the detect_every parameter.
+
+ Args:
+ detect_every: Number of steps between coarse detection calls
+ )pbdoc",
+ py::arg("detect_every"))
+ .def("resolve", &DefaultCollisionSystem::resolve,
+ R"pbdoc(
+ Resolve collisions for a BlockRodSystem.
+
+ This method performs the full collision detection and resolution pipeline:
+ 1. Data extraction from Block
+ 2. Coarse collision detection (HashGrid)
+ 3. Fine collision detection (MaxContacts)
+ 4. Contact batching (UnionFind)
+ 5. Contact resolution (LinearSpringDashpot)
+ 6. Force application to ExternalForces
+
+ Args:
+ system: The BlockRodSystem to resolve collisions for
+ )pbdoc",
+ py::arg("system"));
+}
diff --git a/backend/src/environment/api.h b/backend/src/environment/api.h
new file mode 100644
index 000000000..24f356cf7
--- /dev/null
+++ b/backend/src/environment/api.h
@@ -0,0 +1,16 @@
+#pragma once
+
+#include "../api.h" // For BlockRodSystem
+#include "collision/collision_system.h"
+#include "collision/course_detection/hash_grid.h"
+#include "collision/fine_detection/max_contacts.h"
+#include "collision/batching/union_find.h"
+
+namespace elasticapp::environment {
+
+// Re-export CollisionSystem from collision namespace for convenience
+namespace collision {
+ using DefaultCollisionSystem = CollisionSystem;
+}
+
+} // namespace elasticapp::environment
diff --git a/backend/src/environment/collision/batching/union_find.cpp b/backend/src/environment/collision/batching/union_find.cpp
new file mode 100644
index 000000000..702cec507
--- /dev/null
+++ b/backend/src/environment/collision/batching/union_find.cpp
@@ -0,0 +1,105 @@
+#include "union_find.h"
+#include