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

    !this module defines an object, which comprises all parameters, which define the solution
    !procedure e.g. if staggered or monolithic algorithm will be used, convergence parameters
    !etc.

    USE ABQinterface
    
    PRIVATE::To_upper
    PUBLIC

    TYPE analysisparameters
    
        !-------------------analysis parameters--------------------------------
        !bool variable: .TRUE. IF the monolithic scheme is enforced
        LOGICAL:: monolithic
        !bool variable: .TRUE. IF in monolithic CASE factorized matrix is stored
        !for next increment
        LOGICAL:: save_soe
        !.TRUE. IF the stiffnessmatrix is assumed to be indefinite, .FALSE. IF
        !it is assumed to be positive definite
        LOGICAL:: indefinite_matrix
        !true IF stiffness matrix is symmetric
        LOGICAL:: symmetric_matrix
        !states if the simulation is generating training data for ROM
        LOGICAL:: training
        !.TRUE. IF full simulation, .FALSE. IF reduced ROM projection is applied
        LOGICAL:: ROM_projection
        !.TRUE. IF full integration, .FALSE. IF the problem is element-hyper integration
        LOGICAL:: hyperintegration
        !max number of NR steps in staggered scheme
        INTEGER:: max_iters
        !maximum number of macro gausspoints
        INTEGER:: max_n_macro_GP=27
        INTEGER:: max_n_GP_macro_total=1000000
        INTEGER(KIND=4):: ID_ints=1
        !multiplication factor for NR abort criterion in staggered scheme
        REAL(KIND=AbqRK):: convergence_ratio
        INTEGER(KIND=4):: ID_reals=2
        !array with information -> number of existing RVE definitions
        INTEGER,POINTER:: active_RVE_definitions
        INTEGER(KIND=4):: ID_active_RVE_definitions=3
        !array with information -> macro GP exists and the RVE definition label
        !for the i-th macro GP is found in active_gausspoints(i)
        INTEGER, DIMENSION(:),POINTER:: active_gausspoints
        INTEGER(KIND=4):: ID_active_gausspoints=4
        !array with the information at which times to dump data for training
        REAL(KIND=AbqRK),DIMENSION(:),POINTER:: data_dump
        INTEGER(KIND=4):: ID_data_dump=5
        INTEGER:: n_data_dumps
        
        CONTAINS
        
            PROCEDURE:: read_data => read_data_procedure
            PROCEDURE:: get_pointer => get_pointer_procedure
            PROCEDURE:: deallocate_data => deallocate_data_procedure

    END TYPE analysisparameters
    
    PRIVATE::read_data_procedure,get_pointer_procedure,deallocate_data_procedure
    
CONTAINS

