!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
!
!=============================================================================

INCLUDE 'mkl_spblas.f90'           !type which bundles sparse matrix in CSR form
INCLUDE 'UEL.f'                    !includes the user element library
INCLUDE 'SMAAspUserArraysFortran.f'!interface for converting Abaqus Allocatable
                                   !Cray Pointers to Fortran Pointers
INCLUDE 'type_analysisparameters.f'!includes the type "analysisparameters"
INCLUDE 'type_meshparameters.f'    !includes the type "meshparameters"
INCLUDE 'type_macro_GP_DATA.f'     !includes the type "macro_GP_DATA"
INCLUDE 'type_systemmatrices.f'    !includes the type "systemmatrices"
INCLUDE 'Solver.f'                 !Solver Routine interface
INCLUDE 'FE2MODULE.f'              !main FE2 program
INCLUDE 'manage_data.f'            !UEXTERNERLAB, readdata, manage_STATEV_data
INCLUDE 'utility_transform_stress_stiffness.f' !delivers STRESS & DDSDDE in the
                                   !way Abaqus is expecting it


    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 FE2MODULE
    USE ABQSMA
    USE type_macro_GP_DATA
    USE type_analysisparameters
    USE utiliy_transfer_stress_stiffness
                          
    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 (UGLOBAL,STATEV,...)
    TYPE(macro_GP_DATA)::GP_DATA
    !displacement gradient of the current increment
    REAL(KIND=AbqRK),DIMENSION(:),ALLOCATABLE:: Grad_U_k_1
    !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
    !stress in thickness direction (plane strain)
    REAL(KIND=AbqRK):: STRESS33
    !determinant of deformation gradient
    REAL(KIND=AbqRK):: det_DFGRD1
    !--------------------------------------------------------------------------
    
    !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 (STATEV,displacements, etc...)
        CALL GP_DATA%get_pointer(para,analysis)

    END IF
    
    !--------------------------------------------------------------------------
    
    !allocate displacement gradient of current NR increment
    ALLOCATE(Grad_U_k_1(para%ndof_macro))
    !allocate material tangent to be inserted into the main program
    ALLOCATE(DDSDDE_main(para%ndof_macro,para%ndof_macro))
    !allocate 1. Piola Kirchhoff Stress
    ALLOCATE(PK1(para%ndof_macro))
    
    !----------get displacement gradient of current NR increment k+1-----------
    
    IF (JSTEP(3)==0) THEN !small deformation theory
    
        IF (NTENS==para%ndof_macro) THEN !3D macro - 3D micro; 2D macro - 2D micro
        
            Grad_U_k_1=STRAN+DSTRAN
            
        ELSE IF (para%dimens==2) THEN !plane strain macro - 2D micro
        
            Grad_U_k_1=[STRAN(1)+DSTRAN(1),STRAN(2)+DSTRAN(2),(STRAN(4)+DSTRAN(4))]
        
        ELSE  !plane strain/axisymmetric macro - 3D micro
        
            Grad_U_k_1(1:NTENS)=STRAN+DSTRAN
            Grad_U_k_1(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%ndof_macro
            i=voigt_notation(1,k,para%dimens)
            j=voigt_notation(2,k,para%dimens)
            IF (i==j) THEN
                Grad_U_k_1(k)=STRETCH(i,j)-1.0_AbqRK
            ELSE
                Grad_U_k_1(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%UGLOBAL_t_1=GP_DATA%UGLOBAL_t+(DTIME/GP_DATA%DTIME_old_increment(1))*&
                                (GP_DATA%UGLOBAL_t-GP_DATA%UGLOBAL_old)
            GP_DATA%Grad_U_k=GP_DATA%Grad_U_t+(DTIME/GP_DATA%DTIME_old_increment(1))*&
                            (GP_DATA%Grad_U_t-GP_DATA%Grad_U_old)
            
        ELSE
    
            !only initialize displacements in staggered scheme without extrapolation
            !when KINC<3
            
            GP_DATA%UGLOBAL_t_1=GP_DATA%UGLOBAL_t
            GP_DATA%Grad_U_k=GP_DATA%Grad_U_t
            
        END IF
    END IF
    
    !--------------------------call the main program---------------------------
    
    IF (analysis%monolithic) THEN !monolithic
    
        CALL main_program_monolithic(PK1,DDSDDE_main,Grad_U_k_1,GP_DATA,&
                                    DTIME,TIME,para,PNEWDT,JSTEP(3),JSTEP(1),&
                                    KINC,analysis,STRESS33)
    
    ELSE !staggered 
    
        CALL main_program_staggered(PK1,DDSDDE_main,Grad_U_k_1,GP_DATA,&
                                    PNEWDT,DTIME,TIME,para,JSTEP(3),JSTEP(1),&
                                    KINC,analysis,STRESS33)
    
    END IF
    
    !write displacement gradient of NR LOOP K+1 to LOOP K
    GP_DATA%Grad_U_k=Grad_U_k_1
    
    !--------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,STRESS33,ROTATION,&
                                                         DFGRD1_inv,det_DFGRD1)
    END IF
    
    IF (JSTEP(3)==1) DEALLOCATE(STRETCH,ROTATION,DFGRD1_inv) 
    
    !-----write strain resp. displ. gradient to STATEV for Postprocessing------
    
    IF (NSTATV>=para%ndof_macro) STATEV(1:para%ndof_macro)=Grad_U_k_1
    
    !-----------------------update step indicator------------------------------
    IF (GP_DATA%step_indicator(1)==-1) THEN
        GP_DATA%step_indicator(1)=2
    ELSE
        GP_DATA%step_indicator(1)=GP_DATA%step_indicator(1)+1
    END IF
    
    !----------------------deallocate allocated arrays-------------------------
        
    DEALLOCATE(Grad_U_k_1)
    DEALLOCATE(DDSDDE_main)
    DEALLOCATE(PK1)
    
    END SUBROUTINE
