Motivation

Following theoretical arguments in formal textbooks is hard. Engineers are practical people and when dealing with a new subject I know that the best way for me to learn is understand enough to be able to try out myself and experiment with the math. And there is no better way to experiment than implement the equations is a simple script and continue elaborate on that script till you realize that it’s 5am in the morning and you didn’t even had dinner. Bear with me, the code I want to talk about is the product of such a full immersion!

The theory of plasticity

The first time I had to approach a FEM simulation with a plastic material model I was quite startled. I had never done such a thing and had no idea why it took so long for the computer to tell me every setting was wrong and the results were meaningless. And more importantly, how to get this thing converging??? What is to converge exactly???

Of course I had classes on non-linear equations and the Newton-Raphson method during college, but those thing don’t really stick into my mind till I try out myself. And that’s how I end up studying “Computational Inelasticity” (Simo Hughes 1998) and “Computational Methods for Plasticity” (de Souza Neto Peri Owen 2008). After studing those… I still had just a vague idea of what was going on, I needed to implement the stuff myself!

The jax library

There are a lot of derivatives in mathematical models. And even more derivatives in numerical solutions of mathematical models. This is actually a pretty nasty problems for people developing computer programs from general algorithms written in textbooks cause usually the task of actually doing the math is left to the readed… not ideal. Fortunately, together with global pandemics and climate change the 21st century provided us with automatic differentiation. The power of having derivatives computed by the computer with reasonable accuracy just by asking for it. Even better there is a python library that perform derivatives of python functions, the JAX library, a stunning amazing tool which I made sure to overuse.

My code in 3 classes

I settled down one specific goal with this script: set an equilibrium problem with non-linearities in the material model, which has to be a plastic one, and solve this problem on a 3D cubic domain constituted by just one finite elemen, hence the name: FEM1elem.

For all the derivatives to be computed in the script I used the functionalities of the JAX library, from the computation of the derivatives of the shape function to the computation of the tangent stiffness matrixes for the famous Newton-Raphson method.

The code is divided in 3 python classes:

• PDE_problem takes as inputs a Material object and a Finite_element object (defined using the classes below) together with mesh, bare conditions, initial conditions and prescribed displacements matrices defined at each node, and an initial state variables matrix with values prescribed at the integration points.
  • Method Solutor_NL_statics_displ_control this method elaborates the data structures provided as inputs and sets a Newton-Raphson iteration algorithm to solve the non-linear problem of global equilibrium for a continuum solid. A global residual is calculated leveraging the Material and Finite_element methods and derivatives of this global residual function are calculated with the automatic differentiation method provided by the Jax library upon vectorization of the input objects.
• Finite_Element class contains all the methods needed to perform computations with the finite element method and the relative numerical integration procedures with the Gauss method.
  • Method shape_function shape function value for a specific node (identified with an identification number) calculated at the position provided in natural coordinates of the element.
  • Method coord_trans transformation from natural coordinates to real ones considering the element real coordinates.
  • Method Jacobian_matr Jacobian matrix calculated transforming the coord_trans method with the Jax jacfwd method.
  • Method B_strain_disp_matr strain-displacement matrix calculated using the inverse of the Jacobian of the coordinate transformation, and the derivative of the shape function (obtained leveraging the Jax automatic differentiation method).
  • Method Gauss_integration_int function integration using the Gauss method. The method also contains the function definitions for the internal force vector integrand. The implementation of this method contains some logical switches which are used to speed up the computation avoiding unnecessary calculations in the computations of the global residual and its derivatives.
• Material class takes as input the functional definition of an elastic potential, a yield function and an hardening law together with the respective parameters (whose number is variable to account for different formulations).
  • Method Voigt_stress calculation of stress increment from the strain increment input. This method contains the return mapping algorithm and the update of the state variables for the material point.

That’s it! (Considering the amount of time I spent writing it) I hope that the code could be useful to anyone interested in understanding the FEM method in its inner guts. The code is available at my Github toghether with the instructions needed to run it.