2011-02-03

KABIRA2011-02-03

ナナメ方向への拡散
- こちらでナナメの拡散の係数を1/2倍にすればよいことを書いた
ナナメ具合（＝線維方向）があるときどのように計算するか
- セミナーでいわれた通りタテヨコナナメ４方向平等に扱う必要がある
計算（セミナーでやったもの）
- 線維方向
  - $\vec{l} = \begin{pmatrix}\cos\theta \\ \sin \theta \end{pmatrix}$
- 線維方向の拡散の定数： $k_l$
- 線維に垂直な方向の拡散の定数： $k_t$
- 行列 $R$ を以下のようにおく
- $R=(\vec{l}\text{ } \vec{t} ) = \begin{pmatrix} \cos\theta && -\sin\theta \\ \sin\theta && \cos\theta \end{pmatrix}$
- 行列 $G$ は以下のようになる
- $G=R\begin{pmatrix} k_l && 0 \\ 0 && k_t \end{pmatrix} R^T = \begin{pmatrix} k_l\cos^2\theta+k_t\sin^2\theta && (k_l-k_t)\sin\theta\cos\theta \\ (k_l-k_t)\sin\theta\cos\theta && k_l\sin^2\theta+k_t\cos^2\theta \end{pmatrix}$

$\begin{pmatrix}g_l\cos^2\theta\cos^2\phi + g_t(1-\cos^2\theta\cos^2\phi) & (g_l - g_t)\cos\theta \sin \theta \cos^2 \phi & (g_l - g_t) \cos \theta \cos \phi \sin \phi \\ (g_l - g_t)\cos\theta \sin \theta \cos^2 \phi & g_l \sin^2 \theta \cos^2 \phi +g_t (1-\sin^2 \theta \cos^2 \phi) & (g_l -g_t) \sin\theta \cos \phi \sin \phi \\ (g_l - g_t) \cos \theta \cos \phi \sin \phi & (g_l -g_t) \sin \theta \cos \phi \sin \phi & g_l \sin^2 \phi + g_t \cos^2 \phi \end{pmatrix}$

こちらのソースコードをもとにしてナナメの係数k3,k4を与える
- k1,..,k4 は $\theta$ によって与えられている
- 厳密にはk3やk4に代入するのかが疑問がのこるが...
$\theta$ の値を変えて確かめられる

Nx<-30
Ny<-30
Nt<-100
U<-tempU<-matrix(0,Nx,Ny)
U[10:20,10:20]<-1

#これがl方向とt方向の拡散の定数
kl<-0.05
kt<-0.25

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

k1<-kl*cos(theta)^2+kt*sin(theta)^2
k2<-kl*sin(theta)^2+kt*cos(theta)^2
k3<-(kl-kt)*sin(theta)*cos(theta)
k4<-(kl-kt)*sin(theta)*cos(theta)

#1/2に補正
k3<-k3/2
k4<-k4/2


#k1<-matrix(-0.1,Nx-1,Ny)
#k2<-matrix(0.1,Nx,Ny-1)
#k3<-matrix(0.05,Nx-1,Ny-1)
#k4<-matrix(0.05,Nx-1,Ny-1)


for (t in 1:Nt){
	dUx<-U[2:Nx,]-U[1:(Nx-1),]
	dUy<-U[,2:Ny]-U[,1:(Ny-1)]
	dUxy<-U[2:Nx,2:Ny]-U[1:Nx-1,1:Ny-1]
	dUyx<-U[1:Nx-1,2:Ny]-U[2:Nx,1:Ny-1]
	
	U[1:(Nx-1),]<-U[1:(Nx-1),]+k1*dUx
	U[2:Nx,]<-U[2:Nx,]-k1*dUx
	U[,1:(Ny-1)]<-U[,1:(Ny-1)]+k2*dUy
	U[,2:Ny]<-U[,2:Ny]-k2*dUy
	
	U[1:Nx-1,1:Ny-1]<-U[1:Nx-1,1:Ny-1]+k3*dUxy
	U[2:Nx,2:Ny]<-U[2:Nx,2:Ny]-k3*dUxy
#	U[1:Nx-1,2:Ny]<-U[1:Nx-1,2:Ny]-k4*dUyx
#	U[2:Nx,1:Ny-1]<-U[2:Nx,1:Ny-1]+k4*dUyx

	U[1:Nx-1,2:Ny]<-U[1:Nx-1,2:Ny]+k4*dUyx
	U[2:Nx,1:Ny-1]<-U[2:Nx,1:Ny-1]-k4*dUyx
	
#	persp(U,col=heat.colors(Nx*Ny),theta=30,phi=30,zlim=c(0,1))
	image(U,col=topo.colors(100))
#	print(sum(U))
}

螺旋状に線維を張って、線維にそって拡散させる

