オブジェクト指向プログラミング

自分が考えたクラス継承は以下のような感じです。

まず3次元格子を表すクラスとしてGridクラスがあります。これは3次元空間に物理量を1つだけ持ったものです。

Gridを継承してMatrixBaseクラスがあります。これはGridが保持している値に対して差分を計算する機能を付加したものです。

次のMatrixVは、MatrixBaseを速度を扱う様に特化したものです。MatrixPは圧力、MatrixTは温度を扱う様に特化したものです。

さらにMatrixVは、これを3つ使ってとし、新しいクラスMatrixV3Dのメンバーとしています。

そして、速度のMatrixV3D、圧力のMatrixP、温度のMatrixTをSolver3Dのメンバーとしています。

最初はこの様な構成ではなく、速度も圧力も必要な物理量を全て1つのクラスに固めて1セルとし、それを3次元に展開していこうと思っていました。しかし微分を計算するのに各物理量に対して同じコードを何回も書かねばなず、どうしてもおかしい、という事になってしまいました。そこで必要な物理量を全て1つのクラスに固める、という方針を捨て、3次元空間に1つの物理量しかない、というクラス(Grid)を基本として、それを必要な物理量だけ並行して存在させる、という事にしました。3方向の速度も、その一つ一つは圧力のようなスカラー量と変わらない、という所がミソですね。こんな感じでコーディングしたところ、から新しいを計算する部分など、ほとんど数式と同じに書くことができて、ちょっとびっくりです。


に対応するMatrixV3D内のメソッドは以下の様です。

private double calc( MatrixP p, int component ){
 return ( v(component)
         +m_dt * ( -( uRep( component ) * v[ component ].centralDifference1(_x)
                     +vRep( component ) * v[ component ].centralDifference1(_y)
                     +wRep( component ) * v[ component ].centralDifference1(_z) )
                   + mReInverse * ( v[ component ].centralDifference2(_x)
                                   +v[ component ].centralDifference2(_y)
                                   +v[ component ].centralDifference2(_z) )
                   - p.backwardDifference( component ) 
      )
   );
}

もうちょっとがんばればナブラ()自体もコード化できるかも知れません。

しかしながらMatrixVの内部では、生成の際に自分がなのかなのかなのか指定するようになっています。内部でも方向を完全に区別してコーディングしています。我ながらこの辺は固いコーディングですね。本当ならこの辺も、などと決め打ちにせず、自分に対して水平、垂直といった考え方で抽象化できたハズですが、自分の頭の限界で途中で挫折してしまいました。まあ、空間次元はこれ以上増えも減りもしないので良しとしました。

今回のソルバーでは温度もブジネスク近似で考慮していますが、これもMatrixBaseから派生させたMatrixTで計算します。このような構成であれば密度を変数とした場合も割と簡単に組み込めると思います。
また微分に関してはMatrixBase内だけで計算しています。全てのコード中で微分を計算しているのはここ以外にありません。”同じコードを2度書くな”がオブジェクト指向の鉄則ですからね。なので現在は計算精度が1次ですが、このMatrixBaseの計算制度を上げてやれば、全ての計算の精度が上がることになります。こうゆうスマートさが、自分がオブジェクト指向を気に入っている点です。

おまけで説明すると、物理量へのアクセス機能はGridクラスにあるのですが、基本的にアクセスは着目位置(i、j、k)にオフセット+1などとして行います。そして着目位置は変数m_i、m_j、m_kに記憶しているのですが、この変数をスタティックとしています。なのでこのスタティック変数を動かせば、Gridから派生したクラスならすべて着目位置が動く事になります。同期し忘れ、と言う事は起こりません。

下に演算に関係ある部分のコードリストを示します。ご自分でご利用なさる限り自由に使って頂いて構いません。


/******************************************************************************
  
  Copyright(C) 2007 DeepDigital Co.,Ltd. All Rights Reserved.
  
 ******************************************************************************/

package math;

import javax.vecmath.Vector3d;

public class ValueVector3d extends Vector3d {

    public ValueVector3d( double x1, double x2, double x3 ){
        super( x1,x2,x3 );
    }

    public ValueVector3d(){
        super();
    }

