!DEC$ FREEFORM
!=============================================================================
! 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 FE2MODULE

USE ABQINTERFACE
USE ABQINTERFACE
USE class_Solver
USE type_meshparameters
USE type_macro_GP_DATA
USE type_analysisparameters
USE type_systemmatrices
USE mkl_spblas

CONTAINS

    SUBROUTINE main_program_staggered(PK1,DDSDDE,Grad_U_k_1,GP_DATA,PNEWDT,&
                                      DTIME,TIME,para,nl_geom,&
                                      KSTEP,KINC,analysis,STRESS33)
                            
    !This subroutine is the actual main program for the staggered algorithm.
    !It contains the NR loop and calls the assembling routine, which loops
    !over all elements to compute the overall stiffness, right hand side and
    !macro stress. The displacements are updated until the convergence criterion
    !is fullfilled. At the end of the NR loop the static condensation routine
    !is beeing called and computes the macro stiffness DDSDDE.
        
    IMPLICIT NONE
    
    !this type bundles all macro GP DATA (UGLOBAL,STATEV,...)
    TYPE(macro_GP_DATA),INTENT(INOUT)::GP_DATA
    !this type bundles all mesh parameters
    TYPE(meshparameters),INTENT(IN)::para
    !this type bundles all convergence/analysis parameters
    TYPE (analysisparameters)::analysis
    !macro stress (1. Piola Kirchhoff)        
    REAL(KIND=AbqRK),DIMENSION(:),INTENT(INOUT):: PK1
    !macro stiffness
    REAL(KIND=AbqRK),DIMENSION(:,:),INTENT(OUT):: DDSDDE
    !displacement gradient of this increment
    REAL(KIND=AbqRK),DIMENSION(:),INTENT(IN)::Grad_U_k_1
    !this type bundles all FEM matrices
    TYPE(systemmatrices),TARGET:: matrices
    !If in the staggered scheme the increment fails, this variable is to be set
    !smaller then 1
    REAL(KIND=AbqRK),INTENT(OUT):: PNEWDT
    !loop indices
    INTEGER:: p,iters,i
    !solver type declared in module Solver
    TYPE(solve_soe):: solver
    !error code of solver calls
    INTEGER:: error
    !==0 small deformation theory ==1 large deformation theory
    INTEGER:: nl_geom
    !increment of time
    REAL(KIND=AbqRK),INTENT(IN):: DTIME
    !TIME(1)=current step time, TIME(2)=current total time
    REAL(KIND=AbqRK),DIMENSION(2),INTENT(IN):: TIME
    !number of the step resp. the increment
    INTEGER,INTENT(IN):: KSTEP,KINC
    !maximum force, needed for convergence criterion
    REAL(KIND=AbqRK):: max_force_i
    !maximum norm of residual
    REAL(KIND=AbqRK):: norm_rhs
    !maxval(abs(DELTAU)) of whole time step (stp) or current increment number (i)
    REAL(KIND=AbqRK):: DELTAU_max_stp,DELTAU_max_i
    !stress in thickness direction (plane strain)
    REAL(KIND=AbqRK):: STRESS33

    !------------------allocate the systemmatrices object----------------------

    CALL matrices%allocation(analysis,para,GP_DATA)

    !-----------initialize displacements when in ROM mode----------------------
    
    IF (analysis%solving_process==1) CALL initialize_displacements(para,GP_DATA,matrices)

    !----------------------------initialize solver-----------------------------

    CALL solver%initialize(para,analysis,GP_DATA)

    !---------------------------enforce contraint------------------------------

    CALL enforce_constraint(para,Grad_U_k_1,matrices)

    !---------------------------------NR loop----------------------------------

    iters=0_AbqIK
    DELTAU_max_stp=0.0_AbqRK
    
    DO WHILE (iters<=analysis%max_iters) !only iterate until max_iters is reached

      !-------------------------assemble returns:------------------------------
      !-stiffnessmatrix,righthandside,stress,d_PK1_d_u,d_PK1_d_F,d_r_d_F-

      CALL assemble(para,matrices,PK1,GP_DATA,PNEWDT,nl_geom,&
                    DTIME,TIME,KSTEP,KINC,max_force_i,analysis,STRESS33)

      IF (PNEWDT<1.0_AbqRK) EXIT

      !-------------------factor the matrix by calling the solver---------------

      CALL solver%factor(para,analysis,matrices%k_matrix_values,matrices%k_matrix)

      !--------------check if convergence criterion is reached-----------------

      !get maximum norm of the right hand side
      norm_rhs=MAXVAL(ABS(matrices%rhs_total))
      
      IF (convergence_check(para,analysis,norm_rhs,max_force_i,DELTAU_max_stp,&
                            DELTAU_max_i,iters)) THEN

            IF (iters==analysis%max_iters) THEN
                PNEWDT=0.5_AbqRK
                EXIT
            END IF

       !------------------compute the solution  Kt*deltau=rhs------------------

            CALL solver%solve(para,analysis,matrices%k_matrix_values,&
                              matrices%rhs_total,matrices%DELTAU)

       !get the maximum pos. value of DELTAU of all increments of the time step
            DELTAU_max_i=MAXVAL(ABS(matrices%DELTAU))
            IF (DELTAU_max_i>DELTAU_max_stp) DELTAU_max_stp=DELTAU_max_i
            
       !---------------update UGlOBAL: UGlOBAL=UGlOBAL+DELTAU------------------

            CALL update_displacements(para,GP_DATA,matrices,analysis,Grad_U_k_1)

      ELSE
      
        IF (iters<=analysis%max_iters*0.8) THEN
        
            PNEWDT=1.5_AbqRK
            
        ELSE
      
            PNEWDT=1.0_AbqRK !set PNEWDT to 1.0, because convergence is reached
            
        END IF
                
        EXIT !exit DO LOOP if convergence criterion is reached
                
      END IF

      iters=iters+1
                      
    END DO
    
    !------------call static condensation -> get stiffnessmatrix---------------
    
    IF (PNEWDT>=1.0_AbqRK) CALL static_condensation(para,DDSDDE,PK1,matrices,solver,analysis,STRESS33)

    !---------------------deallocate solver storage----------------------------
                
    CALL solver%finish(para,analysis,GP_DATA)

    !-----------------deallocate the systemmatrices object---------------------

    CALL matrices%deallocation(analysis)
        
    END SUBROUTINE main_program_staggered
    