Nx<-30
Ny<-30
Nt<-200
U<-tempU<-matrix(0,Nx,Ny)
U[15:10,10:5]<-1

#これがl方向とt方向の拡散の定数
kl<-0.25
kt<-0.01

#線維走向の傾き
#theta<-pi/6
A<-matrix(1,Nx,Ny)
theta<-matrix(0,Nx,Ny)
for(i in 1:Nx){
for(j in 1:Ny){
	theta[i,j]<-atan2(j-Ny/2,i-Nx/2)+pi/2+0.2
}
}

k1<-kl*cos(theta)^2+kt*sin(theta)^2
k2<-kl*sin(theta)^2+kt*cos(theta)^2
k3<-(kl-kt)*sin(theta)*cos(theta)
k4<-(kl-kt)*sin(theta)*cos(theta)

#1/2に補正
k3<-k3/2
k4<-k4/2

k1<-k1[-Nx,]
k2<-k2[,-Ny]
k3<-k3[-Nx,-Ny]
k4<-k4[-Nx,-Ny]

#k1<-matrix(-0.1,Nx-1,Ny)
#k2<-matrix(0.1,Nx,Ny-1)
#k3<-matrix(0.05,Nx-1,Ny-1)
#k4<-matrix(0.05,Nx-1,Ny-1)


for (t in 1:Nt){
	dUx<-U[2:Nx,]-U[1:(Nx-1),]
	dUy<-U[,2:Ny]-U[,1:(Ny-1)]
	dUxy<-U[2:Nx,2:Ny]-U[1:Nx-1,1:Ny-1]
	dUyx<-U[2:Nx,1:Ny-1]-U[1:Nx-1,2:Ny]
	
	U[1:(Nx-1),]<-U[1:(Nx-1),]+k1*dUx
	U[2:Nx,]<-U[2:Nx,]-k1*dUx
	U[,1:(Ny-1)]<-U[,1:(Ny-1)]+k2*dUy
	U[,2:Ny]<-U[,2:Ny]-k2*dUy
	
	U[1:Nx-1,1:Ny-1]<-U[1:Nx-1,1:Ny-1]+k3*dUxy
	U[2:Nx,2:Ny]<-U[2:Nx,2:Ny]-k3*dUxy
	U[1:Nx-1,2:Ny]<-U[1:Nx-1,2:Ny]-k4*dUyx
	U[2:Nx,1:Ny-1]<-U[2:Nx,1:Ny-1]+k4*dUyx

#	U[1:Nx-1,2:Ny]<-U[1:Nx-1,2:Ny]+k4*dUyx
#	U[2:Nx,1:Ny-1]<-U[2:Nx,1:Ny-1]-k4*dUyx
	
	par(mfrow=c(1,2))
	layout(matrix(c(1,0,2,1),2,2))
	persp(U,col=topo.colors(Nx*Ny),theta=30,phi=30,zlim=c(0,1))
	image(U,col=topo.colors(100))
#	print(sum(U))
}

$\theta$ についてのメモ

e1_theta<-(theta[1:(Nx-1),]+theta[2:Nx,])/2
e2_theta<-(theta[,1:(Ny-1)]+theta[,2:Ny])/2
e3_theta<-(theta[2:Nx,2:Ny]+theta[1:Nx-1,1:Ny-1])/2
e4_theta<-(theta[1:Nx-1,2:Ny]+theta[2:Nx,1:Ny-1])/2

k1<-kl*cos(e1_theta)^2+kt*sin(e1_theta)^2
k2<-kl*sin(e2_theta)^2+kt*cos(e2_theta)^2
k3<-(kl-kt)*sin(e3_theta)*cos(e3_theta)
k4<-(kl-kt)*sin(e4_theta)*cos(e4_theta)


#1/2に補正
k3<-k3/2
k4<-k4/2
library(rgl)
Nx<-30
Ny<-30
Nz<-30
Nt<-10

R<-10
U<-tempU<-array(0,c(Nx,Ny,Nz))


#kl<-0.01
#kt<-0.01
#kd<-0.01

#theta<-array(pi/4)
#phi<-array(acos(sqrt(6)/3),c(Nx,Ny,Nz))

A<-array(0,c(Nx,Ny,Nz))
for(x in 1:Nx){
for(y in 1:Nx){
for(z in 1:Nz){
	if((x-Nx/2)^2+(y-Ny/2)^2+(z-Nz/2)^2 <= R^2){A[x,y,z]<-1} 
}}}

U[1:ceiling(Nx/2),1:ceiling(Ny/2),1:ceiling(Nz/2)]<-1

U<-U*A

