Perform a polar decomposition of a positive definite 3x3 matrix to obtain the stretch tensor which satisfies C = U^2

depending on the axis of rotation.

See https://en.wikipedia.org/wiki/Rotation_matrix and

Add an option to print tensor components with a simple function call.

Hi @adtzlr,
Thanks for this, it is going to be very useful!

In order to make your library easier to interface with existing codes, would you consider to do a CMakeLists.txt or a Makefile?

Does the lib support 2D calculations?

Add an assertion to both functions in libpower.f that throws an error if the exponent is neither 0, -1, nor a positive integer. This should catch "pow_2(tensor,-2)". If somebody tries to be smart and calls pow_2 with exponent 1/2 to get the square root, we still might get some problems. Even though pow_2(tensor,0.5) is caught by the compiler, pow_2(stress,1/2) is regarded as zero exponent. It is already explained what Power does in the documentation https://adtzlr.github.io/ttb/functions/power.html, but maybe we can add some further notes (iis integer) and protection (i>=-1).

Originally posted by @jfriedlein in #35 (comment) (slightly modified by @adtzlr)

in libassignscalar.f,

       subroutine assignscalar_4(T,w)
        implicit none
        
        type(Tensor4), intent(inout) :: T
        real(kind=8), intent(in) :: w
        
        T%abcd = T%abcd*w
        
       end subroutine assignscalar_4

the assignment must be replaced by

        T%abcd = w

for all assignments.

What is the best way to perform Tensor (3,3,3,3) to Voigt (6,6) contraction?

A symmetric fourth order identity tensor:
Eye4(i,j,k,l) = 1/2 * (Eye(i,k)*Eye(j,l)+Eye(i,l)*Eye(j,k))

is written as Eye4 = identity4(Eye) or Eye = Eye.cdya.Eye - but both methods result in different tensors!
A "true" fourth order identity tensor should be a zero matrix with diagonal elements equal to 1.0:
Eye4 = 0.0 and Eye4(i,i) = 1.0
but only the hardcoded identity4 matrix is true. I think I have to include a sum over the first dual shear indices (4) = (1,2)+(2,1), (5) = (2,3) + (3,1), ...

Hello, I am very interested in your fortran tensor calculation! Are there any known unfixed bugs in your code? Thank you for your generous answer!

with more examples for the functions and also replace some old svg-equations with mathjax.

Pages to be edited

• https://adtzlr.github.io/ttb/api/functions/asarray.html
• https://adtzlr.github.io/ttb/api/functions/astensor.html
• https://adtzlr.github.io/ttb/api/functions/asvoigt.html
• https://adtzlr.github.io/ttb/api/functions/inverse.html
• https://adtzlr.github.io/ttb/api/functions/unimodular.html

Add isotropic tensor functions like sqrt(T), exp(T), etc. Source

Or Eigenvalue/Eigenvector Calculation with the help of LAPACK

First error it gives at the line 'implicit none'. If I comment this line the code does not pop any error in Abaqus, but it gives NaN stress which should not be the case.

It also gives this error in Log file.
"Begin Linking Abaqus/Standard User Subroutines
Creating library standardU.lib and object standardU.exp
libirc.lib(fast_mem_ops.c.obj) : warning LNK4210: .CRT section exists; there may be unhandled static initializers or terminators
End Linking Abaqus/Standard User Subroutines"

What could be the problem?

Thank you for the valuable tensor library. Is Einstein summation included in it?

Discussed in #53

@Cyan2077 to be sure: you mean you have a deformation gradient $\boldsymbol{F}$ at hand and you would like to perform $\boldsymbol{F} = \boldsymbol{R} \boldsymbol{U}$, where $\boldsymbol{R}$ is the rotation matrix and $\boldsymbol{U}$ is the right stretch tensor? Are you interested in the left-polar decomposition $\boldsymbol{F} = \boldsymbol{V} \boldsymbol{R}$ too?

According to the priority of operators in Fortran it is necessary to use brackets for Tensor operators. For example the power operator ** which is overridden by the double dot product .ddot.has the highest priority. This is wrong, especially for a double contraction where one tensor is a result of a crossed dyadic product: A**B.cdya.C - the crossdyadic product has to be evaluated first. So the correct way of defining the equation is: A**(B.cdya.C)

Improve this behaviour...

as these functions are equal.

@adtzlr
I currently implement user material models in LS-Dyna with Fortran. Luckily, I found your tensor toolbox to be able to implement our tensor-based models.
Firstly, thank you so much for giving us such a well-written toolbox with, most importantly, a great documentation.
I found some minor issues when trying to use your toolbox in the Fortran version of LS-Dyna, as described below. Nevertheless, I use it now successfully. I also added some additional functions (volumetric part, deviatoric fourth order tensor, ...).
I was wondering whether I could contribute to your repository, so we get a toolbox for Abaqus, Marc and LS-Dyna.

