************************************************************************
!
!
! Solution variables (or nodal variables) are the displacements, with the
!  referential magnetic field applied using a predefined field.  DOF's 1,
!  2, and 3 are the displacements.
!
! Material behavior is Neo-Hookean rubber elasticity.
!
! This subroutine is for a three-dimensional 8 node isoparametric
!  element as shown below with 1pt (reduced) or 8pt (full) Gauss
!  integration.
!
! In order to avoid locking for the fully-integrated element, we
!  use the F-bar method of de Souza Neto (1996).
!
!  Mechanical, traction- and pressure-type boundary conditions 
!   may be applied to the dummy mesh using the Abaqus built-in 
!   commands *Dload or *Dsload.
!
!
!  8-node     8-----------7
!  brick     /|          /|       zeta
!           / |         / |       
!          5-----------6  |       |     eta
!          |  |        |  |       |   /
!          |  |        |  |       |  /
!          |  4--------|--3       | /
!          | /         | /        |/
!          |/          |/         O--------- xi
!          1-----------2        origin at cube center
!
!     Face numbering follows:
!       Face 1 = nodes 1,2,3,4
!       Face 2 = nodes 5,8,7,6
!       Face 3 = nodes 1,5,6,2
!       Face 4 = nodes 2,6,7,3
!       Face 5 = nodes 3,7,8,4
!       Face 6 = nodes 4,8,5,1
!
! Provided by Shawn A. Chester August 2017
! Modified by Liu Wang Sep 2019
***********************************************************************
!
! User element statement in the input file:
!  *User Element,Nodes=8,Type=U3,Iproperties=<nJProps>,Properties=<nProps>,Coordinates=3,Unsymm
!  1,2,3
!
!     Material Properties Vector
!     --------------------------------------------------------------
!     Gshear   = props(1) ! Shear modulus
!     Kbulk    = props(2) ! Bulk modulus
!     refMagMoment(1,1) = props(3) !
!     refMagMoment(2,1) = props(4) ! reference magnetic moment vector
!     refMagMoment(3,1) = props(5) !

***********************************************************************
!
      SUBROUTINE UEL(RHS,AMATRX,SVARS,ENERGY,NDOFEL,NRHS,NSVARS,
     +     PROPS,NPROPS,coords,MCRD,NNODE,Uall,DUall,Vel,Accn,JTYPE,
     +     TIME,DTIME,KSTEP,KINC,JELEM,PARAMS,NDLOAD,JDLTYP,ADLMAG,
     +     PREDEF,NPREDF,LFLAGS,MLVARX,DDLMAG,MDLOAD,PNEWDT,JPROPS,
     +     NJPROP,PERIOD)
!
      IMPLICIT NONE
!
!*     VARIABLES DEFINED IN UEL, PASSED BACK TO ABAQUS
!
      REAL(8) :: RHS,AMATRX,SVARS,ENERGY
!
!*     VARIABLES PASSED INTO UEL 
!
      REAL(8) :: PROPS,coords,centroidcoords,Uall,DUall,Vel,Accn,TIME,
     1  DTIME,PARAMS,ADLMAG,PREDEF,DDLMAG,PNEWDT,PERIOD
      INTEGER :: NDOFEL,NRHS,NSVARS,NPROPS,MCRD,NNODE,JTYPE,KSTEP,KINC,
     1  JELEM,NDLOAD,JDLTYP,NPREDF,LFLAGS,MLVARX,MDLOAD,JPROPS,NJPROP
!
      DIMENSION RHS(MLVARX,*),AMATRX(NDOFEL,NDOFEL),PROPS(*),
     1 SVARS(*),ENERGY(8),coords(MCRD,NNODE),centroidcoords(3),Uall(NDOFEL),
     2 DUall(MLVARX,*),Vel(NDOFEL),Accn(NDOFEL),TIME(2),PARAMS(*),
     3 JDLTYP(MDLOAD,*),ADLMAG(MDLOAD,*),DDLMAG(MDLOAD,*),
     4 PREDEF(2,NPREDF,NNODE),LFLAGS(*),JPROPS(*)
!
!
      ! Internally defined variables
      !
      integer nDim,nInt,nIntS


      ! Obtain the number of integration points
      !
      nInt = jprops(1)
!
!     Call user element
      if(jtype.eq.3) then
         !
         ! This is a 3D analysis
         !
         nDim = 3
         call U3D8(RHS,AMATRX,SVARS,ENERGY,NDOFEL,NRHS,NSVARS,
     +        PROPS,NPROPS,coords,MCRD,NNODE,Uall,DUall,Vel,Accn,JTYPE,
     +        TIME,DTIME,KSTEP,KINC,JELEM,PARAMS,NDLOAD,JDLTYP,ADLMAG,
     +        PREDEF,NPREDF,LFLAGS,MLVARX,DDLMAG,MDLOAD,PNEWDT,JPROPS,
     +        NJPROP,PERIOD,
     +        nDim,nInt,nIntS)
!       
      else
         write(*,*) 'Only element U3 allowed, currently jtype=',jtype
         call xit
!
      endif
      !
      ! Done with this element, RHS and AMATRX already returned
      !  as output from the specific element routine called
      !
      return
      end subroutine uel

************************************************************************
************************************************************************

      subroutine U3D8(RHS,AMATRX,SVARS,ENERGY,NDOFEL,NRHS,NSVARS,
     +     PROPS,NPROPS,coords,MCRD,NNODE,Uall,DUall,Vel,Accn,JTYPE,
     +     TIME,DTIME,KSTEP,KINC,JELEM,PARAMS,NDLOAD,JDLTYP,ADLMAG,
     +     PREDEF,NPREDF,LFLAGS,MLVARX,DDLMAG,MDLOAD,PNEWDT,JPROPS,
     +     NJPROP,PERIOD,
     +     nDim,nInt,nIntS)
!
      IMPLICIT NONE
!
!     VARIABLES DEFINED IN UEL, PASSED BACK TO ABAQUS
!
      REAL(8) :: RHS,AMATRX,SVARS,ENERGY
!
!     VARIABLES PASSED INTO UEL 
!
      REAL(8) :: PROPS,coords,Uall,DUall,Vel,Accn,TIME,
     1  DTIME,PARAMS,ADLMAG,PREDEF,DDLMAG,PNEWDT,PERIOD
      INTEGER :: NDOFEL,NRHS,NSVARS,NPROPS,MCRD,NNODE,JTYPE,KSTEP,KINC,
     1  JELEM,NDLOAD,JDLTYP,NPREDF,LFLAGS,MLVARX,MDLOAD,JPROPS,NJPROP
