Object Oriented Programming

The picture below shows how the classes inheritance.

Grid class has a 3D grid and take charge of one physical quantity.

MatrixBase inherited Grid class. It calculate difference of physical quantity which the Grid class has.

MatrixV inherited MatrixBase and it's specialized in handling velocity. And MatrixP is specialized in handling pressure, and MatrixT is for temperature.

And here, MatrixV3D has three of MatrixV as , and as its members.

And Solver3D has MatrixV3D for velocity, MatrixP for pressure and MatrixT for temperature as its members.

But it wasn't like this at the beginning. Originally, I thought like having every physical quantity I need into one cell and expand it three dimensionally. But when I tried to calculate difference, I had to write similar codes for every quantities. That made me realized something was very wrong. So I abandoned the idea of having every thing into one cell, and had only single quantity in a three dimensional world, and put them parallel as much as I need. Even vector value like velocity, each three component of velocity can be handled as a scalar value like pressure. And it was surprising for myself when I wrote a code in this way that the code which calculate from just look like the equation itself.


The code for this equation is,

    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 ) 
                 )

It might be possible to implement even nabla itself.

But inside MatrixV, I had to declare it is or or . Also I had to distinguish directions to implement. I must admit it is a hard coding. I guess there must be a way to implement it with ideas like it's parallel or perpendicular to the other, or something like that, instead of implementing things for each directions separately.

I implemented temperature by Boussinesq approximation with MatrixT which inherited MatrixBase. In the same way, I think it's not difficult to implement density as a parameter.
Talking about calculating difference, It's executed only in MatrixBase. It's the only section which do it in the entire code, since one of the cardinal rule of OOP is, "Don't write the same thing twice". Currently the accuracy of the difference of this code is pretty poor, but I can improve it by improving only MatrixBase. This is cool and that's why I like OOP.

One more thing, Grid class has a access to the physical quantity. Accessing is done with current address(i, j, k) plus offset(0, 0, +1). And the current address is memorized as static variables. So, as long as a class inherited this class, the current address of all the classes are synchronized when I move the static variables. I can't mistake the address when I calculate things.

Here I show the code related to the calculation.


/******************************************************************************
 
  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){
                            //cavity velocity
                            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();
    }
}

Here is the package diagram.

Home > Fluid Solver

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