!==============================================================================
        
    SUBROUTINE main_program_monolithic(PK1,DDSDDE,Grad_U_k_1,GP_DATA,&
                                       DTIME,TIME,para,PNEWDT,nl_geom,&
                                       KSTEP,KINC,analysis,STRESS33)
                            
    !This subroutine enforces the monolithic algorithm. The system of equations
    !for updating the displacements is only called once; the assembling routine
    !is only called once, (or twice at the beginning of the step, since the
    !displacements are extrapolated->the stiffness matricies might have been
    !changed). The parameter analysis%save_soe controls, if the factorized matrix
    !is stored (i.e. allocated space not freed) or not. At the end the static
    !condensation routine computes the macro stiffness DDSDDE.
        
    IMPLICIT NONE
    
    !this type bundles all macro GP DATA (UGLOBAL,STATEV,...)
    TYPE(macro_GP_DATA),INTENT(INOUT)::GP_DATA
    !this type bundles all mesh parameters
    TYPE(meshparameters),INTENT(IN)::para
    !this type bundles all convergence/analysis parameters
    TYPE (analysisparameters)::analysis
    !macro stress (1. Piola Kirchhoff)        
    REAL(KIND=AbqRK),DIMENSION(:),INTENT(INOUT):: PK1
    !macro stiffness
    REAL(KIND=AbqRK),DIMENSION(:,:),INTENT(OUT):: DDSDDE
    !displacement gradient of this increment
    REAL(KIND=AbqRK),DIMENSION(:),INTENT(IN)::Grad_U_k_1
    !this type bundles all FEM matrices
    TYPE(systemmatrices):: matrices
    !loop indices
    INTEGER:: p,iters,i,k
    !solver type declared in module Solver
    TYPE(solve_soe):: solver
    !error code of solver calls
    INTEGER:: error
    !PNEWDT<1.0 smaller timestep required, PNEWDT>1.0 timestep can be increased
    REAL(KIND=AbqRK),INTENT(OUT):: PNEWDT
    !==0 small deformation theory ==1 large deformation theory
    INTEGER:: nl_geom
    !increment of time
    REAL(KIND=AbqRK),INTENT(IN):: DTIME
    !TIME(1)=current step time, TIME(2)=current total time
    REAL(KIND=AbqRK),DIMENSION(2),INTENT(IN):: TIME
    !number of the step resp. the increment
    INTEGER,INTENT(IN):: KSTEP,KINC
    !maximum force needed for convergence criterion
    REAL(KIND=AbqRK):: max_force_i
    !maxval(abs(DELTAU)) of current increment number (i)
    REAL(KIND=AbqRK):: DELTAU_max_i
    !stress in thickness direction (plane strain)
    REAL(KIND=AbqRK):: STRESS33

    !------------------allocate the systemmatrices object----------------------
    
    CALL matrices%allocation(analysis,para,GP_DATA)
    
    !-----------initialize displacements when in ROM mode----------------------
    
    IF (analysis%solving_process==1) CALL initialize_displacements(para,GP_DATA,matrices)

    !---------------------------initialize solver------------------------------

    CALL solver%initialize(para,analysis,GP_DATA)
            
    !--------> beginning of the step call assemble first (necessary------------
    !----------because of extrapolation of the nodal displacements)------------

    IF (GP_DATA%step_indicator(1)<2) THEN
    
      !-------------------------assemble returns:------------------------------
      !----stiffnessmatrix,righthandside,stress,d_PK1_d_u,d_PK1_d_F,d_r_d_F---

      CALL assemble(para,matrices,PK1,GP_DATA,PNEWDT,nl_geom,&
                    DTIME,TIME,KSTEP,KINC,max_force_i,analysis,STRESS33)
                                        
      IF (PNEWDT<1.0_AbqRK) GO TO 20
                    
       !--------------------factor the stiffnessmatrix-------------------------

       CALL solver%factor(para,analysis,matrices%k_matrix_values,matrices%k_matrix)

    ELSE IF ((.NOT. analysis%save_soe) .AND. (analysis%solving_process==0)) THEN
    
    !--middle of timestep, factor stiffnessmatrix if factorization isnt saved--
        
        CALL solver%factor(para,analysis,GP_DATA%k_matrix_values_k)
        
    END IF

    
    !detect linear step, if convergence already reached (through extrapolation
    !from last increment) when a time increment starts

    IF (((GP_DATA%step_indicator(1)==1) .AND. &
    ((MAXVAL(ABS(Grad_U_k_1-GP_DATA%Grad_U_k))>0.005*MAXVAL(ABS(Grad_U_k_1))) &
     .OR. (MAXVAL(ABS(matrices%rhs_total))>0.005*max_force_i))) .OR. &
       (GP_DATA%step_indicator(1)/=1)) THEN
    
    !----------------compute the solution  Kt*deltau=rhs-----------------------
    !------------------rhs has additional contribution!------------------------

      IF (GP_DATA%step_indicator(1)>1) THEN

        matrices%rhs_total(:,1)=GP_DATA%rhs_total_k-MATMUL(GP_DATA%d_r_d_F_k,(Grad_U_k_1-GP_DATA%Grad_U_k))

        CALL solver%solve(para,analysis,GP_DATA%k_matrix_values_k,matrices%rhs_total,matrices%DELTAU)

      ELSE

        matrices%rhs_total(:,1)=matrices%rhs_total(:,1)-MATMUL(matrices%d_r_d_F,(Grad_U_k_1-GP_DATA%Grad_U_k))

        CALL solver%solve(para,analysis,matrices%k_matrix_values,matrices%rhs_total,matrices%DELTAU)

      END IF
                
      !-----update UGlOBAL: UGlOBAL=UGlOBAL+DELTAU and enforce contraint-------

      CALL update_displacements(para,GP_DATA,matrices,analysis,Grad_U_k_1)
      
     !-------------------------assemble returns:------------------------------
     !-stiffnessmatrix,righthandside,stress,d_PK1_d_u,d_PK1_d_F,d_r_d_F-

     CALL assemble(para,matrices,PK1,GP_DATA,PNEWDT,nl_geom,&
                   DTIME,TIME,KSTEP,KINC,max_force_i,analysis,STRESS33)
                                      
     IF (PNEWDT<1.0_AbqRK) GO TO 20
     
     !-------------------factor the matrix by calling the solver----------------

     CALL solver%factor(para,analysis,matrices%k_matrix_values,matrices%k_matrix)
        
    END IF

    !----------save additional information for next iteration step-------------

    GP_DATA%rhs_total_k=matrices%rhs_total(:,1)
    GP_DATA%d_r_d_F_k=matrices%d_r_d_F
    IF (analysis%solving_process==0) GP_DATA%k_matrix_values_k=matrices%k_matrix_values

    !------------call static condensation -> get stiffnessmatrix---------------

    CALL static_condensation(para,DDSDDE,PK1,matrices,solver,analysis,STRESS33)
    
    20  CONTINUE !continue from this point if PNEWDT<1.0_AbqRK
    
    !---------------------deallocate solver storage----------------------------
    
    IF (.NOT. analysis%save_soe) CALL solver%finish(para,analysis,GP_DATA)
     
    !------------------allocate the systemmatrices object----------------------

    CALL matrices%deallocation(analysis)

        
    END SUBROUTINE main_program_monolithic