!
      DIMENSION RHS(MLVARX,*),AMATRX(NDOFEL,NDOFEL),PROPS(*),
     1 SVARS(*),ENERGY(8),coords(MCRD,NNODE),Uall(NDOFEL),
     2 DUall(MLVARX,*),Vel(NDOFEL),Accn(NDOFEL),TIME(2),PARAMS(*),
     3 JDLTYP(MDLOAD,*),ADLMAG(MDLOAD,*),DDLMAG(MDLOAD,*),
     4 PREDEF(2,NPREDF,NNODE),LFLAGS(*),JPROPS(*)
!
!     Internally defined variables
      real*8 u(nNode,3),du(nNode,ndofel),phi(nNode),coordsC(mcrd,nNode)
      real*8 uOld(nNode,ndofel),u_t(nNode,ndofel)
      real*8 dPhi(nNode)

      integer i,j,k,l,m,n,nIntPt,nDim,intpt,a1,b1,a11,b11,face
      integer nInt,ii,a12,b12,pe,nSdv,stat,q,nIntV,nIntPtV,p,nIntS
      integer nlSdv,kk,IntS,faceFlag

      real*8 Iden(3,3),Le,Ru(3*nNode,1),tmp,centroidcoords(3)
      real*8 Kuu(3*nNode,3*nNode),sh0(nNode)
      real*8 dshxi(nNode,3),dsh0(nNode,3),dshC0(nNode,3),detMapJ0C
      real*8 Fc_tau(3,3),Fc_t(3,3),detFc_tau,detFc_t,w(nInt),D_tau(3,1)
      real*8 sh(nNode),detMapJ,dsh(nNode,3),detMapJC,Limit,umeror
      real*8 dshC(nNode,3),F_tau(3,3)
      real*8 F_t(3,3),detF_tau,xi(nInt,3),detF,T_tau(3,3),bodyForce(3,1)
      real*8 SpUUMod(3,3,3,3),bodyCharge,Gmat(9,3*nNode),Amat(9,9)
      real*8 BodyForceRes(3*nNode,1),flux
      real*8 G0mat(9,3*nNode),Qmat(9,9),dA,Nvec(1,nNode),AmatPhiU(3,9)
      real*8 xLocal(nIntS),yLocal(nIntS),zLocal(nIntS),wS(nIntS),detF_t
      real*8 Kuphi(3*nNode,nNode),Kphiu(nNode,3*nNode),AmatUPhi(9,3)
      real*8 SpPhiUMod(3,3,3),SpPhiPhiMod(3,3),SpUPhiMod(3,3,3),detMapJ0
      real*8 ER(3,1),magFieldVec(3,1),magFieldMag,Svec(9,1)
	  real*8 gradB(3,3),refMagMoment(3,1)

      real*8 zero,one,two,half,Pi,three,third
      parameter(zero=0.d0,one=1.d0,two=2.d0,half=0.5d0,Pi=3.141592653d0,
     +     three=3.d0,third=1.d0/3.d0)
	 

      ! Identity tensor
      !
      call onem(Iden)


      ! Initialize the residual and tangent matrices to zero.
      !
      Ru  = zero
      Kuu = zero
      Energy = zero


      ! Obtain nodal displacements and magnetic field information
      !
      k = 0
      do i=1,nNode
         do j=1,nDim
            k = k + 1
            u(i,j) = Uall(k)
            du(i,j) = DUall(k,1)
            uOld(i,j) = u(i,j) - du(i,j)
         enddo
      enddo



      ! Obtain current nodal coordinates
      !
	  
      do i=1,nNode
         do j=1,nDim
            coordsC(j,i) = coords(j,i) + u(i,j)
         enddo
      enddo
      !
	  !Obtain the centroid coordinates
	  
	  centroidcoords=Zero
	    do i=1,nDim
         do j=1,nNode
            centroidcoords(i)= coordsC(i,j) +centroidcoords(i)
         enddo
      enddo
	  
	  centroidcoords=centroidcoords/nNode

	  ! Obtain the magnetic field vector using a predefined field
      !  that is implemented through the input file.  Specifically,
      !  "*Field, variable=1" provides the magnetic field amplitude.
      magFieldMag = predef(1,1,1)
	  !write(*,*),'magFieldMag='
	  !write(*,*), magFieldMag
	  !write(*,*),'centroidcoords='
	  !write(*,*), centroidcoords
	  
	  call GetMagField(centroidcoords,magFieldMag,magFieldVec)
      !write(*,*),'magFieldVec='
	  !write(*,*), magFieldVec

         !get gradB by each cloumn
		gradB=zero
        call GetgradB(centroidcoords, magFieldMag, magFieldVec,gradB(:,1),1)
		call GetgradB(centroidcoords, magFieldMag, magFieldVec,gradB(:,2),2)
		call GetgradB(centroidcoords, magFieldMag, magFieldVec,gradB(:,3),3)

	  !write(*,*),'gradB='
	  !write(*,*), gradB
	    refMagMoment(1,1) = props(3)
        refMagMoment(2,1) = props(4)
        refMagMoment(3,1) = props(5)
	  
      ! Impose any time-stepping changes on the increments of
      !  displacement if you want (based on element diagonal)
      !
      Le = dsqrt(((coordsC(1,1)-coordsC(1,7))**two) + 
     +     ((coordsC(2,1)-coordsC(2,7))**two) +
     +     ((coordsC(3,1)-coordsC(3,7))**two))
      !
      do i=1,nNode
         do j=1,nDim
            if(dabs(du(i,j)).gt.10.d0*Le) then
               pnewdt = 0.5
               return
            endif
         enddo
      enddo

      !----------------------------------------------------------------
      ! 
      ! Take this opportunity to perform calculations at the element
      !  centroid.  Here, check for hourglass stabilization and get
      !  the deformation gradient for use in the `F-bar' method.
      !
      ! Reference for the F-bar method:
      !  de Souza Neto, E.A., Peric, D., Dutko, M., Owen, D.R.J., 1996.
      !  Design of simple low order finite elements for large strain
      !  analysis of nearly incompressible solids. International Journal
      !  of Solids and Structures, 33, 3277-3296.
      !
      !
      ! Obtain shape functions and their local gradients at the element
      !  centroid, that means xi=eta=zeta=0.0, and nIntPt=1
      !
      if(nNode.eq.8) then
         call calcShape3DLinear(1,zero,1,sh0,dshxi)
      else
         write(*,*) 'Incorrect number of nodes: nNode.ne.8'
         call xit
      endif


      ! Map shape functions from local to global reference coordinate system
      !
      call mapShape3D(nNode,dshxi,coords,dsh0,detMapJ0,stat)
      if(stat.eq.0) then
         pnewdt = 0.5
         return
      endif



      ! Map shape functions from local to global current coordinate system
      !
      call mapShape3D(nNode,dshxi,coordsC,dshC0,detMapJ0C,stat)
      if(stat.eq.0) then
         pnewdt = 0.5
         return
      endif


      ! Calculate the deformation gradient at the element centroid
      !  at the the beginning and end of the increment for use in 
      !  the `F-bar' method
      !
      Fc_tau = Iden
      Fc_t = Iden
      do i=1,nDim
         do j=1,nDim
            do k=1,nNode
               Fc_tau(i,j) = Fc_tau(i,j) + dsh0(k,j)*u(k,i)
               Fc_t(i,j) = Fc_t(i,j) + dsh0(k,j)*uOld(k,i)
            enddo
         enddo
      enddo
      call mdet(Fc_tau,detFc_tau)
      call mdet(Fc_t,detFc_t)
      !
      ! With the deformation gradient known at the element centroid
      !  we are now able to implement the `F-bar' method later
      !
      !----------------------------------------------------------------




      !----------------------------------------------------------------
      ! Begin the loop over integration points
      !
      ! Obtain integration point local coordinates and weights
      !
      if(nInt.eq.1) then
         call xint3D1pt(xi,w,nIntPt) ! 1-pt integration, nInt=1
      elseif(nInt.eq.8) then
         call xint3D8pt(xi,w,nIntPt) ! 8-pt integration, nInt=8
      else
         write(*,*) 'Invalid number of int points, nInt=',nInt
         call xit
      endif


      ! Loop over integration points
      !
      do intpt=1,nIntPt

         ! Obtain shape functions and their local gradients
         !
         if(nNode.eq.8) then
            call calcShape3DLinear(nIntPt,xi,intpt,sh,dshxi)
         else
            write(*,*) 'Incorrect number of nodes: nNode.ne.8'
            call xit
         endif


         ! Map shape functions from local to global reference coordinate system
         !
         call mapShape3D(nNode,dshxi,coords,dsh,detMapJ,stat)
         if(stat.eq.0) then
            pnewdt = 0.5
            return
         endif


         ! Map shape functions from local to global current coordinate system
         !
         call mapShape3D(nNode,dshxi,coordsC,dshC,detMapJC,stat)
         if(stat.eq.0) then
            pnewdt = 0.5
            return
         endif

         ! Obtain, and modify the deformation gradient at this integration
         !  point.  Modify the deformation gradient for use in the `F-bar'
         !  method if needed.
         !
         F_tau = Iden
         F_t = Iden
         do i=1,nDim
            do j=1,nDim
               do k=1,nNode
                  F_tau(i,j) = F_tau(i,j) + dsh(k,j)*u(k,i)
                  F_t(i,j) = F_t(i,j) + dsh(k,j)*uOld(k,i)
               enddo
            enddo
         enddo
         !
         ! Modify the deformation gradient for the `F-bar' method
         !  only when using the 8 node fully integrated linear
         !  element, do not use the `F-bar' method for any other element
         !
         if((nNode.eq.8).and.(nInt.eq.8)) then
            call mdet(F_tau,detF_tau)
            call mdet(F_t,detF_t)
            F_tau = ((detFc_tau/detF_tau)**third)*F_tau
            F_t = ((detFc_tau/detF_tau)**third)*F_t
         endif
         call mdet(F_tau,detF)



         ! Perform the time integration at this integ. point to compute
         !  all the specific forms and parameters needed for the solution
         !
         call integ(props,nprops,dtime,
     +        F_tau,magFieldVec,
     +        T_tau,
     +        SpUUMod,stat)
         !
         if(stat.eq.0) then
            pnewdt = 0.5
            return
         endif



         ! Compute the matrices that contain the shape functions
         ! and derivatives of the shape functions
         !
         do kk=1,nNode
            Nvec(1,kk) = sh(kk)
         enddo
         !
         !
         Gmat = zero
         do kk=1,nNode
            Gmat(1,1+nDim*(kk-1)) = dshC(kk,1)
            Gmat(2,2+nDim*(kk-1)) = dshC(kk,1)
            Gmat(3,3+nDim*(kk-1)) = dshC(kk,1)
            Gmat(4,1+nDim*(kk-1)) = dshC(kk,2)
            Gmat(5,2+nDim*(kk-1)) = dshC(kk,2)
            Gmat(6,3+nDim*(kk-1)) = dshC(kk,2)
            Gmat(7,1+nDim*(kk-1)) = dshC(kk,3)
            Gmat(8,2+nDim*(kk-1)) = dshC(kk,3)
            Gmat(9,3+nDim*(kk-1)) = dshC(kk,3)
         enddo
         !
         G0mat = zero
         do kk=1,nNode
            G0mat(1,1+nDim*(kk-1)) = dshC0(kk,1)
            G0mat(2,2+nDim*(kk-1)) = dshC0(kk,1)
            G0mat(3,3+nDim*(kk-1)) = dshC0(kk,1)
            G0mat(4,1+nDim*(kk-1)) = dshC0(kk,2)
            G0mat(5,2+nDim*(kk-1)) = dshC0(kk,2)
            G0mat(6,3+nDim*(kk-1)) = dshC0(kk,2)
            G0mat(7,1+nDim*(kk-1)) = dshC0(kk,3)
            G0mat(8,2+nDim*(kk-1)) = dshC0(kk,3)
            G0mat(9,3+nDim*(kk-1)) = dshC0(kk,3)
         enddo




         ! Compute/update the displacement residual vector
         !
         Svec(1,1) = T_tau(1,1)
         Svec(2,1) = T_tau(2,1)
         Svec(3,1) = T_tau(3,1)
         Svec(4,1) = T_tau(1,2)
         Svec(5,1) = T_tau(2,2)
         Svec(6,1) = T_tau(3,2)
         Svec(7,1) = T_tau(1,3)
         Svec(8,1) = T_tau(2,3)
         Svec(9,1) = T_tau(3,3)

		 ! the body force vector may be specified here

         bodyForce =(one/detF)*matmul(gradB,matmul(F_tau, refMagMoment))*0
		 !bodyForce(3,1)=0
		 !bodyForce(2,1)=0
		 
	
		 !write(*,*),'bodyForce='
		 !write(*,*), bodyForce
         BodyForceRes = zero
		 do kk=1,nNode
            BodyForceRes(1+nDim*(kk-1),1) = sh(kk)*bodyForce(1,1)*0
            BodyForceRes(2+nDim*(kk-1),1) = sh(kk)*bodyForce(2,1)*0
            BodyForceRes(3+nDim*(kk-1),1) = sh(kk)*bodyForce(3,1)*0
         enddo
         !
         Ru = Ru + detmapJC*w(intpt)*
     +        (
     +        -matmul(transpose(Gmat),Svec)
     +        + BodyForceRes
     +        )



         ! Compute/update the displacement tangent matrix
         !
         Amat = zero
         Amat(1,1) = SpUUMod(1,1,1,1)
         Amat(1,2) = SpUUMod(1,1,2,1)
         Amat(1,3) = SpUUMod(1,1,3,1)
         Amat(1,4) = SpUUMod(1,1,1,2)
         Amat(1,5) = SpUUMod(1,1,2,2)
         Amat(1,6) = SpUUMod(1,1,3,2)
         Amat(1,7) = SpUUMod(1,1,1,3)
         Amat(1,8) = SpUUMod(1,1,2,3)
         Amat(1,9) = SpUUMod(1,1,3,3)
         Amat(2,1) = SpUUMod(2,1,1,1)
         Amat(2,2) = SpUUMod(2,1,2,1)
         Amat(2,3) = SpUUMod(2,1,3,1)
         Amat(2,4) = SpUUMod(2,1,1,2)
         Amat(2,5) = SpUUMod(2,1,2,2)
         Amat(2,6) = SpUUMod(2,1,3,2)
         Amat(2,7) = SpUUMod(2,1,1,3)
         Amat(2,8) = SpUUMod(2,1,2,3)
         Amat(2,9) = SpUUMod(2,1,3,3)
         Amat(3,1) = SpUUMod(3,1,1,1)
         Amat(3,2) = SpUUMod(3,1,2,1)
         Amat(3,3) = SpUUMod(3,1,3,1)
         Amat(3,4) = SpUUMod(3,1,1,2)
         Amat(3,5) = SpUUMod(3,1,2,2)
         Amat(3,6) = SpUUMod(3,1,3,2)
         Amat(3,7) = SpUUMod(3,1,1,3)
         Amat(3,8) = SpUUMod(3,1,2,3)
         Amat(3,9) = SpUUMod(3,1,3,3)
         Amat(4,1) = SpUUMod(1,2,1,1)
         Amat(4,2) = SpUUMod(1,2,2,1)
         Amat(4,3) = SpUUMod(1,2,3,1)
         Amat(4,4) = SpUUMod(1,2,1,2)
         Amat(4,5) = SpUUMod(1,2,2,2)
         Amat(4,6) = SpUUMod(1,2,3,2)
         Amat(4,7) = SpUUMod(1,2,1,3)
         Amat(4,8) = SpUUMod(1,2,2,3)
         Amat(4,9) = SpUUMod(1,2,3,3)
         Amat(5,1) = SpUUMod(2,2,1,1)
         Amat(5,2) = SpUUMod(2,2,2,1)
         Amat(5,3) = SpUUMod(2,2,3,1)
         Amat(5,4) = SpUUMod(2,2,1,2)
         Amat(5,5) = SpUUMod(2,2,2,2)
         Amat(5,6) = SpUUMod(2,2,3,2)
         Amat(5,7) = SpUUMod(2,2,1,3)
         Amat(5,8) = SpUUMod(2,2,2,3)
         Amat(5,9) = SpUUMod(2,2,3,3)
         Amat(6,1) = SpUUMod(3,2,1,1)
         Amat(6,2) = SpUUMod(3,2,2,1)
         Amat(6,3) = SpUUMod(3,2,3,1)
         Amat(6,4) = SpUUMod(3,2,1,2)
         Amat(6,5) = SpUUMod(3,2,2,2)
         Amat(6,6) = SpUUMod(3,2,3,2)
         Amat(6,7) = SpUUMod(3,2,1,3)
         Amat(6,8) = SpUUMod(3,2,2,3)
         Amat(6,9) = SpUUMod(3,2,3,3)
         Amat(7,1) = SpUUMod(1,3,1,1)
         Amat(7,2) = SpUUMod(1,3,2,1)
         Amat(7,3) = SpUUMod(1,3,3,1)
         Amat(7,4) = SpUUMod(1,3,1,2)
         Amat(7,5) = SpUUMod(1,3,2,2)
         Amat(7,6) = SpUUMod(1,3,3,2)
         Amat(7,7) = SpUUMod(1,3,1,3)
         Amat(7,8) = SpUUMod(1,3,2,3)
         Amat(7,9) = SpUUMod(1,3,3,3)
         Amat(8,1) = SpUUMod(2,3,1,1)
         Amat(8,2) = SpUUMod(2,3,2,1)
         Amat(8,3) = SpUUMod(2,3,3,1)
         Amat(8,4) = SpUUMod(2,3,1,2)
         Amat(8,5) = SpUUMod(2,3,2,2)
         Amat(8,6) = SpUUMod(2,3,3,2)
         Amat(8,7) = SpUUMod(2,3,1,3)
         Amat(8,8) = SpUUMod(2,3,2,3)
         Amat(8,9) = SpUUMod(2,3,3,3)
         Amat(9,1) = SpUUMod(3,3,1,1)
         Amat(9,2) = SpUUMod(3,3,2,1)
         Amat(9,3) = SpUUMod(3,3,3,1)
         Amat(9,4) = SpUUMod(3,3,1,2)
         Amat(9,5) = SpUUMod(3,3,2,2)
         Amat(9,6) = SpUUMod(3,3,3,2)
         Amat(9,7) = SpUUMod(3,3,1,3)
         Amat(9,8) = SpUUMod(3,3,2,3)
         Amat(9,9) = SpUUMod(3,3,3,3)


         Qmat = zero
         Qmat(1,1) = third*(Amat(1,1)+Amat(1,5)+Amat(1,9)) 
     +        - (two/three)*T_tau(1,1)
         Qmat(2,1) = third*(Amat(2,1)+Amat(2,5)+Amat(2,9))
     +        - (two/three)*T_tau(2,1)
         Qmat(3,1) = third*(Amat(3,1)+Amat(3,5)+Amat(3,9))
     +        - (two/three)*T_tau(3,1)
         Qmat(4,1) = third*(Amat(4,1)+Amat(4,5)+Amat(4,9))
     +        - (two/three)*T_tau(1,2)
         Qmat(5,1) = third*(Amat(5,1)+Amat(5,5)+Amat(5,9))
     +        - (two/three)*T_tau(2,2)
         Qmat(6,1) = third*(Amat(6,1)+Amat(6,5)+Amat(6,9))
     +        - (two/three)*T_tau(3,2)
         Qmat(7,1) = third*(Amat(7,1)+Amat(7,5)+Amat(7,9))
     +        - (two/three)*T_tau(1,3)
         Qmat(8,1) = third*(Amat(8,1)+Amat(8,5)+Amat(8,9))
     +        - (two/three)*T_tau(2,3)
         Qmat(9,1) = third*(Amat(9,1)+Amat(9,5)+Amat(9,9))
     +        - (two/three)*T_tau(3,3)
         Qmat(1,5) = Qmat(1,1)
         Qmat(2,5) = Qmat(2,1)
         Qmat(3,5) = Qmat(3,1)
         Qmat(4,5) = Qmat(4,1)
         Qmat(5,5) = Qmat(5,1)
         Qmat(6,5) = Qmat(6,1)
         Qmat(7,5) = Qmat(7,1)
         Qmat(8,5) = Qmat(8,1)
         Qmat(9,5) = Qmat(9,1)
         Qmat(1,9) = Qmat(1,1)
         Qmat(2,9) = Qmat(2,1)
         Qmat(3,9) = Qmat(3,1)
         Qmat(4,9) = Qmat(4,1)
         Qmat(5,9) = Qmat(5,1)
         Qmat(6,9) = Qmat(6,1)
         Qmat(7,9) = Qmat(7,1)
         Qmat(8,9) = Qmat(8,1)
         Qmat(9,9) = Qmat(9,1)
         

         if((nNode.eq.8).and.(nInt.eq.8)) then
            !
            ! This is the tangent using the F-bar method with the
            !  8 node fully integrated linear element
            !
            Kuu = Kuu + detMapJC*w(intpt)*
     +           (
     +           matmul(matmul(transpose(Gmat),Amat),Gmat)
     +           + matmul(transpose(Gmat),matmul(Qmat,(G0mat-Gmat)))
     +           )
         else
            !
            ! This is the tangent NOT using the F-bar method with all
            !  other elements
            !
            Kuu = Kuu + detMapJC*w(intpt)*
     +           (
     +           matmul(matmul(transpose(Gmat),Amat),Gmat)
     +           )
         endif


      enddo
      !
      ! End the loop over integration points
      !----------------------------------------------------------------

      ! Return the RHS vector and the Stiffness matrix.
      !
      rhs(1:24,1) = Ru(:,1)
      amatrx = Kuu
      !
      ! End return of RHS and AMATRX
      !----------------------------------------------------------------

      return
      end subroutine U3D8