SUBROUTINE read_data_procedure(analysis)

    IMPLICIT NONE
    
    CLASS(analysisparameters)::analysis
    !Name of directory in which Abaqus is stared
    CHARACTER(LEN=256)::  OUTDIR
    !lenght of name of this directory
    INTEGER:: LENOUTDIR
    INTEGER(KIND=4), DIMENSION(:),POINTER:: array_int
    REAL(KIND=AbqRK), DIMENSION(:),POINTER:: array_real
    INTEGER:: io_error, error , read_error, unit_numer
    CHARACTER(len=500),DIMENSION(2):: line
    CHARACTER(len=100),DIMENSION(10):: Keywords
    CHARACTER(len=100),DIMENSION(4),PARAMETER:: NecessaryKeywords=['*ALGORITHM','*SYSTEMMATRIX_DEFINITENESS','*SYSTEMMATRIX_SYMMETRY','*SOLVING_PROCESS']
    CHARACTER(len=100),DIMENSION(1),PARAMETER:: NecessaryKeywords_monolithic=['*FACTORIZATION']
    CHARACTER(len=100),DIMENSION(2),PARAMETER:: NecessaryKeywords_staggered=['*EQUILIBRIUM_ITERATIONS','*CONVERGENCE_RATIO']
    CHARACTER(len=100),DIMENSION(1),PARAMETER:: NecessaryKeywords_training=['*DATA_DUMP']
    INTEGER:: i,k
    
    !create array active_gausspoints
    CALL AbaqusArrayCreate(analysis%active_gausspoints,analysis%ID_active_gausspoints,[analysis%max_n_GP_macro_total],0)    
    !create array active_RVE_definitions
    CALL AbaqusArrayCreate(array_int,analysis%ID_active_RVE_definitions,[1],0)
    analysis%active_RVE_definitions=>array_int(1)
    !Integer Analysisparameters
    CALL AbaqusArrayCreate(array_int,analysis%ID_ints,[9],0)
    !Real Analysisparameters
    CALL AbaqusArrayCreate(array_real,analysis%ID_reals,[1],0.0_AbqRK)
    
    !initialize
    analysis%training=.FALSE.
    analysis%ROM_projection=.FALSE.
    
    !open File named FE2_Analysisparameters.cfg with analyisis parameters e.g.
    !convergence_ratio, IF not supplied, THEN default parameters are enforced
    
    CALL GETOUTDIR( OUTDIR, LENOUTDIR )!get directory name and length
    unit_numer=101
    OPEN(unit=unit_numer,file=OUTDIR(1:LENOUTDIR)//'/'//'FE2_Analysisparameters.cfg',&
             status='old',action='READ',iostat = io_error) !open the file
    
    Keywords=''
    
    IF (io_error==0) THEN !read parameters from file
    
        read_error=0 !reports whether there was any error during the READ process
        
        k=1 !running index
        
        DO WHILE (.TRUE.)

            READ(unit_numer,*,iostat=error) line
            
            IF (error/=0) THEN
                EXIT
            ELSE
                CALL To_upper(line(1)) !convert to upper cases
                
                IF (line(1)(1:1)=='*' .AND. line(1)(2:2)/='*') THEN
                                
                    Keywords(k)=trim(line(1))
                    
                    SELECT CASE(Keywords(k)) !select the current statement
                    CASE('*ALGORITHM')
                        CALL To_upper(line(2))
                        IF (line(2)=='STAGGERED') THEN
                            analysis%monolithic=.FALSE.
                        ELSE IF (line(2)=='MONOLITHIC') THEN
                            analysis%monolithic=.TRUE.
                        ELSE
                            read_error=1
                        END IF
                    CASE('*FACTORIZATION')
                        CALL To_upper(line(2))
                        IF (line(2)=='SAVED') THEN
                            analysis%save_soe=.TRUE.
                        ELSE IF (line(2)=='NOTSAVED') THEN
                            analysis%save_soe=.FALSE.
                        ELSE
                            read_error=1
                        END IF
                    CASE('*EQUILIBRIUM_ITERATIONS')
                        READ(line(2),*) analysis%max_iters
                        IF (analysis%max_iters<1)  read_error=1
                    CASE('*CONVERGENCE_RATIO')
                        READ(line(2),*) analysis%convergence_ratio
                        IF ((analysis%convergence_ratio>1.0) .OR. (analysis%convergence_ratio<=0.0) ) read_error=1
                    CASE('*SYSTEMMATRIX_DEFINITENESS')
                        CALL To_upper(line(2))
                        IF (line(2)=='INDEFINITE') THEN
                            analysis%indefinite_matrix=.TRUE.
                        ELSE IF (line(2)=='POSITIV') THEN
                            analysis%indefinite_matrix=.FALSE.
                        ELSE
                            read_error=1
                        END IF
                    CASE('*SYSTEMMATRIX_SYMMETRY')
                        CALL To_upper(line(2))
                        IF (line(2)=='SYMM') THEN
                            analysis%symmetric_matrix=.TRUE.
                        ELSE IF (line(2)=='UNSYMM') THEN
                            analysis%symmetric_matrix=.FALSE.
                        ELSE
                            read_error=1
                        END IF
                    CASE('*TRAINING_DATA')
                        CALL To_upper(line(2))
                        IF (line(2)=='ROM') THEN
                            analysis%training=.TRUE.
                        ELSE IF (line(2)=='NO') THEN
                            analysis%training=.FALSE.
                        ELSE
                            read_error=1
                        END IF
                    CASE('*SOLVING_PROCESS')
                        CALL To_upper(line(2))
                        IF (line(2)=='FULL') THEN
                            analysis%ROM_projection=.FALSE.
                            analysis%hyperintegration=.FALSE.
                        ELSE IF (line(2)=='REDUCED') THEN
                            analysis%ROM_projection=1
                            analysis%hyperintegration=.FALSE.
                        ELSE IF (line(2)=='HYPERREDUCED') THEN
                            analysis%ROM_projection=1
                            analysis%hyperintegration=1
                        ELSE
                            read_error=1
                        END IF
                    CASE('*DATA_DUMP')
                        IF (analysis%training) THEN
                            READ(line(2)(3:),*) analysis%n_data_dumps
                            CALL AbaqusArrayCreate(analysis%data_dump,analysis%ID_data_dump,[analysis%n_data_dumps],0.0_AbqRK)
                            READ(unit_numer,*,iostat=error) analysis%data_dump
                            DO i=2,analysis%n_data_dumps
                                IF (analysis%data_dump(i)<analysis%data_dump(i-1)) THEN
                                    CALL STDB_ABQERR(-3_AbqIK,  'Data dump timesteps must be in ascending order!', 0_AbqIK , 0.0_AbqRK, ' ')
                                END IF
                            END DO
                        ELSE
                            READ(unit_numer,*) !skip this line
                        END IF
                    END SELECT
                    k=k+1
                END IF
            END IF
        END DO
        
        !check IF there has been an error in reading the file
        DO i=1,SIZE(NecessaryKeywords)
            IF (NOT(ANY(Keywords==NecessaryKeywords(i)))) read_error=1
        END DO
        IF (analysis%monolithic==0) THEN
            DO i=1,SIZE(NecessaryKeywords_staggered)
                IF (NOT(ANY(Keywords==NecessaryKeywords_staggered(i)))) read_error=1
            END DO
        ELSE
            DO i=1,SIZE(NecessaryKeywords_monolithic)
                IF (NOT(ANY(Keywords==NecessaryKeywords_monolithic(i)))) read_error=1
            END DO
        END IF
        IF (analysis%training) THEN
            DO i=1,SIZE(NecessaryKeywords_training)
                IF (NOT(ANY(Keywords==NecessaryKeywords_training(i)))) read_error=1
            END DO
        END IF
        
        !display an error and stop the simulation IF an error has been reported
        IF (read_error==1) CALL STDB_ABQERR(-3_AbqIK,  'An error in the FE2_Analysisparameters.cfg'//&
                                      ' file has been found.', 0_AbqIK , 0.0_AbqRK, ' ')
        
            
    ELSE !set default parameters

        !default parameters: monolithic scheme without saving the factorization,
        !indefinite matrix,symmetric matrix
        analysis%monolithic=.TRUE.
        analysis%save_soe=.FALSE.
        analysis%indefinite_matrix=.TRUE.
        analysis%symmetric_matrix=.TRUE.
        analysis%training=.FALSE.
        analysis%ROM_projection=.FALSE.
        analysis%hyperintegration=.FALSE.
    
    END IF
    
    CLOSE (unit=unit_numer) !close file
    
    !Always save the factorization in monolithic ROM simulations, because the memory use is neglegible
    IF (analysis%ROM_projection) analysis%save_soe=.TRUE.

    !assigne values
    array_int(1)  =  analysis%monolithic
    array_int(2)  =  analysis%save_soe
    array_int(3)  =  analysis%max_iters
    array_int(4)  =  analysis%indefinite_matrix
    array_int(5)  =  analysis%symmetric_matrix
    array_int(6)  =  analysis%training
    array_int(7)  =  analysis%ROM_projection
    array_int(8)  =  analysis%hyperintegration
    array_int(9)  =  analysis%n_data_dumps
    array_real(1) =  analysis%convergence_ratio
    
    !print information to message file
    CALL STDB_ABQERR(1_AbqIK,  'The following analysis parameters are being used:', 0_AbqIK , 0.0_AbqRK, ' ')
    IF (analysis%monolithic .AND. analysis%save_soe) THEN
        CALL STDB_ABQERR(1_AbqIK,  'Monolithic solution scheme with saving the sytem matrix factorization.', 0_AbqIK , 0.0_AbqRK, ' ')
    ELSE IF (analysis%monolithic .AND. (.NOT. analysis%save_soe)) THEN
        CALL STDB_ABQERR(1_AbqIK,  'Monolithic solution scheme without saving the sytem matrix factorization.', 0_AbqIK , 0.0_AbqRK, ' ')
    ELSE
        CALL STDB_ABQERR(1_AbqIK,  'Staggered solution scheme with max. %I iterations and a residual force criterion of of %R.', analysis%max_iters , analysis%convergence_ratio, ' ')
    END IF
    IF (analysis%indefinite_matrix .AND. analysis%symmetric_matrix) THEN
        CALL STDB_ABQERR(1_AbqIK,  'The system matrix is assumed to be symmetric and indefinite.', 0_AbqIK , 0.0_AbqRK, ' ')
    ELSE IF (analysis%indefinite_matrix .AND. (.NOT. analysis%symmetric_matrix)) THEN
        CALL STDB_ABQERR(1_AbqIK,  'The system matrix is assumed to be unsymmetric and indefinite.', 0_AbqIK , 0.0_AbqRK, ' ')
    ELSE IF ((.NOT. analysis%indefinite_matrix) .AND. analysis%symmetric_matrix) THEN
        CALL STDB_ABQERR(1_AbqIK,  'The system matrix is assumed to be symmetric and positiv definite.', 0_AbqIK , 0.0_AbqRK, ' ')
    ELSE
        CALL STDB_ABQERR(1_AbqIK,  'The system matrix is assumed to be unsymmetric and positiv definite.', 0_AbqIK , 0.0_AbqRK, ' ')
    END IF
    IF (analysis%ROM_projection) THEN
        CALL STDB_ABQERR(1_AbqIK,  'The micro problems are solved using the reduced order modeling (ROM) technique.', 0_AbqIK , 0.0_AbqRK, ' ')
        IF (analysis%hyperintegration) THEN
            CALL STDB_ABQERR(1_AbqIK,  'The micro problems are element-wise hyperintegrated.', 0_AbqIK , 0.0_AbqRK, ' ')
        ELSE
            CALL STDB_ABQERR(1_AbqIK,  'The micro problems are integrated using all GPs.', 0_AbqIK , 0.0_AbqRK, ' ')
        END IF
    END IF
    IF (analysis%training) THEN
        CALL STDB_ABQERR(1_AbqIK,  'Training data is being generated and outputed.', 0_AbqIK , 0.0_AbqRK, ' ')
    END IF
        
