Tensor Toolbox for Modern Fortran (ttb)
Commercial FEM software packages offer interfaces (user subroutines written in Fortran) for custom defined user materials like UMAT in Abaqus or HYPELA2 in MSC.Marc. In comparison to other scientific programming languages like MATLAB or Python Fortran is not as comfortable to use when dealing with high level programming features of tensor manipulation. On the other hand it's super fast  so why not combine the handy features from MATLAB or Python's NumPy/Scipy with the speed of Fortran? That's the reason why I started working on a simple but effective module called Tensor Toolbox for Modern Fortran. I adopted the idea to my needs from Naumann, C. (2016).
The full documentation is available at https://adtzlr.github.io/ttb. This project is licensed under the terms of the MIT license.
Overview
It provides the following basic operations for tensor calculus (all written in double precision real(kind=8)
):
 Dot Product
C(i,j) = A(i,k) B(k,j)
written asC = A*B
 Double Dot Product
C = A(i,j) B(i,j)
written asC = A**B
 Dyadic Product
C(i,j,k,l) = A(i,j) B(k,l)
written asC = A.dya.B
 Crossed Dyadic Product
C(i,j,k,l) = (A(i,k) B(j,l) + A(i,l) B(j,k))/2
written asC = A.cdya.B
 Addition
C(i,j) = A(i,j) + B(i,j)
written asC = A+B
 Subtraction
C(i,j) = A(i,j)  B(i,j)
written asC = AB
 Multiplication and Division by a Scalar
 Deviatoric Part of Tensor
dev(C) = C  tr(C)/3 * Eye
written asdev(C)
 Transpose and Permutation of indices written as
B = permute(A,1,3,2,4)
 Rank 2 Identity tensor of input type
Eye = identity2(Eye)
withC = Eye*C
 Rank 4 Identity tensor (symmetric variant) of input type
I4 = identity4(Eye)
orI4 = Eye.cdya.Eye
withC = I4**C
orinv(C) = identitiy4(inv(C))**C
 Square Root of a positive definite rank 2 tensor
U = sqrt(C)
 Natural logarithm (1) of a rank 2 tensor
lnC = ln(C)
 Exponential function (1) of a rank 2 tensor
expC = exp(C)
 Assigment of a realvalued Scalar to all components of a Tensor
A = 0.0
orA = 0.d0
 Assigment of a realvalued Array to a Tensor with matching dimensions
