!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 FE² problems",
!   DOI: https://doi.org/10.1016/j.cma.2021.113886
!   N. Lange, G. Huetter, B. Kiefer: "A monolithic hyper ROM FE² method with
!                                     clustered training at finite deformations"
!   DOI: https://doi.org/10.1016/j.cma.2023.116522
!
! Further information on the implementation, structure of the source code,
! examples and tutorials can be found in the file doc/documentation.pdf
!
!=============================================================================

INCLUDE 'include.f' !include all the necessary files with Modules/Func./Subr.

    SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,RPL,DDSDDT,&
    DRPLDE,DRPLDT,STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,&
    CMNAME,NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,&
    CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,JSTEP,KINC)
    !the SUBROUTINE UMAT is the actual materialroutine interface which is called
    !from abaqus at each gausspoint; this subroutine itself only calls
    !the main program and acesses the meshdata, analysis and GP_DATA objects

    USE type_meshparameters
    USE MonolithFE2
    USE type_macro_GP_DATA
    USE type_analysisparameters
    
    IMPLICIT NONE
    
    !---------------------variables from Abaqus interface----------------------
    CHARACTER(LEN=80):: CMNAME
    INTEGER(KIND=AbqIK)::NOEL,NPT,LAYER,KSPT,KINC,NDI,NSHR,NTENS,NSTATV,NPROPS
    REAL(KIND=AbqRK),DIMENSION(NTENS):: STRESS,STRAN,DSTRAN
    REAL(KIND=AbqRK),DIMENSION(NSTATV):: STATEV
    REAL(KIND=AbqRK),DIMENSION(NTENS,NTENS):: DDSDDE
    REAL(KIND=AbqRK),DIMENSION(NTENS)::DDSDDT
    REAL(KIND=AbqRK),DIMENSION(NTENS):: DRPLDE
    REAL(KIND=AbqRK),DIMENSION(2):: TIME
    REAL(KIND=AbqRK),DIMENSION(1):: PREDEF,DPRED
    REAL(KIND=AbqRK),DIMENSION(NPROPS):: PROPS
    REAL(KIND=AbqRK),DIMENSION(3):: COORDS
    REAL(KIND=AbqRK),DIMENSION(3,3):: DROT,DFGRD0,DFGRD1
    INTEGER(KIND=AbqIK),DIMENSION(4):: JSTEP
    REAL(KIND=AbqRK):: DTIME,TEMP,DTEMP,SSE,SPD,SCD,CELENT,DRPLDT,PNEWDT,RPL
    
    !--------------------------------------------------------------------------
    !this type bundles all mesh parameters
    TYPE(meshparameters)::para
    !this type bundles all convergence/analysis parameters
    TYPE (analysisparameters)::analysis
    !this type bundles all macro GP DATA (nodal_sol,SVARS,...)
    TYPE(macro_GP_DATA)::GP_DATA
    !displacement gradient of the current increment
    REAL(KIND=AbqRK),DIMENSION(:),ALLOCATABLE:: Grad_U
    REAL(KIND=AbqRK),DIMENSION(:),ALLOCATABLE:: macro_measures
    !DDSDDE as computed by the main program
    REAL(KIND=AbqRK),DIMENSION(:,:),ALLOCATABLE:: DDSDDE_main
    !1. Piola Kirchhoff Stress
    REAL(KIND=AbqRK),DIMENSION(:),ALLOCATABLE:: PK1
    !rotation tensor
    REAL(KIND=AbqRK),DIMENSION(:,:),ALLOCATABLE::ROTATION !rotation of DFGRD1
    REAL(KIND=AbqRK),DIMENSION(:,:),ALLOCATABLE::STRETCH !right stretch tensor
    REAL(KIND=AbqRK),DIMENSION(:,:),ALLOCATABLE::DFGRD1_inv !inverse of DFGRD1
    !running indice
    INTEGER:: i,j,k
    !determinant of deformation gradient
    REAL(KIND=AbqRK):: det_DFGRD1
    !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])
    !--------------------------------------------------------------------------
    
    !ID of the called GP
    GP_DATA%GP_ID=analysis%max_n_macro_GP*NOEL+NPT

    !get pointer to analysis parameters
    CALL analysis%get_pointer()
    
    !get pointers to the mesh information
    para%rve_number=int(PROPS(1)) !get RVE label of the GP
    CALL para%get_pointer(analysis)
    
    !------------------------------ANALYSIS BEGIN------------------------------
    IF (analysis%active_gausspoints(GP_DATA%GP_ID)==0) THEN !beginning of analysis
    
        !-----------------check the dimension of the problem-------------------
    
       IF (NTENS==3 .AND. para%dimens==3) THEN
       
         CALL STDB_ABQERR(-3_AbqIK, 'Call 3D microproblems only from 3D, plane'//&
                          'strain or axisymmetric macro simulations!',0, 0.0_AbqRK, ' ')
           
       ELSE IF (NTENS==6 .AND. para%dimens==2) THEN
       
          CALL STDB_ABQERR(-3_AbqIK, 'A 3D macro GP called a 2D microproblem!',&
                           0, 0.0_AbqRK, ' ')

       END IF
        
       !---------------------allocate macro GP data----------------------------
       
       IF (GP_DATA%GP_ID>analysis%max_n_GP_macro_total) THEN
        CALL STDB_ABQERR(-3_AbqIK, 'Number of used macro elements higher then'//&
        'currently allowed. Solution: set higher max_n_GP_macro_total!',0, 0.0_AbqRK, ' ')
       END IF
       
       !allocate memory for the GP DATA (UGLOBAL,STATEV,etc.)
        CALL GP_DATA%allocate_data(para,analysis)
        
       !mark GP as active and write the RVE label
       analysis%active_gausspoints(GP_DATA%GP_ID)=para%rve_number
    
    ELSE
    
        !get pointer to macro GP data (nodal_sol,SVARS, etc...)
        CALL GP_DATA%get_pointer(para,analysis)

    END IF
    
    !--------------------------------------------------------------------------
    
    !allocate displacement gradient of current NR increment
    ALLOCATE(Grad_U(para%n_Reaction_force_dof))
    !allocate material tangent to be inserted into the main program
    ALLOCATE(DDSDDE_main(para%n_Reaction_force_dof,para%n_Reaction_force_dof))
    !allocate PK1 (==Cauchy stress in small deformations, symmetric Biot stress in finite deformations)
    ALLOCATE(PK1(para%n_Reaction_force_dof))
    
    !----------get displacement gradient of current NR increment k+1-----------
    
    IF (JSTEP(3)==0) THEN !small deformation theory
    
        IF (NTENS==para%n_Reaction_force_dof) THEN !3D macro - 3D micro; 2D macro - 2D micro
        
            Grad_U=STRAN+DSTRAN
            
        ELSE IF (para%dimens==2) THEN !plane strain macro - 2D micro
        
            Grad_U=[STRAN(1)+DSTRAN(1),STRAN(2)+DSTRAN(2),(STRAN(4)+DSTRAN(4))]
        
        ELSE  !plane strain/axisymmetric macro - 3D micro
        
            Grad_U(1:NTENS)=STRAN+DSTRAN
            Grad_U(NTENS+1:NTENS+2)=0.0_AbqRK
        
        END IF
        
    ELSE !large deformation theory
    
        ALLOCATE(STRETCH(para%dimens,para%dimens),&
                ROTATION(para%dimens,para%dimens),&
                DFGRD1_inv(para%dimens,para%dimens))
        !get polar decomposition of the deformation gradient
        CALL polar_decomposition(para,DFGRD1,ROTATION,STRETCH,DFGRD1_inv,det_DFGRD1)
        
        DO k=1,para%n_Reaction_force_dof
            i=voigt_notation(1,k,para%dimens)
            j=voigt_notation(2,k,para%dimens)
            IF (i==j) THEN
                Grad_U(k)=STRETCH(i,j)-1.0_AbqRK
            ELSE
                Grad_U(k)=STRETCH(i,j)*2.0_AbqRK
            END IF
        END DO
        
    END IF
    
    !-------initialize and extrapolate displacements (staggered scheme)--------
    
    !Only extrapolate displacements if staggered scheme is chosen and increment
    !number is higher or equal 3
    
    IF (.NOT. analysis%monolithic) THEN !staggered scheme
        IF (KINC>2) THEN !extrapolation data available
        
            !u(t+1)=u(t)+(DTIME/DTIME_old_increment)*(u(t)-u(t-1))
            GP_DATA%nodal_sol_t_1=GP_DATA%nodal_sol_t+(DTIME/GP_DATA%DTIME_old_increment)*&
                                (GP_DATA%nodal_sol_t-GP_DATA%nodal_sol_old)
            
        ELSE
    
            !only initialize displacements in staggered scheme without extrapolation
            !when KINC<3
            
            GP_DATA%nodal_sol_t_1=GP_DATA%nodal_sol_t
            
        END IF
    END IF
    
    !--------------------------call the main program---------------------------
    
    ALLOCATE(macro_measures(para%n_Reaction_force_dof+para%n_additional_dof))
    macro_measures=0.0_AbqRK
    macro_measures(1:para%n_Reaction_force_dof)=Grad_U
    
    IF (analysis%monolithic) THEN !monolithic
    
        CALL main_program_monolithic(PK1,DDSDDE_main,macro_measures,GP_DATA,&
                                    DTIME,TIME,para,PNEWDT,JSTEP(3),JSTEP(1),&
                                    KINC,analysis)
    
    ELSE !staggered 
    
        CALL main_program_staggered(PK1,DDSDDE_main,macro_measures,GP_DATA,&
                                    PNEWDT,DTIME,TIME,para,JSTEP(3),JSTEP(1),&
                                    KINC,analysis)
    
    END IF
    
    !--------return stiffness DDSDDE and cauchy stress STRESS------------------
    
    IF (PNEWDT>=1.0_AbqRK) THEN
        CALL get_abaqus_stress_stiffness(STRESS,DDSDDE,PK1,DDSDDE_main,DFGRD1,&
                                         JSTEP(3),para,NTENS,ROTATION,DFGRD1_inv,det_DFGRD1)
    END IF
    
    IF (JSTEP(3)==1) DEALLOCATE(STRETCH,ROTATION,DFGRD1_inv) 
    
    !-----write strain resp. stretch to STATEV for Postprocessing------
    
    IF (NSTATV>=SIZE(Grad_U)) STATEV(1:SIZE(Grad_U))=Grad_U
    
    !-----------------------update step indicator------------------------------
    IF (GP_DATA%step_indicator==-1) THEN
        GP_DATA%step_indicator=2
    ELSE
        GP_DATA%step_indicator=GP_DATA%step_indicator+1
    END IF
    
    CONTAINS
    