!==============================================================================
    
    FUNCTION convergence_check(para,analysis,residual,max_force_i,DELTAU_max_stp,&
                               DELTAU_max_i,iteration)
    !this function returns a bool, if convergence is reached or not
    
    IMPLICIT NONE
    
    !.TRUE. if convergence is not reached
    LOGICAL:: convergence_check
    !this type bundles all mesh,convergence & monolithic/staggered parameters
    TYPE(meshparameters),INTENT(IN)::para
    !this type bundles all convergence/analysis parameters
    TYPE (analysisparameters)::analysis
    !maximum residual force current increment
    REAL(KIND=AbqRK),INTENT(IN):: residual
    !maximum force in increment i
    REAL(KIND=AbqRK),INTENT(IN):: max_force_i
    !maxval(abs(DELTAU)) of whole time step (stp) or current increment number (i)
    REAL(KIND=AbqRK),INTENT(IN):: DELTAU_max_stp,DELTAU_max_i
    !number of the current iteration
    INTEGER(KIND=AbqIK),INTENT(IN):: iteration
    
    
    IF ((residual>=analysis%convergence_ratio*max_force_i) &
         .AND. (.NOT. max_force_i==0.0_AbqRK)) THEN
            
       convergence_check=.TRUE.

    ELSE
    
        convergence_check=.FALSE.
        
    END IF
        
    END FUNCTION convergence_check

!==============================================================================
    
    FUNCTION divergence_check(residual_max_3,iteration)
    !this function returns a bool, if the time step appears to be diverging
    !and sets PNEWDT; currently it isn't used but might become useful in the
    !future
    
    IMPLICIT NONE
    
    !.TRUE. if solution diverges
    LOGICAL:: divergence_check
    !number of the current iteration
    INTEGER(KIND=AbqIK),INTENT(IN):: iteration
    !maximum residual force of i-2,i-1,i (i=current increment)
    REAL(KIND=AbqRK),DIMENSION(3),INTENT(IN):: residual_max_3
    
    divergence_check=.FALSE.
    
    IF (iteration>3) THEN
        IF (MINVAL(residual_max_3(2:3))>residual_max_3(1)) THEN
                    
            divergence_check=.TRUE.
                     
        END IF     
    END IF
    
    END FUNCTION divergence_check

!==============================================================================
     
    SUBROUTINE enforce_constraint(para,Grad_U,matrices)
    !this subroutine enforces the constraine:
    !u_slave_i = u_master_i + Eij * deltax_j
    
    IMPLICIT NONE
    
    !this type bundles all FEM matrices
    TYPE(systemmatrices):: matrices
    !this type bundles all mesh parameters
    TYPE(meshparameters),INTENT(IN)::para
    !macro displacement gradient
    REAL(KIND=AbqRK),DIMENSION(:),INTENT(IN)::Grad_U
    !running indice
    INTEGER::p
    
    DO p=1,para%nEquations
        !u_slave_i_ = u_master_i + Eij * deltax_j
        matrices%UGLOBAL_t_1(int(para%Equations(1,p)))=matrices%UGLOBAL_t_1(int(para%Equations(2,p)))+DOT_PRODUCT(para%Equations(3:,p),Grad_U)
    END DO
    
    END SUBROUTINE enforce_constraint

!==============================================================================

    SUBROUTINE update_displacements(para,GP_DATA,matrices,analysis,Grad_U)
    !this subroutine updates the displacements: UGLOBAL_k_1=UGLOBAL_k+DELTAU
    
    IMPLICIT NONE
    
    !this type bundles all macro GP DATA (UGLOBAL,STATEV,...)
    TYPE(macro_GP_DATA)::GP_DATA
    !this type bundles all mesh parameters
    TYPE(meshparameters),INTENT(IN)::para
    !this type bundles all FEM matrices
    TYPE(systemmatrices):: matrices
    !this type bundles all convergence/analysis parameters
    TYPE (analysisparameters)::analysis
    !macro displacement gradient
    REAL(KIND=AbqRK),DIMENSION(:),INTENT(IN)::Grad_U
    INTEGER::p
    
    IF (analysis%solving_process==1) THEN
        !in ROM case, GP_DATA%UGLOBAL contains the amplitudes of the modes
        GP_DATA%UGLOBAL_t_1=GP_DATA%UGLOBAL_t_1+matrices%DELTAU(:,1)
        !compute the actual displacement increment of all nodes
        matrices%DELTAU_FULL=MATMUL(para%ROM_modes,matrices%DELTAU)
    END IF
    
    DO p=1,para%nRows_reduced  !update UGlOBAL UGlOBAL=UGlOBAL+DELTAU
                               !for all nodes except slaves
        matrices%UGLOBAL_t_1(para%master_reduced_to_global(p))=&
        matrices%UGLOBAL_t_1(para%master_reduced_to_global(p))+matrices%DELTAU_FULL(p,1)
    
    END DO
    
    !update the displacements by enforcing the constraint
    CALL enforce_constraint(para,Grad_U,matrices)
    
    END SUBROUTINE update_displacements
    