END SUBROUTINE read_data_procedure


SUBROUTINE get_pointer_procedure(analysis)

    IMPLICIT NONE

    CLASS(analysisparameters)::analysis
    INTEGER(KIND=4), DIMENSION(:),POINTER:: array_int
    REAL(KIND=AbqRK), DIMENSION(:),POINTER:: array_real
    
    !get the pointers to the previously allocated data
    
    !get integer analysis parameters
    CALL AbaqusArrayAccess(array_int,analysis%ID_ints,[9])
    analysis%monolithic         =   array_int(1)
    analysis%save_soe           =   array_int(2)
    analysis%max_iters          =   array_int(3)
    analysis%indefinite_matrix  =   array_int(4)
    analysis%symmetric_matrix   =   array_int(5)
    analysis%training           =   array_int(6)
    analysis%ROM_projection     =   array_int(7)
    analysis%hyperintegration   =   array_int(8)
    analysis%n_data_dumps       =   array_int(9)
    
    
    !get real analysis parameters
    CALL AbaqusArrayAccess(array_real,analysis%ID_reals,[1])
    analysis%convergence_ratio  =   array_real(1)
    
    CALL AbaqusArrayAccess(analysis%active_gausspoints,analysis%ID_active_gausspoints,[analysis%max_n_GP_macro_total])
    CALL AbaqusArrayAccess(array_int,analysis%ID_active_RVE_definitions,[1])
    analysis%active_RVE_definitions=>array_int(1)
    IF (analysis%training) THEN
        CALL AbaqusArrayAccess(analysis%data_dump,analysis%ID_data_dump,[analysis%n_data_dumps])
    END IF
    
