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

USE ABQinterface
USE type_solver
USE type_meshparameters
USE type_macro_GP_DATA
USE type_analysisparameters
USE type_systemmatrices
USE mkl_spblas

CONTAINS

    SUBROUTINE main_program_staggered(macro_response,DDSDDE,macro_measure,&
                                      GP_DATA,PNEWDT,DTIME,TIME,para,nl_geom,&
                                      KSTEP,KINC,analysis)
    
    !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 system matrix of derivatives and 
    !right hand sides. The nodal solution is 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 response and its derivative
    !w.r.t. the macro measures
        
    IMPLICIT NONE
    
    !this type bundles all macro GP DATA (nodal_sol,SVARS,...)
    TYPE(macro_GP_DATA),INTENT(INOUT)::GP_DATA
    !this type bundles all mesh parameters
    TYPE(meshparameters),INTENT(INOUT)::para
    !this type bundles all convergence/analysis parameters
    TYPE (analysisparameters)::analysis
    !macro response (e.g. 1. Piola Kirchhoff)
    REAL(KIND=AbqRK),DIMENSION(:),INTENT(INOUT):: macro_response
    !derivative of macro repsonse w.r.t macro measures
    REAL(KIND=AbqRK),DIMENSION(:,:),INTENT(OUT):: DDSDDE
    !macro measure (e.g. displacement gradient of this increment)
    REAL(KIND=AbqRK),DIMENSION(:),INTENT(IN)::macro_measure
    !this type bundles all FEM matrices
    TYPE(systemmatrices):: matrices
    !If in the staggered scheme no convergence is reache, this variable is to be
    !set smaller than 1
    REAL(KIND=AbqRK),INTENT(OUT):: PNEWDT
    !loop indices
    INTEGER:: p,iters,i
    !solver type declared in the sover module
    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 element reaction, needed for convergence criterion
    REAL(KIND=AbqRK),DIMENSION(para%ndof_n_max):: max_reaction_elem

    !------------------allocate the systemmatrices object----------------------
    
    CALL matrices%allocation(analysis,para,GP_DATA)
    
    !----------------------------initialize solver-----------------------------
    
    CALL solver%initialize(para,analysis,GP_DATA)
    
    !-------put the macroscopic measure (e.g. displacement gradient) into------
    !------------the nodal solution vector & enforce the constraints-----------
    
    IF (analysis%ROM_projection) THEN
        GP_DATA%nodal_sol_t_1(para%n_ROM_modes+1:)=macro_measure
        !---------------initialize solution when in ROM mode-------------------
        CALL initialize_solution_vector(para,GP_DATA,matrices)
    ELSE
        FORALL(p=1:SIZE(macro_measure)) GP_DATA%nodal_sol_t_1(para%macro_measure_to_micro_DOF(p))=macro_measure(p)
        !---------------------------enforce contraint--------------------------
        CALL enforce_constraint(para,GP_DATA%nodal_sol_t_1)
    END IF
    
    !---------------------------------NR loop----------------------------------

    iters=0_AbqIK
    
    DO WHILE (iters<=analysis%max_iters) !only iterate until max_iters is reached
    
      !-----------assemble returns: matrices of derivatives, RHS's-------------
    
      CALL assemble(para,matrices,GP_DATA,analysis,PNEWDT,nl_geom,DTIME,&
                        TIME,KSTEP,KINC,max_reaction_elem)
    
      IF (PNEWDT<1.0_AbqRK) EXIT

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

      CALL solver%factor(para,analysis,matrices%derivatives(1,1)%matrix)
    
      !--------------check if convergence criterion is reached-----------------
    
      IF (convergence_check(para,analysis,matrices,max_reaction_elem)) 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%derivatives(1,1)%matrix,&
                              matrices%rhs(1)%matrix,matrices%nodal_sol_incr)
        
       !------update nodal solultion values: value_k_1=value_k+incr------------
            
            CALL update_nodal_sol(para,GP_DATA,matrices,analysis)
            
      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 DDSDDE-------------------
    
    IF (PNEWDT>=1.0_AbqRK) CALL static_condensation(para,DDSDDE,macro_response,&
                                               matrices,solver,analysis)
    
    !---------------------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(macro_response,DDSDDE,macro_measure,&
                                       GP_DATA,DTIME,TIME,para,PNEWDT,nl_geom,&
                                       KSTEP,KINC,analysis)
                            
    !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
    !nodal solution is extrapolated). At the end the static condensation
    !routine is beeing called and computes the macro response and its derivative
    !w.r.t. the macro measures
        
    IMPLICIT NONE
    
    !this type bundles all macro GP DATA (nodal_sol,SVARS,...)
    TYPE(macro_GP_DATA),INTENT(INOUT)::GP_DATA
    !this type bundles all mesh parameters
    TYPE(meshparameters),INTENT(INOUT)::para
    !this type bundles all convergence/analysis parameters
    TYPE (analysisparameters)::analysis
    !macro response (e.g. 1. Piola Kirchhoff)
    REAL(KIND=AbqRK),DIMENSION(:),INTENT(INOUT):: macro_response
    !derivative of macro repsonse w.r.t macro measures
    REAL(KIND=AbqRK),DIMENSION(:,:),INTENT(OUT):: DDSDDE
    !macro measure (e.g. displacement gradient of this increment)
    REAL(KIND=AbqRK),DIMENSION(:),INTENT(IN):: macro_measure
    !macro measure increment from NR Step k to k+1
    REAL(KIND=AbqRK),DIMENSION(para%n_Reaction_force_dof):: macro_measure_incr
    !this type bundles all FEM matrices
    TYPE(systemmatrices):: matrices
    !loop indices
    INTEGER:: p,iters,i,k
    !solver type declared in the solver module
    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 element reaction, needed for convergence criterion
    REAL(KIND=AbqRK),DIMENSION(para%ndof_n_max):: max_reaction_elem

    !------------------allocate the systemmatrices object----------------------
    
    CALL matrices%allocation(analysis,para,GP_DATA)
    
    !-----------initialize displacements when in ROM mode----------------------
    
    IF (analysis%ROM_projection) CALL initialize_solution_vector(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 solution)---------------

    IF (GP_DATA%step_indicator<2) THEN
    
      !-----------assemble returns: matrices of derivatives, RHS's-------------

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

       CALL solver%factor(para,analysis,matrices%derivatives(1,1)%matrix)

    ELSE IF (.NOT. analysis%save_soe) THEN
    
    !-middle of timestep, factor the system matrix 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) .AND. &
         convergence_check(para,analysis,matrices,max_reaction_elem)) .OR. &
         (GP_DATA%step_indicator/=1)) THEN
    
    !----------compute the solution  Kt*nodal_value_incr=rhs-------------------
    !---------------rhs has additional contribution!---------------------------
      
      !get macro measure increment
      DO p=1,para%n_Reaction_force_dof
          !temporarely store the macro measure of the last NR step
          macro_measure_incr(p)=macro_measure(p)-matrices%nodal_sol_t_1(para%macro_measure_to_micro_DOF(p))
      END DO
      
      IF (GP_DATA%step_indicator>1) THEN

        matrices%rhs(1)%matrix(:,1)=GP_DATA%rhs_k(:,1)-MATMUL(GP_DATA%d_rhs_d_macro_k,macro_measure_incr)
        
        CALL solver%solve(para,analysis,GP_DATA%k_matrix_values_k,matrices%rhs(1)%matrix,matrices%nodal_sol_incr)

      ELSE

        matrices%rhs(1)%matrix(:,1)=matrices%rhs(1)%matrix(:,1)-MATMUL(matrices%derivatives(1,2)%matrix,macro_measure_incr)
        
        CALL solver%solve(para,analysis,matrices%derivatives(1,1)%matrix,matrices%rhs(1)%matrix,matrices%nodal_sol_incr)

      END IF
      
      !-------put the macroscopic measure (e.g. displacement gradient) into------
      !---------------------the nodal solution vector----------------------------
      
      FORALL(p=1:SIZE(macro_measure)) matrices%nodal_sol_t_1(para%macro_measure_to_micro_DOF(p))=macro_measure(p)
      IF (analysis%ROM_projection) GP_DATA%nodal_sol_t_1(para%n_ROM_modes+1:)=macro_measure
      
      !------update nodal solultion values: value_k_1=value_k+incr------------

      CALL update_nodal_sol(para,GP_DATA,matrices,analysis)
      
      !-----------assemble returns: matrices of derivatives, RHS's-------------

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

     CALL solver%factor(para,analysis,matrices%derivatives(1,1)%matrix)
        
    END IF

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

    GP_DATA%rhs_k=matrices%rhs(1)%matrix
    GP_DATA%d_rhs_d_macro_k=matrices%derivatives(1,2)%matrix
    IF (.NOT. analysis%ROM_projection) GP_DATA%k_matrix_values_k=matrices%derivatives(1,1)%matrix

    !-----------------call static condensation -> get DDSDDE-------------------

    CALL static_condensation(para,DDSDDE,macro_response,matrices,solver,analysis)
    
    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,matrices,max_reaction_elem)
    !this function returns a bool, if convergence is reached or not
    !at the moment it is only suitable for mechanical problems!! if the
    !problem is e.g. thermomechanical the residuals should be assessed
    !for mechanical and thermal forces separately!
    
    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
    !this type bundles all FEM matrices
    TYPE(systemmatrices):: matrices
    !maximum force of all elements
    REAL(KIND=AbqRK),INTENT(IN),DIMENSION(para%ndof_n_max):: max_reaction_elem
    
    IF ((MAXVAL(ABS(matrices%rhs(1)%matrix(:,1)))>analysis%convergence_ratio*MAXVAL(max_reaction_elem))&
        .AND. (.NOT. MAXVAL(max_reaction_elem)==0.0_AbqRK)) THEN
        convergence_check=.TRUE.
    ELSE
        convergence_check=.FALSE.
    END IF
    
    END FUNCTION convergence_check

