!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
!
!=============================================================================
    SUBROUTINE get_abaqus_stress_stiffness(STRESS,DDSDDE,PK1,DDSDDE_main,DFGRD1,nl_geom,para,NTENS,STRESS33)
    !this routine outputs the stress and tangent stiffness as from Abaqus expected
    
    !ordering:
    !STRESS -> 11,22,33,12,13,23 (Abaqus convention)
    !PK1 -> 11,22,33,12,21,23,32,13,31
    !STRESS_cauchy -> 11,22,33,12,23,13

    
    USE type_meshparameters
    USE ABQINTERFACE
    
    IMPLICIT NONE
    
    INTEGER(KIND=AbqIK):: NTENS,nl_geom
    TYPE(meshparameters)::para !this type bundles all mesh parameters
    REAL(KIND=ABQrk),DIMENSION(9,9):: nominal_tangent
    REAL(KIND=ABQrk),DIMENSION(9):: PK1_in
    REAL(KIND=AbqRK),DIMENSION(para%ndof_macro,para%ndof_macro),INTENT(IN):: DDSDDE_main
    REAL(KIND=AbqRK),DIMENSION(para%ndof_macro),INTENT(IN):: PK1
    REAL(KIND=AbqRK),DIMENSION(3,3),INTENT(IN):: DFGRD1
    REAL(KIND=AbqRK),DIMENSION(6):: STRESS_cauchy
    REAL(KIND=AbqRK),DIMENSION(6,6):: tangent
    REAL(KIND=AbqRK),DIMENSION(NTENS),INTENT(OUT):: STRESS
    REAL(KIND=AbqRK),DIMENSION(NTENS,NTENS),INTENT(OUT):: DDSDDE
    REAL(KIND=AbqRK),INTENT(IN):: STRESS33 !stress in thickness direction (plane strain)

    
    SELECT CASE (nl_geom)
    CASE(0) !small deformation theory PK1==STRESS (cauchy)
    
        IF (para%dimens==3 .AND. NTENS==6) THEN !3D: makro/mikro
        
            STRESS=[PK1(1),PK1(2),PK1(3),PK1(4),PK1(8),PK1(6)]
        
            DDSDDE(1:4,1:4)=DDSDDE_main(1:4,1:4)
            DDSDDE(1,5)=DDSDDE_main(1,8);DDSDDE(1,6)=DDSDDE_main(1,6)
            DDSDDE(2,5)=DDSDDE_main(2,8);DDSDDE(2,6)=DDSDDE_main(2,6)
            DDSDDE(3,5)=DDSDDE_main(3,8);DDSDDE(3,6)=DDSDDE_main(3,6)
            DDSDDE(4,5)=DDSDDE_main(4,8);DDSDDE(4,6)=DDSDDE_main(4,6)
            DDSDDE(5,1:4)=DDSDDE_main(8,1:4);DDSDDE(5,5)=DDSDDE_main(8,8)
            DDSDDE(5,6)=DDSDDE_main(8,6)
            DDSDDE(6,1:4)=DDSDDE_main(6,1:4);DDSDDE(6,5)=DDSDDE_main(6,8)
            DDSDDE(6,6)=DDSDDE_main(6,6)
                   
        ELSE IF (para%dimens==2 .AND. NTENS==3) THEN !plane stress
        
            STRESS(1:3)=PK1(1:3)
            
            DDSDDE=DDSDDE_main(1:3,1:3)
            
        ELSE IF (para%dimens==2 .AND. NTENS==4) THEN !plane strain
        
            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 IF (para%dimens==3 .AND. NTENS==4) THEN !makro:plane strain/axisymmetric
                                                    !mikro:3D
           STRESS=PK1(1:4)
           
           DDSDDE=DDSDDE_main(1:4,1:4)
                   
        END IF
    
    CASE(1) !large deformation theory
    
        IF (para%dimens==3 .AND. NTENS==6) THEN !3D: makro/mikro
            
            CALL trans_stress_stiffness(STRESS_cauchy,tangent,DFGRD1,PK1,DDSDDE_main)

            STRESS(1:4)=STRESS_cauchy(1:4);STRESS(5)=STRESS_cauchy(6)
            STRESS(6)=STRESS_cauchy(5)
            
            DDSDDE(1:4,1:4)=tangent(1:4,1:4)
            DDSDDE(5,1:4)=tangent(6,1:4)
            DDSDDE(6,1:4)=tangent(5,1:4)
            DDSDDE(1:4,6)=tangent(1:4,5)
            DDSDDE(6,6)=tangent(5,5)
            DDSDDE(1:4,5)=tangent(1:4,6)
            DDSDDE(5,1:4)=tangent(6,1:4)
            DDSDDE(5,6)=tangent(6,5)
            DDSDDE(6,5)=tangent(5,6)
            DDSDDE(5,5)=tangent(6,6)
                   
        ELSE IF (para%dimens==2 .AND. NTENS==3) THEN !plane stress

            !rearrange DDSDDE_main
            nominal_tangent=0.0_AbqRK
            nominal_tangent(1:2,1:2)=DDSDDE_main(1:2,1:2)
            nominal_tangent(4:5,1:2)=DDSDDE_main(3:4,1:2)
            nominal_tangent(1:2,4:5)=DDSDDE_main(1:2,3:4)
            nominal_tangent(4:5,4:5)=DDSDDE_main(3:4,3:4)
            
            PK1_in=[PK1(1),PK1(2),0.0_AbqRK,PK1(3),PK1(4),0.0_AbqRK,0.0_AbqRK,0.0_AbqRK,0.0_AbqRK]
            
            !PK1=>Cauchy DDSDDE_main=d_PK1_d_defgrad=>DDSDDE=d_STRESS_d_STRAN
            CALL trans_stress_stiffness(STRESS_cauchy,tangent,DFGRD1,PK1_in,&
                 nominal_tangent)
            
            STRESS(1:2)=STRESS_cauchy(1:2)
            STRESS(3)=STRESS_cauchy(4)
            
            DDSDDE(1:2,1:2)=tangent(1:2,1:2)
            DDSDDE(3,1:2)=tangent(4,1:2)
            DDSDDE(1:2,3)=tangent(1:2,4)
            DDSDDE(3,3)=tangent(4,4)
            
        ELSE IF (para%dimens==2 .AND. NTENS==4) THEN !plane strain
            
            !rearrange DDSDDE_main
            nominal_tangent=0.0_AbqRK
            nominal_tangent(1:2,1:2)=DDSDDE_main(1:2,1:2)
            nominal_tangent(4:5,1:2)=DDSDDE_main(3:4,1:2)
            nominal_tangent(1:2,4:5)=DDSDDE_main(1:2,3:4)
            nominal_tangent(4:5,4:5)=DDSDDE_main(3:4,3:4)
            
            PK1_in=[PK1(1),PK1(2),0.0_AbqRK,PK1(3),PK1(4),0.0_AbqRK,0.0_AbqRK,0.0_AbqRK,0.0_AbqRK]
            
            !PK1=>Cauchy DDSDDE_main=d_PK1_d_defgrad=>DDSDDE=d_STRESS_d_STRAN
            CALL trans_stress_stiffness(STRESS_cauchy,tangent,DFGRD1,PK1_in,nominal_tangent)
                         
            STRESS(1:4)=STRESS_cauchy(1:4)
            STRESS(3)=STRESS33
            
            DDSDDE=tangent(1:4,1:4)
                                                       
        ELSE IF (para%dimens==3 .AND. NTENS==4) THEN !makro:plane strain/axisymmetric
                                                    !mikro:3D

            !PK1=>Cauchy DDSDDE_main=d_PK1_d_defgrad=>DDSDDE=d_STRESS_d_STRAN
            CALL trans_stress_stiffness(STRESS_cauchy,tangent,DFGRD1,PK1,DDSDDE_main)
                                                
            STRESS(1:4)=STRESS_cauchy(1:4)
            
            DDSDDE=tangent(1:4,1:4)
            
        END IF
    
    END SELECT
    
    END SUBROUTINE get_abaqus_stress_stiffness

    SUBROUTINE trans_stress_stiffness(STRESS,DDSDDE,DFGRD1,PK1,DDSDDE_nominal)
    !this routine takes the 1. Piola Kirchhoff stress tensor and nominal tangent
    !and transforms it to the Cauchy stress tensor and spatial tangent

    USE ABQINTERFACE
    USE MATH

    ! note: storage order: PK1(11,22,33,12,21,23,32,13,31)
    !                      STRESS(11,22,33,12,23,13)

    IMPLICIT NONE
    !
    REAL(KIND=ABQrk),DIMENSION(6),INTENT(OUT)::STRESS !cauchy stress
    REAL(KIND=ABQrk),DIMENSION(9),INTENT(IN)::PK1 !1.PK stress
    REAL(KIND=ABQrk),DIMENSION(3,3),INTENT(IN)::DFGRD1 !deformation gradient
    REAL(KIND=ABQrk),DIMENSION(9,9),INTENT(IN)::DDSDDE_nominal !nominal tangent
    REAL(KIND=ABQrk),DIMENSION(9,9)::DDSDDE_eulerian !eulerian tangent
    REAL(KIND=ABQrk),DIMENSION(6,6),INTENT(OUT)::DDSDDE ! -"- in voigt notation
    INTEGER(KIND=AbqRK),DIMENSION(9,2),PARAMETER::INDEX_NINE_FEAP=RESHAPE((/1,2,3,1,2,2,3,1,3,1,2,3,2,1,3,2,3,1/),(/9,2/))
    INTEGER(KIND=AbqRK),DIMENSION(6,2),PARAMETER::INDEX_SIX_FEAP=RESHAPE((/1,2,3,1,2,1,1,2,3,2,3,3/),(/6,2/))
    INTEGER(KIND=AbqRK),DIMENSION(3,3),PARAMETER::INDEX_TENS_VEC_FEAP=RESHAPE((/1,5,9,4,2,7,8,6,3/),(/3,3/))
    REAL(KIND=AbqRK),DIMENSION(3,3),PARAMETER::IDENT=RESHAPE((/1.0_AbqRK,0.0_AbqRK,0.0_AbqRK,0.0_AbqRK,1.0_AbqRK,0.0_AbqRK,0.0_AbqRK,0.0_AbqRK,1.0_AbqRK/),(/3,3/))
    INTEGER(KIND=ABQik)::i,j,k,l,n,m,l1,l2,i1,i2
    REAL(KIND=ABQrk):: det_DFGRD !determinant of deformation gradient

    !     transform Nominal stresses/moduli to Eulerian stresses/moduli
    !     STRESS^ab     <= PK1^aB DFGRD1^b_B
    !     DDSDDE_eulerian^abcd <= DDSDDE_nominal^aBcD DFGRD1^b_B DFGRD1^d_D
    !     DDSDDE_eulerian_voigt^abcd <= DDSDDE_eulerian^abcd - delta^ac STRESS^bd

    !initialize these arrays
    STRESS=0.0_AbqRK; DDSDDE=0.0_AbqRK; DDSDDE_eulerian=0.0_AbqRK

    !transform 1.Piola-Kirchhoff stress (PK1) to Kirchhoff stress (STRESS)
    DO n = 1,6
    i = INDEX_SIX_FEAP(n,1)
    j = INDEX_SIX_FEAP(n,2)
    DO l1 = 1,3
        i1 = INDEX_TENS_VEC_FEAP(i,l1)
        STRESS(n) = STRESS(n) + PK1(i1)*DFGRD1(j,l1)
    END DO ! l1
    END DO ! n

    !composition of 2. and 4. slot of tangent with deformation gradient
    DO n = 1,9
    i = INDEX_NINE_FEAP(n,1)
    j = INDEX_NINE_FEAP(n,2)
    DO m = 1,9
        k = INDEX_NINE_FEAP(m,1)
        l = INDEX_NINE_FEAP(m,2)
        DO l1 = 1,3
            DO l2 = 1,3
                i1 = INDEX_TENS_VEC_FEAP(i,l1)
                i2 = INDEX_TENS_VEC_FEAP(k,l2)
                DDSDDE_eulerian(n,m) = DDSDDE_eulerian(n,m) + DDSDDE_nominal(i1,i2)*DFGRD1(j,l1)*DFGRD1(l,l2)
            END DO ! l2
        END DO ! l1
    END DO ! m
    END DO ! n

    !get tangent in Voigt Notation
    DO n = 1,6
    i = INDEX_SIX_FEAP(n,1)
    j = INDEX_SIX_FEAP(n,2)
    i1 = INDEX_TENS_VEC_FEAP(i,j)
    DO m = 1,6
        k = INDEX_SIX_FEAP(m,1)
        l = INDEX_SIX_FEAP(m,2)
        i2 = INDEX_TENS_VEC_FEAP(k,l)
        DDSDDE(n,m) = DDSDDE_eulerian(i1,i2)
    END DO ! m
    END DO ! n
    
    !divide by J=determinant(deformation gradient): Kirchhoff -> Cauchy
    det_DFGRD=FindDet(DFGRD1) !get determinant
    DDSDDE=DDSDDE/det_DFGRD
    STRESS=STRESS/det_DFGRD

    END SUBROUTINE trans_stress_stiffness