************************************************************************

      subroutine integ(props,nprops,dtime,
     +        F_tau,magFieldVec,
     +        T_tau,
     +        SpUUMod,stat)


      implicit none

      integer i,j,k,l,m,n,stat,nprops

      real*8 Iden(3,3),F_tau(3,3),FT(3,3),FinvT(3,3),C_tau(3,3),detC
      real*8 Cinv(3,3),trC,I1bar,detF,Finv(3,3),magFieldVec(3,1)
      real*8 TR_tau(3,3),T_tau(3,3),dTRdF(3,3,3,3),Gshear,Kbulk,G0,trAC
      real*8 SpUUmod(3,3,3,3),lamL,aux
      real*8 props(nprops),dtime,AC(3,3),lamBar
      real*8 Fbar(3,3),AFt(3,3),trFAFt,FAFt(3,3)
      real*8 FA(3,3),dGdF(3,3),refMagMoment(3,1),TR_eq(3,3),TR_mag(3,3)
      real*8 MB(3,3),MBFinvT(3,3),FMBFinvT(3,3)


      real*8 zero,one,two,three,third,half,nine
      parameter(zero=0.d0,one=1.d0,two=2.d0,three=3.d0,third=1.d0/3.d0,
     +     half=1.d0/2.d0,nine=9.d0)


      ! Identity matrix
      !
      call onem(Iden)

      ! Obtain material parameters
      !
      Gshear = props(1)
      Kbulk  = props(2)
      refMagMoment(1,1) = props(3)
      refMagMoment(2,1) = props(4)
      refMagMoment(3,1) = props(5)

      ! Compute the determinant, the inverse, the transpose, 
      !  and the inverse transpose of the deformation gradient
      !
      call matInv3D(F_tau,Finv,detF,stat)
      FT = transpose(F_tau)
      FinvT = transpose(Finv)
      Fbar = (detF**(-two/three))*F_tau


      ! Compute the right Cauchy-Green tensor, its inverse, 
      !  its determinant, and its trace
      !
      C_tau = matmul(transpose(F_tau),F_tau)
      call matInv3D(C_tau,Cinv,detC,stat)
      trC = C_tau(1,1) + C_tau(2,2) + C_tau(3,3)
 
      ! Compute the equilibrium 1st Piola stress, including the 
      !  magnetic contribution
      !
      TR_tau = (detF**(-two/three))*Gshear*(F_tau - third*trC*FinvT)
     +     + Kbulk*detF*(detF - one)*FinvT
     +     -matmul(magFieldVec,transpose(refMagMoment))

      ! Compute the total Cauchy stress
      !
      T_tau = (one/detF)*matmul(TR_tau,transpose(F_tau))

      ! Calculate the material tangent modulus
      !
      dTRdF = zero
      do i=1,3
         do j=1,3
            do k=1,3
               do l=1,3
                  dTRdF(i,j,k,l) = dTRdF(i,j,k,l)
     +                 + (detF**(-two/three))*Gshear*
     +                 (
     +                 (-two/three)*F_tau(i,j)*Finv(l,k)
     +                 + (two/9.d0)*trC*Finv(j,i)*Finv(l,k)
     +                 + Iden(i,k)*Iden(j,l)
     +                 + third*trC*Finv(l,i)*Finv(j,k)
     +                 - (two/three)*Finv(j,i)*F_tau(k,l)
     +                 )
     +                 + detF*Kbulk*
     +                 (
     +                 (detF-one)*Finv(j,i)*Finv(l,k)
     +                 + detF*Finv(j,i)*Finv(l,k)
     +                 - (detF-one)*Finv(l,i)*Finv(j,k)
     +                 )
               enddo
            enddo
         enddo
      enddo


      ! Calculate the material spatial tangent modulus
      !
      SpUUMod = zero
      do i=1,3
         do j=1,3
            do k=1,3
               do l=1,3
                  do m=1,3
                     do n=1,3
                        SpUUMod(i,j,k,l) = SpUUMod(i,j,k,l) +
     +                      (dTRdF(i,m,k,n)*F_tau(j,m)*F_tau(l,n))/detF
                     enddo
                  enddo
               enddo
            enddo
         enddo
      enddo


      return
      end subroutine integ

