
MODULE linear_interpolation

    INTERFACE linear_interpol
        MODULE PROCEDURE linear_interpol_1D
        MODULE PROCEDURE linear_interpol_2D
    END INTERFACE linear_interpol
    
    CONTAINS

        SUBROUTINE linear_interpol_1D(t,timesteps,interpolant,res)
            !interpols linearily within a interpolant matrix
            
            USE ABQinterface
            
            IMPLICIT NONE
            
            REAL(KIND=AbqRK),DIMENSION(:),INTENT(IN):: timesteps
            REAL(KIND=AbqRK),DIMENSION(:,:),INTENT(IN):: interpolant
            REAL(KIND=AbqRK),DIMENSION(:),INTENT(OUT):: res
            REAL(KIND=AbqRK),INTENT(IN):: t
            INTEGER:: i
            
            IF (t<timesteps(1)) THEN
                res=t/timesteps(1)*interpolant(:,1)
            ELSE IF (SIZE(timesteps)==1) THEN
                res=interpolant(:,1)
            ELSE
                DO i=2,size(timesteps)
                    IF (t<=timesteps(i)) THEN
                        res=interpolant(:,i-1)+(interpolant(:,i)-interpolant(:,i-1))*(t-timesteps(i-1))/(timesteps(i)-timesteps(i-1))
                        EXIT
                    ELSE IF ((i==size(timesteps)) .AND. (t>timesteps(i))) THEN
                        res=interpolant(:,i)
                    END IF
                END DO
            END IF
        
        END SUBROUTINE linear_interpol_1D
        
        SUBROUTINE linear_interpol_2D(t,timesteps,interpolant,res)
            !interpols linearily within a interpolant matrix
            
            USE ABQinterface
            
            IMPLICIT NONE
            
            REAL(KIND=AbqRK),DIMENSION(:),INTENT(IN):: timesteps
            REAL(KIND=AbqRK),DIMENSION(:,:,:),INTENT(IN):: interpolant
            REAL(KIND=AbqRK),DIMENSION(:,:),INTENT(OUT):: res
            REAL(KIND=AbqRK),INTENT(IN):: t
            INTEGER:: i
            
            IF (t<timesteps(1)) THEN
                res=t/timesteps(1)*interpolant(:,:,1)
            ELSE IF (SIZE(timesteps)==1) THEN
                res=interpolant(:,:,1)
            ELSE
                DO i=2,size(timesteps)
                    IF (t<=timesteps(i)) THEN
                        res=interpolant(:,:,i-1)+(interpolant(:,:,i)-interpolant(:,:,i-1))*(t-timesteps(i-1))/(timesteps(i)-timesteps(i-1))
                        EXIT
                    ELSE IF ((i==size(timesteps)) .AND. (t>timesteps(i))) THEN
                        res=interpolant(:,:,i)
                    END IF
                END DO
            END IF
        
        END SUBROUTINE linear_interpol_2D

END MODULE linear_interpolation

MODULE DSTRAN_DFGRD1_exponential_map
!computing the current deformation gradient (DFGRD1) form the rate of deformation
!increment (DSTRAN) and the deformation gradient of the last step (DFGRD0)

USE ABQinterface

CONTAINS

    SUBROUTINE DSTRAN_DFGRD1_map(DSTRAN,DFGRD0,DFGRD1,NDI,NSHR)
    
    INTEGER,DIMENSION(2,3):: Notation_shear=reshape([1,2,1,3,2,3],[2,3])
    INTEGER,INTENT(IN):: NDI,NSHR
    REAL(KIND=AbqRK),DIMENSION(3,3),INTENT(IN):: DFGRD0
    REAL(KIND=AbqRK),DIMENSION(3,3),INTENT(OUT):: DFGRD1
    REAL(KIND=AbqRK),DIMENSION(NDI+NSHR),INTENT(IN):: DSTRAN
    REAL(KIND=AbqRK),DIMENSION(:,:),ALLOCATABLE:: temp_matrix1,temp_matrix2
    REAL(KIND=AbqRK),DIMENSION(:),ALLOCATABLE:: eigen_values
    REAL(KIND=AbqRK),DIMENSION(48)::work
    INTEGER:: i,l,k

    ALLOCATE(temp_matrix1(NDI,NDI),temp_matrix2(NDI,NDI),eigen_values(NDI))
    !build up DSTRAN in full notation (NDI,NDI)
    temp_matrix1=0.0_AbqRK
    FORALL(i=1:NDI) temp_matrix1(i,i)=DSTRAN(i)
    DO i=1,NSHR
        k=Notation_shear(1,i) ; l=Notation_shear(2,i)
        temp_matrix1(k,l)=0.5_AbqRK*DSTRAN(NDI+i) ; temp_matrix1(l,k)=0.5_AbqRK*DSTRAN(NDI+i)
    END DO
    !compute eigenvectors and -values of DSTRAN
    CALL dsyev('V','U',NDI,temp_matrix1,NDI,eigen_values,work,size(work),info)
    !compute exp(DSTRAN) and store it in temp_matrix2
    temp_matrix2=0.0_AbqRK
    DO i=1,NDI
        temp_matrix2=temp_matrix2+EXP(eigen_values(i))*MATMUL(RESHAPE(temp_matrix1(:,i),[NDI,1]),RESHAPE(temp_matrix1(:,i),[1,NDI]))
    END DO
    !finally compute DFGRD1
    DFGRD1(1:NDI,1:NDI)=MATMUL(temp_matrix2,DFGRD0(1:NDI,1:NDI))
    DEALLOCATE(temp_matrix1,temp_matrix2,eigen_values)
    
    END SUBROUTINE DSTRAN_DFGRD1_map
    
END MODULE DSTRAN_DFGRD1_exponential_map

MODULE SMAArrays
    !Module to store the Pointer Adresses globally

    INTEGER(KIND=8),DIMENSION(-1000000:1000000):: PointerAdressesInt
    INTEGER(KIND=8),DIMENSION(-1000000:1000000):: PointerAdressesReal
    
END MODULE SMAArrays
