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

    !this module defines an object which comprises all global matrices of derivatives,
    !RHS's & solution increments needed temporarily during one global NR step

    USE ABQINTERFACE
    USE type_analysisparameters
    USE type_meshparameters
    USE type_macro_GP_DATA
    
    PUBLIC
    
    TYPE matrix2D
        REAL(KIND=AbqRK),DIMENSION(:,:),POINTER:: matrix
    END TYPE matrix2D
    
    TYPE systemmatrices

        !Matrix of derivatives 1=micro 2=macro -> (1,1) derivative of micro rhs
        !w.r.t micro measures, (1,2) derivative of micro rhs w.r.t macro measures
        !(2,1) derivative of macro response w.r.t micro measures, (2,2) derivative
        !of macro response w.r.t macro measures
        TYPE(matrix2D),DIMENSION(2,2):: derivatives,derivatives_full
        !global rhs -> (1) micro global rhs, (2) macro response
        TYPE(matrix2D),DIMENSION(2):: rhs,rhs_full
        !nodal solution increment per NR loop
        REAL(KIND=AbqRK),DIMENSION(:,:),POINTER:: nodal_sol_incr,nodal_sol_incr_full
        !nodal solution of current and past time increment
        REAL(KIND=AbqRK),DIMENSION(:),POINTER:: nodal_sol_t,nodal_sol_t_1
        
        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 (nodal_sol,SVARS,...)
    TYPE(macro_GP_DATA),INTENT(IN)::GP_DATA
    
    IF (.NOT. analysis%ROM_projection) THEN
        matrices%nodal_sol_t=>GP_DATA%nodal_sol_t
        matrices%nodal_sol_t_1=>GP_DATA%nodal_sol_t_1
        ALLOCATE(matrices%derivatives(1,1)%matrix(para%nNonZero_reduced,1))
        ALLOCATE(matrices%derivatives(2,1)%matrix(para%n_Reaction_force_dof,para%nRows_reduced))
        ALLOCATE(matrices%derivatives(1,2)%matrix(para%nRows_reduced,para%n_Reaction_force_dof))
        ALLOCATE(matrices%derivatives(2,2)%matrix(para%n_Reaction_force_dof,para%n_Reaction_force_dof))
        matrices%derivatives_full(1,1)%matrix=>matrices%derivatives(1,1)%matrix
        matrices%derivatives_full(1,2)%matrix=>matrices%derivatives(1,2)%matrix
        matrices%derivatives_full(2,1)%matrix=>matrices%derivatives(2,1)%matrix
        matrices%derivatives_full(2,2)%matrix=>matrices%derivatives(2,2)%matrix
        ALLOCATE(matrices%rhs(1)%matrix(para%nRows_reduced,1))
        ALLOCATE(matrices%rhs(2)%matrix(para%n_Reaction_force_dof,1))
        matrices%rhs_full(1)%matrix=>matrices%rhs(1)%matrix
        matrices%rhs_full(2)%matrix=>matrices%rhs(2)%matrix
        ALLOCATE(matrices%nodal_sol_incr(para%nRows_reduced,1))
        matrices%nodal_sol_incr_full=>matrices%nodal_sol_incr
    ELSE
        ALLOCATE(matrices%nodal_sol_t(para%ndof_n_max*para%nNodes))
        ALLOCATE(matrices%nodal_sol_t_1(para%ndof_n_max*para%nNodes))
        ALLOCATE(matrices%derivatives(1,1)%matrix(para%n_ROM_modes,para%n_ROM_modes))
        ALLOCATE(matrices%derivatives(2,1)%matrix(para%n_Reaction_force_dof,para%n_ROM_modes))
        ALLOCATE(matrices%derivatives(1,2)%matrix(para%n_ROM_modes,para%n_Reaction_force_dof))
        ALLOCATE(matrices%derivatives(2,2)%matrix(para%n_Reaction_force_dof,para%n_Reaction_force_dof))
        ALLOCATE(matrices%nodal_sol_incr(para%n_ROM_modes,1))
        ALLOCATE(matrices%nodal_sol_incr_full(para%nRows_reduced,1))
        ALLOCATE(matrices%rhs(1)%matrix(para%n_ROM_modes,1))
        ALLOCATE(matrices%rhs(2)%matrix(para%n_Reaction_force_dof,1))
        IF (.NOT. analysis%hyperintegration) THEN
            ALLOCATE(matrices%rhs_full(1)%matrix(para%nRows_reduced,1))
            matrices%rhs_full(2)%matrix=>matrices%rhs(2)%matrix
            ALLOCATE(matrices%derivatives_full(1,1)%matrix(para%nNonZero_reduced,1))
            ALLOCATE(matrices%derivatives_full(2,1)%matrix(para%n_Reaction_force_dof,para%nRows_reduced))
            ALLOCATE(matrices%derivatives_full(1,2)%matrix(para%nRows_reduced,para%n_Reaction_force_dof))
            matrices%derivatives_full(2,2)%matrix=>matrices%derivatives(2,2)%matrix
        END IF
    END IF
    
END SUBROUTINE allocation_procedure

SUBROUTINE deallocation_procedure(matrices,analysis)

    IMPLICIT NONE
    
    CLASS(systemmatrices)::matrices
    TYPE(analysisparameters)::analysis
    
    IF (.NOT. analysis%ROM_projection) THEN
        DEALLOCATE(matrices%derivatives(1,1)%matrix)
        DEALLOCATE(matrices%derivatives(2,1)%matrix)
        DEALLOCATE(matrices%derivatives(1,2)%matrix)
        DEALLOCATE(matrices%derivatives(2,2)%matrix)
        DEALLOCATE(matrices%rhs(1)%matrix)
        DEALLOCATE(matrices%rhs(2)%matrix)
        DEALLOCATE(matrices%nodal_sol_incr)
    ELSE
        DEALLOCATE(matrices%nodal_sol_t)
        DEALLOCATE(matrices%nodal_sol_t_1)
        DEALLOCATE(matrices%derivatives(1,1)%matrix)
        DEALLOCATE(matrices%derivatives(2,1)%matrix)
        DEALLOCATE(matrices%derivatives(1,2)%matrix)
        DEALLOCATE(matrices%derivatives(2,2)%matrix)
        DEALLOCATE(matrices%nodal_sol_incr)
        DEALLOCATE(matrices%nodal_sol_incr_full)
        DEALLOCATE(matrices%rhs(1)%matrix)
        DEALLOCATE(matrices%rhs(2)%matrix)
        IF (.NOT. analysis%hyperintegration) THEN
            DEALLOCATE(matrices%rhs_full(1)%matrix)
            DEALLOCATE(matrices%derivatives_full(1,1)%matrix)
            DEALLOCATE(matrices%derivatives_full(2,1)%matrix)
            DEALLOCATE(matrices%derivatives_full(1,2)%matrix)
        END IF
    END IF
    
END SUBROUTINE deallocation_procedure

END MODULE type_systemmatrices