!****************************************************************************
!     Element subroutines
!****************************************************************************

      subroutine xint3D1pt(xi,w,nIntPt)
      
      ! This subroutine will get the integration point locations
      !  and corresponding gauss quadrature weights for 3D elements
      !  using a 2 gauss points for integration
      !
      !  xi(nIntPt,3): xi,eta,zeta coordinates for the integration pts
      !  w(nIntPt):    corresponding integration weights
      
      implicit none

      integer nIntPt,nDim

      real*8 xi(1,3),w(1)


      ! Init
      !
      w = 0.d0
      xi = 0.d0


      ! Number of Gauss points
      !
      nIntPt = 1


      ! Gauss weights
      !
      w(1) = 8.d0
      

      ! Gauss pt locations in master element
      !
      xi(1,1) = 0.d0
      xi(1,2) = 0.d0
      xi(1,3) = 0.d0

      return
      end subroutine xint3D1pt
     
************************************************************************

      subroutine xint3D8pt(xi,w,nIntPt)
      
      ! This subroutine will get the integration point locations
      !  and corresponding gauss quadrature weights for 3D elements
      !  using 8 gauss points for integration
      !
      !  xi(nIntPt,3): xi,eta,zeta coordinates for the integration pts
      !  w(nIntPt):    corresponding integration weights
      
      implicit none

      integer nIntPt,nDim

      real*8 xi(8,3),w(8)


      ! Init
      !
      w = 0.d0
      xi = 0.d0


      ! Number of Gauss points
      !
      nIntPt = 8


      ! Gauss weights
      !
      w(1) = 1.d0
      w(2) = 1.d0
      w(3) = 1.d0
      w(4) = 1.d0
      w(5) = 1.d0
      w(6) = 1.d0
      w(7) = 1.d0
      w(8) = 1.d0
      

      ! Gauss pt locations in master element
      !
      xi(1,1) = -dsqrt(1.d0/3.d0)
      xi(1,2) = -dsqrt(1.d0/3.d0)
      xi(1,3) = -dsqrt(1.d0/3.d0)
      xi(2,1) = dsqrt(1.d0/3.d0)
      xi(2,2) = -dsqrt(1.d0/3.d0)
      xi(2,3) = -dsqrt(1.d0/3.d0)
      xi(3,1) = -dsqrt(1.d0/3.d0)
      xi(3,2) = dsqrt(1.d0/3.d0)
      xi(3,3) = -dsqrt(1.d0/3.d0)
      xi(4,1) = dsqrt(1.d0/3.d0)
      xi(4,2) = dsqrt(1.d0/3.d0)
      xi(4,3) = -dsqrt(1.d0/3.d0)
      xi(5,1) = -dsqrt(1.d0/3.d0)
      xi(5,2) = -dsqrt(1.d0/3.d0)
      xi(5,3) = dsqrt(1.d0/3.d0)
      xi(6,1) = dsqrt(1.d0/3.d0)
      xi(6,2) = -dsqrt(1.d0/3.d0)
      xi(6,3) = dsqrt(1.d0/3.d0)
      xi(7,1) = -dsqrt(1.d0/3.d0)
      xi(7,2) = dsqrt(1.d0/3.d0)
      xi(7,3) = dsqrt(1.d0/3.d0)
      xi(8,1) = dsqrt(1.d0/3.d0)
      xi(8,2) = dsqrt(1.d0/3.d0)
      xi(8,3) = dsqrt(1.d0/3.d0)


      return
      end subroutine xint3D8pt