ttb in LS-Dyna:
LS-Dyna probably uses an older Fortran version than Abaqus and Marc. When I include the toolbox, the following error results

.\ttb/libdiv.f(57): error #5286: Ambiguous generic interface OPERATOR(/): previously declared specific procedure DIV_10 is not distinguishable from this declaration. [DIV_10_R4] function div_10_r4(T, w)
... (for all functions ending with r4)

By simply deleting all the module procedure *r4 in ttb_library.f all errors disappear and we can use the toolbox (of course without the r4 variants). Is this because this fortran version cannot yet distinguish between data type-"templated" functions?

proposols:
Maybe we use a new ttb_library.f with the above modification or can incorporate this into a single file. The additional functions (vol, I_dev, ...) could be added to existing and new files.

If you can make time and give me an answer, I would appreciate it.

Finally it became clear what the problem in Issue #10 was, where in LS-Dyna the toolbox gave compiler errors .\ttb/libdiv.f(57): error #5286: Ambiguous generic interface OPERATOR(/): previously declared specific procedure DIV_10 is not distinguishable from this declaration. [DIV_10_R4] function div_10_r4(T, w) ... for all functions ending with r4.

LS-Dyna settings:
In the LS-Dyna makefile, the settings for Windows are:
"SMPD = -DOPENMP -Qopenmp -DAUTODOUBLE -4R8 -4I8"
and similarly in the Linux Makefile
"FF= ... -r8 ..."
"R8"/"r8" imply that by default double precision is used. This means that for instance whenever we write "real, ..." it is interpreted as "real(kind=8), ...". In Abaqus the default seems to be single precision, so "real, ..." is read as "real(kind=4), ...". I only noticed these compilation flags, when moving to Linux, so sorry for not mentioning this in the initial Issue #10.

Consequences for ttb:
In the toolbox all functions ending with "r4" use "real, ...", whereas in the remaining functions "real(kind=8), ..." is specified. If the default precision is r4, then there is a difference between "real, ...", which is read as "real(kind=4)" in the r4-functions, and "real(kind=8), ..." in the other functions. However, when the default precision is r8, then the compiler cannot distinguish between a function with "real, ..." (read as "real(kind=8), ...") and the actual "real(kind=8), ..." functions, thus throwing the above error.

Easy fix:
The above errors can be resolved by simply always specifying what precision shall be used and never relying on the default value. Hence, I added "(kind=4)" to all unspecified "real, ..." (as done in the appended four files based on current release v1.0.2, replaced all "real," by "real(kind=4),"), which eliminated the error. If this also works under Abaqus etc., we could directly return to the old framework (no precompilation flags, classical "include", lower case "ttb.f", ...).

libassignarray.f.txt
libassignscalar.f.txt
libdiv.f.txt
libdot.f.txt

While trying the new branch, I again noticed a change I did a while ago to my local implementation.

When defining the tensor types, I initialised all tensor components to zero, so that an uninitialised tensor is zero and its components are not undefined (numerically chaotic values). I was used to this from deal.ii and it might also protect other naive and innocent people like me. To initialise the components to zero, simply add "= 0." as shown below in the file ttb_library.f:

    type Tensor1
    ! tensor of rank 1
    real(kind=8), dimension(3) :: a = 0.
   end type Tensor1
   
   type Tensor2
    ! tensor of rank 2
    real(kind=8), dimension(3,3) :: ab = 0.
   end type Tensor2
   
   type Tensor2S
    ! symmetric tensor of rank 2 stored as vector
    real(kind=8), dimension(6) :: a6 = 0.
   end type Tensor2S
   
   type Tensor4
    ! tensor of rank 4
    real(kind=8), dimension(3,3,3,3) :: abcd = 0.
   end type Tensor4
   
   type Tensor4S
    ! symmetric tensor of rank 4 stored as 6x6 matrix
    real(kind=8), dimension(6, 6) :: a6b6 = 0.
   end type Tensor4S

This might be something to add, in case there is no reason against it?

I think the docs could be improved.

• Switch to the Just the Docs theme: https://just-the-docs.com
• Activate LaTeX math
• Mention the Abaqus Umat Examples now they are actually usable 🧮

Use the diataxis approach: https://diataxis.fr/: Add Tutorials, How-to (later)
Read the Docs would be great for versionated docs but I'm not sure if the complexity of Sphinx-Docs is worth here.

• ...?

Hi,
I am a new user of fortran trying to use the toolbox to help writing a umat in abaqus. I was trying to start running the example script_umat.f but it couldn't compile in gfortran 6.3.0. I got many errors saying "Derived type 'Tensor1' is being used before it is defined." I wonder if you can give some suggestions on how to fix these errors.
Thank you!

as already pointed out in #37 ,

the result of the crossed-dyadic product of two symmetric tensors $\boldsymbol{A}$ and $\boldsymbol{B}$

