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
.