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

    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_SVARS_data(LOP,DTIME,TIME,KINC,LRESTART)
    END IF

    RETURN
    END SUBROUTINE UEXTERNALDB

    
    SUBROUTINE readdata
    !Read the RVE informations from the respective 1-n inputfiles. This
    !subroutine is only called once at the start of the macro analysis to get the
    !structure and information for type_meshparameters and type_analysisparameters.
    
    USE ABQinterface
    USE type_meshparameters
    USE type_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,element_ID,i,j
    
    !use 1 thread for each internal parallel structure in MKL functions
    CALL mkl_set_num_threads(1)
    
    CALL STDB_ABQERR(1_AbqIK, 'MonolithFE2 Version 3.0 reports:',0, 0.0_AbqRK, ' ')
    
    !---------------------------read analysis parameters-----------------------
    
    CALL analysis%read_data()
    
    !-----------------------read micro mesh information------------------------
    rve_number=1
    element_ID=-1000 !first element ID
    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==rve_number) THEN
            !compute permutation matrix at the beginning of the analysis
            IF (.NOT. analysis%ROM_projection) THEN
                CALL get_permutation_matrix(para(1),analysis)
            END IF
            !give every element a unique ID, so that the UEL can autonomously
            !allocate memory
            DO i=1,para(1)%nElem_types
                DO j=1,para(1)%Elements(i)%NELEM
                    para(1)%Elements(i)%JELEM(j)=element_ID
                    element_ID=element_ID-1
                END DO
            END DO
            IF (analysis%training) EXIT !only one RVE definition in training mode
            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==0) THEN
        CALL STDB_ABQERR(-3_AbqIK, 'No Inputfile supplied!',0, 0.0_AbqRK, ' ')
    END IF
    
    END SUBROUTINE readdata


    SUBROUTINE manage_SVARS_data(LOP,DTIME,TIME,KINC,LRESTART)
    !This Subroutine is called if the timestep converged (LOP==2). Then the
    !SVARS from the SVARS t+1 arrays are to be written to the t
    !arrays; If the step didnt converged, or the next step begings, the SVARs
    !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 type_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 (nodal_sol,SVARS,...)
    TYPE(macro_GP_DATA)::GP_DATA
    !solver type declared in module type_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,elem,rve_number,error
    !name of the current abaqus JOB
    CHARACTER(LEN=256):: JOBNAME
    !lenght of that name
    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
    REAL(KIND=AbqRK),DIMENSION(:),ALLOCATABLE::training_data_integrated,Norm
    
    !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))
    DO rve_number=1,analysis%active_RVE_definitions
        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 nodal solution values

                    GP_DATA%nodal_sol_t_1 = GP_DATA%nodal_sol_t
                
                ELSE
                  !extrapolate nodal solution values
                  
                  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)
                    
                END IF
                
                !set value of step_indicator to 1 ->beginnig of increment
                GP_DATA%step_indicator=1

            ELSE IF (LOP==2) THEN   !end of an (converged!) increment
                !save nodal solution values of time t to time t-1
                GP_DATA%nodal_sol_old=GP_DATA%nodal_sol_t
                !save nodal solution values of time t+1 to time t
                GP_DATA%nodal_sol_t=GP_DATA%nodal_sol_t_1
                !save SVARs of time t+1 to time t
                GP_DATA%SVARS_t=GP_DATA%SVARS_t_1
                !save current time increment DTIME to DTIME_old_increment
                GP_DATA%DTIME_old_increment=DTIME
                
                !----------check in training mode if data has to be outputed---------------

                IF (analysis%training) 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 (.NOT. analysis%ROM_projection) 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%nodal_sol_t_1(para(1)%reduced_dof_to_global_dof(1,k))
                                    IF (k<para(1)%nRows_reduced) WRITE(unit_number,fmt='(A1)',advance='no') ','
                                END DO
                                CLOSE(unit=unit_number)
                            ELSE IF (analysis%ROM_projection) THEN
                                
                                !"integrate" by summing over all element contributions
                                ALLOCATE(training_data_integrated(SIZE(GP_DATA%training_matrix,2)))
                                FORALL(k=1:SIZE(training_data_integrated)) training_data_integrated(k)=SUM(GP_DATA%training_matrix(:,k))
                                
                                DO elem=1,para(1)%n_active_elements
                                    GP_DATA%training_matrix(elem,:)=GP_DATA%training_matrix(elem,:)-training_data_integrated/para(1)%n_active_elements
                                END DO

                                ALLOCATE(Norm(SIZE(GP_DATA%training_matrix,2)))
                                DO k=1,SIZE(GP_DATA%training_matrix,2)
                                    Norm(k)=NORM2(GP_DATA%training_matrix(:,k))
                                    IF (Norm(k)>0.000001) THEN
                                        GP_DATA%training_matrix(:,k)=GP_DATA%training_matrix(:,k)/Norm(k)
                                    END IF
                                END DO

                                !output the training data
                                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(GP_DATA%training_matrix,2)
                                    IF (Norm(j)>0.000001) THEN
                                        IF (j>1 .AND. ANY(Norm(1:j-1)>0.000001)) THEN
                                            WRITE(unit_number,fmt='(A1)',advance='no') NEW_LINE('A')
                                        END IF
                                        DO k=1,size(GP_DATA%training_matrix,1)
                                            WRITE (unit_number ,fmt='(E20.13)',advance='no') GP_DATA%training_matrix(k,j)
                                            IF (k<size(GP_DATA%training_matrix,1)) WRITE(unit_number,fmt='(A1)',advance='no') ','
                                        END DO
                                    END IF
                                END DO
                                DEALLOCATE(Norm)
                                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) .AND. (analysis%ROM_projection)) THEN
        OPEN(unit=101,file=OUTDIR(1:LENOUTDIR)//'/Element_Volume.txt',action='WRITE',status='new')
        DO k=1,para(1)%n_active_elements
            WRITE (101 ,fmt='(E16.9)',advance='no') para(1)%ref_Volume_elements(k)
            IF (k<para(1)%n_active_elements) WRITE(101,fmt='(A1)',advance='no') ','
        END DO
        CLOSE(unit=101)
    END IF

    END SUBROUTINE manage_SVARS_data