    public double getXYRad(){
        return Math.atan2( y,x );
    }

    public double getXZRad(){
        return -Math.atan2( z,getXYLength() );
    }

    private double getXYLength(){
        return Math.pow( ( Math.pow(x, 2.0) + Math.pow(y, 2.0) ), 0.5 );
    }
}



/******************************************************************************
  
  Copyright(C) 2007 DeepDigital Co.,Ltd. All Rights Reserved.
  
 ******************************************************************************/

package FDM;

public interface Global {

    public final int    _x=0, _y=1, _z=2, _cellCenter=3;
    public final int    _u=0, _v=1, _w=2;

}

/******************************************************************************
  
  Copyright(C) 2007 DeepDigital Co.,Ltd. All Rights Reserved.
  
 ******************************************************************************/

package FDM;

public class BasicInfo implements Global{

    protected static int        mMaxX,mMaxY,mMaxZ;
    
    protected static double[]   m_hInverse,m_h2Inverse;
    protected static double     mReInverse;
    protected static double     mRePrInverce;
    protected static double     mGrPerRe2;
    protected static double     m_dt, m_dtInverse;
    protected static int        mCaseNo;
    

    public static void setBasicInfo( int maxX, int maxY, int maxZ,
                                        double re, double pr, double gr, int caseNo ) {
        mMaxX = maxX;
        mMaxY = maxY;
        mMaxZ = maxZ;
        m_hInverse = new double[3];
        m_h2Inverse = new double[3];
        m_hInverse[_x] = (mMaxX-2);
        m_hInverse[_y] = (mMaxY-2);
        m_hInverse[_z] = (mMaxZ-2);
        for( int i=0; i<=2; i++)
            m_h2Inverse[i] = Math.pow( m_hInverse[i],2.0 );
        
        mReInverse      = 1.0/re;
        double max1     = Math.max(m_hInverse[_x],m_hInverse[_y]);
        double max2     = Math.max(max1,m_hInverse[_z]);
        m_dt            = 1.0/max2;
        m_dt            = 0.5/max2;//temp
        m_dtInverse     = 1.0/m_dt;

        mRePrInverce    = mReInverse/pr;
        mGrPerRe2       = gr*mReInverse*mReInverse;
        
        mCaseNo         = caseNo;
    }

}

/******************************************************************************
  
  Copyright(C) 2007 DeepDigital Co.,Ltd. All Rights Reserved.
  
 ******************************************************************************/

package FDM;

import math.ValueVector3d;

public class Grid extends BasicInfo{

    protected float[][][]       q;
    protected static int        m_i,m_j,m_k;
    
    public Grid() {
        q = new float[mMaxX+1][mMaxY+1][mMaxZ+1];
    }
    public Grid( float value ) {
        q = new float[mMaxX+1][mMaxY+1][mMaxZ+1];

        for( int k=0; k <= mMaxZ; k++ ){
            for( int j=0; j <= mMaxY; j++ ){
                for( int i=0; i <= mMaxX; i++ ){
                    q[i][j][k] = value;
                }
            }
        }
    }

    public static void setPosition( int i, int j, int k ){
        m_i = i;
        m_j = j;
        m_k = k;
    }

    public void setValue( double value ){
        q[m_i][m_j][m_k] = (float)value;
    }

    public double q( int i, int j, int k ){
        return q[i][j][k];
    }
    protected double q(){
        return q[m_i][m_j][m_k];
    }
    protected double q( int iDirection, int iStep ){
        switch( iDirection ){
            case _x:
                return q[m_i+iStep][m_j][m_k];
            case _y:
                return q[m_i][m_j+iStep][m_k];
            case _z:
                return q[m_i][m_j][m_k+iStep];
        }
        return 9999;
    }

    public ValueVector3d getValue( int i, int j, int k ){
        return new ValueVector3d( 0,0,q[i][j][k] );
    }
    
}

/******************************************************************************
  
  Copyright(C) 2007 DeepDigital Co.,Ltd. All Rights Reserved.
  
 ******************************************************************************/

package FDM;

public class MatrixBase extends Grid {

