c*********************** problem name: circle ************************ c----------------------------------------------------------------------- c c piecewise lagrange triangle multi grid package c c edition 11.0 - - - june, 2012 c c----------------------------------------------------------------------- subroutine a1xy(x,y,u,ux,uy,rl,itag,values) cx use mthdef implicit real(kind=rknd) (a-h,o-z) implicit integer(kind=iknd) (i-n) real(kind=rknd), dimension(*) :: values character(len=80) :: su common /atest2/iu(100),ru(100),su(100) common /val0/k0,ku,kx,ky,kl,kuu,kux,kuy,kul, + kxu,kxx,kxy,kxl,kyu,kyx,kyy,kyl,klu,klx,kly,kll cy values(k0)=ux*ru(itag) values(kx)=ru(itag) return end c----------------------------------------------------------------------- c c piecewise lagrange triangle multi grid package c c edition 11.0 - - - june, 2012 c c----------------------------------------------------------------------- subroutine a2xy(x,y,u,ux,uy,rl,itag,values) cx use mthdef implicit real(kind=rknd) (a-h,o-z) implicit integer(kind=iknd) (i-n) real(kind=rknd), dimension(*) :: values character(len=80) :: su common /atest2/iu(100),ru(100),su(100) common /val0/k0,ku,kx,ky,kl,kuu,kux,kuy,kul, + kxu,kxx,kxy,kxl,kyu,kyx,kyy,kyl,klu,klx,kly,kll cy values(k0)=uy*ru(itag) values(ky)=ru(itag) return end c----------------------------------------------------------------------- c c piecewise lagrange triangle multi grid package c c edition 11.0 - - - june, 2012 c c----------------------------------------------------------------------- subroutine fxy(x,y,u,ux,uy,rl,itag,values) cx use mthdef implicit real(kind=rknd) (a-h,o-z) implicit integer(kind=iknd) (i-n) real(kind=rknd), dimension(*) :: values common /val0/k0,ku,kx,ky,kl,kuu,kux,kuy,kul, + kxu,kxx,kxy,kxl,kyu,kyx,kyy,kyl,klu,klx,kly,kll cy return end c----------------------------------------------------------------------- c c piecewise lagrange triangle multi grid package c c edition 11.0 - - - june, 2012 c c----------------------------------------------------------------------- subroutine gnxy(x,y,u,rl,itag,values) cx use mthdef implicit real(kind=rknd) (a-h,o-z) implicit integer(kind=iknd) (i-n) real(kind=rknd), dimension(*) :: values character(len=80) :: su common /val1/k0,ku,kl,kuu,kul,klu,kll common /atest2/iu(100),ru(100),su(100) cy if(itag>iu(1)) return c c (x,y) is unit outward normal c if(itag<1) return call uexact(x,y,itag,r,rx,ry,rxx,ryy,rxy) if(itag<=iu(1)) then values(k0)=(x*rx+y*ry)*ru(itag) else values(k0)=ry*ru(itag-1) endif return end c----------------------------------------------------------------------- c c piecewise lagrange triangle multi grid package c c edition 11.0 - - - june, 2012 c c----------------------------------------------------------------------- subroutine gdxy(x,y,rl,itag,values) cx use mthdef implicit real(kind=rknd) (a-h,o-z) implicit integer(kind=iknd) (i-n) real(kind=rknd), dimension(*) :: values character(len=80) :: su common /val2/k0,kl,kll,klb,kub,kic,kim,kil common /atest2/iu(100),ru(100),su(100) cy if(itag<1.or.itag>iu(1)) return call uexact(x,y,itag,r,rx,ry,rxx,ryy,rxy) values(k0)=r return end c----------------------------------------------------------------------- c c piecewise lagrange triangle multi grid package c c edition 11.0 - - - june, 2012 c c----------------------------------------------------------------------- subroutine p1xy(x,y,u,ux,uy,rl,itag,values) cx use mthdef implicit real(kind=rknd) (a-h,o-z) implicit integer(kind=iknd) (i-n) real(kind=rknd), dimension(*) :: values common /val0/k0,ku,kx,ky,kl,kuu,kux,kuy,kul, + kxu,kxx,kxy,kxl,kyu,kyx,kyy,kyl,klu,klx,kly,kll cy xx=0.5e0_rknd yy=0.5e0_rknd e=1.0e4_rknd*exp(-5.0e0_rknd*((x-xx)**2+(y-yy)**2)) values(k0)=e*u**2 values(ku)=e*2.0e0_rknd*u values(kuu)=e*2.0e0_rknd return end c----------------------------------------------------------------------- c c piecewise lagrange triangle multi grid package c c edition 11.0 - - - june, 2012 c c----------------------------------------------------------------------- subroutine p2xy(x,y,dx,dy,u,ux,uy,rl,itag,jtag,values) cx use mthdef implicit real(kind=rknd) (a-h,o-z) implicit integer(kind=iknd) (i-n) real(kind=rknd), dimension(*) :: values common /val0/k0,ku,kx,ky,kl,kuu,kux,kuy,kul, + kxu,kxx,kxy,kxl,kyu,kyx,kyy,kyl,klu,klx,kly,kll cy return end c----------------------------------------------------------------------- c c piecewise lagrange triangle multi grid package c c edition 11.0 - - - june, 2012 c c----------------------------------------------------------------------- subroutine qxy(x,y,u,ux,uy,rl,itag,values) cx use mthdef implicit real(kind=rknd) (a-h,o-z) implicit integer(kind=iknd) (i-n) real(kind=rknd), dimension(*) :: values character(len=80) :: su common /val3/kf,kf1,kf2,kad common /atest2/iu(100),ru(100),su(100) cy call uexact(x,y,itag,r,rx,ry,rxx,ryy,rxy) values(kf)=r-u values(kf1)=rx-ux values(kf2)=ry-uy values(kad)=r return end c----------------------------------------------------------------------- c c piecewise lagrange triangle multi grid package c c edition 11.0 - - - june, 2012 c c----------------------------------------------------------------------- subroutine sxy(rl,s,itag,values) cx use mthdef implicit real(kind=rknd) (a-h,o-z) implicit integer(kind=iknd) (i-n) real(kind=rknd), dimension(*) :: values common /val4/jx,jy,jxs,jys,jxl,jyl,jxss,jyss,jxll,jyll, + jxsl,jysl,jxls,jyls cy values(jx)=cos(s) values(jy)=sin(s) values(jxs)=-values(jy) values(jys)=values(jx) values(jxss)=-values(jx) values(jyss)=-values(jy) return end c----------------------------------------------------------------------- c c piecewise lagrange triangle multi grid package c c edition 11.0 - - - june, 2012 c c----------------------------------------------------------------------- subroutine uexact(x,y,itag,u,ux,uy,uxx,uyy,uxy) cx use mthdef implicit real(kind=rknd) (a-h,o-z) implicit integer(kind=iknd) (i-n) character(len=80) :: su common /atest2/iu(100),ru(100),su(100) cy u=0.0e0_rknd ux=0.0e0_rknd uy=0.0e0_rknd uxx=0.0e0_rknd uxy=0.0e0_rknd uyy=0.0e0_rknd r=sqrt(x**2+y**2) if(r<=0.0e0_rknd) return c pi=3.141592653589793e0_rknd al=ru(25) arg=min(x/r,1.0e0_rknd) arg=max(-1.0e0_rknd,arg) theta=acos(arg) if(itag>=5) theta=2.0e0_rknd*pi-theta c s=r**al*sin(theta*al) al1=al-1.0e0_rknd al2=al-2.0e0_rknd sx=al*r**al1*sin(theta*al1) sy=al*r**al1*cos(theta*al1) sxx=al*al1*r**al2*sin(theta*al2) syy=-sxx sxy=al*al1*r**al2*cos(theta*al2) c c=r**al*cos(theta*al) cx=sy cy=-sx cxx=sxy cyy=-sxy cxy=syy c cf=ru(itag+8) sf=ru(itag+16) u=cf*c+sf*s ux=cf*cx+sf*sx uy=cf*cy+sf*sy uxx=cf*cxx+sf*sxx uyy=cf*cyy+sf*syy uxy=cf*cxy+sf*sxy c return end c----------------------------------------------------------------------- c c piecewise lagrange triangle multi grid package c c edition 11.0 - - - june, 2012 c c----------------------------------------------------------------------- subroutine cnu(ntri,af,cf,sf,vv) cx use mthdef implicit real(kind=rknd) (a-h,o-z) implicit integer(kind=iknd) (i-n) real(kind=rknd), dimension(8) :: cf,sf,af,cd,sd real(kind=rknd), dimension(1000) :: val cy eps=1.0e1_rknd*epsilon(1.0e0_rknd) cf(1)=0.0e0_rknd sf(1)=1.0e0_rknd cd(1)=0.0e0_rknd sd(1)=0.0e0_rknd vmax=2.0e0_rknd n=1000 do i=1,n vv=real(i,rknd)*vmax/real(n,rknd) call veval(ntri,cf,sf,af,vv,val(i),dp) enddo do i=2,n if(val(i-1)>0.0e0_rknd.and.val(i)<=0.0e0_rknd) then vmin=real(i-1,rknd)*vmax/real(n,rknd) vmax=real(i,rknd)*vmax/real(n,rknd) go to 10 endif enddo vv=vmax return c 10 itmax=100 vv=(vmin+vmax)/2.0e0_rknd do k=1,itmax call veval(ntri,cf,sf,af,vv,dd,dp) if(dd>0.0e0_rknd) then vmin=vv else vmax=vv endif vv=vv-dd/dp if(vv<=vmin.or.vv>=vmax) + vv=(vmin+vmax)/2.0e0_rknd if(abs(dd)<=eps) return enddo return end c----------------------------------------------------------------------- c c piecewise lagrange triangle multi grid package c c edition 11.0 - - - june, 2012 c c----------------------------------------------------------------------- subroutine veval(ntri,cf,sf,af,vv,dd,dp) cx use mthdef implicit real(kind=rknd) (a-h,o-z) implicit integer(kind=iknd) (i-n) real(kind=rknd), dimension(8) :: cf,sf,af,cd,sd cy pi=3.141592653589793e0_rknd cf(1)=0.0e0_rknd sf(1)=1.0e0_rknd cd(1)=0.0e0_rknd sd(1)=0.0e0_rknd c do i=1,ntri-1 if(af(i+1)==af(i)) then cf(i+1)=cf(i) sf(i+1)=sf(i) cd(i+1)=cd(i) sd(i+1)=sd(i) else theta=real(i,rknd)*pi/4.0e0_rknd c=cos(theta*vv) s=sin(theta*vv) dc=-theta*s ds=theta*c c ff=c*cf(i)+s*sf(i) dd=-s*cf(i)+c*sf(i) df=c*cd(i)+s*sd(i)+dc*cf(i)+ds*sf(i) dp=-s*cd(i)+c*sd(i)-ds*cf(i)+dc*sf(i) c dd=dd*af(i)/af(i+1) dp=dp*af(i)/af(i+1) c cf(i+1)=c*ff-s*dd sf(i+1)=s*ff+c*dd cd(i+1)=dc*ff-ds*dd+c*df-s*dp sd(i+1)=ds*ff+dc*dd+s*df+c*dp endif enddo theta=real(ntri,rknd)*pi/4.0e0_rknd c=cos(theta*vv) s=sin(theta*vv) dc=-theta*s ds=theta*c c dd=-s*cf(ntri)+c*sf(ntri) dp=-ds*cf(ntri)+dc*sf(ntri)-s*cd(ntri)+c*sd(ntri) return end c----------------------------------------------------------------------- c c piecewise lagrange triangle multi grid package c c edition 11.0 - - - june, 2012 c c----------------------------------------------------------------------- subroutine usrcmd(vx,vy,sf,itnode,ibndry,ip,rp,sp,iu,ru,su) cx use mthdef implicit real(kind=rknd) (a-h,o-z) implicit integer(kind=iknd) (i-n) integer(kind=iknd), dimension(5,*) :: itnode integer(kind=iknd), dimension(7,*) :: ibndry integer(kind=iknd), dimension(100) :: ip,iu integer(kind=iknd), save :: len real(kind=rknd), dimension(*) :: vx,vy real(kind=rknd), dimension(2,*) :: sf real(kind=rknd), dimension(100) :: rp,ru real(kind=rknd), dimension(8) :: psave character(len=80), dimension(100) :: sp,su character(len=80), save, dimension(12) :: file cy data len/12/ data (file(i),i=1,10)/ + 'n i=1,n=a1,a=a1,t=r', 1 'n i=2,n=a2,a=a2,t=r', 2 'n i=3,n=a3,a=a3,t=r', 3 'n i=4,n=a4,a=a4,t=r', 4 'n i=5,n=a5,a=a5,t=r', 5 'n i=6,n=a6,a=a6,t=r', 6 'n i=7,n=a7,a=a7,t=r', 7 'n i=8,n=a8,a=a8,t=r', 8 'n i=1,n=ntri,a=nt,t=i', 9 'n i=2,n=ibc,a=bc,t=i'/ data (file(i),i=11,12)/ + 's n=ibc,v=1,l="neumann"', 1 's n=ibc,v=2,l="dirichlet"'/ c ntf=iu(1) ibc=iu(2) do i=1,8 psave(i)=ru(i) enddo c call usrset(file,len,iu,ru,su) c iu(1)=min(8,iu(1)) iu(1)=max(1,iu(1)) if(iu(2)/=1) iu(2)=2 do i=1,8 if(ru(i)<=0.0e0_rknd) ru(i)=1.0e0_rknd enddo c if(ntf/=iu(1).or.ibc/=iu(2)) ip(41)=-1 do i=1,8 if(psave(i)/=ru(i)) ip(41)=-1 enddo return end c----------------------------------------------------------------------- c c piecewise lagrange triangle multi grid package c c edition 11.0 - - - june, 2012 c c----------------------------------------------------------------------- subroutine gdata(vx,vy,sf,itnode,ibndry,ip,rp,sp,iu,ru,su,sxy) cx use mthdef implicit real(kind=rknd) (a-h,o-z) implicit integer(kind=iknd) (i-n) integer(kind=iknd), dimension(5,*) :: itnode integer(kind=iknd), dimension(7,*) :: ibndry integer(kind=iknd), dimension(100) :: ip,iu integer(kind=iknd), save :: ispd=1,iprob=1,itask=0,ifirst=1 real(kind=rknd), dimension(*) :: vx,vy real(kind=rknd), dimension(2,*) :: sf real(kind=rknd), dimension(100) :: rp,ru character(len=80), dimension(100) :: sp,su cy external sxy c c unit circle divided into 8 equal parts. c this problem can be run with ntf=1,2,...8, c nvf=ntf+2,nbf=ntf+2 c if(ip(41)==1) then sp(1)='circle' sp(2)='circle' sp(3)='circle' sp(6)='circle_mpixxx.rw' sp(7)='circle.jnl' sp(9)='circle_mpixxx.out' c iu(1)=8 iu(2)=2 do i=1,8 ru(i)=1.0e0_rknd enddo endif c ntf=iu(1) ibc=iu(2) nvf=ntf+2 nbf=ntf+2 ip(1)=ntf ip(2)=nvf ip(3)=nbf ip(5)=max(ip(5),ifirst) ip(22)=100000 ip(6)=iprob ip(8)=ispd ip(7)=itask c pi=3.141592653589793e0_rknd do i=1,ntf itnode(1,i)=1 itnode(2,i)=i+1 itnode(3,i)=i+2 itnode(4,i)=0 itnode(5,i)=i enddo vx(1)=0.0e0_rknd vy(1)=0.0e0_rknd do i=2,nvf arg=pi*real(i-2,rknd)/4.0e0_rknd vx(i)=cos(arg) vy(i)=sin(arg) enddo do i=1,nbf ibndry(1,i)=i ibndry(2,i)=i+1 ibndry(3,i)=1 ibndry(4,i)=ibc ibndry(5,i)=0 ibndry(6,i)=0 ibndry(7,i)=i-1 sf(1,i)=0.0e0_rknd sf(2,i)=0.0e0_rknd enddo ibndry(2,nbf)=1 ibndry(3,1)=0 ibndry(3,nbf)=0 ibndry(4,1)=2 ibndry(4,nbf)=1 c pi4=pi/4.0e0_rknd do i=2,nbf-1 sf(1,i)=real(i-2,rknd)*pi4 sf(2,i)=real(i-1,rknd)*pi4 ibndry(3,i)=-i enddo c c call cnu(ntf,ru(1),ru(9),ru(17),ru(25)) c sp(1)='alpha = ' call sreal(sp(1)(9:9),nn,ru(25),3,1) sp(2)=sp(1) sp(3)=sp(1) return end