C The following routine allows the user to specify flux on an element edge
C based on contact pressure in the previous increment.
C A factor needs to be specified to convert the pressure to flux
C *** Works only in deformable-deformable or deformable-rigid geometry contact.
C *** Does not work if the body on which the flux to be applied is a rigid heat transfer body
C         (The rigid body is meshed for the purpose of heat transfer)
c
c------------------------------------------------------
      subroutine flux(f2,temflu,mibody,time)
C      subroutine forcem(press,th1,th2,nn,n)
      implicit real*8 (a-h,o-z)                                               dp
      dimension mibody(7)
c* * * * * *
c
c     user subroutine for non-uniform flux input.
c
c     f2           flux value
c     temflu(1)    estimated temperature
c     temflu(2)    previous volumetric flux
c     temflu(3)    temperature at beginning of increment
c     temflu(4,5,6)integration point coordinates
c     mibody(1)    element number
c     mibody(2)    flux type
c     mibody(3)    integration point number
c     time         time
c
c* * * * * *
      include '../common/form'
      include '../common/array4'
      include '../common/concom'
      include '../common/space'
      include '../common/blnk'

      include '../common/dimen'
      include '../common/array5'
      include '../common/peturb'
      parameter (MAX_ELS=1000, NODES_PER_EL=4,ND_PER_EDGE=2)
      dimension nod(NODES_PER_EL)
      dimension ni(ND_PER_EDGE)
      equivalence (nod(1),lm(1))

C The subroutine is limited to 1000 contact elements on which flux is applied.
C This may be increased per user requirements
	dimension ncont(1000),cf(1000,2),df(1000,2),ncn(4)
	common/mycont/nc,ncont,cf,df

c-----------------------------------------------------

c     The following codes may change from one element type to another
C  Please check Element library documentation for appropriate codes
c

        if (mibody(2) .eq. 3) then
          ni(1) = 1
          ni(2) = 2
        else if (mibody(2) .eq. 6) then
          ni(1) = 2
          ni(2) = 3
        else if (mibody(2) .eq. 9) then
          ni(1) = 3
          ni(2) = 4
        else if (mibody(2) .eq. 11) then
          ni(1) = 4
          ni(2) = 1
        endif


      ichk = 0

C ----------------------------- Begin contact check
C Find how many nodes are in contact 
C if 2 nodes are in contact then the edge is in contact (ichk=2)
      do 80 k=1,ND_PER_EDGE
	nnod = nod(ni(k))
	do 70 i=1,nc
	  if(nnod.eq.ncont(i)) then
	    ichk = ichk+1
	    ncn(k) = i
       write(65,*)"inc,El,Nd,ichk,nc",inc,mibody(1),nod(ni(k)),ichk,nc
	    go to 80
	  endif
 70     continue
 80   continue

c     End of contact check
c------------------------------------------------------
c
C If ichk=2, 2 nodes of the edge are in contact (So apply flux)

      if(ichk.eq.2) then
         factor = 1.0
C         write(66,1001) inc,mibody(1),nn,factor,mibody(2),ichk
C Compute the edge length
 	dfx = df(ncn(2),1)-df(ncn(1),1)
	dfy = df(ncn(2),2)-df(ncn(1),2)
	edgelen = sqrt(dfx*dfx + dfy*dfy)

C Compute the contact pressure
	cfx = 0.
	cfy = 0.
	do 90 i=1,ND_PER_EDGE
	  cfx = cfx + cf(ncn(i),1)
	  cfy = cfy + cf(ncn(i),2)
C	  write(57,*)"Load ",inc,nod(ni(i)),cf(ncn(i),1),cf(ncn(i),2)
 90	continue
          avgcf = sqrt((0.5*cfx)**2+(0.5*cfy)**2)
        cntpres = avgcf/edgelen

C	write(57,*) "Tot load: ",mibody(1),cfx,cfy,avgcf,edgelen
C	write(57,*) "Cont pres ",inc,mibody(1),avgcf,cntpres,factor
      else 
	 factor = 0.
      end if
c
C      f2 = factor * f2
C **** User will have to specify the factor used to convert the contact pressure to flux
	f2 = cntpres * factor
      return
      end
c
c

C ===============================================================
C Routine to extract internal/external node corresponding to
C an external/internal node.

	subroutine ndinex(lint,lext,it)
	implicit real*8 (a-h,o-z)

C lint	  internal node number
C lext	  external node number
C it = 1  convert internal node number to external (user) node
C    = 2  convert external node to internal node

	include '../common/arrays'
	include '../common/develp'
	include '../common/space'
	include '../common/dimen'
	include '../common/prepro'

	lint1=lint
	lext1=lext
	if(it.eq.1) then