!==================================================================================
    
    SUBROUTINE get_abaqus_stress_stiffness(STRESS,DDSDDE,PK1,DDSDDE_main,DFGRD1,&
                                           nl_geom,para,NTENS,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):: det_DFGRD1 !determinant of deformation gradient
    
    !transform stress and stiffness in geometrically nonlinear case
    IF (nl_geom==1) THEN
        CALL trans_stress_stiffness(para,PK1,DDSDDE_main,DFGRD1,ROTATION,DFGRD1_inv,det_DFGRD1)
    END IF
    
    IF (NTENS==para%n_Reaction_force_dof) 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(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
    
    !this type bundles all mesh parameters
    TYPE(meshparameters),INTENT(IN)::para
    !symmetric biot stress on entry, cauchy on exit
    REAL(KIND=AbqRK),DIMENSION(:),INTENT(INOUT)::PK1
    !nominal tangent on entry, spatial on exit
    REAL(KIND=AbqRK),DIMENSION(:,:),INTENT(INOUT)::DDSDDE
    !deformation gradient
    REAL(KIND=AbqRK),DIMENSION(:,:),INTENT(IN)::DFGRD1
    !inverse of DFGRD1
    REAL(KIND=AbqRK),DIMENSION(:,:),INTENT(IN)::DFGRD1_inv
    !(total) rotation matrix
    REAL(KIND=AbqRK),DIMENSION(:,:),INTENT(IN)::ROTATION
    REAL(KIND=AbqRK),DIMENSION(para%n_Reaction_force_dof,para%n_Reaction_force_dof)::&
                                 transformation_matrix_left,transformation_matrix_right
                                 !pivot indices for LU factorization of transformation_matrix_left
    INTEGER,DIMENSION(para%n_Reaction_force_dof)::ipiv
    !determinant of deformation gradient
    REAL(KIND=AbqRK),INTENT(IN):: det_DFGRD1 
    INTEGER:: m,n,i,j,k,l,error
    
    !get the left and right transformation matrices
    DO n=1,para%n_Reaction_force_dof
        k=voigt_notation(1,n,para%dimens)
        l=voigt_notation(2,n,para%dimens)
        IF (k==l) THEN
            DO m=1,para%n_Reaction_force_dof
                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%n_Reaction_force_dof
                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%n_Reaction_force_dof,1,transformation_matrix_left,&
               para%n_Reaction_force_dof,ipiv,PK1,para%n_Reaction_force_dof,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%n_Reaction_force_dof,para%n_Reaction_force_dof,transformation_matrix_left,&
                para%n_Reaction_force_dof,ipiv,DDSDDE,para%n_Reaction_force_dof,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
    
    !this type bundles all mesh parameters
    TYPE(meshparameters),INTENT(IN)::para
    !deformation gradient
    REAL(KIND=AbqRK),DIMENSION(:,:),INTENT(IN)::DFGRD1
    !rotation tensor
    REAL(KIND=AbqRK),DIMENSION(:,:),INTENT(OUT)::ROTATION
    !right stretch tensor
    REAL(KIND=AbqRK),DIMENSION(:,:),INTENT(OUT)::STRETCH
    !inverse of DFGRD1
    REAL(KIND=AbqRK),DIMENSION(:,:),INTENT(OUT)::DFGRD1_inv
    !determinant of deformation gradient
    REAL(KIND=AbqRK),INTENT(OUT):: det_DFGRD1
    !left Cauchy Green tensor
    REAL(KIND=AbqRK),DIMENSION(para%dimens,para%dimens)::C
    REAL(KIND=AbqRK),DIMENSION(para%dimens,para%dimens)::temp
    !eigenvalues of left Cauchy Green tensor
    REAL(KIND=AbqRK),DIMENSION(para%dimens)::eigen_values
    INTEGER::i,error
    !permutation matrix
    INTEGER,DIMENSION(para%dimens)::ipiv
    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 SUBROUTINE
