拡散の方向性

• ナナメ方向への拡散をあらわす
  • こちらの前半にある式を正確に書き表すことにした
  • 基本的には拡散項に関してこのときのような計算方法
  • ナナメの近傍への移動はない
  • つまり、k3やk4は出てこない
• といった条件では形がいびつになることがある

Nx<-30
Ny<-30
Nt<-100
U<-matrix(0,Nx,Ny)
U[11:20,11:20]<-1
#これがl方向とt方向の拡散の定数
kl<-0.20
kt<-0.05

#線維走向の傾き
theta<-pi/4

kxx<-kl*cos(theta)^2+kt*sin(theta)^2
kyy<-kl*sin(theta)^2+kt*cos(theta)^2
kxy<-(kl-kt)*sin(theta)*cos(theta)
kyx<-(kl-kt)*sin(theta)*cos(theta)

for (t in 1:Nt){
	dUx<-U[-1,]-U[-Nx,]
	dUy<-U[,-1]-U[,-Ny]
	
	Ex<--(rbind(dUx,dUx[Nx-1,])+rbind(dUx[1,],dUx))/2
	Ey<--(cbind(dUy,dUy[,Ny-1])+cbind(dUy[,1],dUy))/2

#	Ixx<-Ex*kxx
#	Iyy<-Ey*kyy
	Ixy<-Ey*kxy
	Iyx<-Ex*kyx
	
#	Ixxp<-(Ixx>0)*Ixx
#	Ixxm<-(Ixx<0)*Ixx
#	Iyyp<-(Iyy>0)*Iyy
#	Iyym<-(Iyy<0)*Iyy
	Ixyp<-(Ixy>0)*Ixy
	Ixym<-(Ixy<0)*Ixy
	Iyxp<-(Iyx>0)*Iyx
	Iyxm<-(Iyx<0)*Iyx
	
#	Ix<-(Ixxp+Ixyp)[-Nx,]+(Ixxm+Ixym)[-1,]
#	Iy<-(Iyyp+Iyxp)[,-Ny]+(Iyym+Iyxm)[,-1]

	Ix<--dUx*kxx+Ixyp[-Nx,]+Ixym[-1,]
	Iy<--dUy*kyy+Iyxp[,-Ny]+Iyxm[,-1]
	
		
	U[-1,]<-U[-1,]+Ix
	U[-Nx,]<-U[-Nx,]-Ix
	U[,-1]<-U[,-1]+Iy
	U[,-Ny]<-U[,-Ny]-Iy

	
#	persp(U,col=heat.colors(Nx*Ny),theta=30,phi=30,zlim=c(0,1))
	image(U,col=topo.colors(100))
#	image(Iy)
#	persp(Ix,theta =30,phi=30,col=2)
#	print(sum(U))
}

Nx<-30
Ny<-30
Nt<-100
U<-matrix(0,Nx,Ny)
U[11:20,11:20]<-1
#これがl方向とt方向の拡散の定数
kl<-0.25
kt<-0.05

#線維走向の傾き
theta<-pi/4

kxx<-kl*cos(theta)^2+kt*sin(theta)^2
kyy<-kl*sin(theta)^2+kt*cos(theta)^2
kxy<-(kl-kt)*sin(theta)*cos(theta)
kyx<-(kl-kt)*sin(theta)*cos(theta)

for (t in 1:Nt){
	dUx<-U[-1,]-U[-Nx,]
	dUy<-U[,-1]-U[,-Ny]
	
	Ex<--(rbind(dUx,dUx[Nx-1,])+rbind(dUx[1,],dUx))/2
	Ey<--(cbind(dUy,dUy[,Ny-1])+cbind(dUy[,1],dUy))/2

	Ixx<-Ex*kxx
	Iyy<-Ey*kyy
	Ixy<-Ey*kxy
	Iyx<-Ex*kyx
	
##	ここは使わない	
##	Ixxp<-(Ixx>0)*Ixx
##	Ixxm<-(Ixx<0)*Ixx
##	Iyyp<-(Iyy>0)*Iyy
##	Iyym<-(Iyy<0)*Iyy
##	Ixyp<-(Ixy>0)*Ixy
##	Ixym<-(Ixy<0)*Ixy
##	Iyxp<-(Iyx>0)*Iyx
##	Iyxm<-(Iyx<0)*Iyx

	Ix<--dUx*kxx+(Ixy[-Nx,]+Ixy[-1,])/2
	Iy<--dUy*kyy+(Iyx[,-Ny]+Iyx[,-1])/2

#	Ix<-Ixx+Ixy	
#	Iy<-Iyx+Iyy
	
#	Ix<-(Ix[-1,]+Ix[-Nx,])/2
#	Iy<-(Iy[,-1]+Iy[,-Ny])/2
	
	U[-1,]<-U[-1,]+Ix
	U[-Nx,]<-U[-Nx,]-Ix
	U[,-1]<-U[,-1]+Iy
	U[,-Ny]<-U[,-Ny]-Iy

	
#	persp(U,col=heat.colors(Nx*Ny),theta=30,phi=30,zlim=c(0,1))
	image(U,col=topo.colors(100))
#	image(Iy)
#	persp(Ix,theta =30,phi=30,col=2)
#	print(sum(U))
}