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

    !this module defines an object which comprises all stiffnessmatrices, RHS's & solution
    !increments needed temporarily during one or more NR steps

    USE ABQINTERFACE
    USE type_analysisparameters
    USE type_meshparameters
    USE type_macro_GP_DATA
    
    PUBLIC
    
    TYPE systemmatrices

        !Stiffnessmatrix values
        REAL(KIND=AbqRK),DIMENSION(:),ALLOCATABLE:: k_matrix_values
        REAL(KIND=AbqRK),DIMENSION(:,:),ALLOCATABLE:: k_matrix
        !total rhs array
        REAL(KIND=AbqRK),DIMENSION(:,:),POINTER:: rhs_total,rhs_total_full
        !partial derivative of macro stress (1.PK) w.r.t micro displacements
        REAL(KIND=AbqRK),DIMENSION(:,:),POINTER:: d_PK1_d_u,d_PK1_d_u_full
        !partial derivative of reaction force w.r.t macro deformation gradient
        REAL(KIND=AbqRK),DIMENSION(:,:),POINTER:: d_r_d_F,d_r_d_F_full
        !partial derivative of macro stress (1. PK) w.r.t macro deformation gradient
        REAL(KIND=AbqRK),DIMENSION(:,:),ALLOCATABLE:: d_PK1_d_F
        !displacement of all nodes of step t and t+1
        REAL(KIND=AbqRK),DIMENSION(:),POINTER:: UGLOBAL_t,UGLOBAL_t_1
        !displacement increment per NR loop
        REAL(KIND=AbqRK),DIMENSION(:,:),POINTER:: DELTAU
        REAL(KIND=AbqRK),DIMENSION(:,:),POINTER:: DELTAU_FULL
        
        CONTAINS
        
            PROCEDURE:: allocation => allocation_procedure
            PROCEDURE:: deallocation => deallocation_procedure

    END TYPE systemmatrices
    
    PRIVATE:: allocation_procedure,deallocation_procedure
    
CONTAINS

SUBROUTINE allocation_procedure(matrices,analysis,para,GP_DATA)

    IMPLICIT NONE
    
    CLASS(systemmatrices)::matrices
    TYPE(analysisparameters)::analysis
    TYPE(meshparameters)::para
    !this type bundles all macro GP DATA (UGLOBAL,STATEV,...)
    TYPE(macro_GP_DATA),INTENT(IN)::GP_DATA
    
    IF (analysis%solving_process==0) THEN
        ALLOCATE(matrices%k_matrix_values(para%nNonZero_reduced))
        matrices%UGLOBAL_t=>GP_DATA%UGLOBAL_t
        matrices%UGLOBAL_t_1=>GP_DATA%UGLOBAL_t_1
        ALLOCATE(matrices%DELTAU(para%nRows_reduced,1))
        matrices%DELTAU_FULL=>matrices%DELTAU
        ALLOCATE(matrices%rhs_total(para%nRows_reduced,1))
        matrices%rhs_total_full=>matrices%rhs_total
        ALLOCATE(matrices%d_PK1_d_u(para%ndof_macro,para%nRows_reduced))
        ALLOCATE(matrices%d_r_d_F(para%nRows_reduced,para%ndof_macro))
        ALLOCATE(matrices%d_PK1_d_F(para%ndof_macro,para%ndof_macro))
        matrices%d_PK1_d_u_full=>matrices%d_PK1_d_u
        matrices%d_r_d_F_full=>matrices%d_r_d_F
    ELSE IF (analysis%solving_process==1) THEN
        ALLOCATE(matrices%k_matrix(para%n_ROM_modes,para%n_ROM_modes))
        ALLOCATE(matrices%UGLOBAL_t(para%ndof_n*para%nNodes))
        ALLOCATE(matrices%UGLOBAL_t_1(para%ndof_n*para%nNodes))
        ALLOCATE(matrices%DELTAU(para%n_ROM_modes,1))
        ALLOCATE(matrices%DELTAU_FULL(para%nRows_reduced,1))
        ALLOCATE(matrices%rhs_total(para%n_ROM_modes,1))
        ALLOCATE(matrices%d_PK1_d_u(para%ndof_macro,para%n_ROM_modes))
        ALLOCATE(matrices%d_r_d_F(para%n_ROM_modes,para%ndof_macro))
        ALLOCATE(matrices%d_PK1_d_F(para%ndof_macro,para%ndof_macro))
        IF (.NOT. analysis%hyperintegration) THEN
            ALLOCATE(matrices%d_PK1_d_u_full(para%ndof_macro,para%nRows_reduced))
            ALLOCATE(matrices%d_r_d_F_full(para%nRows_reduced,para%ndof_macro))
            ALLOCATE(matrices%rhs_total_full(para%nRows_reduced,1))
            ALLOCATE(matrices%k_matrix_values(para%nNonZero_reduced))
        END IF
    END IF
    
END SUBROUTINE allocation_procedure

SUBROUTINE deallocation_procedure(matrices,analysis)

    IMPLICIT NONE
    
    CLASS(systemmatrices)::matrices
    TYPE(analysisparameters)::analysis
    
    IF (analysis%solving_process==0) THEN
       DEALLOCATE(matrices%DELTAU)
       DEALLOCATE(matrices%k_matrix_values)
    ELSE IF (analysis%solving_process==1) THEN
       DEALLOCATE(matrices%k_matrix)
       DEALLOCATE(matrices%UGLOBAL_t)
       DEALLOCATE(matrices%UGLOBAL_t_1)
       DEALLOCATE(matrices%DELTAU_FULL)
       DEALLOCATE(matrices%DELTAU)
       IF (.NOT. analysis%hyperintegration) THEN
            DEALLOCATE(matrices%rhs_total_full)
            DEALLOCATE(matrices%k_matrix_values)
            DEALLOCATE(matrices%d_PK1_d_u_full)
            DEALLOCATE(matrices%d_r_d_F_full)
        END IF
    END IF
    DEALLOCATE(matrices%d_PK1_d_u)
    DEALLOCATE(matrices%d_r_d_F)
    DEALLOCATE(matrices%d_PK1_d_F)
    DEALLOCATE(matrices%rhs_total)
    
END SUBROUTINE deallocation_procedure

END MODULE type_systemmatrices