************************************************************************
************************************************************************

      subroutine calcShape3DLinear(nIntPt,xi_int,intpt,sh,dshxi)
      !
      !
      ! Calculate the shape functions and their derivatives at the
      ! given integration point in the master element
      !
      ! This subroutine uses a 8-node linear 3D element as shown
      !
      !      8-----------7
      !     /|          /|       zeta
      !    / |         / |       
      !   5-----------6  |       |     eta
      !   |  |        |  |       |   /
      !   |  |        |  |       |  /
      !   |  4--------|--3       | /
      !   | /         | /        |/
      !   |/          |/         O--------- xi
      !   1-----------2        origin at cube center
      !
      !
      ! sh(i) = shape function of node i at the intpt.
      ! dshxi(i,j) = derivative wrt j direction of shape fn of node i
      ! d2shxi(i,j,k) = derivatives wrt j and k of shape fn of node i

      implicit none

      integer intpt,nDim,nIntPt,i,j

      real*8 xi_int(nIntPt,3),sh(8),dshxi(8,3)
      real*8 d2shxi(8,3,3),xi,eta,zeta

      real*8 zero,one,two,half,fourth,eighth
      parameter(zero=0.d0,one=1.d0,two=2.d0,half=0.5d0,fourth=0.25d0,
     +     eighth=1.d0/8.d0)
      

      ! Location in the master element
      !
      xi = xi_int(intpt,1)
      eta = xi_int(intpt,2)
      zeta = xi_int(intpt,3)
      !
      ! The shape functions
      !
      sh(1) = eighth*(one - xi)*(one - eta)*(one - zeta)
      sh(2) = eighth*(one + xi)*(one - eta)*(one - zeta)
      sh(3) = eighth*(one + xi)*(one + eta)*(one - zeta)
      sh(4) = eighth*(one - xi)*(one + eta)*(one - zeta)
      sh(5) = eighth*(one - xi)*(one - eta)*(one + zeta)
      sh(6) = eighth*(one + xi)*(one - eta)*(one + zeta)
      sh(7) = eighth*(one + xi)*(one + eta)*(one + zeta)
      sh(8) = eighth*(one - xi)*(one + eta)*(one + zeta)
      !
      ! The first derivatives
      !
      dshxi(1,1) = -eighth*(one - eta)*(one - zeta)
      dshxi(1,2) = -eighth*(one - xi)*(one - zeta)
      dshxi(1,3) = -eighth*(one - xi)*(one - eta)
      dshxi(2,1) = eighth*(one - eta)*(one - zeta)
      dshxi(2,2) = -eighth*(one + xi)*(one - zeta)
      dshxi(2,3) = -eighth*(one + xi)*(one - eta)
      dshxi(3,1) = eighth*(one + eta)*(one - zeta)
      dshxi(3,2) = eighth*(one + xi)*(one - zeta)
      dshxi(3,3) = -eighth*(one + xi)*(one + eta)
      dshxi(4,1) = -eighth*(one + eta)*(one - zeta)
      dshxi(4,2) = eighth*(one - xi)*(one - zeta)
      dshxi(4,3) = -eighth*(one - xi)*(one + eta)
      dshxi(5,1) = -eighth*(one - eta)*(one + zeta)
      dshxi(5,2) = -eighth*(one - xi)*(one + zeta)
      dshxi(5,3) = eighth*(one - xi)*(one - eta)
      dshxi(6,1) = eighth*(one - eta)*(one + zeta)
      dshxi(6,2) = -eighth*(one + xi)*(one + zeta)
      dshxi(6,3) = eighth*(one + xi)*(one - eta)
      dshxi(7,1) = eighth*(one + eta)*(one + zeta)
      dshxi(7,2) = eighth*(one + xi)*(one + zeta)
      dshxi(7,3) = eighth*(one + xi)*(one + eta)
      dshxi(8,1) = -eighth*(one + eta)*(one + zeta)
      dshxi(8,2) = eighth*(one - xi)*(one + zeta)
      dshxi(8,3) = eighth*(one - xi)*(one + eta)
      !
      ! The second derivatives
      !
      d2shxi = zero
      d2shxi(1,1,2) = eighth*(one - zeta)
      d2shxi(1,2,1) = d2shxi(1,1,2)
      d2shxi(1,1,3) = eighth*(one - eta)
      d2shxi(1,3,1) = d2shxi(1,1,3)
      d2shxi(1,2,3) = eighth*(one - xi)
      d2shxi(1,3,2) = d2shxi(1,2,3)
      d2shxi(2,1,2) = -eighth*(one - zeta)
      d2shxi(2,2,1) = d2shxi(2,1,2)
      d2shxi(2,1,3) = -eighth*(one - eta)
      d2shxi(2,3,1) = d2shxi(2,1,3)
      d2shxi(2,2,3) = eighth*(one + xi)
      d2shxi(2,3,2) = d2shxi(2,2,3)
      d2shxi(3,1,2) = eighth*(one - zeta)
      d2shxi(3,2,1) = d2shxi(2,1,2)
      d2shxi(3,1,3) = -eighth*(one + eta)
      d2shxi(3,3,1) = d2shxi(2,1,3)
      d2shxi(3,2,3) = -eighth*(one + xi)
      d2shxi(3,3,2) = d2shxi(2,2,3)
      d2shxi(4,1,2) = -eighth*(one - zeta)
      d2shxi(4,2,1) = d2shxi(2,1,2)
      d2shxi(4,1,3) = eighth*(one + eta)
      d2shxi(4,3,1) = d2shxi(2,1,3)
      d2shxi(4,2,3) = -eighth*(one - xi)
      d2shxi(4,3,2) = d2shxi(2,2,3)
      d2shxi(5,1,2) = eighth*(one + zeta)
      d2shxi(5,2,1) = d2shxi(2,1,2)
      d2shxi(5,1,3) = -eighth*(one - eta)
      d2shxi(5,3,1) = d2shxi(2,1,3)
      d2shxi(5,2,3) = -eighth*(one - xi)
      d2shxi(5,3,2) = d2shxi(2,2,3)
      d2shxi(6,1,2) = eighth*(one + zeta)
      d2shxi(6,2,1) = d2shxi(2,1,2)
      d2shxi(6,1,3) = eighth*(one - eta)
      d2shxi(6,3,1) = d2shxi(2,1,3)
      d2shxi(6,2,3) = -eighth*(one + xi)
      d2shxi(6,3,2) = d2shxi(2,2,3)
      d2shxi(7,1,2) = eighth*(one + zeta)
      d2shxi(7,2,1) = d2shxi(2,1,2)
      d2shxi(7,1,3) = eighth*(one + eta)
      d2shxi(7,3,1) = d2shxi(2,1,3)
      d2shxi(7,2,3) = eighth*(one + xi)
      d2shxi(7,3,2) = d2shxi(2,2,3)
      d2shxi(8,1,2) = -eighth*(one + zeta)
      d2shxi(8,2,1) = d2shxi(2,1,2)
      d2shxi(8,1,3) = -eighth*(one + eta)
      d2shxi(8,3,1) = d2shxi(2,1,3)
      d2shxi(8,2,3) = eighth*(one - xi)
      d2shxi(8,3,2) = d2shxi(2,2,3)
      
      return
      end subroutine calcShape3DLinear
