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

MODULE type_macro_GP_DATA
!define a type named macro_GP_DATA which bundels all the state variables 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 freed by "deallocate_GP_DATA"

    USE ABQINTERFACE
    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:: SVARS_t
        INTEGER(KIND=4):: ID_SVARS_t=0
        !internal state variables for all elements at time t+1
        REAL(KIND=AbqRK),DIMENSION(:),POINTER:: SVARS_t_1
        INTEGER(KIND=4):: ID_SVARS_t_1=1
        !micro nodal solution at time t
        REAL(KIND=AbqRK),DIMENSION(:),POINTER:: nodal_sol_t
        INTEGER(KIND=4):: ID_nodal_sol_t=2
        !micro nodal solution at time t+1
        REAL(KIND=AbqRK),DIMENSION(:),POINTER:: nodal_sol_t_1
        INTEGER(KIND=4):: ID_nodal_sol_t_1=3
        !micro nodal solution at time t-1 for extrapolation
        REAL(KIND=AbqRK),DIMENSION(:),POINTER:: nodal_sol_old
        INTEGER(KIND=4):: ID_nodal_sol_old=4
        
        !----------only to be saved in monolithic case-------------------------
        !right hand side of last increment
        REAL(KIND=AbqRK),DIMENSION(:,:),POINTER::rhs_k
        INTEGER(KIND=4):: ID_rhs_k=5
        !partial derivative of reduced micro rhs w.r.t macro measures of last increment
        REAL(KIND=AbqRK),DIMENSION(:,:),POINTER::d_rhs_d_macro_k
        INTEGER(KIND=4):: ID_d_rhs_d_macro_k=6
        !values in CSR Format of micro derivatives of reduced rhs w.r.t. micro measures of last increment
        REAL(KIND=AbqRK),DIMENSION(:,:),POINTER::k_matrix_values_k
        INTEGER(KIND=4):: ID_k_matrix_values_k=7
        !-1 signals that the analysis began;1 signals the beginning of a timestep;
        !2-i signals the middle of a timestept
        INTEGER,POINTER::step_indicator
        INTEGER(KIND=4):: ID_step_indicator=8
        !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=9,ID_solver_real=10
        !element wise training matrix for the hyperintegration
        REAL(KIND=AbqRK),DIMENSION(:,:),POINTER::training_matrix
        INTEGER(KIND=4):: ID_training_matrix=11
        !---------------------------------------------------------------------
        
        !integer with DTIME of old increment (needed for displ. extrapolation)
        REAL(KIND=AbqRK),POINTER:: DTIME_old_increment
        INTEGER(KIND=4):: ID_DTIME_old_increment=12
        
        !------------------------------------------------------------------------
        !number of already outputed snapshots in training
        INTEGER,POINTER::i_data_dump
        INTEGER(KIND=4):: ID_data_dump=13

        
        INTEGER(KIND=4)::ID_total=14
        
        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 (nodal_sol, SVARS...)
    CLASS(macro_GP_DATA):: GP_DATA
    !this type bundles all mesh parameters
    TYPE(meshparameters)::para
    !this type bundles all convergence/analysis parameters
    TYPE (analysisparameters)::analysis
    INTEGER,DIMENSION(:),POINTER::int_pointer
    REAL(KIND=AbqRK),DIMENSION(:),POINTER:: real_pointer
    INTEGER:: n_DOF
    
    IF (analysis%ROM_projection) THEN
        n_DOF=para%n_ROM_modes+para%n_Reaction_force_dof+para%n_additional_dof
    ELSE
        n_DOF=para%ndof_n_max*para%nNodes
    END IF
    
    CALL AbaqusArrayCreate(GP_DATA%SVARS_t,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_SVARS_t,[para%NSVARS_total],0.0_AbqRK)
                                    
    CALL AbaqusArrayCreate(GP_DATA%SVARS_t_1,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_SVARS_t_1,[para%NSVARS_total],0.0_AbqRK)

    CALL AbaqusArrayCreate(GP_DATA%nodal_sol_t,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_nodal_sol_t,[n_DOF],0.0_AbqRK)
                                          
    CALL AbaqusArrayCreate(GP_DATA%nodal_sol_t_1,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_nodal_sol_t_1,[n_DOF],0.0_AbqRK)
   
    CALL AbaqusArrayCreate(GP_DATA%nodal_sol_old,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_nodal_sol_old,[n_DOF],0.0_AbqRK)
    
    CALL AbaqusArrayCreate(int_pointer,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_step_indicator,[1],-1)
    GP_DATA%step_indicator=>int_pointer(1)
    
    IF (analysis%monolithic) THEN !in monolithic case
        
        IF (.NOT. analysis%ROM_projection) THEN
            CALL AbaqusArrayCreate(GP_DATA%rhs_k,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_rhs_k,[para%nRows_reduced,1],0.0_AbqRK)
            CALL AbaqusArrayCreate(GP_DATA%d_rhs_d_macro_k,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_d_rhs_d_macro_k,[para%nRows_reduced,para%n_Reaction_force_dof],0.0_AbqRK)
            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,1],0.0_AbqRK)
        ELSE
            CALL AbaqusArrayCreate(GP_DATA%rhs_k,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_rhs_k,[para%n_ROM_modes,1],0.0_AbqRK)
            CALL AbaqusArrayCreate(GP_DATA%d_rhs_d_macro_k,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_d_rhs_d_macro_k,[para%n_ROM_modes,para%n_Reaction_force_dof],0.0_AbqRK)
        END IF
        
    END IF
    
    IF (analysis%training) THEN
        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)
        IF (analysis%ROM_projection) THEN
            CALL AbaqusArrayCreate(GP_DATA%training_matrix,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_training_matrix,[para%n_active_elements,para%n_ROM_modes+para%n_Reaction_force_dof+1+para%n_additional_hyper_outputs],0.0_AbqRK)
        END IF
    END IF
    
    !Create Array which saves DTIME of the old increment for extrapolation
    CALL AbaqusArrayCreate(real_pointer,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_DTIME_old_increment,[1],0.0_AbqRK)
    GP_DATA%DTIME_old_increment=>real_pointer(1)
        
    END SUBROUTINE allocate_data_procedure
    
    
    SUBROUTINE allocate_solver_data_procedure(GP_DATA,n_reals,n_ints)
    
    !bundeld macro GP DATA (nodal_sol, SVARS...)
    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 (nodal_sol, SVARS...)
    CLASS(macro_GP_DATA):: GP_DATA
    !this type bundles all mesh parameters
    TYPE(meshparameters)::para
    !this type bundles all convergence/analysis parameters
    TYPE (analysisparameters)::analysis
    INTEGER,DIMENSION(:),POINTER::int_pointer
    REAL(KIND=AbqRK),DIMENSION(:),POINTER:: real_pointer
    INTEGER:: n_DOF
    
    IF (analysis%ROM_projection) THEN
        n_DOF=para%n_ROM_modes+para%n_Reaction_force_dof+para%n_additional_dof
    ELSE
        n_DOF=para%ndof_n_max*para%nNodes
    END IF
    
    CALL AbaqusArrayAccess(GP_DATA%SVARS_t,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_SVARS_t,[para%NSVARS_total])
                                    
    CALL AbaqusArrayAccess(GP_DATA%SVARS_t_1,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_SVARS_t_1,[para%NSVARS_total])

    CALL AbaqusArrayAccess(GP_DATA%nodal_sol_t,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_nodal_sol_t,[n_DOF])
                                          
    CALL AbaqusArrayAccess(GP_DATA%nodal_sol_t_1,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_nodal_sol_t_1,[n_DOF])
   
    CALL AbaqusArrayAccess(GP_DATA%nodal_sol_old,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_nodal_sol_old,[n_DOF])
    
    CALL AbaqusArrayAccess(int_pointer,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_step_indicator,[1])
    GP_DATA%step_indicator=>int_pointer(1)
    
    IF (analysis%monolithic) THEN !in monolithic case
        
        IF (.NOT. analysis%ROM_projection) THEN !only in a full simulation save the stiffnessmatrix values
            CALL AbaqusArrayAccess(GP_DATA%rhs_k,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_rhs_k,[para%nRows_reduced,1])
            CALL AbaqusArrayAccess(GP_DATA%d_rhs_d_macro_k,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_d_rhs_d_macro_k,[para%nRows_reduced,para%n_Reaction_force_dof])
            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,1])
        ELSE
            CALL AbaqusArrayAccess(GP_DATA%rhs_k,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_rhs_k,[para%n_ROM_modes,1])
            CALL AbaqusArrayAccess(GP_DATA%d_rhs_d_macro_k,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_d_rhs_d_macro_k,[para%n_ROM_modes,para%n_Reaction_force_dof])
        END IF
        
        IF (analysis%save_soe) 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) THEN
        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)
        IF (analysis%ROM_projection) THEN
            CALL AbaqusArrayAccess(GP_DATA%training_matrix,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_training_matrix,[para%n_active_elements,para%n_ROM_modes+para%n_Reaction_force_dof+1+para%n_additional_hyper_outputs])
        END IF
    END IF
    
    CALL AbaqusArrayAccess(real_pointer,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_DTIME_old_increment,[1])
    GP_DATA%DTIME_old_increment=>real_pointer(1)
        
    END SUBROUTINE get_pointer_procedure
    
    
    SUBROUTINE deallocate_data_procedure(GP_DATA,analysis)
    
    !bundeld macro GP DATA (nodal_sol, SVARS...)
    CLASS(macro_GP_DATA):: GP_DATA
    !this type bundles all convergence/analysis parameters
    TYPE (analysisparameters)::analysis
    INTEGER,DIMENSION(:),POINTER::int_pointer
    REAL(KIND=AbqRK),DIMENSION(:),POINTER:: real_pointer
    
    CALL AbaqusArrayDelete(GP_DATA%SVARS_t,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_SVARS_t)
    CALL AbaqusArrayDelete(GP_DATA%SVARS_t,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_SVARS_t_1)
    CALL AbaqusArrayDelete(GP_DATA%nodal_sol_t,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_nodal_sol_t)
    CALL AbaqusArrayDelete(GP_DATA%nodal_sol_t_1,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_nodal_sol_t_1)
    CALL AbaqusArrayDelete(GP_DATA%nodal_sol_old,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_nodal_sol_old)
    CALL AbaqusArrayDelete(int_pointer,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_step_indicator)
    
    IF (analysis%monolithic) THEN !in monolithic case
        CALL AbaqusArrayDelete(GP_DATA%rhs_k,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_rhs_k)
        CALL AbaqusArrayDelete(GP_DATA%d_rhs_d_macro_k,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_d_rhs_d_macro_k)
        
        IF (.NOT. analysis%ROM_projection) 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) 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
    
    IF (analysis%training) THEN
        CALL AbaqusArrayDelete(int_pointer,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_data_dump)
        IF (analysis%ROM_projection) THEN
            CALL AbaqusArrayDelete(GP_DATA%training_matrix,GP_DATA%GP_ID*GP_DATA%ID_total-GP_DATA%ID_training_matrix)
        END IF
    END IF
    
    CALL AbaqusArrayDelete(real_pointer,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
