!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 ABQinterface
    
    IMPLICIT NONE
    
    INTEGER(KIND=AbqIK),INTENT(IN)::LOP,LRESTART,KSTEP,KINC
    REAL(KIND=AbqRK),DIMENSION(2),INTENT(IN)::TIME
    REAL(KIND=AbqRK),INTENT(IN)::DTIME
    
    IF (LOP==0) THEN
        CALL readdata !get data that define the micro mesh
    END IF

    IF ((LOP==1) .OR. (LOP==2) .OR. (LOP==3)) THEN
        CALL manage_STATEV_data(LOP,DTIME,TIME,KINC,LRESTART)
    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(:),ALLOCATABLE::para
    !this type bundles all convergence/analysis parameters
    TYPE (analysisparameters)::analysis
    INTEGER::rve_number
    
    CALL mkl_set_num_threads(1) !use 1 thread for each internal parallel
                                !structure in MKL functions 
    
    CALL STDB_ABQERR(1_AbqIK, 'MonolithFE2 Version 2.0 reports:',0, 0.0_AbqRK, ' ')
    
    !---------------------------read analysis parameters-----------------------
    
    CALL analysis%read_data()
    
    !-----------------------read micro mesh information------------------------
    rve_number=1
    DO WHILE (.TRUE.)
        IF (ALLOCATED(para)) DEALLOCATE(para); ALLOCATE(para(1))
        para(1)%rve_number=rve_number
        CALL para(1)%read_data(analysis) !read the information
        IF (analysis%active_RVE_definitions(1)==rve_number) THEN
            !compute permutation matrix at the beginning of the analysis
            IF (analysis%solving_process==0) THEN
                CALL get_permutation_matrix(para(1),analysis)
            END IF
            rve_number=rve_number+1
        ELSE
            EXIT
        END IF
    END DO
    
    !---------------------------output possible error--------------------------
    
    !Stop the Analysis, if not at least 1 mesh information file is supplied
    IF (analysis%active_RVE_definitions(1)==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,TIME,KINC,LRESTART)
    !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(:),ALLOCATABLE::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,INTENT(IN):: KINC
    !DTIME of current increment
    REAL(KIND=AbqRK),INTENT(IN):: DTIME
    !TIME of current increment
    REAL(KIND=AbqRK),DIMENSION(2),INTENT(IN)::TIME
    !running indices
    INTEGER::i,j,k,l,m,e,rve_number,error
    !name of the current abaqus JOB
    CHARACTER(LEN=256):: JOBNAME
    !lenght of thatname
    INTEGER:: LENJOBNAME
    INTEGER,INTENT(IN):: LOP
    !a variable that is of no interest when called from Abaqus, but when called
    !from UMAT_Driver then this variable contains the Macro GP)
    INTEGER(KIND=AbqIK),INTENT(IN):: LRESTART
    !Name of directory in which UMAT_Driver is stared
    CHARACTER(LEN=256)::  OUTDIR
    !lenght of name of this directory
    INTEGER:: LENOUTDIR
    INTEGER:: unit_number
    !get the force snapshots when being in hyper ROM training mode
    REAL(KIND=AbqRK),DIMENSION(:,:),ALLOCATABLE:: F_Snaps
    REAL(KIND=AbqRK),DIMENSION(:),ALLOCATABLE:: F_int
    REAL(KIND=AbqRK):: Volume
    
    !regularly free unused memory allocated by MKL
    CALL mkl_free_buffers
    
    !get Abaqus directory
    CALL GETOUTDIR(OUTDIR,LENOUTDIR)
    
    !-------------------get pointer to analysis parameters---------------------

    CALL analysis%get_pointer()
    
    ALLOCATE(para(analysis%active_RVE_definitions(1)))
    DO rve_number=1,analysis%active_RVE_definitions(1)
        para(rve_number)%rve_number=rve_number
        CALL para(rve_number)%get_pointer(analysis)
    END DO
    
    !---------loop over all (possibly existing) macro gausspoints--------------

    DO i=1,analysis%max_n_GP_macro_total
    
        IF (LRESTART<0) THEN
            m=-LRESTART
        ELSE
            m=i
        END IF

        IF (analysis%active_gausspoints(m)/=0) THEN!check if gausspoint is active

            rve_number=analysis%active_gausspoints(m)
            GP_DATA%GP_ID=m
                            
            !get pointer to GP data
            CALL GP_DATA%get_pointer(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) THEN
                    !initialize displacements and strain

                    GP_DATA%UGLOBAL_t_1 =GP_DATA%UGLOBAL_t
                    GP_DATA%Grad_U_k     =GP_DATA%Grad_U_t
                
                ELSE
                  !extrapolate displacements and strain:
                  
                  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)
                                        
                END IF
                
                !set value of step_indicator to 1 ->beginnig of increment
                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
                !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
                !save current time increment DTIME to DTIME_old_increment
                GP_DATA%DTIME_old_increment(1)=DTIME
                
                !----------check in training mode if data has to be outputed---------------
                
                IF (analysis%training==1) THEN
                    IF (GP_DATA%i_data_dump<analysis%n_data_dumps) THEN
                        IF (TIME(1)>=analysis%data_dump(GP_DATA%i_data_dump+1) .AND. (TIME(1)>0.0_AbqRK)) THEN
                        
                            GP_DATA%i_data_dump=GP_DATA%i_data_dump+1
                            DO WHILE (GP_DATA%i_data_dump<=analysis%n_data_dumps)
                                IF (GP_DATA%i_data_dump==analysis%n_data_dumps) EXIT
                                IF (TIME(1)>=analysis%data_dump(GP_DATA%i_data_dump+1)) THEN
                                    GP_DATA%i_data_dump=GP_DATA%i_data_dump+1
                                ELSE
                                    EXIT
                                END IF
                            END DO
                            
                            !get the Jobname
                            WRITE(JOBNAME,*) m-analysis%max_n_macro_GP
                            
                            unit_number=100+m
                            
                            IF (analysis%solving_process==0) THEN
                                OPEN(unit=unit_number,file=OUTDIR(1:LENOUTDIR)//'/training_data-u-'//TRIM(ADJUSTL(JOBNAME))//'.txt',action='WRITE',position='append',status='unknown')
                                DO k=1,para(1)%nRows_reduced
                                    WRITE (unit_number ,fmt='(E20.13)',advance='no') GP_DATA%UGLOBAL_t_1(para(1)%master_reduced_to_global(k))
                                    IF (k<para(1)%nRows_reduced) WRITE(unit_number,fmt='(A1)',advance='no') ','
                                END DO
                                CLOSE(unit=unit_number)
                            ELSE IF (analysis%solving_process==1) THEN
                                
                                !F_Snaps contains the Snapshots of the forcecontribution to all Modes (including macroscopic strain) of all IPs for one timestep
                                !F_int contains the integrated forces
                                ALLOCATE(F_Snaps(para(1)%nElem*para(1)%NGP_local,para(1)%n_ROM_modes+para(1)%ndof_macro),F_int(para(1)%n_ROM_modes+para(1)%ndof_macro))
                                
                                !compute contribution of one integration point to the force of all modes for one timestep
                                DO e=1,para(1)%nElem
                                   DO k=1,para(1)%NGP_local
                                      F_Snaps((e-1)*para(1)%NGP_local+k,:para(1)%n_ROM_modes)=MATMUL(TRANSPOSE(para(1)%ROM_modes_elem(:,:,e)),GP_DATA%F_int_ip((e-1)*para(1)%NGP_local+k,:))
                                      F_Snaps((e-1)*para(1)%NGP_local+k,para(1)%n_ROM_modes+1:)=MATMUL(TRANSPOSE(para(1)%equations_elem(:,:,e)),GP_DATA%F_int_ip((e-1)*para(1)%NGP_local+k,:))
                                   END DO
                                END DO
                                
                                !substract the volume average from the force snapshots
                                F_int=MATMUL(TRANSPOSE(F_Snaps),para(1)%IP_Weights)
                                Volume=SUM(para(1)%IP_Weights) !compute the RVE Volume (without holes)
                                FORALL(i=1:size(F_int)) F_Snaps(:,i)=F_Snaps(:,i)-(F_int(i)/Volume)
                                
                                OPEN(unit=unit_number,file=OUTDIR(1:LENOUTDIR)//'/training_data-f-'//TRIM(ADJUSTL(JOBNAME))//'.txt',action='WRITE',position='append',status='unknown')
                                DO j=1,size(F_Snaps,2)
                                    DO k=1,size(F_Snaps,1)
                                        WRITE (unit_number ,fmt='(E20.13)',advance='no') F_Snaps(k,j)
                                        IF (k<size(F_Snaps,1)) WRITE(unit_number,fmt='(A1)',advance='no') ','
                                    END DO
                                    IF (j<size(F_Snaps,2)) WRITE(unit_number,fmt='(A1)',advance='no') NEW_LINE('A')
                                END DO
                                
                                CLOSE(unit=unit_number)
                                
                            END IF
                        END IF
                    END IF
                END IF

            ELSE IF ((LOP==3) .AND. (LRESTART<0)) THEN !at the simulation end deallocate data
                
                CALL GP_DATA%deallocate_data(analysis)    
                
            END IF

        END IF
        
        IF (LRESTART<0) EXIT
    
    END DO
    
    !------------output integration point weights when being in training mode-----------
    IF ((LOP==3) .AND. (LRESTART>=0) .AND. (analysis%training==1) .AND. (analysis%solving_process==1)) THEN
        OPEN(unit=101,file=OUTDIR(1:LENOUTDIR)//'/IP_Weights.txt',action='WRITE',status='new')
        DO k=1,para(1)%NGP_local*para(1)%nElem
            WRITE (101 ,fmt='(E16.9)',advance='no') para(1)%IP_Weights(k)
            IF (k<para(1)%NGP_local*para(1)%nElem) WRITE(101,fmt='(A1)',advance='no') ','
        END DO
        CLOSE(unit=101)
    END IF
        
    END SUBROUTINE manage_STATEV_data