!==============================================================================

   SUBROUTINE initialize_displacements(para,GP_DATA,matrices)
   !in case of ROM simulation initialize full displacement
   
   IMPLICIT NONE
   
   !this type bundles all macro GP DATA (UGLOBAL,STATEV,...)
    TYPE(macro_GP_DATA)::GP_DATA
    !this type bundles all mesh parameters
    TYPE(meshparameters),INTENT(IN)::para
    !this type bundles all FEM matrices
    TYPE(systemmatrices):: matrices
    REAL(KIND=AbqRk),DIMENSION(para%nRows_reduced):: U_red
    REAL(KIND=AbqRk),DIMENSION(:),POINTER:: U_full
    REAL(KIND=AbqRk),DIMENSION(:),POINTER:: U
    REAL(KIND=AbqRk),DIMENSION(:),POINTER:: Grad_U
    INTEGER::i,p
    
    DO i=1,2 !do it for last time step and current time step
        IF (i==1) THEN
            U=>GP_DATA%UGLOBAL_t_1
            U_full=>matrices%UGLOBAL_t_1
            Grad_U=>GP_DATA%Grad_U_k
        ELSE
            U=>GP_DATA%UGLOBAL_t
            U_full=>matrices%UGLOBAL_t
            Grad_U=>GP_DATA%Grad_U_t
        END IF
        U_red=matmul(para%ROM_modes,U)
        U_full=0.0_AbqRK
        DO p=1,para%nRows_reduced  !initialize masters
            U_full(para%master_reduced_to_global(p))=U_red(p)
        END DO
        DO p=1,para%nEquations !initialize slaves
            !u_slave_i_* = u_master_i + Eij * deltax_j
            U_full(int(para%Equations(1,p)))=U_full(int(para%Equations(2,p)))+DOT_PRODUCT(para%Equations(3:,p),Grad_U)
        END DO
    END DO
    
   END SUBROUTINE initialize_displacements
    
