加上迭代步后投影（按照SIMPLE想法）

/*the Sum and the DATE:2012-8-18*/
//本程序先测试实现纯流体的流动，
//边界条件：左边界为一抛物线，右边界也是和左边界相同的抛物线
//边界条件：上下边界是固定边界u=v=0
//内部节点初始值全部为0
#include<iostream>
#include<math.h> 
#include <iomanip>
#include <fstream>

#include "mathoperation.h"
#include "Solve-ustar.h"
#include "Solve-vstar.h"
#include "Solver-Pressure.h"
#include "SUPie.h"
#include "SVPie.h"
#include "eta.h"
#include"SVelPlusOne.h"
using namespace std;

void (*Sol_Press)(int ,int ,double  ,double ,double ,double ** ,double ** ,double ** ,double **,double **); //定义一个函数指针变量 Sol_Press

void (*Sol_ustar)(double ** ,double ** ,double ** ,int , int ); //定义一个函数指针变量 Sol_ustar
void (*Sol_vstar)(double ** ,double ** ,double ** ,int , int ); //定义一个函数指针变量 Sol_vstar

void (*Sol_upie)(double **ustar,double **upie,double **Pres,int N,int M ); //定义一个函数指针变量 Sol_upie
void (*Sol_vpie)(double **ustar,double **upie,double **Pres,int N,int M ); //定义一个函数指针变量 Sol_vpie