    public MatrixBase(){
    }
    public MatrixBase( float value ){
        super( value );
    }

    
    protected double centralDifference1( int direction ){
        return ( q( direction, +1 ) - q( direction, -1 ) )
                * m_hInverse[ direction ] * 0.5;
    }
    protected double centralDifference2( int direction ){
        return ( q( direction, +1 ) - q()*2.0 + q( direction, -1 ) )
                * m_h2Inverse[ direction ];
    }
    
    protected double forwardDifference( int direction ){
        return ( q( direction, +1 ) - q() ) * m_hInverse[ direction ];
    }
    protected double backwardDifference( int direction ){
        return ( q() - q( direction, -1 ) ) * m_hInverse[ direction ];
    }

}

/******************************************************************************
  
  Copyright(C) 2007 DeepDigital Co.,Ltd. All Rights Reserved.
  
 ******************************************************************************/

package FDM;

public class MatrixV extends MatrixBase {

    private final int           mComponent;
    
    public MatrixV( int component ){
        mComponent = component;
    }

    protected double represent( int currentComponent ){
        return represent( currentComponent, 0 );
    }
    protected double represent( int currentComponent, int step ){
        switch( mComponent ){
            case _u:
                switch( currentComponent ){
                    case _u:
                        return q[m_i+step][m_j][m_k];
                    case _v:
                        return ( q[m_i+step][m_j][m_k]   + q[m_i+step+1][m_j][m_k]
                                + q[m_i+step][m_j-1][m_k] + q[m_i+step+1][m_j-1][m_k] ) * 0.25;
                    case _w:
                        return ( q[m_i+step][m_j][m_k]   + q[m_i+step+1][m_j][m_k]
                                + q[m_i+step][m_j][m_k-1] + q[m_i+step+1][m_j][m_k-1] ) * 0.25;
                    case _cellCenter:
                        return ( q[m_i+step][m_j][m_k] + q[m_i+step+1][m_j][m_k] ) * 0.5;
                }
                break;
            case _v:
                switch( currentComponent ){
                    case _u:
                        return ( q[m_i][m_j+step][m_k] + q[m_i][m_j+step+1][m_k]
                                + q[m_i-1][m_j+step][m_k] + q[m_i-1][m_j+step+1][m_k] ) * 0.25;
                    case _v:
                        return q[m_i][m_j+step][m_k];
                    case _w:
                        return ( q[m_i][m_j+step][m_k]   + q[m_i][m_j+step+1][m_k]
                                + q[m_i][m_j+step][m_k-1] + q[m_i][m_j+step+1][m_k-1] ) * 0.25;
                    case _cellCenter:
                        return ( q[m_i][m_j+step][m_k] + q[m_i][m_j+step+1][m_k] ) * 0.5;
                }
                break;
            case _w:
                switch( currentComponent ){
                    case _u:
                        return ( q[m_i][m_j][m_k+step] + q[m_i][m_j][m_k+step+1]
                                + q[m_i-1][m_j][m_k+step] + q[m_i-1][m_j][m_k+step+1] ) * 0.25;
                    case _v:
                        return ( q[m_i][m_j][m_k+step]   + q[m_i][m_j][m_k+step+1]
                                + q[m_i][m_j-1][m_k+step] + q[m_i][m_j-1][m_k+step+1] ) * 0.25;
                    case _w:
                        return q[m_i][m_j][m_k+step];
                    case _cellCenter:
                        return ( q[m_i][m_j][m_k+step] + q[m_i][m_j][m_k+step+1] ) * 0.5;
                }
                break;
        }
        return 9999;
    }
    
