Linear elasticity in C++

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 C++, you need to install a C++ compiler and the build tool cmake. Additionally, you need a library for linear algebra, and a library to generate Python bindings. We recommend Eigen for linear algebra and either pybind11 or nanobind for the Python bindings. You can install it into your current conda-environment by running the following command:

mamba install cmake, clang, pybind11, eigen -c conda-forge

After that you can create a new C++ project with the following commands:

mkdir elasticity_cpp
mkdir elasticity_cpp/src
touch elasticity_cpp/CMakeLists.txt
touch src/main.cpp

This will create a new folder called elasticity_cpp with the following structure:

elasticity_cpp
├── CMakeLists.txt
└── src
    └── main.cpp

In order to write a constitutive model which is then linked to Python, you need to define the dependencies in the CMakeLists.txt file. For linear algebra on small matrices, the Eigen library is used. For interfacing with Python, the pybind11 library is used:

cmake_minimum_required(VERSION 3.15...3.26)
project(elasticity_cpp LANGUAGES CXX)

set(PYBIND11_FINDPYTHON ON)
find_package(pybind11 CONFIG REQUIRED)
find_package(Eigen3 REQUIRED NO_MODULE)

pybind11_add_module(elasticity_cpp MODULE src/main.cpp)
#install(TARGETS _core DESTINATION ${SKBUILD_PROJECT_NAME})
target_link_libraries(elasticity_cpp PRIVATE pybind11::module Eigen3::Eigen)

Writing the model

An example for a full source code of an elasticity model in C++ is shown below:

#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
#include <pybind11/eigen.h>
#include <eigen3/Eigen/Core>
#include <iostream>

inline Eigen::Vector<double, 6> strain_from_grad_u(const Eigen::Matrix3d &grad_u)
{
    constexpr double FACTOR = 0.707106781186547524400844362104849039;
    Eigen::Vector<double, 6> strain;
    strain(0) = grad_u(0, 0);
    strain(1) = grad_u(1, 1);
    strain(2) = grad_u(2, 2);
    strain(3) = FACTOR * (grad_u(0, 1) + grad_u(1, 0));
    strain(4) = FACTOR * (grad_u(0, 2) + grad_u(2, 0));
    strain(5) = FACTOR * (grad_u(1, 2) + grad_u(2, 1));
    return strain;
}

struct Elasticity3D
{
    Eigen::Matrix<double, 6, 6> D_;
    Elasticity3D(double E, double nu)
    {
        double lambda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu));
        double mu = E / (2.0 * (1.0 + nu));
        D_ << lambda + 2.0 * mu, lambda, lambda, 0.0, 0.0, 0.0,
            lambda, lambda + 2.0 * mu, lambda, 0.0, 0.0, 0.0,
            lambda, lambda, lambda + 2.0 * mu, 0.0, 0.0, 0.0,
            0.0, 0.0, 0.0, 2.0 * mu, 0.0, 0.0,
            0.0, 0.0, 0.0, 0.0, 2.0 * mu, 0.0,
            0.0, 0.0, 0.0, 0.0, 0.0, 2.0 * mu;
    }

    void evaluate(
        double t,
        double del_t,
        const Eigen::Ref<const Eigen::VectorXd> &del_grad_u,
        Eigen::Ref<Eigen::VectorXd> &stress,
        Eigen::Ref<Eigen::VectorXd> &tangent,
        std::map<std::string, Eigen::Ref<Eigen::VectorXd>> &history)
    {
        int n_ip = del_grad_u.size() / 9;
        for (int ip = 0; ip < n_ip; ip++)
        {
            Eigen::Vector<double, 9> grad_u_flat = del_grad_u.segment<9>(9 * ip);
            stress.segment<6>(6 * ip) += D_ * strain_from_grad_u(Eigen::Map<Eigen::Matrix3d>(grad_u_flat.data()));
            tangent.segment<36>(36 * ip) = Eigen::Map<Eigen::Vector<double, 36>>(D_.data());
        }
    }

    std::map<std::string, int> history_dim() { 
        return {};
    }

};

namespace py = pybind11;

PYBIND11_MODULE(elasticity_cpp, m)
{
    m.doc() = "";
    py::class_<Elasticity3D> elasticity_3d(m, "Elasticity3D");
    elasticity_3d.def(py::init<double,double>(), py::arg("E"), py::arg("nu"));
    elasticity_3d.def("evaluate", &Elasticity3D::evaluate);
    elasticity_3d.def("history_dim", &Elasticity3D::history_dim);
}

Compilation

To compile the C++ code, you need to create a build directory and run CMake:

cmake -DCMAKE_CXX_COMPILER=clang -S src/main.cpp -B build
cmake --build build

This creates the shared library build/elasticity_cpp.platform_tag.so where platform_tag is the platform specific tag of your system.

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:

import elastictity_cpp

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 C++ 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 C++ extension. The full source code of the model is shown below:

from elastictity_cpp 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