Linear elasticity in Rust
All files for this example can be found in the GitHub repository.
Setting up the project
In order to build a linear elasticity model in Rust, you need to install the rust compiler which comes with the build tool Cargo. You can install it into your current conda-environment by running the following command:
After that you can create a new rust project using Cargo:
This will create a new folder called elasticity_rs with the following structure:
In order to write a constitutive model which is then linked to rust, you need to define the dependencies in the Cargo.toml file. For linear algebra on small matrices, the nalgebra crate is used. For interfacing with numpy.ndarray, the numpy crate is used. The PyO3 crate is used to create a Python module from the rust code. The dependencies are defined as follows:
[package]
name = "elasticity_rs"
version = "0.1.0"
edition = "2021"
# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
[lib]
name = "elasticity_rs"
crate-type = ["cdylib"]
[dependencies]
numpy = {version = "*", features=["nalgebra"]}
nalgebra = "*"
[dependencies.pyo3]
version = "*"
# "extension-module" tells pyo3 we want to build an extension module (skips linking against libpython.so)
# "abi3-py38" tells pyo3 (and maturin) to build using the stable ABI with minimum Python version 3.8
features = ["extension-module", "abi3-py38"]
Note that crate-type = ["cdylib"] is used to create a dynamic library which can be imported into Python.
Writing the model
An example for a full source code of an elasticity model in Rust is shown below:
use std::collections::HashMap;
use nalgebra::{Const, DVectorView, Dyn, SMatrix, SVector, SVectorView};
use numpy::{PyReadonlyArray1, PyReadonlyArray2, PyReadwriteArray1};
use pyo3::{exceptions::PyRuntimeError, prelude::*};
pub fn strain_from_grad_u(grad_u: &SMatrix<f64, 3, 3>) -> SVector<f64, 6> {
// creates the strain rate mandel vector directly from the velocity gradient L instead of the rate
// of deformation tensor D. Therefore, the factor is 1/sqrt(2) instead of sqrt(2)
const FACTOR: f64 = std::f64::consts::FRAC_1_SQRT_2;
SVector::<f64, 6>::new(
grad_u.m11,
grad_u.m22,
grad_u.m33,
FACTOR * (grad_u.m12 + grad_u.m21),
FACTOR * (grad_u.m13 + grad_u.m31),
FACTOR * (grad_u.m23 + grad_u.m32),
)
}
#[pyclass]
struct Elasticity3D {
D: SMatrix<f64, 6, 6>,
}
#[pymethods]
impl Elasticity3D {
#[new]
fn new(E: f64, nu: f64) -> PyResult<Self> {
let mut D = SMatrix::<f64, 6, 6>::zeros();
let lambda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu));
let mu = E / (2.0 * (1.0 + nu));
D[(0, 0)] = lambda + 2.0 * mu;
D[(0, 1)] = lambda;
D[(0, 2)] = lambda;
D[(1, 0)] = lambda;
D[(1, 1)] = lambda + 2.0 * mu;
D[(1, 2)] = lambda;
D[(2, 0)] = lambda;
D[(2, 1)] = lambda;
D[(2, 2)] = lambda + 2.0 * mu;
D[(3, 3)] = 2.0 * mu;
D[(4, 4)] = 2.0 * mu;
D[(5, 5)] = 2.0 * mu;
Ok(Self { D })
}
fn evaluate(
&self,
t: f64,
del_t: f64,
del_grad_u: PyReadonlyArray1<f64>,
stress: PyReadwriteArray1<f64>,
tangent: PyReadwriteArray1<f64>,
history: HashMap<String, PyReadwriteArray1<f64>>,
) -> PyResult<()> {
//convert numpy arrays to nalgebra matrices
let del_grad_u = del_grad_u
.try_as_matrix::<Dyn, Const<1>, Const<1>, Dyn>()
.unwrap();
let mut stress = stress
.try_as_matrix_mut::<Dyn, Const<1>, Const<1>, Dyn>()
.unwrap();
let mut tangent = tangent
.try_as_matrix_mut::<Dyn, Const<1>, Const<1>, Dyn>()
.unwrap();
let n_ip = del_grad_u.nrows() / 9;
for ip in 0..n_ip {
let mut view_stress = stress.fixed_view_mut::<6, 1>(ip * 6, 0);
let mut view_tangent = tangent.fixed_view_mut::<36, 1>(ip * 36, 0);
//create matrix from column slice. This is actually false because fenics is row-wise and
// nalgebra column-wise, but it makes no difference for the strain calculation
let del_grad_u_full = SMatrix::<f64, 3, 3>::from_column_slice(
del_grad_u.fixed_view::<9, 1>(ip * 9, 0).as_slice(),
);
view_stress += self.D * strain_from_grad_u(&del_grad_u_full);
view_tangent.copy_from_slice(&self.D.as_slice());
}
Ok(())
}
fn history_dim(&self) -> HashMap<String, i32> {
HashMap::new()
}
}
#[pymodule]
fn elasticity_rs(_py: Python, m: &PyModule) -> PyResult<()> {
m.add_class::<Elasticity3D>()?;
Ok(())
}
Compilation
The rust code can be compiled using Cargo:
This creates the shared library target/release/libelasticity_rs.so. This shared library cannot be imported directly, since the name of the shared library is not the same as the name of the module. To fix this, you can create a link to the shared library with the following command:
Importing the shared library into Python
If the created shared library is located within the same directory as your Python script, you can import it directly with the following code:
You may also create a soft or hard link to the shared library in your Python script directory, or import the script using the path to the shared library with
import importlib.util
from pathlib import Path
path_as_string = "..."
module_path = Path(path_as_string)
print(module_path.stem)
spec = importlib.util.spec_from_file_location(
module_path.stem.split(".")[0], module_path
)
module = importlib.util.module_from_spec(spec)
spec.loader.exec_module(module)
model = module.Elasticity3D(...)
Create an IncrSmallStrainModel
So far, we have created a class that can compute the stress and tangent stiffness for a given strain, however, since it is a Rust extension, Python cannot recognize it as a IncrSmallStrainModel. This is important because the IncrSmallStrainProblem requires the method constraint to determine the modeling assumptions about the strains and stresses that the model is implemented in. To do this, we need to create a Python class that inherits from IncrSmallStrainModel and calls the Rust extension. The full source code of the model is shown below:
from elastictity_rs import Elasticity3D
from fenics_constitutive import IncrSmallStrainModel, StressStrainConstraint
class PyElasticity3D(IncrSmallStrainModel):
def __init__(self, E, nu):
self._model = Elasticity3D(E, nu)
def evaluate(
self,
t: float,
del_t: float,
grad_del_u: np.ndarray,
stress: np.ndarray,
tangent: np.ndarray,
history: dict[str, np.ndarray] | None,
) -> None:
self._model.evaluate(t, del_t, grad_del_u, stress, tangent, history)
@property
def history_dim(self):
return self._model.history_dim()
@property
def constraint(self):
return StressStrainConstraint.FULL