subroutine uel(rhs,amatrx,svars,energy,ndofel,nrhs,nsvars,props 1 ,nprops,coords,mcrd,nnode,u,du,v,a,jtype,time,dtime,kstep, 2 kinc,jelem,params,ndload,jdltyp,adlmag,predef,npredf,lflags 3 ,mlvarx,ddlmag,mdload,pnewdt,jprops,njprop,period) c include 'aba_param.inc' c dimension rhs(mlvarx,*),amatrx(ndofel,ndofel),svars(nsvars), 1 energy(8),props(*),coords(mcrd,nnode), 2 u(ndofel),du(mlvarx,*),v(ndofel),a(ndofel),time(2), 3 params(3),jdltyp(mdload,*),adlmag(mdload,*), 4 ddlmag(mdload,*),predef(2,npredf,nnode),lflags(*),jprops(*) c dimension gpx(9),gwei(9),phi(8),phix(8),phic(8) c c print *,u(1),u(2),du(1,nhrs),du(2,nhrs),time(1),lflags(3) c level set calculation dpos=0.01*time c c material property definition rho = 1. spec = 1. conduc = 1. c initialization (nrhs=1) do k1=1,ndofel rhs(k1,nrhs)=0. do k2=1,ndofel amatrx(k2,k1)=0. enddo enddo if (lflags(3).eq.4) return c transient analysis if (lflags(1).eq.33) then c determine gauss point locations gpx(1)=-1./sqrt(3.) gpx(2)=1./sqrt(3.) gwei(1)=1. gwei(2)=1. c determine node level sets rnl1=abs(coords(1,1)-dpos) rnl2=abs(coords(1,2)-dpos) c get jacobian dxdc=abs(coords(1,2)-coords(1,1))/2. ajacob=dxdc c assemble amatrx and rhs do k=1,4 c loop through gauss pts c=gpx(k) c get shape functions and derivatives phi(1)=(1.-c)/2. phi(2)=(1.+c)/2. c get ip position pos=coord(1,1)+ajacob(1.+c) phi(3)=abs(c) phic(1)=-0.5 phic(2)=0.5 phix(1)=phic(1)*(1./ajacob) phix(2)=phic(2)*(1./ajacob) c interpolate temperatures to int points dtdx=u(1)*phix(1)+u(2)*phix(2) t=u(1)*phi(1)+u(2)*phi(2) told=(u(1)-du(1,nrhs))*phi(1)+(u(2)-du(2,nrhs))*phi(2) c other housekeeping cond=conduc dtdt=(t-told)/dtime we=gwei(k)*ajacob c Assemble Element Stiffness Matrix and Add to Global do ki=1,2 c loop over nodes rhs(ki,nrhs)=rhs(ki,nrhs)-we*(phi(ki)*rho*spec*dtdt 1 + cond*(phix(ki)*dtdx)) do kj=1,2 amatrx(ki,kj)=amatrx(ki,kj)+we*(phi(ki)*phi(kj) 1 *rho*spec/dtime+cond*(phix(ki)*phix(kj))) end do end do enddo end if return end