$$C_{ijkl} = \frac{1}{2} \left(A_{ik} \cdot B_{jl} +A_{il} \cdot B_{jk} \right) $$

is not symmetric and hence, can't be stored as a 6x6 tensor.

The implementation is wrong. The current implementation is only valid if $\boldsymbol{A} = \boldsymbol{B}$!

Hi,

I know that the library has not been tested in ubuntu but I would like to ask if anyone knows how to read .F extensions in Abaqus in Ubuntu.
The library is not working because the pre-processor is not called.

abaqus14 job=rectangle_fibers user=hgo_crimpingCopy.F
Abaqus Error: The following file(s) could not be located: hgo_crimpingCopy.F.f, hgo_crimpingCopy.F.for, hgo_crimpingCopy.F.c, hgo_crimpingCopy.F.cpp, hgo_crimpingCopy.F.C, hgo_crimpingCopy.F.c++, hgo_crimpingCopy.F.o
Abaqus/Analysis exited with error(s).

If I rename the file with a .f, there is a pre-processor violation and *Use tensor does not work.

Thanks in advance,

Iulen,

Voigt Notation is in ordering (11,22,33,12,23,31). Abaqus uses different ordering. Add option to allow both orderings in the future.

Hi @adtzlr ,
When I use the toolbox to write Abaqus subroutine, T%ab(i,j) gives me T(j,i) instead of T(i,j). This doesn't not happen when I test the subroutine outside of the Abaqus framework, (e.g. in Intel Fortran). And I couldn't debug what is the cause of this issue. Have you experienced similar issue or could you give some suggests how I could fix this?
Thank you,
Yiling

Voigt notation is possible but marked as 'experimental' because some functions could be wrong and some are still missing.

The result of the inner crossed-dyadic product of two symmetric tensors $\boldsymbol{A}$ and $\boldsymbol{B}$

$$C_{ijkl} = A_{ik} \cdot B{jl} $$

is unsymmetric and hence, can't be stored as a 6x6 tensor.

$$C_{ijkl} \ne C_{jikl} \qquad \text{and} \qquad C_{ijkl} \ne C_{ijlk} $$

The implementation is wrong.

...remove it 🗑️

Hi Andreas @adtzlr ,
I don't want to hold you up, so the following is absolutely not a priority (just nice to have). However, would you mind giving me a Go/Nogo answer soon?

For an upcoming presentation I would like to refer to the ttb (and advertise it). Just for the fun I played around with a logo some time ago, because I think we're still missing one for the toolbox (always nicer than showing a written URL).
This is what I came up with (editable .svg file appended):

I just wanted to check in with you, whether I may use this logo for the time being. Of course, it's your toolbox, you decide and we can alter the logo to your liking.

Another question that arises in this context, is how the toolbox could be cited (others for instance: https://www.tensortoolbox.org/, https://sites.google.com/site/aenader/umat-workshop). But that does not need to be elaborated in the near future.

Thanks and once again, lots of appreciation for the work you've put into the toolbox, it facilitates and improves my work majorly,
Johannes

PNG and SVG:
ttb logo.zip

Hi, My work involves using St. Venant-Kirchoff material. How can I use this subroutine in my CalculiX code?

The implementation of function piola4s is very slow (even slower than full tensor version piola4) because first the fourth order tensor is converted to tensor storage, then the piola operation is performed and lastly it is converted back to Voigt notation.

Hi, many thanks for your efforts as I was trying to heavily use your library in the development of UMAT. However, due to software limitation, only the .f file is compatible with my current settings. I am wondering if it is easy to integrate your library by calling a .f file instead of .F file?

Thanks!

Hi @adtzlr ,
I got unexpected results when I tried to use dot_11 to calculate dot product of two vectors.
I think it should be
dot_11 = dot_11 + T1%a(i)*T2%a(i)
instead of
dot_11 = dot_11 + T1%a(i)*T2%a(j)

Please correct me if I am wrong.

The output should be Tensor2 instead of Tensor2s.

ideally within a CI/CD pipeline using Github Actions.

Hi, Andreas:

I have a couple questions for this umat example using ttb.

First, have you tested running of this umat file on Abaqus? I am assuming not, because apparently the updated stress (STRESS) and constitutive tangent (DDSDDE) is not defined, despite they were assigned by s and d in your code.

Secondly, since I was running the umat developed on my own but it always have convergent issues. I noticed you have an additional step to push forward to Jaumann tangent rate of Cauchy stress for Abaqus. In the past, I always understand the DDSDDE required by Abaqus can be approximated by dSdE, where S is the 2nd PK stress. However, it seems that should be dsigma/dD, where sigma is Cauchy stress and D is the displacement, Would you mind expand a little bit about how did you get the equation on line 60 of your example code.

Best,

due to a recent switch of the Jekyll theme and a switch to GitHub Actions.