END SUBROUTINE get_pointer_procedure


SUBROUTINE deallocate_data_procedure(analysis)

    IMPLICIT NONE

    CLASS(analysisparameters)::analysis
    INTEGER(KIND=4), DIMENSION(:),POINTER:: array_int
    REAL(KIND=AbqRK), DIMENSION(:),POINTER:: array_real
    
    !deallocate the previously allocated data
    
    CALL AbaqusArrayDelete(analysis%active_gausspoints,analysis%ID_active_gausspoints)
    CALL AbaqusArrayDelete(array_int,analysis%ID_active_RVE_definitions)
    CALL AbaqusArrayDelete(array_int,analysis%ID_ints)
    CALL AbaqusArrayDelete(array_real,analysis%ID_reals)
    IF (analysis%training) THEN
        CALL AbaqusArrayDelete(analysis%data_dump,analysis%ID_data_dump)
    END IF

END SUBROUTINE deallocate_data_procedure

SUBROUTINE To_upper(str)
   !this subroutine converts a string with upper and lower cases to a purely
   !upper case expression
   
   IMPLICIT NONE
   
   CHARACTER(*), INTENT(IN OUT) :: str
   INTEGER :: i
    
   DO i = 1, LEN(str)
   SELECT CASE(str(i:i))
       CASE("a":"z")
       str(i:i) = ACHAR(IACHAR(str(i:i))-32)
   END SELECT
   END DO
   
END SUBROUTINE To_upper

END MODULE type_analysisparameters