!==============================================================================
     
    SUBROUTINE enforce_constraint(para,nodal_sol)
    !this subroutine enforces the constraints:
    !x_1=b*x_2+c*x_3+... (assuming a=-1)
    
    IMPLICIT NONE
    
    !nodal solution
    REAL(KIND=AbqRK),DIMENSION(:),INTENT(INOUT):: nodal_sol
    !this type bundles all mesh parameters
    TYPE(meshparameters),INTENT(IN)::para
    !running indice
    INTEGER::i,j
    
    DO i=1,para%nEquations
        nodal_sol(para%Equation_members(para%Equation_n_members(i)))=0.0_AbqRK
        DO j=para%Equation_n_members(i)+1,para%Equation_n_members(i+1)-1
            nodal_sol(para%Equation_members(para%Equation_n_members(i)))=&
            nodal_sol(para%Equation_members(para%Equation_n_members(i)))+&
            nodal_sol(para%Equation_members(j))*para%Equation_factors(j)
        END DO
    END DO
    
    END SUBROUTINE enforce_constraint
    
!==============================================================================

   SUBROUTINE initialize_solution_vector(para,GP_DATA,matrices)
   !in case of ROM simulation initialize full solution vector
   
   IMPLICIT NONE
   
   !this type bundles all macro GP DATA (nodal_sol,STATEV,...)
   TYPE(macro_GP_DATA)::GP_DATA
   !this type bundles all FEM matrices
   TYPE(systemmatrices):: matrices
   !this type bundles all mesh parameters
   TYPE(meshparameters),INTENT(IN)::para
   REAL(KIND=AbqRk),DIMENSION(:),POINTER:: nodal_sol,nodal_sol_full
   INTEGER::i,p
   
   DO i=1,2 !do it for last time step and current time step
       IF (i==1) THEN !current timestep
           nodal_sol=>GP_DATA%nodal_sol_t_1
           nodal_sol_full=>matrices%nodal_sol_t_1
       ELSE !last timestep
           nodal_sol=>GP_DATA%nodal_sol_t
           nodal_sol_full=>matrices%nodal_sol_t
       END IF
       !solution of all independent DOFs
       matrices%nodal_sol_incr_full(:,1)=MATMUL(para%ROM_modes,nodal_sol(1:para%n_ROM_modes))
       !insert full solution values
       DO p=1,para%nRows_reduced
           nodal_sol_full(para%reduced_dof_to_global_dof(1,p))=matrices%nodal_sol_incr_full(p,1)
       END DO
       !insert macro values
       DO p=1,para%n_Reaction_force_dof+para%n_additional_dof
           nodal_sol_full(para%macro_measure_to_micro_DOF(p))=nodal_sol(para%n_ROM_modes+p)
       END DO
       !enforce constraint
       CALL enforce_constraint(para,nodal_sol_full)
   END DO
    
   END SUBROUTINE initialize_solution_vector
    