!************************************************************************

      subroutine mapShape3D(nNode,dshxi,coords,dsh,detMapJ,stat)
      !
      ! Map derivatives of shape fns from xi-eta-zeta domain
      !  to x-y-z domain.  This subroutine works for both 8-node
      !  linear and 20-node quadratic 3D elements.
      !
      implicit none

      integer i,j,k,nNode,ieror,stat

      real*8 dshxi(nNode,3),dsh(nNode,3),coords(3,nNode)
      real*8 mapJ(3,3),mapJ_inv(3,3),detmapJ

      real*8 zero,one,two,half,fourth,eighth
      parameter(zero=0.d0,one=1.d0,two=2.d0,half=0.5d0,fourth=0.25d0,
     +     eighth=1.d0/8.d0)
      

      ! Calculate the mapping Jacobian matrix:
      !
      mapJ = zero
      do i=1,3
        do j=1,3
          do k=1,nNode
              mapJ(i,j) = mapJ(i,j) + dshxi(k,i)*coords(j,k)
          end do
        end do
      end do


      ! Calculate the inverse and the determinant of Jacobian
      !
      call matInv3D(mapJ,mapJ_inv,detMapJ,stat)


      ! Calculate first derivatives wrt x, y, z
      !
      dsh = transpose(matmul(mapJ_inv,transpose(dshxi)))


      ! The second derivatives may be calculated.
      !

      return
      end subroutine mapShape3D

