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

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,extension)
     use ABQINTERFACE
     use MATH
     use ExtendedUELModule
     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
                           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)
                          RHS_ip(NDOFEL), & ! Right hand sinde of an integration point
                          Determinant, Determinant0, Determinant_rel ! Determinant of the JacobiMatrix for numerical integration.
     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
     type(UELextension),optional:: extension
     logical,dimension(NGP):: ip_loop ! States if an integration point is active or not
     logical:: training ! If training data (internal force of IP's) is to be outputed or not
     logical:: BMAT_VOLGP_saved ! If B Matrix and IP Weight's are stored

     !---------Check wether the extension is present, unfortunately due to FORTRAN 77 assumed----------------
     !---------------------size arrays, this cannot be done using present------------------------------------
     if (NJPROP>0) then
       if (JPROPS(1)==1) then
          BMAT_VOLGP_saved=.true.
          training=extension%training
          ip_loop=.false.
           do i=1,NGP
             if (extension%integration_points(i,JELEM)/=0) ip_loop(i)=.true.
           end do
         else
           ip_loop=.true.
           BMAT_VOLGP_saved=.false.
          training=.false.
         end if
      else
        ip_loop=.true.
        BMAT_VOLGP_saved=.false.
       training=.false.
      end if
     !--------------------------------------------------------------------------------------------------------
     
     ! Set PNEWDT_aba to a high value (return the lowest of all encountered in the UMAT)
     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

     ! Check if the total (not per integration point) number of state variables in the input file is big enough.
     if (TIME(2)==0.0_AbqRK) 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) ! Check the material properties
     end if

     ! Allocate array STATEV, if NSTATV /= 0
     if (NSTATV/=0) then
        allocate(STATEV(NSTATV)) ! Store NSTATV state variables @ each integration point.
     end if
     
     ! Initialize RHS (right hand side) and AMATRX (stiffness matix)
     RHS=0.0_AbqRK
     AMATRX=0.0_AbqRK
     
     ! Put DU (increment of displacement) in requiered formmat DU=[DU11,DU12,DU21,DU22,DU31...]
     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 with last timestep configuration as reference!
        COORDS=COORDS0+reshape(U-DU,(/NDIM,NNODES/))
     end select

     ! Loop over Gauss-points.
     do IntegrPoint=1,NGP ! NGP - number of Gauss-points. See UEL_lib.*: InteGETNOSTATEVgr... module. 
    
       if (ip_loop(IntegrPoint)) then ! In reduced cubature mode only compute the active integration points
       
          ! Restore stresses from previous time increment
          SVARSindex=NSTATVTOTPIP*(IntegrPoint-1)
          STRESS=SVARS( SVARSindex+1:SVARSindex+NTENS )
          SVARSindex=SVARSindex+NTENS 
          STATEV=SVARS( SVARSindex+1:SVARSindex+NSTATV)

         if ((.NOT. BMAT_VOLGP_saved) .or. (LFLAGS(2)==1) .or. (TIME(2)==0.0_AbqRK)) then
          ! Call 'ShapeFuncDeriv' function (See UEL_lib.*) to calculate derivative of the shape function.
          ! GPPos - array with Gauss-points' coordinates (See UEL_lib.*: GAUSS quadrature module).
          ShapeFunctionDerivative=ShapeFuncDeriv(GPPos(IntegrPoint,:))
          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 (LFLAGS(2)==1) then
            Determinant0=FindDet(Jacobi(ShapeFunctionDerivative,COORDS0))
            Determinant_rel=Determinant/Determinant0
          else
            Determinant0=Determinant
            Determinant_rel=1.0_AbqRK
          end if
          
          ! Volume of current Gauss point
          if (BMAT_VOLGP_saved) then
             if (.NOT. (extension%IP_Weights(extension%integration_points(IntegrPoint,JELEM))>0.0_AbqRK)) then
                 extension%IP_Weights(extension%integration_points(IntegrPoint,JELEM))=GPWeight(IntegrPoint)*Determinant0
             end if
             VolGP=extension%IP_Weights(extension%integration_points(IntegrPoint,JELEM))
          else
             VolGP=GPWeight(IntegrPoint)*Determinant0
          end if

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

          if (BMAT_VOLGP_saved .and. (LFLAGS(2)==0)) then
            forall(i=1:NDOFEL) extension%B_Matrices(extension%integration_points(IntegrPoint,JELEM),NTENS*(i-1)+1:NTENS*i)=BmatStress(:,i)
          end if

          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 initialized as in abaqus
            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/))
            
            ! In small deformations compute the strains of the old increment
            STRAN=matmul(BMatU,U-DU)

          case (1) ! Geometrically nonlinear analysis, updated lagrangian, with last timestep as reference configuration
          
            ! 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=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
            
            ! Finally compute the rotation increment
            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
            
            ! Rotate stresses
            STRESSTEMP=STRESS
            call ROTSIG(STRESSTEMP,DROT,STRESS,1_AbqIK,NDI,NSHR)
            
            STRAN=0.0_AbqRK ! In geometrically nonlinear analysis, STRAN can't be used and will be always set to zero!

        end select

        else ! When reduced integration is used and B-Matrix and integration point weights are given
            forall(i=1:NDOFEL) BMatU(:,i)=extension%B_Matrices(extension%integration_points(IntegrPoint,JELEM),NTENS*(i-1)+1:NTENS*i)
            BMatStress=BMatU
            VolGP=extension%IP_Weights(extension%integration_points(IntegrPoint,JELEM))
            Determinant_rel=1.0_AbqRK
            STRAN=matmul(BMatU,U-DU)
        end if

        ! 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,JELEM,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*Determinant_rel*matmul( transpose(BMatStress), matmul(DDSDDE(1:NTENS,1:NTENS),BMatStress))

        ! Assemble RHS vector: [He] = summ_over_integr_points: W_i*det|Jacobian|*[B_transpose]*[Stress]
        RHS_ip=matmul(transpose(BMatStress),STRESS)*Determinant_rel
        RHS=RHS-VolGP*RHS_ip

        ! Energies
        energyIP=0.0_AbqRK
        energyIP(2)=SSE
        energyIP(3)=SCD
        energyIP(4)=SPD

        ! Compute average stress(3,3) of element as volume integral (required for MonolithFEsqr)
        if (NTENS==4 .AND. NDI==3 .AND. BMAT_VOLGP_saved) extension%STRESS33=extension%STRESS33+VolGP*STRESS(3)

        ! Update the Energy
        ENERGY=ENERGY+VolGP*energyIP

        ! Update SVARS array. NOTE: SVARS is not passed to UMAT
        SVARSindex=NSTATVTOTPIP*(IntegrPoint-1)
        SVARS(SVARSindex+1:SVARSindex+NTENS)=STRESS
        SVARSindex=SVARSindex+NTENS
        SVARS(SVARSindex+1:SVARSindex+NSTATV)=STATEV

        ! Output the internal force of all integration points when being in training mode
        if (training) extension%F_int_ip(NGP*(JELEM-1)+IntegrPoint,:)=RHS_ip
        
      end if
        
     end do
     
     RHS_ABQ(1:NDOFEL,1)=RHS(:)
     
     ! Deallocate STATEV array
     if (allocated(STATEV)) deallocate(STATEV)
     
end subroutine