C convert the internal node number to user node number
	  lext1=lint1
	  la2=inpnum+lint1-1
	  if(joptit.ne.0) lext1=igetsh(ints(la2),1)
	  if(nnoids.ne.0) lext1=ints(inoids+lext1-1)
	  lext=lext1
	else if(it.eq.2) then
C convert the external node to internal node number
	  lint1=ibsrch(lext1,ints(inoids),numnp,1)
	  la2=inpnum+lint1-1
	  if(joptit.ne.0)lint1=igetsh(ints(la2),0)
	  lint=lint1
	endif
	
	return
	end

C ===============================================================
C This subroutine returns the coordinates/displacements of a node.
	subroutine ncoord(lint,ncflag,cdnode)
	implicit real*8 (a-h,o-z)

C   lint	internal node number
C   ncflag = 0  Original coordinates of a node
C          = 1  Displacements of a node
C          = 2  Deformed coordinates of a node
C  cdnode 	Returned array of coordinates/displacements based on ncflag

	include '../common/dimen'
	include '../common/array2'
	include '../common/space'
	dimension cdnode(12),disp(12)
	jrdpre=0

C Original coordinates
	if (ncflag .eq. 0) then
	  call vecftc(cdnode,vars(ixord),ncrdmx,ncrd,lint,jrdpre,2,1)
C Displacements
	else if (ncflag .eq. 1) then
	  call vecftc(cdnode,vars(idsxt),ndegmx,ncrd,lint,jrdpre,2,5)
C Current coordinates
	else if (ncflag .eq. 2) then
	  call vecftc(cdnode,vars(ixord),ncrdmx,ncrd,lint,jrdpre,2,1)
	  call vecftc(disp,vars(idsxt),ndegmx,ncrd,lint,jrdpre,2,5)
	  do 10 i=1,ncrdmx
	    cdnode(i) = cdnode(i) + disp(i)
 10	  continue
	endif
	
	return
	end



C ===============================================================
C The following routine checks for the nodes in contact and puts together
C arrays of displacements and contact forces.

      subroutine uedinc(inc1,incsub1)
      implicit real*8 (a-h,o-z)
c     dummy user subroutine that would get called at the end of each
c     increment

c
c     for each node 
c        icntyp = 0  no contact
c               = 1  node is contacting rigid surface only
c               = 2  node is a tied node
c               = 3  node is a retained node
c               >10  node is a retained node in multiple ties
c
c	contf	contact force components
c
C ndiflx	No. of flexible bodies
C nbct		No. of nodes on the boundary of contact body
C nbcn		Upper bound to the no. of nodes on the boundary of contact body
C setwr2	Internal routine, copies contact information to iwrkn2
C vecftc	Internal routine, used to copy nodal information from global vector
C inod		internal node number (that is in contact)
C contf		Array containing contact forces for a given node
C i2or3		2d or 3D problem
C lext		external (user) node number
C cf		contact force array
C df		deformed coordinate array
C nc		Total number of nodes in contact

      include '../common/array4'
      include '../common/array5'
      include '../common/arrays'
      include '../common/develp'
      include '../common/space'
      include '../common/peturb'
      include '../common/dimen'
      include '../common/form'
      include '../common/form1'


	dimension contf(10),ncont(1000),cf(1000,2),df(1000,2)
	dimension cdnode(3)
	save incold

	common/mycont/nc,ncont,cf,df
	
	if(incold.ne.inc1) then
	  incold=inc1
	endif
	
	nc = 0
c
	do 100 md=1,ndiflx
	 nbct=ints(inbct+md-1)
	 do 200 i=1,nbct
	  inod=ints(inf+(md-1)*nbcn+i-1)
	  call vecftc(contf,vars(ipload),ndegmx,ndeg,inod,jrdpre,0,15)

C convert the internal node number to user node number
C  Need to modify for 3D analysis
C	  call ndinex(inod,lext,1)
          call setwr2(ints(iwrkn2))
	  icntyp = ints(iwrkn2+inod-1)
	  if(icntyp.gt.0) then
	   if(i2or3.eq.2) then
	    nc = nc + 1
	    ncont(nc) = inod
	    cf(nc,1) = contf(1)
	    cf(nc,2) = contf(2)
	    call ncoord(inod,2,cdnode)
	    df(nc,1) = cdnode(1)
	    df(nc,2) = cdnode(2)
111	    format(3i7,5X,2e14.5)
	   else if(i2or3.eq.3) then
112	    format(3i7,5X,3e14.5)
	   endif
	  endif


200	continue
100	continue
c
      return
      end