!==============================================================================


    SUBROUTINE update_nodal_sol(para,GP_DATA,matrices,analysis)
    !this subroutine updates the nodal solution: 
    
    IMPLICIT NONE
    
    !this type bundles all macro GP DATA (nodal_sol,SVARS...)
    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
    INTEGER::p
    
    IF (analysis%ROM_projection) THEN
        !in ROM case, GP_DATA%nodal_sol contains the amplitudes of the modes
        GP_DATA%nodal_sol_t_1(1:para%n_ROM_modes)=GP_DATA%nodal_sol_t_1(1:para%n_ROM_modes)&
                                                  +matrices%nodal_sol_incr(:,1)
        !compute the actual nodal value increment from the amplitudal increment
        matrices%nodal_sol_incr_full=MATMUL(para%ROM_modes,matrices%nodal_sol_incr)
    END IF
    
    DO p=1,para%nRows_reduced  !update the nodal solution for all independent DOF
        matrices%nodal_sol_t_1(para%reduced_dof_to_global_dof(1,p))=&
        matrices%nodal_sol_t_1(para%reduced_dof_to_global_dof(1,p))+&
        matrices%nodal_sol_incr_full(p,1)
    END DO
    
    !enforce the constraint
    CALL enforce_constraint(para,matrices%nodal_sol_t_1)
    
    END SUBROUTINE update_nodal_sol
    
