c*********************** problem name: jcn ************************ 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 common /val0/k0,ku,kx,ky,kl,kuu,kux,kuy,kul, + kxu,kxx,kxy,kxl,kyu,kyx,kyy,kyl,klu,klx,kly,kll cy call cbeta(x,y,betax,betay,itag) values(k0)=ux+u*betax values(ku)=betax values(kx)=1.0e0_rknd 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 common /val0/k0,ku,kx,ky,kl,kuu,kux,kuy,kul, + kxu,kxx,kxy,kxl,kyu,kyx,kyy,kyl,klu,klx,kly,kll cy call cbeta(x,y,betax,betay,itag) values(k0)=uy+u*betay values(ku)=betay values(ky)=1.0e0_rknd 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 common /val1/k0,ku,kl,kuu,kul,klu,kll cy 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/ival,iob,iu(98),top,bottom,du, + angle,ct,st,xt(7),yt(7),ru(80),su(100) cy values(k0)=top if(itag==1) values(k0)=bottom values(kic)=values(k0) 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 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 common /val3/kf,kf1,kf2,kad cy za=0.5e0_rknd*abs(u) if(za>1.0e5_rknd) then as=2.0e0_rknd*za else as=za+sqrt(1.0e0_rknd+za*za) endif c if(u<0.0e0_rknd) then v=-log10(as) else v=log10(as) endif c values(kf)=v 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 return end c----------------------------------------------------------------------- c c piecewise lagrange triangle multi grid package c c edition 11.0 - - - june, 2012 c c----------------------------------------------------------------------- subroutine cbeta(x,y,betax,betay,itag) cx use mthdef implicit real(kind=rknd) (a-h,o-z) implicit integer(kind=iknd) (i-n) character(len=80) :: su common /atest2/ival,iob,iu(98),top,bottom,du, + angle,ct,st,xt(7),yt(7),ru(80),su(100) cy betax=du*(ct*xt(itag)+st*yt(itag)) betay=du*(ct*yt(itag)-st*xt(itag)) 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 character(len=80), dimension(100) :: sp,su character(len=80), save, dimension(10) :: file cy data len/7/ data (file(i),i= 1, 7)/ + 'n i= 1,n=obtuse,a= o,t=i', 1 'n i= 1,n=top ,a= t,t=r', 2 'n i= 2,n=bottom,a= b,t=r', 3 'n i= 3,n=du ,a= d,t=r', 4 'n i= 4,n=angle ,a= a,t=r', 5 's n=obtuse,v=0,l="do not trisect elements"', 6 's n=obtuse,v=1,l="trisect all elements"'/ c aa=ru(4) c call usrset(file,len,iu,ru,su) c if(iu(1)/=0) ip(41)=-1 if(ru(4)/=aa) then pi=3.141592653589793e0_rknd theta=ru(4)*pi/180.0e0_rknd ru(5)=cos(theta) ru(6)=sin(theta) endif 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, dimension(23) :: b1,b2,b4 integer(kind=iknd), save, dimension(7) :: ib1,ib2 integer(kind=iknd), save, dimension(11) :: ix,iy integer(kind=iknd), save :: ntf,nvf,nbf,ispd,iprob integer(kind=iknd), save :: ifirst real(kind=rknd), dimension(*) :: vx,vy real(kind=rknd), dimension(2,*) :: sf real(kind=rknd), dimension(100) :: rp,ru real(kind=rknd), save :: hmax,grade character(len=80), dimension(100) :: sp,su cy external sxy c data b1/ 1, 2, 3, 4, 1, 6, 7, 6, 8, 9, + 8,10,11,12,13,14,15,16,17, 7, 1 12,13,14/ data b2/ 2, 3, 4, 5, 6, 7,12, 8, 9,15, + 10,11, 5,13,14, 3,16,17, 4, 9, 1 15,16,17/ data b4/ 2, 1, 1, 1, 1, 0, 0, 1, 0, 0, + 1, 2, 1, 0, 0, 0, 0, 0, 0, 0, 1 0, 0, 0/ c data ib1/ 8,6, 7,12,14,14,3/ data ib2/11,8,20,21,15,23,2/ c data ix/175,200,225,250,275,175,200,175,200,175,275/ data iy/0,0,0,0,0,-25,-25,-50,-50,-100,-100/ data ntf,nvf,nbf,ispd,iprob,ifirst/7,17,23,0,1,1/ data hmax,grade/0.1e0_rknd,1.5e0_rknd/ c c if(ip(41)==1) then iu(1)=0 ru(1)=1.0e10_rknd ru(2)=1.0e-5_rknd ru(3)=(40.0e0_rknd+15.0e0_rknd*log(10.0e0_rknd)) ru(4)=0.0e0_rknd ru(5)=1.0e0_rknd ru(6)=0.0e0_rknd endif c if(iu(1)/=0) go to 100 c c set up regions c sp(2)='continuity equation' sp(1)='continuity equation' sp(3)='continuity equation' sp(6)='jcn.rw' sp(7)='jcn.jnl' sp(9)='jcn.out' c ip(5)=max(ip(5),ifirst) ip(6)=iprob ip(8)=ispd ip(19)=1 ip(20)=5 rp(15)=hmax rp(16)=grade c ss=2.99563e-2_rknd pi=3.141592653589793e0_rknd c s1=ss*1.0e-2_rknd do i=1,11 vx(i)=real(ix(i),rknd)*s1 vy(i)=real(iy(i),rknd)*s1 enddo do i=1,3 arg=(12.0e0_rknd+real(i,rknd))*pi/8.0e0_rknd c=cos(arg)*ss s=sin(arg)*ss vx(11+i)=vx(2)+0.25e0_rknd*c vy(11+i)=0.25e0_rknd*s vx(14+i)=vx(2)+0.5e0_rknd*c vy(14+i)=0.5e0_rknd*s enddo c c constants for a1xy and a2xy c ru(7)=0.0e0_rknd ru(14)=0.0e0_rknd c ru(8)=0.0e0_rknd ru(15)=1.0e0_rknd/(vy(9)-vy(7)) c xx=(vx(9)-vx(7)+vx(15)-vx(12))/2.0e0_rknd yy=(vy(9)-vy(7)+vy(15)-vy(12))/2.0e0_rknd dd=1.0e0_rknd/(xx*xx+yy*yy) ru(9)=xx*dd ru(16)=yy*dd c xx=(vx(15)-vx(12)+vx(16)-vx(13))/2.0e0_rknd yy=(vy(15)-vy(12)+vy(16)-vy(13))/2.0e0_rknd dd=1.0e0_rknd/(xx*xx+yy*yy) ru(10)=xx*dd ru(17)=yy*dd c xx=(vx(16)-vx(13)+vx(17)-vx(14))/2.0e0_rknd yy=(vy(16)-vy(13)+vy(17)-vy(14))/2.0e0_rknd dd=1.0e0_rknd/(xx*xx+yy*yy) ru(11)=xx*dd ru(18)=yy*dd c xx=(vx(17)-vx(14)+vx(4)-vx(3))/2.0e0_rknd yy=(vy(17)-vy(14)+vy(4)-vy(3))/2.0e0_rknd dd=1.0e0_rknd/(xx*xx+yy*yy) ru(12)=xx*dd ru(19)=yy*dd c ru(13)=0.0e0_rknd ru(20)=0.0e0_rknd c do i=1,nbf ibndry(1,i)=b1(i) ibndry(2,i)=b2(i) ibndry(3,i)=0 ibndry(4,i)=b4(i) ibndry(5,i)=0 ibndry(6,i)=0 ibndry(7,i)=0 sf(1,i)=0.0e0_rknd sf(2,i)=0.0e0_rknd enddo ibndry(7,12)=1 ibndry(7,1)=3 c do i=1,ntf itnode(1,i)=ib1(i) itnode(2,i)=ib2(i) itnode(3,i)=0 itnode(4,i)=0 itnode(5,i)=i enddo c ip(1)=ntf ip(2)=nvf ip(3)=nbf return c c make obtuse angles c 100 iu(1)=0 ntf=ip(1) nt0=ntf nvf=ip(2) nbf=ip(3) maxt=ip(83) maxv=ip(84) if (ntf*3>maxt) return if (nvf+ntf>maxv) return do it=1,nt0 iv1=itnode(1,it) iv2=itnode(2,it) iv3=itnode(3,it) iv4=itnode(4,it) iv5=itnode(5,it) c nvf=nvf+1 vx(nvf)=(vx(iv1)+vx(iv2)+vx(iv3))/3.0e0_rknd vy(nvf)=(vy(iv1)+vy(iv2)+vy(iv3))/3.0e0_rknd c itnode(1,ntf+1)=nvf itnode(2,ntf+1)=iv2 itnode(3,ntf+1)=iv3 itnode(4,ntf+1)=iv4 itnode(5,ntf+1)=iv5 itnode(1,ntf+2)=iv1 itnode(2,ntf+2)=nvf itnode(3,ntf+2)=iv3 itnode(4,ntf+2)=iv4 itnode(5,ntf+2)=iv5 itnode(3,it)=nvf ntf=ntf+2 enddo ip(1)=ntf ip(2)=nvf return end