    protected void setBoundaryCondition(){
        switch( mComponent ){
            case _u:
                //x section
                for( int j=0; j <= mMaxY; j++ ){
                    for( int k=0; k <= mMaxZ; k++ ){
                        //max
                        q[mMaxX  ][j][k]    =  q[mMaxX-2][j][k];
                        //max-1
                        q[mMaxX-1][j][k]    =  0.0f;
                        //zero
                        q[   0   ][j][k]    =  q[   2   ][j][k];
                        //one
                        q[   1   ][j][k]    =  0.0f;
                    }
                }

                if(mCaseNo==2){
                    //air conditioning IN
                    for( int j=2*mMaxY/6; j <= 3*mMaxY/6; j++ ){
                        //max-1
                        q[mMaxX  ][j][mMaxZ-3]  =  0.5f;
                        q[mMaxX-1][j][mMaxZ-3]  =  0.5f;
                        q[mMaxX  ][j][mMaxZ-4]  =  0.5f;
                        q[mMaxX-1][j][mMaxZ-4]  =  0.5f;
                    }
                    //air conditioning OUT
                    for( int j=2*mMaxY/6; j <= 3*mMaxY/6; j++ ){
                        //max-1
                        q[mMaxX  ][j][mMaxZ-5]  =  -1.0f;
                        q[mMaxX-1][j][mMaxZ-5]  =  -1.0f;
                    }
                }
                if(mCaseNo==4){
                    //air conditioning OUT
                    for( int j=1; j <= 3; j++ ){
                        for( int k=mMaxZ-4; k <= mMaxZ-2; k++ ){
                            q[mMaxX  ][j][k]    = -1.0f;
                            q[mMaxX-1][j][k]    = -1.0f;
                        }
                    }
                }
                //y section
                for( int k=0; k <= mMaxZ; k++ ){
                    for( int i=0; i <= mMaxX; i++ ){
                        //max
                        q[i][mMaxY  ][k]    =  0.0f;
                        //max-1
                        q[i][mMaxY-1][k]    = -q[i][mMaxY-2][k];
                        //zero
                        q[i][   0   ][k]    = -q[i][   1   ][k];
                    }
                }
                //z section
                for( int i=0; i <= mMaxX; i++ ){
                    for( int j=0; j <= mMaxY; j++ ){
                        //max
                        q[i][j][mMaxZ  ]    =  0.0f;
                        //max-1
                        q[i][j][mMaxZ-1]    = -q[i][j][mMaxZ-2];
                        if(mCaseNo==1){
                            //キャビティ速度
                            q[i][j][mMaxZ-1]    =  1.0f;
                        }
                        //zero
                        q[i][j][   0   ]    = -q[i][j][   1   ];
                    }
                }
                break;
            case _v:
                //x section
                for( int j=0; j <= mMaxY; j++ ){
                    for( int k=0; k <= mMaxZ; k++ ){
                        //max
                        q[mMaxX  ][j][k]    =  0.0f;
                        //max-1
                        q[mMaxX-1][j][k]    = -q[mMaxX-2][j][k];
                        //zero
                        q[   0   ][j][k]    = -q[   1   ][j][k];
                    }
                }
                //y section
                for( int k=0; k <= mMaxZ; k++ ){
                    for( int i=0; i <= mMaxX; i++ ){
                        //max
                        q[i][mMaxY  ][k]    =  q[i][mMaxY-2][k];
                        //max-1
                        q[i][mMaxY-1][k]    = 0.0f;
                        //zero
                        q[i][   0   ][k]    =  q[i][   2   ][k];
                        //one
                        q[i][   1   ][k]    =  0.0f;
                    }
                }
                //z section
                for( int i=0; i <= mMaxX; i++ ){
                    for( int j=0; j <= mMaxY; j++ ){
                        //max
                        q[i][j][mMaxZ  ]    =  0.0f;
                        //max-1
                        q[i][j][mMaxZ-1]    = -q[i][j][mMaxZ-2];
                        //zero
                        q[i][j][   0   ]    = -q[i][j][   1   ];
                    }
                }
                break;
            case _w:
                //x section
                for( int j=0; j <= mMaxY; j++ ){
                    for( int k=0; k <= mMaxZ; k++ ){
                        //max
                        q[mMaxX  ][j][k]    =  0.0f;
                        //max-1
                        q[mMaxX-1][j][k]    = -q[mMaxX-2][j][k];
                        //zero
                        q[   0   ][j][k]    = -q[   1   ][j][k];
                    }
                }
                //y section
                for( int k=0; k <= mMaxZ; k++ ){
                    for( int i=0; i <= mMaxX; i++ ){
                        //max
                        q[i][mMaxY  ][k]    =  0.0f;
                        //max-1
                        q[i][mMaxY-1][k]    = -q[i][mMaxY-2][k];
                        //zero
                        q[i][   0   ][k]    = -q[i][   1   ][k];
                    }
                }
                //z section
                for( int i=0; i <= mMaxX; i++ ){
                    for( int j=0; j <= mMaxY; j++ ){
                        //max
                        q[i][j][mMaxZ  ]    =  q[i][j][mMaxZ-2];
                        //max-1
                        q[i][j][mMaxZ-1]    =  0.0f;
                        //zero
                        q[i][j][   0   ]    =  q[i][j][   2   ];
                        //one
                        q[i][j][   1   ]    =  0.0f;
                    }
                }
                if(mCaseNo==4){
                    //air conditioning IN
                    //max-1
                    for( int i=(mMaxX-2)/2-1; i <= (mMaxX-2)/2+1; i++ ){
                        for( int j=(mMaxY-2)/2-1; j <= (mMaxY-2)/2+1; j++ ){
                            q[i][j][mMaxZ  ]    =  1.0f;
                            q[i][j][mMaxZ-1]    =  1.0f;
                        }
                    }
                }
                break;
        }
    }
}