!==============================================================================
    
    SUBROUTINE assemble(para,matrices,PK1,GP_DATA,PNEWDT,nl_geom,&
                        DTIME,TIME,KSTEP,KINC,max_force_i,analysis,STRESS33)
    
    !this subroutine assembles k_matrix_values,rhs_total,STRESS,
    !d_PK1_d_u,d_PK1_d_F,d_r_d_F,all_STATEV_t_1 by looping over all
    !elements and calling the UEL (which itself calls the UMAT at each GP)
    
    USE ExtendedUELModule !needed for an optional UEL argument when calling UEL
                          !in ROM mode with reduced integration
    
    IMPLICIT NONE
    
    !this type bundles all mesh,convergence & monolithic/staggered parameters
    TYPE(meshparameters),INTENT(IN)::para
    !this type bundles all macro GP DATA (UGLOBAL,STATEV,...)
    TYPE(macro_GP_DATA),INTENT(INOUT)::GP_DATA
    !this type bundles all convergence/analysis parameters
    TYPE (analysisparameters)::analysis
    !this type bundles all FEM matrices
    TYPE(systemmatrices):: matrices
    
    !----see in ABAQUS user manual for definition of UEL variabels------
    !right hand side
    REAL(KIND=AbqRK), DIMENSION(para%ndof_e,1):: RHS
    !increment of the displacements
    REAL(KIND=AbqRK), DIMENSION(para%ndof_e,1):: DU
    !element stiffness matrix
    REAL(KIND=AbqRK), DIMENSION(para%ndof_e,para%ndof_e):: AMATRX
    !local/element coordinates
    REAL(KIND=AbqRK),DIMENSION(para%dimens,para%nnodes_e):: COORDS
    !displacement/Velocity/Acceleration
    REAL(KIND=AbqRK), DIMENSION(para%ndof_e):: U,V,A
    REAL(KIND=AbqRK):: PARAMS(0),ADLMAG(0,0),SVARS(para%n_STATEV_elem),&
                       ENERGY(8),PREDEF(2,0,3),DDLMAG(0,0),PERIOD,PNEWDT_in
    INTEGER(KIND=AbqIK):: NPREDF,LFLAGS(5),MDLOAD,NDLOAD,NJPROP,MCRD,NPROPS
    INTEGER(KIND=AbqIK):: JPROPS(1), JDLTYP(0,0)
    INTEGER(KIND=AbqIK),INTENT(IN):: KSTEP,KINC
    REAL(KIND=AbqRK),INTENT(OUT):: PNEWDT
    REAL(KIND=AbqRK),INTENT(IN):: DTIME
    REAL(KIND=AbqRK), DIMENSION(2),INTENT(IN):: TIME
    !--------------------------------------------------------------------
    !temporary matrices, so they doen't have to be allocated every loop in
    !hyperintegration mode
    REAL(KIND=AbqRK),DIMENSION(:,:),ALLOCATABLE:: temp_matrix1,temp_matrix2
    
    !-------------------------------------------------------------------
    !running indices
    INTEGER:: l,i,j,e,elem
    !macro stress (1. Piola Kirchhoff)
    REAL(KIND=AbqRK), DIMENSION(:),INTENT(OUT):: PK1
    !==0 small deformation theory ==1 large deformation theory
    INTEGER:: nl_geom
    !maximum force needed for convergence criterion
    REAL(KIND=AbqRK):: max_force_i
    !stress in thickness direction (plane strain)
    REAL(KIND=AbqRK):: STRESS33
    
    !-------------------------------------------------------------------
    !optional argument for UEL, which cannot be used in UEL when being called from Abaqus
    TYPE(UELextension):: extension
    
    JPROPS(1)=1 !mark the extension to be present
    
    extension%B_Matrices=>para%B_Matrices
    extension%IP_Weights=>para%IP_Weights
    extension%integration_points=>para%integration_points
    
    IF (analysis%training==1) THEN
        extension%training=.TRUE.
        extension%F_int_ip=>GP_DATA%F_int_ip
    ELSE
        extension%training=.FALSE.
    END IF
    
    !---------------------initialize these variables---------------------------

    PNEWDT=2.0_AbqRK
    LFLAGS=[1_AbqIK,nl_geom,1_AbqIK,0_AbqIK,0_AbqIK]
    MCRD=para%dimens
    NJPROP=SIZE(JPROPS)

    !----set all arrays and matrices to zero before building them up again-----
    
    IF (analysis%hyperintegration) THEN
        matrices%k_matrix=0.0_AbqRK
        matrices%rhs_total=0.0_AbqRK
        matrices%d_PK1_d_u=0.0_AbqRK
        IF (.NOT. analysis%symmetric_matrix) matrices%d_r_d_F=0.0_AbqRK
        matrices%d_PK1_d_F=0.0_AbqRK
        ALLOCATE(temp_matrix1(para%ndof_e,para%n_ROM_modes),temp_matrix2(para%ndof_e,para%ndof_macro))
    ELSE
        matrices%k_matrix_values=0.0_AbqRK
        matrices%rhs_total_full=0.0_AbqRK
        matrices%d_PK1_d_u_full=0.0_AbqRK
        IF (.NOT. analysis%symmetric_matrix) matrices%d_r_d_F_full=0.0_AbqRK
        matrices%d_PK1_d_F=0.0_AbqRK
    END IF
    PK1=0.0_AbqRK
    max_force_i=0.0_AbqRK
    extension%STRESS33=0.0_AbqRK
    
    !-------loop over all elements to build k_matrix_values,rhs_total,---------
    !-----------STRESS,d_PK1_d_u,d_PK1_d_F,all_STATEV_t_1----------------
    
    DO elem=1,para%n_elem_to_hyper !loop over all elements containing active integration points
    
        e=para%elem_to_hyper(elem)
        
        !------------read element coordinates form global COORDSGLOBAL-----------
        
        l=1
        DO i=1,para%nnodes_e
            DO j=1,para%dimens
            COORDS(j,i)=para%COORDSGLOBAL(para%element_to_global_dof(l,e))
            l=l+1
            END DO
        END DO
    
        !---------read the element displacement vector U and increment DU--------
        
        DO i=1,para%ndof_e
            DU(i,1)=matrices%UGLOBAL_t_1(para%element_to_global_dof(i,e))&
                    -matrices%UGLOBAL_t(para%element_to_global_dof(i,e))
            U(i)=matrices%UGLOBAL_t_1(para%element_to_global_dof(i,e))
        END DO
    
        !------------------call the element routine UEL--------------------------
        
        SVARS=GP_DATA%STATEV_t(:,elem) !read internal state variables
        RHS=0.0_AbqRK
        AMATRX=0.0_AbqRK
        
        !number of Properties of the elements material
        NPROPS=int(para%PROPS(para%element_to_material(e),1))
            
        !call user element library with the optional argument extension
        CALL UEL(RHS,AMATRX,SVARS,ENERGY,para%ndof_e,1,para%n_STATEV_elem,&
            para%PROPS(para%element_to_material(e),2:NPROPS+1),NPROPS,&
            COORDS,MCRD,para%nnodes_e,U,DU,V,A,para%JTYPE,TIME,DTIME,KSTEP,&
            KINC,e,PARAMS,NDLOAD,JDLTYP,ADLMAG,PREDEF,NPREDF,LFLAGS,&
            para%ndof_e,DDLMAG,MDLOAD,PNEWDT_in,JPROPS,NJPROP,PERIOD,extension)
        
        IF (PNEWDT_in<PNEWDT) PNEWDT=PNEWDT_in
                        
        GP_DATA%STATEV_t_1(:,elem)=SVARS !save SVARS for the next timestep t+1
        
        
        !-----------sort AMATRX values into global k_matrix_values array---------
        !-------------sort RHS values into global rhs_total array----------------
        
        !take symmetric part of AMATRX, if symmetric_matrix==.TRUE.
        IF (analysis%symmetric_matrix) AMATRX=0.5_AbqRK*(AMATRX+TRANSPOSE(AMATRX))
        
        IF (analysis%hyperintegration) THEN
            CALL assemble_ROM_stiffness_and_RHS(para,analysis,matrices,RHS,AMATRX,PK1,max_force_i,e,elem,temp_matrix1,temp_matrix2)
        ELSE
            CALL assemble_full_stiffness_and_RHS(para,analysis,matrices,RHS,AMATRX,PK1,max_force_i,e)
        END IF
        
    END DO
    
    !when being in ROM mode, but the problem is fully integrated, do the ROM projection with
    !the unreduced Matrices by multiplying with the ROM modes
    IF (analysis%solving_process==1) THEN
        IF (.NOT. analysis%hyperintegration) THEN
            CALL ROM_projection(matrices,para,analysis)
        ELSE
            DEALLOCATE(temp_matrix1,temp_matrix2)
        END IF
    END IF
    
    IF (analysis%symmetric_matrix) THEN !use the fact that the stiffness matrix is symmetric
        matrices%d_r_d_F=TRANSPOSE(matrices%d_PK1_d_u)
    END IF
    
    STRESS33=extension%STRESS33
    
    END SUBROUTINE assemble
    
    !==============================================================================
    
    SUBROUTINE assemble_full_stiffness_and_RHS(para,analysis,matrices,RHS,AMATRX,PK1,max_force_i,e)
    !this subroutine assembles the full stiffness matrices and RHS's
    
    IMPLICIT NONE
    
    !this type bundles all mesh,convergence & monolithic/staggered parameters
    TYPE(meshparameters),INTENT(IN)::para
    !this type bundles all FEM matrices
    TYPE(systemmatrices),INTENT(INOUT):: matrices
    !this type bundles all convergence/analysis parameters
    TYPE (analysisparameters),INTENT(IN)::analysis
    !right hand side
    REAL(KIND=AbqRK), DIMENSION(para%ndof_e,1),INTENT(INOUT):: RHS
    !element stiffness matrix
    REAL(KIND=AbqRK), DIMENSION(para%ndof_e,para%ndof_e):: AMATRX
    !macro stress (1. Piola Kirchhoff)
    REAL(KIND=AbqRK), DIMENSION(:),INTENT(OUT):: PK1
    !maximum force needed for convergence criterion
    REAL(KIND=AbqRK),INTENT(INOUT):: max_force_i
    REAL(KIND=AbqRK):: max_force_i_elem
    INTEGER:: l,i,j,i1,i2,j1,j2
    !element number
    INTEGER,INTENT(IN)::e
    
    l=0
    DO i=1,para%ndof_e
        DO j=1,para%ndof_e
            l=l+1
            IF (para%values_to_global_reduced(l,e)/=-1) THEN!==-1 indicates entry not needed
                matrices%k_matrix_values(para%values_to_global_reduced(l,e))=&
                matrices%k_matrix_values(para%values_to_global_reduced(l,e))+AMATRX(i,j)
            END IF
        END DO
        
        i1=ABS(para%element_to_reduced_dof(i,e))
        IF (i1/=0) THEN    !==0 indicates boundary condition
            matrices%rhs_total_full(i1,1)=matrices%rhs_total_full(i1,1)+RHS(i,1)
        END IF
    END DO
    
    !-sort slave AMATRX and RHS values into d_PK1_d_F,d_PK1_d_u,d_r_d_F and PK1-
    IF (ANY (para%element_to_reduced_dof(:,e)<=0)) THEN !only elements with slaves contribute
    
        IF (.NOT. analysis%symmetric_matrix==1) THEN !avoid inefficient if queries in the loop over i and j
        
            DO i=1,para%ndof_e !loop over all entries of AMATRX
                i1=para%element_to_reduced_dof(i,e)
                i2=para%global_node_to_slave_node(para%element_to_global_dof(i,e))
                IF (i2>0) THEN
                    DO j=1,para%ndof_e
                        j1=para%element_to_reduced_dof(j,e)
                        j2=para%global_node_to_slave_node(para%element_to_global_dof(j,e))
                        
                        IF (j2>0) THEN !!Kii!!
                            matrices%d_PK1_d_F=matrices%d_PK1_d_F+AMATRX(i,j)*MATMUL(RESHAPE(para%Equations(3:,i2),[para%ndof_macro,1]),RESHAPE(para%Equations(3:,j2),[1,para%ndof_macro]))
                            IF (ABS(j1)>0) matrices%d_PK1_d_u_full(:,ABS(j1))=matrices%d_PK1_d_u_full(:,ABS(j1))+AMATRX(i,j)*para%Equations(3:,i2)
                            IF (ABS(i1)>0) matrices%d_r_d_F_full(ABS(i1),:)=matrices%d_r_d_F_full(ABS(i1),:)+AMATRX(i,j)*para%Equations(3:,j2)
                        ELSE IF (j1>0) THEN !!Kid!!
                            matrices%d_PK1_d_u_full(:,j1)=matrices%d_PK1_d_u_full(:,j1)+AMATRX(i,j)*para%Equations(3:,i2)
                        END IF
                    END DO
                        
                    !build macro STRESS (1.PK) from slave node rhs contribution
                    PK1=PK1-RHS(i,1)*para%Equations(3:,i2)
                    
                    !get maxval of RHS
                    max_force_i_elem=ABS(RHS(i,1))
                    IF (max_force_i_elem>max_force_i) max_force_i=max_force_i_elem
                    
                ELSE IF (i1>0) THEN !!Kdi!!
                    DO j=1,para%ndof_e
                        j2=para%global_node_to_slave_node(para%element_to_global_dof(j,e))
                        IF (j2>0) THEN
                            matrices%d_r_d_F_full(i1,:)=matrices%d_r_d_F_full(i1,:)+AMATRX(i,j)*para%Equations(3:,j2)  !contribution to d_r_d_F
                        END IF
                    END DO
                END IF
            END DO
        
        ELSE
        
            DO i=1,para%ndof_e !loop over all entries of AMATRX
                i1=para%element_to_reduced_dof(i,e)
                i2=para%global_node_to_slave_node(para%element_to_global_dof(i,e))
                IF (i2>0) THEN
                    DO j=i,para%ndof_e
                        j1=para%element_to_reduced_dof(j,e)
                        j2=para%global_node_to_slave_node(para%element_to_global_dof(j,e))
                        IF (j2>0) THEN !!Kii!!
                            matrices%d_PK1_d_F=matrices%d_PK1_d_F+AMATRX(i,j)*MATMUL(RESHAPE(para%Equations(3:,i2),[para%ndof_macro,1]),RESHAPE(para%Equations(3:,j2),[1,para%ndof_macro]))
                            IF (ABS(j1)>0) matrices%d_PK1_d_u_full(:,ABS(j1))=matrices%d_PK1_d_u_full(:,ABS(j1))+AMATRX(i,j)*para%Equations(3:,i2)
                        ELSE IF (j1>0) THEN !!Kid!!
                            matrices%d_PK1_d_u_full(:,j1)=matrices%d_PK1_d_u_full(:,j1)+AMATRX(i,j)*para%Equations(3:,i2)
                        END IF
                    END DO
                        
                    !build macro STRESS (1.PK) from slave node rhs contribution
                    PK1=PK1-RHS(i,1)*para%Equations(3:,i2)
                
                    !get maxval of RHS
                    max_force_i_elem=ABS(RHS(i,1))
                    IF (max_force_i_elem>max_force_i) max_force_i=max_force_i_elem
                
                END IF
            END DO
        END IF
    END IF
    
    END SUBROUTINE assemble_full_stiffness_and_RHS
    
    !==============================================================================
    
    SUBROUTINE assemble_ROM_stiffness_and_RHS(para,analysis,matrices,RHS,AMATRX,PK1,max_force_i,e,elem,temp_matrix1,temp_matrix2)
    !this subroutine assembles the ROM stiffness matrices and RHS's using the ROM Galerkin projection modes, by directly sorting
    !the AMATRX entries to the respectiv ROM entries
    
    IMPLICIT NONE
    
    !this type bundles all mesh,convergence & monolithic/staggered parameters
    TYPE(meshparameters),INTENT(IN)::para
    !this type bundles all FEM matrices
    TYPE(systemmatrices),INTENT(INOUT):: matrices
    !this type bundles all convergence/analysis parameters
    TYPE (analysisparameters),INTENT(IN)::analysis
    !right hand side
    REAL(KIND=AbqRK), DIMENSION(:,:),INTENT(IN):: RHS
    !element stiffness matrix
    REAL(KIND=AbqRK), DIMENSION(:,:),INTENT(IN):: AMATRX
    REAL(KIND=AbqRK), DIMENSION(:,:),ALLOCATABLE:: C
    !macro stress (1. Piola Kirchhoff)
    REAL(KIND=AbqRK), DIMENSION(:),INTENT(OUT):: PK1
    REAL(KIND=AbqRK), DIMENSION(:,:),INTENT(INOUT):: temp_matrix1,temp_matrix2
    !maximum force needed for convergence criterion
    REAL(KIND=AbqRK),INTENT(INOUT):: max_force_i
    REAL(KIND=AbqRK):: max_force_i_elem
    !element number
    INTEGER,INTENT(IN)::e,elem
    
    !build up K_dd (ROM reduced)
    !-> at first compute K_elem*Phi_elem and store it in a temporary matrix since it is needed again
    IF (analysis%symmetric_matrix) THEN !symmetric matrix multiplication
        CALL dsymm('L','U',para%ndof_e,para%n_ROM_modes,1.0_AbqRK,AMATRX,para%ndof_e,para%ROM_modes_elem(:,:,elem),para%ndof_e,0.0_AbqRK,temp_matrix1,para%ndof_e)
    ELSE
        temp_matrix1=MATMUL(AMATRX,para%ROM_modes_elem(:,:,elem))
    END IF
    
    !-> now compute Phi_elem^T*(K_elem*Phi_elem)
    matrices%k_matrix=matrices%k_matrix+MATMUL(TRANSPOSE(para%ROM_modes_elem(:,:,elem)),temp_matrix1)
    
    !build up rhs (ROM reduced)
    matrices%rhs_total=matrices%rhs_total+MATMUL(TRANSPOSE(para%ROM_modes_elem(:,:,elem)),RHS)
    
    IF (ANY (para%element_to_reduced_dof(:,e)<=0)) THEN
        
        !-> now compute A_elem^T*(K_elem*Phi_elem) using the temporary matrix to compute d_PK1_d_u
        matrices%d_PK1_d_u=matrices%d_PK1_d_u+MATMUL(TRANSPOSE(para%equations_elem(:,:,elem)),temp_matrix1)
        
        !-> at first compute K_elem*A_elem and store it in a temporary matrix
        IF (analysis%symmetric_matrix) THEN
            CALL dsymm('L','U',para%ndof_e,para%ndof_macro,1.0_AbqRK,AMATRX,para%ndof_e,para%equations_elem(:,:,elem),para%ndof_e,0.0_AbqRK,temp_matrix2,para%ndof_e)
        ELSE
            temp_matrix2=MATMUL(AMATRX,para%equations_elem(:,:,elem))
            !-> now compute A_elem^T*(K_elem*Phi_elem) to compute d_PK1_d_F when being in unsymmetric mode
            matrices%d_r_d_F=matrices%d_r_d_F+MATMUL(TRANSPOSE(para%ROM_modes_elem(:,:,elem)),temp_matrix2)
        END IF
        
        !-> now compute A_elem^T*(K_elem*A_elem) to compute d_r_d_F
        matrices%d_PK1_d_F=matrices%d_PK1_d_F+MATMUL(TRANSPOSE(para%equations_elem(:,:,elem)),temp_matrix2)
        
        !compute macro stress
        PK1=PK1-MATMUL(TRANSPOSE(para%equations_elem(:,:,elem)),RHS(:,1))
        
        !get the maximum RHS value of all elements
        max_force_i_elem=MAXVAL(ABS(RHS))
        IF (max_force_i_elem>max_force_i) max_force_i=max_force_i_elem
        
    END IF
    
    END SUBROUTINE assemble_ROM_stiffness_and_RHS