#KX<-kl*cos(theta)^2*cos(phi)^2+kt*(1-cos(theta)^2*cos(phi)^2)
#KY<-kl*sin(theta)^2*cos(phi)^2+kt*(1-sin(theta)^2*cos(phi)^2)
#KZ<-kl*sin(theta)^2+kt*cos(phi)^2
#KXY<-KxY<-(kl-kt)*sin(theta)*cos(theta)*cos(phi)^2
#KYZ<-KyZ<-(kl-kt)*sin(theta)*cos(phi)*sin(phi)
#KZX<-KzX<-(kl-kt)*cos(theta)*cos(phi)*sin(phi)

#線維を考えずに定数にする場合
KX<-KY<-KZ<-array(0.01,c(Nx,Ny,Nz))
KXY<-KYZ<-KZX<-array(0.01,c(Nx,Ny,Nz))
KxY<-KyZ<-KzX<-array(0.01,c(Nx,Ny,Nz))


kX<-(KX[-Nx,,]+KX[-1,,])/2*A[-Nx,,]*A[-1,,]
kY<-(KY[,-Ny,]+KY[,-1,])/2*A[,-Ny,]*A[,-1,]
kZ<-(KZ[,,-Nz]+KZ[,,-1])/2*A[,,-Nz]*A[,,-1]

kXY<-(KXY[-1,-1,]+KXY[-Nx,-Ny,])/2*A[-Nx,-Ny,]*A[-1,-1,]
kYZ<-(KYZ[,-1,-1]+KYZ[,-Ny,-Nz])/2*A[,-Ny,-Nz]*A[,-1,-1]
kZX<-(KZX[-1,,-1]+KZX[-Nx,,-Nz])/2*A[-Nx,,-Nz]*A[-1,,-1]

kxY<-(KxY[-Nx,-1,]+KxY[-1,-Ny,])/2*A[-Nx,-Ny,]*A[-1,-1,]
kyZ<-(KyZ[,-Ny,-1]+KyZ[,-1,-Nz])/2*A[,-Ny,-Nz]*A[,-1,-1]
kzX<-(KzX[-1,,-Nz]+KzX[-Nx,,-1])/2*A[-Nx,,-Nz]*A[-1,,-1]

kXY<-kXY/2
kYZ<-kYZ/2
kZX<-kZX/2

kxY<-kxY/2
kyZ<-kyZ/2
kzX<-kzX/2

for (t in 1:Nt){
	kdUX<-(U[-1,,]-U[-Nx,,])*kX
	kdUY<-(U[,-1,]-U[,-Ny,])*kY
	kdUZ<-(U[,,-1]-U[,,-Nz])*kZ
	
	kdUXY<-(U[-1,-1,]-U[-Nx,-Ny,])*kXY
	kdUYZ<-(U[,-1,-1]-U[,-Ny,-Nz])*kYZ
	kdUZX<-(U[-1,,-1]-U[-Nx,,-Nz])*kZX
	
	kdUxY<-(U[-Nx,-1,]-U[-1,-Ny,])*kxY
	kdUyZ<-(U[,-Ny,-1]-U[,-1,-Nz])*kyZ
	kdUzX<-(U[-1,,-Nz]-U[-Nx,,-1])*kzX
	
		
	U[-Nx,,]<-U[-Nx,,]+kdUX
	U[-1,,]<-U[-1,,]-kdUX
	U[,-Ny,]<-U[,-Ny,]+kdUY
	U[,-1,]<-U[,-1,]-kdUY
	U[,,-Nz]<-U[,,-Nz]+kdUZ
	U[,,-1]<-U[,,-1]-kdUZ
	
	U[-Nx,-Ny,]<-U[-Nx,-Ny,]+kdUXY
	U[-1,-1,]<-U[-1,-1,]-kdUXY
	U[,-Ny,-Nz]<-U[,-Ny,-Nz]+kdUYZ
	U[,-1,-1]<-U[,-1,-1]-kdUYZ
	U[-Nx,,-Nz]<-U[-Nx,,-Nz]+kdUZX
	U[-1,,-1]<-U[-1,,-1]-kdUZX
	
	U[-Nx,-1,]<-U[-Nx,-1,]-kdUxY
	U[-1,-Ny,]<-U[-1,-Ny,]+kdUxY
	U[,-Ny,-1]<-U[,-Ny,-1]-kdUyZ
	U[,-1,-Nz]<-U[,-1,-Nz]+kdUyZ
	U[-1,,-Nz]<-U[-1,,-Nz]-kdUzX
	U[-Nx,,-1]<-U[-Nx,,-1]+kdUzX
	
	plot3d(slice.index(U,1),slice.index(U,2),slice.index(U,3),col=rainbow(1000)[ifelse(U < 0.001,1,U*1000)],alpha=ifelse(U < 0.001,0,1))
	rgl.viewpoint(theta=300,phi=30)

#	print(sum(U));print(max(U))
}
plot3d(slice.index(U,1),slice.index(U,2),slice.index(U,3),col=rainbow(1000)[ifelse(U < 0.001, 1,U*1000)],alpha=ifelse(U < 0.001,0,1))