A = B
where B is an Array and A a Tensor  Assigment of a Tensor in Voigt notation to a Tensor in tensorial notation and vice versa
(1) series expansion of isotropic tensor function
The idea is to create derived data types for rank 1, rank 2 and rank 4 tensors (and it's symmetric variants). In a next step the operators are defined in a way that Fortran calls different functions based on the input types of the operator: performing a dot product between a vector and a rank 2 tensor or a rank 2 and a rank 2 tensor is a different function. Best of it: you don't have to take care of that.
Basic Usage
The most basic example on how to use this module is to download the module, put the 'ttb'Folder in your working directory and add two lines of code:
include 'ttb/ttb_library.f'
program script101_ttb
use Tensor
implicit none
! user code
end program script101_ttb
The include 'ttb/ttb_library.f'
statement replaces the line with the content of the ttbmodule. The first line in a program or subroutine is now a use Tensor
statement. That's it  now you're ready to go.
Tensor or Voigt Notation
It depends on your preferences: either you store all tensors in full tensor dimension(3,3)
or in voigt dimension(6)
notation. The equations remain (nearly) the same. Dot Product, Double Dot Product  every function is implemented in both full tensor and voigt notation. Look for the voigtcomments in an example of a user subroutine for MSC.Marc.
Access Tensor components by Array
Tensor components may be accessed by a conventional array with the name of the tensor variable T
followed by a percent operator %
and a keyword as follows:

Tensor of rank 1 components as array:
T%a
. ith component of T:T%a(i)

Tensor of rank 2 components as array:
T%ab
. i,j component of T:T%ab(i,j)

Tensor of rank 4 components as array:
T%abcd
. i,j,k,l component of T:T%abcd(i,j,k,l)

Symmetric Tensor of rank 2 (Voigt) components as array:
T%a6
. ith component of T:T%a6(i)

Symmetric Tensor of rank 4 (Voigt) components as array:
T%a6b6
. i,j component of T:T%a6b6(i,j)
(at least minor symmetric)
Warning: Output as array
It is not possible to access tensor components of a tensor valued function in a direct way s = symstore(S1)%a6
 unfortunately this is a limitation of Fortran. To avoid the creation of an extra variable it is possible to use the asarray(T,i_max[,j_max,k_max,l_max])
function to access tensor components. i_max,j_max,k_max,l_max
is not the single component, instead a slice T%abcd(1:i_max,1:j_max,1:k_max,1:l_max)
is returned. This can be useful when dealing with mixed formulation or variation principles where the last entry/entries of stress and strain voigt vectors are used for the pressure boundary. To export a full stress tensor S1
to voigt notation use:
s(1:ndim) = asarray( voigt(S1), ndim )
d(1:ndim,1:ndim) = asarray( voigt(C4), ndim, ndim )
Abaqus Users: Output as abqarray
To export a stress tensor to Abaqus Voigt notation use asabqarray
which reorders the storage indices to 11,22,33,12,13,23
. This function is available for Tensor2s
and Tensor4s
data types.
s(1:ndim) = asabqarray( symstore(S1), ndim )
ddsdde(1:ndim,1:ndim) = asabqarray( symstore(C4), ndim, ndim )
A note on the Permutation of Indices
The permutation function reorders indices in the given order for a fourth order tensor of data type Tensor4
. Example: (i,j,k,l) > (i,k,j,l)
with permute(C4,1,3,2,4)
.
External Libraries
This library is not using any LAPACK functions. Instead you can compile your subroutine with o3
flag where at least the Intel compiler translates nested forloops to Intel MPI functions. However commercial FEM packages have predefined optimization flags which should not or only be changed with caution. I'm open for ideas how to use LAPACK and Intel MPI in the future  please let me know.
NeoHookean Material
With the help of the Tensor module the Second PiolaKirchhoff stress tensor S
of a nearlyincompressible NeoHookean material model is basically a oneliner:
Second Piola Kirchhoff Stress Tensor
S = mu*det(C)**(1./3.)*dev(C)*inv(C)+p*det(C)**(1./2.)*inv(C)
While this is of course not the fastest way of calculating the stress tensor it is extremely short and readable. Also the second order tensor variables S, C
and scalar quantities mu, p
have to be created at the beginning of the program. A minimal working example for a very simple umat user subroutine can be found in script_umat.f. The program is just an example where umat is called and an output information is printed. It is shown that the tensor toolbox is only used inside the material user subroutine umat.
Material Elasticity Tensor
Isochoric part of the material elasticity tensor C4_iso
of a nearlyincompressible NeoHookean material model:
C4_iso = det(F)**(2./3.) * 2./3.* (
* tr(C) * identity4(inv(C))
*  (Eye.dya.inv(C))  (inv(C).dya.Eye)
* + tr(C)/3. * (inv(C).dya.inv(C)) )
Example of MSC.Marc HYPELA2
Here you can find an example of a nearlyincompressible version of a NeoHookean material for MSC.Marc. It works only in Total Lagrange (no push forward implemented). Updated Lagrange is implemented with a push forward of both stress tensor and tangent matrix. Herrmann Elements are automatically detected. As HYPELA2 is called twice per iteration the stiffness calculation is only active during stage lovl == 4
. One of the best things is the supersimple switch from tensor to voigt notation: Change data types of all symmetric tensors and save the right CauchyGreen deformation tensor in voigt notation. See commented lines for details.
Download HYPELA2: NeoHooke, MSC.Marc, Total Lagrange, Tensor Toolbox
Credits
Naumann, C.: Chemischmechanisch gekoppelte Modellierung und Simulation oxidativer Alterungsvorgänge in Gummibauteilen (German). PhD thesis. Fakultät für Maschinenbau der Technischen Universität Chemnitz, 2016.