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

MODULE type_macro_GP_DATA
!define a type named macro_GP_DATA which bundels all the STATEV data of
!a macro GP; the data is stored in abaqus allocatable arrays, which are
!created by the routine "allocate_macro_GP_DATA"; the data is accessed by the
!routine "get_pointer_macro_GP_DATA" and deleted by "deallocate_GP_DATA"

    USE ABQINTERFACE
    USE ABQSMA
    USE type_meshparameters
    
    IMPLICIT NONE
    
    PUBLIC
    
    TYPE macro_GP_DATA
        
        INTEGER:: GP_ID !identification number of the macro GP
        
        !--------------------saved in all cases--------------------------------
        !internal state variables for all elements at time t
        REAL(KIND=AbqRK),DIMENSION(:,:),POINTER:: STATEV_t
        INTEGER(KIND=4):: ID_STATEV_t=0
        !internal state variables for all elements at time t+1
        REAL(KIND=AbqRK),DIMENSION(:,:),POINTER:: STATEV_t_1
        INTEGER(KIND=4):: ID_STATEV_t_1=1
        !micro displacements at time t
        REAL(KIND=AbqRK),DIMENSION(:),POINTER:: UGLOBAL_t
        INTEGER(KIND=4):: ID_UGLOBAL_t=2
        !micro displacements at time t+1
        REAL(KIND=AbqRK),DIMENSION(:),POINTER:: UGLOBAL_t_1
        INTEGER(KIND=4):: ID_UGLOBAL_t_1=3
        !micro displacements at time t-1
        REAL(KIND=AbqRK),DIMENSION(:),POINTER:: UGLOBAL_old
        INTEGER(KIND=4):: ID_UGLOBAL_old=4
        !macro displacement gradient of the last iteration step
        REAL(KIND=AbqRK),DIMENSION(:),POINTER::Grad_U_k
        INTEGER(KIND=4):: ID_Grad_U_k=5
        !macro displacement gradient of the last (converged) increment at time t
        REAL(KIND=AbqRK),DIMENSION(:),POINTER::Grad_U_t
        INTEGER(KIND=4):: ID_Grad_U_t=6
        !macro displacement gradient at time t-1
        REAL(KIND=AbqRK),DIMENSION(:),POINTER::Grad_U_old
        INTEGER(KIND=4):: ID_Grad_U_old=7
        
        !----------only to be saved in monolithic case-------------------------
        !right hand side of last increment
        REAL(KIND=AbqRK),DIMENSION(:),POINTER::rhs_total_k
        INTEGER(KIND=4):: ID_rhs_total_k=8
        !partial derivative of r w.r.t macro deformation gradient of last increment
        REAL(KIND=AbqRK),DIMENSION(:,:),POINTER::d_r_d_F_k
        INTEGER(KIND=4):: ID_d_r_d_F_k=9
        !values in CSR Format of Stiffnessmatrix of last increment
        REAL(KIND=AbqRK),DIMENSION(:),POINTER::k_matrix_values_k
        INTEGER(KIND=4):: ID_k_matrix_values_k=10
        !-1 signals that the analysis began;1 signals the beginning of a timestep;
        !2-i signals the middle of a timestept
        INTEGER,DIMENSION(:),POINTER::step_indicator
        INTEGER(KIND=4):: ID_step_indicator=11
        !integer and real arrays that can be used by a solver
        INTEGER(KIND=8),DIMENSION(:),POINTER::solver_int
        REAL(KIND=AbqRK),DIMENSION(:),POINTER::solver_real
        INTEGER(KIND=4):: ID_solver_int=12,ID_solver_real=13
        !internal force vectors of all integration points
        REAL(KIND=AbqRK),DIMENSION(:,:),POINTER::F_int_ip
        INTEGER(KIND=4):: ID_F_int_ip=14
        !---------------------------------------------------------------------
        
        !integer with DTIME of old increment (needed for displ. extrapolation)
        REAL(KIND=AbqRK), DIMENSION(:),POINTER:: DTIME_old_increment
        INTEGER(KIND=4):: ID_DTIME_old_increment=15
        
        !------------------------------------------------------------------------
        !number of already outputed snapshots in training
        INTEGER,POINTER::i_data_dump
        INTEGER(KIND=4):: ID_data_dump=16

        
        INTEGER(KIND=4)::ID_total=17
        
        CONTAINS
        
            PROCEDURE:: allocate_data => allocate_data_procedure
            PROCEDURE:: allocate_solver_data => allocate_solver_data_procedure
            PROCEDURE:: get_pointer => get_pointer_procedure
            PROCEDURE:: deallocate_data => deallocate_data_procedure
        
    END TYPE macro_GP_DATA
    
    PRIVATE:: allocate_data_procedure,allocate_solver_data_procedure,&
              get_pointer_procedure,deallocate_data_procedure
    
