ソースコードを簡潔に直したいのですが…
有限積分法(参考 http://www.akita-nct.jp/yamamoto/study/thesis/2005/thesis_namekawa.pdf)
を用いて、格子点 ( i, j, k )の磁場成分 Bx, By, Bz をその周辺を囲む電場成分 Ex, Ey, Ez で逐次計算するコードは次のようになります。
for ( i = 0; i < i_max; i++ ) {
for ( j = 0; j < j_max; j++ ) {
for ( k = 0; k < k_max; k++ ) {
Bx [ i ][ j ][ k ] = Ey [ i + 1 ][ j ][ k ] - Ey [ i + 1 ][ j ][ k + 1 ] + Ez [ i + 1 ][ j + 1 ][ k ] - Ez [ i + 1 ][ j ][ k ];
By [ i ][ j ][ k ] = Ez [ i ][ j + 1 ][ k ] - Ez [ i + 1 ][ j + 1 ][ k ] + Ex [ i ][ j + 1 ][ k + 1 ] - Ex [ i ][ j + 1 ][ k ];
Bz [ i ][ j ][ k ] = Ex [ i ][ j ][ k + 1 ] - Ex [ i ][ j + 1 ][ k + 1 ] + Ey [ i + 1 ][ j ][ k + 1 ] - Ey [ i ][ j ][ k + 1 ];
}
}
}
電場成分の係数行列 C ( 0, 1, -1 のいずれかをもつ ) を求める必要が生じた、すなわち
CCC E B
CCC * E = B
CCC E B
のような行列計算に変更するため、次のように書き直しましたが係数行列 C を求める部分を簡潔にできませんでした。
for ( i = 0; i < i_max; i++ ) {
for ( j = 0; j < j_max; j++ ) {
for ( k = 0; k < k_max; k++ ) {
for ( x = 0; x < 3; x++ ) {
row = index ( x, i, j, k, i_max, j_max, k_max );
y = ( x + 1 ) % 3;
z = ( x + 2 ) % 3;
if ( x == 0 ) {
C[ row ][ index ( y, i + 1, j, k, i_max, j_max, k_max )] = 1;
C[ row ][ index ( y, i + 1, j, k + 1, i_max, j_max, k_max )] = -1;
C[ row ][ index ( z, i + 1, j + 1, k, i_max, j_max, k_max )] = 1;
C[ row ][ index ( z, i + 1, j, k, i_max, j_max, k_max )] = -1;
}
if ( x == 1 ) {
C[ row ][ index ( y, i, j + 1, k, i_max, j_max, k_max )] = 1;
C[ row ][ index ( y, i + 1, j + 1, k, i_max, j_max, k_max )] = -1;
C[ row ][ index ( z, i, j + 1, k + 1, i_max, j_max, k_max )] = 1;
C[ row ][ index ( z, i, j + 1, k, i_max, j_max, k_max )] = -1;
}
if ( x == 2 ) {
C[ row ][ index ( y, i, j, k + 1, i_max, j_max, k_max )] = 1;
C[ row ][ index ( y, i, j + 1, k + 1, i_max, j_max, k_max )] = -1;
C[ row ][ index ( z, i + 1, j, k + 1, i_max, j_max, k_max )] = 1;
C[ row ][ index ( z, i, j, k + 1, i_max, j_max, k_max )] = -1;
}
}
}
}
}
int index ( const int d, const int i, const int j, const int k, const int i_max, const int j_max, const int k_max )
{
return ( ( d * i_max + i ) * j_max + j ) * k_max + k;
}
void Solver ( const int i_max, const int j_max, const int k_max)
{
int row, col,
cell_num = i_max * j_max * k_max;
for ( row = 0; row < cell_num; row++ ) {
for ( col = 0; col < cell_num; col++ ) {
B [ row ] += C [ row ][ col ] * E [ col ];
}
}
}
添え字計算が多く、非常に複雑なソースコードですが、規則性があるので、
皆様の力をお借りして簡潔に表現したいです。
よろしくお願いします。