void (*Sol_VelPlusOne)(double **u,double **upie,double **v,double **vpie,double **eta,int N,int M);
//double Length;
//double Wide;
//double delta_t;
int main()
{
    int N;
    cout<<"Please the row"<<endl;//输入行
    cin>>N;
    int M;
    cout<<"Please the lie"<<endl;//输入列
    cin>>M;
    double Length=30.0;   //求解区域的长
    //double Length;
    double Wide=15.0;      //求解区域的宽
    //double Wide;
    double delta_x=Length/N;
    double delta_y=Wide/M;
    double delta_t=0.001;
    //double delta_t;
    double Re=1;

    double **Pres=new double *[N+1];//Pres代表压力 
    for(int k=0;k<N+1;k++)
        Pres[k]= new double [M+1];
    for(int i=0;i<N+1;i++)
        for(int j=0;j<M+1;j++)
            Pres[i][j]=0;

    double **ustar=new double *[N+1];
    for(int k=0;k<N+1;k++)
        ustar[k]=new double [M+1];
    for(int i=0;i<N+1;i++)
        for(int j=0;j<M+1;j++)
            ustar[i][j]=0;
    cout<<endl;
    double **vstar=new double *[N+1];
    for(int k=0;k<N+1;k++)
        vstar[k]=new double [M+1];
    for(int i=0;i<N+1;i++)
        for(int j=0;j<M+1;j++)
            vstar[i][j]=0;
    cout<<endl;
    double **u=new double *[N+1];
    for(int k=0;k<N+1;k++)
        u[k]=new double [M+1];
    for(int i=0;i<N+1;i++)
        for(int j=0;j<M+1;j++)
            u[i][j]=0;
    /*---------初值条件---边界--------------*/
    for(int i=0;i<1;i++)
    {
        for(int j=0;j<=M;j++)
        {
            u[i][j]=-(delta_y*j-15)*(delta_y*j-0)*8.0/225;
        }
    }
    for(int i=N;i<N+1;i++)
    {
        for(int j=0;j<=M;j++)
        {
            u[i][j]=-(delta_y*j-15)*(delta_y*j-0)*8.0/225;//抛物线最大值为2
        }
    }
    for(int i=0;i<N+1;i++)
    {
        for(int j=0;j<1;j++)
        {
            u[i][j]=0;
        }
    }
    for(int i=0;i<N+1;i++)
    {
        for(int j=M;j<M+1;j++)
        {
            u[i][j]=0;
        }
    }
    /*------------------------------------------------------------------------------------------------------*/
    cout<<"-----Output the u ------"<<endl;
    for(int i=0;i<N+1;i++)
    {
        for(int j=0;j<M+1;j++)
        {
            cout<<" "<<u[i][j];
        }
    cout<<endl;
    }
    /*---------------------给v分配空间并赋初值---------------------------*/
    double **v=new double *[N+1];
    for(int k=0;k<N+1;k++)
        v[k]=new double [M+1];
    for(int i=0;i<N+1;i++)
        for(int j=0;j<M+1;j++)
            v[i][j]=0;
    cout<<endl;
    /*给u'开辟空间并赋初值*/
    double **upie=new double *[N+1];
    for(int k=0;k<N+1;k++)
        upie[k]= new double [M+1];
    for(int i=0;i<N+1;i++)
        for(int j=0;j<M+1;j++)
            upie[i][j]=0;

    /*给v'开辟空间并赋初值*/
    double **vpie=new double *[N+1];
    for(int k=0;k<N+1;k++)
        vpie[k]= new double [M+1];
    for(int i=0;i<N+1;i++)
        for(int j=0;j<M+1;j++)
            vpie[i][j]=0;
/*---------初值条件---边界--------u'------*/
    /*for(int i=0;i<1;i++)
    {
        for(int j=0;j<=M;j++)
        {
            upie[i][j]=2;
        }
    }
    for(int i=N;i<N+1;i++)
    {
        for(int j=0;j<=M;j++)
        {
            upie[i][j]=2;
        }
    }
    for(int i=0;i<N+1;i++)
    {
        for(int j=0;j<1;j++)
        {
            upie[i][j]=2;
        }
    }
    for(int i=0;i<N+1;i++)
    {
        for(int j=M;j<M+1;j++)
        {
            upie[i][j]=2;
        }
    }*/
        //std::ofstream outvel;
        //outvel.open("F:\\fluid\\cC-V-IBM\\vel.txt");    
        //cout<<"the N+1 Setp  velocity!"<<endl;
        //for(int j=M;j>=0;j--)
        //{
        //    for(int  i=0;i<N+1;i++)
        //    {
        //        //cout<<" "<<u[i][j]; 
        //        outvel<<i*delta_x<<" "<<j*delta_y<<" "<<u[i][j]<<" "<<v[i][j]<<" "<<Pres[i][j];
        //        outvel<<endl;
        //    }
        //    cout<<endl;
        //    outvel<<endl;
        //}
        //outvel.close();
    double **eta=new double *[N];
    for(int k=0;k<N;k++)
        eta[k]= new double [M];
    for(int i=0;i<N;i++)
        for(int j=0;j<M;j++)
            eta[i][j]=10;//本来eta[i][j]在0~1之间，这里取10为了防止错误，便于检查
    SolEta((double **)eta, N, M);//更新eta
    for(int k=1;k<2;k++)
    {
        /*---------初值条件---边界--------------*/
        for(int i=0;i<1;i++)
        {
            for(int j=0;j<=M;j++)
            {
                upie[i][j]=-(delta_y*j-15)*(delta_y*j-0)*8.0/225;
            }
        }
        for(int i=N;i<N+1;i++)
        {
            for(int j=0;j<=M;j++)
            {
                upie[i][j]=-(delta_y*j-15)*(delta_y*j-0)*8.0/225;//抛物线最大值为2
            }
        }
        for(int i=0;i<N+1;i++)
        {
            for(int j=0;j<1;j++)
            {
                upie[i][j]=0;
            }
        }
        for(int i=0;i<N+1;i++)
        {
            for(int j=M;j<M+1;j++)
            {
                upie[i][j]=0;
            }
        }
        cout<<"---the "<<k<<"   step---"<<endl;
        Sol_ustar=CUEqution; //将函数指针指向函数CU;
        (*Sol_ustar)((double **) ustar,(double **) u,(double **) v,N,M); //调用CU函数
        
        Sol_vstar=CVEqution; //将函数指针指向函数CV;
        (*Sol_vstar)((double **)vstar,(double **)u,(double **)v,N,M); //调用CV函数
        double epsilon=0.01;//速度修正的精度
        int NUM=0;//压力修正累积的步数
        int Iterstep=20;//压力修正过程中最大的迭代步数
        do{
            Sol_Press=Sol_P; //将函数指针指向函数压力;
            (*Sol_Press)(N,M,delta_x,delta_y,delta_t,(double **) ustar,(double **) vstar,(double **) upie,(double **) vpie,(double **) Pres); //调用CV函数

            Sol_upie=SolUPie; //将函数指针指向函数SolUPie;
            (*Sol_upie)((double **)ustar,(double **)upie,(double **)Pres, N, M ); //调用CU函数

            Sol_vpie=SolVPie; //将函数指针指向函数SolVPie;
            (*Sol_vpie)((double **)vstar,(double **)vpie,(double **)Pres, N, M ); //调用CV函数

            //Sol_VelPlusOne=VelPlusOne; //将函数指针指向函数Sol-Vel在下一步的值;
            //(*Sol_VelPlusOne)((double **)u,(double **)upie,(double **)v,(double **)vpie,(double **)eta, N, M);
            NUM++;
            std::ofstream out;
            out.open("F:\\fluid\\cC-V-IBM\\vel.txt");    
            cout<<"the N+1 Setp  velocity!"<<endl;
            for(int j=M;j>=0;j--)
            {
                for(int  i=0;i<N+1;i++)
                {
                    cout<<" "<<u[i][j]; 
                    out<<i*delta_x<<" "<<j*delta_y<<" "<<u[i][j]<<" "<<v[i][j]<<" "<<upie[i][j]<<" "<<vpie[i][j]<<" "<<Pres[i][j];
                    out<<endl;
                }
                cout<<endl;
                out<<endl;
            }
            out.close();
            cout<<"Matnorm2(u,ustar,N+1,M+1)="<<Matnorm2(u,ustar,N+1,M+1)<<endl;
            cout<<"Matnorm2(v,vstar,N+1,M+1)="<<Matnorm2(v,vstar,N+1,M+1)<<endl;
        }while((Matnorm2(u,ustar,N+1,M+1)>epsilon||Matnorm2(v,vstar,N+1,M+1)>epsilon)&&NUM<Iterstep);
        //当压力修正步数不够设置值且速度二范数大于精度要求，做修正
        cout<<"NUM="<<NUM;
    }
    /*std::ofstream out;
    out.open("F:\\fluid\\cC-V-IBM\\vel.txt");    
    cout<<"the N+1 Setp  velocity!"<<endl;
    for(int j=M;j>=0;j--)
    {
        for(int  i=0;i<N+1;i++)
        {
            cout<<" "<<u[i][j]; 
            out<<i*delta_x<<" "<<j*delta_y<<" "<<u[i][j]<<" "<<v[i][j]<<" "<<upie[i][j]<<" "<<vpie[i][j]<<" "<<Pres[i][j];
            out<<endl;
        }
        cout<<endl;
        out<<endl;
    }
    out.close();*/
    
    for(int t=0;t<N;t++)
        delete [] eta[t];
    delete [] eta;

    for(int t=0;t<N+1;t++)
        delete [] Pres[t];
    delete [] Pres;

    for(int t=0;t<N+1;t++)
        delete []u[t];
    delete []u;
    for(int t=0;t<N+1;t++)
        delete []v[t];
    delete []v;
    for(int t=0;t<N+1;t++)
        delete []ustar[t];
    delete []ustar;
    for(int t=0;t<N+1;t++)
        delete []vstar[t];
    delete []vstar;

    /*撤销u'空间*/
    for(int t=0;t<N+1;t++)
        delete [] upie[t];
    delete [] upie;

    /*撤销v'空间*/
    for(int t=0;t<N+1;t++)
        delete [] vpie[t];
    delete [] vpie;

    return 0;
}