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

    USE ABQINTERFACE
    USE ABQSMA
    
    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
        !max number of NR steps in staggered scheme
        INTEGER:: max_iters
        !multiplication factor for NR abort criterion in staggered scheme
        REAL(KIND=AbqRK):: convergence_ratio
        !maximum number of macro gausspoints
        INTEGER:: max_n_macro_GP=27
        INTEGER:: max_n_GP_macro_total=600000
        !array with information -> RVE definition existing
        INTEGER, DIMENSION(:),POINTER:: active_RVE_definitions
        !array with information -> macro GP exists
        INTEGER, DIMENSION(:),POINTER:: active_gausspoints
        !integer with DTIME of old increment (needed for displ. extrapolation)
        REAL(KIND=AbqRK), DIMENSION(:),POINTER:: DTIME_old_increment
        
        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*256  OUTDIR
    !lenght of name of this directory
    INTEGER:: LENOUTDIR
    INTEGER, DIMENSION(:),POINTER:: array_int
    REAL(KIND=AbqRK), DIMENSION(:),POINTER:: array_real
    INTEGER:: io_error, error , read_error
    CHARACTER(len=500),DIMENSION(2):: line
    CHARACTER(len=100),DIMENSION(10):: Keywords
    CHARACTER(len=100),DIMENSION(3),PARAMETER:: NecessaryKeywords=['*ALGORITHM','*SYSTEMMATRIX_DEFINITENESS','*SYSTEMMATRIX_SYMMETRY']
    CHARACTER(len=100),DIMENSION(1),PARAMETER:: NecessaryKeywords_monolithic=['*FACTORIZATION']
    CHARACTER(len=100),DIMENSION(2),PARAMETER:: NecessaryKeywords_staggered=['*EQUILIBRIUM_ITERATIONS','*CONVERGENCE_RATIO']
    INTEGER:: i,k
    
    !create array active_gausspoints
    CALL SMAIntArrayCreateFortran(analysis%active_gausspoints,1,analysis%max_n_GP_macro_total,0)    
    !create array active_RVE_definitions
    CALL SMAIntArrayCreateFortran(analysis%active_RVE_definitions,2,5,0)
    !Create Array which saves DTIME of the old increment for extrapolation
    CALL SMAFloatArrayCreateFortran(analysis%DTIME_old_increment,3,1,0.0_AbqRK)
        
    !Integer Analysisparameters
    CALL SMAIntArrayCreateFortran(array_int,4,5,0)
    !Float Analysisparameters
    CALL SMAFloatArrayCreateFortran(array_real,5,1,0.0_AbqRK)
    
    !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
    OPEN(unit=26,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(26,*,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=0
                        ELSE IF (line(2)=='MONOLITHIC') THEN
                            analysis%monolithic=1
                        ELSE
                            read_error=1
                        END IF
                    CASE('*FACTORIZATION')
                        CALL To_upper(line(2))
                        IF (line(2)=='SAVED') THEN
                            analysis%save_soe=1
                        ELSE IF (line(2)=='NOTSAVED') THEN
                            analysis%save_soe=0
                        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=1
                        ELSE IF (line(2)=='POSITIV') THEN
                            analysis%indefinite_matrix=0
                        ELSE
                            read_error=1
                        END IF
                    CASE('*SYSTEMMATRIX_SYMMETRY')
                        CALL To_upper(line(2))
                        IF (line(2)=='SYMM') THEN
                            analysis%symmetric_matrix=1
                        ELSE IF (line(2)=='UNSYMM') THEN
                            analysis%symmetric_matrix=0
                        ELSE
                            read_error=1
                        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
        
        !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=1
        analysis%save_soe=0
        analysis%max_iters=8
        analysis%indefinite_matrix=1
        analysis%symmetric_matrix=1
        analysis%convergence_ratio=0.00001
    
    END IF
        
    CLOSE (unit=26) !close file

    !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_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
    
END SUBROUTINE read_data_procedure


SUBROUTINE get_pointer_procedure(analysis)

    IMPLICIT NONE

    CLASS(analysisparameters)::analysis
    INTEGER, DIMENSION(:),POINTER:: array_int
    REAL(KIND=AbqRK), DIMENSION(:),POINTER:: array_real
    
    !get the pointers to the previously allocated data
    
    !get integer analysis parameters
    CALL SMAIntArrayAccessFortran1D(array_int,4,5)
    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)
    
    !get real analysis parameters
    CALL SMAFloatArrayAccessFortran1D(array_real,5,1)
    analysis%convergence_ratio=array_real(1)
    
    CALL SMAIntArrayAccessFortran1D(analysis%active_gausspoints,1,analysis%max_n_GP_macro_total)
    CALL SMAIntArrayAccessFortran1D(analysis%active_RVE_definitions,2,5)
    CALL SMAFloatArrayAccessFortran1D(analysis%DTIME_old_increment,3,1)
    
    
END SUBROUTINE get_pointer_procedure


SUBROUTINE deallocate_data_procedure(analysis)

    CLASS(analysisparameters)::analysis
    INTEGER:: rve_number
    
    !deallocate the previously allocated data
    
    CALL SMAIntArrayDeleteFortran(1)   !active_gausspoints
    CALL SMAIntArrayDeleteFortran(2)   !active_RVE_definitions
    CALL SMAFloatArrayDeleteFortran(3) !DTIME_old_increment
    CALL SMAIntArrayDeleteFortran(4)   !integer Analysisparameters
    CALL SMAFloatArrayDeleteFortran(5) !real Analysisparameters
    

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