!****************************************************************************
!     Utility subroutines
!****************************************************************************

      subroutine matInv3D(A,A_inv,det_A,istat)
      !
      ! Returns A_inv, the inverse and det_A, the determinant
      ! Note that the det is of the original matrix, not the
      ! inverse
      !
      implicit none
      !
      integer istat
      !
      real*8 A(3,3),A_inv(3,3),det_A,det_A_inv


      istat = 1
      
      det_A = A(1,1)*(A(2,2)*A(3,3) - A(3,2)*A(2,3)) -
     +        A(2,1)*(A(1,2)*A(3,3) - A(3,2)*A(1,3)) +
     +        A(3,1)*(A(1,2)*A(2,3) - A(2,2)*A(1,3))
      
      if (det_A .le. 0.d0) then
        write(*,*) 'WARNING: subroutine matInv3D:'
        write(*,*) 'WARNING: det of mat=',det_A
        istat = 0
        return
      end if
          
      det_A_inv = 1.d0/det_A
        
      A_inv(1,1) = det_A_inv*(A(2,2)*A(3,3)-A(3,2)*A(2,3))
      A_inv(1,2) = det_A_inv*(A(3,2)*A(1,3)-A(1,2)*A(3,3))
      A_inv(1,3) = det_A_inv*(A(1,2)*A(2,3)-A(2,2)*A(1,3))
      A_inv(2,1) = det_A_inv*(A(3,1)*A(2,3)-A(2,1)*A(3,3))
      A_inv(2,2) = det_A_inv*(A(1,1)*A(3,3)-A(3,1)*A(1,3))
      A_inv(2,3) = det_A_inv*(A(2,1)*A(1,3)-A(1,1)*A(2,3))
      A_inv(3,1) = det_A_inv*(A(2,1)*A(3,2)-A(3,1)*A(2,2))
      A_inv(3,2) = det_A_inv*(A(3,1)*A(1,2)-A(1,1)*A(3,2))
      A_inv(3,3) = det_A_inv*(A(1,1)*A(2,2)-A(2,1)*A(1,2))
      

      return
      end subroutine matInv3D

!****************************************************************************

      subroutine mdet(A,det)
      !
      ! This subroutine calculates the determinant
      ! of a 3 by 3 matrix [A]
      !
      implicit none
      !
      real*8  A(3,3),det


      det = A(1,1)*A(2,2)*A(3,3) 
     +	  + A(1,2)*A(2,3)*A(3,1)
     +	  + A(1,3)*A(2,1)*A(3,2)
     +	  - A(3,1)*A(2,2)*A(1,3)
     +	  - A(3,2)*A(2,3)*A(1,1)
     +	  - A(3,3)*A(2,1)*A(1,2)


      return
      end subroutine mdet
	
