!DEC$ FREEFORM
!-----------------------------------------------------------------------------------
!-----------------------------------------------------------------------------------
! actual Element definition
! Geralf Hütter,
! Rostyslav Skrypnyk, 28.02.2014
! Stephan Roth, 18.10.2017 (axi- and spheri-symmetric elements)
!-----------------------------------------------------------------------------------
!-----------------------------------------------------------------------------------
SUBROUTINE UEL(RHS_ABQ,AMATRX,SVARS,ENERGY,NDOFEL,PROPS,NPROPS,COORDS,MCRD,NNODE,U,DU_ABQ,V,A,JTYPE,TIME,DTIME,KSTEP,KINC,JELEM, &
               PARAMS,NDLOAD,JDLTYP,ADLMAG,PREDEF,NPREDF,LFLAGS,MLVARX,DDLMAG,MDLOAD,PNEWDT,JPROPS,NJPROP,PERIOD)
     USE ABQINTERFACE
     USE MATH
     IMPLICIT NONE
     ! INCLUDE 'ABA_PARAM.INC' - comment out implicit declaration
     ! 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,PERIOD
     REAL(KIND=AbqRK):: RHS_ABQ(:,:),AMATRX(:,:),PROPS(:), SVARS(:),ENERGY(8),COORDS(:,:),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), & ! NNODES - number of nodes (See UEL_lib.*: module Ansatz...).
                          JacobiMatrix(NDIM, NDIM), & ! NDIM - number of dimensions (See UEL_lib.*: module Ansatz...).
                          JacobiMatrixInverse(NDIM,NDIM), &
                          BMatrix(NTENS, NDOFEL), &
                          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 ! Determinant of the JacobiMatrix for numerical integration.
     INTEGER(KIND=AbqIK) :: IntegrPoint ! Integration point number.
     REAL(KIND=AbqRK) :: radius, ShapeFunction(NNODES), energyIP(8)

     ! Call 'GetNSTATV' function to get a number of solution-dependent state variables. (See UXMAT1 module for more details).
     NSTATV=GetNSTATV(NTENS,NDIM)
     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_AbqIK,"A bigger number of solution-dependent state variables is required. " // &
                          "Please note that this is a total number of SDV's, not per integration point.", &
                          0_AbqIK, 0.0_AbqRK, " ")
       END IF
       CALL CheckMaterialParameters(PROPS)
     END IF

     ! Initialize RHS and LHS.
     RHS = 0.0_AbqRK
     AMATRX = 0.0_AbqRK
     DU = DU_ABQ(1:NDOFEL,1)
     ! Loop over Gauss-points.
     DO IntegrPoint=1, NGP ! NGP - number of Gauss-points. See UEL_lib.*: InteGETNOSTATEVgr... module.
          !
          ! 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).
          ! 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 axi- or spheri-symmetric element
          IF (SYM) THEN
            ShapeFunction = ShapeFuncDispl(GPPos(IntegrPoint,:))
            radius        = DOT_PRODUCT(COORDS(1,:), ShapeFunction(:))
            Determinant   = Determinant * PreFactor(radius)
          ELSE
            radius = 0.0_AbqRK
            ShapeFunction = 0.0_AbqRK
          END IF
          !
          ! Call 'BMatStrain' function to get B-matrix for strain calculation. See UEL_simple...f* file to find out more about BMatStrain.
          BMatrix = BMatStrain(JacobiMatrixInverse, ShapeFunctionDerivative, ShapeFunction, radius)
          !
          ! Compute strain at time 't' for UMAT (IMPORTANT: array U stores displacements at 't+dt')
          STRAN = MATMUL(BMatrix, (U - DU) )
          ! Compute increment of strain
          DSTRAN = MATMUL(BMatrix, DU )
          ! Recover state variables and stress from last increment
          STRESS=SVARS( 1+(IntegrPoint-1)*NSTATVTOTPIP:IntegrPoint*NSTATVTOTPIP-NSTATV-NTENS )
          ! Call UMAT subroutine to get material tangent & stress.
          CALL UXMAT (STRESS,SVARS(IntegrPoint*NSTATVTOTPIP-NSTATV+1:IntegrPoint*NSTATVTOTPIP),DDSDDE,SSE,SPD,SCD,RPL,DDSDDT,DRPLDE,DRPLDT,STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF_UMAT,DPRED,CMNAME, &
                      NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,GPPos(IntegrPoint,:),DROT,PNEWDT,CELENT,DFGRD0,DFGRD1,NOEL,IntegrPoint,LAYER,KSPT,KSTEP,KINC)
          ! Assemble stiffness matrix: [Ke] = summ_over_integr_points: W_i*det|Jacobian|*[B_transpose]*[DDSDDE]*[B].
          AMATRX = AMATRX + GPWeight(IntegrPoint)*Determinant*MATMUL( TRANSPOSE(BMatrix),MATMUL(DDSDDE,BMatrix) ) ! GPWeight() - See UEL_lib.*: GAUSS quadrature module.
          ! Assemble RHS vector: [He] = summ_over_integr_points: W_i*det|Jacobian|*[B_transpose]*[Stress]
          RHS = RHS - GPWeight(IntegrPoint)*Determinant*MATMUL( TRANSPOSE(BMatrix),STRESS )
          !
          ! energies
          energyIP = 0.0_AbqRK
          energyIP(2) = SSE
          energyIP(3) = SCD
          energyIP(4) = SPD
          ! update
          ENERGY = ENERGY+GPWeight(IntegrPoint)*Determinant*energyIP
          !
          ! Update SVARS array (store stress)
          SVARS( 1+(IntegrPoint-1)*NSTATVTOTPIP:IntegrPoint*NSTATVTOTPIP-NSTATV-NTENS ) = STRESS
          SVARS( 1+(IntegrPoint-1)*NSTATVTOTPIP+NTENS:IntegrPoint*NSTATVTOTPIP-NSTATV ) = STRAN+DSTRAN
     END DO
     RHS_ABQ(1:NDOFEL,1) = RHS(:)
     RETURN
END SUBROUTINE
