!=============================================================================
! 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
    
    PUBLIC
    
    TYPE macro_GP_DATA
    
        !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
        !pointer to handle (where factorized stiffnessmatrix is being stored)
        INTEGER(KIND=8),DIMENSION(:),POINTER::handle_nbr
        INTEGER(KIND=4):: ID_handle_nbr=11
        !-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=12
        !----------------------------------------------------------------------
        INTEGER(KIND=4)::ID_total=14
        
        CONTAINS
        
            PROCEDURE:: allocate_data => allocate_data_procedure
            PROCEDURE:: get_pointer => get_pointer_procedure
            PROCEDURE:: deallocate_data => deallocate_data_procedure
        
    END TYPE macro_GP_DATA
    
    PRIVATE:: read_data_procedure,get_pointer_procedure,deallocate_data_procedure
    
CONTAINS
    
    SUBROUTINE allocate_data_procedure(GP_DATA,GP_ID,para,analysis)
    
    !bundeld macro GP DATA (UGLOBAL, STATEV...)
    CLASS(macro_GP_DATA):: GP_DATA
    !ID of the macro Gausspoint
    INTEGER:: GP_ID
    !this type bundles all mesh,convergence & monolithic/staggered parameters
    TYPE(meshparameters)::para
    !this type bundles all convergence/analysis parameters
    TYPE (analysisparameters)::analysis
    
    CALL SMAFloatArrayCreateFortran2D(GP_DATA%STATEV_t,GP_ID*GP_DATA%ID_total-ID_STATEV_t,[para%n_STATEV_elem,para%nElem],0.0_AbqRK)
                                    
    CALL SMAFloatArrayCreateFortran2D(GP_DATA%STATEV_t_1,GP_ID*GP_DATA%ID_total-GP_DATA%ID_STATEV_t_1,[para%n_STATEV_elem,para%nElem],0.0_AbqRK)

    CALL SMAFloatArrayCreateFortran(GP_DATA%UGLOBAL_t,GP_ID*GP_DATA%ID_total-GP_DATA%ID_UGLOBAL_t,para%dimens*para%nNodes,0.0_AbqRK)
                                          
    CALL SMAFloatArrayCreateFortran(GP_DATA%UGLOBAL_t_1,GP_ID*GP_DATA%ID_total-GP_DATA%ID_UGLOBAL_t_1,para%dimens*para%nNodes,0.0_AbqRK)
   
    CALL SMAFloatArrayCreateFortran(GP_DATA%UGLOBAL_old,GP_ID*GP_DATA%ID_total-GP_DATA%ID_UGLOBAL_old,para%dimens*para%nNodes,0.0_AbqRK)
                                          
    IF (analysis%monolithic) THEN !in monolithic case
    
        CALL SMAFloatArrayCreateFortran(GP_DATA%Grad_U_k,GP_ID*GP_DATA%ID_total-GP_DATA%ID_Grad_U_k,para%dimens**2,0.0_AbqRK)
                                                       
        CALL SMAFloatArrayCreateFortran(GP_DATA%Grad_U_t,GP_ID*GP_DATA%ID_total-GP_DATA%ID_Grad_U_t,para%ndof_macro,0.0_AbqRK)
                                                       
        CALL SMAFloatArrayCreateFortran(GP_DATA%Grad_U_old,GP_ID*GP_DATA%ID_total-GP_DATA%ID_Grad_U_old,para%ndof_macro,0.0_AbqRK)
                                                               
        CALL SMAFloatArrayCreateFortran(GP_DATA%rhs_total_k,GP_ID*GP_DATA%ID_total-GP_DATA%ID_rhs_total_k,para%nRows_reduced,0.0_AbqRK)
        
        CALL SMAFloatArrayCreateFortran2D(GP_DATA%d_r_d_F_k,GP_ID*GP_DATA%ID_total-GP_DATA%ID_d_r_d_F_k,[para%nRows_reduced,para%ndof_macro],0.0_AbqRK)
                                  
        CALL SMAIntArrayCreateFortran(GP_DATA%step_indicator,GP_ID*GP_DATA%ID_total-GP_DATA%ID_step_indicator,1,-1)
                
        IF (analysis%save_soe) THEN !if stiffnessmatrix factorization is being stored
        
            CALL SMAIntArrayCreateFortran_dp(GP_DATA%handle_nbr,GP_ID*GP_DATA%ID_total-GP_DATA%ID_handle_nbr,64,0)
        
        END IF 
        
        CALL SMAFloatArrayCreateFortran(GP_DATA%k_matrix_values_k,GP_ID*GP_DATA%ID_total-GP_DATA%ID_k_matrix_values_k,para%nNonZero_reduced,0.0_AbqRK)
    
    END IF
    
    END SUBROUTINE allocate_data_procedure
    
    
    SUBROUTINE get_pointer_procedure(GP_DATA,GP_ID,para,analysis)
    
    !bundeld macro GP DATA (UGLOBAL, STATEV...)
    CLASS(macro_GP_DATA):: GP_DATA
    !ID of the macro Gausspoint
    INTEGER:: GP_ID
    !this type bundles all mesh,convergence & monolithic/staggered parameters
    TYPE(meshparameters)::para
    !this type bundles all convergence/analysis parameters
    TYPE (analysisparameters)::analysis
    
    CALL SMAFloatArrayAccessFortran2D(GP_DATA%STATEV_t,GP_ID*GP_DATA%ID_total-ID_STATEV_t,[para%n_STATEV_elem,para%nElem])
                                    
    CALL SMAFloatArrayAccessFortran2D(GP_DATA%STATEV_t_1,GP_ID*GP_DATA%ID_total-GP_DATA%ID_STATEV_t_1,[para%n_STATEV_elem,para%nElem])

    CALL SMAFloatArrayAccessFortran1D(GP_DATA%UGLOBAL_t,GP_ID*GP_DATA%ID_total-GP_DATA%ID_UGLOBAL_t,para%dimens*para%nNodes)
                                          
    CALL SMAFloatArrayAccessFortran1D(GP_DATA%UGLOBAL_t_1,GP_ID*GP_DATA%ID_total-GP_DATA%ID_UGLOBAL_t_1,para%dimens*para%nNodes)
   
    CALL SMAFloatArrayAccessFortran1D(GP_DATA%UGLOBAL_old,GP_ID*GP_DATA%ID_total-GP_DATA%ID_UGLOBAL_old,para%dimens*para%nNodes)
                                          
    IF (analysis%monolithic) THEN !in monolithic case
    
        CALL SMAFloatArrayAccessFortran1D(GP_DATA%Grad_U_k,GP_ID*GP_DATA%ID_total-GP_DATA%ID_Grad_U_k,para%ndof_macro)
                                                       
        CALL SMAFloatArrayAccessFortran1D(GP_DATA%Grad_U_t,GP_ID*GP_DATA%ID_total-GP_DATA%ID_Grad_U_t,para%ndof_macro)
                                                       
        CALL SMAFloatArrayAccessFortran1D(GP_DATA%Grad_U_old,GP_ID*GP_DATA%ID_total-GP_DATA%ID_Grad_U_old,para%ndof_macro)
                                                               
        CALL SMAFloatArrayAccessFortran1D(GP_DATA%rhs_total_k,GP_ID*GP_DATA%ID_total-GP_DATA%ID_rhs_total_k,para%nRows_reduced)
        
        CALL SMAFloatArrayAccessFortran2D(GP_DATA%d_r_d_F_k,GP_ID*GP_DATA%ID_total-GP_DATA%ID_d_r_d_F_k,[para%nRows_reduced,para%ndof_macro])
                                                
        CALL SMAIntArrayAccessFortran1D(GP_DATA%step_indicator,GP_ID*GP_DATA%ID_total-GP_DATA%ID_step_indicator,1)
               
        IF (analysis%save_soe) THEN
        
            CALL SMAIntArrayAccessFortran1D_dp(GP_DATA%handle_nbr,GP_ID*GP_DATA%ID_total-GP_DATA%ID_handle_nbr,64)
        
        END IF
        
        CALL SMAFloatArrayAccessFortran1D(GP_DATA%k_matrix_values_k,GP_ID*GP_DATA%ID_total-GP_DATA%ID_k_matrix_values_k,para%nNonZero_reduced)
        

    END IF
    
    END SUBROUTINE get_pointer_procedure
    
    
    SUBROUTINE deallocate_data_procedure(GP_DATA,analysis,GP_ID)
    
    !bundeld macro GP DATA (UGLOBAL, STATEV...)
    CLASS(macro_GP_DATA):: GP_DATA
    !this type bundles all convergence/analysis parameters
    TYPE (analysisparameters)::analysis
    !ID of the macro Gausspoint
    INTEGER:: GP_ID
    
    CALL SMAFloatArrayDeleteFortran(GP_ID*GP_DATA%ID_total-ID_STATEV_t)
    CALL SMAFloatArrayDeleteFortran(GP_ID*GP_DATA%ID_total-GP_DATA%ID_STATEV_t_1)
    CALL SMAFloatArrayDeleteFortran(GP_ID*GP_DATA%ID_total-GP_DATA%ID_UGLOBAL_t)
    CALL SMAFloatArrayDeleteFortran(GP_ID*GP_DATA%ID_total-GP_DATA%ID_UGLOBAL_t_1)
    CALL SMAFloatArrayDeleteFortran(GP_ID*GP_DATA%ID_total-GP_DATA%ID_UGLOBAL_old)
    
    IF (analysis%monolithic) THEN !in monolithic case
        CALL SMAFloatArrayDeleteFortran(GP_ID*GP_DATA%ID_total-GP_DATA%ID_Grad_U_k)   
        CALL SMAFloatArrayDeleteFortran(GP_ID*GP_DATA%ID_total-GP_DATA%ID_Grad_U_t)
        CALL SMAFloatArrayDeleteFortran(GP_ID*GP_DATA%ID_total-GP_DATA%ID_Grad_U_old)
        CALL SMAFloatArrayDeleteFortran(GP_ID*GP_DATA%ID_total-GP_DATA%ID_rhs_total_k)
        CALL SMAFloatArrayDeleteFortran(GP_ID*GP_DATA%ID_total-GP_DATA%ID_d_r_d_F_k)
        CALL SMAIntArrayDeleteFortran(GP_ID*GP_DATA%ID_total-GP_DATA%ID_step_indicator)
            
        IF (analysis%save_soe) THEN
            CALL SMAIntArrayDeleteFortran_dp(GP_ID*GP_DATA%ID_total-GP_DATA%ID_handle_nbr)
        END IF
        
        CALL SMAFloatArrayDeleteFortran(GP_ID*GP_DATA%ID_total-GP_DATA%ID_k_matrix_values_k)

    END IF
    
    END SUBROUTINE deallocate_data_procedure

END MODULE type_macro_GP_DATA