!****************************************************************************

      subroutine onem(A)
      !
      ! This subroutine stores the identity matrix in the
      ! 3 by 3 matrix [A]
      !
      implicit none
      !
      integer i,j
      !
      real*8 A(3,3)


      do i=1,3
         do J=1,3
	    if (i .eq. j) then
              A(i,j) = 1.0
            else
              A(i,j) = 0.0
            end if
         end do
      end do


      return
      end subroutine onem
!
!*****************************************************************************
       subroutine GetgradB(centroidcoords,magFieldMag,magFieldVec,gradB,n)
	  !gradB of n cloumn, e.g. n=1, gradB(1,1),gradB(2,1),gradB(3,1), i.e. partial Bx partial x,partial By partial x, partial Bz partial x
	  !gradB of n cloumn, e.g. n=2, gradB(1,2),gradB(2,2),gradB(2,3), i.e. partial Bx partial y,partial By partial y, partial Bz partial y
	  !gradB of n cloumn, e.g. n=3, gradB(1,3),gradB(2,3),gradB(3,3), i.e. partial Bx partial z,partial By partial z, partial Bz partial z
	  implicit none
	  REAL(8):: gradB,IncCentroid,IncMag,inc,centroidcoords,magFieldVec,magFieldMag
	  integer n
      DIMENSION gradB(3,1),IncCentroid(3),IncMag(3,1),centroidcoords(3),magFieldVec(3,1)
	  
	  !numerically obtain the gradient by setting infinitesmal increment inc=0.00001 (should be smaller than the element size)
	  inc=0.0010/1000
	  IncCentroid=centroidcoords
	  
	  !only change the centroid coordinate in n-direction, n=1=x, n=2=y, n=3=z directions
	  IncCentroid(n)=centroidcoords(n)+inc
	  
	  call GetMagField(IncCentroid,magFieldMag,IncMag)
	  
	  ! change from 1D row to 1D column
	
      gradB(1,1)=(IncMag(1,1)-magFieldVec(1,1))/inc 
      gradB(2,1)=(IncMag(2,1)-magFieldVec(2,1))/inc 
	  gradB(3,1)=(IncMag(3,1)-magFieldVec(3,1))/inc 
	  !write(*,*), 'check magnetic field'
	  !write(*,*),  magFieldVec
	  !write(*,*), 'check inc magnetic field'
	  !write(*,*),  IncMag
   	  end subroutine GetgradB  
	  
	  
!***********************************************************************	  
	  subroutine GetMagField(centroidcoords,magFieldMag,magFieldVec)
	  
	  implicit none
	  REAL(8):: centroidcoords,magFieldVec,shiftPosition,globalcentroid,a,b,c,L,x,y,z,gap,tempB, ansx,ansy,ansz,coef,magFieldMag
      DIMENSION centroidcoords(3),globalcentroid(3),magFieldVec(3,1),shiftPosition(3), ansx(8),ansy(4),ansz(4)
	  !define magnet and bar size, here dimesions are normalized by L (L=1 is the size of abauqs model)
	  L=60.0/1000
	  !unit [m]
	  a=L
	  b=L
	  c=L
	  gap=20.0/1000
	
	  
	  !magFieldMag is the flux density on the surface of the rigid magnet  B0
	  !Calculate the mu*M/pi according to Eq. (11) in reference "Alternative method to calculate the magnetic ﬁeld of permanent magnets with azimuthal symmetry"
	  coef=3.141592654/2-atan(c*b/(2*a*sqrt(c**2+b**2+(2*a)**2)))
	  tempB=magFieldMag/coef

      !ShiftPosition defines the position vector from the center of the magnet to CAE model (0,0,0) 
	  shiftPosition(1)=L+gap
      shiftPosition(2)=0.0
	  shiftPosition(3)=0.0

	  globalcentroid=shiftPosition+centroidcoords
	 
	  x=globalcentroid(1)
	  y=globalcentroid(2)
	  z=globalcentroid(3)
	  
	  !if set the z-symmetric plane at z=0.
	  !write(*,*),'check centroidcoords'
	  !write(*,*), centroidcoords
	  !write(*,*),'check golcentroid'
	  !write(*,*), globalcentroid
	  !Please define Bz conponent here
	
	  call function2(a,b,c,-x,-z,y,ansz(1))
	  call function2(a,b,c,x,z,y,ansz(2))
	  call function2(a,b,c,-x,z,y,ansz(3))
	  call function2(a,b,c,x,-z,y,ansz(4))
	  magFieldVec(3,1)=tempB/4*log(ansz(1)*ansz(2)/(ansz(3)*ansz(4)))
	  
	!Please define By conponent here  
	  call function2(a,b,c,-x,-y,z,ansy(1))
	  call function2(a,b,c,x,y,z,ansy(2))
	  call function2(a,b,c,-x,y,z,ansy(3))
	  call function2(a,b,c,x,-y,z,ansy(4))
	  magFieldVec(2,1)=tempB/4*log(ansy(1)*ansy(2)/(ansy(3)*ansy(4)))
	  
	!Please define Bx conponent here	
	  call function1(a,b,c,-x,y,z,ansx(1))
	  call function1(a,b,c,-x,y,-z,ansx(2))
	  call function1(a,b,c,-x,-y,z,ansx(3))
	  call function1(a,b,c,-x,-y,-z,ansx(4))
	  call function1(a,b,c,x,y,z,ansx(5))
	  call function1(a,b,c,x,y,-z,ansx(6))
	  call function1(a,b,c,x,-y,z,ansx(7))
	  call function1(a,b,c,x,-y,-z,ansx(8))
	  
	 ! dont forget the negtive sign here
	 magFieldVec(1,1)=-tempB/4*sum(ansx)
	 !
	 !write(*,*), 'check magnetic field'
	 !write(*,*),  magFieldVec
	 
	  end subroutine GetMagField
!******************************************************************************	  	
	 
	 subroutine function1(a,b,c,x,y,z,ans)
	 implicit none
	 REAL(8):: a,b,c,x,y,z,ans
	 ans=atan ((z+c)*(y+b)/((x+a)*sqrt((x+a)**2+(y+b)**2+(z+c)**2)))
	 end subroutine function1
	 
	 subroutine function2(a,b,c,x,y,z,ans)
	 implicit none
	 REAL(8):: a,b,c,x,y,z,ans
	 ans=(sqrt((x+a)**2+(z-c)**2+(y+b)**2)+c-z)/(sqrt((x+a)**2+(y+b)**2+(z+c)**2)-c-z)
	 
	 end subroutine function2
!*******************************************************************************
	  

	  
	  
	  
	  
	  