/******************************************************************************
  
  Copyright(C) 2007 DeepDigital Co.,Ltd. All Rights Reserved.
  
 ******************************************************************************/

package FDM;

public class MatrixT extends MatrixBase {

    float WALL_TEMP;
    
    public MatrixT() {
        super(0.5f);
    }

    public void calcAll( MatrixV3D V ){
        setBoundaryCondition();
        for( int k=1; k <= mMaxZ-2; k++ ){
            for( int j=1; j <= mMaxY-2; j++ ){
                for( int i=1; i <= mMaxX-2; i++ ){
                    setPosition( i,j,k );
                    setValue( calc( V ) );
                }
            }
        }
    }
    
    private double calc( MatrixV3D V ){
        return q() + m_dt * ( -( V.uRep(_cellCenter) * centralDifference1(_x)
                                + V.vRep(_cellCenter) * centralDifference1(_y)
                                + V.wRep(_cellCenter) * centralDifference1(_z) )
                        + mRePrInverce * ( centralDifference2(_x)
                                            +centralDifference2(_y)
                                            +centralDifference2(_z)));
    }
    
    private void setBoundaryCondition(){
        WALL_TEMP=0.5f;
        //x section
        for( int j=0; j <= mMaxY; j++ ){
            for( int k=0; k <= mMaxZ; k++ ){
                //max
                q[mMaxX  ][j][k]    =  0.0f;
                //max-1
                //q[mMaxX-1][j][k]  =  q[mMaxX-2][j][k];
                q[mMaxX-1][j][k]    =  WALL_TEMP;
                //zero
                //q[   0   ][j][k]  =  q[   1   ][j][k];
                q[   0   ][j][k]    =  WALL_TEMP;
            }
        }

        if(mCaseNo==2){
            //air conditioning OUT
            for( int j=2*mMaxY/6; j <= 3*mMaxY/6; j++ ){
                //max-2
                q[mMaxX-1][j][mMaxZ-5]  =  0.0f;
            }
        }
        if(mCaseNo==4){
            //air conditioning OUT
            for( int j=1; j <= 3; j++ ){
                for( int k=mMaxZ-4; k <= mMaxZ-2; k++ ){
                    q[mMaxX-2][j][k]    = 1.0f;
                }
            }
        }
        //y section
        for( int i=0; i <= mMaxX; i++ ){
            for( int k=0; k <= mMaxZ; k++ ){
                //max
                q[i][mMaxY  ][k]    =  0.0f;
                //max-1
                //q[i][mMaxY-1][k]  =  q[i][mMaxY-2][k];
                q[i][mMaxY-1][k]    =  WALL_TEMP;
                //zero
                //q[i][   0   ][k]  =  q[i][   1   ][k];
                q[i][   0   ][k]    =  WALL_TEMP;
            }
        }
        //z section
        for( int i=0; i <= mMaxX; i++ ){
            for( int j=0; j <= mMaxY; j++ ){
                //max
                q[i][j][mMaxZ  ]    =  0.0f;
                //max-1
                //q[i][j][mMaxZ-1]  =  q[i][j][mMaxZ-2];
                q[i][j][mMaxZ-1]    =  WALL_TEMP;
                //zero
                //q[i][j][   0   ]  =  q[i][j][   1   ];
                q[i][j][   0   ]    =  WALL_TEMP;
            }
        }

        if(mCaseNo==2){
            //z section heating
            for( int i=2*mMaxX/3; i <= mMaxX-3; i++ ){
                for( int j=mMaxY/3; j <= 2*mMaxY/3; j++ ){
                    //zero
                    q[i][j][   0   ]    =  1.0f;
                }
            }
        }
        if(mCaseNo==3){
            //z section heating
            for( int i=(mMaxX-1)/2-1; i <= (mMaxX-1)/2+1; i++ ){
                for( int j=(mMaxY-1)/2-1; j <= (mMaxX-1)/2+1; j++ ){
                    //zero
                    q[i][j][   0   ]    =  1.0f;
                }
            }
        }
    }
}

