Linear elasticity as an Abaqus UMAT
All files for this example can be found in the GitHub repository.
Setting up the project
In the preceding examples, we have implemented a linear elasticity model using one language and looped over all quadrature points. However, an Abaqus UMAT is only written for one quadrature point. Therefore, we need to write code implementing a for loop in which the UMAT is called for each quadrature point. This is done in a C++ class that is then compiled into a shared library. The shared library is then imported into Python as in the previous examples.
In addition to the dependencies from the C++ example, we need to install a Fortran compiler. You can install it into your current conda-environment by running the following command:
After that you can create a new C++ project with the following commands:
mkdir umat
mkdir umat/src
touch umat/CMakeLists.txt
touch src/main.cpp
touch src/umat_linear_elastic.f
This will create a new folder called umat with the following structure:
In order to write a constitutive model with 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(umat LANGUAGES CXX)
set(PYBIND11_FINDPYTHON ON)
find_package(pybind11 CONFIG REQUIRED)
find_package(Eigen3 REQUIRED NO_MODULE)
pybind11_add_module(umat MODULE src/main.cpp)
#install(TARGETS _core DESTINATION ${SKBUILD_PROJECT_NAME})
target_link_libraries(umat PRIVATE pybind11::module Eigen3::Eigen)
Writing the model
An example for a full source code of an elasticity model as an Abaqus UMAT is shown below:
And the C++ code that calls the UMAT is shown below:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 | |
Compilation
To compile the C++ code, you need to create a build directory and run CMake:
This creates the shared library build/umat.platform_tag.so where platform_tag is the platform specific tag of your system.
The Fortran code can be compiled using 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 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