!==============================================================================
    
    SUBROUTINE assemble(para,matrices,GP_DATA,analysis,PNEWDT,nl_geom,DTIME,&
                        TIME,KSTEP,KINC,max_reaction_elem)
    
    !this subroutine assembles the global matrices of derivatives, global rhs's,
    !and all element state variables by looping over all elements and
    !calling the UEL (which itself calls the UMAT at each GP)
    
    IMPLICIT NONE
    
    !this type bundles all mesh parameters
    TYPE(meshparameters),INTENT(INOUT)::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------
    REAL(KIND=AbqRK):: PARAMS(0),ADLMAG(0,0),ENERGY(8),PREDEF(2,0,3),DDLMAG(0,0),&
                      PERIOD,PNEWDT_in
    INTEGER(KIND=AbqIK):: NPREDF,LFLAGS(5),MDLOAD,NDLOAD,MCRD
    INTEGER(KIND=AbqIK):: 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
    
    !-------------------------------------------------------------------
    !running indices
    INTEGER:: k,l,i,j,e,elem
    !==0 small deformation theory ==1 large deformation theory
    INTEGER,INTENT(IN):: nl_geom
    !maximum reaction needed for convergence criterion
    REAL(KIND=AbqRK),INTENT(OUT),DIMENSION(para%ndof_n_max):: max_reaction_elem
    
    !---------------------initialize these variables---------------------------
    
    PNEWDT=2.0_AbqRK
    LFLAGS=[1_AbqIK,nl_geom,1_AbqIK,0_AbqIK,0_AbqIK]
    MCRD=para%dimens
    
    !----set all arrays and matrices to zero before building them up again-----
    
    IF (analysis%ROM_projection .AND. (.NOT. analysis%hyperintegration)) THEN
        forall(i=1:2,j=1:2) matrices%derivatives_full(i,j)%matrix=0.0_AbqRK
        forall(i=1:2) matrices%rhs_full(i)%matrix=0.0_AbqRK
    ELSE
        forall(i=1:2,j=1:2) matrices%derivatives(i,j)%matrix=0.0_AbqRK
        forall(i=1:2) matrices%rhs(i)%matrix=0.0_AbqRK
    END IF
    
    max_reaction_elem=0.0_AbqRK
    
    !--loop over all elements to build global matrices of derviatives and rhs--
    
    DO elem=1,para%n_active_elements !loop over all active elements
        
        k=para%active_elements(1,elem)
        e=para%active_elements(2,elem)
        
        !------read element coordinates from global Coordinates_global---------
        
        DO i=1,para%Elements(k)%NNODE
            para%Elements(k)%COORDS(:,i)=&
            para%Coordinates_global(:,para%Elements(k)%element_to_node(i,e))
        END DO
        
        !------read the element nodal solution vector and its increment--------
        
        DO i=1,para%Elements(k)%NDOFEL
            para%Elements(k)%U(i)=matrices%nodal_sol_t_1(para%Elements(k)%element_dof_to_global_dof(i,e))
            para%Elements(k)%DU(i,1)=para%Elements(k)%U(i)-matrices%nodal_sol_t(para%Elements(k)%element_dof_to_global_dof(i,e))
        END DO
        
        !------------------call the element routine UEL--------------------------
        
        !read internal state variables of the last time step
        para%Elements(k)%SVARS=GP_DATA%SVARS_t(para%Elements(k)%local_to_global_SVARS(e):para%Elements(k)%local_to_global_SVARS(e)+para%Elements(k)%NSVARS-1)
        para%Elements(k)%RHS=0.0_AbqRK !initialize RHS
        para%Elements(k)%AMATRX=0.0_AbqRK !initialize AMATRX
        
        !call user element which has to be compiled into MonolithFE2
        CALL UEL(para%Elements(k)%RHS,para%Elements(k)%AMATRX,para%Elements(k)%SVARS,ENERGY,para%Elements(k)%NDOFEL,1,&
                 para%Elements(k)%NSVARS,para%Elements(k)%PROPS,para%Elements(k)%NPROPS,para%Elements(k)%COORDS,MCRD,&
                 para%Elements(k)%NNODE,para%Elements(k)%U,para%Elements(k)%DU,para%Elements(k)%V,para%Elements(k)%A,&
                 para%Elements(k)%JTYPE,TIME,DTIME,KSTEP,KINC,para%Elements(k)%JELEM(e),PARAMS,NDLOAD,JDLTYP,ADLMAG,PREDEF,NPREDF,LFLAGS,&
                 para%Elements(k)%NDOFEL,DDLMAG,MDLOAD,PNEWDT_in,para%Elements(k)%JPROPS,para%Elements(k)%NJPROP,PERIOD)
        
        !return the lowest encountered value of PNEWDT
        IF (PNEWDT_in<PNEWDT) PNEWDT=PNEWDT_in
        
        !in the first timestep save the reference volume if in hyper ROM training mode
        IF (analysis%training) THEN
            IF (analysis%ROM_projection) THEN
                CALL store_hyper_intergation_data(para,GP_DATA,TIME,elem,e,k)
            END IF
        END IF
        
        !save state variables for the next timestep t+1
        GP_DATA%SVARS_t_1(para%Elements(k)%local_to_global_SVARS(e):para%Elements(k)%local_to_global_SVARS(e)+para%Elements(k)%NSVARS-1)=para%Elements(k)%SVARS
        
        !get the absolut largest rhs entry (to be compared to the largest variable of the residual) -> convergence check
        DO i=1,para%Elements(k)%NDOFEL
            max_reaction_elem(para%Elements(k)%element_dof(2,i))=MAX(max_reaction_elem(para%Elements(k)%element_dof(2,i)),&
                                                                     ABS(para%Elements(k)%RHS(i,1)))
        END DO
        
        !------------sort AMATRX values into global derivatives array----------
        !------------------sort RHS values into global rhs array---------------
        
        !take symmetric part of AMATRX, if symmetric_matrix==.TRUE.
        IF (analysis%symmetric_matrix) para%Elements(k)%AMATRX=0.5_AbqRK*(para%Elements(k)%AMATRX+TRANSPOSE(para%Elements(k)%AMATRX))
        
        IF (analysis%hyperintegration) THEN
            CALL assemble_ROM_derivative_matrix_and_RHS(para,analysis,matrices,k,e)
        ELSE
            CALL assemble_full_derivative_matrix_and_RHS(para,analysis,matrices,k,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%ROM_projection .AND. (.NOT. analysis%hyperintegration)) THEN
        CALL ROM_projection(matrices,para,analysis)
    END IF
    
    IF (analysis%symmetric_matrix) THEN !use the fact that the stiffness matrix is symmetric
        matrices%derivatives(1,2)%matrix=TRANSPOSE(matrices%derivatives(2,1)%matrix)
    END IF
    
    END SUBROUTINE assemble
    
    !==============================================================================
    
    SUBROUTINE store_hyper_intergation_data(para,GP_DATA,TIME,elem,e,k)
    !this subroutine stores the hyper integration training data in the hyper ROM training
    !simulation containing of modal internal forces and additional measures (e.g. stress
    !power) stored in the state variables array
    !this type bundles all mesh parameters
    TYPE(meshparameters),INTENT(IN)::para
    !this type bundles all macro GP DATA (UGLOBAL,STATEV,...)
    TYPE(macro_GP_DATA),INTENT(INOUT)::GP_DATA
    REAL(KIND=AbqRK), DIMENSION(2),INTENT(IN):: TIME
    INTEGER:: elem,e,k
    
    IF (TIME(2)==0.0_AbqRK) THEN
        para%ref_Volume_elements(elem)=para%Elements(k)%SVARS(para%Elements(k)%NSVARS-para%n_additional_hyper_outputs)
    ELSE
        GP_DATA%training_matrix(elem,1:para%n_ROM_modes)=MATMUL(TRANSPOSE(para%Elements(k)%ROM_modes(:,:,e)),para%Elements(k)%RHS(:,1))
        GP_DATA%training_matrix(elem,para%n_ROM_modes+1:para%n_ROM_modes+para%n_Reaction_force_dof)=MATMUL(TRANSPOSE(para%Elements(k)%Equations(:,:,e)),para%Elements(k)%RHS(:,1))
        GP_DATA%training_matrix(elem,para%n_ROM_modes+para%n_Reaction_force_dof+1:)=para%Elements(k)%SVARS(para%Elements(k)%NSVARS-para%n_additional_hyper_outputs:)
    END IF
    
    END SUBROUTINE store_hyper_intergation_data
    
    !==============================================================================
    
    SUBROUTINE assemble_full_derivative_matrix_and_RHS(para,analysis,matrices,k,e)
    !this subroutine assembles the full matrices of derivatives and RHS's, by sorting the local values
    !into the global arrays
    
    IMPLICIT NONE
    
    !this type bundles all mesh 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
    !element type and element number
    INTEGER,INTENT(IN)::k,e
    !running index
    INTEGER::i
    
    !sort in the RHS values
    DO i=para%Elements(k)%rhs_to_global_column(e),para%Elements(k)%rhs_to_global_column(e+1)-1
        matrices%rhs_full(para%Elements(k)%rhs_to_global(2,i))%matrix(para%Elements(k)%rhs_to_global(3,i),1)=&
        matrices%rhs_full(para%Elements(k)%rhs_to_global(2,i))%matrix(para%Elements(k)%rhs_to_global(3,i),1)+&
        para%Elements(k)%rhs_to_global_factor(i)*para%Elements(k)%RHS(para%Elements(k)%rhs_to_global(1,i),1)
    END DO
    
    !sort in the AMATRX values
    DO i=para%Elements(k)%AMATRX_to_global_column(e),para%Elements(k)%AMATRX_to_global_column(e+1)-1
        matrices%derivatives_full(para%Elements(k)%AMATRX_to_global(3,i),para%Elements(k)%AMATRX_to_global(4,i))%matrix(para%Elements(k)%AMATRX_to_global(5,i),para%Elements(k)%AMATRX_to_global(6,i))=&
        matrices%derivatives_full(para%Elements(k)%AMATRX_to_global(3,i),para%Elements(k)%AMATRX_to_global(4,i))%matrix(para%Elements(k)%AMATRX_to_global(5,i),para%Elements(k)%AMATRX_to_global(6,i))+&
        +para%Elements(k)%AMATRX_to_global_factor(i)*para%Elements(k)%AMATRX(para%Elements(k)%AMATRX_to_global(1,i),para%Elements(k)%AMATRX_to_global(2,i))
    END DO
    
    END SUBROUTINE assemble_full_derivative_matrix_and_RHS
    
    !==============================================================================
    
    SUBROUTINE assemble_ROM_derivative_matrix_and_RHS(para,analysis,matrices,k,e)
    !this subroutine assembles the ROM matrices of derivatives and RHS's using the ROM modes and Galerkin projection, by directly sorting
    !the AMATRX entries to the respectiv ROM entries
    
    IMPLICIT NONE
    
    !this type bundles all mesh parameters
    TYPE(meshparameters),INTENT(INOUT)::para
    !this type bundles all FEM matrices
    TYPE(systemmatrices),INTENT(INOUT):: matrices
    !this type bundles all convergence/analysis parameters
    TYPE (analysisparameters),INTENT(IN)::analysis
    !!element type and element number
    INTEGER,INTENT(IN)::k,e
    
    !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%Elements(k)%NDOFEL,para%n_ROM_modes,1.0_AbqRK,para%Elements(k)%AMATRX,para%Elements(k)%NDOFEL,&
                   para%Elements(k)%ROM_modes(:,:,e),para%Elements(k)%NDOFEL,0.0_AbqRK,para%Elements(k)%temp_matrix1,para%Elements(k)%NDOFEL)
    ELSE
        para%Elements(k)%temp_matrix1=MATMUL(para%Elements(k)%AMATRX,para%Elements(k)%ROM_modes(:,:,e))
    END IF
    
    !now compute Phi_elem^T*(K_elem*Phi_elem)
    matrices%derivatives(1,1)%matrix=matrices%derivatives(1,1)%matrix+MATMUL(TRANSPOSE(para%Elements(k)%ROM_modes(:,:,e)),para%Elements(k)%temp_matrix1)*para%Elements(k)%multiplication_factor(e)
    
    !build up global rhs (ROM reduced)
    matrices%rhs(1)%matrix=matrices%rhs(1)%matrix+MATMUL(TRANSPOSE(para%Elements(k)%ROM_modes(:,:,e)),para%Elements(k)%RHS)*para%Elements(k)%multiplication_factor(e)
    
    IF (para%Elements(k)%constraint_element_dof(e)==1) THEN !if the element is constrained to macro dof
        
        !now compute A_elem^T*(K_elem*Phi_elem) using the temporary matrix
        matrices%derivatives(2,1)%matrix=matrices%derivatives(2,1)%matrix+MATMUL(TRANSPOSE(para%Elements(k)%Equations(:,:,e)),para%Elements(k)%temp_matrix1)*para%Elements(k)%multiplication_factor(e)
        
        IF (analysis%symmetric_matrix) THEN
            !compute K_elem*A_elem and store it in a temporary matrix
            CALL dsymm('L','U',para%Elements(k)%NDOFEL,para%n_Reaction_force_dof,1.0_AbqRK,para%Elements(k)%AMATRX,para%Elements(k)%NDOFEL,&
                       para%Elements(k)%Equations(:,:,e),para%Elements(k)%NDOFEL,0.0_AbqRK,para%Elements(k)%temp_matrix2,para%Elements(k)%NDOFEL)
        ELSE
            !compute K_elem*A_elem and store it in a temporary matrix
            para%Elements(k)%temp_matrix2=MATMUL(para%Elements(k)%AMATRX,para%Elements(k)%Equations(:,:,e))
            !Phi_elem^T*(K_elem*A_elem)
            matrices%derivatives(1,2)%matrix=matrices%derivatives(1,2)%matrix+MATMUL(TRANSPOSE(para%Elements(k)%ROM_modes(:,:,e)),para%Elements(k)%temp_matrix2)*para%Elements(k)%multiplication_factor(e)
        END IF
        
        !now compute A_elem^T*(K_elem*A_elem)
        matrices%derivatives(2,2)%matrix=matrices%derivatives(2,2)%matrix+MATMUL(TRANSPOSE(para%Elements(k)%Equations(:,:,e)),para%Elements(k)%temp_matrix2)*para%Elements(k)%multiplication_factor(e)
        
        !compute macro response
        matrices%rhs(2)%matrix=matrices%rhs(2)%matrix+MATMUL(TRANSPOSE(para%Elements(k)%Equations(:,:,e)),para%Elements(k)%RHS)*para%Elements(k)%multiplication_factor(e)
        
    END IF
    
    END SUBROUTINE assemble_ROM_derivative_matrix_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 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(1)%matrix=MATMUL(TRANSPOSE(para%ROM_modes),matrices%rhs_full(1)%matrix)
    matrices%derivatives(2,1)%matrix=MATMUL(matrices%derivatives_full(2,1)%matrix,para%ROM_modes)
    IF (.NOT. analysis%symmetric_matrix) matrices%derivatives(1,2)%matrix=MATMUL(TRANSPOSE(para%ROM_modes),matrices%derivatives_full(1,2)%matrix)
    
    !create a handle containing the unprojected matrix of derivatives (micro-micro) 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%derivatives_full(1,1)%matrix)
    
    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%derivatives(1,1)%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,macro_response,matrices,solver,analysis)
    !this subroutine performs the static condensation by computing DDSDDE
    !(matrix of derivative of macro response w.r.t macro measures)
    
    IMPLICIT NONE
    
    !this type bundles all mesh parameters
    TYPE(meshparameters),INTENT(IN)::para
    !this type bundles all convergence/analysis parameters
    TYPE (analysisparameters)::analysis
    !macro response(e.g. 1.Piola Kirchhoff Stress)
    REAL(KIND=AbqRK), DIMENSION(:),INTENT(INOUT):: macro_response
    !macro tangent "stiffness"
    REAL(KIND=AbqRK), DIMENSION(:,:),INTENT(OUT):: DDSDDE
    !this type bundles all FEM matrices
    TYPE(systemmatrices):: matrices
    !partial derivative of micro nodal unknowns w.r.t macro deformation gradient
    REAL(KIND=AbqRK), DIMENSION(:,:),ALLOCATABLE:: d_micro_nodal_sol_d_macro_measure
    !partial derivative of macro stress w.r.t micro reaction force
    REAL(KIND=AbqRK), DIMENSION(:,:),ALLOCATABLE::d_macro_response_d_micro_residual
    !running indice
    INTEGER:: i,j
    TYPE(solve_soe):: solver
    
    !-----------the macro response was previously stored in rhs(2)-------------
    
    macro_response=-matrices%rhs(2)%matrix(:,1) !negative sign because of definition of element RHS
    
    !-------------compute d_micro_nodal_sol_d_macro_measure--------------------
    
    ALLOCATE(d_micro_nodal_sol_d_macro_measure(SIZE(matrices%rhs(1)%matrix,1),para%n_Reaction_force_dof))
    d_micro_nodal_sol_d_macro_measure=0.0_AbqRK
    CALL solver%solve(para,analysis,matrices%derivatives(1,1)%matrix,&
            matrices%derivatives(1,2)%matrix,d_micro_nodal_sol_d_macro_measure)

    !-------------------------compute DDSDDE-----------------------------------    
        
    DDSDDE=matrices%derivatives(2,2)%matrix-&
      MATMUL(matrices%derivatives(2,1)%matrix,d_micro_nodal_sol_d_macro_measure)

    !---compute algorithmically consistent macro response in monolithic case---
    
    IF (analysis%monolithic) THEN
        
        IF (analysis%symmetric_matrix) THEN !matrices%derivatives(2,1)==matrices%derivatives(1,2)
            macro_response=macro_response+&
            MATMUL(TRANSPOSE(d_micro_nodal_sol_d_macro_measure),matrices%rhs(1)%matrix(:,1))
        ELSE
            !Allocate Matrix of derivative of macro repsonse w.r.t. micro residual
            ALLOCATE(d_macro_response_d_micro_residual(SIZE(matrices%rhs(1)%matrix,1),para%n_Reaction_force_dof))
            
            CALL solver%solve(para,analysis,matrices%derivatives(1,1)%matrix,&
                    TRANSPOSE(matrices%derivatives(2,1)%matrix),d_macro_response_d_micro_residual)
            
            macro_response=macro_response+&
            MATMUL(TRANSPOSE(d_macro_response_d_micro_residual),matrices%rhs(1)%matrix(:,1))
            
        END IF
        
    END IF
    
    !----------------divide by the RVEs (Area=2D, Volume=3D)-------------------
    
    DDSDDE=DDSDDE/para%RVE_Volume
    macro_response=macro_response/para%RVE_Volume
        
END SUBROUTINE static_condensation

END MODULE MonolithFE2
