!DEC$ FREEFORM
!===============================================================
! Mechanical Element for Abaqqus interface UEL
! using shape functions etc. from UELlib
! based on example element in UELlib by R. Skypnyk, G. Huetter and S. Roth
!
! 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
!
! features: 
! - small and large deformation analyses
! - Update Lagrange Formulation 
! - Stress rotation algorithm of Hughes and Winget
!
!-------------------------------------------------------------------
! Version history:
!
! Geralf Hütter, 2020-11-10: first version
! Nils Lange, 24.11.2020
! - AMATRX corrected (VolGP)
! - return the symmetric part of AMATRX
! - return lowest value of PNEWDT from all GP UMAT calls
! - geometrical stiffness added
! Nils Lange, 10.03.2021
! - do not return symmetric part of AMATRX, as the symmetrization
!   is done by the main program
! - all geometrical terms removed
!===============================================================

SUBROUTINE UEL(RHS_ABQ,AMATRX,SVARS,ENERGY,NDOFEL,PROPS,NPROPS,COORDS0,MCRD,NNODE,U,DU_ABQ,V,A,JTYPE,TIME,DTIME,KSTEP,KINC,JELEM, &
               PARAMS,NDLOAD,JDLTYP,ADLMAG,PREDEF,NPREDF,LFLAGS,MLVARX,DDLMAG,MDLOAD,PNEWDT_ABQ,JPROPS,NJPROP,PERIOD)
     USE ABQINTERFACE
     USE MATH
     IMPLICIT NONE
     ! Declaration of interface of UEL (See paragraph 1.1.27 of Abaqus User Subroutines Reference Manual).
     INTEGER(KIND=AbqIK):: NDOFEL,NPROPS,NNODE,JTYPE,KSTEP,KINC,JELEM,NDLOAD,NJPROP,MCRD,NPREDF,MLVARX,MDLOAD,JDLTYP(:,:),LFLAGS(*),JPROPS(:)
     REAL(KIND=AbqRK):: DTIME,PNEWDT,PNEWDT_ABQ,PERIOD
     REAL(KIND=AbqRK):: RHS_ABQ(:,:),AMATRX(:,:),PROPS(:), SVARS(:),ENERGY(8),COORDS0(:,:),U(:),&
                        DU_ABQ(:,:),V(:),A(:),TIME(2),PARAMS(*),ADLMAG(MDLOAD,*),DDLMAG(MDLOAD,*), &
                        PREDEF(2,NPREDF,NNODE)
     ! Declaration of interface of UMAT to call it afterwards (See paragraph 1.1.40 of Abaqus User Subroutines Reference Manual).
     ! Note that some arguments are already declared in UEL declaration section.
     INTEGER(KIND=AbqIK),PARAMETER:: NTENS = NDI+NSHR ! Size of the stress or strain component array.
     INTEGER(KIND=AbqIK):: NSTATV,NSTATVTOTPIP, & ! Number of solution-dependent state variables that are associated with UMAT
                           NOEL, & ! Element number.
                           LAYER, & ! Layer number (for composite shells and layered solids).
                           KSPT ! Section point number within the current layer.
     REAL(KIND=AbqRK):: SSE, SPD, SCD, & ! Specific elastic strain energy, plastic dissipation, and “creep” dissipation, respectively.
                        RPL, & ! Volumetric heat generation per unit time.
                        DRPLDT, & ! Variation of RPL with respect to the temperature.
                        TEMP, & ! Temperature at the start of the increment.
                        DTEMP, & ! Increment of temperature.
                        CELENT ! Characteristic element length.
     REAL(KIND=AbqRK):: STRESS(NTENS), & ! Stress array, must be updated in UMAT.
                        STRAN(NTENS), & ! Array containing the total strains at the beginning of the increment.
                        DSTRAN(NTENS), & ! Array of strain increments.
                        DDSDDE(NTENS,NTENS), & ! Jacobian matrix of the constitutive model / material tangent
                        DDSDDT(NTENS), & ! Variation of the stress increments with respect to the temperature.
                        DRPLDE(NTENS), & ! Variation of RPL with respect to the strain increments.
                        PREDEF_UMAT(1), & ! Array of interpolated values of predefined field variables at this point.
                                          ! _UMAT added to avoid collision with PREDEF array of UEL subroutine.
                        DPRED(1), & ! Array of increments of predefined field variables.
                        DROT(3,3), & ! Rotation increment matrix.
                        DFGRD0(3,3), & ! Array containing the deformation gradient at the beginning of the increment.
                        DFGRD1(3,3) ! Array containing the deformation gradient at the end of the increment.
     REAL(KIND=AbqRK),ALLOCATABLE:: STATEV(:) ! Array containing the solution-dependent state variables.
     CHARACTER*80 CMNAME ! User-defined material name, left justified.
     !
     ! Declaration of local variables.
     REAL (KIND=AbqRK) :: ShapeFunctionDerivative(NDIM, NNODES),ShapeFunction(NNODES), & ! NNODES - number of nodes (See UEL_lib.*: module Ansatz...).
                          JacobiMatrix(NDIM, NDIM), & ! NDIM - number of dimensions (See UEL_lib.*: module Ansatz...).
                          JacobiMatrixInverse(NDIM,NDIM),JacobiMatrixInverse_STRESS(NDIM,NDIM),AMATRX_GEOM(NNODES,NNODES), &
                          BmatU(NTENS, NDOFEL), BmatStress(NTENS, NDOFEL), BmatUComp(NDIM, NNODES),BmatStressComp(NDIM, NNODES),&
                          STRESSTEMP(NTENS),&
                          COORDS(NDIM,NNODES),COORDSGP(3), & ! Node coordinates in current configuration
                          DU(NDOFEL), &
                          RHS(NDOFEL), & ! Additional array used to avoid confusion with the size of RHS_ABQ array:
                                         ! RHS_ABQ corresponds to RHS(MLVARX,NRHS) from paragraph 1.1.27 of Abaqus User Subroutines Reference Manual, and
                                         ! value of MLVARX is different from expected one. For this reason, it was decided to create a new one, RHS(NDOFEL)
                          Determinant, Determinant0,& ! Determinant of the JacobiMatrix for numerical integration.
                          STRESS33 !stress in thickness direction in plane strain simulations
     INTEGER(KIND=AbqIK) :: IntegrPoint, J,I,N,K,L,& ! Integration point number.
                            SVARSindex ! Coefficient to work with SVARS array
     REAL(KIND=AbqRK) :: radius, energyIP(8),VolGP
         
     !set PNEWDT_aba to a high value
     PNEWDT_ABQ=20.0_AbqRK
     ! Call 'GetNSTATV' function to get a number of solution-dependent state variables
     NSTATV=GetNSTATV(NTENS,NDIM,PROPS)
     NSTATVTOTPIP=NSTATV+NTENS+NTENS
     ! Check if the total (not per integration point) number of state variables in the input file is big enough.
     if ( TIME(2).eq.0 ) then! Do it if current value of total time is zero.
        if ( size(SVARS)<NGP*NSTATVTOTPIP ) then
            CALL STDB_ABQERR(-3,'A larger number of solution-dependent state variables is required. Please note that this is a total number of SDVs, not per integration point.',JELEM,TIME,CMNAME)
        elseif (NNODE.ne.NNODES) then
            CALL STDB_ABQERR(-3,'Element type U%I requires %I nodes but %I were assigner to element %I.',(/JTYPE, NNODES, NNODE,JELEM/),TIME,CMNAME)
        elseif (NDOFEL.ne.NDIM*NNODES) then
            CALL STDB_ABQERR(-3,'Element type U%I requires %I DOFs but %I were assigned to element %I.',(/JTYPE, NDIM*NNODE, NDOFEL,JELEM/),TIME,CMNAME)
        endif
        call CheckMaterialParameters(PROPS)
     end if
     ! Allocate array STATEV, if NSTATV /= 0
     if (NSTATV.ne.0) then
        allocate(STATEV(NSTATV)) ! Store NSTATV state variables @ each integration point.
     end if
     ! Initialize RHS and LHS.
     RHS = 0.0_AbqRK
     AMATRX = 0.0_AbqRK
     DU = DU_ABQ(1:NDOFEL,1)
     ! Compute nodal coordinates for initial or current configuration
     select case (LFLAGS(2))
      case (0) !GEOMETRICALLY LINEAR ANALYSIS
        COORDS=COORDS0
      case (1) ! GEOMETRICALLY NONLINEAR ANALYSIS, UPDATED LAGRANGIAN
        COORDS=COORDS0+reshape(U,(/NDIM,NNODES/))
     end select
     ! Loop over Gauss-points.
     DO IntegrPoint=1, NGP ! NGP - number of Gauss-points. See UEL_lib.*: InteGETNOSTATEVgr... module.
          ! Restore stresses from previous time increment
          SVARSindex = NSTATVTOTPIP*(IntegrPoint-1)
          STRESS = SVARS( SVARSindex+1:SVARSindex+NTENS )
          SVARSindex = SVARSindex+NTENS 
          STRAN = SVARS( SVARSindex+1:SVARSindex+NTENS )
          SVARSindex = SVARSindex+NTENS 
          STATEV = SVARS( SVARSindex+1:SVARSindex+NSTATV)	

          ! Call 'ShapeFuncDeriv' function (See UEL_lib.*) to calculate derivative of the shape function.
          ShapeFunctionDerivative = ShapeFuncDeriv(GPPos(IntegrPoint,:)) ! GPPos - array with Gauss-points' coordinates (See UEL_lib.*: GAUSS quadrature module).
          ShapeFunction=ShapeFunc(GPPos(IntegrPoint,:))
          ! Call 'Jacobi' function (See Math.*) to calculate Jacobi matrix.
          JacobiMatrix = Jacobi(ShapeFunctionDerivative, COORDS)
          ! Call 'Inverse' subroutine (See Math.*) to calculate Jacobi matrix inverse and determinant.
          CALL Inverse(JacobiMatrix, JacobiMatrixInverse, Determinant)
          !
          if (NTENS==4 .AND. NDI==3) then !get Determinant0 for integrating STRESS_33 over Element
            if (LFLAGS(2)==0) then
				Determinant0=Determinant
            else
				Determinant0=FindDet(Jacobi(ShapeFunctionDerivative, COORDS0))
            end if
          end if
          !
          !Volume of current Gauss point
          VolGP=GPWeight(IntegrPoint)*Determinant
          ! (Current) coordinates of GP
          COORDSGP=0.0_AbqRK
          COORDSGP(1:NDIM)=MATMUL(COORDS,ShapeFunction)
          ! if axi- or spheri-symmetric element
          IF (SYM) THEN
            radius        = DOT_PRODUCT(COORDS(1,:), ShapeFunction(:))
            VolGP   = VolGP * PreFactor(radius)
          ELSE
            radius = 0.0_AbqRK
          END IF
          !
          ! Call 'BMatTensSym' function to get B-matrix for strain calculation. See UEL_simple...f* file to find out more about BMatTensSym.
		  BmatStressComp = BMatScal(JacobiMatrixInverse, ShapeFunctionDerivative,ShapeFunction,radius) !gradient of shape functions dN_I/dx_i
          BmatStress=BMatTensSym_SFG(BmatStressComp,ShapeFunction,radius) ! put these quantities into conventional B-matrix format
          !BmatStress = BMatTensSym(JacobiMatrixInverse_STRESS, ShapeFunctionDerivative, ShapeFunction, radius)
          select case (LFLAGS(2))
          case (0) !GEOMETRICALLY LINEAR ANALYSIS
            ! B-matrices in weak form and for strains coincide
            BMatU=BMatStress
            ! no rotation			
            DROT=reshape((/1,0,0,0,1,0,0,0,1/),(/3,3/))
            ! deformation gradients
            DFGRD0=reshape((/1,0,0,0,1,0,0,0,1/),(/3,3/))
            DFGRD1=reshape((/1,0,0,0,1,0,0,0,1/),(/3,3/))
          case (1) ! GEOMETRICALLY NONLINEAR ANALYSIS, UPDATED LAGRANGIAN
            ! integration after Hughes & Winget (Abaqus Theory Manual, 3.2.2)
            JacobiMatrix = Jacobi(ShapeFunctionDerivative, COORDS-reshape(0.5*DU,(/NDIM,NNODES/))) !all derivatives with respect to half increment
            JacobiMatrixInverse=Inverse(JacobiMatrix)
            ! B-matrices w.r.t. half-increment after Hughes & Winget
            BmatUComp = BMatScal(JacobiMatrixInverse, ShapeFunctionDerivative,ShapeFunction,radius) !for components of grad(u) corresponding to dN_I/dx_i
            BMatU=BMatTensSym_SFG(BmatUComp,ShapeFunction,radius) ! put these quantities into conventional B-matrix format
            ! gradient of displacement increment, temporarily stored in DFGRD1
            DFGRD1(NDIM+1:,:)=0.0_AbqRK
            DFGRD1(1:NDIM,NDIM+1:)=0.0_AbqRK
            DFGRD1(1:NDIM,1:NDIM)=matmul(reshape(DU,(/NDIM,NNODES/)) , transpose(BmatUComp))
            ! Compute its antimetric part as rotation increment dW, stored temporarily in DROT
            DROT(1,2)=(DFGRD1(1,2)-DFGRD1(2,1))/2
            DROT(1,3)=(DFGRD1(1,3)-DFGRD1(3,1))/2
            DROT(2,3)=(DFGRD1(2,3)-DFGRD1(3,2))/2
            DROT(2,1)=-DROT(1,2)
            DROT(3,1)=-DROT(1,3)
            DROT(3,2)=-DROT(2,3)
            FORALL (J=1:3) DROT(J,J)=2 ! Store 2*tident+dW temporarily in DROT before computing it actually
            DROT=matmul(inverse(Transpose(DROT)),DROT)
            !compute derivatives of shape functions w.r.t. to initial configuration for deformation gradient for both configurations
            JacobiMatrix = Jacobi(ShapeFunctionDerivative, COORDS0)
            CALL Inverse(JacobiMatrix, JacobiMatrixInverse, Determinant)
            IF (SYM) radius = DOT_PRODUCT(COORDS0(1,:), ShapeFunction(:))
            BmatUComp = BMatScal(JacobiMatrixInverse, ShapeFunctionDerivative,ShapeFunction,radius)
            DFGRD1=0.0_AbqRK
            DFGRD1(1:NDIM,1:NDIM)=matmul(reshape(U,(/NDIM,NNODES/)) , transpose(BmatUComp))
            FORALL (J=1:3) DFGRD1(J,J)=DFGRD1(J,J)+1
            DFGRD0=0.0_AbqRK
            DFGRD0(1:NDIM,1:NDIM)=matmul(reshape(U-DU,(/NDIM,NNODES/)) , transpose(BmatUComp))
            FORALL (J=1:3) DFGRD0(J,J)=DFGRD0(J,J)+1
