! This Subroutine Implements the Level Set Method
	! J. Grogan - 07/10/13
	subroutine uexternaldb(lop,lrestart,time,dtime,kstep,kinc)
		include 'aba_param.inc'
		dimension 	time(2)
		integer	  	Element (100,4)
		real		Node(100,2),LSet(100)
	!	
		if(lop==0)then
	!		Get Mesh Data
			call getMesh(Element,Node,numElem,numNodes)
	!		Get Initial Level Set
			call initialLSet(LSet,Node,numNodes)
		elseif(lop==1)then
	!		Update Level Set
			call updateLSet(Node,numNodes,Element,numElem,dtime,LSet)
		endif	
    return
    end	
	! This subroutine returns the finite element mesh connectivity data
	subroutine getMesh(Element,Node,numElem,numNodes)
		include 'aba_param.inc'
		integer	  	Element (100,4)
		real		Node(100,2)
		character(256) outdir,jobname,input
		call getoutdir(outdir,lenoutdir)
		call getjobname(jobname,lenjobname)
		filename=trim(outdir)//trim(jobname)//'.inp'
		open(unit=107,file=filename,status='old')
		read(107,*)input
		do while (index(input,'*Node')==0)
			read(107,*)input
		enddo
		ierr=0
		numNodes=0
		do while (ierr==0)
			read(107,*)nodeNum,xcor,ycor,zcor
			if(ierr==0)then
				Node(nodeNum,1)=xcor
				Node(nodeNum,2)=ycor
				numNodes=numNodes+1
			endif
		enddo
		do while (index(input,'*Element')==0)
			read(107,*)input
		enddo	
		numElem=0
		do while (ierr==0)
			read(107,*)elNum,n1,n2,n3,n4
			if(ierr==0)then
				Element(elNum,1)=n1
				Element(elNum,2)=n2
				Element(elNum,3)=n3
				Element(elNum,4)=n4
				numElem=numElem+1
			endif
		enddo	
		close(107)
	end subroutine
	!	This subroutine calculates the initial Level Set
	subroutine initialLSet(LSet,Node,numNodes)
		include 'aba_param.inc'
		real		Node(100,2),LSet(100)	
		!centx=4.
		!centy=4.
		!rad=2.1
		!do i=1,numNodes
		!    dist=sqrt((Node(i,1)-centx)*(Node(i,1)-centx)+(Node(i,2)-centy)*(Node(i,2)-centy))
		!    LSet(i)=dist-rad
		!enddo
		do i=1,numNodes
		   dist=Node(i,1)-0.1
		   LSet(i)=dist
		enddo	
	end subroutine
	! This subroutine updates the Level Set
	subroutine updateLSet(Node,numNodes,Element,numElem,dtime,LSet)
		include 'aba_param.inc'
		integer	  	Element(100,4),NBelem(100,4)
		real		Node(100,2),NGlobal(100,2),NLocal(100,2)
		real		LSet(100),LSetLocal(100),F(100),Fred(100)
		real		A(100,100),Ared(100,100),RHS(100),RHSred(100)
		real 		M(100,100),MGL(100,100),f1(100),f2(100),f3(100)
		! parameters
		bandWidth=10.
		av_Length=1.
		h2=0.00001*av_Length
		visc=0.0005
		dt=0.01
		! explicit update of Level Set
		do istep=1:floor(dtime/dt)
			! Identify Narrow Band Elements and Get Local Level Set
			call getNarrowBand(NBelem,NBElems,NGlobal,NLocal,NBNodes,LSetLocal, &
				& bandWidth,LSet,Element,numElem,numNodes)				
			! Identify Scalar Velocity on Nodes Crossed By Interface - F
			call getF(F,LSetLocal,NBElems,NBNodes,NBelem)
			! Get 'Stiffness' Matrix - A
			call getA(A,Node,NLocal,NBelem,NBNodes,NBElems,LSetLocal)
			! Apply BCs
			RHS=-matmul(A,F)
			iindex=0
			do i=1,NBNodes
				if (F(i)==0.)then
					iindex=iindex+1
					RHSred(iindex)=RHS(i)
					jindex=0
					do j=1,NBNodes
						if (F(j)==0.)then
							jindex=jindex+1;
							Ared(iindex,jindex)=A(i,j)
						endif
					enddo
				endif
			enddo
			! Solve for Fred
			Fred=(Ared^-1)*RHSred'
			! Get F
			iindex=0
			do i=1,NBNodes
				if (F(i)==0.)then
					iindex=iindex+1
					F(i)=Fred(iindex)
				endif
			enddo
			! Get Level Set Equation Terms
			call getTerms(M,MGLS,f1,f2,f3,Node,NLocal,NBelem,NBNodes,NBElems,&
				& LSetLocal,visc,h2,F)
			LSetLocal=LSetLocal-((((M+MGLS)^-1)*dt)*(f1+f2+f3))'
			! Reinitialize LS
			!call fastMarch(LSetLocal,NBelem,Node,NLocal,NBNodes)
			do i=1,NBNodes
				LSet(NLocal(i))=LSetLocal(i)
			end  
		enddo
	end	subroutine	
	!	This subroutine identifies elements in the narrow band
	subroutine getNarrowBand(NBelem,NBElems,NGlobal,NLocal,NBNodes,LSetLocal,&
				& bandWidth,LSet,Element,numElem,numNodes)
		include 'aba_param.inc'
		integer	  	Element(100,4),NBelem(100,4)
		real		Node(100,2),NGlobal(100,2),NLocal(100,2)
		real		LSet(100),LSetLocal(100)
		! Identify Narrow Band Elements
		NBElems=0
		NBNodes=0
		NGlobal=0.
		do i=1,numElem
			check=0
			do iNd=1,4
				 if (abs(LSet(Element(i,iNd)))<=bandWidth*(delX+delY)/2.)then
					 check=1
				 endif
			enddo
			! If an element is in the narrow band split it into triangles
			if (check==1)then
				 for j=1,4
					 if (NGlobal(Element(i,j))==0)then
						 NBNodes=NBNodes+1
						 NGlobal(Element(i,j))=NBNodes
						 NLocal(NBNodes)=Element(i,j)                 
					 endif
				 endddo
				 NBElems=NBElems+1
				 NBelem(NBElems,1)=NGlobal(Element(i,1))
				 NBelem(NBElems,2)=NGlobal(Element(i,2))
				 NBelem(NBElems,3)=NGlobal(Element(i,3))                     
				 NBElems=NBElems+1
				 NBelem(NBElems,1)=NGlobal(Element(i,1))
				 NBelem(NBElems,2)=NGlobal(Element(i,3))
				 NBelem(NBElems,3)=NGlobal(Element(i,4)) 
			endif
		enddo
		! Get local Level Set
		do i=1,NBNodes
			LSetLocal(i)=LSet(NLocal(i))
		enddo  		
	end subroutine	
	!	This subroutine extends the interface velocity throughout the computational domain
	subroutine getF(F,LSetLocal,NBElems,NBNodes,NBelem)  
		include 'aba_param.inc'
		real		LSetLocal(100),F(100),L(3)
		F=0.
		do i=1,NBElems 
			do j=1,3
				L(j)=sign(1.,LSetLocal(NBelem(i,j)))
			enddo
			if (L(1) /= L(2) .or. L(1) /= L(3))then
				do j=1,3
					F(NBelem(i,j))= 1.
				enddo
			endif
		enddo
	end	subroutine
	!	This subroutine gets the 'stiffness' matrix - A
	subroutine getA(A,Node,NLocal,NBelem,NBNodes,NBElems,LSetLocal)
		include 'aba_param.inc'
		integer	  	NBelem(100,4)
		real		Node(100,2),NLocal(100,2)
		real		LSetLocal(100),A(100,100),AfL(3,3),AfLGLS(3,3)
		real		gx(3),hx(3),phi(3),phig(3),phih(3),phix(3),phiy(3)
		A=0.
		do i=1,NBElems            
			gx(1)=2./3.
			gx(2)=1./6.
			gx(3)=1./6.
			hx(1)=1./6.
			hx(2)=1./6.
			hx(3)=2./3.
			AfL=0.
			AfLGLS=0.
			x1=Node(NLocal(NBelem(i,1)),1)
			y1=Node(NLocal(NBelem(i,1)),2)
			x2=Node(NLocal(NBelem(i,2)),1)
			y2=Node(NLocal(NBelem(i,2)),2)
			x3=Node(NLocal(NBelem(i,3)),1)
			y3=Node(NLocal(NBelem(i,3)),2)
			do j=1,3
				g=gx(j)
				h=hx(j)            
				phi(1)=1.-g-h
				phi(2)=g
				phi(3)=h
				phig(1)=-1.
				phig(2)=1.
				phig(3)=0.   
				phih(1)=-1.
				phih(2)=0.
				phih(3)=1.            
				djac=2.*abs(x1*(y2-y3)+x2*(y3-y1)+x3*(y1-y2))
				do k=1,3
					phix(k)=(1./djac)*((-y1+y3)*phig(k)+(y1-y2)*phih(k))
					phiy(k)=(1./djac)*((x1-x3)*phig(k)+(-x1+x2)*phih(k))
				enddo
				delphi=[phix;phiy]
				nodalLset=[LSetLocal(NBelem(i,1));LSetLocal(NBelem(i,2));LSetLocal(NBelem(i,3))]
				set=phi*nodalLset
				delset=delphi*nodalLset
				AfL=AfL+(phi'*sign(set))*(delset'*delphi)/3.
				AfLGLS=AfLGLS+(delphi'*delset)*(1./norm(delset))*(delset'*delphi)/3.       
			enddo
			do k=1,3
				do j=1,3
				   A(NBelem(i,j),NBelem(i,k))=A(NBelem(i,j),NBelem(i,k))+AfL(j,k)+AfLGLS(j,k)
				enddo
			enddo
		enddo    
	end subroutine
	! This subroutine gets the neccessary terms for the level set equation
	subroutine getTerms(M,MGLS,f1,f2,f3,Node,NLocal,NBelem,NBNodes,NBElems,&
				& LSetLocal,visc,h2,F)
		include 'aba_param.inc'
		integer	  	NBelem(100,4)
		real		Node(100,2),NLocal(100,2),LSetLocal(100)	
		real 		M(100,100),MGLS(100,100),f1(100),f2(100),f3(100)	
		real 		ML(3,3),MGLSL(3,3),f1L(3),f2L(3),f3L(3)	
		real		gx(3),hx(3),phi(3),phig(3),phih(3),phix(3),phiy(3)		
		
		M=0.
		MGLS=0.
		f1=0.
		f2=0.
		f3=0.
		do i=1,NBElems  
			ML=0.
			MGLSL=0.
			f1L=0.
			f2L=0.
			f3L=0.
			gx(1)=2./3.
			gx(2)=1./6.
			gx(3)=1./6.
			hx(1)=1./6.
			hx(2)=1./6.
			hx(3)=2./3.
			x1=Node(NLocal(NBelem(i,1)),1)
			y1=Node(NLocal(NBelem(i,1)),2)
			x2=Node(NLocal(NBelem(i,2)),1)
			y2=Node(NLocal(NBelem(i,2)),2)
			x3=Node(NLocal(NBelem(i,3)),1)
			y3=Node(NLocal(NBelem(i,3)),2)
			do j=1,3
				g=gx(j)
				h=hx(j)            
				phi(1)=1.-g-h
				phi(2)=g
				phi(3)=h
				phig(1)=-1.
				phig(2)=1.
				phig(3)=0.    
				phih(1)=-1.
				phih(2)=0.
				phih(3)=1.            
				djac=abs(x1*(y2-y3)+x2*(y3-y1)+x3*(y1-y2))
				do k=1,3
					phix(k)=(1./djac)*((-y1+y3)*phig(k)+(y1-y2)*phih(k))
					phiy(k)=(1./djac)*((x1-x3)*phig(k)+(-x1+x2)*phih(k))
				end 
				delphi=[phix;phiy]
				nodalLset=[LSetLocal(NBelem(i,1));LSetLocal(NBelem(i,2));LSetLocal(NBelem(i,3))]
				nodalF=[F(NBelem(i,1));F(NBelem(i,2));F(NBelem(i,3))]
				delset=delphi*nodalLset;                             
				Floc=phi*nodalF
				ML=ML+(phi'*phi)/3.
				MGLSL=MGLSL+((delphi'*(delset/norm(delset)))*Floc*(h2/abs(Floc)))*phi/3.
				f1L=f1L+phi'*Floc*norm(delset)/3.
				f2L=f2L+(delphi'*(delset/norm(delset))*Floc)*(h2/abs(Floc))*Floc*norm(delset)/3.
				vs=h2*((abs(visc+Floc*norm(delset)))/(norm(Floc*delset)+h2))
				f3L=f3L+vs*delphi'*delset/3.
			enddo
			do k=1,3
				do j=1,3
				   M(NBelem(i,j),NBelem(i,k))=M(NBelem(i,j),NBelem(i,k))+ML(j,k)
				   MGLS(NBelem(i,j),NBelem(i,k))=MGLS(NBelem(i,j),NBelem(i,k))+MGLSL(j,k)          
				enddo
				f1(NBelem(i,k))=f1(NBelem(i,k))+f1L(k)
				f2(NBelem(i,k))=f2(NBelem(i,k))+f2L(k)
				f3(NBelem(i,k))=f3(NBelem(i,k))+f3L(k)      
			enddo
		enddo   
	end subroutine	
	!	This subroutine uses the fast marching method to re-initialize the level set
	call fastMarch(LSetLocal,NBelem,Node,NLocal,NBNodes)
		include 'aba_param.inc'
		integer	  	NBelem(100,4),nstat(100)
		real		Node(100,2),NLocal(100,2),LSetLocal(100),newlSet(100),L(3)	
		newlSet=LSetLocal
		! Reinitialize LS
		nstat=0
		do i=1,NBElems
			do j=1,3
				L(j)=sign(1.,lSetLocal(NBelem(i,j)))
			enddo
			if (L(1) /= L(2) .or. L(1) /= L(3))then      
				do j=1,3
					nstat(NBelem(i,j))=1
				enddo
			endif
		enddo  
		maincheck=0
		do while(maincheck==0)
			lmin=1000.
			avlmin=1000.
			eindex=0
			nindex=0
			maincheck=1
			do i=1,NBElems
				if (nstat(NBelem(i,1))+nstat(NBelem(i,2))+nstat(NBelem(i,3))==2)then
					maincheck=0
					check=0
					ltot=0.
					do j=1,3
						if (nstat(NBelem(i,j))==0)then
							if (abs(lSetLocal(NBelem(i,j)))<=lmin)then
								check=1
								tempindex=j
							endif
						endif
						ltot=ltot+abs(lSetLocal(NBelem(i,j)))
					enddo
					if (check==1 .and. ltot/3.<=avlmin)then
						eindex=i
						nindex=tempindex
						lmin=lSetLocal(NBelem(eindex,nindex))
						avlmin=ltot/3.
					endif
				endif
			enddo
			if (maincheck==0)then
				! Find New LS for point
				xp=Node(NLocal(NBelem(eindex,nindex)),1)
				yp=Node(NLocal(NBelem(eindex,nindex)),2)
				count=0
				do i=1,3
					if (i/=nindex)then
						icount=icount+1
						x(icount)=Node(NLocal(NBelem(eindex,i)),1)
						y(icount)=Node(NLocal(NBelem(eindex,i)),2)
						lloc(icount)=newlSet(NBelem(eindex,i))
					endif
				enddo
				delxa=x(1)-xp
				delya=y(1)-yp
				delxb=x(2)-xp
				delyb=y(2)-yp
				N=[delxa delya; delxb delyb]
				M=N^-1
				A=(M(1)*M(1)+M(2)*M(2))
				B=(M(3)*M(3)+M(4)*M(4))
				C=2.*(M(1)*M(3)+M(2)*M(4))
				a=A+B+C
				b=-2.*lloc(1)*A-2.*lloc(2)*B-C*(lloc(1)+lloc(2))
				c=lloc(1)*lloc(1)*A+lloc(2)*lloc(2)*B+lloc(1)*lloc(2)*C-1.
				templ1=(-b+sqrt(b*b-4.*a*c))/(2.*a)
				templ2=(-b-sqrt(b*b-4.*a*c))/(2.*a)
				if (abs(templ1)>abs(templ2))then
					newlSet(NBelem(eindex,nindex))=templ1
				else
					newlSet(NBelem(eindex,nindex))=templ2
				endif
				nstat(NBelem(eindex,nindex))=1 
			endif
		enddo
		LSetLocal=newlSet
	end	subroutine