!==============================================================================

    SUBROUTINE ROM_projection(matrices,para,analysis)
    !this routine performs the rom projection phi^T*f_int=f_int_reduced,
    !phi^T*K*phi=K_reduced etc. which is more efficient in case of ROM without
    !hyperintegration
        
    IMPLICIT NONE
    
    !this type bundles all FEM matrices
    TYPE(systemmatrices),INTENT(INOUT):: matrices
    !this type bundles all mesh,convergence & monolithic/staggered parameters
    TYPE(meshparameters),INTENT(IN)::para
    !this type bundles all convergence/analysis parameters
    TYPE (analysisparameters)::analysis
    TYPE(SPARSE_MATRIX_T)::k_matrix_handle
    TYPE(MATRIX_DESCR)::description
    INTEGER:: error
    REAL(KIND=AbqRK),DIMENSION(para%nRows_reduced,para%n_ROM_modes)::C !temporary matrix
    
    !do the ROM projections
    matrices%rhs_total=MATMUL(TRANSPOSE(para%ROM_modes),matrices%rhs_total_full)
    matrices%d_PK1_d_u=MATMUL(matrices%d_PK1_d_u_full,para%ROM_modes)
    IF (.NOT. analysis%symmetric_matrix==1) matrices%d_r_d_F=MATMUL(TRANSPOSE(para%ROM_modes),matrices%d_r_d_F_full)
    
    !create a handle containing the full stiffnessmatrix and its sparsity structure
    error=mkl_sparse_d_create_csr(k_matrix_handle,SPARSE_INDEX_BASE_ONE,para%nRows_reduced,&
                                   para%nRows_reduced,para%rowIndex_reduced(1:para%nRows_reduced),&
                                   para%rowIndex_reduced(2:para%nRows_reduced+1),para%columnIndex_reduced,&
                                   matrices%k_matrix_values)
    IF (analysis%symmetric_matrix) THEN
        description%type=SPARSE_MATRIX_TYPE_SYMMETRIC
        description%mode=SPARSE_FILL_MODE_UPPER
        description%diag=SPARSE_DIAG_NON_UNIT
    ELSE
        description%type=SPARSE_MATRIX_TYPE_GENERAL
    END IF
    C=0.0_AbqRK
    !perform a sparse multiplication k_matrix_full*phi
    error=mkl_sparse_d_mm(SPARSE_OPERATION_NON_TRANSPOSE,1.0_AbqRK,k_matrix_handle,description,&
                           SPARSE_LAYOUT_COLUMN_MAJOR,para%ROM_modes,para%n_ROM_modes,para%nRows_reduced,&
                           0.0_AbqRK,C,para%nRows_reduced)
    !perform a dense multiplication phi^T*(k_matrix_full*phi)
    matrices%k_matrix=MATMUL(TRANSPOSE(para%ROM_modes),C)
    error=mkl_sparse_destroy(k_matrix_handle) !free the allocated memory for the matrix handle
    
    END SUBROUTINE ROM_projection

