!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 UEXTERNALDB(LOP,LRESTART,TIME,DTIME,KSTEP,KINC)
    !This Abaqus interface is called several times during the analysis.
    !LOP==0 indicates that the analysis is started. Then readdata is
    !called and reads the mesh/analysis paraemeter information which is stored
    !then in Abaqus allocatable arrays; LOP==2 indicates, that the analysis
    !converged in the time increment and the STATEVs of the time t+1
    !are to be written to the STATEVs of the time t; LOP==1 indicates
    !that a time increment starts, then the arrays of time t are to be
    !written to the arrays with time t+1. LOP==3 indicates that the
    !analysis is finished. Then all the arrays are to be deallocated.
    
    USE class_Solver

    INCLUDE 'ABA_PARAM.INC'

    DIMENSION TIME(2)
                    
    IF (LOP==0) THEN

        CALL readdata !get data that define the micro mesh

        CALL mkl_set_num_threads(1) !use 1 thread for each internal
                                    !parallel structure in MKL functions 
        
    END IF
    
    IF ((LOP==1) .OR. (LOP==2) .OR. (LOP==3)) THEN
        CALL manage_STATEV_data(LOP,DTIME,KINC)
    END IF
        
    RETURN
    END SUBROUTINE UEXTERNALDB

    
    SUBROUTINE readdata
    !Read the informations for the RVE meshes, coming from python script
    !FE2 Inputgenerator. This subroutine is only called once at the start
    !of the macro analysis. The arrays are saved in abaqus allocatable global
    !arrays (threadsafe; important for multithreading). Besides analysis
    !parameters are read from the file FE2_Analysisparameters.cfg if it exists.

    USE ABQINTERFACE
    USE type_meshparameters
    USE ABQSMA
    USE class_Solver
    
    IMPLICIT NONE
    
    !this type bundles all mesh parameters
    TYPE(meshparameters),DIMENSION(5)::para
    !this type bundles all convergence/analysis parameters
    TYPE (analysisparameters)::analysis
    INTEGER::rve_number,error
    
    CALL STDB_ABQERR(1_AbqIK, 'MonolithFE2 Version 1.02 reports:',0, 0.0_AbqRK, ' ')
    
    !---------------------------read analysis parameters-----------------------
    
    CALL analysis%read_data()
    
    !-----------------------read micro mesh information------------------------
    
    DO rve_number=1,5
        CALL para(rve_number)%read_data(analysis,rve_number)
        IF (analysis%active_RVE_definitions(rve_number)/=0) THEN
            !compute permutation matrix at the beginning of the analysis
            error=get_permutation_matrix(para(rve_number),analysis)
        END IF
    END DO
    
    !---------------------------output possible error--------------------------
    
    !Stop the Analysis, if not at least 1 mesh information file is supplied
    IF (NOT(ANY(analysis%active_RVE_definitions/=0))) THEN
        CALL STDB_ABQERR(-3_AbqIK, 'No Inputfile supplied!',0, 0.0_AbqRK, ' ')
    END IF
    
    END SUBROUTINE readdata


    SUBROUTINE manage_STATEV_data(LOP,DTIME,KINC)
    !This Subroutine is called if the timestep converged (LOP==2). Then the
    !STATEVs from the STATEV t+1 arrays are to be written to the t
    !arrays; If the step didnt converged, or the next step begings, the STATEVs
    !from the t arrays are to be written to the t+1 array (LOP=1); If the
    !Analysis is finished (LOP==3), all allocated arrays are to be deallocated
    
    USE ABQINTERFACE
    USE type_meshparameters
    USE type_macro_GP_DATA
    USE ABQSMA
    USE class_Solver
    USE type_analysisparameters
    
    IMPLICIT NONE

    !this type bundles all mesh parameters
    TYPE(meshparameters),DIMENSION(5)::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
    !solver type declared in module class_Solver
    TYPE(solve_soe):: solver
    !current increment number
    INTEGER:: KINC
    !DTIME of current increment
    REAL(KIND=AbqRK):: DTIME
    !running indices
    INTEGER::i,j,rve_number,error
    INTEGER,INTENT(IN):: LOP
        
    !-------------------get pointer to analysis parameters---------------------
    CALL analysis%get_pointer()
    
    DO rve_number=1,5
        IF (analysis%active_RVE_definitions(rve_number)/=0) THEN
            CALL para(rve_number)%get_pointer(rve_number)
        END IF
    END DO
    
    !---------loop over all (possibly existing) macro gausspoints--------------

    DO i=1,analysis%max_n_GP_macro_total
                                         
        IF (analysis%active_gausspoints(i)/=0) THEN!check if gausspoint is active
        
            rve_number=analysis%active_gausspoints(i)
                            
            !get pointer to GP data
            CALL GP_DATA%get_pointer(i,para(rve_number),analysis)
        
            IF (LOP==1) THEN    !start of an increment
                !write time t array to time t+1 array (necessary if convergence
                !fails and smaller increment begins in monolithic analysis)
                !==> initialization
                
                IF ((KINC<3) .AND. (analysis%monolithic)) THEN
                    !initialize displacements and strain in monolithic case:

                    GP_DATA%UGLOBAL_t_1 =GP_DATA%UGLOBAL_t
                    GP_DATA%Grad_U_k     =GP_DATA%Grad_U_t
                
                ELSE IF ((KINC>=3) .AND. (analysis%monolithic)) THEN
                  !extrapolate displacements and strain in monolithic case:
                  
                  GP_DATA%UGLOBAL_t_1=GP_DATA%UGLOBAL_t+&
                  (DTIME/analysis%DTIME_old_increment(1))*&
                  (GP_DATA%UGLOBAL_t-GP_DATA%UGLOBAL_old)
                  
                  GP_DATA%Grad_U_k=GP_DATA%Grad_U_t+(DTIME/analysis%DTIME_old_increment(1))*&
                  (GP_DATA%Grad_U_t-GP_DATA%Grad_U_old)
                                        
                END IF
                
                !set value of step_indicator to 1 ->beginnig of increment
                IF (analysis%monolithic) GP_DATA%step_indicator(1)=1
                
            ELSE IF (LOP==2) THEN   !end of an (converged!) increment
                        
                !save node displacement values of time t to time t-1
                GP_DATA%UGLOBAL_old=GP_DATA%UGLOBAL_t
                !save node displacement values of time t+1 to time t
                GP_DATA%UGLOBAL_t=GP_DATA%UGLOBAL_t_1
                !save STATEV'2 of time t+1 to time t
                GP_DATA%STATEV_t=GP_DATA%STATEV_t_1
                
                IF (analysis%monolithic) THEN

                    !save strain of time t to time t-1
                    GP_DATA%Grad_U_old=GP_DATA%Grad_U_t
                    !save strain of time t+1 to time t
                    GP_DATA%Grad_U_t=GP_DATA%Grad_U_k
                    
                END IF
                
            ELSE IF (LOP==3) THEN    !at end of anaylsis deallocate all GP_DATA
                                
                !deallocate solver data if para%save_soe==.TRUE.
                IF (analysis%monolithic .AND. analysis%save_soe) THEN
                    
                    !get pointer to solver storage
                    DO j=1,64
                        solver%handle(j)%dummy=GP_DATA%handle_nbr(j)
                    END DO
                    !deallocate data
                    error=solver%finish() 
                
                END IF
                
                CALL GP_DATA%deallocate_data(analysis,i)    
                
            END IF
        
        END IF
    
    END DO
      
    IF (LOP==2) THEN
        !save current time increment DTIME to DTIME_old_increment
        analysis%DTIME_old_increment(1)=DTIME

        !free unused memory allocated by MKL
        CALL mkl_free_buffers

    END IF
        
    IF (LOP==3) THEN   !at end of anaylsis delete arrays with micromesh data
    
        !deallocate mesh information data
        DO rve_number=1,5
            IF (analysis%active_RVE_definitions(rve_number)/=0) THEN
                CALL para(rve_number)%deallocate_data(rve_number)
            END IF
        END DO
        
        !deallocate analysis parameter data
        CALL analysis%deallocate_data()
                        
    END IF
    
    END SUBROUTINE manage_STATEV_data
