!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 '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 'Solver.f'                 !Direct Sparse Solver Routine interface
INCLUDE 'FE2MODULE.f'              !main FE2 program
INCLUDE 'manage_data.f'            !UEXTERNERLAB, readdata, manage_STATEV_data
INCLUDE 'utility_transform_stress_stiffness.f'


    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
                          
    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
    !definition number of the meshdata
    INTEGER::rve_number
    !ID of the macro Gausspoint
    INTEGER:: GP_ID
    INTEGER::n1,n2
    !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
    !running indice
    INTEGER:: i
    !stress in thickness direction (plane strain)
    REAL(KIND=AbqRK):: STRESS33
    !--------------------------------------------------------------------------
    
                   
    !get RVE label of the GP
    rve_number=int(PROPS(1))
    
    !ID of the called GP
    GP_ID=analysis%max_n_macro_GP*NOEL+NPT
    
    !get pointer to analysis parameters
    CALL analysis%get_pointer()
        
    !get pointers to the mesh information
    CALL para%get_pointer(rve_number)
    
    !------------------------------ANALYSIS BEGIN------------------------------
    IF (analysis%active_gausspoints(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
       
       !------------------------get n_STATEV_elem------------------------------
       
        !get n_STATEV_elem (number of sol. dep. state variables of one micro element)
        !n_STATEV_elem depends on the used Elementtype (JTYPE) and the PROPS
        !-> GET_n_STATEV_elem(n_STATEV_elem,JTYPE,PROPS,NPROPS)
       
        CALL GET_n_STATEV_elem(n1,para%JTYPE,PROPS(4),int(PROPS(2)))
        
        IF (int(PROPS(3))>0) THEN !2Materials are being used
            CALL GET_n_STATEV_elem(n2,para%JTYPE,PROPS(4+int(PROPS(2))),int(PROPS(3)))
        ELSE
            n2=0
        END IF
        
        IF (n1>n2) THEN !assigne the higer amount of STATEV_elem
            para%n_STATEV_elem=n1
        ELSE
            para%n_STATEV_elem=n2
        END IF
        
       !---------------------allocate macro GP data----------------------------
       
       IF (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(GP_ID,para,analysis)
        
       !mark GP as active and write the RVE label
       analysis%active_gausspoints(GP_ID)=rve_number

    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 pointer to macro GP data (STATEV,displacements, etc...)
    CALL GP_DATA%get_pointer(GP_ID,para,analysis)

    !----------get displacement gradient of current NR increment k+1-----------
    
    !ordering Grad_U -> 2D: 11,22,12,21 ;  3D: 11,22,33,12,21,23,32,13,31
    
    IF (JSTEP(3)==0) THEN !small deformation theory
    
        IF (para%dimens==3 .AND. NTENS==6) THEN
        
           Grad_U_k_1=[STRAN(1)+DSTRAN(1),STRAN(2)+DSTRAN(2),STRAN(3)+DSTRAN(3),&
                  (STRAN(4)+DSTRAN(4))*0.5_AbqRK,(STRAN(4)+DSTRAN(4))*0.5_AbqRK,&
                  (STRAN(6)+DSTRAN(6))*0.5_AbqRK,(STRAN(6)+DSTRAN(6))*0.5_AbqRK,&
                  (STRAN(5)+DSTRAN(5))*0.5_AbqRK,(STRAN(5)+DSTRAN(5))*0.5_AbqRK]
                   
        ELSE IF (para%dimens==2 .AND. NTENS==3) THEN
        
           Grad_U_k_1=[STRAN(1)+DSTRAN(1),STRAN(2)+DSTRAN(2),&
                  (STRAN(3)+DSTRAN(3))*0.5_AbqRK,(STRAN(3)+DSTRAN(3))*0.5_AbqRK]
                  
        ELSE IF (para%dimens==2 .AND. NTENS==4) THEN
        
           Grad_U_k_1=[STRAN(1)+DSTRAN(1),STRAN(2)+DSTRAN(2),&
                  (STRAN(4)+DSTRAN(4))*0.5_AbqRK,(STRAN(4)+DSTRAN(4))*0.5_AbqRK]
                   
        ELSE IF (para%dimens==3 .AND. NTENS==4) THEN
        
           Grad_U_k_1=[STRAN(1)+DSTRAN(1),STRAN(2)+DSTRAN(2),STRAN(3)+DSTRAN(3),&
                   (STRAN(4)+DSTRAN(4))*0.5_AbqRK,(STRAN(4)+DSTRAN(4))*0.5_AbqRK&
                    ,0.0_AbqRK,0.0_AbqRK,0.0_AbqRK,0.0_AbqRK]
                   
        END IF
                   
    ELSE !large deformation theory
    
        IF (para%dimens==3) THEN
        
           Grad_U_k_1=[DFGRD1(1,1)-1.0_AbqRK,DFGRD1(2,2)-1.0_AbqRK,&
                       DFGRD1(3,3)-1.0_AbqRK,DFGRD1(1,2),DFGRD1(2,1),&
                       DFGRD1(2,3),DFGRD1(3,2),DFGRD1(1,3),DFGRD1(3,1)]
                   
        ELSE IF (para%dimens==2) THEN
        
           Grad_U_k_1=[DFGRD1(1,1)-1.0_AbqRK,DFGRD1(2,2)-1.0_AbqRK,&
                       DFGRD1(1,2),DFGRD1(2,1)]
                                              
        END IF
    
    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 ((KINC>=3) .AND. (.NOT. analysis%monolithic)) THEN
    
        !u(t+1)=u(t)+(DTIME/DTIME_old_increment)*(u(t)-u(t-1))
        GP_DATA%UGLOBAL_t_1=GP_DATA%UGLOBAL_t+(DTIME/analysis%DTIME_old_increment(1))*&
                            (GP_DATA%UGLOBAL_t-GP_DATA%UGLOBAL_old)
        
    ELSE IF ((KINC<3) .AND. (.NOT. analysis%monolithic)) THEN

        !only initialize displacements in staggered scheme without extrapolation
        !when KINC<3
        
        GP_DATA%UGLOBAL_t_1=GP_DATA%UGLOBAL_t
        
    END IF
            
    !--------------------------call the main program---------------------------

    IF (analysis%monolithic) THEN !monolithic
    
        CALL main_program_monolithic(PK1,DDSDDE_main,Grad_U_k_1,GP_DATA,&
                                     PROPS,NPROPS,DTIME,TIME,para,PNEWDT,&
                                     JSTEP(3),JSTEP(1),KINC,analysis,STRESS33)
                                                                          
        !write displacement gradient of NR LOOP K+1 to LOOP K
        GP_DATA%Grad_U_k=Grad_U_k_1
    
    ELSE !staggered 
    
        CALL main_program_staggered(PK1,DDSDDE_main,Grad_U_k_1,GP_DATA,&
                                    PNEWDT,PROPS,NPROPS,DTIME,TIME,para,&
                                    JSTEP(3),JSTEP(1),KINC,analysis,STRESS33)
        
    END IF
      
    !--------return stiffness DDSDDE and cauchy stress STRESS------------------
    
    CALL get_abaqus_stress_stiffness(STRESS,DDSDDE,PK1,DDSDDE_main,DFGRD1,&
                                     JSTEP(3),para,NTENS,STRESS33)
    
    !-----write strain resp. displ. gradient to STATEV for Postprocessing------
    
    IF (NSTATV>=9) THEN
        IF (para%dimens==3) THEN
            STATEV(1:9)=Grad_U_k_1(1:9)
        ELSE
            STATEV(1:2)=Grad_U_k_1(1:2)
            STATEV(3)=0.0_AbqRK
            STATEV(4:5)=Grad_U_k_1(3:4)
            STATEV(6:9)=0.0_AbqRK
        END IF
    END IF
            
    !----------------------deallocate allocated arrays-------------------------
    
    DEALLOCATE(Grad_U_k_1)
    DEALLOCATE(DDSDDE_main)
    DEALLOCATE(PK1)
                         
    END SUBROUTINE