!dec$ if defined(DEBUG)
        ! Additional output if compiled in debug mode
        write(7,*) "DRot:" !det: ", JacobiMatDet, ":
        write(7,'(3(E13.6,1X))') DROT(1,:)
        write(7,'(3(E13.6,1X))') DROT(2,:)
        write(7,'(3(E13.6,1X))') DROT(3,:)
!dec$ endif
            ! Rotate stresses
            STRESSTEMP=STRESS
            CALL ROTSIG(STRESSTEMP,DROT,STRESS,1_AbqIK,NDI,NSHR)
            STRESSTEMP=STRAN
            CALL ROTSIG(STRESSTEMP,DROT,STRAN,2_AbqIK,NDI,NSHR)
        end select
        ! Compute increment of strains
        DSTRAN = matmul(BMatU, DU )
        !set PNEWDT to 1.0, the UMAT corrects this value if necessary
        PNEWDT=1.0_AbqRK
        ! Call UMAT subroutine to get material tangent & stress.
        CALL UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,RPL,DDSDDT,DRPLDE,DRPLDT,STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF_UMAT,DPRED,CMNAME,&
                  NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDSGP,DROT,PNEWDT,CELENT,DFGRD0,DFGRD1,NOEL,IntegrPoint,LAYER,KSPT,[KSTEP,0,LFLAGS(2),0],KINC)
                          
        !return the lowest PNEWDT value
        if (PNEWDT<PNEWDT_ABQ) PNEWDT_ABQ=PNEWDT
                
        ! Assemble stiffness matrix: [Ke] = summ_over_integr_points: W_i*det|Jacobian|*[B_transpose]*[DDSDDE]*[B].
        AMATRX = AMATRX + VolGP*matmul( transpose(BMatStress), matmul(DDSDDE(1:NTENS,1:NTENS),BMatStress) ) ! GPWeight() - See UEL_lib.*: GAUSS quadrature module.
        ! Assemble RHS vector: [He] = summ_over_integr_points: W_i*det|Jacobian|*[B_transpose]*[Stress]
        RHS = RHS - VolGP*matmul( transpose(BMatStress),STRESS)
        !
        ! energies
        energyIP = 0.0_AbqRK
        energyIP(2) = SSE
        energyIP(3) = SCD
        energyIP(4) = SPD
        !
        if (NTENS==4 .AND. NDI==3) then !compute average stress(3,3) of element as volume integral
            STRESS33=STRESS33+GPWeight(IntegrPoint)*Determinant0*STRESS(3)
        end if
        
        ! update
        ENERGY = ENERGY+VolGP*energyIP
        !
        ! Update SVARS array. NOTE: 1)SVARS is not passed to UMAT;
        SVARSindex = NSTATVTOTPIP*(IntegrPoint-1)
        SVARS( SVARSindex+1:SVARSindex+NTENS )=STRESS
	    SVARSindex = SVARSindex+NTENS
        SVARS( SVARSindex+1:SVARSindex+NTENS )=STRAN+DSTRAN
        SVARSindex = SVARSindex+NTENS
        SVARS( SVARSindex+1:SVARSindex+NSTATV)=STATEV
        
     END DO
     
     !
     if (NTENS==4 .AND. NDI==3) then !compute average stress(3,3) of element as volume integral (for MonolithFEsqr)
            ENERGY(8)=STRESS33
        end if
          
     RHS_ABQ(1:NDOFEL,1) = RHS(:)
     
     !deallocate STATEV array
     if (allocated(STATEV)) deallocate(STATEV)
          
     RETURN
END SUBROUTINE