!==============================================================================
    
    SUBROUTINE static_condensation(para,DDSDDE,PK1,matrices,solver,analysis,STRESS33)
    !this subroutine performs the static condensation by computing DDSDDE
    !(macro stiffness matrix) [essentially call the solver and compute
    !k_matrix*d_u_d_E=(-1)*d_r_d_F and DDSDDE=d_PK1_d_F+d_PK1_d_u*d_u_d_E
    
    IMPLICIT NONE
    
    !this type bundles all mesh,convergence & monolithic/staggered parameters
    TYPE(meshparameters),INTENT(IN)::para
    !this type bundles all convergence/analysis parameters
    TYPE (analysisparameters)::analysis
    !macro stress (1. Piola Kirchhoff)
    REAL(KIND=AbqRK), DIMENSION(:),INTENT(INOUT):: PK1
    !macro tangentstiffness
    REAL(KIND=AbqRK), DIMENSION(:,:),INTENT(OUT):: DDSDDE
    !this type bundles all FEM matrices
    TYPE(systemmatrices):: matrices
    !partial derivative of micro displacement w.r.t macro deformation gradient
    REAL(KIND=AbqRK), DIMENSION(:,:),ALLOCATABLE:: d_u_d_F
    !partial derivative of macro stress w.r.t micro reaction force
    REAL(KIND=AbqRK), DIMENSION(:,:),ALLOCATABLE::d_PK1_d_r
    !running indice
    INTEGER:: i,j
    TYPE(solve_soe):: solver
    !stress in thickness direction (plane strain)
    REAL(KIND=AbqRK):: STRESS33
        
    !--------------------initialize Arrays to be computed----------------------
    
    !Allocate Matrix of derivatives with respect to the full or reduced (ROM) residual
    IF (analysis%solving_process==1) THEN
        ALLOCATE(d_u_d_F(para%n_ROM_modes,para%ndof_macro))
    ELSE IF (analysis%solving_process==0) THEN
        ALLOCATE(d_u_d_F(para%nRows_reduced,para%ndof_macro))
    END IF
    
    d_u_d_F=0.0_AbqRK

    !----------------------------compute d_u_d_F-------------------------------
    
    CALL solver%solve(para,analysis,matrices%k_matrix_values,matrices%d_r_d_F,d_u_d_F)

    !-------------------------compute DDSDDE-----------------------------------    
        
    DDSDDE=matrices%d_PK1_d_F-MATMUL(matrices%d_PK1_d_u,d_u_d_F)

    !-----------compute algorithmically consistent stress in monolithic case------------      
    
    IF (analysis%monolithic) THEN
        
        !Allocate Matrix of derivative of PK1 with respect to the full or reduced (ROM) residual
        IF (analysis%solving_process==0) THEN
            ALLOCATE(d_PK1_d_r(para%nRows_reduced,para%ndof_macro))
        ELSE IF (analysis%solving_process==1) THEN
            ALLOCATE(d_PK1_d_r(para%n_ROM_modes,para%ndof_macro))
        END IF

        CALL solver%solve(para,analysis,matrices%k_matrix_values,TRANSPOSE(matrices%d_PK1_d_u),d_PK1_d_r)

        !in monolithic case STRESS has additional contributions
        PK1=PK1-MATMUL(TRANSPOSE(d_PK1_d_r),matrices%rhs_total(:,1))
        
        DEALLOCATE(d_PK1_d_r)
        
    END IF

    !----------------divide by the RVEs (Area=2D, Volume=3D)-------------------
    
    DDSDDE=DDSDDE/para%RVE_Volume(1)
    PK1=PK1/para%RVE_Volume(1)
    STRESS33=STRESS33/para%RVE_Volume(1)
        
END SUBROUTINE static_condensation

END MODULE FE2MODULE