CONTAINS
    
    SUBROUTINE allocate_data_procedure(GP_DATA,para,analysis)
    
    !bundeld macro GP DATA (UGLOBAL, STATEV...)
    CLASS(macro_GP_DATA):: GP_DATA
    !this type bundles all mesh,convergence & monolithic/staggered parameters
    TYPE(meshparameters)::para
    !this type bundles all convergence/analysis parameters
    TYPE (analysisparameters)::analysis
    INTEGER::nElem,nDof,nDof_red
    INTEGER,DIMENSION(:),POINTER::int_pointer
    
    IF (analysis%hyperintegration) THEN
        nElem=para%n_elem_to_hyper
    ELSE
        nElem=para%nElem
    END IF
    IF (analysis%solving_process==0) THEN
        nDof=para%ndof_n*para%nNodes
        nDof_red=para%nRows_reduced
    ELSE IF (analysis%solving_process==1) THEN
        nDof=para%n_ROM_modes
        nDof_red=para%n_ROM_modes
    END IF
    
    CALL AbaqusArrayCreate(GP_DATA%STATEV_t,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_STATEV_t,[para%n_STATEV_elem,nElem],0.0_AbqRK)
                                    
    CALL AbaqusArrayCreate(GP_DATA%STATEV_t_1,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_STATEV_t_1,[para%n_STATEV_elem,nElem],0.0_AbqRK)

    CALL AbaqusArrayCreate(GP_DATA%UGLOBAL_t,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_UGLOBAL_t,[nDof],0.0_AbqRK)
                                          
    CALL AbaqusArrayCreate(GP_DATA%UGLOBAL_t_1,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_UGLOBAL_t_1,[nDof],0.0_AbqRK)
   
    CALL AbaqusArrayCreate(GP_DATA%UGLOBAL_old,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_UGLOBAL_old,[nDof],0.0_AbqRK)
    
    CALL AbaqusArrayCreate(GP_DATA%step_indicator,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_step_indicator,[1],-1)
    
    CALL AbaqusArrayCreate(GP_DATA%Grad_U_k,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_Grad_U_k,[para%ndof_macro],0.0_AbqRK)
                                                       
    CALL AbaqusArrayCreate(GP_DATA%Grad_U_t,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_Grad_U_t,[para%ndof_macro],0.0_AbqRK)

    CALL AbaqusArrayCreate(GP_DATA%Grad_U_old,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_Grad_U_old,[para%ndof_macro],0.0_AbqRK)
                                          
    IF (analysis%monolithic) THEN !in monolithic case
                                                               
        CALL AbaqusArrayCreate(GP_DATA%rhs_total_k,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_rhs_total_k,[nDof_red],0.0_AbqRK)
        
        CALL AbaqusArrayCreate(GP_DATA%d_r_d_F_k,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_d_r_d_F_k,[nDof_red,para%ndof_macro],0.0_AbqRK)
        
        IF (analysis%solving_process==0) THEN !only in a full simulation save the stiffnessmatrix values
            CALL AbaqusArrayCreate(GP_DATA%k_matrix_values_k,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_k_matrix_values_k,[para%nNonZero_reduced],0.0_AbqRK)
        END IF
    
    END IF
    
    IF (analysis%training==1) THEN
        CALL AbaqusArrayCreate(GP_DATA%F_int_ip,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_F_int_ip,[para%nElem*para%NGP_local,para%ndof_e],0.0_AbqRK)
        CALL AbaqusArrayCreate(int_pointer,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_data_dump,[1],0)
        GP_DATA%i_data_dump=>int_pointer(1)
    END IF
    
    !Create Array which saves DTIME of the old increment for extrapolation
    CALL AbaqusArrayCreate(GP_DATA%DTIME_old_increment,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_DTIME_old_increment,[1],0.0_AbqRK)
        
    END SUBROUTINE allocate_data_procedure
    
    
    SUBROUTINE allocate_solver_data_procedure(GP_DATA,n_reals,n_ints)
    
    !bundeld macro GP DATA (UGLOBAL, STATEV...)
    CLASS(macro_GP_DATA):: GP_DATA
    !number of reals and ints to be allocated
    INTEGER(KIND=4),INTENT(IN):: n_reals,n_ints
    INTEGER(KIND=4):: n_r,n_i
    INTEGER(KIND=8)::INITVAL
    
    INITVAL=0
    
    IF (n_reals<1) THEN
        n_r=1
    ELSE
        n_r=n_reals
    END IF
    IF (n_ints<1) THEN
        n_i=1
    ELSE
        n_i=n_ints
    END IF

    CALL AbaqusArrayCreate(GP_DATA%solver_int,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_solver_int,[n_i],INITVAL)
    CALL AbaqusArrayCreate(GP_DATA%solver_real,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_solver_real,[n_r],0.0_AbqRK)

    END SUBROUTINE allocate_solver_data_procedure


    SUBROUTINE get_pointer_procedure(GP_DATA,para,analysis)
    
    !bundeld macro GP DATA (UGLOBAL, STATEV...)
    CLASS(macro_GP_DATA):: GP_DATA
    !this type bundles all mesh,convergence & monolithic/staggered parameters
    TYPE(meshparameters)::para
    !this type bundles all convergence/analysis parameters
    TYPE (analysisparameters)::analysis
    INTEGER::nElem,nDof,nDof_red
    INTEGER,DIMENSION(:),POINTER::int_pointer
    
    IF (analysis%hyperintegration) THEN
        nElem=para%n_elem_to_hyper
    ELSE
        nElem=para%nElem
    END IF
    IF (analysis%solving_process==0) THEN
        nDof=para%ndof_n*para%nNodes
        nDof_red=para%nRows_reduced
    ELSE IF (analysis%solving_process==1) THEN
        nDof=para%n_ROM_modes
        nDof_red=para%n_ROM_modes
    END IF
    
    CALL AbaqusArrayAccess(GP_DATA%STATEV_t,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_STATEV_t,[para%n_STATEV_elem,nElem])
                                    
    CALL AbaqusArrayAccess(GP_DATA%STATEV_t_1,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_STATEV_t_1,[para%n_STATEV_elem,nElem])

    CALL AbaqusArrayAccess(GP_DATA%UGLOBAL_t,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_UGLOBAL_t,[nDof])
                                          
    CALL AbaqusArrayAccess(GP_DATA%UGLOBAL_t_1,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_UGLOBAL_t_1,[nDof])
   
    CALL AbaqusArrayAccess(GP_DATA%UGLOBAL_old,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_UGLOBAL_old,[nDof])
    
    CALL AbaqusArrayAccess(GP_DATA%step_indicator,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_step_indicator,[1])
    
    CALL AbaqusArrayAccess(GP_DATA%Grad_U_k,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_Grad_U_k,[para%ndof_macro])
                                                       
    CALL AbaqusArrayAccess(GP_DATA%Grad_U_t,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_Grad_U_t,[para%ndof_macro])

    CALL AbaqusArrayAccess(GP_DATA%Grad_U_old,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_Grad_U_old,[para%ndof_macro])
                                          
    IF (analysis%monolithic) THEN !in monolithic case
                                                               
        CALL AbaqusArrayAccess(GP_DATA%rhs_total_k,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_rhs_total_k,[nDof_red])
        
        CALL AbaqusArrayAccess(GP_DATA%d_r_d_F_k,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_d_r_d_F_k,[nDof_red,para%ndof_macro])
        
        IF (analysis%solving_process==0) THEN !only in a full simulation save the stiffnessmatrix values
            CALL AbaqusArrayAccess(GP_DATA%k_matrix_values_k,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_k_matrix_values_k,[para%nNonZero_reduced])
        END IF
        
        IF (analysis%save_soe .OR. (analysis%solving_process==1)) THEN !if the factorization is to be saved
            CALL AbaqusArrayAccess(GP_DATA%solver_int,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_solver_int,[1])
            CALL AbaqusArrayAccess(GP_DATA%solver_real,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_solver_real,[1])
        END IF        

    END IF
    
    IF (analysis%training==1) THEN
        CALL AbaqusArrayAccess(GP_DATA%F_int_ip,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_F_int_ip,[para%nElem*para%NGP_local,para%ndof_e])
        CALL AbaqusArrayAccess(int_pointer,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_data_dump,[1])
        GP_DATA%i_data_dump=>int_pointer(1)
    END IF
    
    CALL AbaqusArrayAccess(GP_DATA%DTIME_old_increment,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_DTIME_old_increment,[1])
        
    END SUBROUTINE get_pointer_procedure
    
    
    SUBROUTINE deallocate_data_procedure(GP_DATA,analysis)
    
    !bundeld macro GP DATA (UGLOBAL, STATEV...)
    CLASS(macro_GP_DATA):: GP_DATA
    !this type bundles all convergence/analysis parameters
    TYPE (analysisparameters)::analysis
    INTEGER,DIMENSION(:),POINTER::int_pointer
    
    CALL AbaqusArrayDelete(GP_DATA%STATEV_t,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_STATEV_t)
    CALL AbaqusArrayDelete(GP_DATA%STATEV_t,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_STATEV_t_1)
    CALL AbaqusArrayDelete(GP_DATA%UGLOBAL_t,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_UGLOBAL_t)
    CALL AbaqusArrayDelete(GP_DATA%UGLOBAL_t_1,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_UGLOBAL_t_1)
    CALL AbaqusArrayDelete(GP_DATA%UGLOBAL_old,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_UGLOBAL_old)
    CALL AbaqusArrayDelete(GP_DATA%step_indicator,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_step_indicator)
    CALL AbaqusArrayDelete(GP_DATA%Grad_U_k,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_Grad_U_k)   
    CALL AbaqusArrayDelete(GP_DATA%Grad_U_t,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_Grad_U_t)
    CALL AbaqusArrayDelete(GP_DATA%Grad_U_old,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_Grad_U_old)
    
    IF (analysis%monolithic) THEN !in monolithic case
        CALL AbaqusArrayDelete(GP_DATA%rhs_total_k,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_rhs_total_k)
        CALL AbaqusArrayDelete(GP_DATA%d_r_d_F_k,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_d_r_d_F_k)
        
        IF (analysis%solving_process==0) THEN
            CALL AbaqusArrayDelete(GP_DATA%k_matrix_values_k,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_k_matrix_values_k)
        END IF
        
        IF (analysis%save_soe .OR. (analysis%solving_process==1)) THEN
            CALL AbaqusArrayDelete(GP_DATA%solver_int,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_solver_int)
            CALL AbaqusArrayDelete(GP_DATA%solver_real,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_solver_real)
        END IF

    END IF
    
    CALL AbaqusArrayDelete(GP_DATA%DTIME_old_increment,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_DTIME_old_increment)
    
    END SUBROUTINE deallocate_data_procedure

END MODULE type_macro_GP_DATA