/******************************************************************************
  
  Copyright(C) 2007 DeepDigital Co.,Ltd. All Rights Reserved.
  
 ******************************************************************************/

package FDM;

public class MatrixP extends MatrixBase {

    protected final double  _A;
    protected final double  _MAX_ITERATION;

    protected Grid              Q;
    
    public MatrixP() {
        
        _A = 1.0 / ( 2.0*( m_h2Inverse[_x] + m_h2Inverse[_y] + m_h2Inverse[_z] ) );
        _MAX_ITERATION=500;

        Q =     new Grid();
    }

    public void calc( MatrixV3D V ){
        double preError=-1.0, error=0.0;
        double pD=0.0;
        
        resetQ( V );
        for( int iteration=0; iteration<_MAX_ITERATION; iteration++ ){
            error=0.0;
            setBoundaryCondition( V );
            for( int i=1; i <= mMaxX-2; i++ ){
                for( int j=1; j <= mMaxY-2; j++ ){
                    for( int k=1; k <= mMaxZ-2; k++ ){
                        setPosition( i,j,k );
                        pD = calc();
                        error+=Math.abs( q()-pD );
                        setValue( pD );
                    }
                }
            }
            System.out.println( "error=" + error );
            if( error / preError > 0.9999 ) return;
            preError = error;
        }
        System.out.println( "exceeded iteration loop" );
    }

    private double calc(){
        return _A * ( ppDifference(_x) + ppDifference(_y) + ppDifference(_z) - Q.q() );
    }
    
    private double ppDifference( int iDirection ){
        return ( q( iDirection, +1 ) + q( iDirection, -1 ) ) * m_h2Inverse[ iDirection ];
    }

    private void resetQ( MatrixV3D V ){
        for( int i=1; i <= mMaxX-2; i++ ){
            for( int j=1; j <= mMaxY-2; j++ ){
                for( int k=1; k <= mMaxZ-2; k++ ){
                    setPosition( i,j,k );
                    Q.setValue( calcQ( V ) );
                }
            }
        }
    }
    private double calcQ( MatrixV3D V ){
        return ( (  V.v[_u].forwardDifference(_x)
                    +V.v[_v].forwardDifference(_y)
                    +V.v[_w].forwardDifference(_z) ) * m_dtInverse
                - (  Math.pow( V.v[_u].forwardDifference(_x), 2.0 )
                    +Math.pow( V.v[_v].forwardDifference(_y), 2.0 )
                    +Math.pow( V.v[_w].forwardDifference(_z), 2.0 )
                + ( V.vRepForwardDifference(_x) * V.uRepForwardDifference(_y)
                   +V.wRepForwardDifference(_x) * V.uRepForwardDifference(_z)
                   +V.wRepForwardDifference(_y) * V.vRepForwardDifference(_z) ) * 2.0 ) );
    }
    
