!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 type_meshparameters
USE ABQINTERFACE
USE class_Solver
USE type_macro_GP_DATA
USE type_analysisparameters

CONTAINS

    SUBROUTINE main_program_staggered(PK1,DDSDDE,Grad_U_k_1,GP_DATA,PNEWDT,&
                                      PROPS,NPROPS,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(para%ndof_macro),INTENT(INOUT):: PK1
    !macro stiffness
    REAL(KIND=AbqRK),DIMENSION(para%ndof_macro,para%ndof_macro),INTENT(OUT):: DDSDDE
    !displacement gradient of this increment
    REAL(KIND=AbqRK),DIMENSION(para%ndof_macro),INTENT(IN)::Grad_U_k_1
    !mechanical properties
    INTEGER,INTENT(IN):: NPROPS
    REAL(KIND=AbqRK),DIMENSION(NPROPS),INTENT(IN):: PROPS
    !If in the staggered scheme the increment fails, this variable is to be set
    !smaller then 1
    REAL(KIND=AbqRK),INTENT(OUT):: PNEWDT
    !displacement increment per NR loop
    REAL(KIND=AbqRK),DIMENSION(para%nRows_reduced):: DELTAU
    !Stiffnessmatrix values
    REAL(KIND=AbqRK),DIMENSION(para%nNonZero_reduced):: k_matrix_values
    !total rhs array
    REAL(KIND=AbqRK),DIMENSION(para%nRows_reduced):: rhs_total
    !partial derivative of macro stress (1.PK) w.r.t micro displacements
    REAL(KIND=AbqRK),DIMENSION(para%ndof_macro,para%nRows_reduced):: d_PK1_d_u
    !partial derivative of reaction force w.r.t macro deformation gradient
    REAL(KIND=AbqRK),DIMENSION(para%nRows_reduced,para%ndof_macro):: d_r_d_F
    !partial derivative of macro stress (1. PK) w.r.t macro deformation gradient
    REAL(KIND=AbqRK),DIMENSION(para%ndof_macro,para%ndof_macro):: d_PK1_d_F
    !loop indices
    INTEGER:: p,iters,i
    !solver type declared in module class_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
    !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
        
    !----------------------------initialize solver-----------------------------

    error=solver%initialize(para,analysis)
                                          
    !---------------------------enforce contraint------------------------------
        
    CALL enforce_constraint(para,Grad_U_k_1,GP_DATA)
        
    !---------------------------------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,k_matrix_values,rhs_total,PK1,d_PK1_d_u,d_r_d_F,&
                    d_PK1_d_F,PROPS,NPROPS,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---------------

      error=solver%factor(para,analysis,k_matrix_values)

      !--------------check if convergence criterion is reached-----------------
      
      IF (convergence_check(para,analysis,MAXVAL(ABS(rhs_total)),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------------------
            error=solver%solve(para,analysis,k_matrix_values,1,rhs_total,DELTAU)

       !get the maximum pos. value of DELTAU of all increments of the time step
            DELTAU_max_i=maxval(abs(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,DELTAU)
            
      ELSE
      
        IF (iters<=analysis%max_iters/2 .AND. PNEWDT>=0.999_AbqRK) THEN
        
            PNEWDT=1.5_AbqRK
            
        ELSE IF (iters>analysis%max_iters/2 .AND. PNEWDT>=0.999_AbqRK) THEN
      
            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) THEN
        
       CALL static_condensation(para,DDSDDE,PK1,k_matrix_values,rhs_total,&
                                d_PK1_d_u,d_r_d_F,d_PK1_d_F,solver,analysis)
                                
     END IF
                             
    !---------------------deallocate solver storage----------------------------
                
    error=solver%finish()
        
    END SUBROUTINE main_program_staggered
    
!==============================================================================
        
    SUBROUTINE main_program_monolithic(PK1,DDSDDE,Grad_U_k_1,GP_DATA,PROPS,&
									   NPROPS,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(para%ndof_macro),INTENT(INOUT):: PK1
    !macro stiffness
    REAL(KIND=AbqRK),DIMENSION(para%ndof_macro,para%ndof_macro),INTENT(OUT):: DDSDDE
    !displacement gradient of this increment
    REAL(KIND=AbqRK),DIMENSION(para%ndof_macro),INTENT(IN)::Grad_U_k_1
    !mechanical properties
    INTEGER,INTENT(IN):: NPROPS
    REAL(KIND=AbqRK),DIMENSION(NPROPS),INTENT(IN):: PROPS
    !displacement increment per NR loop
    REAL(KIND=AbqRK),DIMENSION(para%nRows_reduced):: DELTAU
    !Stiffnessmatrix values
    REAL(KIND=AbqRK),DIMENSION(para%nNonZero_reduced):: k_matrix_values
    !total rhs array
    REAL(KIND=AbqRK),DIMENSION(para%nRows_reduced):: rhs_total
    !partial derivative of macro stress (1.PK) w.r.t micro displacements
    REAL(KIND=AbqRK),DIMENSION(para%ndof_macro,para%nRows_reduced):: d_PK1_d_u
    !partial derivative of reaction force w.r.t macro deformation gradient
    REAL(KIND=AbqRK),DIMENSION(para%nRows_reduced,para%ndof_macro):: d_r_d_F
    !partial derivative of macro stress w.r.t macro strain
    REAL(KIND=AbqRK),DIMENSION(para%ndof_macro,para%ndof_macro):: d_PK1_d_F
    !loop indices
    INTEGER:: p,iters,i,k
    !solver type declared in module class_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

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

    IF ((analysis%save_soe .AND. (GP_DATA%step_indicator(1)==-1)) .OR. (.NOT. analysis%save_soe)) THEN

        error=solver%initialize(para,analysis)
                         
        !save pointer to solver storage
        IF (analysis%save_soe) THEN
            DO k=1,64
                GP_DATA%handle_nbr(k)=solver%handle(k)%dummy
            END DO
        END IF
    
    ELSE
    
        !assigne pointer to solver storage
        DO k=1,64
            solver%handle(k)%dummy=GP_DATA%handle_nbr(k)
        END DO

    END IF
            
    !--------> 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,k_matrix_values,rhs_total,PK1,d_PK1_d_u,d_r_d_F,&
                    d_PK1_d_F,PROPS,NPROPS,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-------------------------

       error=solver%factor(para,analysis,k_matrix_values)

    END IF
    
    !--middle of timestep, factor stiffnessmatrix if factorization isnt saved--
    
    IF (GP_DATA%step_indicator(1)>1 .AND. (.NOT. analysis%save_soe)) THEN

        error=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(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

        error=solver%solve(para,analysis,GP_DATA%k_matrix_values_k,1,GP_DATA%rhs_total_k-MATMUL(GP_DATA%d_r_d_F_k,(Grad_U_k_1-GP_DATA%Grad_U_k)),DELTAU)

      ELSE

        error=solver%solve(para,analysis,k_matrix_values,1,rhs_total-MATMUL(d_r_d_F,(Grad_U_k_1-GP_DATA%Grad_U_k)),DELTAU)

      END IF
          
      !---------------update UGlOBAL: UGlOBAL=UGlOBAL+DELTAU------------------
       
      CALL update_displacements(para,GP_DATA,DELTAU)
       
     !---------------------------enforce contraint-----------------------------
        
      CALL enforce_constraint(para,Grad_U_k_1,GP_DATA)
    
     !-------------------------assemble returns:------------------------------
     !-stiffnessmatrix,righthandside,stress,d_PK1_d_u,d_PK1_d_F,d_r_d_F-
      
     CALL assemble(para,k_matrix_values,rhs_total,PK1,d_PK1_d_u,d_r_d_F,&
                   d_PK1_d_F,PROPS,NPROPS,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----------------

     error=solver%factor(para,analysis,k_matrix_values)
        
    END IF

    !----------save additional information for next iteration step-------------
                    
    GP_DATA%rhs_total_k=rhs_total
    GP_DATA%d_r_d_F_k=d_r_d_F
    GP_DATA%k_matrix_values_k=k_matrix_values
    
    
    !------------call static condensation -> get stiffnessmatrix---------------

    CALL static_condensation(para,DDSDDE,PK1,k_matrix_values,rhs_total,&
                             d_PK1_d_u,d_r_d_F,d_PK1_d_F,solver,analysis)
    
    !-----------------------update step indicator------------------------------
    IF (GP_DATA%step_indicator(1)==-1) THEN
        GP_DATA%step_indicator(1)=2
    ELSE
        GP_DATA%step_indicator(1)=GP_DATA%step_indicator(1)+1
    END IF
    
    20  CONTINUE !continue from this point if PNEWDT<1.0_AbqRK
    
   !---------------------deallocate solver storage----------------------------

    IF (.NOT. analysis%save_soe) error=solver%finish()  
        
    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,GP_DATA)
    !this subroutine enforces the constraine:
    !u_slave_i = u_master_i + Eij * deltax_j
    
    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
    !macro displacement gradient
    REAL(KIND=AbqRK),DIMENSION(para%ndof_macro),INTENT(IN)::Grad_U
    !running indice
    INTEGER::p,i
    
    
    DO p=1,para%nEquations
            
      !u_slave_i_* = u_master_i
      GP_DATA%UGLOBAL_t_1(int(para%Equations(1,p)))=&
      GP_DATA%UGLOBAL_t_1(int(para%Equations(2,p)))
                
      DO i=1,para%ndof_macro
       !u_slave_i_ = u_slave_i_* + Eij * deltax_j
       GP_DATA%UGLOBAL_t_1(int(para%Equations(1,p)))=&
       GP_DATA%UGLOBAL_t_1(int(para%Equations(1,p)))&
       +para%Equations(2+i,p)*Grad_U(i)
       
      END DO
            
    END DO
    
    END SUBROUTINE enforce_constraint

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

    SUBROUTINE update_displacements(para,GP_DATA,DELTAU)
    !this subroutine updates the displacements: UGLOBAL_k_1=UGLOBAL_k+DELTAU
    
    !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
    !displacement increment per NR loop
    REAL(KIND=AbqRK),DIMENSION(para%nRows_reduced):: DELTAU
    
    
       DO p=1,para%nRows_reduced  !update UGlOBAL UGlOBAL=UGlOBAL+DELTAU
                                      !for all nodes except slaves
         GP_DATA%UGLOBAL_t_1(para%master_reduced_to_global(p))=&
         GP_DATA%UGLOBAL_t_1(para%master_reduced_to_global(p))+DELTAU(p)
                    
       END DO
                    
       DO p=1,(para%nreduced_master_slave) !update slave displacements
         GP_DATA%UGLOBAL_t_1(para%global_master_reduced_to_slave(1,p))=&
         GP_DATA%UGLOBAL_t_1(para%global_master_reduced_to_slave(1,p))+&
         DELTAU(para%global_master_reduced_to_slave(2,p))
         
       END DO
    
    
    END SUBROUTINE update_displacements
    
!==============================================================================

    SUBROUTINE assemble(para,k_matrix_values,rhs_total,PK1,d_PK1_d_u,d_r_d_F,&
                        d_PK1_d_F,PROPS,NPROPS,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)
    
    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
    
    !----see in ABAQUS user manual for definition of UEL variabels------
    !number of properties
    INTEGER,INTENT(IN)::NPROPS
    !mechanical properties
    REAL(KIND=AbqRK),DIMENSION(NPROPS),INTENT(IN):: PROPS
    !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_micro
    INTEGER(KIND=AbqIK):: JPROPS(0), 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:: l,i,j,e,n,k,i1,i2,i3,j1,j2,j3
    !boundarys indices for PROPS of Material of the Element
    INTEGER:: lower_boundary_PROPS,upper_boundary_PROPS
    !macro stress (1. Piola Kirchhoff)
    REAL(KIND=AbqRK), DIMENSION(para%ndof_macro),INTENT(OUT):: PK1
    !Stiffnessmatrix values
    REAL(KIND=AbqRK), DIMENSION(para%nNonZero_reduced),INTENT(OUT):: k_matrix_values
    !total rhs array
    REAL(KIND=AbqRK), DIMENSION(para%nRows_reduced),INTENT(OUT):: rhs_total
    !partial derivative of macro stress w.r.t micro displacements
    REAL(KIND=AbqRK), DIMENSION(para%ndof_macro,para%nRows_reduced),INTENT(OUT):: d_PK1_d_u
    !partial derivative of reaction force w.r.t macro strain
    REAL(KIND=AbqRK), DIMENSION(para%nRows_reduced,para%ndof_macro),INTENT(OUT):: d_r_d_F
    !partial derivative of macro stress w.r.t macro strain
    REAL(KIND=AbqRK), DIMENSION(para%ndof_macro,para%ndof_macro),INTENT(OUT):: d_PK1_d_F
    !==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
    
    !---------------------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-----
    
    k_matrix_values=0.0_AbqRK
    d_PK1_d_u=0.0_AbqRK
    d_r_d_F=0.0_AbqRK
    d_PK1_d_F=0.0_AbqRK
    rhs_total=0.0_AbqRK
    PK1=0.0_AbqRK
    max_force_i=0.0_AbqRK
    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 e=1,para%nElem !loop over all elements
         
      !------------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)=GP_DATA%UGLOBAL_t_1(para%element_to_global_dof(i,e))&
                -GP_DATA%UGLOBAL_t(para%element_to_global_dof(i,e))
                
        U(i)=GP_DATA%UGLOBAL_t_1(para%element_to_global_dof(i,e))
      END DO
                  
      !------------------call the element routine UEL--------------------------
      
      SVARS=GP_DATA%STATEV_t(:,e) !read internal state variables
      RHS=0.0_AbqRK
      AMATRX=0.0_AbqRK
        
      IF (para%element_to_material(e)==1) THEN   !get material PROPS of element
        lower_boundary_PROPS=4
        NPROPS_micro=int(PROPS(2))
        upper_boundary_PROPS=3+NPROPS_micro
      ELSE IF (para%element_to_material(e)==2) THEN
        lower_boundary_PROPS=4+int(PROPS(2))
        NPROPS_micro=int(PROPS(3))
        upper_boundary_PROPS=NPROPS
      END IF
              
      !call user element library
            
      CALL UEL(RHS,AMATRX,SVARS,ENERGY,para%ndof_e,1,para%n_STATEV_elem,&
           PROPS(lower_boundary_PROPS:upper_boundary_PROPS),NPROPS_micro,&
           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)
           
      STRESS33=STRESS33+ENERGY(8)
           
      IF (PNEWDT_in<PNEWDT) PNEWDT=PNEWDT_in
                      
      GP_DATA%STATEV_t_1(:,e)=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))

      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
            k_matrix_values(para%values_to_global_reduced(l,e))=&
            k_matrix_values(para%values_to_global_reduced(l,e))+AMATRX(i,j)
          END IF
        END DO
        
        IF (para%element_to_reduced_dof(i,e)/=0) THEN    !==0 indicates boundary
          rhs_total(abs(para%element_to_reduced_dof(i,e)))=&          !condition
          rhs_total(abs(para%element_to_reduced_dof(i,e)))+RHS(i,1)
        END IF
      END DO

      !-sort slave AMATRX and RHS values into d_PK1_d_F,d_PK1_d_u,STRESS-
                        
      IF (ANY (para%element_to_reduced_dof(:,e)<=0)) THEN !only elements with
                                                          !slaves contribute
        DO i=1,para%ndof_e !loop over all entries of AMATRX
          DO j=1,para%ndof_e
            i1=para%element_to_reduced_dof(i,e)
            j1=para%element_to_reduced_dof(j,e)
            i2=para%element_to_global_dof(i,e)
            j2=para%element_to_global_dof(j,e)
            i3=para%global_node_to_slave_node(i2)
            j3=para%global_node_to_slave_node(j2)
           
            IF ((i1<=0) .AND. (j1<=0)) THEN !!Kdd!!

               DO k=1,para%ndof_macro !contribution to d_PK1_d_F
                 DO l=1,para%ndof_macro
                   d_PK1_d_F(k,l)=d_PK1_d_F(k,l)+&
                   AMATRX(i,j)*para%Equations(2+k,i3)*para%Equations(2+l,j3)
                 END DO
                        
                 IF (abs(j1)>0) THEN !contribution to d_PK1_d_u
                   d_PK1_d_u(k,abs(j1))=d_PK1_d_u(k,abs(j1))+&
                   AMATRX(i,j)*para%Equations(2+k,i3)
                   
                 END IF
                   
                 IF (abs(i1)>0) THEN !contribution to d_r_d_F
                   d_r_d_F(abs(i1),k)=d_r_d_F(abs(i1),k)+&
                   AMATRX(i,j)*para%Equations(2+k,j3)
                 END IF
                        
               END DO
           
            ELSE IF ((i1<=0 .AND. .NOT. i3==0) .AND. (j1>0)) THEN !!Kdi!!
              
              DO k=1,para%ndof_macro !contribution to d_PK1_d_u
                d_PK1_d_u(k,j1)=d_PK1_d_u(k,j1)+&
                AMATRX(i,j)*para%Equations(2+k,i3)
              END DO
              
            ELSE IF ((j1<=0 .AND. .NOT. j3==0) .AND. (i1>0)) THEN !!Kid!!
            
              DO k=1,para%ndof_macro !contribution to d_r_d_F
                d_r_d_F(i1,k)=d_r_d_F(i1,k)+&
                AMATRX(i,j)*para%Equations(2+k,j3)
              END DO
            
            END IF
                
          END DO
            
          !build macro STRESS (1.PK) from slave node rhs contribution
          IF (i1<=0 .AND. .NOT. i3==0) THEN  !detect slave
            DO j=1,para%ndof_macro
              PK1(j)=PK1(j)-&
              RHS(i,1)*para%Equations(2+j,i3)
            END DO
            !get maxval of RHS
            IF (ABS(RHS(i,1))>max_force_i) max_force_i=ABS(RHS(i,1))
          END IF
            
        END DO
        
      END IF
    
    END DO
    
    !divde the integrated thickness stress (plane strain) by the area
    DO i=1,para%dimens
        STRESS33=STRESS33/para%delta_x(i)
    END DO
    
    END SUBROUTINE assemble

!==============================================================================
    
    SUBROUTINE static_condensation(para,DDSDDE,PK1,k_matrix_values,rhs_total,&
                                   d_PK1_d_u,d_r_d_F,d_PK1_d_F,solver,analysis)
    !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(para%ndof_macro),INTENT(INOUT):: PK1
    !global reduced stiffnessmatrix in CSR Format
    REAL(KIND=AbqRK), DIMENSION(para%nNonZero_reduced),INTENT(IN):: k_matrix_values
    !macro tangentstiffness
    REAL(KIND=AbqRK), DIMENSION(para%ndof_macro,para%ndof_macro),INTENT(OUT):: DDSDDE
    !partial derivative of macro stress w.r.t micro displacements
    REAL(KIND=AbqRK), DIMENSION(para%ndof_macro,para%nRows_reduced),INTENT(IN):: d_PK1_d_u
    !partial derivative of reaction force w.r.t macro deformation
    REAL(KIND=AbqRK), DIMENSION(para%nRows_reduced,para%ndof_macro),INTENT(IN):: d_r_d_F
    !partial derivative of macro stress w.r.t maro deformation gradient
    REAL(KIND=AbqRK), DIMENSION(para%ndof_macro,para%ndof_macro),INTENT(IN):: d_PK1_d_F
    !partial derivative of micro displacement w.r.t macro deformation gradient
    REAL(KIND=AbqRK), DIMENSION(para%nRows_reduced,para%ndof_macro):: d_u_d_F
    !partial derivative of macro stress w.r.t micro reaction force
    REAL(KIND=AbqRK), DIMENSION(para%nRows_reduced,para%ndof_macro)::d_PK1_d_r
    !total rhs array
    REAL(KIND=AbqRK), DIMENSION(para%nRows_reduced),INTENT(IN):: rhs_total
    !running indice
    INTEGER:: i
    TYPE(solve_soe):: solver
    !error code of solver calls
    INTEGER:: error
        
    !--------------------initialize Arrays to be computed----------------------
    
    DDSDDE=0.0_AbqRK
    d_u_d_F=0.0_AbqRK
    d_PK1_d_r=0.0_AbqRK
    
    !--------------------compute d_u_d_F---------------------
	error=solver%solve(para,analysis,k_matrix_values,para%ndof_macro,d_r_d_F,d_u_d_F)
    !-----------compute algorithmically consistent stress in monolithic case------------      
    IF (analysis%monolithic) THEN 
        error=solver%solve(para,analysis,k_matrix_values,para%ndof_macro,TRANSPOSE(d_PK1_d_u),d_PK1_d_r)
        !in monolithic case STRESS has additional contributions
        PK1=PK1-MATMUL(TRANSPOSE(d_PK1_d_r),rhs_total)    
    END IF
    
    DO i=1,para%dimens !stress has to be divided by area (resp. volume in 3D)
        PK1=PK1/para%delta_x(i)
    END DO
                
    !-------------------------compute DDSDDE-----------------------------------
    
    DDSDDE=d_PK1_d_F-MATMUL(d_PK1_d_u,d_u_d_F)
    
    !------------divide DDSDDE by the RVEs (Area=2D, Volume=3D)----------------
    
    DO i=1,para%dimens
      DDSDDE=DDSDDE/para%delta_x(i)
    END DO
        
END SUBROUTINE static_condensation


END MODULE FE2MODULE
