!DEC$ FREEFORM
!=============================================================================
! Monolithic FE^2
! Nils Lange, Geralf Huetter, Bjoern Kiefer
!   Nils.Lange@imfd.tu-freiberg.de, Geralf.Huetter@imfd.tu-freiberg.de, 
!   Bjoern.Kiefer@imfd.tu-freiberg.de
! distributed under CC BY-NC-SA 4.0 license
! (https://creativecommons.org/licenses/by-nc-sa/4.0/)
! Reference: 
!   N. Lange, G. Huetter, B. Kiefer: "An efficient monolithic solution scheme for FE2 problems",
!   https://arxiv.org/abs/2101.01802
!
! Further information on the implementation, structure of the source code,
! examples and tutorials can be found in the file doc/documentation.pdf
!
!=============================================================================

MODULE utiliy_transfer_stress_stiffness
!module with utility routines to transfer stress and stiffness from the internal used
!conventions of MonolithFE2 into the conventions of Abaqus

USE type_meshparameters
USE ABQINTERFACE

IMPLICIT NONE

!Notation as defined in Abaqus Convention
INTEGER,DIMENSION(2,6,3),PARAMETER:: voigt_notation=reshape([1,1,0,0,0,0,0,0,0,0,0,0,&
                                                             1,1,2,2,1,2,0,0,0,0,0,0,&
                                                             1,1,2,2,3,3,1,2,1,3,2,3],[2,6,3])
                                                                
CONTAINS
    
    SUBROUTINE get_abaqus_stress_stiffness(STRESS,DDSDDE,PK1,DDSDDE_main,DFGRD1,&
                                           nl_geom,para,NTENS,STRESS33,ROTATION,&
                                           DFGRD1_inv,det_DFGRD1)
    !this routine outputs the stress and tangent stiffness as expected from Abaqus
    
    IMPLICIT NONE
    
    INTEGER(KIND=AbqIK):: NTENS,nl_geom
    TYPE(meshparameters)::para !this type bundles all mesh parameters
    REAL(KIND=AbqRK),DIMENSION(:,:),INTENT(INOUT):: DDSDDE_main
    REAL(KIND=AbqRK),DIMENSION(:),INTENT(INOUT):: PK1
    REAL(KIND=AbqRK),DIMENSION(:,:),INTENT(IN):: DFGRD1
    REAL(KIND=AbqRK),DIMENSION(:,:),INTENT(IN):: ROTATION
    REAL(KIND=AbqRK),DIMENSION(:,:),INTENT(IN)::DFGRD1_inv
    REAL(KIND=AbqRK),DIMENSION(:),INTENT(OUT):: STRESS
    REAL(KIND=AbqRK),DIMENSION(:,:),INTENT(OUT):: DDSDDE
    REAL(KIND=AbqRK),INTENT(IN):: STRESS33 !stress in thickness direction (plane strain)
    REAL(KIND=AbqRK),INTENT(IN):: det_DFGRD1 !determinant of deformation gradient
    
    !transform stress and stiffness in geometrically nonlinear case
    IF (nl_geom==1) CALL trans_stress_stiffness(para,PK1,DDSDDE_main,DFGRD1,ROTATION,DFGRD1_inv,det_DFGRD1)

    IF (NTENS==para%ndof_macro) THEN !3D macro - 3D micro; 2D macro - 2D micro
    
        DDSDDE=DDSDDE_main
        STRESS=PK1
    
    ELSE IF (para%dimens==2) THEN !plane strain macro - 2D micro
    
        STRESS(1:2)=PK1(1:2)
        STRESS(3)=STRESS33
        STRESS(4)=PK1(3)
    
        DDSDDE(1:2,1:2)=DDSDDE_main(1:2,1:2)
        DDSDDE(4,1:2)=DDSDDE_main(3,1:2)
        DDSDDE(1:2,4)=DDSDDE_main(1:2,3)
        DDSDDE(4,4)=DDSDDE_main(3,3)
    
    ELSE  !plane strain/axisymmetric macro - 3D micro
    
        STRESS=PK1(1:NTENS)
        DDSDDE=DDSDDE_main(1:NTENS,1:NTENS)
    
    END IF

    END SUBROUTINE get_abaqus_stress_stiffness
    
!==================================================================================

    SUBROUTINE trans_stress_stiffness(para,PK1,DDSDDE,DFGRD1,ROTATION,DFGRD1_inv,det_DFGRD1)
    !this routine takes the (symmetric) biot stress tensor and tangent biot stress after stretch
    !and transforms it to the Cauchy stress tensor and spatial tangent

    IMPLICIT NONE
    
    TYPE(meshparameters),INTENT(IN)::para !this type bundles all mesh parameters
    REAL(KIND=AbqRK),DIMENSION(:),INTENT(INOUT)::PK1 !biot stress on entry, cauchy on exit
    REAL(KIND=AbqRK),DIMENSION(:,:),INTENT(INOUT)::DDSDDE !nominal tangent on entry, spatial on exit
    REAL(KIND=AbqRK),DIMENSION(:,:),INTENT(IN)::DFGRD1 !deformation gradient
    REAL(KIND=AbqRK),DIMENSION(:,:),INTENT(IN)::DFGRD1_inv !inverse of DFGRD1
    REAL(KIND=AbqRK),DIMENSION(:,:),INTENT(IN)::ROTATION !(total) rotation matrix
    REAL(KIND=AbqRK),DIMENSION(para%ndof_macro,para%ndof_macro)::transformation_matrix_left,transformation_matrix_right
    INTEGER,DIMENSION(para%ndof_macro)::ipiv !pivot indices for LU factorization of transformation_matrix_left
    REAL(KIND=AbqRK),INTENT(IN):: det_DFGRD1 !determinant of deformation gradient
    INTEGER:: m,n,i,j,k,l,error
    
    !get the left and right transformation matrices
    DO n=1,para%ndof_macro
        k=voigt_notation(1,n,para%dimens)
        l=voigt_notation(2,n,para%dimens)
        IF (k==l) THEN
            DO m=1,para%ndof_macro
                i=voigt_notation(1,m,para%dimens)
                j=voigt_notation(2,m,para%dimens)
                transformation_matrix_left(m,n)=ROTATION(k,i)*DFGRD1_inv(j,l)+ROTATION(l,j)*DFGRD1_inv(i,k)
                IF (i==j) THEN
                    transformation_matrix_right(n,m)=ROTATION(i,k)*DFGRD1(j,l)
                ELSE
                    transformation_matrix_right(n,m)=0.5*(ROTATION(i,k)*DFGRD1(j,l)+ROTATION(j,k)*DFGRD1(i,l))
                END IF
            END DO
        ELSE
            DO m=1,para%ndof_macro
                i=voigt_notation(1,m,para%dimens)
                j=voigt_notation(2,m,para%dimens)
                transformation_matrix_left(m,n)=ROTATION(k,i)*DFGRD1_inv(j,l)+ROTATION(l,j)*DFGRD1_inv(i,k)&
                                                +ROTATION(l,i)*DFGRD1_inv(j,k)+ROTATION(k,j)*DFGRD1_inv(i,l)
                IF (i==j) THEN
                    transformation_matrix_right(n,m)=ROTATION(i,k)*DFGRD1(j,l)+ROTATION(i,l)*DFGRD1(j,k)
                ELSE
                    transformation_matrix_right(n,m)=0.5*(ROTATION(i,k)*DFGRD1(j,l)+ROTATION(i,l)*DFGRD1(j,k)+ROTATION(j,k)*DFGRD1(i,l)+ROTATION(j,l)*DFGRD1(i,k))
                END IF
            END DO
        END IF
    END DO
    
    transformation_matrix_left=transformation_matrix_left*det_DFGRD1/2.0_AbqRK
    
    !transform stress stress_cauchy=transformation_matrix_left^-1*biot_stress
    CALL dgesv(para%ndof_macro,1,transformation_matrix_left,para%ndof_macro,ipiv,PK1,para%ndof_macro,error)
    
    !transform stiffness DDSDDE=transformation_matrix_left^-1*DDSDDE_nominal*transformation_matrix_right; use dgetrs because dgesv already outputs the LU factorization
    CALL dgetrs('N',para%ndof_macro,para%ndof_macro,transformation_matrix_left,para%ndof_macro,ipiv,DDSDDE,para%ndof_macro,error)
    DDSDDE=MATMUL(DDSDDE,transformation_matrix_right)
    
    END SUBROUTINE trans_stress_stiffness

!==================================================================================
    
    SUBROUTINE polar_decomposition(para,DFGRD1,ROTATION,STRETCH,DFGRD1_inv,det_DFGRD1)
    !this routine performs the polar decomposition of the deformation gradient
    !into the right stretch and rotation tensor
        
    IMPLICIT NONE
    
    TYPE(meshparameters),INTENT(IN)::para !this type bundles all mesh parameters
    REAL(KIND=AbqRK),DIMENSION(:,:),INTENT(IN)::DFGRD1 !deformation gradient
    REAL(KIND=AbqRK),DIMENSION(:,:),INTENT(OUT)::ROTATION !rotation tensor
    REAL(KIND=AbqRK),DIMENSION(:,:),INTENT(OUT)::STRETCH !right stretch tensor
    REAL(KIND=AbqRK),DIMENSION(:,:),INTENT(OUT)::DFGRD1_inv !inverse of DFGRD1
    REAL(KIND=AbqRK),INTENT(OUT):: det_DFGRD1 !determinant of deformation gradient
    REAL(KIND=AbqRK),DIMENSION(para%dimens,para%dimens)::C !left Cauchy Green tensor
    REAL(KIND=AbqRK),DIMENSION(para%dimens,para%dimens)::temp
    REAL(KIND=AbqRK),DIMENSION(para%dimens)::eigen_values !eigenvalues of left Cauchy Green tensor
    INTEGER::i,error
    INTEGER,DIMENSION(para%dimens)::ipiv !permutation matrix
    INTEGER::lwork
    REAL(KIND=AbqRK),DIMENSION(48)::work
    
    !compute eigenvectors and -values of the left Cauchy Green tensor
    !the eigenvectors are stored in C; symmetric solver
    C=MATMUL(TRANSPOSE(DFGRD1(1:para%dimens,1:para%dimens)),DFGRD1(1:para%dimens,1:para%dimens))
    !call the eigenvalue/-vector solver -> C is overwritten with the eigenvectors
    lwork=SIZE(work)
    CALL dsyev('V','U',para%dimens,C,para%dimens,eigen_values,work,lwork,error)
    eigen_values=SQRT(eigen_values)
    STRETCH=0.0_AbqRK; DFGRD1_inv=0.0_AbqRK; det_DFGRD1=1.0_AbqRK
    DO i=1,para%dimens
        temp=MATMUL(RESHAPE(C(:,i),[para%dimens,1]),RESHAPE(C(:,i),[1,para%dimens]))
        STRETCH=STRETCH+temp*eigen_values(i)
        !temporary store the inverse of STRETCH in DFGRD1_inv
        DFGRD1_inv=DFGRD1_inv+temp/eigen_values(i)
        !get the determinant of the deformation gradient
        det_DFGRD1=det_DFGRD1*eigen_values(i)
    END DO
    !finally compute the rotation: F*U^-1=R
    ROTATION=MATMUL(DFGRD1(1:para%dimens,1:para%dimens),DFGRD1_inv)
    !compute the inverse of F -> F^-1=U^-1*R^T
    DFGRD1_inv=MATMUL(DFGRD1_inv,TRANSPOSE(ROTATION))
    
    END SUBROUTINE polar_decomposition

END MODULE utiliy_transfer_stress_stiffness