    private void setBoundaryCondition( MatrixV3D V ){
        //x section
        for( int j=0; j <= mMaxY; j++ ){
            for( int k=0; k <= mMaxZ; k++ ){
                //max
                q[mMaxX  ][j][k]    =  (float)0.0;
                //max-1
                q[mMaxX-1][j][k]    =  (float)(q[mMaxX-2][j][k] - mReInverse*2.0*m_hInverse[_x] * V.v[_u].q[mMaxX-2][j][k]);
                //zero
                q[   0   ][j][k]    =  (float)(q[   1   ][j][k] - mReInverse*2.0*m_hInverse[_x] * V.v[_u].q[   2   ][j][k]);
            }
        }
        //case 1
        //air conditioning IN
        for( int j=2*mMaxY/6; j <= 3*mMaxY/6; j++ ){
            //max-1
            q[mMaxX-1][j][mMaxZ-3]  =  q[mMaxX-2][j][mMaxZ-3];
            q[mMaxX-1][j][mMaxZ-4]  =  q[mMaxX-2][j][mMaxZ-4];
            q[mMaxX-1][j][mMaxZ-5]  =  q[mMaxX-2][j][mMaxZ-5];
        }
        
        //y section
        for( int i=0; i <= mMaxX; i++ ){
            for( int k=0; k <= mMaxZ; k++ ){
                //max
                q[i][mMaxY  ][k]    =  (float)0.0;
                //max-1
                q[i][mMaxY-1][k]    =  (float)(q[i][mMaxY-2][k] - mReInverse*2.0*m_hInverse[_y] * V.v[_v].q[i][mMaxY-2][k]);
                //zero
                q[i][   0   ][k]    =  (float)(q[i][   1   ][k] - mReInverse*2.0*m_hInverse[_y] * V.v[_v].q[i][   2   ][k]);
            }
        }
        //z section
        for( int i=0; i <= mMaxX; i++ ){
            for( int j=0; j <= mMaxY; j++ ){
                //max
                q[i][j][mMaxZ  ]    =  (float)0.0;
                //max-1
                q[i][j][mMaxZ-1]    =  (float)(q[i][j][mMaxZ-2] - mReInverse*2.0*m_hInverse[_z] * V.v[_w].q[i][j][mMaxZ-2]);
                //zero
                q[i][j][   0   ]    =  (float)(q[i][j][   1   ] - mReInverse*2.0*m_hInverse[_z] * V.v[_w].q[i][j][   2   ]);
            }
        }
    }
}

/******************************************************************************
  
  Copyright(C) 2007 DeepDigital Co.,Ltd. All Rights Reserved.
  
 ******************************************************************************/

package FDM;

import math.ValueVector3d;

public class MatrixV3D extends BasicInfo {

    MatrixV[]           v;
    
    public MatrixV3D() {
        
        v = new MatrixV[3];
        v[_u] = new MatrixV( _u );
        v[_v] = new MatrixV( _v );
        v[_w] = new MatrixV( _w );
        
    }

    private void setBoundaryCondition(){
        for( int vComponent=_x; vComponent<=_z; vComponent++ ){
            v[vComponent].setBoundaryCondition();
        }
    }
    
    public void setPosition( int i, int j, int k ){
        Grid.setPosition( i,j,k );
    }
    
    public double v( int component ){
        return v[component].q();
    }

    public double uRep( int currentComponent ){
        return v[_u].represent( currentComponent );
    }
    public double vRep( int currentComponent ){
        return v[_v].represent( currentComponent );
    }
    public double wRep( int currentComponent ){
        return v[_w].represent( currentComponent );
    }

    public double uRepForwardDifference( int direction ){
        return ( v[_u].represent( direction, +1 ) - v[_u].represent( direction, 0 ) )
                * m_hInverse[ direction ] * 0.5;
    }
    public double vRepForwardDifference( int direction ){
        return ( v[_v].represent( direction, +1 ) - v[_v].represent( direction, 0 ) )
                * m_hInverse[ direction ] * 0.5;
    }
    public double wRepForwardDifference( int direction ){
        return ( v[_w].represent( direction, +1 ) - v[_w].represent( direction, 0 ) )
                * m_hInverse[ direction ] * 0.5;
    }

    public void calc( MatrixP p, MatrixT T ){
        setBoundaryCondition();
        for( int i=1; i <= mMaxX-2; i++ ){
            for( int j=1; j <= mMaxY-2; j++ ){
                    for( int k=1; k <= mMaxZ-2; k++ ){
                    setPosition( i,j,k );
                    v[_u].setValue( calc( p, _u ) );
                    v[_v].setValue( calc( p, _v ) );
                    v[_w].setValue( calc( p, _w ) + m_dt * mGrPerRe2 * T.q() );
                }
            }
        }
    }

    private double calc( MatrixP p, int component ){
        return ( v(component)
                 +  m_dt * ( -( uRep( component ) * v[ component ].centralDifference1(_x)
                               +vRep( component ) * v[ component ].centralDifference1(_y)
                               +wRep( component ) * v[ component ].centralDifference1(_z) )
                             + mReInverse * ( v[ component ].centralDifference2(_x)
                                             +v[ component ].centralDifference2(_y)
                                             +v[ component ].centralDifference2(_z) )
                             - p.backwardDifference( component )  
                            )
                );
    }

    public ValueVector3d getValue( int i, int j, int k ){
        return new ValueVector3d( v[_u].q[i][j][k],v[_v].q[i][j][k],v[_w].q[i][j][k] );
    }
}

/******************************************************************************
  
  Copyright(C) 2007 DeepDigital Co.,Ltd. All Rights Reserved.
  
 ******************************************************************************/

package FDM;

import graphic3D.RenderObjectForFlow;

public class Solver3D {

    MatrixV3D           V;
    MatrixP             p;
    MatrixT             T;

    RenderObjectForFlow     mRenderObject;
    
    final int   RENDERING_INTERVAL;
    final int   RENDERING_WAIT;
    
    public Solver3D(int maxX, int maxY, int maxZ,
                   double re, double pr, double gr,
                   RenderObjectForFlow renderObject,
                   int caseNo ) {

        mRenderObject = renderObject;
        BasicInfo.setBasicInfo( maxX,maxY,maxZ,re,pr,gr,caseNo );
        
        V = new MatrixV3D();
        p = new MatrixP();
        T = new MatrixT();
        
        RENDERING_INTERVAL = 2;
        RENDERING_WAIT = (int)Math.floor((maxX)*(maxY)*(maxZ)*0.1);
    }

    public void calc(){
        Thread mainThread = Thread.currentThread();

        System.out.println( "start"  );
        mRenderObject.zeroClear();
        for( int i=0; i<=100000000; i++ ){
            long startMilisec = System.currentTimeMillis();
            System.out.println( "loop" + i );

            p.calc( V );
            V.calc( p,T );
            T.calcAll( V );

            if( Math.IEEEremainder( i, RENDERING_INTERVAL ) == 0 ){
                mRenderObject.rendering(V,T);
                try{
                    mainThread.sleep(RENDERING_WAIT);
                } catch( InterruptedException e ){}
            }
            System.out.println( "time" +  ( System.currentTimeMillis() -startMilisec ) );
        }
    }
}



/******************************************************************************
  
  Copyright(C) 2007 DeepDigital Co.,Ltd. All Rights Reserved.
  
 ******************************************************************************/

import FDM.Solver3D;
import graphic3D.Applet3D;
import graphic3D.RenderObjectForFlow;

public class Flow {

    public static void main(String[] args) {

        int caseNo = 2;
//      final int maxNo = 22;
        final int maxNo = 14;
        Applet3D app3D = new Applet3D( new RenderObjectForFlow( maxNo,maxNo,maxNo,caseNo ) );

        Solver3D solver;
        switch( caseNo ){
        case 1:
            solver = new Solver3D( maxNo,maxNo,maxNo,
                    100,1.0,0.0,
                    (RenderObjectForFlow)app3D.getRenderObject(), caseNo );
            break;
        case 2:
            solver = new Solver3D( maxNo,maxNo,maxNo,
                    350,0.7,10000.0,
                    (RenderObjectForFlow)app3D.getRenderObject(), caseNo );
            break;
        case 3:
            solver = new Solver3D( maxNo,maxNo,maxNo,
                    500,7.0,10000.0,
                    (RenderObjectForFlow)app3D.getRenderObject(), caseNo );
            break;
        default:
            solver = new Solver3D( maxNo,maxNo,maxNo,
                    100,1.0,0.0,
                    (RenderObjectForFlow)app3D.getRenderObject(), caseNo );
            break;
        }
        
        solver.calc();
    }
}

またパッケージ構成は以下の様になっています。

HOME > 流体解析

Copyright(C) Since 2007 DeepDigital Co.,Ltd. All Rights Reserved.